CC BY
16
Russian Journal of Earth Sciences
Estimation of Groundwater Recharge from Rainfall for Arid Coastal Plain of Ninh Thuan Province, Vietnam
Nguyen Van Hoang1, Trinh Hoai Thu2, Renat Shakirov3, Tran Thi Thuy Huong2, and
Nadezhda Syrbu*'3
1 Institute of Geological Sciences VAST, Hanoi, Vietnam 2Institute ofMarine Geology and Geophysics VAST, Hanoi, Vietnam 3Il'ichev Pacific Oceanological Institute FEB RAS, Vladivosstok, Russia
Received 18 June 2020; accepted 25 August 2021; published 17 February 2022.
Estimation of groundwater recharge from rainfall is a key factor for determining groundwater resources in water development and management for supporting sustainable socio-economic development, especially for arid areas. The paper presents finite element modeling in the simulation of moisture transfer in unsaturated soils through the relationship between soil moisture, soil suction, unsaturated permeability and soil-moisture dispersivity. Those parameters required for soil moisture transfer are derived from the soil-water characteristic curve functions. Element sizes and time steps used in the modelling have been selected based on a detailed analysis of numerical simulation errors. The methodology had been applied to arid coastal plain area of Ninh Thuan province, Vietnam. Five subsurface soil types in the study area have been collected and analyzed, for saturated permeability, porosity, saturated soil water content, field moisture content, etc. Hourly rainfall data of the years 2014-2018 have been analyzed and grouped into different-duration rainfall events (1-hour, 2-hour, 3-hour and so on). The different rainfall durations and depths of rainfall events and temporal infiltration determined by the moisture transfer modelling have allowed determining the groundwater recharge from the rainfall data. The results show that during the rainy months from May to December 2014-2018, the groundwater recharge from the rainfall is very varying through the modeled soil profiles, from 0.280 m (silty clay) to 0.470 m (sandy silt), which is equivalent to 33.3%-55.2% of the rainfall depth during May-December. Lower infiltration in silty clay is due to low permeability and in the sand is due to low suction, and higher infiltration in silt and sandy silt is thanks to their higher moisture dispersivity. On average, in terms of annual rainfall and soil properties, the average infiltration during May-December is 0.380 m which is equivalent to 44.9% of the rainfall depth, which is about 289 x 106 m3 of rainwater infiltrated into the Quaternary aquifer over 760 km2 of the coastal plain of Ninh Thuan province. The results would be very useful for effective water resources development and management in a given specific hydrogeological condition for such a severe drought area where water is extremely essential.
Keywords: tropical savannah climate, drought, moisture transfer, finite element (FE), higher-order element function
Citation: Nguyen Van Hoang, Trinh Hoai Thu, Renat Shakirov, Tran Thi Thuy Huong, and Nadezhda Syrbu, (2022), Estimation of Groundwater Recharge from Rainfall for Arid Coastal Plain of Ninh Thuan Province, Vietnam, Russ. J. Earth. Sci., Vol. 22, ES1001, 10.2205/2022ES000775.
1 Introduction
The South central region of Vietnam, includ-
ing Ninh Thuan province, has been heavily suffer-
ing from drought and salinization. In accordance with the Koppen-Geiger climate classification sys-
tem and a 1-kilometer resolution Koppen-Geiger climate map Beck et al. [2018] (Figure 1), Ninh Thuan province has a tropical savanna climate, i.e.,
severe dry season, and prevailing drought conditions during the year. Ninh Thuan province is one of the provinces which are most affected by
'Corresponding author: syrbu@poi.dvo.ru
the drought. Severe drought in recent years (2012 to 2016) led to reduced water availability for irrigated agriculture, especially for higher-value crops such as black pepper, coffee, dragon and grapefruits. In 2016, there was 5770 ha of agricultural land in Ninh Thuan province, which is about 30% of agricultural land, could not be planted due to water shortage and 5500 households were suffering from a shortage of domestic water [VNCSC-NDP&C, 2016]. In the context of the drought situation, the Asian development bank had financed a project to improve supply and water efficiency in drought-affected provinces in the central coast
Figure 1: Map of the study area in Köppen-Geiger climate classification map of Vietnam
and central highlands of Vietnam, including Ninh Thuan province, one main objective of which is a qualitative groundwater resource assessment since the increase of the frequency and intensity of extreme weather conditions had resulted in increasing water scarcity, and the farmers in these regions are heavily relying on groundwater for domestic use and irrigating the high-value crops [Jacobs, 2017]. which could contribute to the sustainable socio-economic development of the regions.
The coastal plain area of Ninh Thuan province has annual rainfall (P) of 750 mm - 900 mm and annual potential evapotranspiration (ETP) of 1840 mm [Hai, 2015]. Sam and Vuong [2008] stated that the analysis and interpretation of different drought indices and drought frequencies for Ninh Thuan province had shown that the aridity index (IA) proposed by UNEP [1997] for susceptibility analysis of desertification, the ratio between P and ETP, is the most suitable index for aridity study of the coastal plain area of the province. The analysis results allowed the authors to conclude that the coastal plain area of Ninh Thuan province belongs to the arid area (from low aridity to high aridity) with IA from 1 to less than 0.25. The annual average aridity index for the coastal plain area of Ninh Thuan province is 0.45 which is corresponding to susceptibility to desertification as classified by UNEP [1997]. Climate change exacerbates the desertification of Ninh Thuan province and makes its conditions more severe: the area of desertification is 41.02 ha which is 12.21% of the province
area) [Loan, 2018]. In the context of the aridity of the coastal plain area of Ninh Thuan province, effective water resources development and management are very important for the area. It is well known that in arid and semi-arid environments, groundwater recharge is the most critical to the sustainable use of water. However, so far a well-based groundwater recharge assessment for the Southern coastal plain region of Vietnam in general and the coastal plain area of Ninh Thuan province, in particular, has not been implemented.
Groundwater recharge from precipitation needs to be determined. This study attempts to estimate the groundwater recharge from the precipitation in about 760 km2 of the coastal plain area of Ninh Thuan province (Figure 1).
Many different methods can be applied to access groundwater recharge dependently upon given hy-drogeological and climatic situations. The most common methods applicable for arid and semiarid lands have been presented and reviewed by Kinzelbach et al. [2002]. The authors described and classified the most important methods, evaluated the advantages and disadvantages, and judged the accuracy and reliability of the methods. Four categories of the methods have been grouped as 1) Direct measurements; 2) Water balance methods; 3) Darcyan methods, and 4) Tracer methods. The accuracy rating of recharge estimate by the methods has been classed as 1) Class 1: within a factor of 2; 2) Class 2: within a factor of 5; and 3) Class 3: within a factor of 10 or more [Kinzel-
bach et al., 2002]. Since different recharge estimate methods have different accuracy and reliability, the authors have suggested that a successful estimation of groundwater recharge would have resulted from the utilization of several independent methods. Within this study, an unsaturated moisture transfer by Richards' equation [Richards, 1931] which requires a more complex model with even more unknown parameters is applied. Despite a low accuracy rating of recharge estimate by unsaturated moisture transfer by Richards' equation as concluded by Kinzelbach et al. [2002], a sophisticated and well-based theoretical unsatu-rated moisture transfer with soil-water characteristic curve (SWCC) parameters is used in this work, the results of which would prove that the estimate accuracy highly depends upon the accuracy and reliability of the input data of the model.
2 Hydrogeological conditions of Ninh Thuan coastal plain
The coastal plain of Ninh Thuan province has an area of about 760 km2, including about 200 km2 of coastal sand dune [Tu et al., 2016] (Figure 1) and consists of Quaternary deposits, which are divided into two aquifers, the Holocene aquifer (qh) and the Pleistocene aquifer (qp), and an undivided Quaternary aquifer. The Holocene aquifer (qh) and the Pleistocene aquifer (qp) are the main aquifers having a potential water resource for socio-economic use, the features of which may be referred to in the work of Tu et al. [2016].
The Holocene aquifer is the first aquifer from the ground surface and consists of silt and sand with grit and gravel. The aquifer has a very thin thickness from less than a meter to about 10 meters, on average 5 meters. The aquifer is unconfined and has a water level depth from few decimeters to 6 meters, on average 2.5 meters. About 57% of the Holocene aquifer area has high total dissolved solids (TDS) more than 1.02 g/l with the water chemical form of Chloride-Sodium or ChlorideSulfate or Sodium-Calcium water. For the area of freshwater distribution, in general, the aquifer water has an adequate quality for domestic use. However, in some places Nitrate concentration is much higher than 15 mg/l and Nitrogen and Sulfate concentration is higher than 400 mg/l as specified by the Vietnam national technical regulation on groundwater quality [MONRE, 2015].
The Pleistocene aquifer has a total distributed area of about 362 km2. There are 152 km2 of the aquifer, underlying Holocene aquifer, while the remaining 210 km2 are exposed to the ground surface where the Holocene aquifer is absent. The Pleistocene aquifer consists of silt and sand with calcareous grit in the upper part and of medium
to coarse quartz sand with gravel in the lower part. The aquifer has a variable thickness from less than a meter to about 44 meters, on average 7 meters. The groundwater level is from few decimeters above the ground surface to 6 meters below the ground surface, on average 2.45 meters below the ground surface. About 52% of the aquifer area has high total dissolved solids (TDS) from 1.42 g/l to 19 g/l with the water chemical form of Chloride-Sodium or Chloride-Bicarbonate-Sodium water. For the area of freshwater distribution, in general, the aquifer water has an adequate quality for domestic use. However, in some places, nitrate concentration is much higher than 15 mg/l of Nitrogen as specified by the Vietnam national technical regulation on ground-water quality [MONRE, 2015].
There is no impermeable layer in between the Holocene and Pleistocene aquifers. Therefore, the two aquifers in combination form one hydraulic aquifer with a single water level.
In the study area (plains of Ninh Thuan province), underneath the Pleistocene aquifer in the Upper Cretaceous effusive formation consisting of dacite, ryodacite, felsite andesitodacite and their tuff. The drilling data showed its maximal thickness of about 67 meters. The aquifer is low permeable: two pumping tests in the aquifer gave transmissivity values of 1.18 m2/d and 1.28 m2/d. The hydrogeological map of the study area is presented in Figure 2 and a typical hydrogeological cross-section in the middle of the plain along line AB from north-west to south-east is presented in Figure 3.
3 Materials and methods
Soil profile from the surface to the top of the aquifer in the study area had been determined from various previous studies and projects. The hydraulic conductivity parameter of the soil had been collected and analyzed to classify into different representative ranges. The number of the hydraulic conductivity parameter data is 52. As the SWCC study is not common in the field of hy-drogeology and groundwater science in Vietnam, the SWCC data of the soil are not available for the study area. The SWCC of the soils is used as provided by Geo-Slope [2015] for silty clay, silt, silty sand and sand. An additional sandy silt SWCC between silt and silty sand. The hourly monitoring rainfall data in the study area in the period 2014-2019 have been provided by the Hydro-Meteorological Information and Data Center under the Vietnam Meteorological and Hydrologi-cal. The methods used in the study mainly are the statistical analysis for obtaining the representative parameter values, the parameter identification by
Figure 2: Hydrogeological map
Figure 3: Hydrogeological cross-section along with line AB
the best fitting algorithm, soil moisture transport parameter calculation by various formulas of the SWCC theory, Galerkin FEM with higher-order element functions for modelling the soil moisture transfer with the calculated soil SWCC-based parameters.
4 SWCC-based modelling of unsaturated moisture parameter
4.1 Equation of soil moisture transfer in unsat-urated soil
The equation describing the moisture transfer in unsaturated soil with an assumption that air does not move is as follows Polubarinova-Kochina [1977]
dux + duy + dUz = _ dËw
dx dy dz dt
VX = -K (6W ) -g ;
vy = -K (eu ) ^ ; dh
dy dZ
Vz = -K e ) ^ = -K (ew )!-£ + !
Equations (1) and (2) give:
ddw dt
d Ik(eu) dh
dx\ y dx
d Ik(eu) dh
dy\ y dy
+ d_ / KiM dh
dz\ y dz
dew dt
d (k(ew)dM + ±Ik(M+
dx\ y dx ] dy\ y dy
+ 1(K (ew )
dz\ y
dty dz dz dz
Which is the Richards equation in terms of suction pressure [Philip, 1969].
With water density equal to one and suction pressure y as a function of VWWC eu (3) is calculated in vertical direction z:
de„
dt dz
= Î\K (eu )
dy deu dew dz
(1)
+ K (eu ) . (4)
The component D(6w) = K(6w) is called soil-moisture diffusivity [Philip, 1957] (as cited by [Gilding, 1991]) with unit L2T-1, and then (4) has the following form:
where: eu is volumetric water content (VWC), t is a time, vx; Vy and vz are moisture transfer rates in x, y and z direction, respectively:
deu dt
d_ dz
d (eu )
deu dz
dK (eu ) dz .
(5)
(2)
Initial condition: moistrure content is known at time t = 10:
eu = e u (z).
Boundary conditions:
Dirichlette condition:
eu = eu
Neumann condition:
on ru
de^
dn
= g on ru
(3)
where: h = y + z; h is a total head; ^ is soil suction; y is a water density, and K(6w) is hydraulic conductivity of unsaturated soil.
4.2 Finite element modelling of soil moisture transfer in unsaturated soil
Finite element (FE) method [Zienkiewics and Morgan] can be applied in solving (5) with the
+
use of quadratic element shape functions for a more accurate numerical solution than the use of linear shape functions. For simplicity let us use Dz instead of Dz(0w). Applying Green lemma to (5) with the approximation of water content 0W « 0W = Lm=1 9w,mNm gives:
L (D# )
=-J. (
dK (0W ) dz
Dzdtdw)dz + Jr DdïW^* T+
d0w
dK(0w) d0w
dz dt
The boundary integration jr is present only for the elements having sides on boundaries Tgw:
Ja (
Dz^dz + f (o^WN.ldr*
J, (
D9Nm9Wl0 +
D? —--—0m \dz +
z z
f (qcWiNz)dr+
q
dK(0w) 90,
dz dt
K = I w/Sm W)dz;
JL( '
z z
f =
dK (0w )
dz
+ f (qWlNz) dr. Jra
d0 dt
i at„
+ AK\0n +
- zn+(1 -A)K
Dn+1
+ (1 - X)fn + Xfn+1,
ty was proposed by Gardner [1964] [Fredlund et al., 2011]
0w = —s—. w 1+atyn
Where: 0w is soil moisture content; 0g is saturated soil moisture content; ty is soil suction; a and n are fitting soil parameters associated with the soil-water characteristic curve (SWCC).
Since that time, there are numerous bestfit equations that have been proposed for the SWCC of unsaturated soil [Brooks and Corey, 1964; Brutsaert, 1966; Fredlund and Xing, 1994; van Genuchten, 1980; McKee and Bumb, 1984, 1987]. Leong and Rahardjo [1997] showed that the equations proposed by van Genuchten [1980] and Fredlund and Xing [1994] give more flexibility to best-fit the measured data. Fredlund and Xing [1994] proposed an equation expressed as:
dK (0w) d0w
+ ~dzr = -ft. (6)
0w = C (v)
0S
[ln{ e + ( V/a)n\X
Putting 0w ~ 0W = Lm=1 0w,mNm into (6) results in:
where C(ty) = 1 - ln(1+ty/tyr) is a correction fac-
ln(1+106/tyr)
tor, tyr is the suction corresponding to the soil residual water content; a is the fitting parameter related to the air-entry value of the soil (kPa), n is the fitting parameter related to the slope of the SWCC, m is the fitting parameter related to the soil residual water content, e = 2.71828, ty is soil suction (kPa). Leong and Rahardjo [1997] concluded that C(ty) can be assumed to be unity without affecting the initial portion of the SWCC, which helps to reduce the number of parameters in the equation, i.e.:
K0 +^= f. (7)
0W =
0S
Applying finite element method to (7) for the time domain, the following is [Zienkiewics and Morgan]:
[ln{ e + (ty/a)"}]m
The unsaturated permeability K(0w) of the soil can be determined through the relative coefficient of permeability kr(0) [Fredlund et al., 1994]:
K(0^) = K x kr(0);
kr (0) =
where A = 0 ^ 1.
Zienkiewics and Morgan in very detail showed that the schemes with A > 0.5 are always unconditionally stable for any values of time step At.
4.3 Soil moisture transfer parameters in accordance with the soil-water characteristic curve
The first equation describing the relationship between soil moisture content 0w and soil suction
IM) dcIL
0S 0S - Z
V2(Z)
dZ, (8)
where: K is the hydraulic conductivity of saturated soil; 0 is the effective VWC defined as 0 = 0w - 0r; 0w is the VWC; 0L is the lowest VWC on experimental SWCC; 0r is residual VWC; 0s is the saturated VWC; Z is a dummy integration variable; Se = (0w - 0r) / (0S - 0r) is the effective saturation.
In accordance to Kunze et al. [1968] [Fredlund et al., 1994] an improved prediction accuracy would have resulted if the complete SWCC is used, namely with the availability of residual VWC 0r where
+
1
L
kr = 0. In this case, the integration in (8) should be carried out in the interval from 0 r to 0
S-
kr ( e) =
Lr I 2(o dzIL,
S es - Z
)er y(Z) i jer f2(z) The soil-moisture diffusivity is:
dy
dZ. (9)
d ( e u) = K ( 0u)
d e u
4.4 Procedure of FE modelling of soil moisture transfer with SWCC-based parameters
Fortran programming code of groundwater modelling by finite element method using linear and higher-order element functions in the study NCCB-DHUD.2012-G/04 [Hoang et al., 2018] was modified for moisture transfer modelling in this work by the flow chart in Figure 4. Several subroutines were programmed to determine all the necessary SWCC-based parameters (relationships between the soil suction, relative permeability, soil-moisture diffusivity and VWC), which are required in the soil moisture modelling by the above-described theoretical fundamentals. The flow chart of the soil moisture transfer FE modelling is presented in Figure 4 along with the required parameters derived from SWCC function.
5 Site condition and estimate of groundwater recharge from precipitation
5.1 Soil physical parameters
Since the groundwater recharge from rainfall occurs through the infiltration via the surface soil layer, the surface soil layer physical parameters need to be identified. The topsoil layer in the plain area of Ninh Thuan province consists of alluvial deposits and has a thickness of more than one meter. This topsoil layer serves as agricultural soil, which has very specific characteristics, one of which is its light texture [Tu et al., 2016]. The pedologic studies of agricultural soils in Ninh Thuan province [DARD, 2017] have shown that the topsoils in the plain area of Ninh Thuan province are from light to medium texture, which is loamy sand and sandy loam with 10-20% of clay, 10-24% of silt and sand, and 10-13% moisture range.
Hydraulic conductivity (K) of the Holocene and Pleistocene was determined by field pumping tests. The results of 38 pumping tests in the Holocene aquifer and 17 pumping tests in the Pleistocene aquifer gave the average, minimal and maximal values of horizontal hydraulic conductivity of 1.86 m/d, 0.17 m/d and 5.15 m/d for Holocene aquifer, and 1.61 m/d, 0.26 m/d
C^^Start
Selection of model domain: 1D vertical z direction or 2D in horizontal x and vertical z
Inputting initial moisture content and moisture dispersion coefficient at step n = 0
Selection of element sizes Axj < 10Dx( 0), Ayi < 10Dy( 0), time step and number of run times
Compilation of FEM mesh: coordinates of each node, numbering of nodes and elements
n = 0
Specifying boundary conditions for each time step n (two vertical sides have a gradient of moisture equal to 0)
Determination of relation between moisture transfer coefficient and moisture. Updating moisture transfer coefficient for all elements for a given moisture contents at time step n:
Kr ( e) = srUs dzii;
■ es es-Z
yW)dz\e yz)dz
u - er; K(eu) = K x kr ( e);
d (ew) = K(eu) df
Building up system of equations for moisture variables at time step n +1 with known moisture at time step n
Solving the system of equations for moisture variable; Calculation of soil water pressure.
II-
n <N
n> N
End
Figure 4: Flow chart of moisture transfer FE modelling with SWCC - based parameters
and 4.79 m/d for Pleistocene aquifer [Tu et al., 2016]. There is no clear distinction between the hydraulic conductivity of the Holocene and Pleistocene aquifers. Statistical analysis has shown that the hydraulic conductivity of the Holocene and Pleistocene together follows a normal distribution
Table 1: SWCC Parameters and Soil Physical Properties of the 5 Modeled Soils
Soii SWCC parameters Saturated VWC Initial VWC Hyd. Conductivity
a n m es в 0 K (m/d)
Silty clay 30 1.68 1.03 0.470 0.13 0.013
Silt 17 1.65 1.25 0.430 0.12 0.058
Sandy silt 13 1.63 1.35 0.405 0.11 0.102
Silty sand 9 1.61 1.45 0.380 0.10 0.146
Sand 2.2 2.00 1.34 0.350 0.09 0.200
Figure 5: The normal distribution of horizontal hydraulic conductivity of Holocene-Pleistocene aquifer
with the mean of 1.55 m/d and standard deviation of 1.10 m/d. Most of the hydraulic conductivity data are less than 2 m/d with an average of 1.02 m/d and a standard deviation of 0.44 m/d (Figure 5). Therefore, the moisture transfer analysis is carried out for five cases of typical hydraulic conductivity values: minimal (0.13 m/d), average minus standard deviation (0.58 m/d), average (1.02 m/d), average plus standard deviation (1.46 m/d) and maximal (2.00 m/d). Since there is no information on vertical hydraulic conductivity, a reference ratio between the horizontal (Kh) and vertical hydraulic conductivity (Kv) of 10 is used to get Kv. Therefore, the values of Kv of the five modelling cases are: 0.013 m/d, 0.058 m/d, 0.102m/d, 0.146 m/d and 0.200 m/d.
Since there are no SWCC data existing for the soils under consideration, the SWCC of the soils is used as provided by Geo-Slope [2015] for silty clay, silt, silty sand and sand. An additional sandy silt SWCC between silt and silty sand was interpolated for the case of average K. The a, n and m parameters of the SWCCs and soil physical properties are presented in Table 1. The SWCCs are presented in Figure 6 and the relative coefficient of permeability determined by Eq. (9) is presented in Figure 7.
5.2 Selection of element size in FE modelling
Finite element solutions have three types of error: 1) discretization error due to the element size:
Figure 6: SWCCs of the modeled soils
the smaller element size the less error; 2) roundoff error due to truncation in computer calculation process; and 3) approximation of the exact mathematical model by FE method [Zienkiewics and Morgan]. The modelling applicants are most interested in discretization error, i.e., in the selection of the element size. For estimating rational element size which ensures a small enough discretization error in our modelling, let us consider the steady-state moisture transfer:
d_ dz
d (вв. ) '-t
+ dK (ву )
+ dz
To select rational element size, preliminary soil-
Figure 7: Relative permeability coefficient of the modeled soils
Table 2: Preliminary Soil-Moisture Diffusivity of Different Modeled Soil Cases
Modeled soils
Saturated hydraulic conductivity (m/d)
Porosity/Saturated Natural VWC VWC
D (6w)(m2/d)
Max Min
Silty clay 0.013 0.470 0.13 0.3431 0.00004
Silt 0.058 0.430 0.12 0.8296 0.00024
Sandy silt 0.102 0.405 0.11 1.1152 0.00035
Silty sand 0.146 0.380 0.10 1.1125 0.00038
Sand 0.200 0.330 0.09 0.1089 0.00080
moisture diffusivity had been determined for the modelling cases. Element size of 0.02 m was used in these simulation cases and soil-moisture diffu-sivity values of the different elements were determined. The minimal and maximal values of soil-moisture diffusivity are given in Table 2.
The less value of soil-moisture diffusivity, the smaller element size is needed to be used as the error estimator for a quadratic element [Zienkiewics and Morgan]:
(he)2 rh\ 12 Jc J
^R d x.
where: Oe is element domain; Rqs is an error within element e; he is element size.
Let us consider a hypothetical case with the ratio between permeability gradient and moisture coefficient value equal to 0.01, i.e., the following equation:
520w =0.1.
dz2
If a quadratic element function is used with element size he, then the error estimate is:
PI\ce =
(he)
he
12 1
(he)2 rhe J0
^R d x —
12
0.12dx — 0.01 ^h
12
xih. (10)
The error estimate by Eq. (10) for different element sizes he and absolute relative maximal VWC error estimate with average minimal VWC of 0.11 are listed in Table 3. Therefore, the element size from 0.01 m to 0.05 m would well meet our modelling accuracy requirement. Regarding the time step, Zienkiewics and Morgan in very details showed for our moisture transfer Eq. (5), the most accurate scheme is the Crank-Nicolson scheme (the error order is O(at2)), and for the forward and backward difference schemes, the time step needs to be:
forward difference scheme : at < (he)2/[6D (0w)]; backward difference scheme : at < (he)2/[2D (0w)],
where the partial differential equation used by Zienkiewics and Morgan is the same as our moisture transfer equation where D(9w) = 1. In our work, even the time Galerkin scheme is used as recommended by Huyakorn and Pinder [1983] with A = 2/3 for convection-diffusion equations. The time step of 6 minutes would be appropriate.
The element size of 0.01 m used in this work is within the minimal range of grid spacing or element size used by other authors, which is from 0.01 m to 0.25 m, where the grid spacing or element size is constant over the model domain length of from 1 m to 10 m [Carrera-Hernández et al., 2012]. Carrera-Hernández et al. [2012] used different elements sizes, both constant and variable, over domain lengths of 2 m, 4 m, 6 m and 12 m for studying the sensitivity of groundwater recharge estimate to boundary conditions and spatial discretization in unsaturated flow modelling. The authors concluded that the variable spatial discretization from 0.001 m with increment by a factor of 1.1 to a maximum value of 0.10 m provided the fastest simulation times and the better accuracy of the numerical modelling.
5.3 Initial profile of soil volumetric water content
Initial VWC of the modeled soil profiles in the model for moisture transfer prediction needs to be determined in the present study. The initial time of the modelling is the beginning of the rainy season, i.e., the profiles of soil VWC need to be determined. The VWC profiles of the unsaturated soil layers may be determined by simulation of the models with specified saturated VWC at the bottom and natural field VWC at the ground surface. The results showed that the water content had reached a quasi-steady state, i.e., the soil water content has not changed after several months. The determined initial VWC of the modeled soil profiles are presented in Figure 8, which are well corresponding to coarse, fine and very fine soils described by Bear [2000] in Figure 9.
Table 3: Error Estimate for Different Element Sizes
he (m)
0.01 0.02 0.025 0.03 0.04 0.05
_abs \\E|_
0.00003 0.00008 0.00011 0.00015 0.00023 0.00032
Absolute relative maximal VWC error estimate (100 x abs ||E||/MinVWC) (%)
0.027 0.073 0.100 0.136 0.209 0.291
Figure 8: Initial field water content and suction profiles for the modeled soils
5.4 Boundary conditions
The first type of boundary condition is assigned to the upper model domain of the ground surface since the moisture transfer model is to be carried out during excessive rainfall events which last continuously for 1 hr, 2 hrs, 3 hrs, etc. That is, the soil moisture content at the ground surface is corresponding to the saturated water content.
The lower model domain is also set up as the first-type boundary condition since it directly contacts the water table. This is also in accordance with many other authors who assigned the boundary as either fixed water table, free drainage (unit gradient), or head-dependent [Carrera-Hernández et al., 2012]. Carrera-Hernández et al. [2012] basing on their simulation results had recommended
a fixed head for the lower boundary as the long-term simulations provide some more realistic soil moisture dynamics.
5.5 Selection of modelling depth of the soil
Suppose that at the initial time t = 0 the initial soil VWC profile is corresponding to the hydrostatic condition (Figure 10a). As water input is applied at the surface, after time t1, t2 and t the soil moisture propagates deeper (Figure 10a). With the assumption that the soil water content changes shapely [Tarboton, 2003], i.e., after time 11, t2 and t the soil moisture is equal to saturation (n) in the depth interval 0-L1, 0 - L2 and 0 - L, and is equal to initial moisture in the depth below L1, L2 and L, respectively (Figure 10b). L is called wetting
Figure 9: Typical water content profiles for different soils
front depth. Several models have been proposed for the estimation of the wetting front depth.
The most popular models are Green and Ampt [1911], Horton [1939] Philip [1957], which fundamentally have no advantages of one over the other [Tarboton, 2003]. By Green and Ampt model, the wetting front depth formulae is [Nie et al., 2017]:
Kt i L
L = -KB (h0 + ^)ln(1 + fe^
where: t is the time from the infiltration beginning; A0 is the difference between saturated VWC and initial VWC; h0 is the ponding depth; tym is the suction head at the wetting front.
The monitoring rainfall data in the analysis period (2014-2019) have shown that the longest rainfall duration is 33 hours in 2018, which would result in a maximal wetting front depth of from 0.741 m (for silty clay, the lowest permeable soil under consideration) to 1.088 m (for sand, the highest permeable soil under consideration). Therefore, the modeled length is selected to be from a maximal value of 1.1 m (for sand) to a minimal value of 0.5 m (for silty clay).
Figure 10: Illustration of sharp wetting front: a-Real water content distribution; b-Assumed shape water content distribution
5.6 Simulation results
The simulation results for the modeled soils, water content profiles for the initial field condition are presented in paragraph 4.3, and soil VWC and suctions along the soil columns are presented in Figure 8. The temporal soil VWC along with the depth for the five soils are presented in Figure 11, and the accumulative rainwater infiltration is presented in Figure 12 and Table 4
5.7 Estimation of potential groundwater recharge due to rainwater infiltration during May-November 2014-2018
The rainy season in the study area is from May to December with rainfall, concentrated during September-November. The average total rainfall is 0.811 m, with 95% during May to December and 58% during September-November for 1980-2014 [Son, 2016] (Figure 13). Estimation of groundwa-ter recharge due to rainwater infiltration is to be determined for the 2014-2018 rainy seasons from May to December.
The analysis of hourly rainfall data from May to December in 2014-2018 had provided rainfall events of different rainfall durations and depths (1 hr, 2 hrs, 3 hrs). The rainwater infiltration during each rainfall event equals infiltration determined by the moisture transfer modelling if the rainfall minus evaporation is greater than the infiltration amount determined by the moisture transfer modelling, otherwise equals rainfall minus evaporation. The average observed evaporation during May-December is 100 mm/month [Son, 2016] which is equivalent to 0.135 mm/hr. The accumulative rainwater infiltration for the rainfall durations of from 1 hr to 36 hrs for the five selected soil is given in Table 4. The estimated total rainwater infiltration during the 2014-2018 rainy seasons in terms of total infiltration depth for the five selected soil is given in Table 5.
The analysis gave varying rainwater infiltration during rainy seasons through the modeled soil profiles, from 0.280 m (silty clay) to 0.470 m (sandy silt), and from 33.3% to 55.2% of the rainfall depth during May-December. Low infiltration results in silty clay are due to low permeability and in sand are due to low suction, and high infiltration results in silt and sandy silt are caused by their higher moisture transfer. Since the evaporation is small since the rainfall duration is not long, the runoff is about from 66.7% (in the area of distribution of silty clay) to 24.8% (in the area of distribution of sandy silt). On average of total rainfall and soil properties, the average infiltration during May-December is 0.380 m, which is equivalent to 44.9% of the rainfall depth. The total volume of 289x106 m3 of rainwater may infiltrate into
Figure 11: Distribution of VWC versus depth
Figure 12: Accumulative rainwater infiltration
the Quaternary aquifer over 760 km2 of the coastal plain area of Ninh Thuan province.
Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
Figure 13: 1980-2014 average monthly rainfall
6 Discussions and concluding remarks
The estimated groundwater recharge via the rainwater infiltration through unsaturated soil was determined by the FE modelling using higherorder element functions with SWCC-based parameters in the combination of duration and rainfall depth of rainfall events in the period of 20142018.
The methodology proposed in this work seems to be generalized in terms that the unsaturated soil moisture transfer is carried out separately as an independent study. Then the rainfall data are analyzed to be grouped in continuous rainfall events of different durations. The rainfall depths of different-duration rainfall events and hourly evaporation in combination with temporal accumulative rainwater infiltration, determined by the unsaturated soil moisture transfer, shall provide the groundwater recharge during all the rainfall events within the study period. The unsaturated soil moisture modelling in this work had been carried out strictly by the given SWCC, expressing a relationship between soil moisture and suction as well as the unsaturated permeability and soil-moisture diffusivity.
As the groundwater recharge through all soil types and all depth profiles has been determined, then the groundwater recharge analysis over any
Table 4: Accumulative Rainwater Infiltration
Time (hr) Accumulative rainwater infiltration (m)
K = 0.013 m/d K = 0.058 m/d K = 0.102 m/d K = 0.146 m/d K = 0.200 m/d
1 0.00330 0.00596 0.00637 0.00576 0.00367
2 0.00395 0.00815 0.00964 0.00741 0.00419
3 0.00472 0.00994 0.01137 0.00859 0.00471
4 0.00537 0.01159 0.01260 0.00962 0.00539
5 0.00593 0.01257 0.01364 0.01049 0.00614
6 0.00643 0.01363 0.01457 0.01129 0.00690
7 0.00688 0.01441 0.01541 0.01205 0.00760
8 0.00731 0.01516 0.01619 0.01276 0.00820
9 0.00770 0.01584 0.01692 0.01344 0.00878
10 0.00807 0.01652 0.01762 0.01409 0.00934
11 0.00842 0.01715 0.01827 0.01470 0.00986
12 0.00875 0.01774 0.01890 0.01528 0.01038
13 0.00907 0.01831 0.01950 0.01584 0.01089
14 0.00937 0.01886 0.02008 0.01639 0.01142
15 0.00966 0.01938 0.02064 0.01692 0.01195
16 0.00994 0.01989 0.02118 0.01744 0.01245
17 0.01021 0.02038 0.02170 0.01794 0.01293
18 0.01048 0.02085 0.02221 0.01842 0.01339
19 0.01073 0.02131 0.02270 0.01889 0.01386
20 0.01098 0.02176 0.02318 0.01935 0.01431
21 0.01122 0.02220 0.02365 0.01980 0.01474
22 0.01145 0.02263 0.02411 0.02024 0.01513
23 0.01168 0.02304 0.02455 0.02067 0.01551
24 0.01191 0.02345 0.02499 0.02108 0.01588
25 0.01212 0.02385 0.02542 0.02150 0.01625
26 0.01246 0.02423 0.02584 0.02190 0.01663
27 0.01277 0.02462 0.02625 0.02229 0.01701
28 0.01305 0.02499 0.02665 0.02268 0.01740
29 0.01313 0.02536 0.02705 0.02307 0.01779
30 0.01340 0.02571 0.02744 0.02345 0.01819
31 0.01365 0.02607 0.02782 0.02382 0.01859
32 0.01376 0.02641 0.02820 0.02419 0.01899
33 0.01399 0.02675 0.02857 0.02455 0.01939
34 0.01423 0.02709 0.02894 0.02491 0.01979
35 0.01436 0.02742 0.02930 0.02526 0.02019
36 0.01458 0.02774 0.02966 0.02561 0.02058
area under concern may be carried out in an inexpensive work with the mapping of the distribution of soil types and unsaturated depths, and the analysis and interpretation of the rainfall data in the study area. An important issue is that the rainwater infiltration rate depends on the initial soil moisture content, which would require some more labor and time-consuming work.
Regarding the accuracy of the estimate of groundwater recharge from rainwater infiltration through unsaturated zone by solving Richards' equation which would be within a factor of 10 as [Kinzelbach et al., 2002] pointed out, the results of this work would show a much smaller factor. The accuracy of the groundwater recharge estimate would depend on the accuracy and appropriateness of input parameters (unsaturated soil
parameters), the accuracy of boundary and initial conditions, the selected element size and time step, and the way of the use of precipitation data in combination with the infiltration rates determined by the numerical modelling.
Finally, the results of the estimate of rainwater infiltration during rainy seasons to recharge the groundwater in the arid area of the coastal plain of Ninh Thuan province have shown that a large volume of rainwater would percolate into the aquifers. The groundwater resources are derived from the rainwater infiltration would be effectively managed for dealing with the water shortage in the area. The following two engineering measures would be worthwhile addressed. One measure is the construction of an underground dike to have twofold purposes of avoiding the groundwater dis-
Table 5: Total Rainwater Infiltration During May-December 2014-2018
Year Rainfall depth (m) Infiltration K = 0.013 m/d K = 0.058 m/d K = 0.102 m/d K = 0.146 m/d K = 0.200 m/d
2014 0.503 Ra (m) 0.193 0.292 0.304 0.272 0.204
Pb (%) 38.3 58.0 60.5 54.0 40.6
2015 0.805 R (m) 0.228 0.356 0.381 0.324 0.238
P (%) 28.3 44.2 47.3 40.3 29.6
2016 1.285 R (m) 0.380 0.641 0.673 0.583 0.409
P (%) 29.5 49.9 52.4 45.3 31.8
2017 0.847 R (m) 0.303 0.482 0.511 0.439 0.316
P (%) 35.7 56.9 60.2 51.8 37.3
2018 0.827 R (m) 0.287 0.437 0.458 0.408 0.308
P (%) 34.8 52.8 55.4 49.4 37.3
Average 0.853 R (m) 0.280 0.440 0.470 0.410 0.300
P (%) 33.3 52.4 55.2 48.2 35.3
R (m) 0.380
P (%) 44.9
charge into the sea and salinization from the sea [Hoang, 2005]. Another measure is the desaliniza-tion of existing salty groundwater by abstracting groundwater wells with an enhanced fresh surface water recharge to the groundwater [Hoang et al., 2018].
7 Acknowledgment:
This paper has been completed within the implementation of the grant research project KHCBTD.01/19-21: "Groundwater quality and quantity assessment for the coastal area of Ninh Thuan province for groundwater resources management under the context of drought and climate change" and with the financial support of the grant FEB RAS 18-1-008 (QTRU02.01/18-19). A number of registration AAAA-A18-118121090011-9, Grant of the FEB RAS "Abnormal gas-geochemical fields and petrochemical features as indicators of a hydrocarbon fluid, deep permeable zones, active tectonics and geoecology of Central Vietnam and the adjacent shelf" (No. 20-BAHT-010) and project VAST-FEB RAS code QTRU02.02/21-22. The research was carried out as part os State Program for basic scientific research AAAA-A19-119122090009-2 and 121021500055-0.
Conflict of Interest
None.
References
Bear, J., Computer-Mediated, Distance Learning Course on Modeling Groundwater Flow and Contaminant Transport, Topic D: Modeling Flow in the Unsaturated Zone, WebPage, 2000.
Beck, H. E., N. E. Zimmermann, T. R. McVicar, N. Vergopolan, A. Berg, and E. F. Wood, Present and Future Koppen-Geiger Climate Classification Maps at 1-km Resolution, Sci Data, 5, 18,021, doi:10.1038/sdata.2018.214, 2018.
Brooks, R. H., and A. T. Corey, Hydraulic Properties of the Porous Medium, Colorado State University (Fort Collins), Hydrology Paper, Nr., 3,1964.
Brutsaert, W., Probability Laws for Pore-Size Distributions, Soil Science, 101(2), 85-92,1966.
Carrera-Hernández, J. J., B. D. Smerdon, and C. A. Mendoza, Estimating Groundwater Recharge Through Unsaturated Flow Modelling: Sensitivity to Boundary Conditions and Vertical Discretization, Journal of Hydrology, 452-453, 90-101, 2012.
DARD, The Effect of the Pedologic Features of Agricultural Soils on the Quality of Ninh Thuan Grape, Ninh Thuan Department of Agriculture and Rural Development, 2017.
Fredlund, D. G., and A. Xing, Equations for the Soil-Water Characteristic Curve, Canadian Geotechnical Journal, 31(3), 521-532, 1994.
Fredlund, D. G., A. Xing, and S. Huang, Predicting the Permeability Function for Unsaturated Soils Using the Soil-Water Characteristic Curve, Canadian Geotechni-cal Journal, 31, 533-546, 1994.
Fredlund, D. G., S. Daichao, and Z. Jidong, Estimation of soil suction from the soil-water characteristic curve, Canadian Geotechnical Journal, 48,186-198, 2011.
Gardner, W. R., Water Movement Below the Root Zone, in Proc 8th Int Congr Soil Sci, RompresFilatelia, Bucharest, 1964.
van Genuchten, M. T., A Closed-form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils, Soil Science Society of America Journal, 44(5), 892-898, doi:10.2136/sssaj1980.03615995004400050002x, 1980.
Geo-Slope, Seepage Modeling with SEEP/W, an Engineering Methodology, 2015.
Gilding, B. H., Qualitative Mathematical Analysis of the Richards Equation, Transport in Porous Media, 6, 651666, 1991.
Green, W. H., and G. Ampt, Study on Soil Physics-1 the Flow of Air and Water Through Soils, The Journal of Agricultural Science, 4(1), 1-24, doi:10.1017/S0021859600001441, 1911.
Hai, T. T., Study and Assessment of Modern Geo-Tectonics of Vietnam Central Region and its Role in Natural Hazards for Predicting and Preventing Hazards in the Context of Climate Change, National Target Program KHCN-BDKH/11-15, 2015.
Hoang, N. V., On the Economics of Underground Dike to Protect Groundwater Salinization into Groundwa-ter Abstraction Facilities in the Coastal Areas, Vietnam Journal of Agriculture and Rural Development, pp. 54-55, 2005.
Hoang, N. V., T. N. Thanh, N. D. Roi, T. D. Huy, and T. T. Tung, The Potential of Desalination of Brackish Groundwater Aquifer Thanks to Salt-Intrusion Prevention River Gates in the Red River Delta, Vietnam, pp. 2747-2771, doi:10.1007/s10668-017-0014-x, 2018.
Horton, R. E., Analysis of Runoff Plat Experiments With Varying Infiltration Capacity, Eos, Transactions American Geophysical Union, 20(4), 693-711, doi:10.1029/TR020i004p00693, 1939.
Huyakorn, P. S., and G. F. Pinder, Computational Method in Subsurface Flow, Elsevier, Academic Press Inc., doi:10.1016/c2012-0-01564-5, 473 pages, 1983.
Jacobs, Water Efficiency Improvement in Drought Affected Provinces in the Central Coast and Central Highlands of Vietnam, Tech. rep., 2017.
Kinzelbach, W., W. Aeschbach, C. Alberich, I. B. Goni, U. Beyerle, P. Brunner, W.-H. Chiang, J. Rueedi, and K. Zoellmann, A Survey of Methods for Groundwater Recharge in Arid and Semi-Arid Regions. Early Warning and Assessment Report Series, UNEP/DEWA/RS.02-2, United Nations Environment Programme, Nairobi, Kenya, 2002.
Kunze, R. J., G. Uehara, and K. Graham, Factors Important in the Calculation of Hydraulic Conductivity, Soil Science Society of America Journal, 32(6), 760-765, doi:10.2136/sssaj1968.03615995003200060020x,
1968.
Leong, E. C., and H. Rahardjo, Permeability Functions for Unsaturated Soils, Journal of Geotechnical and Geoenvironmental Engineering, 123(12), 1118-1126, doi:10.1061/(ASCE)1090-0241(1997)123:12(1118), 1997.
Loan, T. T., Strategic Environmental Assessment for the Reviewing and Amending the Up to the Year of 2025 and 2030 - Vision Agricultural and Rural Planning of the Southern Central Region of Vietnam in the Context of Climate Change, National Institute of Agricultural Planning and Projection - MARD, 2018.
McKee, C. R., and A. C. Bumb, The Importance of Un-saturated Flow Parameters in Designing a Monitoring System for Hazardous Wastes and Environmental Emergencies, in Proceedings of the Hazardous Materials Control Research Institute, National Conference, pp. 50-58,1984.
McKee, C. R., and A. C. Bumb, Flow-Testing Coalbed Methane Production Wells in The Presence of Water and Gas, SPE formation Evaluation, 2(04), 599-608, 1987.
MONRE, QCVN09-MT:2015/BTNMT: National Technical Regulation on Groundwater Quality, Tech. rep., Ministry of Natural Resources and Environment, 2015.
Nie, W.-B., Y.-B. Li, L.-J. Fei, and X.-Y. Ma, Approximate Explicit Solution to the Green-Ampt Infiltration Model for Estimating Wetting Front Depth, Water, 9(8), doi:10.3390/w9080609, 2017.
Philip, J. R., The Theory of Infiltration. 4. Sorptivity and Algebraic Infiltration Equations, Soil Science, 84, 257264, 1957.
Philip, J. R., Theory of Infiltration, pp. 215-296, Elsevier, doi:10.1016/B978-1-4831-9936-8.50010-6,
1969.
Polubarinova-Kochina, P. Y., Theory of Groundwater Movement, Moscow, U.S.S.R.: Nauka. (In Russian), 1977.
Richards, L. A., Capillary Conduction of Liquids Through Porous Mediums, Physics, 1(5), 318-333, doi:10.1063/1.1745010, 1931.
Sam, L., and N. D. Vuong, The Selection to Research Formula of Drought Index and Applying to Calculate Droughty Frequency in Ninh Thuan Province, in Proceedings of 2008 Southern Vietnam Water Resources Academy Scientific and Technology Works, 2008.
Son, H. T., Study of Water Resources of the Desertification Area of Ninh Thuan Province in The Context of Climate Change and Proposal of Adaptive Measures, Ph.D. thesis, Vietnam Academy of Science and Technology, 2016.
Tarboton, D. G., Rainfall-Runoff Processes, Utah State VNCSCNDP&C, Report on Drought, Salinization and Pre-
University, 2003. vention Measures for UN Meeting, Vietnam Central
Steering Committee for Natural Disaster Prevention
Tu, N. T., et al., Hydrogeological Mapping at 1:50000 and Control 2016 Scale of Ninh Thuan and Binh Thuan Province, Tech. rep., Division for Water Resources Planning and Investigation for the Central Region of Vietnam, 2016.
UNEP, World Atlas of Desertification, 2nd ed., p. 182, Zienkiewics, O. C., and K. Morgan, Finite Elements and
London, UNEP, 1997. Approximation, 283-289 pp., John Willey & Sons.
