The rainfall/snowmelt and interception processes result in a net amount of rainfall reaching the soil surface. The infiltration process results in water at the surface, because the rainfall intensity is higher than the infiltration rate or, if the soil has a closed lower boundary, when the storage amount is exceeded.
There are two categories of infiltration processes available:
- A simple 1 or 2 layer solution of the Darcy equations: following a Green & Ampt (Kutilek and Nielsen, 2004), Smith & Parlange model (Morgan et al., 1998) or simple subtraction of ksat. These options are coupled to a simple drainage function.
- A full multilayered soil water balance, using the SWATRE model (Bastiaansen et al., 1996), using the Richards equation over a number of nodes in the soil profile.
Regardless of the method the profile can be open or closed at the lower boundary. A soil has a given depth with a storage determined by the difference between the porosity and the initial moisture content. The user can select to make the soil profile impermeable at the bottom. In other words the soildepth map is not automatically limiting, only when the infiltration option “impermeable lower soil boundary” is selected. If not, there will be percolation from the bottom of the profile (see below).
Simple 1 and 2 layer infiltration
The simple infiltration model assumes a wetting front parallel to the soil surface, and a soil that consists of one or two homogeneous soil layers. The driving forces for infiltration are gravity and the matrix suction at the wetting front. Above the wetting front the soil is assumed to reach saturation immediately, below the wetting front the initial moisture content prevails. Despite these simplifications, the infiltration will behave well in most situations. For each timestep a potential infiltration rate (fpot) is calculated ad compared to the rainfall intensity. If fpot is larger than the rainfall rate the actual infiltration rate (fact) equals the rainfall rate, else it equals the fpot. Because if e small timestep in openLISEM (generally in the order of 10 sec), it is not necessary to have an iterative procedure which estimates the time to ponding. A direct explicit solution works well.
Green & Ampt and Smith & Parlange
Both Green and Ampt (REF) and Smith and Parlange (REF) use a simplification of the Darcy equation for vertical water flow:
f = -Ks•(dh/dz + 1)
where: f is the infiltration rate (mm/h), Ks is the saturated hydraulic conductivity (mm/h), dh is the sum of the pressure by the water layer on the surface and the matrix suction at the wetting front (mm) and dz is the infiltration depth (mm). Following hydraulic conventions, the negative sign in these equations denotes a flow downward. Because dz is the cumulative infiltration F divided by the available (empty) pore space, it follows that:
dz = F/(θs-θi)
where: θs is the porosity (-) and θi is the initial moisture content below the wetting front (-).
Following a solution by Kutilek and Nielsen (2004), Green and Ampt simply combine these two equations:
f = fpot = -Ks•(dh•(θs-θi)/F + 1)
The negative sign indicates the direction of flow (downward). F will increase during the infiltration process until the bottom of the soil is reached.
The figure on the right shows the potential infiltration rate curve (fpot) for a 1:2 year event in Kampala, Uganda.
Following a solution form the EUROSEM model by (Morgan et al. 1998), Smith and Parlange calculate the fpot in a slightly different way:
C = exp(F/(dh•(θs-θi))
f = fpot = -Ks•C/(C-1)
The figure below shows the runoff generated from a 1 layer Green and Ampt solution in openLISEM.
Finally the option “subtract Ksat” in the openLISEM interface simply assumes:
f = fpot = -Ks
This method is only used for quick and dirty assessments.
Soils can be compacted and/or crusted. Two additional maps are needed for this. First a map specifying how much of the gridcell is crusted/compacted, called crustfrc.map or compfrc.map. Second a map specifybing the Ksat of crusted and/or compacted areas ksatcrst.map or ksatcomp.map. In the infiltration module an effective Ksat is calculated:
ksateff = Ksat1*(1-fcrust-fcomp) + Ksatcrust*fcrust + Ksatcomp*fcomp
Crusting or compaction can be switched on or off in the interface, it is not automatically done when the fraction and Ksat files are present. If a 2-layer infiltration method is chosen, the crusting/compaction is only applied to the first layer (Ksat1).
Swatre (Soil watr balance for terrestrial ecosystems) is created by Bastiaansen et al. (1996) at the Wageningen University (see here). Infiltration and soil water transport in soils are simulated by a solution of the well-known Richards equation, which combines the Darcy equation and the continuity equation:
K = the hydraulic conductivity (m/s);
h = the pressure (matric) potential (m);
θ = the volumetric water content (m3/m3);
z = the gravitational potential or height above a reference level (m);
t = time (s)
Using the soil water capacity C(h)=d(θ(h))/dh (the slope of the soil water retention curve θ(h), the unsaturated flow equation is derived:
This equation is solved for every soil layer using a so called explicit linearization using a Thomas algorithm for tri-diagonal matrices. For each timestep in openLISEM, the system of equations is solved using the smallest necessary timestep within the SWATRE module. This has proven stable enough for most situations compared to fully implicit solutions and is faster.
openLISEM uses tables with values for θ, h(θ) and k(h) instead of for instance Mualem/Van Genuchten or Brooks-Corey equations. In this way, users can apply their own measured data directly. Of course these tales can also be created with water retention and hydraulic conductivity equations.
The database to run SWATRE requires a relatively elaborate setup. It enbales a full 3D definition of the soil, using a classified map with soil units, called profile.map. The user has to define this classified map him/herself. Depending on what influences infiltration (soil texture and/or soil structure), some areas are better defined with soil units, others with land use types or a combination.
The numbers in the unit map are coupled to a text file profile.inp describing the soil profiles. This file first has the number of nodes and their depth, then for each unique profile number a description of the soil layers. Each layer has a name (such as “crust1” or “A_horizon”). These are names of tables (crust.txt or A_horizon.txt) that contains the values for the hydraulic conductivity curves and the water retention curves (θ, h(θ) and K(h) values). The figure below explains this setup.