Moisture Retention Curves in Aquifer Evaporation

INFLUENCE OF MOISTURE RETENTION CURVES ON THE ESTIMATION OF EVAPORATION OF SHALLOW AQUIFERS

[hr height=”20″]

The objective of this article is to determine the influence of unsaturated soil characteristic curves in determining evaporation in a column of bare soil. The VS2DTI code (Lappala et al., 1987) is used to determine evaporation as a function of groundwater depth and a sensitivity analysis of the different parameters of the characteristic curves is carried out. Additionally, the effect of using the curve h (θ) obtained in the laboratory is compared with those proposed from the granulometric curve.

[hr height=”20″ style=”dots” line=”default” themecolor=”1″]

THEORETICAL FRAMEWORK

Flow equation in the unsaturated zone

The equation that governs the flow in the unsaturated zone (Richards, 1931), is obtained by combining Darcy’s law with the continuity equation, which for a vertical flow in one dimension is written as:

[image src=”” size=”” width=”” height=”” align=”center” stretch=”0″ border=”0″ margin_top=”” margin_bottom=”” link_image=”” link=”” target=”” hover=”” alt=”” caption=”” greyscale=”” animate=””]

Where C (h) is the specific moisture capacity (L-1); K is the hydraulic conductivity (LT-1); z is the vertical coordinate (L), positive downward; h is the pressure head (L), a function of the volumetric water content θ (L3L-3). To obtain the solution of the vertical flow equation in one dimension, it is required to specify the appropriate initial and boundary conditions. The initial condition is given by the pressure head h throughout the profile and the boundary conditions that can be applied are the potential condition or flow condition.

Numerical solution of the Richards equation

There are several codes that allow solving the flow equation in the unsaturated zone. The spatial and temporal discretization of the domain is carried out through finite difference or finite element schemes. By means of finite differences it is possible to obtain adequate solutions for most of the problems, and it presents limitations when the geometry of the domain is complex and when the properties of the elements do not have a homogeneous behavior in all its directions. The finite element method allows adjusting to geometries of all kinds, for which it requires implementing appropriate discretization methods.

The VS2DTI (Variably Saturated 2 Dimensional Transport Interface) code solves the Richards equation for the flow in the unsaturated zone (Healy, 1990), using a finite difference scheme that allows to analyze cases of one or two dimensions with Cartesian or radial coordinates. . The suction and hydraulic conductivity curves are represented by the relationships proposed by van Genuchten (1980), Brooks and Corey (1964), Haverkamp et al. (1977) or experimental points. The code allows to establish boundary boundary conditions on hydraulic head (H), pressure (h), flow, infiltration, evaporation, transpiration of plants and filtration surfaces.

The evaporation condition can be imposed in two stages. The first stage occurs when the soil near the surface is moist and the water leaves the system at a rate equivalent to the evaporative demand of the atmosphere, known as the potential rate. This rate is maintained as long as the soil is able to satisfy this demand. As the soil dries up, the second stage begins in which the evaporative rate falls rapidly due to the decrease in the value of unsaturated hydraulic conductivity (Tyler et al., 2006).

The evaporative flux (EVPa) is established due to the pressure gradient existing between the atmosphere (imaginary cell) and the first cell of the domain represented by cell 1, as presented in Figure 1. (Lappala et al., 1987; Liu et al., 2005) and is written as:

where EVPa is the actual evaporation at the surface [LT-1]; EVP (t) is the potential evaporation from the soil [LT-1]; Ks is the saturated hydraulic conductivity [LT-1]; Kr is the relative hydraulic conductivity (-); SRES is the surface resistance that corresponds to the reciprocal of the distance from the node to the surface [L-1]; HA is the pressure potential of the atmosphere [L]; and h is the pressure potential at the first node on the surface of the modeled volume [L].

Figure 1. Schematic for the evaporation edge condition

The SRES parameter is used to model the existence of a less permeable surface crust than the soil (e.g. saline crust). HA is calculated using the Kelvin equation (Lappala et al., 1987) as:

where R is the gas constant (8.31 g m2 s − 2 K − 1 mol − 1); T is the air temperature (K); Mw is the molecular weight of water (0.018 Kg mol-1); g is the acceleration of gravity (9.81 m s − 2); and Rh is the relative humidity of the atmosphere (-).

[hr height=”20″]

METHODOLOGY

To achieve the proposed objectives, the characteristic curves were estimated from laboratory tests and the granulometric characteristics of a fine sand sample. Evaporation was determined by the numerical code VS2DTI.

Estimation of Characteristic Curves

The hydrodynamic characteristic curves of the soil were determined by two methods: 1) from laboratory tests, and 2) from empirical models or pedotransference functions (PFT). In the laboratory, the suction curve h (θ) was determined by the hanging column (ASTM, 2003) and pressure cooker (ASTM, 2003) methods, which consisted of applying different suctions to a soil sample and recording the gravimetric moisture content. First, the hanging column method was performed, saturating the sample with water and applying suctions between 0 and 80 kpa, maintaining each suction for one day and then the sample was transferred to the pressure cooker and suctions between 80 and 1500 kpa were applied. . The samples were subsequently dried in an oven at 110 +/- 5 ° C. The saturated hydraulic conductivity Ks was measured with a constant head permeameter, and the unsaturated hydraulic conductivity K (h) was determined by adjusting the parameters of the equations of van Genuchten (1980), which are written as:

Where Se is the effective saturation [-]; θh, θr and θs are the volumetric, residual and saturated moisture content respectively [L3L-3]; Ks is the saturated hydraulic conductivity [LT-1]; m and n are constants [-]; and m = 1-1 / n.

The estimation of the characteristic curves by means of the PFT was carried out with the ROSETTA model (Schaap et al., 2001) that estimates the parameters of the fit proposed by van Genuchten (1980), using neural networks (Schaap and Leij, 1998) by means of 5 hierarchical models that depend on the input data (USDA soil classification, texture, soil density, field capacity or moisture content at 330 cm suction and wilting point or moisture content at suction pressure 15 bar). In this work the models called H2 and H3 were used. The H2 model uses the percentage of sand, silt and clay as input data, and the H3 model incorporates, in addition to the texture, the density of dry soil as a predictor. The granulometric characteristics of the soil were determined by sieving and sedimentation tests.

Spatial discretization and boundary conditions

For the simulation of the evaporative flow, the VS2DTI code was used with a discretization in radial coordinates for the case of a circular column of 15.3 cm in diameter and 150 cm in height, and a variable refinement with finer cells in the upper area where the produces the evaporative flow. A total of 230 cells were used. The top cell has a thickness of 10-4 cm and the following cells grow in geometric progression at a rate of 5% more than the thickness of the previous one.

As the upper border condition of the domain, an evaporation condition was used that allows setting a value for the evaporation potential (EVP) and for the pressure potential of the atmosphere (HA). It was indicated that there was no crust on the surface that limited the evaporative flow and therefore the SRES variable was taken as 2 / Dz, considering the value of Dz of the upper cell. For the lower edge of the domain, a constant load edge condition (H = cte) was set to define the depth of the water table, which corresponds to the case in which the soil suction pressure is zero (h = 0). For the side edges a waterproof edge condition was imposed.

[hr height=”20″]

RESULTS AND DISCUSSION

Estimation of the characteristic curves h (θ) and K (θ) The model was applied to estimate the evaporation from groundwater in a column of soil whose granulometry, density, moisture content and hydraulic conductivity are known. In Figure 2 the granulometric curve is presented, which corresponds to a typical curve of sand with fine material less than 10%, and according to the classification of the USDA (United States Department of Agriculture) the soil studied has 95% of sand, 2% silt, and 3% clay.

The experimentally determined suction curve points were adequate to describe its shape in the suction pressure range between 0 and -10 m.c.a. The hanging column method is effective in defining the first section of the curve (0 to 1 m.c.a.) as is the pressure cooker method for higher suctions. It is recommended that both methods be applied in an intermediate range of pressures to verify that the values measured by both methods agree.

Table 1 shows the van Genuchten parameters of equations (4) and (5) for h (θ) and K (h). The suction curve adjusted to the experimental data of the fine sand sample is capable of maintaining its pores with a moisture content close to saturation (q »qs) in a wide capillary strip close to 60 cm. For a water depth of less than 60 cm, the humidity profile remains saturated and the evaporative flow will be dominated by the capacity of the atmosphere to induce it. The soil will then be able, according to the VS2DTI code, to provide the flow demanded at all times within this depth range of the water table.

Figure 3 shows the suction h (q) and hydraulic conductivity K (q) curves adjusted with van Genuchten (1980). The curves h (θ) estimated from the granulometry show clear differences mainly in the area close to saturation (Figure 3a). The saturated moisture content (qs) obtained through the H2 and H3 models were higher than that measured in the laboratory. The curve measured in the laboratory shows a capillary band that extends to a suction close to 50 cm, while for the adjustments using models H2 and H3 this band does not exceed 30 cm. The residual moisture content (qr) obtained with the H2 and H3 models is very similar to that measured in the laboratory.

In Figure 3b it can be observed that the differences between the curve fitted from the Ks value measured in the laboratory and those estimated by the H2 and H3 models from the granulometry are important. The domain for which the curve K (q) is defined is different, because the upper end (qs) is different for each of the settings. These differences are greater near the saturation condition and decrease towards the residual moisture content (q »qr). The slope of the K (q) curve associated with each adjustment is also modified due to these differences already mentioned, but also due to the variation of the parameter n in the van Genuchten adjustment that affects the value of m in equation (5) .

The greatest differences are found in the parameters p, n and Ks. Models H2 and H3 overestimate the values of p and Ks, and underestimate the values of n. An overestimation of the parameter p produces an underestimation of the bubbling or air inlet pressure (hb = 1 / p) and therefore a decrease in the height of the capillary fringe, while an increase (decrease) in the parameter n produces an increase (decrease) in the height of the capillary fringe (Figure 4). The differences between the estimates of the suction and experimental hydraulic conductivity curves and the curves fitted from the granulometry are mainly explained by the methodology used in the empirical models. The neural network method used by Schaap and Leij (1998) strongly depends on the data used in the training process of said networks, mainly on the quantity and textural distribution of the data, and despite the fact that the training for this model was carried out with about 10,000 h (θ) points and 1500 K (h) points, the uncertainty lies in the use for soil types other than those used in the training of neural networks.

Estimation of daily evaporation as a function of the depth of the water table

Figure 5 shows the evaporative flow EV (mm / day) estimated by solving equation (1) with the VS2DTI code, where depths are observed at which the calculated evaporation corresponds to the potential evaporation, called decoupling depth, between 60 and 80 cm, with an abrupt decay phase and a negligible evaporation flow under 140 cm depth. No great differences are observed in the evaporation values obtained.

The influence of the parameters p and n on the evaporation flux that is established for different depths of the water table is directly reflected with the depth of decoupling. The greater the thickness of the capillary strip, the greater the depth to which the soil is capable of providing the flow demanded by the atmosphere. This is reflected in the evaporation curves obtained by means of the numerical model where the evaporation curve associated with the hydrodynamic characteristics adjusted from the experimental data presents a decoupling depth greater than that obtained for the curves associated with the adjustments using the models H2 and H3. The capillary fringe of the suction curve adjusted from the experimental data is of the order of 65 cm, while the thickness of the capillary fringe of the suction curves obtained from the granulometry using models H2 and H3 are of the order 30 cm.

Sensitivity analysis of VG parameters in evaporation flow

As a base case of evaporation flow, the one calculated with the curves h (q) and K (q) adjusted to the experimental measurements was used and a variation of 5% was applied on each parameter and the evaporative flow was obtained for different depths of the nappa: 80, 100 and 140 cm deep. Table 2 shows the percentage variation of the evaporative flow obtained according to each disturbed parameter, where negative values represent a percentage decrease in the evaporative flow with respect to the base situation considered.

The parameters qr and qs do not affect the behavior of the model on the evaporative flow calculated when the layer is over 100 cm deep, although it shows a variation of the order of 0.2% when it is at a depth of 140 cm. The sensitivity associated with the HA parameter does not seem to be relevant. In all cases, a variation of less than 0.1% was obtained, so that its influence is considered negligible in the face of variations of that magnitude (± 5%).

The saturated hydraulic conductivity has a similar sensitivity for the different depths of the water table considered. Likewise, a variation in the calculated evaporative flow is observed that is equivalent to that imposed on the parameter Ks, since when the value of Ks is increased or decreased there is an increase or decrease of the same percentage magnitude in the calculated evaporative flow., This It occurs when the actual evaporation is less than the potential evaporation since the evaporative flow is proportional to the saturated hydraulic conductivity and HA does not vary since the temperature is constant, as it appears in equation 2.

The model is very sensitive to changes in the coefficients p and n. The highest percentage variation for these parameters was recorded when the water table is 140 cm deep, with a flow increase of 193.11% for the parameter of p and 119.20% for the case of n with respect to the initial situation. As with the parameter p, when n decreases, the calculated evaporative flux increases and decreases as n increases. This shows that there is a negative relationship between the evaporative flux and the parameters p and n.

[hr height=”20″]

CONCLUSIONS AND RECOMMENDATIONS

Despite the large differences observed in the estimation of the characteristic curves h (q) and k (q) obtained from experimental tests and from the PFTs, these do not generate large differences in the values of the evaporation flux estimated in the modeling .

The greatest differences observed in the characteristic curves and in the estimated evaporation flux are produced by the parameters p, n and Ks. Models H2 and H3 overestimate the values of p and Ks, and underestimate the values of n.

The VS2DTI code is capable of predicting evaporative flux but has some limitations. The edge condition at the surface produces a decoupling depth that appears to be determined by the thickness of the capillary fringe, assuming an unsaturated flow that is not real; The model solves the equations only for water transport, for which it imposes an edge condition as a function of potential evaporation generating a decoupling distance that does not exist in reality and that is due to the movement of liquid water, water vapor water and heat are strongly interrelated. Therefore, simulations and analytical solutions that are based only on the Richards equation are not able to describe the entire profile.

INFLUENCE OF MOISTURE RETENTION CURVES ON THE ESTIMATION OF EVAPORATION OF SHALLOW AQUIFERS

THEORETICAL FRAMEWORK

Flow equation in the unsaturated zone

Numerical solution of the Richards equation

METHODOLOGY

Estimation of Characteristic Curves

Spatial discretization and boundary conditions

RESULTS AND DISCUSSION

Estimation of daily evaporation as a function of the depth of the water table

Sensitivity analysis of VG parameters in evaporation flow

CONCLUSIONS AND RECOMMENDATIONS

Dejar un comentario Cancelar respuesta