7.2 Challenges of Numerical Modeling
The challenge to simulate water movement in the upper layers revolves around the exponential decline in saturated hydraulic conductivity. This is compounded by profound changes in saturation associated with macropores prevalent in the upper layers, especially in Sphagnum dominated systems. Evaporation from the drained moss surface generates extreme soil water pressures (< 10,000 cm H2O) that cannot be confidently parameterized in pressure-moisture-conductivity relationships needed to simulate flow. Moreover, while precipitation is rapidly infiltrated, the role of macropores is key—but requires the use of dual permeability functions to simulate appropriately. While these functions are available in the Hydrus 1D model, for example, they have a greater number of parameters that need to be characterized either through measurement or statistical procedures. There is scant information about the values of these parameters in the peat literature, so the use of these functions modeling of water flow in peat and peatlands has been rare.
Like simulating water flow, numerically representing solute transport in peat and peatlands is difficult. The considerations outlined below apply to one-, two-, and three-dimensional approaches, although their application has been primarily in one-dimensional simulations due to the difficulties associated with parameterization. A primary consideration when modeling solute transport in peat is its dual porosity nature, which creates conditions where a given solute can enter the immobile porosity and be removed from the advective flux. Thus, solute transport is commonly characterized by the mobile—immobile solute transport model (van Genuchten and Wagenet, 1989) to account for the immobile porosity within peat (𝜙mob and 𝜙im are described in Section 3.1, Saturated Zone Properties and Processes). However, there is some evidence that the transfer rate of solutes into immobile porosity can be sufficiently high in some peat that the mobile— immobile solute transport model simplifies to the advection—dispersion model (Simhayov et al., 2018; McCarter et al., 2019) or that peat can be represented by dual permeability models (Liu et al., 2017).
The advection dispersion solute transport model assumes that the mobile porosity is approximately the same as the total porosity, a key assumption that is not generally met when simulating solute transport in peat and peatlands. Conversely, dual permeability models are a group of models that simulate at least two different pore domains with dramatically different hydraulic and/or solute transport properties. For instance, in highly- degraded fen peat, macropores can accelerate the appearance of a solute pulse, but complete breakthrough of the solute plume occurs much later, often exhibiting a multi-modal flushing curve in breakthrough experiments (Liu et al., 2017).
Reactive solutes, such as cations or nutrients, require further parameterization beyond the hydrophysical parameters discussed above (i.e., characterization of biogeochemical parameters of the peat and the solute). Like most peat systems, these parameters are often not well characterized in the literature. Even assuming a conservative tracer such as chloride in peat can result in erroneous results. For instance, fitting of breakthrough curves using one-dimensional modeling has shown that anion adsorption (McCarter et al., 2018) and anion exclusion from narrow throat pores can operate in peat (McCarter et al., 2019), either retarding or accelerating solute breakthrough.
Furthermore, most reactive solutes are either simply modeled with first order decay/production coefficients or adsorption isotherms (often Langmuir, Freundlich, or Linear isotherms), depending on the specific processes in question. Simulating cation transport in peat requires representation of its cation exchange capacity, especially since this is strongly related to organic matter content (Gharedaghloo and Price, 2021). Similarly, first order decay/production coefficients are often used for most reactive solutes, regardless of the complexity of biogeochemical reactions being simulated. Numerical modeling of reactive solutes in peat is often limited by the inability to properly account for the complexity of organic matter composition (peat), mostly due to the lack of a detailed mechanistic understanding of the various biogeochemical processes and associated numerical expressions of such processes in peat.
In spite of the challenges, one -dimensional simulations of vertical flow and transport have been used to illustrate the nature of important hydrological processes in peat and peatlands. One-dimensional models have shown that the pore distribution in mosses controls evaporation from peatlands; certain Sphagnum hummock species better maintain moisture in spite of their elevated position (McCarter and Price, 2014); and, air entry pressure is positively correlated with carbon accumulation (Kettridge et al., 2016). However, simulations of capillary rise of water and solutes have shown that even adjacent and visually similar hummocks can have marked variability in hydraulic properties that strongly affect solute distribution (Balliston and Price, 2020).
While laboratory evaluation of hydraulic properties governing water retention and unsaturated hydraulic conductivity are useful for simulating flow conditions in a laboratory column, calibration of a one-dimensional model representing a field site can generate distinctly different parameter values (Elliott and Price, 2020). Simulation of these processes in multidimensional models is possible, but the difficulty in parameterization is compounded by the added spatial variability inherent in all peatland systems.