waterFlowInfiltrate.f90 File Reference

Quantifies movement of liquid water between soil layers under conditions of infiltration. More...

Functions/Subroutines
subroutine	waterflowinfiltrate (IVEG, THLIQ, THICE, TBARW, BASFLW, TBASFL, RUNOFF, TRUNOF, ZFAV, LZFAV, THLINV, QFG, WLOST, FI, EVAP, R, TR, TPOND, ZPOND, DT, ZMAT, WMOVE, TMOVE, THLIQX, THICEX, TBARWX, DELZX, ZBOTX, FDT, TFDT, PSIF, THLINF, GRKINF, THLMAX, THTEST, ZRMDR, FDUMMY, TDUMMY, THLDUM, THIDUM, TDUMW, TRMDR, ZF, FMAX, TUSED, RDUMMY, ZERO, WEXCES, FDTBND, WADD, TADD, WADJ, TIMPND, DZF, DTFLOW, THLNLZ, THLQLZ, DZDISP, WDISP, WABS, THPOR, THLRET, THLMIN, BI, PSISAT, GRKSAT, THLRAT, THFC, DELZW, ZBOTW, XDRAIN, DELZ, ISAND, IGRN, IGRD, IFILL, IZERO, LZF, NINF, IFIND, ITER, NEND, ISIMP, IGDR, IG, IGP1, IGP2, ILG, IL1, IL2, JL, N)

Detailed Description

Quantifies movement of liquid water between soil layers under conditions of infiltration.

Author

D. Verseghy, M. Lazare, Y. Delage, R. Soulis, J. Melton

Function/Subroutine Documentation

waterflowinfiltrate()

subroutine waterflowinfiltrate	(	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), intent(in)	BASFLW,
		real, dimension(ilg), intent(in)	TBASFL,
		real, dimension(ilg), intent(in)	RUNOFF,
		real, dimension(ilg), intent(in)	TRUNOF,
		real, dimension (ilg), intent(out)	ZFAV,
		integer, dimension (ilg), intent(out)	LZFAV,
		real, dimension(ilg), intent(out)	THLINV,
		real, dimension (ilg), intent(in)	QFG,
		real, dimension (ilg), intent(in)	WLOST,
		real, dimension (ilg), intent(in)	FI,
		real, dimension (ilg), intent(in)	EVAP,
		real, dimension (ilg), intent(in)	R,
		real, dimension (ilg), intent(in)	TR,
		real, dimension (ilg), intent(in)	TPOND,
		real, dimension (ilg), intent(in)	ZPOND,
		real, dimension (ilg), intent(in)	DT,
		real, dimension (ilg,igp2,igp1), intent(in)	ZMAT,
		real, dimension (ilg,igp2), intent(inout)	WMOVE,
		real, dimension(ilg,igp2), intent(inout)	TMOVE,
		real, dimension(ilg,igp1), intent(inout)	THLIQX,
		real, dimension(ilg,igp1), intent(inout)	THICEX,
		real, dimension(ilg,igp1), intent(inout)	TBARWX,
		real, dimension (ilg,igp1), intent(inout)	DELZX,
		real, dimension (ilg,igp1), intent(inout)	ZBOTX,
		real, dimension (ilg,igp1), intent(inout)	FDT,
		real, dimension (ilg,igp1), intent(inout)	TFDT,
		real, dimension (ilg,igp1), intent(inout)	PSIF,
		real, dimension(ilg,igp1), intent(inout)	THLINF,
		real, dimension(ilg,igp1), intent(inout)	GRKINF,
		real, dimension(ilg,ig), intent(inout)	THLMAX,
		real, dimension(ilg,ig), intent(inout)	THTEST,
		real, dimension (ilg,igp1), intent(inout)	ZRMDR,
		real, dimension(ilg,igp1), intent(inout)	FDUMMY,
		real, dimension(ilg,igp1), intent(inout)	TDUMMY,
		real, dimension(ilg,ig), intent(inout)	THLDUM,
		real, dimension(ilg,ig), intent(inout)	THIDUM,
		real, dimension (ilg,ig), intent(inout)	TDUMW,
		real, dimension(ilg), intent(inout)	TRMDR,
		real, dimension(ilg), intent(inout)	ZF,
		real, dimension (ilg), intent(in)	FMAX,
		real, dimension (ilg), intent(in)	TUSED,
		real, dimension(ilg), intent(inout)	RDUMMY,
		real, dimension (ilg), intent(in)	ZERO,
		real, dimension(ilg), intent(in)	WEXCES,
		real, dimension(ilg), intent(in)	FDTBND,
		real, dimension (ilg), intent(in)	WADD,
		real, dimension (ilg), intent(in)	TADD,
		real, dimension (ilg), intent(in)	WADJ,
		real, dimension(ilg), intent(in)	TIMPND,
		real, dimension (ilg), intent(in)	DZF,
		real, dimension(ilg), intent(in)	DTFLOW,
		real, dimension(ilg), intent(in)	THLNLZ,
		real, dimension(ilg), intent(in)	THLQLZ,
		real, dimension(ilg), intent(in)	DZDISP,
		real, dimension (ilg), intent(in)	WDISP,
		real, dimension (ilg), intent(in)	WABS,
		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)	THLRAT,
		real, dimension (ilg,ig), intent(in)	THFC,
		real, dimension (ilg,ig), intent(in)	DELZW,
		real, dimension (ilg,ig), intent(in)	ZBOTW,
		real, dimension(ilg), intent(in)	XDRAIN,
		real, dimension(ig), intent(in)	DELZ,
		integer, dimension(ilg,ig), intent(in)	ISAND,
		integer, dimension (ilg), intent(inout)	IGRN,
		integer, dimension (ilg), intent(in)	IGRD,
		integer, dimension(ilg), intent(inout)	IFILL,
		integer, dimension(ilg), intent(in)	IZERO,
		integer, dimension(ilg), intent(inout)	LZF,
		integer, dimension(ilg), intent(inout)	NINF,
		integer, dimension(ilg), intent(in)	IFIND,
		integer, dimension (ilg), intent(in)	ITER,
		integer, dimension (ilg), intent(in)	NEND,
		integer, dimension (ilg), intent(in)	ISIMP,
		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
	)

Parameters

[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]	basflw	Base flow from bottom of soil column [m]
[in]	tbasfl	Temperature of base flow from bottom of soil column [K]
[in]	runoff	Total runoff from soil column [m]
[in]	trunof	Temperature of total runoff from soil column [K]
[in]	qfg	Evaporation from soil surface (diagnostic) \([kg m^{-2} s^{-1}]\)
[in]	wlost	Residual amount of water that cannot be supplied by surface stores \([kg m^{-2}]\)
[out]	zfav	Average depth of wetting front over current time step [m]
[out]	thlinv	Liquid water content behind the wetting front \([m^3 m^{-3}]\)
[out]	lzfav	Soil layer index in which the average position of the wetting front occurs
[in]	fi	Fractional coverage of subarea in question on modelled area [ ]
[in]	evap	Evaporation rate from ground surface \([m s^{-1}]\)
[in]	r	Rainfall rate at ground surface \([m s^{-1}]\)
[in]	tr	Temperature of rainfall [C]
[in]	tpond	Temperature of ponded water [C]
[in]	zpond	Depth of ponded water on soil surface [m]
[in]	dt	Time period over which water movement takes place [s]
[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]	thlrat	Fractional saturation of soil behind the wetting front \([ ] (f_{inf})\)
[in]	thfc	Field capacity \([m^3 m^{-3}]\)
[in]	delzw	Permeable depth of soil layer \([m] (\Delta z_{g,w})\)
[in]	zbotw	Depth to permeable bottom of soil layer [m]
[in]	xdrain	Drainage index for water flow at bottom of soil profile [ ]
[in]	delz	Overall thickness of soil layer [m]
[in]	isand	Sand content flag
[in,out]	igrn	Flag to indicate whether calculations in this subroutine are to be done
[in]	igrd	Flag to indicate whether calculations in subroutine waterFlowNonInfiltrate are to be done
[in]	igdr	Index of soil layer in which bedrock is encountered

In loop 50, the flag IGRN 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; and that either the rainfall rate is greater than zero or that ponded water exists on the surface. If any of these conditions is not met, IGRN is set to zero.

In loop 100, a series of local arrays THLIQX, THICEX, TBARWX, DELZX and ZBOTX are defined for use in the subsequent infiltration calculations, with a “y” dimension of IG+1, where IG is the number of soil layers. The entries from 1 to IG are set equal to the corresponding values in THLIQ, THICE, TBARW, DELZW and ZBOTW. An effective saturated hydraulic conductivity GRKSATF is calculated for each soil layer, applying an empirical correction for the presence of ice. This ice content factor, \(f_{ice}\), is calculated from Zhao and Gray (1997) as:

\(f_{ice} = [1.0 – min(1.0, \theta_i / \theta_p)]^2\)

To further account for the presence of ice, a modified pore volume THPORF for the soil layer is calculated as the maximum of the available pore volume \(\theta_p\) – \(\theta_i\) (where \(\theta_p\) is the total pore volume and \(\theta_i\) is the ice content of the layer), the actual liquid water content \(\theta_l\), and the minimum residual liquid water content \(\theta_{l, min}\). (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 \(\theta_l\) or \(\theta_{l, min}\) may be greater than \(\theta_p\) – \(\theta_i\).) Finally, the water content THLINF and the hydraulic conductivity GRKINF behind the wetting front are evaluated, following the analysis of Green and Ampt (1911) [41] and treating the change in soil moisture due to infiltration as a downward-propagating square wave. THLINF is calculated as the maximum of \(f_{inf} ( \theta_p – \theta_i)\), \(\theta_l\), and \(\theta_{l, min}\), where \(f_{inf}\) represents the fractional saturation of the soil behind the wetting front, corresponding to a hydraulic conductivity of half the saturated value (this correction is applied in order to account for the fact that as infiltration occurs, a small amount of air is usually trapped in the soil). GRKINF is calculated from GRKSATF, THLINF and THPORF, using the classic Clapp and Hornberger (1978) [21] equation

\(K(z) = K_{sat} ( \theta_l / \theta_p)^{(2b + 3)}\)

where K(z) is the hydraulic conductivity at a depth z, \(K_{sat}\) is the hydraulic conductivity at saturation, \(\theta_p\) is the pore volume and b is an empirical coefficient.

In loop 150, values for the IG+1 entries in each of the above arrays are assigned. If the permeable thickness DELZW of the IG layer is less than its overall thickness DELZ (indicating the presence of bedrock) the IG+1 entries are set to 0; otherwise they are set to the values for the IG layer in the case of THLIQX, THICEX, TBARWX and THLINF, and for DELZX and ZBOTX they are set to large positive numbers. GRKINF is set to the value for the IG layer multiplied by the drainage parameter XDRAIN, which takes a value of 0 if the soil is underlain by an impermeable layer (as in a bog), and a value of 1 otherwise.

In loop 200, the soil water suction across the wetting front psi_f is calculated for each soil layer, using equation 25 in Verseghy (1991):

\(\Psi_f = b[\Psi_{inf} K_{inf} - \Psi(z) K(z)]/[K_{inf} (b+3)]\)

where \(\Psi_{inf}\) and \(K_{inf}\) are the soil moisture suction and the hydraulic conductivity behind the wetting front, and \(\Psi(z)\) and \(K(z)\) are the soil moisture suction and the hydraulic conductivity ahead of the wetting front. The soil moisture suction values are obtained using the classic Clapp and Hornberger (1978) equation:

\(\Psi(z) = \Psi_{sat} (\theta_l / \theta_p)^{(-b)}\)

where \(\Psi_{sat}\) is the soil moisture suction at saturation.

In loop 400, a test is carried out to determine whether saturated conditions are already present in the soil. Generally it is assumed that a period of unsaturated flow occurs initially, so the flag IFILL is set to 1, the depth of the wetting front ZF is set to zero, and the flag LZF indicating the index of the soil layer in which the wetting front occurs is set to 1. If a pond is present on the soil surface or if the hydraulic conductivity if the first soil layer is very small, it is deemed that saturated flow begins immediately, so IFILL is set to zero. Then each soil layer is checked in turn, to ascertain whether the liquid water content is greater than THLINF. If so, it is concluded that saturated flow is occurring in this layer; ZF is set to ZBOTW, the depth of the bottom of the soil layer, LZF is set to the index of the next layer, and IFILL is set to zero. For each layer found to be undergoing saturated flow, its water content is stored in the J+1 level of the water movement matrix WMOVE, and its water temperature is stored in the J+1 level of the matrix TMOVE. The counter NINF is set to the number of soil layers behind and including the one with the wetting front, plus 2.

A series of subroutines is called to complete the infiltration calculations. Subroutine waterInfiltrateUnsat treats the period of unsaturated flow up the point when saturated flow begins. Subroutine waterInfiltrateSat treats the period of saturated flow. Subroutine waterBaseflow reassigns liquid water contents and temperatures at the conclusion of the infiltration period. THLIQ, THICE and TBARW are then set to the corresponding updated values of THLIQX, THICEX and TBARWX. The average depth of the wetting front in the soil layer in which it occurs, the index of the soil layer, and the moisture content behind the wetting front are stored in output variables ZFAV, LZFAV and THLINV respectively. Finally, if there is any time remaining in the current time step following the infiltration period, subroutine waterFlowNonInfiltrate is called to treat the movement of soil water over the remainder of the time step.