- Generated on Mon Feb 3 2020 19:58:04 for CLASSIC by
1.8.13
|
CLASSIC
Canadian Land Surface Scheme including Biogeochemical Cycles
|
Quantifies movement of liquid water between soil layers under non-infiltrating conditions, in response to gravity and tension forces. More...
Functions/Subroutines | |
| subroutine | waterflownoninfiltrate (IVEG, THLIQ, THICE, TBARW, FDT, TFDT, BASFLW, TBASFL, RUNOFF, TRUNOF, QFG, WLOST, FI, EVAP, R, ZPOND, DT, WEXCES, THLMAX, THTEST, THPOR, THLRET, THLMIN, BI, PSISAT, GRKSAT, THFC, DELZW, XDRAIN, ISAND, LZF, IGRN, IGRD, IGDR, IG, IGP1, IGP2, ILG, IL1, IL2, JL, N) |
Quantifies movement of liquid water between soil layers under non-infiltrating conditions, in response to gravity and tension forces.
| subroutine waterflownoninfiltrate | ( | integer, intent(in) | IVEG, |
| real, dimension (ilg,ig), intent(inout) | THLIQ, | ||
| real, dimension (ilg,ig), intent(inout) | THICE, | ||
| real, dimension (ilg,ig), intent(inout) | TBARW, | ||
| real, dimension (ilg,igp1), intent(inout) | FDT, | ||
| real, dimension (ilg,igp1), intent(inout) | TFDT, | ||
| real, dimension(ilg), intent(inout) | BASFLW, | ||
| real, dimension (ilg), intent(inout) | TBASFL, | ||
| real, dimension(ilg), intent(inout) | RUNOFF, | ||
| real, dimension (ilg), intent(inout) | TRUNOF, | ||
| real, dimension (ilg), intent(inout) | QFG, | ||
| real, dimension (ilg), intent(inout) | WLOST, | ||
| real, dimension (ilg), intent(in) | FI, | ||
| real, dimension (ilg), intent(in) | EVAP, | ||
| real, dimension (ilg), intent(in) | R, | ||
| real, dimension (ilg), intent(in) | ZPOND, | ||
| real, dimension (ilg), intent(in) | DT, | ||
| real, dimension(ilg), intent(inout) | WEXCES, | ||
| real, dimension(ilg,ig), intent(inout) | THLMAX, | ||
| real, dimension(ilg,ig), intent(inout) | THTEST, | ||
| real, dimension (ilg,ig), intent(in) | THPOR, | ||
| real, dimension(ilg,ig), intent(in) | THLRET, | ||
| real, dimension(ilg,ig), intent(in) | THLMIN, | ||
| real, dimension (ilg,ig), intent(in) | BI, | ||
| real, dimension(ilg,ig), intent(in) | PSISAT, | ||
| real, dimension(ilg,ig), intent(in) | GRKSAT, | ||
| real, dimension (ilg,ig), intent(in) | THFC, | ||
| real, dimension (ilg,ig), intent(in) | DELZW, | ||
| real, dimension(ilg), intent(in) | XDRAIN, | ||
| integer, dimension (ilg,ig), intent(in) | ISAND, | ||
| integer, dimension (ilg), intent(in) | LZF, | ||
| integer, dimension (ilg), intent(in) | IGRN, | ||
| integer, dimension (ilg), intent(inout) | IGRD, | ||
| integer, dimension (ilg), intent(in) | IGDR, | ||
| 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 | ||
| ) |
| [in] | iveg | Subarea type flag |
| [in,out] | thliq | Volumetric liquid water content of soil layer \([m^3 m^{-3}] (\theta_l)\) |
| [in,out] | thice | Volumetric frozen water content of soil layer \([m^3 m^{-3}] (\theta_i)\) |
| [in,out] | tbarw | Temperature of water in soil layer [C] |
| [in,out] | fdt | Water flow at soil layer interfaces during current time step [m] |
| [in,out] | tfdt | Temperature of water flowing between soil layers [C] |
| [in,out] | basflw | Base flow from bottom of soil column [m] |
| [in,out] | tbasfl | Temperature of base flow from bottom of soil column [K] |
| [in,out] | runoff | Total runoff from soil column [m] |
| [in,out] | trunof | Temperature of total runoff from soil column [K] |
| [in,out] | qfg | Evaporation from soil surface (diagnostic) \([kg m^{-2} s^{-1}]\) |
| [in,out] | wlost | Residual amount of water that cannot be supplied by surface stores \([kg m^{-2}]\) |
| [in] | fi | Fractional coverage of subarea in question on modelled area [ ] \((X_i)\) |
| [in] | evap | Evaporation rate from ground surface \([m s^{-1}]\) |
| [in] | r | Rainfall rate at ground surface \([m s^{-1}]\) |
| [in] | zpond | Depth of ponded water on soil surface [m] |
| [in] | dt | Time period over which water movement takes place [s] |
| [in] | igrn | Flag to indicate whether calculations in subroutine waterFlowInfiltrate are done |
| [in] | lzf | Index of soil layer in which wetting front is located |
| [in,out] | igrd | Flag to indicate whether calculations in this subroutine are to be done |
| [in] | igdr | Index of soil layer in which bedrock is encountered |
| [in] | thpor | Pore volume in soil layer \([m^3 m^{-3}] (\theta_p)\) |
| [in] | thlret | Liquid water retention capacity for organic soil \([m^3 m^{-3} ] (\theta_{l,ret})\) |
| [in] | thlmin | Residual soil liquid water content remaining after freezing or evaporation \([m^3 m^{-3}] (\theta_{l,min})\) |
| [in] | bi | Clapp and Hornberger empirical "b" parameter [ ] (b) |
| [in] | psisat | Soil moisture suction at saturation \([m] (\Psi_{sat})\) |
| [in] | grksat | Hydraulic conductivity of soil at saturation \([m s^{-1}] (K_{sat})\) |
| [in] | thfc | Field capacity \([m^3 m^{-3}]\) |
| [in] | delzw | Permeable depth of soil layer \([m] (\Delta_{zg,w})\) |
| [in] | xdrain | Drainage index for water flow at bottom of soil profile [ ] |
| [in] | isand | Sand content flag |
In loop 50, the flag IGRD is first set to 1 for all grid cells where the calculations in this subroutine are to be performed. The necessary conditions are: that the surface being modelled is not a glacier or ice sheet (ISAND > -4); that the time period DT is greater than zero; that the infiltration calculations in subroutine waterFlowInfiltrate are not simultaneously being performed (IGRN = 0); and that the rainfall rate and the depth of ponded water are both vanishingly small. If any of these conditions is not met, IGRD is set to zero.
In loop 100, if the surface being modelled is not an ice sheet (ISAND = -4) and if the soil layer in question is not completely rock (ISAND = -3), the maximum possible water content of each soil layer, THLMAX, is calculated as the maximum of the available pore volume THPOR – THICE (where THPOR is the total pore volume and THICE is the ice content of the layer), the actual liquid water content THLIQ, and the minimum residual liquid water content THLMIN. The last two conditions are required because in the case of a saturated soil undergoing freezing, since water expands when frozen, the sum of the liquid and frozen volumetric water contents may be greater than the pore volume, and thus THLIQ or THLMIN may be greater than THPOR – THICE. An effective saturated hydraulic conductivity GRKSATF of the soil layer is also defined, applying an empirical correction for the presence of ice. This ice content factor, fice, is calculated from Zhao and Gray (1997) as:
\(f_{ice} = [1.0 – min((\theta_p - \theta_{l, min} )/\theta_p , \theta_i / \theta_p)]^2\)
The pore volume of the soil layer is corrected for the presence of ice by defining the effective porevolume THPORF as equivalent to THLMAX.
In loops 150 and 200, the theoretical amounts of water FDT flowing across the soil layer boundaries are calculated. FDT is evaluated as the product of the flow rate F(z) at the given depth z, and the period of time (DT) over which the flow is occurring. At z=0, the water flux is simply equal to the soil surface evaporation rate. Within the soil, the flow rate F at a given depth z is obtained using equation 21 from Verseghy (1991):
\(F(z) = K(z) [-b \Psi(z)/\theta_l(z) \times d\theta_l / dz + 1]\)
where K(z) is the hydraulic conductivity and psi(z) is the soil moisture suction at depth z, and b is an empirical parameter developed by Clapp and Hornberger (1978) [21]. K(z) and \(\Psi(z)\) are calculated following Clapp and Hornberger as
\(K(z) = K_{sat} (\theta_l/\theta_p)^{(2b + 3)}\)
\(\Psi(z) = \Psi_{sat} (\theta_l/\theta_p)^{(-b)}\)
where \(K_{sat}\) and \(\Psi_{sat}\) are the values of \(K\) and \(\Psi\) respectively at saturation.
At the bottom of the permeable soil depth, in layer IGDR, if the liquid water content of the soil is greater than the field capacity, the vertical flow out of the bottom of the soil profile is calculated using a relation derived from Soulis et al. (2010):
\(F(z_b) = K_{sat} \times min{1, (\theta_l /\theta_p)/[1 – 1/(2b + 3)]^{(2b + 3)}}\)
This flow rate is multiplied by a drainage parameter XDRAIN, which is set to 0 if the soil is underlain by an impermeable layer (as in a bog), and to 1 otherwise.
Between soil layers, values of \(K(z)\) and \(\Psi(z)\) must be determined. This requires the estimation of soil properties at the layer interfaces. For \(\theta_l\), \(\theta_p\) and \(d(\theta_l)/dz\), slightly different approaches are followed if the permeable soil layer thickness DELZW increases with depth (the normal case), or if it decreases between one soil layer and the next (indicating the presence of an impermeable barrier at some depth in the second layer). In the first case, the pore volume THPBND, liquid water content THLBND, and liquid water gradient DTHLDZ are simply calculated as arithmetic averages of the values above and below the interface. This has the effect of weighting the interface values towards the upper layer values, roughly emulating a surfaceward exponential decay curve of liquid water content. In the second case, in order to avoid a spurious weighting towards the lower layer, the gradients of liquid water content and pore volume between the upper and lower layers are calculated as linear relations, and then solved for the values at the interface.
The Clapp and Hornberger b parameter at the interface is calculated as a simple average of the values in the upper and lower soil layers. The saturated hydraulic conductivity \(K_{sat, bnd}\) is calculated as a harmonic mean, and the saturated soil moisture suction \(\Psi_{sat, bnd}\) as a geometric mean, of the top and bottom layer values:
\( K_{sat, bnd} = K_{sat, t} K_{sat, b} ( \Delta z_{g, w, , t} + \Delta z_{g, w, , b} ) / ( K_{sat, t} \Delta z_{g, w, b} + K_{sat, b} \Delta z_{g, w, t} ) \)
\( \Psi_{sat, bnd} = \Psi_{sat, t}^{\Delta zg, w, t/(\Delta zg, w, t + \Delta zg, w, b)} \Psi_{sat, b}^{\Delta zg, w, b/(\Delta zg, w, t + \Delta zg, w, b)} \)
Finally, K(z), \(\Psi(z)\) and F(z) at the interface are calculated from the equations given above. At the end of the loop, if set to zero behind the soil layer LZF containing the wetting front, and the flow rate at the bottom of LZF is constrained to be \(\geq\) 0.
In the next several loops, checks are carried out to ascertain whether the theoretically determined flow rates at the soil layer interfaces can be supported by the water contents in the layers. At the beginning of loop 250, the soil surface evaporation rate is addressed. A trial liquid water content THTEST is calculated for the first soil layer, reflecting the removal of the evaporated water. If THTEST is less than the minimum soil water content, F(0) is set to zero and \(\theta_l\) is set to \(\theta_{l, min}\). The excess surface flux that was not met by the removed liquid water is converted to a frozen water content, and an attempt is made to remove it from the frozen water in the layer. (The energy involved in converting the required frozen water mass to liquid water is obtained by decreasing the water temperature of the layer.) The frozen water content of the layer is adjusted to reflect the removal of the required amount. If the demand exceeds the supply of frozen water available, the frozen water content is set to zero, the diagnostic evaporative flux QFG at the soil surface is adjusted, and the remaining water flux not able to be met by the first soil layer is assigned to the variable WLOST.
In the remainder of the 250 loop, checks are carried out for soil layers with liquid water contents already effectively equal to \(\theta_{l, min}\). If the flux at the top of the layer is upward and that at the bottom of the layer is downward, both are set to zero. If both are downward, only the flux at the bottom is set to zero; if both are upward, only the flux at the top is set to zero. If the soil layer is an organic one and the liquid water content is less than the layer retention capacity THLRET and the flux at the bottom of the layer is downward, it is set to zero.
In the 300 loop, checks are carried out for soil layers with liquid water contents already effectively equal to THLMAX. If the flux at the top of the layer is downward and that at the bottom of the layer is upward, both are set to zero. If both are downward, then if the top flux is greater than the bottom flux, it is set equal to the bottom flux. If both are upward, then if magnitude of the bottom flux is greater than that of the top flux, it is set equal to the top flux.
In the 400 loop, for each soil layer THTEST is recalculated as the liquid water content resulting from the updated boundary fluxes, and is compared with a residual value THLTHR. For organic soil layers deeper than the first layer, THLTHR is set to the minimum of \(\theta_{l, ret}\) and \(\theta_l\), since only evapotranspiration can cause the soil moisture to fall below the retention capacity. (The first layer is excepted because surface evaporation can drive its moisture content below \(\theta_{l, ret}\).) For mineral soils, \(\theta_{l, min}\) is set to THLMIN. If THTEST < THLTHR, then if the flow at the bottom of the soil layer is downward, it is recalculated as the sum of the flow at the top of the layer plus the amount of water that must be removed to make the liquid water content of the layer equal to THLTHR. If the flow at the bottom of the layer is upward, the flow at the top of the layer is recalculated as the sum of the flow at the bottom of the layer minus the amount of water that must be removed to make the liquid water content of the layer equal to THLTHR. THTEST of the current layer is then reset to THLTHR, and THTEST of the overlying and underlying layers (if any) are recalculated using the new interface fluxes.
In the 500 loop, for each soil layer THTEST is compared with THLMAX. If THTEST > THLMAX, two temporary variables are defined: WLIMIT as the amount of water required to raise the liquid water content of the soil layer to THLMAX, and WEXCES as the amount of water representing the excess of THTEST over THLMAX. If the flux at the top of the layer is downward and the flux at the bottom of the layer is upward, then if the negative of the flux at the bottom of the layer is greater than WLIMIT, it is set equal to -WLIMIT and the flux at the top is set to zero; otherwise WEXCES is subtracted from the flux at the top. If both the flux at the top and the flux at the bottom are downward, then WEXCES is subtracted from the flux at the top. If both are upward, then WEXCES is added to the flux at the bottom, and a correction is applied to the flux at the bottom of the underlying layer. Finally, THTEST is recalculated for all the soil layers using the new interface fluxes.
In the first part of the 600 loop, a correction is performed in the unlikely event that as a result of the above adjustments, the flux at the bottom of the last permeable soil layer has become upward. If this occurs, all of the fluxes above are corrected by the inverse of this same amount. However, this will result in a small spurious upward water flux at the surface. Since the liquid water flows are now accounted for, an attempt is made, as in loop 250, to remove the required water from the frozen water store in the first layer. The excess that cannot be supplied from this source is assigned to WLOST.
In the second part of the 600 loop, the temperatures TFDT of the water fluxes at the top and bottom of the layers in the permeable soil profile are set equal to the temperatures of the water in the first and last layers respectively. The flow at the bottom of the profile and the temperature of the flow are used to update BASFLW and TBASFL, the baseflow from the current subarea and its temperature, and RUNOFF and TRUNOF, the total runoff from the grid cell in question and its temperature.
In the 700 loop, for each successive soil layer the temperature of the water flux at the bottom of the layer is first assigned. If the flux is downward, TFDT is set to the temperature of the water in the current layer; if it is upward, TFDT is set to the temperature of the water in the layer below. The temperature of the water in the current layer is then updated using the FDT and TFDT values at the top and bottom of the layer, and the liquid water content of the layer is set to THTEST.
1.8.13