On the optimal level of complexity for the representation of groundwater-dependent wetland systems in land surface models

Wetland systems are among the largest stores of carbon on the planet, the most biologically diverse of all ecosystems, and dominant controls on the hydrologic cycle. However, their representation in land surface models (LSMs), which are the terrestrial lower boundary of Earth system models (ESMs) that inform climate actions, is limited. Here, we explore different possible parameterizations to represent wetland–groundwater–upland interactions with varying levels of system and computational complexity. We perform a series of numerical experiments informed by field observations from a particular type of wetland, called a fen, at the well-instrumented White Gull Creek in Saskatchewan, in the boreal region of North America. In this study, we focus on how modifying the modelling connection between the upland and the wetland affects the system's outcome. We demonstrate that the typical representation of groundwater-dependent wetlands in LSMs, which ignores interactions with groundwater and uplands, can be inadequate. We show that the optimal level of model complexity depends on the land cover, soil type, and the ultimate modelling purpose, being nowcasting and prediction, scenario analysis, or diagnostic learning.

1 Introduction

The Canadian boreal region covers about half of the land area of Canada. Around 85 % of all Canadian wetlands (∼ 1.3×10⁶ km²) are located in the boreal region. Wetlands are vital elements in landscapes, as they can mitigate the effect of floods, store carbon from the atmosphere, improve water quality, absorb pollutants, and provide a habitat for a wide range of endangered wildlife and plants (Mitsch et al., 2013). Types of wetlands are bogs, fens, swamps, marshes, and shallow water. Each type of wetlands differs in terms of hydrology, water level, morphology, vegetation, and biological aspects (Canada Committee on Ecological (Biophysical) Land Classification, 1988). Fens are wetlands that have accumulated peat of over 40 cm, hydrologically interact with the surrounding groundwater (GW) and surface water (SW), and maintain a water level at or above ground level for most of the year (Gingras et al., 2018). Fens constitute around 65 % of the peatland area within the boreal plain ecozone (Goodbrand, 2013). Fens critically depend on GW discharge fluxes to sustain their moisture and water levels. Understanding the lateral hydrological interactions between GW and SW in GW-dependent wetland/fen systems is crucial to improve their representation in land surface models (LSMs) (Rivera, 2014). Such improvements can directly improve the simulation of land's energy and water balance as well as different hydrological-cycle components like evaporation and streamflow (Blyth et al., 2021).

LSMs were originally proposed to estimate the vertical fluxes (energy and water) of the land surface, which is a necessary lower boundary condition for climate models (Manabe, 1969). Over the past decades, these models have been extensively modified to represent different processes such as soil moisture and vegetation dynamics. However, many recent studies have highlighted the deficiencies in the current LSMs and discussed the scientific motivation behind improving their process representations (Clark et al., 2015; Davison et al., 2016; Fan et al., 2019). Lateral water movement, GW dynamics, wetland hydrology, hillslope hydrology, and GW–SW interactions are examples of the elements that are either missing or need more realistic representation. The representation of these processes requires the inclusion of sub-grid heterogeneity in LSMs (Blyth et al., 2021). Capturing the complex interactions and heterogeneity of fine-scale landscape features would enable LSMs to provide more realistic and reliable simulations, facilitating a better understanding of land–atmosphere interactions and their implications for climate dynamics (Fan et al., 2019).

Significant advancements have recently been achieved with respect to representing both vertical and lateral water exchange of GW with SW in LSMs (de Graaf et al., 2017; Fan et al., 2013; Kollet and Maxwell, 2008; Maxwell and Miller, 2005; Maxwell et al., 2015); however, most of these efforts have focused on one-way process coupling (i.e., water moves in one direction from SW to GW as recharge), and any feedback from the GW system to the SW system has been neglected. Multiple studies have implemented two-way coupling between the GW aquifer and soil column within LSMs, such as the Community Land Model (CLM) and the Variable Infiltration Capacity (VIC) model, to enhance GW system simulation (Scheidegger et al., 2021; Zampieri et al., 2012; Zeng et al., 2018). Zampieri et al. (2012) introduced a simple parameterization in CLM v3.5 (Dai et al., 2003) that allowed for river–GW bidirectional flow. This parameterization showed an improvement in the soil moisture and surface temperature simulations. Zeng et al. (2018) developed a module for simulating lateral water exchange between grid cells and integrated it into CLM4.5. This module facilitates lateral exchange between grid cells but does not allow lateral exchange between GW and rivers. Scheidegger et al. (2021) incorporated a GW model into the VIC model that allows for a bidirectional water exchange between the soil and GW aquifer. Their study showed that including the 2D GW model significantly affected evapotranspiration (ET), runoff, and recharge of the different simulated grid cells (computational elements). Despite these recent efforts, a significant research gap remains with respect to effectively representing the two-way lateral exchange between GW and open waterbodies, including lakes and wetlands.

Typically, the coupling between different processes can be represented via one of three approaches: uncoupled, one-way coupled, or two-way coupled. The choice of the appropriate approach depends on the desired complexity level that can predict the variables of interest (Ogden, 2021). The literature generally demonstrates that coupling different models can improve the simulation of various hydrological-cycle components, including runoff, soil moisture, and water table fluctuations. However, such “full” coupling approaches can be both systematically and computationally complex, often rendering their applications impractical. In addition, constraints with respect to data availability, especially concerning subsurface processes, can limit the applicability of complex approaches. Therefore, in practice, more parsimonious approaches to accounting for SW–GW interactions may be “optimal” for specific modelling goals (Blyth et al., 2021; Ogden, 2021). What has remained elusive, however, is a thorough characterization of trade-offs between model complexity and adequacy to represent SW–GW processes for a particular landscape and variables of interest (Yalew et al., 2018). Another factor that affects the selection of optimal model complexity is the ultimate modelling purpose. Models can be used for (1) nowcasting and prediction, which focus on simulating and predicting the expected behaviour of the system of interest in the near future; (2) scenario analysis, wherein the model simulates the system under changing conditions over a longer term; and (3) diagnostic learning, which focuses on understanding how the system behaved during a historic period (Razavi et al., 2022). Ultimately, the selection of the optimal model complexity will depend on the specific modelling task, the available data, and the desired level of accuracy.

Here, to address this gap, we aim to characterize the optimal level of complexity required to represent the interactions between uplands, GW, and wetlands in different landscape configurations. To achieve this, we conduct a series of modelling experiments using four approaches, ranging from a full disconnect to full coupling. We particularly focus on a well-instrumented and well-studied fen system in the boreal region of North America. Our primary objective is to examine how the simulated upland and fen systems respond to changes in modelling complexity. Groundwater table (GWT) observations are utilized to derive a reasonable set of parameters representing the upland–GW–fen system. Additionally, using numerical experiments, we investigate the performance of a SW-fed wetland, providing a contrasting non-fen wetland scenario against the GW-dominated case. Based on insights gained in these experiments, we provide some recommendations on how to improve LSMs with respect to representing critical processes related to wetland types such as fens; these processes are often missing or poorly represented in the current generation of LSMs.

Figure 1 Detailed view of the study area. Panel (a) presents Canada's boreal region (gray shading). Panel (b) provides a focused view of the White Gull Creek basin (WGCB) area and the two flux towers used in the study (Old Jack Pine, OJP, and Fen); the shaded gray area represents the Pine Fen Creek Catchment (PFCC) that includes the Pine Fen (PF). Panel (c) is a focused view of the OJP area, the piezometer at OJP (POJP), and the Pine Fen (PF). Panel (d) gives a cross-section of OJP and the PF (the dashed black line is assumed to be the line of symmetry).

2 Study area and data

The study area is located within the White Gull Creek basin (WGCB), located north of Prince Albert, Saskatchewan (Barr et al., 2012), as shown in Fig. 1. The study area and transect are set using the Canadian Digital Elevation Model (CDEM; Government of Canada, 2022).

The upland is the area around the Old Jack Pine (OJP) flux tower (53.92^∘ N, 104.69^∘ W; Fig. 1), which is part of the Boreal Ecosystem Research and Monitoring Sites (BERMS). The upland transect ends at the boundaries of the WGCB; the OJP site is located roughly in the middle of the transect (Fig. 1d) and is co-located with the piezometer (POJP; Fig. 1c) that is used to calibrate and validate the model performance. The transect length of 3450 m is divided into the upland hillslope (which is 3300 m in length) and the 150 m wide fen (which is half the total width of the fen). At the OJP site, the dominant land cover is Jack pine (Pinus banksiana Lamb.); the soil has a sandy texture, poor nutrition, and high drainage; and the water table is around 5 m below the ground (Barr et al., 2012). The meteorological data at OJP are available each 30 min for 23 years from 1997 to 2018. The GWT observations are available at the POJP location (Fig. 1) from 2003 to 2018.

The lowland part of our transect is a fen known as Pine Fen (PF) which is situated in Pine Fen Creek catchment (PFCC), a tributary of White Gull Creek (39.9 % of PFCC land cover is peatland) (Goodbrand et al., 2019). The PF site is a peatland mosaic surrounded by forests of Jack pine and black spruce. The average peat thickness is 0.65 m with a maximum depth of about 2 m. Unfortunately, we do not have direct meteorological observations from PF; therefore, as a proxy, we use observations from the BERMS fen flux tower site (FEN hereafter). FEN is located just outside the WGCB boundary, about 8 km south of the basin (53.78^∘ N, 104.69^∘ W; Fig. 1b). The FEN and PF sites are similar in terms of the peat soils, vegetation, and topography; therefore, FEN is considered a reasonable proxy for the vertical fluxes at PF. At FEN, forcing data are recorded with a 30 min resolution from 2003 to 2018, and the observed ET rates are from 2004 to 2010 and from 2013 to 2019 at 30 min intervals.

Figure 2 Schematic of the different scenarios representing the connection between the upland, GW, and fen: (a) uncoupled upland–fen and no GW (V0), (b) uncoupled upland–GW–fen (V1), (c) chained model (V2), and (d) coupled model (V3). The vertical fluxes are rain on the ground (R_g), snowmelt (M), ET from the fen (E_f), and soil drainage as recharge into GW (R). The lateral fluxes are upland runoff (R_O) and GW discharge (Q_G).

3 Upland–GW–fen model description

3.1 Conceptual background

Our study is based on a real field site, using real field observations, as described in the previous section. However, this site is not a perfectly constrained hillslope–fen system, i.e., a hillslope with 1D horizontal flow and a no-flow boundary condition at the interfluve. Therefore, we use an abstracted hypothetical hillslope configuration to simulate the vertical and lateral flows of water within the atmosphere–upland–GW–fen system. This configuration is physically realistic, allowing us to test the implications of different hillslope–fen coupling mechanisms in a controlled manner (Fig. 2).

Our model has three distinct components: (1) the upland soil water balance, which generates GW recharge and runoff; (2) a lateral GW flow model beneath the upland, which may discharge water into the fen; and (3) a simple fen water balance model, which receives inputs from rainfall, snowmelt, and runoff and lateral GW discharge (in some cases), and loses water to evaporation and discharge into a stream channel.

The model is driven by site observations of precipitation and other meteorological variables. A typical LSM is used to simulate ET fluxes in the upland and fen, GW recharge, runoff fluxes from the upland, and snowmelt into the fen. In this case, runoff fluxes are much smaller compared with the GW recharge. The water table and GW discharge are simulated by a simple 1D unconfined aquifer model. The water level at the fen and its discharge flux into an adjacent stream are simulated by a simple fen water balance model. The connections between the three models of the upland, GW, and fen are configured in four ways (as shown in Fig. 2), resulting in four different collective models:

• V0.

  Uncoupled upland–fen model. The upland soil water balance and the fen are simulated as being independent of one another, and there is no GW model. Discharge from the upland in the form of surface runoff and soil drainage (baseflow) are combined with the fen discharge and routed into the river. This configuration is representative of many LSMs.

• V1.

  Uncoupled upland–GW–fen model. Soil drainage from the upland recharges the unconfined aquifer. Surface runoff from the upland, GW discharge, and fen discharge are combined and routed into the river. The upland and the fen are again completely independent of one another. This configuration is representative of an LSM that has a GW store with a one-way vertical connection with the soil column, such as the Community Land Model (CLM) and the Variable Infiltration Capacity (VIC) model (Clark et al., 2015).

• V2.

  Chained model. Soil drainage from the upland recharges the unconfined aquifer, and GW discharge contributes to storage into the fen. Discharge from the fen is routed into the river. GW is independent of the fen, but the fen depends on discharge from GW.

• V3.

  Coupled model. Soil drainage from the upland recharges the unconfined aquifer. GW discharge into the fen is determined based on the head gradient between the GW and fen, and two-way water exchange is considered between the GW and fen. Surface runoff from the upland also goes into the fen. Discharge from the fen is routed into the river. The GW and fen are mutually dependent on one another.

3.2 Upland soil water balance model

The upland soil water balance is simulated using the MESH/CLASS land surface model (Canadian Land-surface scheme; Pietroniro et al., 2007; Wheater et al., 2022) with a point-scale set-up using the data from the OJP site, which is referred to as MESH_OJP hereafter. The land cover of the grid cell is represented by one canopy type – evergreen needle leaf – as most of the vegetation in the OJP area consists of Jack pine trees. MESH_OJP is forced using seven metrological components (precipitation, short-wave radiation, long-wave radiation, wind speed, specific humidity, temperature, and atmospheric pressure), which were collected from the OJP flux tower every 30 min from 1997 to 2018. The soil depth is assumed to equal 4.1 m and is divided into three layers (i.e., the default CLASS configuration). We are mainly interested in two outputs from the model: (1) soil drainage (R) is used as either baseflow to the river (V0) or recharge to the GW (V1, V2, and V3) and (2) surface runoff (R_O) is used as SW input to the river (V0 and V1) or input to the fen (V2 and V3). The MESH/CLASS model uses an infiltration excess based on the Green–Ampt equation to calculate the soil infiltration, and the excess ponded water on the top of soil column is used to calculate the overland runoff.

3.3 Upland–GW model

The GW underneath the upland zone is represented as a 1D horizontal unconfined aquifer, bound by a no-flow boundary on the right (the GW divide) and the fen on the left (Fig. 2). The problem is governed by the 1D Boussinesq equation:

$$\begin{array}{}\text{(1)}& S_{\mathrm{y}}{\displaystyle \frac{\partial h}{\partial t}}={\displaystyle \frac{K}{\mathrm{2}}}{\displaystyle \frac{{\partial }^{\mathrm{2}}{h}^{\mathrm{2}}}{\partial x^{\mathrm{2}}}}+R,\end{array}$$

where S_y (unitless) is the specific yield, h (m) is the head in the aquifer (water table level), t (days) is time, K (m d⁻¹) is the lateral hydraulic conductivity, and R (m d⁻¹) is the recharge rate.

Equation 1 is solved numerically using a block-cantered finite-difference solution on a regularly spaced grid ($\mathrm{d}x=\mathrm{0.01}\cdot L)$, integrated in time using the method of lines with the SciPy ODE solver “odeint” (Virtanen et al., 2020), and the solution is coded up in Python. The lateral fluxes are calculated at cell boundaries using Darcy's law. The initial condition is assumed to be a uniform hydraulic head in the aquifer at the same water level as in the fen (h_f). The right-hand boundary condition is considered to be a GW divide and, thus, a no-flow condition ($q(x=\mathrm{0},\mathrm{0}\le t\le t_{\text{max}})=\mathrm{0})$. The flux on the left-hand boundary is determined using either a fixed head boundary (V1 and V2) or based on the head gradient between the GW and the fen (V3). The GW model is driven by the recharge fluxes (R), which are output from the upland soil water balance model and are assumed to be spatially uniform.

3.4 Fen water balance model

The fen is modelled as a simple lumped store, with the water balance equation:

$$\begin{array}{}\text{(2)}& n_{\mathrm{f}}{w}_{\mathrm{f}}{\displaystyle \frac{\partial h_{\mathrm{f}}}{\partial t}}=(R_{\mathrm{g}}+M-E_{\mathrm{f}})w_{\mathrm{f}}+R_{\mathrm{O}}L+Q_{\mathrm{G}}+Q_{\text{in}}-Q_{\text{out}},\end{array}$$

where n_f (unitless) is the porosity of the fen's material; w_f (m) is the width of the fen; h_f (m) is the water level of the fen; L (m) is the aquifer/hillslope length (in the x direction); $R_{\mathrm{g}},M,E_{\mathrm{f}}$, and R_O (m d⁻¹) are the rainfall, snowmelt, evaporation, and runoff from the upland, respectively; Q_G (m² d⁻¹) is the lateral GW inflow; and Q_out (m² d⁻¹) is the outflow from the fen, assuming a unit area of 1 m in the y direction.

For E_f, R_g, and M, these fluxes are simulated by a point-scale MESH model using the data from the FEN site (hereafter referred to as MESH_FEN), which reasonably represent the conditions at the PF site. The forcing data (2003–2018) from the FEN flux tower are used to drive the MESH_FEN model. To simulate the peatland in the MESH/CLASS model, the land cover is assumed to be grass, while the soil type is set to organic soil with three soil layers of changing properties – fibric, hemic, and sapric (Letts et al., 2000). We manually calibrate the minimum stomatal resistance (r_s,min = 100 s m⁻¹) using local observations of the driving meteorological variables. We compared the simulated E_f with both the observed fluxes and with potential evaporation calculated using the Penman–Monteith equation, and we found that all three are consistent, showing that the evaporation from the FEN is unstressed (i.e., not water limited). The values of Q_G are either zero (V0 and V1) or equal to the GW discharge at the right-hand boundary (V2 and V3).

The values of Q_out are generated when a spilling threshold, denoted as h_spill, is exceeded. Q_out is assumed to be a non-linear function of storage above h_spill as given by the following:

$$\begin{array}{}\text{(3)}& Q_{\text{out}}=w_{\mathrm{f}}{c}_{\text{spill}}(h_{\text{spill}}-h_{\mathrm{f}}{)}^n.\end{array}$$

Table 1 Monte Carlo analysis parameter ranges for the uncoupled upland model.

The value of h_spill is assumed to equal the elevation of the fixed head boundary at the uncoupled GW model (Fig. 2b, c). The values of c_spill and n (coefficient and exponent, respectively) are calibration parameters and are somewhat arbitrarily chosen within reasonable ranges. As there are no available outflow observations at the FEN site to allow for calibration, we adopt the recommended ranges provided by Razavi and Gupta (2019) for a fast reservoir with non-linear response (similar function/response to the fen storage). Razavi and Gupta (2019)'s values of k and α (coefficient and exponent, respectively) are utilized to define the c_spill and n values (c_spill=0.1 and n=1.5, respectively). While this approach lacks specific calibration for our site, it is deemed acceptable for achieving our study's objective, which is to compare different model configurations (V0, V1, V2, and V3) using the same parameters.

4 Model analysis and performance evaluation

4.1 Calibration strategy

The performance of the different collective models is evaluated using the GWT observations at the POJP location by calculating the root-mean-square error (RMSE). For the upland water balance model, we use the same calibrated parameters as the study of Nazarbakhsh et al. (2020), who used the CLASS model to assess the controls on ET in a seasonally frozen forest. For the other parameters, Monte Carlo simulations with 15 000 realizations (randomly generated parameters from a uniform distribution bounded by the feasible parameter space in Table 1) are used to run the uncoupled upland–GW (V1) model. The behavioural runs are identified as the realizations with RMSE < 0.08 m (threshold is chosen arbitrarily based on expert judgement and calculated for the period from 2003 to 2009) and are used to perform the uncertainty analysis of the GWT simulation. The parameter set with lowest RMSE is considered the calibrated parameter set and is used to validate the model. This parameter set is also used to run all the other collective models throughout the study.

4.2 Effect of different upland properties

We develop two hypothetical numerical experiments to explore the conditions under which different levels of model complexity may be necessary. Experiment 1 focuses on hillslope geometry by considering five different values of the hillslope length (L, which is the length of the upland). L values of 100, 300, 500, 1000, 2000, and 3000 m are considered, and all other parameters are the same as in the original model set-up. This experiment uses the chained (V2) and the coupled (V3) versions of the model. Experiment 2 focuses on soil properties by comparing the original (sandy soil) set-up with a fine-grained soil representative of mineral soil that can typically be found in the prairie area. This is achieved using an alternative configuration of the MESH model in which the parameters are exchanged with values to represent a grass land cover and a fine-grained soil texture (clayey soil), resulting in different amounts of runoff, infiltration, and soil drainage. In the upland algorithm, the values of K and S_y are set to 1.35 m d⁻¹ and 0.1, respectively, to characterize fine-grained soil. The simulated GWT and GW fluxes using the new MESH/CLASS run (representing fine-grained soil) are compared with the model results when using the original study set-up (upland algorithm forced with MESH_OJP).

5 Results and discussion

We assess the performance of the model in the upland (Sect. 5.1) and fen (Sect. 5.2) independently. Next, in Sect. 5.3, we assess the sensitivity of the simulated outflow from the integrated upland–GW–fen system, which corresponds to the outflow at the grid-cell scale in LSMs (i.e., the bulk system outflow) to the different model configurations. Lastly, in Sect. 5.4, we describe the results of the two numerical experiments exploring the effect of changed upland properties and different wetland functionality.

5.1 Model performance in the upland

In the upland, the performance of the vertical land surface fluxes has been explored by Nazarbakhsh et al. (2020). We are able to assess the upland model's performance with respect to reproducing observations of the water table elevation. As explained earlier, there is no GW in V0 of the model. In the uncoupled (V1) and chained (V2) models, the GW simulations are identical. In the coupled model (V3), the water table simulations may, in principle, differ from those in V1 and V2. Therefore, we compare V1 and V3 separately here.

Figure 3 Panel (a) shows the upland recharge rates into the GW aquifer that are generated using the MESH_OJP model and used to drive the upland component. Panel (b) presents a comparison between the simulated (V1) and observed changes in the GWT at the OJP site. The uncertainty results (GW discharge and GWT) are obtained from the behavioural realizations using Monte Carlo analysis.

5.1.1 Uncoupled model calibration and validation

The upland component in the uncoupled model is driven by the recharge values that are generated using the MESH_OJP model (Fig. 3a). Figure 3b shows a comparison between the simulated and observed GWT at the OJP site using the uncoupled upland model. The behavioural runs (69 runs with an RMSE < 0.08 m) from Monte Carlo analysis are used to generate the uncertainty bounds in the simulation of both GW discharge and GWT (Fig. 3). It can be seen from Fig. 3a that the recharge and the simulated GW discharge responded to each other consistently. Also, the changes in the GWT corresponded to both the recharge and the GW discharge peaks.

For the uncoupled upland model, the best RMSE value is equal to 0.058 m for calibration and 0.28 m for validation. The uncoupled upland model can simulate the GWT in the calibration period with a narrow uncertainty bound. In the validation period, the simulated GWT matched the observations until the spring of 2011, when a discrepancy is noticed, and the GWT is underestimated thereafter. The GWT underestimation is caused by low recharge rates from 2011 to the end of the simulation, which might be caused by either undercatch in the observed precipitation or problems with the MESH/CLASS model with respect to simulating the recharge rates in this period. The MESH/CLASS model problem could be because of the overestimation of ET rates at the upland site, which means that the MESH/CLASS model might need recalibration to a longer period of data. The reason for the overestimation of ET rates at OJP could be the missing representation of the wetland and GW system in the MESH model (see Figs. S1 and S2 in the Supplement). We should note that, although the model showed underestimation of the GWT magnitude (from 2011 to 2018), it captured the same pattern during the same period. Based on its performance, we believe that our model acceptably serves its purpose, as the main aim of this study is to perform numerical experiments to help us understand the upland–GW–fen system dynamics.

Figure 4 Comparison between the simulated GWT using both the uncoupled and coupled models with the observations.

5.1.2 Upland – uncoupled vs. coupled

Figure 4 shows a comparison between the simulated GWT using the uncoupled/chained (V1/V2) and coupled (V3) model versions and using the observations (note that the upland is still uncoupled, as in V1, in the V2 model configuration). The overall performance of the coupled version showed only a slight improvement (RMSE = 0.18 m) over the uncoupled version (RMSE = 0.22 m). A comparison between the simulated uncoupled and coupled systems shows that, in the period of record, considering the effect of the fen system does not affect the simulated GWT underneath the upland in the case of OJP and PF. However, the impact of this coupling might become more profound for other sites with different settings or the same site under different climate conditions.

Figure 5 Comparison between the fen's output (fen outflow and change in the fen water level) for the three modelling scenarios: uncoupled (V1), chained (V2), and coupled (V3).

5.2 Model performance in the fen

For the fen, we are not able to directly test our model performance due to a lack of data; instead, we explore the sensitivity of the fluxes to the change in the modelling configuration (interaction between the upland and fen). The fen's outflow and changes in water level are compared for the three model versions (V1, V2, and V3) in Fig. 5. In the uncoupled (V1) model, when there is no GW inflow to the fen, the estimated outflow and water level changes are unrealistic, as the outflow is almost zero and the water level kept decreasing from one year to another. On the other hand, the chained (V2) and coupled (V3) models had a reasonable simulation of the outflow and water level changes in the fen, although with differences from each other, particularly in terms of the flow rate. The overall trends in the flow rates of V2 and V3 look similar, but the higher-frequency features (e.g., daily flows) show different dynamics from time to time. The flow rates of V3 are affected by the two-way water exchange between the upland and fen, which are based on the variable fen water level, unlike the concept of a constant fen water level that is used in V2. This is apparent in 2004, wherein the V3 model generated negative flow rates (i.e., water flows from the fen into the upland). This comparison shows that the GW inflow from the upland into the fen cannot be ignored when simulating the fen; however, the chained modelling approach might be deemed adequate to capture the fen system dynamics. A coupled configuration is needed to study the short-term impacts and changes in the fen outflow.

Figure 6 Comparison between GW discharge into the fen and total outflow from the grid cell into the river at daily and annual scales for the V0–V3 modelling scenarios (Fig. 2).

5.3 Outflow of the integrated upland–fen system (grid-cell scale)

In this section, we investigated the total outflow from the integrated upland–GW–fen system as a whole unit under different modelling configurations, which are the possible approaches to represent such system in LSMs, to simulate the amount of the total flow that discharges into the river network (streamflow; see the dashed blue arrows in Fig. 2). In LSMs, a grid cell can contain multiple components, such as upland and fen. In this case, the total outflow of the grid cell is the combined outflow from both the upland and fen components after considering the interaction between them based on the used modelling configuration. This is done to assess the optimal level of model complexity that can simulate the streamflow adequately.

In the chained (V2) and coupled (V3) models, the simulated GW discharge into the fen from both models is almost the same at a daily and annual scale (Fig. 6). Also, the total outflow from the grid cell into river showed no significant difference in the two cases. In the case of the uncoupled upland–GW–fen model (V1), there is no GW discharge into the fen, but the daily and annual total outflow into the river is similar to that of V2 and V3. Therefore, at the grid-cell scale, the interaction between the upland and the fen had no effect on the total outflow from the grid cell nor discharge into the river for our model configuration.

In the case of V0, GW storage is not accounted for, and all of the soil drainage (recharge) is considered to be baseflow, which is discharged directly into the river network (the case in most of the current LSMs). Therefore, the simulated outflow (of V0) into the river is significantly different (with overall greater magnitudes) compared with the other three versions that account for GW storage (Fig. 6). This means that consideration of the GW dynamics underneath the upland is essential for a reasonable simulation of the total streamflow into the river network.

5.4 The effect of different upland properties on upland–fen interactions

Here, we run two additional numerical experiments to explore the optimal level of modelling complexity for different upland site properties. The experiments are hypothetical (no observations) and represent other possible sites' conditions.

Figure 7 Simulated upland GWT for both the chained and coupled model versions using different hillslope length (L) values.

Figure 8 Simulated upland GW fluxes for both the chained and coupled model versions using different hillslope length (L) values.

5.4.1 Experiment 1: different hillslope lengths

Figure 7 shows the simulated upland GWT using different upland hillslope lengths (width of the fen is constant) and compares the results in the case of chained (V2) and coupled (V3) upland configurations. In the case of horizontally large aquifers (L > 1000 m), as in our original study set-up, there is no significant difference in the simulated GWT and GW flux when using the two model configurations. In contrast, for small hillslopes (lengths between 100 and 500 m), the chained model is not able to reasonably capture the fluctuations in the upland GWT (GWT is almost constant), due to a fixed head boundary condition on the upland side. Also, the difference can be seen when comparing the simulated GW fluxes in the two model configurations (Fig. 8). For small hillslopes, the coupled model can capture the water amounts that move from the fen into the upland (negative flux values), which are considerable amounts that are frequently present throughout the year.

For large hillslopes, the GW size (extent) is significantly greater than the fen size (140 m); therefore, large amounts of water move from the upland area into the fen. As a result, the upland area controls the dynamics of the whole system. In such cases, the upland can be simulated independently with no consideration of interaction with the fen. In the case of small hillslopes, however, the fen is the dominant contributor to the system, as the water moves continually in both directions. In such cases, the coupling between the upland and fen is essential.

Figure 9 Comparison between two different upland site conditions and the simulated upland GWT in each case: (a) OJP with a coarse-grained soil texture and evergreen needleleaf canopy and (b) St. Denis with a fine-grained soil texture and grass land cover.

5.4.2 Experiment 2: different upland soil properties

Figure 9 illustrates a comparison between two upland cases: (1) coarse-grained soil (high permeability) and forest land cover (dense vegetation) with a long hillslope length, which is our original study set-up (Fig. 9a), and (2) fine-grained soil (low permeability) and grass land cover (Fig. 9b). In MESH_OJP, changing the vegetation type impacts ET rates from the upland and R_O, whereas altering the soil properties affects the infiltrated water and R rates. In this experiment, the first case represents the original model configuration at OJP, whereas the second case is designed to maximize surface runoff from the upland into the wetland. This provides a contrasting scenario with a different wetland type (other than a fen), in contrast to the GW-dependent wetland case. The parameters used to configure the MESH model for the second case are based on a model configuration presented in Ireson and Sanchez-Rodriguez (2021) and are similar to the parameters and properties of the St. Denis site.

In the case of coarse-grained soil (Fig. 9a), the main contribution to the upland GW system is high recharge rates because of the high infiltrability of the soil (coarse-grained/sandy) and the relatively small amounts of surface runoff to the fen. Thus, in this case, the dominant component of the upland–fen system is the GW fluxes from the upland into the fen (through subsurface water movement). Water arrives as precipitation on the upland, infiltrates into the soil and recharges the aquifer, and finally moves laterally in the aquifer to discharge into the fen. To represent these system dynamics, the chained (V2) modelling approach for the upland system seems adequate to simulate the GW system (GWT and GW discharge) of the upland (this might not be applicable to relatively small hillslopes, as demonstrated in Experiment 1 in Sect. 5.4.1), as the difference between the results when using the chained (V2) and coupled (V3) models is relatively small (Fig. 9a). However, the coupled (V3) model can simulate the dynamics of the daily GW flows due to the frequent change in the fen water level. Accordingly, the flow direction is reversed to be from the fen into the upland (negative flow values) from 2004 to 2005 (Fig. 9a).

The system dynamics are different when the soil has a fine-grained texture and the vegetation has a low density (Fig. 9b). In this case, large amounts of surface runoff move directly from the hillslope into the wetland, whereas very limited amounts of water may infiltrate into the upland aquifer. Figure 9b shows that the chained model is not able to simulate any of the GW dynamics underneath the upland, whereas the coupled model captured the system dynamics (GW fluxes in Fig. 9b). In this case, the wetland is mainly fed by the surface runoff fluxes, and the water moves laterally into the upland (especially during snowmelt season). Hence, the main flow path is from the wetland into the upland GW. Therefore, the upland GW dynamics are dominated by the subsurface water fluxes coming from the neighbouring wetlands. It is obvious that considering only the chained approach (one-way exchange between the upland and wetland) in the case of fine-grained upland soil cannot reasonably capture the real dynamics of the system. However, the full coupling between the upland and the wetland (two-way water exchange) allows the model to represent the actual dynamics of the upland aquifer underneath fine-grained soil layers.

6 Conclusions

The insights gained from applying alternative model configurations to the upland–fen system in this study are as follows:

1. We were able to reasonably simulate the GW dynamics underneath the upland using the 1D Boussinesq equation. There are no significant differences between the coupled and uncoupled modelling approaches with respect to simulating the upland water table elevation, as the dominant flow direction is from the upland to the fen.

2. To accurately simulate the water level in the fen, it is important not to ignore the GW input from the upland. However, there is no significant difference between the chained (one-way interaction) and coupled (two-way interaction) approaches in terms of the simulated fen water level and outflow.

3. Including upland–fen interactions had no significant impact on the discharge into the river network. However, the inclusion of GW storage had a major impact on the timing and magnitude of river discharge.

We found that a coupled fen–upland modelling approach becomes essential when the size of the fen is large relative to the upland area. This is due to substantial bidirectional water exchanges between the fen and the upland GW at various times of the year. The coupled modelling approach is also more likely to be necessary when simulating uplands with fine-grained soils. In this case, the wetland receives more surface runoff and less GW input, resulting in a significant loss of water to the GW system.

In general, if the main objective of the model is to simulate streamflow, the coupling (two-way interaction) between the upland and fen (or GW-dependent wetland types) can likely be ignored. However, GW dynamics must be represented in LSMs, as they significantly affect the total outflow (streamflow) from the entire system. Conversely, if the simulation of the storage and fluxes within fen/wetlands are of interest, the chained modelling approach represents the least complex level necessary to account for the contributions of the surrounding upland and GW systems.

This study provides insights into the necessary model complexity for simulating an upland–GW–fen system within LSMs. The outcomes of this study can aid in improving process representation in LSMs and guiding current and future hydrological modelling practices in fen-dominated areas. However, other types of wetlands need more investigation using data/observations including different site conditions. This can result in more accurate simulation of the water cycle in those regions, contributing to enhanced water resource management and allocation as well as improving LSMs' ability to predict the effects of future climate change on wetlands.