A modification of the mixed form of Richards equation and its application in vertically inhomogeneous soils

Recently, new soil data maps were developed, which include vertical soil properties like soil type. Implementing those into a multilayer Soil-Vegetation-Atmosphere-Transfer (SVAT) scheme, discontinuities in the water content occur at the interface between dissimilar soils. Therefore, care must be taken in solving the Richards equation for calculating vertical soil water fluxes. We solve a modified form of the mixed (soil water and soil matric potential based) Richards equation by subtracting the equilibrium state of soil matrix potentialψE from the hydraulic potential ψh. The sensitivity of the modified equation is tested under idealized conditions. The paper will show that the modified equation can handle with discontinuities in soil water content at the interface of layered soils.


Introduction
The Richards equation (RE) (Richards, 1931) is commonly used in many Land Surface Models (LSMs) to describe vertical soil water flow in unsaturated zones.It is a combination of mass conservation and Darcy's law.There are three different expressions for it, (1) by soil water content ("θ-based"), (2) by hydraulic potential ("ψ-based") and (3) by both, hydraulic potential and soil water content ("mixed form") as the dependent variables (e.g., Hillel, 1980;Bohne, 2005).The "ψ-based" RE is often used by soil scientists and has the advantage to be continuous in the dependent variable at soil interfaces, but the numerical solution shows errors in mass balance for any iteration method (Picard, Newton-Raphson, etc.) (Celia et al., 1990).The "θ-based" RE is favored by climate modelers because of its excellent mass balance, but has the disadvantages that it cannot handle water flows in saturated zones and that discontinuities in the dependent variable occur across soil layer boundaries (Hills et al., 1989).Using the "mixed form" of RE with a finite element method for numerical solution, Celia et al. (1990) showed, that mass balance problems were consistently low in a homogeneous soil and that discontinuity problems are better-natured than in the "θ-based" RE.
Correspondence to: F. Kalinka (frank.kalinka@iau.uni-frankfurt.de)Implementing a soil data map that include information about vertical soil properties like soil type, care must be taken in solving RE to model unsaturated soil water flow in layered soils, because discontinuities in the water content occur at the interface between dissimilar soils.Many efforts have been done to handle with this problem (Hills et al., 1989;Romano et al., 1998;Schaudt and Morrill, 2002;Matthews et al., 2004;Barontini et al., 2007).But to a greater or lesser extent, all methods and results show problems in mass conservation.Zeng and Decker (2008) applied a modified mixed form of RE for unsaturated homogeneous soils, where mass conservation is nearly given.We expected, that this combination will also lead to better results in modeling water flows in inhomogeneous soils and herefore, implemented this formulation into the LSM model TERRA-ML, which is included in the non-hydrostatic mesoscale weather forecast model COSMO (http://www.cosmo-model.org/)as well as in the regional climate model COSMO-CLM (http: //www.clm-community.eu/).
In the Theory section we give a short overview over TERRA-ML and discuss the applied modifications.In the Result section we show the application of the modified RE in a sand-loam-sand and a loam-sand-loam soil column under idealized conditions (no infiltration, no gravitational drainage and initialization with the hydrostatic equilibrium state).Finally the results are summarized and discussed within the Conclusion section.

Model description
TERRA-ML is the actually implemented LSM of the nonhydrostatic mesoscale weather forecast model COSMO and of the regional climate model COSMO-CLM.Only a brief description of the model can be given here but more details can be found in Schrodin and Heise (2001), Doms et al. (2005) and Heise et al. (2006).Combining conservation of mass with Darcy's law leads to the mixed form of RE, describing one-dimensional water flow in soils where is the hydraulic conductivity and S is a sink term (e.g.transpiration loss in the root zone).Taking the hydraulic diffusivity (e.g.Kabat et al., 1997) into account, leads to the θ-based RE for soil water fluxes used in TERRA-ML The hydraulic conductivity K and the hydraulic diffusivity D both depend on soil moisture θ and soil properties.Both functions are parameterized according to Rijtema (1969) using exponential functions.Soil matric potential ψ m for unfrozen soil is defined according to Clapp and Hornberger (1978) at the node depth z i of each layer i as where the air entry potential at saturation ψ sat is an exponential function of percentages of sand, and the pore size distribution index B is a linear function of percentages of clay (Cosby and Hornberger, 1984).Celia et al. (1990) showed, that best simulation results with respect to mass conservation are obtained using the mixed form of RE.Additionally, using an inhomogeneous soil type distribution, Eq. (3) will be discontinuous in θ at the interface between two different soil horizons, so that ∂θ/∂z is not solvable using small soil layer thicknesses.However, the hydrological potential ψ is continuous, that is the reason why we test Eq.(1) as possible alternative to Eq. (3).Furthermore Grasselt et al. (2008) showed that, using a linear function of K, rather than the operationally used exponential function used in TERRA, also leads to results closer to observations.That is why we have implemented an empirical linear function of hydraulic conductivity according to Clapp and Hornberger (1978).

Equilibrium state
For idealized case studies in numerical weather prediction as well as in climate simulations, it is common to do model-runs in equilibrium state as initialisation.In our case the hydrostatic equilibrium state of Eq. ( 1) is theoretically reached, if no fluxes occur between the soil layers (∂θ E /∂t = 0), that is as a special solution under assumption of a constant ψ h in saturated zones (C = ψ sat + z WTD ).Here z WTD is the water table depth (WTD) and the subscript E in θ E and ψ m,E denotes the equilibrium state.The layer averaged θE,i can then be derived by integrating over soil layers.

Modified Richards equation
Because the derivative of C with depth vanishes, we can rewrite Eq. ( 1) to The method of subtracting the equilibrium state has the theoretical advantage, that sharp bends in the hydraulic potential ψ h , arising at the interface between two dissimilar soils because of different soil type properties, will be minimized in the numerical solution.Furthermore, small discontinuities in ψ m , occurring at the boundary between layered soils as a result of the numerical solution of Eq. ( 4) in combination with the soil type dependent hydrological parameters given in Table 1, are also reduced.

Model setup
For all simulations we used the mixed form of RE (Eq. 1) first and the modified form (Eq. 6) afterwards.Here, the operationally used "θ-based" RE (Eq. 3) is not used anymore because of its disadvantage that it can not simulate soil water fluxes across soil layer boundaries, as discussed in Sect.2.1.In general, horizontal fluxes are not considered in TERRA-ML, we therefore simulated only one column.The following model setup has been used: -specific soil type properties B, ψ sat , % sand and % clay were used as given in Table 1, -time stepping ∆t = 60 s, The last four assumptions ensure a closed system, where no fluxes should be simulated in a perfect model.

Varying layer thickness
Applying the model setup described above with the old formulation of RE (Eq.1), using a WTD of 3m and a seven layer soil distribution, fluxes up to 0.05 mm per hour occur especially from the sand column into the loam column in both, a loam-sand-loam distribution and a sand-loam-sand distribution (grey lines in Fig. 1).These fluxes can not be removed by using smaller soil layer thicknesses to minimize ∂z and therefore to improve results of numerical solution.Figure 2 shows fluxes occurring while using homogeneous soil layer thicknesses of 10 cm.Fluxes from the sand layer into the loam layer do not get smaller, they rather enlarge up to 0.6 mm/h in a loam-sand-loam distribution and up to 0.3 mm/h in a sand-loam-sand column (grey lines).
Using the modified formulation of RE (Eq. 6) instead of Eq. ( 1), no fluxes occur so that mass is perfectly conserved, independent of soil layer thickness and soil type distribution (black lines in Figs. 1 and 2, pink line in Fig. 2).

Influence of water table depth
We have implemented different WTDs to show its influence on soil water fluxes in an inhomogeneous soil type distribution.If the WTD is below the model domain (3 m and 5 m), fluxes only occur at the interfaces between different soil horizons (Fig. 2).If the WTD is in the model domain (0.5, 1 and 2 m), layers with depth below or at the WTD become saturated, so that θ E = θ sat and ψ E = ψ sat .Figure 3 shows simulated fluxes for different WTDs.It is conspicuous, that below and at saturated layers, huge fluxes occur up to 230 mm h −1 .Again they especially occur at the interfaces between different soil horizons as a result of the influence of different soil properties and this shows that the mixed form of RE is not adequate in simulating water flows in saturated soils.Testing the sensitivity of the modified RE to the WTD, results again indicate perfect mass conservation.The black, red and the dark green lines in Fig. 3 show, that no fluxes occur, independent on WTD and soil type distribution.This demonstrates the efficiency of the modified RE.

Conclusions
This study shows that the "θ-based" and the "mixed form" of RE are not adequate in simulating water transports in nonhomogeneous unsaturated soils.Therefore, we tested a modified form of the "mixed" RE after Zeng and Decker (2008), where the equilibrium state is subtracted of the hydraulic potential.Results show almost perfect mass conservation, independent of soil layer thicknesses used for the numerical solution, of the WTD (if it is below or in the model domain) and of what soil type distribution is adopted.However, all simulations were done under idealized test cases (neglecting the sink term S , prescribing a water table depth, initializing soil moisture with its equilibrium state, no infiltration as upper boundary condition and no gravitational drainage).Further simulations under more realistic conditions still have to be done.

Figure 1 .
Figure 1.(left) Simulated mean water fluxes [mm h −1 ] in a loam-sand-loam soil column after an integration period of 30 days with a WTD of 3 m: the grey line shows water fluxes of the old formulation (Eq.1), the black line of the new formulation (Eq.6).(right) Same as (left) exept with a sand-loam-sand soil column.

Figure 2 .
Figure 2. (left) Simulated mean water fluxes [mm h −1] in a loam-sand-loam soil column after an integration period of 30 days with homogeneous soil layer thicknesses of 10 cm and with different WTDs at 3 m and 5 m in a loam-sand-loam soil column: the grey and red lines show fluxes in mm h −1 of the old formulation (Eq.1), the black and purple lines of the new formulation of RE (Eq.6).(right) Same as (left) except using a sand-loam-sand soil type distribution.

Figure 3 .
Figure 3. (left) Same as Fig. 2a, except using WTDs within the model domain at 0.5, 1 and 2 m (indicated by the blue lines).Results of the new formulation of RE (black, purple and dark green lines) overlap each other.(right) Same as (left) except using a sand-loam-sand soil type distribution.

Table 1 .
Soil type dependent hydrological parameters Soil-type % sand % clay ψ sat [mm] B