CLASSIC
Canadian Land Surface Scheme including Biogeochemical Cycles
waterInfiltrateUnsat.f90 File Reference

Evaluates infiltration of water into soil under unsaturated conditions. More...

Functions/Subroutines

subroutine waterinfiltrateunsat (WMOVE, TMOVE, LZF, NINF, ZF, TRMDR, R, TR, PSIF, GRKINF, THLINF, THLIQX, TBARWX, DELZX, ZBOTX, DZF, TIMPND, WADJ, WADD, IFILL, IFIND, IG, IGP1, IGP2, ILG, IL1, IL2, JL, N)
 

Detailed Description

Evaluates infiltration of water into soil under unsaturated conditions.

Author
D. Verseghy, M. Lazare

Function/Subroutine Documentation

◆ waterinfiltrateunsat()

subroutine waterinfiltrateunsat ( real, dimension (ilg,igp2), intent(inout)  WMOVE,
real, dimension (ilg,igp2), intent(out)  TMOVE,
integer, dimension (ilg), intent(inout)  LZF,
integer, dimension(ilg), intent(out)  NINF,
real, dimension (ilg), intent(inout)  ZF,
real, dimension (ilg), intent(inout)  TRMDR,
real, dimension (ilg), intent(in)  R,
real, dimension (ilg), intent(in)  TR,
real, dimension (ilg,igp1), intent(in)  PSIF,
real, dimension(ilg,igp1), intent(in)  GRKINF,
real, dimension(ilg,igp1), intent(in)  THLINF,
real, dimension(ilg,igp1), intent(in)  THLIQX,
real, dimension(ilg,igp1), intent(in)  TBARWX,
real, dimension (ilg,igp1), intent(in)  DELZX,
real, dimension (ilg,igp1), intent(in)  ZBOTX,
real, dimension (ilg), intent(inout)  DZF,
real, dimension(ilg), intent(inout)  TIMPND,
real, dimension(ilg), intent(inout)  WADJ,
real, dimension(ilg), intent(inout)  WADD,
integer, dimension (ilg), intent(in)  IFILL,
integer, dimension(ilg), intent(inout)  IFIND,
integer, intent(in)  IG,
integer, intent(in)  IGP1,
integer, intent(in)  IGP2,
integer, intent(in)  ILG,
integer, intent(in)  IL1,
integer, intent(in)  IL2,
integer, intent(in)  JL,
integer, intent(in)  N 
)
Parameters
[in,out]wmoveWater movement matrix \([m^3 m^{-2}]\)
[out]tmoveTemperature matrix associated with ground water movement [C}
[in,out]zfDepth of the wetting front [m]
[in,out]trmdrRemainder of time step after unsaturated infiltration ceases [s]
[in,out]lzfIndex of soil layer in which wetting front is located
[out]ninfNumber of levels involved in water movement
[in]rRainfall rate at ground surface \([m s^{-1}]\)
[in]trTemperature of rainfall [C]
[in]psifSoil water suction across the wetting front [m]
[in]grkinfHydraulic conductivity of soil behind the wetting front \([m s^{-1}]\)
[in]thlinfVolumetric liquid water content behind the wetting front \([m^3 m^{-3}]\)
[in]thliqxVolumetric liquid water content of soil layer \([m^3 m^{-3}]\)
[in]tbarwxTemperature of water in soil layer [C]
[in]delzxPermeable depth of soil layer \([m] (\Delta z_{g,w})\)
[in]zbotxDepth of bottom of soil layer [m]
[in]ifillFlag indicating whether unsaturated infiltration is occurring

The infiltration rate \(F_{inf}\) under conditions of a constant water supply can be expressed, e.g. in Mein and Larson (1973), as

\(F_{inf} = \)K_{inf} \( [(\Psi_f + z_f)/ z_f ]\)

where \(K_{inf}\) is the hydraulic conductivity of the soil behind the wetting front, \(\Psi_f\) is the soil moisture suction across the wetting front, and \(z_f\) is the depth of the wetting front. It can be seen that \(F_{inf}\) decreases with increasing \(z_f\) to an asymptotic value of \(K_{inf}\). Thus, if the rainfall rate r is less than \(K_{inf}\), the actual infiltration rate is limited by r, i.e. \(F_{inf} = r\). Otherwise, \(F_{inf}\) will be equal to R until the right-hand side of the above equation becomes less than r, after which point the above equation applies and ponding of excess water begins on the surface. The depth of the wetting front at this time \(t_p\) can be calculated by setting \(F_{inf}\) equal to r in the above equation and solving for \(z_f\). This results in:

\(z_f = \Psi_f /[r/K_{inf} – 1]\)

The amount of water added to the soil up to the time of ponding is \(t_p r\), or \(z_f (\theta_{inf} – \theta_l)\), where \(\theta_l\) and \(\theta_{inf}\) are respectively the liquid water content of the soil before and after the wetting front has passed. Setting these two equal and solving for \(t_p\) results in

\(t_p = z_f (\theta_{inf} – \theta_l)/r\)

In the 100 loop, a check is done for each successive soil layer to compare the infiltration rate in the layer with the rainfall rate. If \(K_{inf} < R\), a test calculation is performed to determine where the depth of the wetting front would theoretically occur at the ponding time \(t_p\). If the calculated value of \(z_f\) is less than the depth of the top of the soil layer, \(z_f\) is set to the depth of the top of the layer; if \(z_f\) falls within the soil layer, that value of \(z_f\) is accepted. In both cases, the index LZF is set to the index of the layer, and the flag IFIND, indicating that \(z_f\), has been successfully located, is set to 1. If the infiltration rate in the soil layer is greater than the rainfall rate, \(z_f\) is provisionally set to the bottom of the current layer, and LZF to the index of the next layer. IFIND remains zero. If the infiltration rate in the layer is vanishingly small, \(z_f\) is set to the depth of the top of the current layer, LZF to the index of the overlying layer, and IFIND to 1.

If LZF is greater than 1, some adjustment to the equation for \(t_p\) above is required to account for the fact that the values of \(\theta_{inf}\) and \(\theta_l\) in the layer containing the wetting front may differ from those in the overlying layers. The equation for \(t_p\) above can be rewritten as

\(t_p = [z_f [\theta_{inf}(z_f) – \theta_l(z_f)] + w_{adj}]/r\)

where \(w_{adj}\) is calculated as

\(LZF-1\)

\(w_{adj} = \Sigma[(\theta_{inf. i} – \theta_{l, i} ) – (\theta_{inf} (z_f) – \theta_l (z_f))] \Delta z_{g, w}\)

\(i=1\)

The adjusting volume WADJ is calculated in loop 200, and the time to ponding TIMPND in loop 250. If TIMPND is greater than the amount of time remaining in the current time step TRMDR, then unsaturated infiltration is deemed to be occurring over the entire time step. In this case, the amount of water infiltrating over the time step is assigned to the first level of the water movement matrix WMOVE and to the accounting variable WADD, and the temperature of the infiltrating water is assigned to the first level of the matrix TMOVE.

In loop 300 WADD is partitioned over the soil profile by comparing in turn the liquid water content of each soil layer with the calculated liquid water content behind the wetting front THLINF, and decrementing WADD layer by layer until a layer is reached in which the remainder of WADD is insufficient to raise the liquid water content to THLINF. If this condition is reached, LZF is set to the index of the soil layer; the depth of the wetting front DZF within the layer, obtained as WADD/(THLINF-THLIQX), is added to the depth of the bottom of the overlying layer to obtain ZF.

In loop 400, the water content in each soil layer J existing above ZF is assigned to the J+1 level of the water movement matrix WMOVE, and the respective water temperatures are assigned to TMOVE.

If TIMPND < TRMDR, the amount of water infiltrating between the start of the time step and TIMPND is again assigned to the first level of the water movement matrix WMOVE, and the temperature of the infiltrating water is assigned to the first level of the matrix TMOVE. The depth DZF of the wetting front within the layer containing it is calculated by subtracting the depth of the bottom of the overlying layer from ZF.

In loop 500, the water content in each soil layer J existing above ZF is assigned to the J+1 level of the water movement matrix WMOVE, and the respective water temperatures are assigned to TMOVE.

Finally, the time remaining in the current time step after the period of unsaturated infiltration is recalculated, and the counter NINF is set to LZF+1.