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 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.
tropical savannah climate, drought, moisture transfer, finite element (FE), higher-order element function
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 Köppen-Geiger climate classification sys- tem and a 1-kilometer resolution Köppen-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 condi- tions during the year. Ninh Thuan province is one of the provinces which are most affected by the drought. Severe drought in recent years (2012 to 2016) led to reduced water availability for irri- gated agriculture, especially for higher-value crops such as black pepper, coffee, dragon and grape- fruits. 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 suffer- ing from a shortage of domestic water [VNCSC- NDP&C, 2016]. In the context of the drought sit- uation, 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 ex- treme weather conditions had resulted in increas- ing water scarcity, and the farmers in these re- gions are heavily relying on groundwater for do- mestic use and irrigating the high-value crops [Ja- cobs, 2017]. which could contribute to the sustain- able 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 analy- sis 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 av- erage aridity index for the coastal plain area of Ninh Thuan province is 0.45 which is correspond- ing 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 desertifica- tion 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 man- agement are very important for the area. It is well known that in arid and semi-arid environ- ments, 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 semi- arid 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 meth- ods; 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 esti- mate methods have different accuracy and relia- bility, the authors have suggested that a successful estimation of groundwater recharge would have resulted from the utilization of several indepen- dent 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. De- spite a low accuracy rating of recharge estimate by unsaturated moisture transfer by Richards’ equa- tion as concluded by Kinzelbach et al. [2002], a sophisticated and well-based theoretical unsatu- rated moisture transfer with soil-water character- istic 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.
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] (Fig- ure 1) and consists of Quaternary deposits, which are divided into two aquifers, the Holocene aquifer (qh) and the Pleistocene aquifer (qp), and an undi- vided 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 thick- ness 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 Chloride- Sulfate 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 Sul- fate concentration is higher than 400 mg/l as spec- ified 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 re- maining 210 km2 are exposed to the ground sur- face 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 aver- age 7 meters. The groundwater level is from few decimeters above the ground surface to 6 meters below the ground surface, on average 2.45 me- ters 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 fresh- water distribution, in general, the aquifer water has an adequate quality for domestic use. How- ever, 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 consist- ing 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 pre- sented 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.
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 dif- ferent 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 pro- vided by Geo-Slope [2015] for silty clay, silt, silty sand and sand. An additional sandy silt SWCC between silt and silty sand. The hourly moni- toring rainfall data in the study area in the pe- riod 2014–2019 have been provided by the Hydro- Meteorological Information and Data Center un- der 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 el- ement functions for modelling the soil moisture transfer with the calculated soil SWCC-based parameter
Equation of soil moisture transfer in unsat- urated soil
The equation describing the moisture transfer in
Which is the Richards equation in terms of suction
pressure [Philip, 1969].
With water density equal to one and suction pressure ψ as a function of VWWC θw (3) is cal- culated in vertical direction z:
unsaturated soil with an assumption that air does
not move is as follows Polubarinova-Kochina [1977]
The component D(θw) = K (θw) dψ
is called
where: θw is volumetric water content (VWC), t is a time, υx; υy and υz are moisture transfer rates in x, y and z direction, respectively:
soil-moisture diffusivity [Philip, 1957] (as cited by [Gilding, 1991]) with unit L2T −1, and then (4) has the following form:
where: h = ψ/γ + z; h is a total head; ψ is soil suc- tion; γ is a water density, and K (θw) is hydraulic conductivity of unsaturated soil.
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(θw). Applying Green lemma
ψ was proposed by Gardner [1964] [Fredlund et al., 2011]
Where: θ is soil moisture content; θ is sat-
urated soil moisture content; ψ is soil suction; a
and n are fitting soil parameters associated with the soil-water characteristic curve (SWCC).
Since that time, there are numerous best- fit equations that have been proposed for the
SWCC of unsaturated soil [Brooks and Corey, 1964; Brutsaert, 1966; Fredlund and Xing, 1994;
The boundary integration Γ
+ ∂K (θw) = ∂θw .
∂z ∂t
is present only for the
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
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
+ ∂K (θw) = ∂θw .
SWCC, m is the fitting parameter related to the soil
∂z
Z !
∂N ∂W
∂t residual water content, e = 2.71828, ψ is soil suc- tion (kPa). Leong and Rahardjo [1997] concluded
that C(ψ) can be assumed to be unity without
K = Dz m l
Ω Z ∂z ∂z
dz;
affecting the initial portion of the SWCC, which helps to reduce the number of parameters in the
f = ∂K (θw) +
Γ
q W N d .
equation, i.e.:
∂z Γq
c l z
θw = θS .
Kθ + dθ = f . (7)
dt
ln e + (ψ/a)n m
Applying finite element method to (7) for the time domain, the following is [Zienkiewics and Morgan]:
The unsaturated permeability K (θw) of the soil can be determined through the relative coefficient of permeability kr (θ) [Fredlund et al., 1994]:
Zienkiewics and Morgan in very detail showed that the schemes with λ 0.5 are always unconditionally stable for any values of time step ∆t.
Soil moisture transfer parameters in ac- cordance with the soil-water characteristic curve
The first equation describing the relationship between soil moisture content θw and soil suction
where: K is the hydraulic conductivity of saturated soil; θ is the effective VWC defined as θ = θw θr ; θw is the VWC; θL is the lowest VWC on experi- mental SWCC; θr is residual VWC; θS is the sat- urated VWC; ζ is a dummy integration variable; Se = (θw θr ) / (θS θr ) is the effective saturation. In accordance to Kunze et al. [1968] [Fredlund et al., 1994] an improved pre- diction accuracy would have resulted if the complete SWCC is used, namely with the availability of residual VWC θr where
kr = 0. In this case, the integration in (8) should be carried out in the interval from θr to θS :
Selection of model domain: 1D vertical
z direction or 2D in horizontal x and vertical z
The soil-moisture diffusivity is:
D(θ ) = K (θ ) dψ .
w w dθw
Inputting initial moisture content and moisture dispersion coefficient at step n = 0
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 (relation- ships between the soil suction, relative permeabil- ity, 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 mod- elling is presented in Figure 4 along with the re- quired parameters derived from SWCC function.
Selection of element sizes ∆xi 10Dx (θ),
≤
≤
∆yi 10Dy (θ), 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 conte,nts at time step n:
K (θ) = R θ θ−ζ dζ R θS θS −ζ dζ;
Site condition and estimate of groundwater recharge from precip- itation
Soil physical parameters
Since the groundwater recharge from rainfall oc- curs 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 al- luvial deposits and has a thickness of more than one meter. This topsoil layer serves as agricul- tural soil, which has very specific characteristics,
r θr ψ2(ζ) θr ψ2(ζ)
θ = θw − θr ; K (θw) = K × kr (θ);
dθw
D(θw) = K (θw) dψ
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.
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 pump- ing 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 con- ductivity 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
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 Pleis- tocene 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
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 (Fig- ure 5). Therefore, the moisture transfer analysis is carried out for five cases of typical hydraulic con- ductivity values: minimal (0.13 m/d), average mi- nus 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 infor-
Figure 6: SWCCs of the modeled soils
the smaller element size the less error; 2) round- off error due to truncation in computer calculation process; and 3) approximation of the exact mathe- matical model by FE method [Zienkiewics and Mor- gan]. 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:
mation on vertical hydraulic conductivity, a refer- ence ratio between the horizontal (K ) and vertical
" #
∂ D(θw) ∂θw +
∂K (θw) = 0.
hydraulic conductivity (K
h
v) of 10 is used to get
Kv.
∂z ∂z ∂z
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 parame- ters of the SWCCs and soil physical properties are presented in Table 1. The SWCCs are presented in Figure 6 and the relative coefficient of permeabil- ity determined by Eq. (9) is presented in Figure 7.
Selection of element size in FE modelling
Finite element solutions have three types of er- ror: 1) discretization error due to the element size:
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
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 deter- mined. 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 er- ror estimator for a quadratic element [Zienkiewics and Morgan]:
(he )2 Z he
where the partial differential equation used by Zienkiewics and Morgan is the same as our mois- ture transfer equation where D(θw) = 1. In our work, even the time Galerkin scheme is used as recommended by Huyakorn and Pinder [1983] with λ = 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 ele- ment size used by other authors, which is from
0.01 m to 0.25 m, where the grid spacing or el-
Ω
Ω
kEk2 e = 12
R2 e dx.
0
ement size is constant over the model domain length of from 1 m to 10 m [Carrera-Hernández
where: Ωe is element domain; RΩe is an error within element e; he is element size.
Let us consider a hypothetical case with the ratio between permeability gradient and moisture coef- ficient value equal to 0.01, i.e., the following equa- tion:
∂2θw = 0
et al., 2012]. Carrera-Hernández et al. [2012] used different elements sizes, both constant and vari- able, 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 spa- tial discretization in unsaturated flow modelling. The authors concluded that the variable spatial
∂z2
.1.
discretization from 0.001 m with increment by a
If a quadratic element function is used with ele- ment size he, then the error estimate is:
(he )2 Z he
factor of 1.1 to a maximum value of 0.10 m pro- vided the fastest simulation times and the better accuracy of the numerical modelling.
Initial profile of soil volumetric water content
Initial VWC of the modeled soil profiles in the
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 or- der is O(∆t2)), and for the forward and backward difference schemes, the time step needs to be:
forward difference scheme : ∆t < (he)2/[6D (θw)]; backward difference scheme : ∆t < (he)2/[2D (θw)],
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 sea- son, i.e., the profiles of soil VWC need to be de- termined. 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 sur- face. 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 de- scribed by Bear [2000] in Figure 9.
Table 3: Error Estimate for Different Element Sizes
Figure 8: Initial field water content and suction profiles for the modeled soils
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 con- tinuously for 1 hr, 2 hrs, 3 hrs, etc. That is, the soil moisture content at the ground surface is cor- responding to the saturated water content.
The lower model domain is also set up as the first-type boundary condition since it directly con- tacts the water table. This is also in accordance with many other authors who assigned the bound- ary as either fixed water table, free drainage (unit gradient), or head-dependent [Carrera-Hernández et al., 2012]. Carrera-Hernández et al. [2012] bas- ing 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.
Selection of modelling depth of the soil
– −
Suppose that at the initial time t = 0 the ini- tial soil VWC profile is corresponding to the hy- drostatic 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 t1, 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
Simulation results
The simulation results for the modeled soils, wa- ter 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 pre- sented in Figure 12 and Table 4
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 funda- mentally 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]:
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
where: t is the time from the infiltration begin- ning; ∆θ is the difference between saturated VWC and initial VWC; h0 is the ponding depth; ψm is the suction head at the wetting front.
The monitoring rainfall data in the analysis pe- riod (2014–2019) have shown that the longest rain- fall 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 min- imal 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
during each rainfall event equals infiltration de- termined 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 mi- nus evaporation. The average observed evapo- ration 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 rain- fall 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 in- filtration results in silt and sandy silt are caused by their higher moisture transfer. Since the evap- oration 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 equiv- alent to 44.9% of the rainfall depth. The total vol- ume of 289×106 m3 of rainwater may infiltrate into
Figure 11: Distribution of VWC versus depth
6 Discussions and concluding re- marks
Figure 12: Accumulative rainwater infiltration
the Quaternary aquifer over 760 km2 of the coastal plain area of Ninh Thuan province.
Figure 13: 1980–2014 average monthly rainfall
The estimated groundwater recharge via the rainwater infiltration through unsaturated soil was determined by the FE modelling using higher- order element functions with SWCC-based param- eters in the combination of duration and rainfall depth of rainfall events in the period of 2014– 2018.
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 accu- mulative 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 car- ried 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 area under concern may be carried out in an inex- pensive work with the mapping of the distribution of soil types and unsaturated depths, and the anal- ysis and interpretation of the rainfall data in the study area. An important issue is that the rain- water 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 esti- mate would depend on the accuracy and appro- priateness 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 com- bination with the infiltration rates determined by the numerical modelling.
Finally, the results of the estimate of rainwa- ter 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
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 man- agement 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 num-
ber 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, ac- tive 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.
1. Bear, J., Computer-Mediated, Distance Learning Course on Modeling Groundwater Flow and Contaminant Transport, Topic D: Modeling Flow in the Unsaturated Zone, WebPage, 2000.
2. Beck, H. E., N. E. Zimmermann, T. R. McVicar, N. Vergopolan, A. Berg, and E. F. Wood, Present and Future Köppen-Geiger Climate Classification Maps at 1-km Resolution, Sci Data, 5, 18,021, doi:10.1038/sdata.2018.214, 2018.
3. Brooks, R. H., and A. T. Corey, Hydraulic Properties of the Porous Medium, Colorado State University (Fort Collins), Hydrology Paper, Nr., 3, 1964.
4. Brutsaert, W., Probability Laws for Pore-Size Distributions, Soil Science, 101(2), 85–92, 1966.
5. 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.
6. 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.
7. Fredlund, D. G., and A. Xing, Equations for the SoilWater Characteristic Curve, Canadian Geotechnical Journal, 31(3), 521–532, 1994.
8. Fredlund, D. G., A. Xing, and S. Huang, Predicting the Permeability Function for Unsaturated Soils Using the Soil-Water Characteristic Curve, Canadian Geotechnical Journal, 31, 533–546, 1994.
9. 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.
10. Gardner, W. R., Water Movement Below the Root Zone, in Proc 8th Int Congr Soil Sci, RompresFilatelia, Bucharest, 1964.
11. 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.
12. Geo-Slope, Seepage Modeling with SEEP/W, an Engineering Methodology, 2015.
13. Gilding, B. H., Qualitative Mathematical Analysis of the Richards Equation, Transport in Porous Media, 6, 651– 666, 1991.
14. Green, W. H., and G. Ampt, Study on Soil Physics-
15. 1 the Flow of Air and Water Through Soils, The Journal of Agricultural Science, 4(1), 1–24, doi:10.1017/S0021859600001441, 1911.
16. 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.
17. Hoang, N. V., On the Economics of Underground Dike to Protect Groundwater Salinization into Groundwater Abstraction Facilities in the Coastal Areas, Vietnam Journal of Agriculture and Rural Development, pp. 54–55, 2005.
18. 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-0014x, 2018.
19. 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.
20. 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.
21. Jacobs, Water Efficiency Improvement in Drought Affected Provinces in the Central Coast and Central Highlands of Vietnam, Tech. rep., 2017.
22. 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.022, United Nations Environment Programme, Nairobi, Kenya, 2002.
23. 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.
24. 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.
25. 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.
26. McKee, C. R., and A. C. Bumb, The Importance of Unsaturated 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.
27. 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.
28. MONRE, QCVN09-MT:2015/BTNMT: National Technical Regulation on Groundwater Quality, Tech. rep., Ministry of Natural Resources and Environment, 2015.
29. 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.
30. Philip, J. R., The Theory of Infiltration. 4. Sorptivity and Algebraic Infiltration Equations, Soil Science, 84, 257– 264, 1957.
31. Philip, J. R., Theory of Infiltration, pp. 215–296, Elsevier, doi:10.1016/B978-1-4831-9936-8.50010-6, 1969.
32. Polubarinova-Kochina, P. Y., Theory of Groundwater Movement, Moscow, U.S.S.R.: Nauka. (In Russian), 1977.
33. Richards, L. A., Capillary Conduction of Liquids Through Porous Mediums, Physics, 1(5), 318–333,
34. doi:10.1063/1.1745010, 1931.
35. 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.
36. 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.
37. Tarboton, D. G., Rainfall-Runoff Processes, Utah State University, 2003.
38. Tu, N. T., et al., Hydrogeological Mapping at 1:50000 Scale of Ninh Thuan and Binh Thuan Province, Tech. rep., Division for Water Resources Planning and Investigation for the Central Region of Vietnam, 2016.
39. UNEP, World Atlas of Desertification, 2nd ed., p. 182, London, UNEP, 1997.
40. VNCSCNDP&C, Report on Drought, Salinization and Prevention Measures for UN Meeting, Vietnam Central Steering Committee for Natural Disaster Prevention and Control, 2016.
41. Zienkiewics, O. C., and K. Morgan, Finite Elements and Approximation, 283–289 pp., John Willey & Sons.
