Revisiting water retention curves for simple hydrological modelling of peat

ABSTRACT

Accurate modelling of peat water contents (θ) is critical for wetland studies. We modified the Campbell and Van Genuchten soil water retention curves (SWRCs) by replacing their empirical parameters with measurable properties. Combining the water table depth (d_WT) into SWRCs, we derived formulae for calculating volumetric θ from d_WT, coded in a simple model to test our hypotheses that d_WT is a reliable predictor of θ for peat of low and high water holding capacity at near-saturation. We compared our simulations with time-domain reflectometry θ measurements at Mer Bleue bog (Ontario, Canada) and the Western Peatland fen (Alberta, Canada). Constraining Campbell SWRC at extreme drying and waterlogging rather than wilting point and field capacity produced superior results. The Van Genuchten SWRC was approximated by hyperbolic and inverse hyperbolic segments. When simplified, its performance reconciles with that of the modified Campbell SWRC. Overall, our formulae performed well with generalized peat parameters.

KEYWORDS:

• peat
• soil water retention curves
• water contents
• water table
• modelling

Editor A. Fiori Associate editor L. Ruiz

1 Introduction

Northern peatlands store large amounts of carbon (C) over time; estimates of the current stock vary between 200 and 450 Pg globally (Gorham 1991, Roulet et al. 2007, Loisel et al. 2014). This vast C store is sensitive to peat hydrology (Clymo 1983, Lafleur et al. 2003). Peat water status is critical for C sequestration in plants and soil (Strack et al. 2006, Sulman et al. 2009), as it affects plant–water relations and productivity (Bubier et al. 2003a, 2003b) and microbial activity along aerobic–anaerobic gradients in peat (Lafleur et al. 2005a, 2005b). Thus, accurate quantification of soil water contents (θ) in peat and better understanding of their dynamics will improve our understanding of hydrological controls on peatland processes, including microbial respiration, bryophyte and vascular productivity and respiration, net ecosystem exchange of C, and peatland–atmosphere energy exchange (Peichl et al. 2014).

However, peatland hydrology is still not well understood largely due to the variability of peat properties and complexity of peat profiles, with sharp transitions between the poorly decomposed top peat with large air-filled porosity that facilitates water infiltration and more decomposed and heavily waterlogged deep peat that impedes drainage (Holden 2009, Weber et al. 2017a, 2017b). Varying thicknesses of well-aerated peat on top and waterlogged peat below also play a role in the complexity of peatland hydrology. Furthermore, peat from different peatlands, or different zones of the same profile, can differ widely in its near-saturation water holding capacity (Mezbahuddin et al. 2016). Peat of low water holding capacity at near-saturation (Cai et al. 2010, Sulman et al. 2010) rapidly loses water with small water table (WT) drawdowns, while peat of high water holding capacity at near-saturation retains large amounts of water with similar WT fluctuations (Parmentier et al. 2009).

1.1 Overview of peat hydrological modelling

Ecohydrological models complement field observations for better understanding of peatland hydrology. Many of them simulate θ by solving the water balance of the ecosystem from incoming precipitation, outgoing evapotranspiration, and incoming and outgoing surface runoff, subsurface lateral and vertical drainage, and recharge. Examples of such models are the Deng and Wang model (Deng and Wang 2017), ECOSYS (Grant 2001), HYDRUS (Schwarzel et al. 2006), FLOCOPS (Kennedy and Price 2004), COUP (Jansson and Karlberg 2004), ECOUL (Garnier et al. 1997), Visual MODFLOW and MODFLOW (McDonald and Harbaugh 1988). Many of these models require data-intensive (hourly, daily) continuous meteorological inputs of precipitation, air temperature, relative humidity, and wind speed (Dimitrov et al. 2014a, 2014b) that often are not taken at the modelled site, and thus subject to uncertainty. Also, these models require many hydrological and soil parameters, such as saturated hydraulic conductivities that may differ in vertical and lateral directions (Reeve et al. 2000), evapotranspiration forcing, aspects and terrain slopes, and largely unknown field boundary conditions (Grant 2001). Such parameter values may vary widely among peatlands or even within the same peatland (Fraser 1999).

All the above make model calibration of one- and three-dimensional (1D and 3D) ecohydrological models time- and effort-consuming. Furthermore, the empirical parameters that shape the soil water retention curves (SWRCs) (Fig. 1), i.e. the curve-fitting parameter in the hyperbolic Campbell SWRC, or the generally unknown pore distribution and shape parameters in the sigmoidal Van Genuchten SWRC, are derived from observations on specific sites or often as aggregated values from pooled records (Letts et al. 2000). These parameters are not measurable quantities, and when applied to other sites they introduce inevitable uncertainties.

Figure 1. Generalized soil water retention curves (SWRC), hyperbolic Campbell SWRC and sigmoidal Van Genuchten SWRC with high and low inflection points. The Campbell SWRC and low-inflection sigmoidal Van Genuchten SWRC are characteristic for peat of low water retention at near-saturation, while the high-inflection sigmoidal Van Genuchten SWRC is characteristic for peat of high water retention at near-saturation. The key soil parameters that determine the SWRC curves are the pairs of {total porosity (θ_p); water potential at saturation (–Ψ_m,sat)}, {water-filled porosity at field capacity (θ_FC); water potential at field capacity (–Ψ_m,FC)}, {water-filled porosity at wilting point (θ_WP); water potential at wilting point (–Ψ_m,WP)}, {residual soil water content (θ_r); hydroscopic potential (–Ψ_m,HC, i.e. the soil matric potential in the air-dry condition)}, and only for Van Genuchten SWRC {water-filled porosity at the inflection point (θ_LX); water potential at the inflection point (–Ψ_m,LX)}

Simple ecohydrological models simulate ecosystem hydrology from site-measured water table depth (d_WT) used as a model input (Kurbatovam et al. 2009, St-Hilaire et al. 2010). Peat θ values in these models are reported to be calculated from d_WT, interpreted to be the soil matric water potential (Ψ_m) expressed as height above the WT using linear functions (Zhang et al. 2002, Barr et al. 2012) or via the Campbell SWRC (Frolking et al. 2002). A shortcoming is that these models do not explicitly simulate the dynamics of θ, as observed in the field. Also, such models rarely account for the macropore fraction in peat, which is reported to reach rather large volumetric values in the near-surface layers (Baird 1997, Schwarzel et al. 2002, Wallage and Holden 2011), thus, ultimately affecting the soil water retention and aeration.

1.2 Rationale for revisiting SWRC for peat

Relying upon a simplifying assumption that water in the unsaturated peat is always in equilibrium between gravitational and matric tension due to the high peat porosity, Frolking et al. (2002) suggested a modified version of the original Campbell SWRC function, i.e.:

${\theta }_{\mathrm{i}}={\theta }_{\mathrm{p}}{\left(\frac{d_{\mathrm{W}\mathrm{T}}-d_{\mathrm{i}}}{{\psi }_{\mathrm{m},\mathrm{s}\mathrm{a}\mathrm{t}}}\right)}^{-\frac{1}{\beta }}$ with  $\left(d_{\mathrm{W}\mathrm{T}}-d_{\mathrm{i}}\right)>{\psi }_{\mathrm{m},\mathrm{s}\mathrm{a}\mathrm{t}}$

(1) ${\theta }_{\mathrm{i}}={\theta }_{\mathrm{p}}\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}\left(d_{\mathrm{W}\mathrm{T}}-d_{\mathrm{i}}\right)\le {\psi }_{\mathrm{m},\mathrm{s}\mathrm{a}\mathrm{t}}$(1)

where the soil matric potential is expressed as the difference between d_WT (m) and depth d_i (m) to a soil layer i, Ψ_m,sat is the air-entry matric potential (m), θ_i is volumetric soil water content (m³ m⁻³), θ_p is the soil total porosity (m³ m⁻³), and β is an empirical curve-fitting parameter (dimensionless). Their approach allows the calculation of θ_i directly from the SWRC, thereby avoiding the complexity of calculating water fluxes driven by water potential gradients, constrained by hydraulic conductivities, and subject to continuity of water flow. However, the disadvantages of Equation (1)(1) ${\theta }_{\mathrm{i}}={\theta }_{\mathrm{p}}\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}\left(d_{\mathrm{W}\mathrm{T}}-d_{\mathrm{i}}\right)\le {\psi }_{\mathrm{m},\mathrm{s}\mathrm{a}\mathrm{t}}$(1) are the empirical parameter β obtained by curve-fitting and the restricted applicability only to peat of low water holding capacity at near-saturation. Also, Equation (1)(1) ${\theta }_{\mathrm{i}}={\theta }_{\mathrm{p}}\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}\left(d_{\mathrm{W}\mathrm{T}}-d_{\mathrm{i}}\right)\le {\psi }_{\mathrm{m},\mathrm{s}\mathrm{a}\mathrm{t}}$(1) does not account for the large peat macropore fractions θ_M (m³ m⁻³) declining in depth.

Some ecohydrological process-based models, such as ECOSYS (Grant 2001), use a log-transformed version of the Campbell SWRC function (Fig. 1), which allows for the substitution of the empirical β with physically measurable parameters. In particular, ECOSYS sets the log-transformed line of Campbell SWRC by solving for the line slope and intercept by water potentials and water contents at field capacity, Ψ_m,FC and θ_FC, and wilting point, Ψ_m,WP and θ_WP. The log-transformed Campbell SWRC worked well for peat of low water holding capacity at near-saturation (Dimitrov et al. 2010a). However, it was much less accurate for peat with high water holding capacity at near-saturation (Mezbahuddin et al. 2016). For this reason, Mezbahuddin et al. (2016) modified ECOSYS to use Van Genuchten SWRC for peat with high water holding capacity while using Campbell SWRC for peat with low water holding capacity. Despite such accommodations, empirical parameters cannot be avoided in the original Van Genuchten SWRC, and setting the log-transformed line of Campbell SWRC between field capacity (Ψ_m,FC; θ_FC) and wilting point (Ψ_m,WP; θ_WP) proved challenging for peat (Dimitrov et al. 2010a). While applicable for mineral soils, where θ varies predominantly between θ_FC and θ_WP, θ in peat often varies outside that range, rapidly shifting between the extremes at near saturation and drying, Thus, calculating peat θ in the segments to the right of θ_FC and to the left of θ_WP (Fig. 1) that are not well defined has greater uncertainty compared to the internal line segment between θ_FC and θ_WP.

1.3 Objectives and hypotheses

Our study had three objectives. The first is to develop formulae for simple calculation of peat θ from d_WT by modifying the Campbell and Van Genuchten SWRCs to avoid their empirical parameters and replace them by quantifiable physical soil properties that could be measured on site or found in the literature as representative for a variety of site conditions. The second is to code these formulae in a newly developed simple model, named the Peat Hydrological Model (PHM), that reduces the complexity of peat hydrology simulations by calculating dynamic peat θ at any user-specified soil depth from dynamic d_WT observations. The third is to test PHM for a variety of peat types and properties, and peatland microtopography forms within the same site and between peatland sites.

The main (conceptual) hypothesis of this study is that d_WT, which is closely related to soil volumetric water content, can be a reliable predictor of peat θ at various soil depths, in peat of low and high water holding capacity at near-saturation. The modelling hypothesis implemented by the PHM model is that log-transformed Campbell and Van Genuchten SWRCs corrected for peat macroporosity and parameterized for peat type (i.e. fibric, hemic, and sapric) will improve simulations of peat θ at different depths predicted from daily d_WT for a variety of peatland conditions. Therefore, we test the general validity of the simplifying assumption that highly porous unsaturated peat reaches equilibrium within a day (Frolking and Crill 1994, Frolking et al. 2002). We used PHM to address five specific research hypotheses:

1. The log-transformed line of the Campbell SWRC in peat of low water holding capacity at near-saturation should be determined by the pairs of soil hydrological parameters that confine its largest hyperbolic segment, i.e. {total porosity (θ_p); water potential at saturation (Ψ_m,sat)} and {residual soil water content (θ_r); hydroscopic potential (Ψ_m,HC – that is, the soil matric potential in the air-dry condition)} (Fig. 1). Our expectation was that {Ψ_m,sat; θ_p} and {Ψ_m,HC; θ_r} would determine the Campbell SWRC reasonably well, as these pairs are well constrained for peat θ variation in the curve segments around saturation and drying, as well as in between.

2. Alternatively, we expected that setting the log-transformed line of Campbell SWRC between field capacity (Ψ_m,FC; θ_FC) and wilting point (Ψ_m,WP; θ_WP) would be less accurate than setting it between {Ψ_m,sat; θ_p} and {Ψ_m,HC; θ_r} for peat subject to extreme drying and waterlogging.

3. The original sigmoidal Van Genuchten SWRC can be approximated by a series of hyperbolic (Campbell-type) segments, starting with an inverse hyperbolic segment confined between {θ_p; Ψ_m,sat} and {the water content θ_LX at the inflection point, i.e. the point at which the change in the direction of curvature occurs; the water potential Ψ_m,LX at the inflection point}, followed by Campbell-type segments between {θ_LX; Ψ_m,LX} and {θ_FC; Ψ_m,FC}, {θ_FC; Ψ_m,FC} and {θ_WP; Ψ_m,WP}, {θ_WP; Ψ_m,WP} and {θ_r; Ψ_m,HC} (Fig. 1, from right to left). Our expectation was that the modified Van Genuchten SWRC with a high inflection point (high-flex) (Fig. 1) would perform well for peat of high water holding capacity at near-saturation.

4. The modified Van Genuchten SWRC can be further simplified for peat with low water holding capacity at near-saturation by approximation with only two hyperbolic segments, i.e. a relatively small inverse segment between {θ_p; Ψ_m,sat} and {θ_LX; Ψ_m,LX} on the right, and a relatively large segment to the left from {θ_LX; Ψ_m,LX} (Fig. 1). It was expected that both modified Campbell and Van Genuchten with low inflection point (low-flex) SWRCs (Fig. 1) would perform well for peat of low water holding capacity at near-saturation.

5. Generalized peat properties of {θ_p; Ψ_m,sat} and {θ_r; Ψ_m,HC} for fibric, hemic, and sapric peat, aggregated for different peatlands and reported in independent studies, provide working parameter sets for both the modified Campbell and Van Genuchten SWRCs, which can be applied to different peatlands, peatland microtopographic forms (hummocks, hollows), and types of peat by degree of decomposition (fibric, hemic, sapric).

2 Materials and methods

2.1 Simulating soil water contents in peat

The newly developed formulae for calculating volumetric θ_i at any position i within the soil profile at a given depth d_i (m) below the soil surface are developed mathematically in the Supplementary material for this paper. Here we report only the relevant final equations involved in PHM simulations, where all θ, Ψ and d values are given in m³ m⁻³, MPa, and m, respectively. Following Frolking et al. (2002) at a daily time step, originally based on Frolking and Crill (1994):

(2) ${\psi }_{\mathrm{m},\mathrm{i}}=\frac{-\left(d_{\mathrm{W}\mathrm{T}}-d_{\mathrm{i}}\right)}{100}$(2)

Then the hyperbolical Campbell SWRC (Fig. 1) was log-transformed and parameterized by {θ_p; Ψ_m,sat} and {θ_r; Ψ_m,HC}, thereby constrained between extreme drying and waterlogging:

(3) ${\theta }_{\mathrm{i}}=\left(\frac{{\theta }_{\mathrm{p}}-{\theta }_{\mathrm{M}}}{{\theta }_{\mathrm{p}}}\right){{\theta }_{\mathrm{r}}}^{\left[1-\frac{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{H}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{i}}\right)}{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{H}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{s}\mathrm{a}\mathrm{t}}\right)}\right]}{{\theta }_{\mathrm{p}}}^{\left[\frac{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{H}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{i}}\right)}{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{H}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{s}\mathrm{a}\mathrm{t}}\right)}\right]}$(3)

which was modified for peat by plugging in Equation (2)(2) ${\psi }_{\mathrm{m},\mathrm{i}}=\frac{-\left(d_{\mathrm{W}\mathrm{T}}-d_{\mathrm{i}}\right)}{100}$(2) . Also, the log-transformed hyperbolical Campbell SWRC was parameterized by {θ_FC; Ψ_m,FC} and {θ_WP; Ψ_m,WP}, and thereby constrained between field capacity and wilting point too:

(4) ${\theta }_{\mathrm{i}}=\left(\frac{{\theta }_{\mathrm{p}}-{\theta }_{\mathrm{M}}}{{\theta }_{\mathrm{p}}}\right){{\theta }_{\mathrm{F}\mathrm{C}}}^{\left[1-\frac{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{F}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{i}}\right)}{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{F}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{W}\mathrm{P}}\right)}\right]}{{\theta }_{\mathrm{W}\mathrm{P}}}^{\left[\frac{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{F}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{i}}\right)}{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{F}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{W}\mathrm{P}}\right)}\right]}$(4)

which was modified for peat by plugging in Equation (2)(2) ${\psi }_{\mathrm{m},\mathrm{i}}=\frac{-\left(d_{\mathrm{W}\mathrm{T}}-d_{\mathrm{i}}\right)}{100}$(2) .

Unlike the hyperbolical Campbell SWRC, the sigmoidal Van Genuchten SWRC cannot be log-transformed directly to avoid the unknown empirical parameters. Instead, we suggest that its curvature could be approximated by several different hyperbolic segments in the outlined solution below, i.e. between {θ_r; Ψ_m,HC} and {θ_WP; Ψ_m,WP}, {θ_WP; Ψ_m,WP} and {θ_FC; Ψ_m,FC}, {θ_FC; Ψ_m,FC} and {θ_LX; Ψ_m,LX}, and an inverse hyperbolic segment between {θ_LX; Ψ_m,LX} and {θ_p; Ψ_m,sat} to the right of the inflection point where the curvature changes (Fig. 1). Therefore:

(5a) ${\theta }_{\mathrm{i}}=\left(\frac{{\theta }_{\mathrm{p}}-{\theta }_{\mathrm{M}}}{{\theta }_{\mathrm{p}}}\right){{\theta }_{\mathrm{r}}}^{\left[1-\frac{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{H}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{i}}\right)}{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{H}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{W}\mathrm{P}}\right)}\right]}{{\theta }_{\mathrm{W}\mathrm{P}}}^{\left[\frac{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{H}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{i}}\right)}{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{H}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{W}\mathrm{P}}\right)}\right]}$(5a)

(5b) ${\theta }_{\mathrm{i}}=\left(\frac{{\theta }_{\mathrm{p}}-{\theta }_{\mathrm{M}}}{{\theta }_{\mathrm{p}}}\right){{\theta }_{\mathrm{W}\mathrm{P}}}^{\left[1-\frac{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{W}\mathrm{P}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{i}}\right)}{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{W}\mathrm{P}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{F}\mathrm{C}}\right)}\right]}{{\theta }_{\mathrm{F}\mathrm{C}}}^{\left[\frac{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{W}\mathrm{P}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{i}}\right)}{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{W}\mathrm{P}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{F}\mathrm{C}}\right)}\right]}$(5b)

(5c) ${\theta }_{\mathrm{i}}=\left(\frac{{\theta }_{\mathrm{p}}-{\theta }_{\mathrm{M}}}{{\theta }_{\mathrm{p}}}\right){{\theta }_{\mathrm{F}\mathrm{C}}}^{\left[1-\frac{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{F}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{i}}\right)}{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{F}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{L}\mathrm{X}}\right)}\right]}{{\theta }_{\mathrm{L}\mathrm{X}}}^{\left[\frac{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{F}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{i}}\right)}{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{F}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{L}\mathrm{X}}\right)}\right]}$(5c)

(5d) $\begin{array}{rl}& {\theta }_{\mathrm{i}}=\left(\frac{{\theta }_{\mathrm{p}}-{\theta }_{\mathrm{M}}}{{\theta }_{\mathrm{p}}}\right)\left\{{\left[\frac{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{L}\mathrm{X}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{s}\mathrm{a}\mathrm{t}}\right)}{\mathrm{l}\mathrm{n}\left({\theta }_{\mathrm{L}\mathrm{X}}\right)-\mathrm{l}\mathrm{n}\left({\theta }_{\mathrm{p}}\right)}\right]}^2\right.\\ & \left.\left\{\frac{1}{{\theta }_{\mathrm{L}\mathrm{X}}}-\frac{{{\theta }_{\mathrm{p}}}^{\left[\frac{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{i}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{L}\mathrm{X}}\right)}{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{L}\mathrm{X}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{s}\mathrm{a}\mathrm{t}}\right)}\right]}}{{{\theta }_{\mathrm{L}\mathrm{X}}}^{\left[\frac{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{i}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{s}\mathrm{a}\mathrm{t}}\right)}{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{L}\mathrm{X}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{s}\mathrm{a}\mathrm{t}}\right)}\right]}}\right\}+{\theta }_{\mathrm{L}\mathrm{X}}\right\}\end{array}$(5d)

when Ψ_i is between {θ_r; Ψ_m,HC} and {θ_WP; Ψ_m,WP}, {θ_WP; Ψ_m,WP} and {θ_FC; Ψ_m,FC}, {θ_FC; Ψ_m,FC} and {θ_LX; Ψ_m,LX}, and {θ_LX; Ψ_m,LX} and {θ_p; Ψ_m,sat}, respectively, which could be further adapted for peat by plugging in Equation (2)(2) ${\psi }_{\mathrm{m},\mathrm{i}}=\frac{-\left(d_{\mathrm{W}\mathrm{T}}-d_{\mathrm{i}}\right)}{100}$(2) . This method provides an opportunity to reconcile the Campbell SWRC and Van Genuchten SWRC for soil substrate, such as peat of low water holding capacity at near-saturation. For this purpose, the segment confined between {θ_r; Ψ_m,HC} and {θ_LX; Ψ_m,LX} of the low-flex Van Genuchten SWRC (Fig. 1) was approximated as a single hyperbolic segment, i.e.:

(6a) ${\theta }_{\mathrm{i}}=\left(\frac{{\theta }_{\mathrm{p}}-{\theta }_{\mathrm{M}}}{{\theta }_{\mathrm{p}}}\right){{\theta }_{\mathrm{r}}}^{\left[1-\frac{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{H}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{i}}\right)}{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{H}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{L}\mathrm{X}}\right)}\right]}{{\theta }_{\mathrm{L}\mathrm{X}}}^{\left[\frac{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{H}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{i}}\right)}{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{H}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{L}\mathrm{X}}\right)}\right]}$(6a)

(6b) $\begin{array}{rl}& {\theta }_{\mathrm{i}}=\left(\frac{{\theta }_{\mathrm{p}}-{\theta }_{\mathrm{M}}}{{\theta }_{\mathrm{p}}}\right)\left\{{\left[\frac{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{L}\mathrm{X}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{s}\mathrm{a}\mathrm{t}}\right)}{\mathrm{l}\mathrm{n}\left({\theta }_{\mathrm{L}\mathrm{X}}\right)-\mathrm{l}\mathrm{n}\left({\theta }_{\mathrm{p}}\right)}\right]}^2\right.\\ & \left.\left\{\frac{1}{{\theta }_{\mathrm{L}\mathrm{X}}}-\frac{{{\theta }_{\mathrm{p}}}^{\left[\frac{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{i}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{L}\mathrm{X}}\right)}{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{L}\mathrm{X}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{s}\mathrm{a}\mathrm{t}}\right)}\right]}}{{{\theta }_{\mathrm{L}\mathrm{X}}}^{\left[\frac{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{i}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{s}\mathrm{a}\mathrm{t}}\right)}{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{L}\mathrm{X}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{s}\mathrm{a}\mathrm{t}}\right)}\right]}}\right\}+{\theta }_{\mathrm{L}\mathrm{X}}\right\}\end{array}$(6b)

when Ψ_i is between {θ_r; Ψ_m,HC} and {θ_LX; Ψ_m,LX}, and {θ_LX; Ψ_m,LX} and {θ_p; Ψ_m,sat}, respectively, which was further adapted for peat by plugging in Equation (2)(2) ${\psi }_{\mathrm{m},\mathrm{i}}=\frac{-\left(d_{\mathrm{W}\mathrm{T}}-d_{\mathrm{i}}\right)}{100}$(2) . Equations (3(3) ${\theta }_{\mathrm{i}}=\left(\frac{{\theta }_{\mathrm{p}}-{\theta }_{\mathrm{M}}}{{\theta }_{\mathrm{p}}}\right){{\theta }_{\mathrm{r}}}^{\left[1-\frac{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{H}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{i}}\right)}{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{H}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{s}\mathrm{a}\mathrm{t}}\right)}\right]}{{\theta }_{\mathrm{p}}}^{\left[\frac{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{H}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{i}}\right)}{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{H}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{s}\mathrm{a}\mathrm{t}}\right)}\right]}$(3) –6) are all constrained by:

(7) ${\theta }_{\mathrm{i}}\in \left[{\theta }_{\mathrm{r}};\left.{\theta }_{\mathrm{p}}\right)\right.with{d}_{\mathrm{W}\mathrm{T}}>d_{\mathrm{i}}and{\theta }_{\mathrm{i}}={\theta }_{\mathrm{p}}with{d}_{\mathrm{W}\mathrm{T}}\le d_{\mathrm{i}}$(7)

Equations (2(2) ${\psi }_{\mathrm{m},\mathrm{i}}=\frac{-\left(d_{\mathrm{W}\mathrm{T}}-d_{\mathrm{i}}\right)}{100}$(2) –7(7) ${\theta }_{\mathrm{i}}\in \left[{\theta }_{\mathrm{r}};\left.{\theta }_{\mathrm{p}}\right)\right.with{d}_{\mathrm{W}\mathrm{T}}>d_{\mathrm{i}}and{\theta }_{\mathrm{i}}={\theta }_{\mathrm{p}}with{d}_{\mathrm{W}\mathrm{T}}\le d_{\mathrm{i}}$(7) ) are the same as Equations (S3), (S4a), (S7a), (S9a), (S9b), (S9c), (S9d), (S10a), (S10b) and (S5) in the Supplementary material, respectively. Applying Equations (2)(2) ${\psi }_{\mathrm{m},\mathrm{i}}=\frac{-\left(d_{\mathrm{W}\mathrm{T}}-d_{\mathrm{i}}\right)}{100}$(2) , (3(3) ${\theta }_{\mathrm{i}}=\left(\frac{{\theta }_{\mathrm{p}}-{\theta }_{\mathrm{M}}}{{\theta }_{\mathrm{p}}}\right){{\theta }_{\mathrm{r}}}^{\left[1-\frac{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{H}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{i}}\right)}{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{H}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{s}\mathrm{a}\mathrm{t}}\right)}\right]}{{\theta }_{\mathrm{p}}}^{\left[\frac{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{H}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{i}}\right)}{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{H}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{s}\mathrm{a}\mathrm{t}}\right)}\right]}$(3) ), (4(4) ${\theta }_{\mathrm{i}}=\left(\frac{{\theta }_{\mathrm{p}}-{\theta }_{\mathrm{M}}}{{\theta }_{\mathrm{p}}}\right){{\theta }_{\mathrm{F}\mathrm{C}}}^{\left[1-\frac{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{F}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{i}}\right)}{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{F}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{W}\mathrm{P}}\right)}\right]}{{\theta }_{\mathrm{W}\mathrm{P}}}^{\left[\frac{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{F}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{i}}\right)}{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{F}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{W}\mathrm{P}}\right)}\right]}$(4) ), (6) and (7(7) ${\theta }_{\mathrm{i}}\in \left[{\theta }_{\mathrm{r}};\left.{\theta }_{\mathrm{p}}\right)\right.with{d}_{\mathrm{W}\mathrm{T}}>d_{\mathrm{i}}and{\theta }_{\mathrm{i}}={\theta }_{\mathrm{p}}with{d}_{\mathrm{W}\mathrm{T}}\le d_{\mathrm{i}}$(7) ) for peat of low water-holding capacity at near-saturation, and Equations (2)(2) ${\psi }_{\mathrm{m},\mathrm{i}}=\frac{-\left(d_{\mathrm{W}\mathrm{T}}-d_{\mathrm{i}}\right)}{100}$(2) , (5) and (7(7) ${\theta }_{\mathrm{i}}\in \left[{\theta }_{\mathrm{r}};\left.{\theta }_{\mathrm{p}}\right)\right.with{d}_{\mathrm{W}\mathrm{T}}>d_{\mathrm{i}}and{\theta }_{\mathrm{i}}={\theta }_{\mathrm{p}}with{d}_{\mathrm{W}\mathrm{T}}\le d_{\mathrm{i}}$(7) ) for peat of high water-holding capacity at near-saturation, θ at any depth is calculated from daily d_WT (m). The calculations are computerized by coding these formulae into the simple model code of PHM, written in C++, which reads daily d_WT values for periods of simulations and uses pairs of {θ_r; Ψ_m,HC}, {θ_WP; Ψ_m,WP}, {θ_FC; Ψ_m,FC}, {θ_p; Ψ_m,sat} and/or {θ_LX; Ψ_m,LX} measured on site or reported in the literature to produce output for daily θ at any user-specified depth d_i. The PHM accounts for θ_M too, if any macropore fraction occurs at the corresponding d_i.

2.2 Sites and measurements

Two sites, Mer Bleue bog, with peat of low water retention at near-saturation (Lafleur et al. 2005a, 2005b), and the Western Peatland fen (WPL), with peat of high water retention at near-saturation (Flanagan and Syed 2011), were chosen for this study. Mer Bleue bog, located near Ottawa, Ontario, Canada, is an ombrotrophic peatland. It is slightly domed, with peat depth varying from ~0.3 m along the periphery up to ~6 m towards the central part, with an average depth of ~4 m across most of the site. The bog surface has a pronounced hummock–hollow microtopography, with hummocks occupying ~70% and hollows ~30% of the surface area (Fig. 2). Based on Von Post degree of humification, the least decomposed and highly microporous, well-drained fibric peat occupies the upper ~0.35 m in hummocks, and the upper ~0.1 m in hollows. Below the fibric peat, a thin layer of transitional hemic peat spreads across hummocks and hollows between ~0.35 and ~0.45 m depth below the hummock surface. Farther below, the most decomposed, waterlogged sapric peat exists down to the contact with mineral soil.

Figure 2. The hummocky microtopography of Mer Bleue bog, Ontario, Canada. The hummocks of ~1 m mean diameter occupy about 70% of the surface area, alternating with the adjacent hollows which occupy the remaining ~30%, with average difference of ~0.25 m between the top of the hummock domes and the hollow surface

Daily WT was measured at Mer Bleue bog during the period 1998–2004 and varied from ~0.7 m to ~0.25 m below the hummock surface, and from ~0 m to ~0.5 m in hollows. The d_WT referenced to the hummock surface was recorded continuously in wells, using a float and counterweight system attached to a potentiometer, and complemented by manual d_WT measurements to verify the potentiometer readings (Lafleur et al. 2005a, 2005b). The θ was measured continuously at 0.1 and 0.2 m (fibric peat, dry), 0.3 m (fibric peat, wet), 0.4 m (hemic peat) and 0.5 m (sapric peat) below the hummock surface, and 0.03 m (fibric peat) and 0.15 m (hemic peat) below the hollow surface using water content probes (model CS615, Campbell Scientific). The calibration and normalization of the probes for Mer Bleue bog are described in Lafleur et al. (2001, 2003, 2005a).

The WPL fen is located in the boreal region of Alberta, Canada, about 80 km northeast of the town of Athabasca. Its richer nutrient status resulted in wide range of mosses, e.g. Sphagnum mosses, feather moss, and brown moss, which resulted in more decomposed and less microporous fibric peat, compared to that at Mer Bleue bog, occupying the top 20–30 cm of the hummock domes of ~2 m diameter. The peat depth around the location of the instruments is ~2 m, and daily WT below the average hummock surface was measured in a well using a custom float and counterweight system attached to a potentiometer. Peat θ at 0.075, 0.010 and 0.125 m below the peat surface were measured continuously with water content probes (model CS616, Campbell Scientific) that were previously calibrated in the lab in containers of commercial peat moss held at a range of water contents (Syed et al. 2006, Flanagan and Syed 2011).

2.3 Model tests and statistics

To address our hypotheses, several model runs were conducted using PHM, each at a daily time step and continuously through the entire time series of available data for Mer Bleue bog (1998–2004) and the WPL fen (2005–2007). A summary of the features of each PHM run is provided in Table 1, including the model run name, the type of SWRC and hydrological parameters to constrain it, the type of peat and the site where applied, and the specific hypotheses that each model run was designed to test.

Table 1. Summary of the features of different peat hydrological model (PHM) runs and the hypotheses that they address. SWRC: soil water retention curve, described by a specific equation; period: the entire time series in years throughout which each model run was accomplished continuously at a daily time step: LWHC: low water holding capacity at near-saturation; HWHC: high water holding capacity at near-saturation; Ψ_m,sat and θ_p: total porosity and water potential at saturation; Ψ_m,HC and θ_r: residual water content and hygroscopic potential; Ψ_m,FC and θ_FC: water content and water potential at field capacity; Ψ_m,WP and θ_WP: water content and water potential at wilting point; and Ψ_m,LX and θ_LX: water content and water potential at inflection point

PHM model run	SWRC, Equation	Hydrological parameters^1,2 (Sources)	Peat type	Site	Period	Hypotheses
CAMP1	Modified Campbell, Equation (3)(3) ${\theta }_{\mathrm{i}}=\left(\frac{{\theta }_{\mathrm{p}}-{\theta }_{\mathrm{M}}}{{\theta }_{\mathrm{p}}}\right){{\theta }_{\mathrm{r}}}^{\left[1-\frac{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{H}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{i}}\right)}{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{H}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{s}\mathrm{a}\mathrm{t}}\right)}\right]}{{\theta }_{\mathrm{p}}}^{\left[\frac{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{H}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{i}}\right)}{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{H}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{s}\mathrm{a}\mathrm{t}}\right)}\right]}$(3)	Ψm,sat (Letts et al. 2000) and θp (Letts et al. 2000), Ψm,HC (Letts et al. 2000) and θr (Thien and Graveel 2003)	LWHC	Mer Bleue bog	1998–2004	1, 5
CAMP2	Modified Campbell, Equation (4)(4) ${\theta }_{\mathrm{i}}=\left(\frac{{\theta }_{\mathrm{p}}-{\theta }_{\mathrm{M}}}{{\theta }_{\mathrm{p}}}\right){{\theta }_{\mathrm{F}\mathrm{C}}}^{\left[1-\frac{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{F}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{i}}\right)}{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{F}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{W}\mathrm{P}}\right)}\right]}{{\theta }_{\mathrm{W}\mathrm{P}}}^{\left[\frac{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{F}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{i}}\right)}{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{F}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{W}\mathrm{P}}\right)}\right]}$(4)	Ψm,FC (Letts et al. 2000) and θFC (Dimitrov et al. 2010a), Ψm,WP (Letts et al. 2000) and θWP (Dimitrov et al. 2010a)	LWHC	Mer Bleue bog	1998–2004	2
VANG1	Modified Van-Genuchten, Equation (6)	Ψm,sat (Letts et al. 2000) and θp (Letts et al. 2000), Ψm,LX (this study) and θLX (this study), Ψm,HC (Letts et al. 2000) and θr (Thien and Graveel 2003)	LWHC	Mer Bleue bog	1998–2004	4, 5
VANG2	Modified Van-Genuchten, Equation (5)	Ψm,sat (Letts et al. 2000) and θp (Letts et al. 2000), Ψm,LX (this study) and θLX (this study), Ψm,FC (Letts et al. 2000) and θFC (Mezbahuddin et al. 2016), Ψm,WP (Letts et al. 2000) and θWP (Mezbahuddin et al. 2016), Ψm,HC (Letts et al. 2000) and θr (Thien and Graveel 2003)	HWHC	Western Peatland fen	2005–2007	3, 5
VANG2,JAS	Modified Van-Genuchten, Equation (5)	Ψm,sat (Letts et al. 2000) and θp (Letts et al. 2000), Ψm,LX (this study) and θLX (this study), Ψm,FC (Letts et al. 2000) and θFC (Mezbahuddin et al. 2016), Ψm,WP (Letts et al. 2000) and θWP (Mezbahuddin et al. 2016), Ψm,HC (Letts et al. 2000) and θr (Thien and Graveel 2003)	HWHC	Western Peatland fen	2005–2007 June, July, August only	3, 5

¹Pairs of hydrological parameters, each consisting of water potential (MPa) and soil water content (m³ m⁻³), to constrain the SWRC.

²Simulated soil water contents account for macroporosity (m³ m⁻³) at Mer Bleue bog (Dimitrov et al. 2010a) and Western Peatland fen (Mezbahuddin et al. 2016).

2.3.1 Simulations for peat of low water holding capacity at near-saturation

The model runs CAMP₁, CAMP₂ and VANG₁, designed to test Hypotheses 1, 2 and 4 for peat of low water holding capacity at near saturation (Table 1), were parameterized, initialized, and run for the peat profiles of hummocks and hollows at Mer Bleue bog (Fig. 3). The hummock and hollow profiles consist of 10 and eight layers, respectively. The peat layers were selected to encompass the profile of the dry fibric peat (the top 25 cm in hummocks and the top 5 cm in hollows), the wet fibric peat (25–35 cm in hummocks, 5–20 cm in hollows), the transitional hemic peat (35–45 cm in hummocks, 10–20 cm in hollows), and the waterlogged sapric peat (below 45 cm in hummocks and below 20 cm in hollows), as well as to capture the zone of WT fluctuations during 1998–2004. For comparison purposes, the peat stratigraphy in Fig. 3 was also set to be consistent with previous modelling studies at the same site (Dimitrov et al. 2010a, 2010b).

Figure 3. Schematic representation of Mer Bleue bog peat profile, with alternating hummocks (hmk) and hollows (hlw). n_hmk and n_hlw: numbers of soil layers in hummocks and hollows; d_i,hmk and d_i,hlw: depths to soil layer bottom (with respect to hummock and hollow surfaces); θ_M: volumetric macropore fraction; Ψ_m,sat: water potential at saturation; θ_p: total porosity; Ψ_m,HC: hydroscopic potential (i.e. the soil matric potential in the air-dry condition); θ_r: residual soil water content; Ψ_m,FC and θ_FC, and Ψ_m,WP and θ_WP: the water potentials and water contents at field capacity and wilting point, respectively; and Ψ_m,LX and θ_LX: water potential and water content at the inflection point (i.e. where change in the direction of curvature of the water retention curve occurs). The peat hydrological parameters are assigned accordingly to the model runs CAMP₁ and CAMP₂ (modified Campbell soil water retention curve), and VANG₁ (modified Van Genuchten soil water retention curve)

The alternating wetter and drier years in the 1998–2004 period at Mer Bleue bog provided an opportunity to examine how well the PHM performed during a sequence of hydrologically diverse years, including the year 2000 with the highest average and minimum annual d_WT of 0.36 and 0.45 m, respectively, then 2002 with the lowest average and minimum annual d_WT of 0.46 and 0.74 m, respectively. The short-term model response to rainfall was tested against measured θ and d_WT following an extreme rain event of 124.5 mm on 9 September 2004. Simulated θ by CAMP₁, CAMP₂ and VANG₁ were compared with observed values at depths of 0.1, 0.2, 0.3, 0.4 and 0.5 m in hummocks, and 0.03 and 0.15 m in hollows. The model performance of VANG₁ and CAMP₂ was compared against that of CAMP₁.

2.3.2 Simulations for peat of high water holding capacity at near-saturation

The model runs VANG₂ and VANG_2,JAS, i.e. the version of VANG₂ only for the months of June, July and August to avoid the effects of spring thaw and/or fall freeze not currently handled by PHM, were designed to test Hypothesis 3 for peat of high water holding capacity at near saturation (Table 1) and were parameterized, initialized, and run for peat (only hummock available) at the WPL fen (Fig. 4). The hummock peat profile (no available data for hollows) consists of 15 layers and was set to encompass the profile of the fibric peat (the top 19 cm), and the waterlogged sapric peat (below 45 cm in hummocks and below 20 cm in hollows), as well as to capture the zone of WT fluctuations during that period. The peat stratigraphy is consistent with previous studies conducted to simulate the hydrology of the same site (Mezbahuddin et al. 2016). The 2005–2007 period at the WPL fen provided an opportunity to examine PHM’s ability to simulate seasonal and inter-annual θ dynamics during a sequence starting with the wet year 2005 with the highest average and minimum annual d_WT (0.34 m and 0.38 m), followed by the average year 2006 with moderate average and minimum annual d_WT (0.42 m and 0.53 m), and then the dry year 2007 with low average and minimum annual d_WT (0.53 m and 0.65 m). Simulated θ by VANG₂ and VANG_2,JAS were compared with observed values at depths of 0.075, 0.010 and 0.0125 m.

Figure 4. Schematic representation of Western Peatland fen peat profile (hummocks). n: number of soil layer; d_i: depth to soil layer bottom (with respect to hummock surface); θ_M: volumetric macropore fraction; Ψ_m,sat: water potential at saturation; θ_p: total porosity; Ψ_m,HC: hydroscopic potential (i.e. the soil matric potential in the air-dry condition); θ_r: residual soil water content; Ψ_m,FC and θ_FC, and Ψ_m,WP and θ_WP: the water potentials and water contents at field capacity and wilting point, respectively; and Ψ_m,LX and θ_LX: water potential and water content at the inflection point (i.e. where change in the direction of the curvature of the water retention curve occurs). The peat hydrological parameters are assigned accordingly to the model runs VANG₂ and VANG_2,JAS (modified Van Genuchten soil water retention curve)

2.3.3 Testing the validity of generalized sets of peat properties

Hypothesis 5 was tested with CAMP₁, VANG₁, VANG₂ and VANG_2,JAS model runs applied for different peatlands, microtopographic forms and types of peat using generalized values (Figs 3 and 4) aggregated from independent research on peat from different sites (Table 1). The θ_p and θ_r values for the fibric, hemic and sapric peat were taken from Letts et al. (2000). The values of Ψ_m,sat were taken from the same literature source; however, they were slightly modified to reconcile Ψ_m,sat and Van Genuchten α values in fibric peat, where hydrological changes and θ measurements were most frequent, given that α = 1/Ψ_m,sat (Mezbahuddin et al. 2016). We paired to θ_r the general value of soil matric potential in the air-dry condition (hydroscopic coefficient) Ψ_m,HC (Thien and Graveel 2003). The other parameters in Table 1 and Figs 3 and 4 were aligned with previous modelling studies on Mer Bleu bog (Dimitrov et al. 2010a, 2010b) and WPL fen (Mezbahuddin et al. 2016).

2.3.4 Statistics of model–data comparisons

To evaluate performance of each model run and to compare runs, linear regressions of simulated θ (y-axis) on observed θ (x-axis) and model accuracy statistics were obtained at depths of 0.1, 0.2, 0.3, 0.4 and 0.5 m in hummocks, and 0.03 and 0.15 m in hollows at Mer Bleue bog, and at depths of 0.075, 0.100 and 0.125 m in hummocks at WPL fen. The overall model accuracy was addressed by the root mean square error (RMSE) between simulations and observations and its systematic (RMSE_s) and unsystematic (RMSE_u) components (Willmott 1981). The RMSE is a measure of the average difference between observed and modelled values, expressed as follows (Ritter and Munoz-Carpena 2013):

(8a) $\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=\sqrt{\sum _{i=1}^n\frac{{\left(P_i-O_i\right)}^2}{n}}$(8a)

where P_i and O_i represent the sample (of size n). The relative discrepancies between observations and simulations were evaluated by Willmott’s index of agreement d, which represents the relative model accuracy with respect to the magnitude of observed values, as follows (Willmott 1982):

(8b) $d=1-\frac{{\sum }_{i=1}^n{\left(P_i-O_i\right)}^2}{{\sum }_{i=1}^n{\left(|P_i^{\prime }\left|+|O_i^{\prime }\right|\right)}^2},with{P}_i^{\prime }=P_i-\bar{O}and{O}_i^{\prime }=O_i-\bar{O}$(8b)

where 0 ≤ d ≤ 1 and $\bar{O}$ is the average of the observed values. In addition, the Nash-Sutcliffe efficiency coefficient was calculated to assess the model’s predictive power, as follows (Krause et al. 2005):

(8c) $\mathrm{N}\mathrm{S}\mathrm{E}=1-\frac{{\sum }_{i=1}^n{\left(P_i-O_i\right)}^2}{{\sum }_{i=1}^n{\left(\bar{O}-O_i\right)}^2}$(8c)

Further to the above, the 95% confidence intervals (CI_95%) were evaluated for RMSE and NSE by applying the statistical hydrological software FITEVAL (Ritter and Munoz-Carpena, 2013, 2020).

Model runs CAMP₁, VANG₁ and CAMP₂ for Mer Bleue bog were evaluated for their individual performance and then compared with each other. Also, the model runs VANG₂ and VANG_2,JAS (only for the months of July, August, September) for the WPL fen were evaluated for their performance and compared with each other, to avoid the complicating effects of spring thawing on θ in the peat of high water holding capacity at near-saturation that are not currently handled by PHM.

3 Results

3.1 Modelling the hydrology of peat of low water holding capacity at near-saturation

3.1.1 Simulations by modified Campbell SWRC constrained at extreme drying and waterlogging (CAMP₁), and simplified Van Genuchten SWRC (VANG₁)

Both CAMP₁ and VANG₁ simulated closely seasonal and inter-annual dynamics of θ in hummocks at depth of 0.1, 0.2, 0.3, 0.4 and 0.5 m below the hummock surface (Figs 5–9) and in hollows at depth of 0.03 and 0.15 m below the hollow surface (Figs 10 and 11), all of them reflecting the WT dynamics. Both runs captured changes in θ at 0.2 m (Fig. 6) and its pronounced dynamics at 0.3 m (Fig. 7), with peaks between June and August 2000, April and May 2001, April, May and July 2002, May and June 2003, and to a lesser extent between April and June 2004 associated with large d_WT fluctuations. Both CAMP₁ and VANG₁ followed closely the observed θ at 0.4 m, thereby capturing the magnitude and dynamics of θ with WT drawdown during the summer (Fig. 8). In addition, the model runs captured the magnitude of summer drop in θ at 0.5 m while maintaining conditions close to saturation in agreement with observations for most of the time series (Fig. 9). In addition, both CAMP₁ and VANG₁ captured θ dynamics in the upper dry hollow peat at the very top 0.03 m with pronounced peaks between June and July 2000, April and May 2001, April, May and July 2002, and May and June 2003 (Fig. 10). They also captured the θ dynamics in the deeper wet peat at 0.15 m with pronounced summer drawdowns (Fig. 11). Both runs performed better compared to previous modelling studies at the same site during the spring thaw and summer drought (see Dimitrov et al. 2010a, fig. 7(d)).

Figure 5. Daily simulated and observed soil water contents (θ) at 10 cm depth in hummocks and daily observed water table (WT) depths referred to the hummock surface at Mer Bleue bog. Simulated θ by CAMP₁, VANG₁, and CAMP₂ model runs were driven by the available observed daily WT records during the period 1998–2004

Figure 6. Daily simulated and observed soil water contents (θ) at 20 cm depth in hummocks and daily observed water table (WT) depths referred to the hummock surface at Mer Bleue bog. Simulated θ by CAMP₁, VANG₁, and CAMP₂ model runs were driven by the available observed daily WT records during the period 1998–2004

Figure 7. Daily simulated and observed soil water contents (θ) at 30 cm depth in hummocks and daily observed water table (WT) depths referred to the hummock surface at Mer Bleue bog. Simulated θ by CAMP₁, VANG₁, and CAMP₂ model runs were driven by the available observed daily WT records during the period 1998–2004

Figure 8. Daily simulated and observed soil water contents (θ) at 40 cm depth in hummocks and daily observed water table (WT) depths referred to the hummock surface at Mer Bleue bog. Simulated θ by CAMP₁, VANG₁, and CAMP₂ model runs were driven by the available observed daily WT records during the period 1998–2004

Figure 9. Daily simulated and observed soil water contents (θ) at 50 cm depth in hummocks and daily observed water table (WT) depths referred to the hummock surface at Mer Bleue bog. Simulated θ by CAMP₁, VANG₁, and CAMP₂ model runs were driven by the available observed daily WT records during the period 1998–2004

Figure 10. Daily simulated and observed soil water contents (θ) at 3 cm depth in hollows and daily observed water table (WT) depths referred to the hollow surface at Mer Bleue bog. Simulated θ by CAMP₁, VANG₁, and CAMP₂ model runs were driven by the available observed daily WT records during the period 1998–2004

Figure 11. Daily simulated and observed soil water contents (θ) at 15 cm depth in hollows and daily observed water table (WT) depths referred to the hollow surface at Mer Bleue bog. Simulated θ by CAMP₁, VANG₁, and CAMP₂ model runs were driven by the available observed daily WT records during the period 1998–2004

The most intensive rain event in the 1998–2004 period produced 124.5 mm over 24 h at Mer Bleue bog on September 9 in 2004. Consistent with observations, CAMP₁ and VANG₁ captured the timing and magnitude of θ changes at depth, i.e. the slight peak in hummock θ at 0.2 m (Fig. 6) caused by the rainfall, and rapid waterlogging to full saturation at 0.3 m during and immediately after the rainfall (Fig. 7), which was coincident with an immediate rise and gradual decline of d_WT. Faster drawdown of simulated θ at 0.3 m compared to observations was possibly due to slightly higher water retention in the field, which might be due to some peat hysteresis effects that PHM does not account for at present.

Various statistics for both CAMP₁ and VANG₁ of observed on modelled daily θ at 0.1 m, 0.2, 0.3, 0.4 and 0.5 m in hummocks, and 0.03 and 0.15 m in hollows are given in Table 2, along with their CI_95% values. All linear regressions were significant at α = 0.05 with p < 2e-16 (Table 2). At 0.1 m depth in hummocks, R² = 28% indicated that 28% of variation in θ observations was explained. The slope and intercept indicated that simulations matched the observations with no pronounced bias. Although the RMSE between 0.02 and 0.03 was low, the low values of θ in highly macroporous dry fibric peat resulted in relatively small differences between the average observed θ and observed and simulated values, which in turn resulted in a relatively smaller denominator with respect to the numerator in Equation (8b)(8b) $d=1-\frac{{\sum }_{i=1}^n{\left(P_i-O_i\right)}^2}{{\sum }_{i=1}^n{\left(|P_i^{\prime }\left|+|O_i^{\prime }\right|\right)}^2},with{P}_i^{\prime }=P_i-\bar{O}and{O}_i^{\prime }=O_i-\bar{O}$(8b) , hence low d (0.31). Systematic underestimation of θ by both CAMP₁ and VANG₁ resulted in a large sum of squares, hence a large numerator and negative NSE values in Equation (8c)(8c) $\mathrm{N}\mathrm{S}\mathrm{E}=1-\frac{{\sum }_{i=1}^n{\left(P_i-O_i\right)}^2}{{\sum }_{i=1}^n{\left(\bar{O}-O_i\right)}^2}$(8c) . However, the relatively small RMSE_u and higher RMSE_s approaching RMSE indicated that model parameterization would need improvement rather than model formulae. At 0.2 m depth in hummocks, slopes and intercepts indicated that CAMP₁ and VANG₁ were unbiased and explained almost half of the variation of θ observations with R² = 0.46. Compared to 0.1 m depth, greater d (0.71) indicated that both runs followed observations even more closely. The NSE become positive with the upper limits of its CI_95% of 0.76, reaching the zone of satisfactory model performance (Moriasi et al. 2007, Ritter and Munoz-Carpena 2013). The relatively small values of RMSE between 0.02 and 0.06 and its systematic and unsystematic components indicated good model accuracy.

Table 2. Linear regressions of simulated values on measured daily volumetric soil water contents (m³ m⁻³), Willmott’s index of agreement (d), root mean square error (RMSE), its systematic and unsystematic components (RMSE_s and RMSE_u) and Nash-Sutcliffe efficiency (NSE) at depths of 0.1, 0.2, 0.3, 0.4 and 0.5 m in hummocks, and 0.03 and 0.15 m in hollows at Mer Bleue bog over the 1998–2004 period. All regressions are significant at α = 0.05 (p < 2e-16). CI_95% indicates the 95% confidence interval. N: number; R²: coefficient of determination

	N	R^2	Slope	Intercept	d	RMSE (CI95%)	RMSEs	RMSEu	NSE (CI95%)
	(-)	(-)	(-)	(m^3 m^−3)	(-)	(m^3 m^−3)	(m^3 m^−3)	(m^3 m^−3)	(-)
Model run: CAMP1
Hummock, 0.1 m	990	0.28	0.97	0.03	0.31	0.03 [0.02; 0.03]	0.03	0.01	−8.23 (−14.09; −4.95]
Hummock, 0.2 m	1462	0.46	0.97	0.02	0.71	0.04 [0.02; 0.06]	0.02	0.03	0.20 [−0.67; 0.76]
Hummock, 0.3 m	1212	0.51	0.71	0.04	0.84	0.17 [0.11; 0.22]	0.07	0.15	0.60 [0.23; 0.78]
Hummock, 0.4 m	1270	0.75	0.78	0.17	0.93	0.10 [0.08; 0.13]	0.04	0.08	0.69 [0.44; 0.85]
Hummock, 0.5 m	969	0.34	1.04	0.01	0.62	0.07 [0.06; 0.11]	0.04	0.05	0.23 [−0.78; 0.67]
Hollow, 0.03 m	1237	0.66	0.81	0.03	0.89	0.09 [0.06; 0.12]	0.03	0.09	0.71 [0.58; 0.80]
Hollow, 0.15 m	720	0.66	0.80	0.19	0.89	0.11 [0.08; 0.14]	0.05	0.10	0.62 [−0.18; 0.82]
Model run: VANG1
Hummock, 0.1 m	990	0.28	0.99	0.03	0.30	0.03 [0.02; 0.04]	0.03	0.01	−8.47 [−14.47; −5.13]
Hummock, 0.2 m	1462	0.46	0.98	0.02	0.71	0.04 [0.02; 0.06]	0.02	0.03	0.19 [−0.72; 0.76]
Hummock, 0.3 m	1212	0.51	0.71	0.04	0.84	0.17 [0.11; 0.23]	0.07	0.15	0.60 [0.23; 0.77]
Hummock, 0.4 m	1270	0.74	0.76	0.19	0.92	0.10 [0.07; 0.13]	0.05	0.09	0.66 [0.42; 0.82]
Hummock, 0.5 m	969	0.32	0.85	0.16	0.64	0.07 [0.05; 0.10]	0.04	0.06	0.15 [−0.94; 0.65]
Hollow, 0.03 m	1237	0.65	0.80	0.04	0.89	0.09 [0.07; 0.13]	0.03	0.09	0.71 [0.57; 0.79]
Hollow, 0.15 m	720	0.67	0.78	0.20	0.89	0.11 [0.08; 0.15]	0.05	0.10	0.61 [−0.29; 0.83]
Model run: CAMP2
Hummock, 0.1 m	990	0.28	18.57	−0.13	0.07	0.09 [0.08; 0.09]	0.09	0.01	−45.93 [−64.84; −33.93]
Hummock, 0.2 m	1462	0.36	1.12	0.09	0.23	0.10 [0.09; 0.11]	0.09	0.03	−4.47 [−19.85; −1.64]
Hummock, 0.3 m	1212	0.50	0.78	0.10	0.80	0.17 [0.13; 0.23]	0.08	0.15	0.35 [−0.05; 0.59]
Hummock, 0.4 m	1270	0.72	0.72	0.24	0.90	0.12 [0.09; 0.15]	0.08	0.09	0.53 [0.20; 0.76]
Hummock, 0.5 m	969	–	< 8e-31	0.29	0.32	0.07 [0.06; 0.09]	0.04	0.08	−0.07 [−0.38; 0.08]
Hollow, 0.03 m	1237	0.63	0.87	0.11	0.80	0.13 [0.11; 0.16]	0.09	0.10	0.16 [−0.58; 0.50]
Hollow, 0.15 m	720	0.63	0.72	0.28	0.85	0.13 [0.10; 0.18]	0.09	0.10	0.41 [−1.12; 0.82]

In the zone of greatest WT fluctuations in wet fibric peat at 0.3 m in hummocks, both CAMP₁ and VANG₁ followed closely the observations explaining >50% of their variation (R² = 0.51). While RMSE variation between 0.11 and 0.22 (CAMP₁) or 0.23 (VANG₁) was greatest at that depth, the d went up to 0.84, indicating improved model performance compared to 0.1 and 0.2 m depths above. Furthermore, NSE = 0.6 with its upper limits of CI_95% of 0.78 (CAMP₁) or 0.77 (VANG₁) were > 0.5, thus indicating satisfactory model performance at that depth (Moriasi et al. 2007). At 0.4 m depth, the PHM performance improved further compared to the depths above, with lower RMSE and RMSE_u, d > 0.9, and NSE > 0.65 with upper limits of CI_95% > 0.8, thereby indicating good model performance (Ritter and Munoz-Carpena 2013). The CAMP₁ and VANG₁ explained ~75% of the variation of θ observations with R² of 0.75 (CAMP₁) or 0.74 (VANG₁). Relatively large RMSE_s for both runs implied the potential to improve the PHM performance at 0.3 m and 0.4 m in hummocks by improving the model parameterization, especially with {θ_LX; Ψ_m,LX} at inflection points (Figs 1 and 3). The smaller WT variation in peat at 0.5 m in hummocks that is predominantly at near saturation resulted in lower R², d, and NSE values, compared to those at 0.3 m and 0.4 m, although the slopes approached unity and intercepts were near zero for both model runs. Both CAMP₁ and VANG₁ performed comparatively well at hollows for fibric peat at 0.03 m and hemic peat at 0.15 m with R² ≥ 0.65, d values of 0.89, and NSE values > 0.6 or 0.7, respectively, indicating satisfactory model performance (Moriasi et al. 2007, Ritter and Munoz-Carpena 2013).

3.1.2 Simulations by modified Campbell SWRC constrained at wilting point and field capacity (CAMP₂)

The PHM run CAMP₂ with Mer Bleue data was conducted to show the effects of parameterizing PHM with Ψ and θ at field capacity and wilting point in peat of low water holding capacity at near-saturation, instead of Ψ and θ at waterlogging and extreme drying used to parameterize PHM for CAMP₁ (Fig. 3). During the entire period 1998–2004, CAMP₂ underestimated θ of dry fibric peat at 0.1 and 0.2 m in hummocks (Figs 5 and 6) and at 0.03 m in hollows (Fig. 10), while maintaining permanently waterlogged conditions at 50 cm depth in hummocks (Fig. 9). However, CAMP₂ performed well in the intermediate zone of wet fibric peat at 0.3 m in hummocks (Fig. 7) and 0.15 m in hollows (Fig. 11), as well as in lower hemic peat at 0.4 m in hummocks (Fig. 8). For predominantly dry fibric peat at 0.1 m in hummocks and 0.03 m in hollows, and for predominantly waterlogged sapric peat at 0.5 m in hummocks, where θ often varied outside of the range between field capacity and wilting point, CAMP₂ performed much worse than CAMP₁. CAMP₁ was also better able to capture the extremely low and high θ in both near-surface and deep peat. At 0.4 m in hummocks and 0.15 m in hollows, where θ usually varied between field capacity and wilting point, CAMP₂ performed similarly to CAMP_1.

Linear regressions suggested that CAMP₂ performed comparably well to CAMP₁ at 0.3 and 0.4 m in hummocks, as well as at 0.03 and 0.15 m in hollows, i.e. within the zone of pronounced WT fluctuations and θ varying greatly between the wilting point and field capacity (Table 2). However, the slopes and intercepts for CAMP₂ showed strong biases for θ in hummocks at 0.1 and 0.2 m (underestimating observations), and at 0.5 m (overestimating observations). The corresponding RMSE and d values indicated low model accuracy at those depths. On the other hand, much higher RMSE_s values compared to RMSE_u (similar to RMSE_u of CAMP₁) indicated that the model–data misfit was mostly due to inappropriate model parameterization by {θ_WP; Ψ_m,WP} and {θ_FC; Ψ_m,FC}, instead of {θ_r; Ψ_m,HC} and {θ_p; Ψ_m,sat}, which worked well in CAMP₁. The latter was confirmed by generally low NSE values that were <0.5, except at 0.4 m depth in hummocks, where NSE = 0.53, indicating unsatisfactory CAMP₂ model performance (Moriasi et al. 2007, Ritter and Munoz-Carpena 2013).

3.2 Modelling hydrology of peat of high water holding capacity at near-saturation (VANG₂ and VANG_2,JAS)

The PHM runs VANG₂ and its June, July and August version VANG_2,JAS closely simulated the observed seasonal and inter-annual dynamics of θ in WPL hummocks at depths of 0.075, 0.1 and 0.125 m below the hummock surface and associated WT fluctuations during 2005–2007 (Fig. 12). Specifics of local peat properties at these depths (Fig. 4) produced slightly more waterlogged conditions at 0.1 m than above (0.0075 m) and below (0.125 m), which was captured in simulations of θ at the corresponding depths. Inconsistencies in VANG₂ performances were mainly associated with periods of spring thaw and excessive evapotranspiration with summer drying in the field, as currently PHM does not simulate freeze–thaw cycles or evapotranspiration. The lack of simulated freeze–thaw cycles was partially compensated by VANG_2,JAS, which is a VANG₂ version confined to June, July and August.

Figure 12. Daily simulated and observed soil water contents (θ) at different depths in hummocks and daily observed water table (WT) depths referred to the hummock surface at Western Peatland fen. (a) Simulated θ by VANG₂ model run were driven by the available observed daily WT records during the period 2005–2007. (b) Simulated θ by VANG_2,JAS model runs were driven by the available observed daily WT records only in June, July and August during the period 2005–2007

The linear regressions for VANG₂ of observed on modelled daily θ at 0.075, 0.1 and 0.125 m were significant at α = 0.05, p < 2e-16 (Table 3). The R², slope, and intercept at 0.0075 and 0.125 m indicated that the simulations were unbiased and explained >33% of variation in θ observations. The 0.1 m depth had a lower R² and higher RMSE_s, possibly a result of local peat stratigraphic inconsistencies (Fig. 4), which produced some uncertainty in model parameterization at that depth compared to the peat above and below. Eliminating the effects of spring thaw in VANG_2.JAS resulted in pronounced improvement of model–observation statistics compared with VANG₂. The d values increased from ~0.70 to ~0.80 and the RMSE and RMSE_u decreased, indicating improved model accuracy. For VANG2, NSE was 0.42 (>0.4) at 0.075 m depth and 0.37 (>0.36) at 0.125 m depth, which according to some sources is acceptable (Ramanarayanan et al. 1997, Motovilov et al. 1999). However, the NSE values for all three soil depths improved for VANG_2,JAS and were above the widely accepted threshold of 0.5 (Moriasi et al. 2007). Furthermore, their upper CI_95% boundaries indicated PHM abilities for good (>0.8) and very good (>0.9) model performance (Ritter and Munoz-Carpena 2013).

Table 3. Linear regressions of simulated values on measured daily volumetric soil water contents (m³ m⁻³), Willmott’s index of agreement (d), root mean square error (RMSE), its systematic and unsystematic components (RMSE_s and RMSE_u) and Nash-Sutcliffe efficiency (NSE), at depths of 0.075, 0.100 and 0.125 m in hummocks at Western Peatland fen over the 2005–2007 period. All regressions are significant at α = 0.05 (p < 2e-16). CI_95% indicates the 95% confidence interval. N: number; R²: coefficient of determination

	N	R^2	Slope	Intercept	d	RMSE (CI95%)	RMSEs	RMSEu	NSE (CI95%)
	(-)	(-)	(-)	(m^3 m^−3)	(-)	(m^3 m^−3)	(m^3 m^−3)	(m^3 m^−3)	(-)
Model run: VANG2
Hummock, 0.075 m	682	0.33	1.01	−0.01	0.71	0.16 [0.12; 0.20]	0.01	0.16	0.42 [0.03; 0.74]
Hummock, 0.100 m	682	0.25	0.72	0.18	0.70	0.17 [0.13; 0.22]	0.05	0.16	0.31 [−0.05; 0.67]
Hummock, 0.125 m	682	0.35	0.95	–0.01	0.73	0.16 [0.13; 0.20]	0.03	0.16	0.37 [−0.10; 0.69]
Model run: VANG2,JAS
Hummock, 0.075 m	341	0.50	1.08	−0.03	0.80	0.14 [0.09; 0.19]	0.01	0.14	0.70 [0.32; 0.92]
Hummock, 0.100 m	341	0.45	0.68	0.22	0.80	0.14 [0.09; 0.18]	0.07	0.11	0.54 [0.14; 0.81]
Hummock, 0.125 m	341	0.50	0.98	−0.02	0.81	0.14 [0.10; 0.18]	0.02	0.14	0.62 [0.25; 0.84]

4 Discussion

4.1 Comparison to previous modelling work

The PHM with modified Campbell and Van Genuchten SWRCs, parameterized at waterlogging and extreme drying (CAMP₁ and VANG₁), performed similarly to or better than ECOSYS’s 3D ecohydrological model previously tested at Mer Bleue bog hummocks for 1998–2004 (Table 2 in this paper vs table 1 in Dimitrov et al. 2010a); there were no reported statistics for the performance of ECOSYS in hollows. ECOSYS parameterized at field capacity and wilting point partially offset the excessive drying effects in near-surface peat by explicit modelling of the evapotranspiration drag, thereby maintaining higher θ within the top peat compared to PHM (CAMP₂). Previous modelling studies on hydrology of the WPL fen employing ECOSYS (Mezbahuddin et al. 2016) did not explicitly report the statistics for model–data fit. Instead, R², d and RMSE values were bar-plotted for each year in the 2005–2007 period, which makes direct comparison to PHM (VANG₂) statistics impossible. However, R², d, RMSE, and RMSE_u values of VANG₂ and especially VANG_2,JAS (Table 3) fit within the range of the corresponding ECOSYS statistics.

4.2 Implications for SWRC of organic and mineral soils

Our newly developed formulae for modified Campbell and Van Genuchten SWRCs avoid the shape parameters β in Equation (3)(3) ${\theta }_{\mathrm{i}}=\left(\frac{{\theta }_{\mathrm{p}}-{\theta }_{\mathrm{M}}}{{\theta }_{\mathrm{p}}}\right){{\theta }_{\mathrm{r}}}^{\left[1-\frac{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{H}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{i}}\right)}{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{H}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{s}\mathrm{a}\mathrm{t}}\right)}\right]}{{\theta }_{\mathrm{p}}}^{\left[\frac{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{H}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{i}}\right)}{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{H}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{s}\mathrm{a}\mathrm{t}}\right)}\right]}$(3) and α and n in Equations (5) and (6) by using pairs of measurable hydrological parameters, i.e. {Ψ_m,sat; θ_p}, {Ψ_m,LX; θ_LX},{Ψ_m,FC; θ_FC}, {Ψ_m,WP; θ_WP}, {Ψ_m,WP; θ_WP} and {Ψ_m,HC; θ_r}. Thus, our formulae avoid complicated curve fitting for other sites that adds uncertainties to θ calculations. The generality of Equation (3)(3) ${\theta }_{\mathrm{i}}=\left(\frac{{\theta }_{\mathrm{p}}-{\theta }_{\mathrm{M}}}{{\theta }_{\mathrm{p}}}\right){{\theta }_{\mathrm{r}}}^{\left[1-\frac{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{H}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{i}}\right)}{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{H}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{s}\mathrm{a}\mathrm{t}}\right)}\right]}{{\theta }_{\mathrm{p}}}^{\left[\frac{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{H}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{i}}\right)}{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{H}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{s}\mathrm{a}\mathrm{t}}\right)}\right]}$(3) for the modified Campbell SWRC, and Equations (5) and (6) for the modified Van Genuchten SWRC and its simplified version, should be applicable for mineral soils too. The validity of Equations (3)(3) ${\theta }_{\mathrm{i}}=\left(\frac{{\theta }_{\mathrm{p}}-{\theta }_{\mathrm{M}}}{{\theta }_{\mathrm{p}}}\right){{\theta }_{\mathrm{r}}}^{\left[1-\frac{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{H}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{i}}\right)}{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{H}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{s}\mathrm{a}\mathrm{t}}\right)}\right]}{{\theta }_{\mathrm{p}}}^{\left[\frac{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{H}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{i}}\right)}{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{H}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{s}\mathrm{a}\mathrm{t}}\right)}\right]}$(3) , (5) and (6) for mineral soils could be tested directly vs {Ψ; θ} observations. Our formulae can further be tested for modelling soil water retention in 1D or 3D hydrological schemes of process-based models. The modified Van Genuchten SWRC, i.e. Equation (5), provides a new perspective for describing more complex SWRCs by introducing more inflection points, thereby creating even smaller sections if field data suggests they are necessary. Alternating hyperbolic and reverse hyperbolic segments would describe a variety of SWRC shapes (Kutilek and Nielsen 1994). This also implies that there is potential to expand application of our modified SWRCs, i.e. Equations (3)(3) ${\theta }_{\mathrm{i}}=\left(\frac{{\theta }_{\mathrm{p}}-{\theta }_{\mathrm{M}}}{{\theta }_{\mathrm{p}}}\right){{\theta }_{\mathrm{r}}}^{\left[1-\frac{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{H}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{i}}\right)}{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{H}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{s}\mathrm{a}\mathrm{t}}\right)}\right]}{{\theta }_{\mathrm{p}}}^{\left[\frac{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{H}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{i}}\right)}{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{H}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{s}\mathrm{a}\mathrm{t}}\right)}\right]}$(3) and (5), to non-homogeneous soil substrates, reclamation substrates that are often designed as mixtures of alternating organic and mineral soil layers (Daly et al. 2012, Wytrykush et al. 2012).

4.3 Advantages and uncertainties for peatland ecosystem modelling

The newly introduced PHM is a simple model for simulating peat water contents at any depth from field WT records, by applying either the Campbell or Van Genuchten SWRC, i.e. Equations (3)(3) ${\theta }_{\mathrm{i}}=\left(\frac{{\theta }_{\mathrm{p}}-{\theta }_{\mathrm{M}}}{{\theta }_{\mathrm{p}}}\right){{\theta }_{\mathrm{r}}}^{\left[1-\frac{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{H}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{i}}\right)}{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{H}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{s}\mathrm{a}\mathrm{t}}\right)}\right]}{{\theta }_{\mathrm{p}}}^{\left[\frac{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{H}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{i}}\right)}{\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{H}\mathrm{C}}\right)-\mathrm{l}\mathrm{n}\left(-{\psi }_{\mathrm{m},\mathrm{s}\mathrm{a}\mathrm{t}}\right)}\right]}$(3) , (5) and (6). The PHM functions well with generalized and site-specific peat hydrological parameters (Figs 3 and 4) and could be a reasonable substitution for more complex 1D and 3D hydrological schemes in ecohydrological modelling of peatlands. Compared to these more complex models, PHM is less data-intensive and is computationally more efficient. These economies are achieved largely by omitting variable parameters, such as saturated hydraulic conductivities, that vary in vertical and horizontal directions; terrain slopes and aspects; and largely unknown lateral and vertical field boundary conditions. In addition, PHM does not require as many meteorological input drivers as the complex process-based models do. Instead, available d_WT records and general peat parameters drive calculations of peat θ. Thus, PHM is applicable at large spatial scales and at many geographical locations, indicative of national and international wetland databases (Zoltai et al. 2000, Loisel et al. 2014) for which continuous (for dynamic simulations) or sporadic (for static estimations) WT records exist.

A limitation of PHM for large-scale spatial application is the lack of WT observations for many peatlands. However, the PHM could also be driven by d_WT estimated from large-scale satellite or airborne imagery (Kalacska et al. 2018) or provided by open-source web-based collections. The simple data requirements of PHM warrant its potential to easily connect to big data sets provided by open application programming interface (API) sources, which makes the model suitable for large-scale scientific simulations that benefit from rapid development of new information technologies and web applications.

Other current limitations of PHM are the lack of explicit processes for evapotranspiration, freeze–thaw, and summer drying. Even though their effects are partially captured in declining and rising WT levels, there is clearly a need for special treatment of these processes, as well as to introduce peat hysteresis, elasticity, and plasticity, in the next stage of the model’s development. The PHM performance for the upper peat of low water holding capacity at near-saturation at Mer Bleue bog was relatively less sensitive to spring thaw and summer droughts, compared to the WPL. This difference could be attributed to the smaller water holding capacity per unit volume of Mer Bleue fibric peat, with its smaller matrix and larger macropore volumetric fractions, compared to the WPL peat (Figs 3 and 4). An inherent uncertainty for all models is estimating peat macroporosity (Dimitrov et al. 2010a). Although peat macroporosities have been studied for open bogs (Blodau and Moore 2002), blanket peatlands (Holden 2009), and forested peatlands (Silins and Rothwell 1998), few studies report the macropore volumetric fractions distribution in fibric, hemic, and sapric peat. It has been shown that quantifying macropore fractions solely by size is rather arbitrary, while quantifying by water potential at preferential gravitational flow might be site-specific and generally uncertain (Skopp 1981, Baird 1997).

5 Conclusion

This study supported the main hypothesis that the WT level, as a focal point of all inputs and outputs of the hydrological balance, can be a reliable predictor of water contents at any depth in peats with low and high water holding capacity at near-saturation. The original Campbell and Van Genuchten SWRCs can be modified mathematically to avoid computationally expensive and site-specific curve-fitted parameters, which can be replaced by measurable hydrological properties reconciled for typical peat profiles. The modified Campbell SWRC constrained at excessive drying and waterlogging provided superior results when applied for peat of low water holding capacity at near-saturation (Hypothesis 1), compared to when constrained at field capacity and wilting point (Hypothesis 2). It was found that the field capacity and wilting point do not adequately constrain the Campbell SWRC hyperbola for the highly porous and dry fibric peat at near-surface and waterlogged sapric peat at depth. The sigmoidal Van Genuchten SWRC approximated by a series of inverse hyperbolic and hyperbolic segments was applicable to and performed well for peat of high water holding capacity at near-saturation (Hypothesis 3). When further simplified by means of only two segments, i.e. one small inverse hyperbolic to saturation and another large hyperbolic segment to excessive drying, the performance of the modified Van Genuchten SWRC for the peat of low water holding capacity at near-saturation was almost indiscernible from that of the Campbell SWRC (Hypothesis 4). New formulae for calculating peat water contents at any depth from WT depth have been developed and coded in a simple model for peat hydrology that works reasonably well for different types of peat and peatland sites. The new formulae performed well against independent field measurements for different peatland sites, when parameterized by generalized peat properties (Hypothesis 5), and may be applicable for large-scale simulations and soil substrate with diverse SWRC specifics.

Acknowledgements

Data were collected with funding from Canadian Carbon Program and Fluxnet Canada Research Network (FCRN) through its major sponsors the Natural Sciences and Engineering Council of Canada (NSERC), Canadian Foundation for Climate and Atmospheric Sciences, and BIOCAP Canada. The authors are greatly thankful to Prof. Lawrence B. Flanagan of the University of Lethbridge for sharing the hydrological data collection at Western Peatland fen. Also, the authors are thankful to Alberta Innovation and Advanced Education to Campus Alberta Innovation Program (CAIP) Research Chair (No.RCP-12-001-BCAIP) Prof. Junye Wang and Dean of Science and Technology at Athabasca University Prof. Lisa Carter for their support of the preliminary stages of this work. Furthermore, the authors wish to thank one anonymous reviewer and the associate editor for their detailed and useful comments and suggestions that resulted in improvement of this manuscript.

Disclosure statement

No potential conflict of interest was reported by the authors.

Supplementary material

Supplemental data for this article can be accessed here.