Modeling water flow and volumetric water content in a degraded peat comparing unimodal with bimodal porosity and flux with pressure head boundary condition

Degraded peatlands release large amounts of greenhouse gases. The development of effective mitigation and management measures requires an understanding of relevant site‐specific biogeochemical and hydraulic processes. However, the simulation of water fluxes and vadose zone state variables of degrading peatlands relies on proper process description, parameterization of hydraulic functions, and representation of the boundary conditions. The objective of this study was to analyze the effects of unimodal versus bimodal soil hydraulic functions and pressure head versus flux‐type lower boundary conditions (LBCs) on the calculated hydraulic characteristics of a degraded peat profile. HYDRUS‐1D was used to study the hydraulic flow dynamics parameterized with data from a weighable groundwater lysimeter for the period from May 1 to December 31, 2019. Simulations comparing uni‐ and bimodal hydraulic functions showed only minor differences. Simulations of soil water pressure at a depth of 30 cm using a flux‐type LBC (RMSE: 27 cm, where RMSE is root mean square error) performed better than simulations using a pressure head LBC (RMSE: 48 cm). The pressure head LBC performed better at simulating volumetric water contents in 30‐cm depth than the flux LBC variant (RMSE: 0.05 vs. 0.09 cm3 cm−3). For specific site conditions with a shallow, fluctuating groundwater table and temporary air entrapment, the choice of LBC was important for a more accurate simulation of soil water fluxes and volumetric water content.

Vadose Zone Journal gases (GHGs), potentially transforming these peatlands from carbon sinks into carbon sources (Tiemeyer et al., 2016).Therefore, the management of peatland drainage is closely linked to GHG emissions and has significant environmental implications.Recognizing the importance of addressing GHG emissions, especially from drained peatlands, the German climate protection strategy (BMUV, 2021) emphasizes the need to reduce emissions from German peatlands by raising the water table (WT).Achieving this goal requires effective water management in the affected areas.As a result, understanding the hydraulic processes that influence the water balance in degraded peat soils and impact GHG emissions is crucial.This knowledge is essential for the development of adaptation measures and management concepts aimed at mitigating GHG emissions from peatlands and restoring their ecological functions.
The hydraulic processes within degraded peatlands are known to be highly complex (Rezanezhad et al., 2016).As a result, there has been a growing interest in the field of hydraulic modeling to better understand these intricate processes.Numerical vadose zone models, such as HYDRUS-1D (Šimůnek et al., 2013), are commonly used for this purpose.These models rely on the Richards equation and require specific parameters related to soil water retention and soil hydraulic conductivity curves.These curves are often described using the van Genuchten-Mualem (vGM) functions (Mualem, 1976;van Genuchten, 1980).To accurately simulate the behavior of water in peatlands, numerical flow models need to account for soil hydraulic parameters(SHPs), including the saturated hydraulic conductivity (K s ), which are essential for describing soil water retention and hydraulic conductivity curves.Determining the effective SHPs typically involves inverse modeling, a process that derives parameter values that best match observations and model simulations (Wöhling et al., 2008).However, standard numerical models based on the Richards equation often do not consider air entrapment and preferential flow in macropores and local nonequilibrium conditions (Dettmann et al., 2014;Ma et al., 2010;Šimůnek et al., 2003).These macropores, which can be as large as 1-3 mm in diameter, play a significant role in determining soil hydraulic properties (Jarvis, 2007), including those of peat soils (Rezanezhad et al., 2016;Wallor, Herrmann et al., 2018).The influence of near-surface macropores on the hydraulic properties of peat soils was first reported in Baird (1997), and similar observations were made in highly decomposed peat soils (Schindler et al., 2003;Wallor, Herrmann et al., 2018).Macroporous soils have often been characterized macroscopically with dual-porosity models that assume two pore domains (i.e., macropore and matrix) with distinct soil hydraulic characteristics (Durner, 1994;Gerke & van Genuchten, 1993;Šimůnek et al., 2003;Vogel et al., 2010).Assuming local equilibrium conditions in the pressure head between the macropore and matrix pore systems, the flow

Core Ideas
• Understanding degraded peatland hydraulics based on groundwater lysimeter data.• Simulation of the peat soil water balance as affected by water table dynamics.• Comparing unimodal with bimodal soil hydraulic models in water balance simulation.• Analyzing the choice of the lower boundary conditions by scenario simulations.• Discussion on missing processes not accounted for to improve simulation results.
can be described using a single-porosity model, where the soil hydraulic properties of both the macro-and micropore systems are represented with a bimodal function (Durner, 1994).This bimodal-type model effectively combines the two soil water retention curves of dual-porosity soil into a single double-S-shaped retention curve (Filipović et al., 2020;Gerke, 2006).
Only a few studies have employed bimodal soil hydraulic models to analyze peat soil data, yielding improved hydraulic modeling (Dettmann et al., 2014;Dimitrov et al., 2010;Jung et al., 2017;Weber et al., 2017).However, these studies examined various peatland types at different degradation stages, impacting hydraulic properties (e.g., Wang et al., 2021).The suitability of bimodal soil water retention models for describing hydraulic processes in degraded peat soils remains unclear.
Studies involving in situ and lysimeter-scale measurements in degraded or aerated peat soils remain limited (Durner et al., 2008;Schwärzel & Bohl, 2003;Schwärzel et al., 2002;Schwärzel, Šimůnek, Stoffregen et al., 2006).Additionally, only one study has applied soil hydraulic parameters in a soil water balance simulation model (Wallor, Rosskopf et al., 2018).Consequently, it remains uncertain whether high-resolution measurement and modeling of transient field processes in heterogeneously degraded peatland sites can be achieved for use as a foundation in geochemical and reactive transport modeling.In terms of determining soil hydraulic parameters from data obtained in field lysimeter experiments, Durner et al. (2008) explored the simultaneous estimation of hydraulic properties for multiple soil layers through inverse modeling, albeit in a theoretical context.
Numerical models based on Richards' equation necessitate well-defined boundary conditions (BCs).In studies focused on shallow groundwater sites, the most commonly utilized and measured lower BC is the soil water pressure head, typically represented by the water table depth (WTD) at such locations (e.g., Akhtar et al., 2013;Naylor et al., 2016).Frequently, modeling efforts are coupled with inverse parameter estimation of soil hydraulic properties (SHPs) (e.g., Forkutsa et al., 2009;Karimov et al., 2014;Kelleners et al., 2005;Schwärzel, Šimůnek, van Genuchten et al., 2006).Various other conditions have been applied at the lower BCs, including (i) no flow-this condition is especially employed in investigations related to actual evapotranspiration (ET a ) or air entrapment (Dettmann et al., 2014;Schwärzel, Šimůnek, Stoffregen et al., 2006;Zhang et al., 2011) and (ii) free drainage-this condition is used to study nutrient leaching (e.g., Bian et al., 2021;Jha et al., 2017).However, we did not come across any shallow groundwater modeling studies where variable flux conditions were applied based on inflow and outflow measurements at the lower boundary.Collecting such data requires significant effort, particularly when utilizing a weighable groundwater lysimeter (Dietrich et al., 2016).Nevertheless, it's worth noting that the choice of a lower BC can substantially impact the parameterization of SHPs, although this issue was not encountered in studies involving the inverse estimation of SHPs for shallow groundwater table sites.
Sound hydraulic modeling relies significantly on the availability of high-resolution and accurate data (Vereecken et al., 2010).Field measurements conducted under transient conditions are often considered more representative than laboratory experiments, which typically operate on a smaller spatial scale (Durner et al., 2008;Kelleners et al., 2005;Schwärzel, Šimůnek, Stoffregen et al., 2006;Vrugt et al., 2008).In field settings, the measurement of various state variables and knowledge of water fluxes at boundaries are both considered desirable (Abbasi et al., 2003).These boundary fluxes can be accurately provided by a lysimeter, making it a valuable tool (Durner et al., 2008;Schrader et al., 2013).Therefore, it is assumed that soil hydraulic modeling, based on lysimeter data characterized by high spatial and temporal resolutions, can contribute significantly to a more detailed exploration of hydraulic processes in degraded peatland soils.
The primary objective of this study was to provide a quantitative description of hydraulic processes within degraded peat soils characterized by a shallow, fluctuating WT.This description would serve as a foundation for subsequent hydrogeochemical process analyses of peatlands.To achieve this goal, our study had two main aims: (1) to perform inverse modeling for the estimation of hydraulic properties, including uni-and bimodal water retention and hydraulic conductivity functions.This estimation was carried out for both flux and pressure head lower boundary conditions (LBCs); (2) to assess and compare the impact of these different hydraulic functions and lower BCs on the simulated water balance.This comparison was conducted by analyzing the model's performance against observed time series data of various state variables and fluxes measured within the lysimeter.
The software package HYDRUS-1D was utilized for the key purposes in this study: parameter fitting to estimate hydraulic properties and simulation of water balance components.These components encompassed fluxes at the lysimeter bottom and state variables such as volumetric water content, soil water pressure head, and WTD.

Experimental site and measurements
The study site is situated within the Spreewald wetland, a region located to the southeast of Berlin, Germany, with coordinates at 51˚52′N latitude and 14˚02′E longitude.This area consists of extensive grasslands that have been artificially drained.The peatland soils in this region have been subjected to degradation and are characterized by a shallow groundwater table, which experiences seasonal fluctuations ranging from 10 to 125 cm below the soil surface.The soil type at this site is classified as Mollic Gleysol according to FAO-ISRIC classification (FAO-ISRIC, 2006).It consists of three distinct soil horizons, as detailed in Table 1.To conduct this study, a weighable groundwater lysimeter was employed.This cylindrical lysimeter had a surface area of 1 m 2 and a volume of 2 m 3 .The lysimeter contained an undisturbed soil monolith of Mollic Gleysol and naturally occurring vegetation that mirrors the native flora of the site.The lysimeter's groundwater table was adjusted to match the depth of the surrounding WT using a system of chambers and pumps, as illustrated in Figure 1.A reference water-level gauge positioned at the wetland site continuously measured the representative WT for the entire site.The WT in the lysimeter was controlled and synchronized with the representative WT of the site through a reservoir chamber located adjacent to the lysimeter.Water could be automatically transferred into or out of the reservoir chamber from a nearby well (about 2-3 m apart), maintaining equilibrium with the lysimeter.This setup allowed for precise adjustments of the lysimeter's WT to mimic the surrounding conditions.For a more comprehensive understanding of the lysimeter station's methodology and technology, readers are referred to Dietrich et al. (2016).The Spreewald region in northeast Germany has experienced alluvial sedimentation processes, resulting in the formation of extensive layers comprising clay and organic matter that are interspersed among typical sandy layers.These particular layers are recognized for their relatively low hydraulic conductivities (Riek et al., 2014).This low hydraulic conductivity is characteristic for the second soil layer of the soil monoliths investigated in this study (Table 1).
The lysimeter station was equipped with a weather station that provided basic meteorological data in hourly resolution (i.e., net radiation, air temperature, relative humidity, and wind speed).To compute the reference potential evapotranspiration (ET p ) for a short, well-watered grass, the methodology T A B L E 1 Basic properties of the soil in the lysimeter as reported in Dietrich and Steidl (2021).Note: The particle size classes are sand (2-0.063mm), silt (0.063-0.002 mm), and clay (<0.002 mm) and C org is organic carbon content; the soil texture was classified according to the German KA5 system (Eckelmann et al., 2005).

Soil layer
F I G U R E 1 Schematic diagram of the groundwater lysimeter station, adapted from Dietrich et al. (2016).WTD is water table depth, ET a is actual evapotranspiration, P is precipitation, R in and R out are in and output fluxes at the lysimeter bottom, ΔS is water storage change of the soil monolith, θ is volumetric water content, here measured with a profile soil moisture probe at eight depths, and h is soil water pressure head measured here at three depths.
outlined in the FAO-56 irrigation and drainage paper by Allen et al. (1998) was employed.This calculation utilized hourly meteorological data collected from May 1 to December 31, 2019, (Table 2) for the purpose of this modeling study.

Numerical model
The one-dimensional vertical flow in a variably saturated medium was described using a modified form of the Richards equation for one-dimensional vertical flow and numerically solved using software HYDRUS-1D, version 4.17 (Šimůnek et al., 2013).
where θ is the volumetric water content (cm 3 cm −3 ), t is the time (h), h is the soil water pressure head (cm), z is the soil depth defined positive upward (cm), K(h) is the unsaturated hydraulic conductivity function (cm h −1 ), and sink is the sink term for root water uptake (cm 3 cm −3 h −1 ).
The water retention and hydraulic conductivity functions were described using two approaches.The first approach was the frequently used vGM (Mualem, 1976;van Genuchten, 1980).
where θ r and θ s are the residual and saturated water contents (cm 3 cm −3 ), respectively; α (cm −1 ), m (-), and n (-) are empirical vGM model parameters; K s is the saturated hydraulic conductivity (cm h −1 ); S e is the effective saturation; and l is a pore connectivity parameter.
The second approach was the bimodal Durner-vGM model, which was applied to the upper two soil layers because of the existence of macroporosity and dual porosity in the degraded fens, (e.g., Dettmann et al., 2014;Wallor, Herrmann et al., 2018).In this model, different pore size ranges are each represented by a vGM retention function that overlap in a single curve in which they are weighted with the factor ω (Durner, 1994).In HYDRUS-1D, two pore size distributions with predominantly larger (macropores) and smaller (micropores) mean pore sizes (as indicated by the inverse of the vGM-α-values) can be considered.A multimodal unsaturated hydraulic conductivity function, which combines the multimodal water retention function (Durner, 1994) with the conductivity representation (Mualem, 1976;van Genuchten, 1980), is given by Priesack and Durner (2006) as follows: where the sum of the weighting factors ω 1 to ω k is equal to 1, and k = 2 for a bimodal model.

Model setup
The observed ET a obtained from the field lysimeter data was split into actual evaporation (E a ) and actual transpiration (T a ) according to Ritchie (1972) and used in HYDRUS-1D as follows: where r Extinct is the canopy extinction coefficient, which was set to 0.75, and LAI is the leaf area index, which was linearly interpolated in time between regular LAI measurements.This procedure resulted in T a values that were, on average, larger than the E a values, which was deemed realistic for a year-round vegetated grassland.HYDRUS-1D was designed to use E p and T p values to calculate ET a .However, here, we assumed splitting observed ET a into E a and T a using LAI to be more accurate than using E a and T a values derived from T A B L E 3 Model parameters and initial and boundary conditions used in the HYDRUS-1D simulations.

Initial conditions
Pressure heads h in three depths measured at the initial time were extrapolated linearly throughout the depth of the soil profile.

Boundary conditions
Upper boundary: Atmospheric boundary conditions, input as P and E a , as well as T a Lower boundary: Either flux (R in ) or soil water pressure head (h).

Water stress function
The Feddes parameters were adapted to wetland plants-assuming no water stress for saturated conditions.P0 = 0, P0pt = 0, P2H = −200 cm, P3 = −800 cm.The critical stress index for water uptake was set to 1.

Root distribution and density
Roots were visible down to a depth of 30 cm from visual inspection.In HYDRUS-1D, the root density was assumed to be at a maximum of 100% in 0-to 30-cm depth, and it was linearly decreasing from 100% in 30 cm down to 0% in 40-cm depth.
Note: The parameters for the Feddes model are P0, which is the value of the pressure head below which root water extraction begins, P0pt is the pressure head below which the maximum root extraction begins, P2H is the value of the pressure head below which roots can no longer extract at maximum rate, and P3 is the value of pressure head at which root extraction ceases.
Abbreviations: E a , actual evaporation; h, soil water pressure head; P is precipitation; R in , water flux at the bottom boundary (+ inflow, − outflow); T a is actual transpiration.
F I G U R E 2 Illustration of steps in parameterization of model variables; θ is volumetric water content, h is soil water pressure head, HYDRUS-1D is the modeling software (Šimůnek et al., 2013), WTD is water table depth, RETC is a curve fitting software (van Genuchten et al., 1991), R in is water fluxes at the lower boundary (+ inflow, − outflow), ROSETTA is a pedotransfer function software (Schaap et al., 2001).Ks, saturated hydraulic conductivity; SHP, soil hydraulic property.
meteorologically based ET p values.Since HYDRUS-1D usually applies E a and T a values as upper boundary condition, we assumed here a relatively small soil water pressure head threshold value of −1.5 × 10 4 cm to ensure a negligibly small reduction of our E a values.The critical stress index for the water uptake was set to 1 (Table 3).The parameterization of soil hydraulic properties (SHPs) and K s values involved multiple steps, as outlined in Figure 2.
Initially, preliminary SHP values were generated for each soil layer utilizing the pedotransfer function software ROSETTA (Schaap et al., 2001).These preliminary estimates were calculated based on the particle size distribution and bulk density of the degraded peat soil.It is important to note that these initial retention curves were used solely as starting approximations in the RETC software (van Genuchten et al., 1991).The RETC software was employed to refine the retention curves by using measured pairs of water content versus soil water pressure head data for soil layers 1 and 3.This calibration process used data from May 1 to August 25, 2019.This iterative approach was essential to adjust the SHP values to represent the soil's hydraulic behavior under field conditions.The initial values of the SHPs obtained with RETC were used for the inverse optimization of the hydraulic parameters using HYDRUS-1D, in which the objective function was minimized using the Levenberg-Marquardt nonlinear minimization method (Marquardt, 1963;Šimůnek et al., 2013).Table 4 lists the data used in the objective function.Data were weighted using the mean ratio option in HYDRUS-1D.
Because the saturated water content θ s was measured, this variable was not optimized during the inverse modeling exercise.Residual water content θ r , and the parameters α, n, and K s were optimized, as they influence the shapes of the retention and hydraulic conductivity curves.
The inversely optimized parameters for the bimodal soil water retention curve and soil hydraulic conductivity curve (Durner, 1994) included the parameters ω2, α2, and n2 of the second pore system.These parameters were optimized for layers 1 and 2 to represent the dual-porosity features of the degraded peat.A ratio of macropore domain, ω2, of approximately 10% was assumed according to observations (e.g., Liu & Lennartz, 2015).As the second pore domain in the range of larger pore sizes and macropores (Šimůnek et al., 2013) should be represented by a retention curve with larger α2 and n2 values than those of the soil matrix, the initial values for the inverse optimization of α2 and n2 were gradually increased until convergence.Sandy layer 3 was assumed to have unimodal porosity owing to its homogeneously course-textured appearance.
Three model variants differing between the applied LBC and retention models were considered (Table 4).The first model variant was characterized by a flux LBC and a unimodal hydraulic function.The second model variant also uses a flux LBC but with a bimodal (dual-porosity) hydraulic function.The third model variant uses a pressure head as the LBC and a bimodal (dual-porosity) hydraulic function.All the model variants had the same upper boundary conditions.The LBC was used as an input to the inverse model in HYDRUS-1D.The retention models applied in HYDRUS-1D were either a single-porosity vGM model (Mualem, 1976;van Genuchten, 1980) or bimodal Durner-type dual-porosity vGM model (Durner, 1994).
The simulations encompassed the period from May 1 to December 31, 2019, with the data from May 1 to August 25, 2019, serving as the calibration period (Figure 3).Notably, this period concluded just before a substantial rainfall event.
It is worth mentioning that the validation period included the occurrence of the largest rainfall event of the year, involving 74 mm of rain within a 3-h timeframe.Throughout the calibration period, soil hydraulic parameters (SHPs) were determined inversely and subsequently applied to the entire simulation period.Model performance during the validation period, spanning from August 26 to December 31, 2019, was assessed using several metrics, including the root mean square error (RMSE), the Nash-Sutcliffe model efficiency coefficient (NSE) as proposed by Nash and Sutcliffe (1970), and the Wilmott index of agreement (IoA) based on Willmott et al. (2012).These metrics were employed to compare and evaluate the agreement between simulated and observed values.where O i and S i are the observed and simulated values, respectively, Ō is the average of the observed values, and n is the number of observed/simulated values.

Lysimeter and weather data
The model input (Table 2) was based on the measured lysimeter data (Figure 3).Throughout the study period, the lysimeter received a cumulative precipitation (P) of 540 mm and provided a total ET a of 640 mm.The largest rainfall event occurred on August 25, 2019 (74 mm within 3 h).The difference between E a and P (Figure 3b) represents the net amount of water that passes through the soil surface, and T a represents the plant root water uptake.The difference between the upper and lower boundary water fluxes is the water extracted from the soil water storage and is equal to the water storage change ΔS.The R in displayed an overall inflow of 25 cm from May to the largest rainfall event at the end of August.After the large rainfall event that occurred on August 25, 2019, the bottom flux, R in , accumulated to an inflow of 30 cm.At the end of the year, the sum of the inflow at the lower boundary (30 cm) was larger than the cumulative difference between ET a and P (27 cm) at the upper boundary.This resulted in an increase in water storage by 3 cm throughout the observation period (Figure 3c).Corresponding to the change in water storage, the WTD decreased from −33 cm down to −127 cm between May 1 and August 25, 2019 (Figure 3d).From the beginning of the time period until just before the largest rainfall event, storage decreased by 12 cm (Figure 3c).During the largest rainfall event, storage increased by 7 cm.The lowest pressure head value measured at a depth of 30 cm was −368 cm before sensor failure at −278 cm on July 16 (Figure 3e).While θ 30 and θ 40 in layer 1 also decreased strongly from 0.65 and 0.54 cm 3 cm −3 down to values of 0.27 and 0.30 cm 3 cm −3 , respectively, θ 50 in layer 2 and θ 60 in layer 3 revealed smaller relative changes (0.45-0.29 cm 3 cm −3 and 0.41-0.22cm 3 cm −3 , respectively, Figure 3f) despite the WTD at −127 cm.The drying out of soil layer 1 was interrupted by the largest rainfall event, during which the measured WTD increased by several decimeters to −49 cm within 1 h and to less than −10 cm below soil surface within 6 h (Figure 3d).After the rainfall event, the values of WTD, h, and θ in soil layer 1 receded again.Around December 1, 2019, the values of WTD, h, and θ increased to values similar to or larger than the values measured at the beginning of the time series in May.The values of θ 60 and θ 90 remained between 0.2 and 0.3 cm 3 cm −3 .T A B L E 6 Soil hydraulic properties estimated from inverse optimization in HYDRUS-1D.

Retention curves of the model variants
The parameter values for the van Genuchten model determined in the first step were obtained using ROSETTA and RETC (Table 5).
The soil hydraulic properties and the corresponding retention and unsaturated hydraulic conductivity curves for each variant and soil layer were obtained by inverse optimization (Table 6, Figure 4).The retention curves all display the same θ s value, as this value was set and not optimized.The curves differ in the value of parameter θ r , in steepness, which is mainly influenced by the van Genuchten parameter n, and in the location of the inflection point at the saturated end of the retention curve.The dual-porosity curves for the two upper horizons of variants 2 and 3 and the overlay with the measured data showed that the retention curves fit the data (Figure 4, Table 6).

Simulated time series
During the period from May to August 2019, it was observed that variant 3 of the model consistently underestimated the inflow at the lower boundary.By the end of this study period, a significant difference of 11 cm was noted, representing an underestimation of approximately 43% when compared to the measured inflow (Figure 5a).Interestingly, variant 3 exhibited a peculiar behavior during the largest rainfall event in August 2019 and the second-largest event in June 2019.It simulated inflow at the lower boundary, which contradicted the actual field lysimeter observations.Following the substantial rainfall event, the model continued to underestimate the overall inflow.It predicted 4 cm of outflow, while the measured inflow amounted to 5 cm between August 25 and the end of the time series on December 31, 2019.The simulation of WTD was conducted using model variants 1 and 2, as it served as an input for variant 3.Both variants 1 and 2 produced a good fit up until the occurrence of the largest rainfall event on August 25, 2019, as indicated by the measure of fit values (variant 1: RMSE 11.88 cm, NSE 0.83; variant 2: RMSE 11.74 cm, NSE 0.84, Table 7).During the calibration period, there was a maximum underestimation of the WTD by 44 cm, which occurred on July 21, 2019.However, the most substantial deviation was observed during the significant rainfall event on August 25, 2019.During and after this event, the simulated WTD deviated significantly from the measured WTD by up to 96 cm, as depicted in Figure 5b.The measured WTD increased from −123 to 8 cm below the soil surface, a rise of over 115 cm within 7 h.In contrast, the simulated WTD increased from −103 to −85 cm after 24 h, F I G U R E 4 Soil water retention curves and unsaturated hydraulic conductivity curves, identified by inverse optimization for the soil layers 1, 2, and 3 for the three model variants 1, 2, and 3, and measurements Obs30 and Obs90 in 30-and 90-cm depth in the calibration period.θ, volumetric water content; h, soil water pressure head; K, hydraulic conductivity.

T A B L E 7
Goodness-of-fit measures for simulated versus observed water table depth (WTD) of variants 1 and 2, cumulative fluxes at lower boundary (R in ) of variant 3, and volumetric water contents (θ) and soil water pressure head (h) in 30-and 60-cm depth, in the calibration and validation period.representing an increase of only 18 cm.In the later part of the validation period, the simulated WTD continued to underestimate the measured values by 45 cm on September 25, 2019, until another rainfall event occurred on October 5, 2019.

Model
The cumulative water storage change, as shown in Figure 5c, was slightly underestimated in variants 1 and 2, with an underestimation of 1 cm until the end of the calibration period in August and 2 cm by the end of the validation period.This underestimation is attributed to a reduction in the input ET a by HYDRUS-1D, as it only accepts the input of ET p and then reduces it to ET a .This reduction in ET a was observed due to dry conditions during the summer.In contrast, variant 3 consistently underestimated the change in water storage during most of the calibration period, with a difference of up to 3.5 cm at the end of the calibration period.However, after the rainfall event on August 25, variant 3 overestimated the measured water storage by up to 3.6 cm.The degree of overestimation decreased during the validation period, reaching an underestimation of 1.6 cm by the end of the validation period and the year.

Model outputs-Simulated vertical water fluxes
The HYDRUS-1D model variants enabled the evaluation of vertical water fluxes and volumetric water contents θ throughout the soil profile at every modeling time step.We selected the values for three time steps under two different conditions.The first three time steps show a situation without P and with water extraction using a high ET a .Water flowed mainly from the bottom to the top (Figure 6).The second three time steps show the situation after a rainfall event, with infiltration and mainly downward flux (Figure 7).
The upward fluxes throughout the soil profile during noon in early summer (high ET a and no P, June 6, 2012, at 12:00) are defined by the input E a values of the model at each time step with values between 0 and 0.02 cm h −1 at the soil surface (Figure 6a-c).Variants 1 and 2 show an upward flux throughout the entire soil profile with a near-uniform value of 0.018 cm h −1 in layer 3 and values between 0.05 and 0.06 cm h −1 in layer 1.In contrast, variant 3 displays a downward flux of −0.01 cm h −1 in layer 3, due to decreasing WTD.Note that the WTD is a BC for variant 3.In soil layers 1 and 2, variant 3 displays a smaller upward flux of 0.01 cm h −1 compared to model variants 1 and 2. The depth of the largest upward flux in all variants corresponded to the deepest simulated rooting depth of 40 cm.The simulated values corresponding to  Figure 7a,c shows the behavior of water fluxes in the soil profile during and after a rainfall event.At 10:00 (Figure 7a), water infiltrated the soil profile at a Darcian velocity of approximately 1.36 cm h −1 .This corresponded to the input data that were identical for all three variants.The flux decreases with increasing depth.The downward flux reached soil layer 2 in variants 1 and 2 within the first hour of the rainfall event.In variant 3, downward flux was observed only in the first 10 cm (Figure 7a).Two hours later (12:00 a.m., Figure 7b), variants 1 and 2 display a downward flux between −0.2 and −0.5 cm h −1 in the vadose zone up to 50-cm depth.In the capillary fringe at 55-cm depth, the values decrease down to −0.08 cm h −1 throughout layer 3, as determined by the R in LBC (Figure 7b).Variant 3 displays smaller downward fluxes of −0.4 cm h −1 in layer 1 up to 25-cm depth.Beyond, variant 3 displays a small upward flux of 0.02 cm h −1 in layer 3 (Figure 7b).In the next time step at 02:00 p.m., variants 1 and 2 display no relevant fluxes (Figure 7c).In contrast, variant 3 displays an upward flux of 0.47 cm h −1 in soil layer 3 (Figure 7f).In this time step, WTD increases in variant 3 because WTD is the LBC and an input for this model variant.
The simulated θ value of 0.5 cm 3 cm −3 in variant 3 is still lower than the θ values in variants 1 and 2 for layers 1 and 2. The progression of the θ values in variant 3, together with the fluxes, displayed a slower but gradually advancing wetting front compared to variants 1 and 2.

DISCUSSION
Shallow WTD conditions can readily supply moisture to the root zone, resulting in larger soil water pressure heads (h) within the unsaturated zone compared to sites with a lower WTD (Brauer et al., 2018).However, during the calibration period, the tensiometer-measured h 30 values reached as low as −368 cm (until sensor failure), indicating a reduced or intermittent water supply to the root zone in the following summer period.Variants 1 and 2 were able to satisfactorily describe the measured h 30 data in soil layer 1 during both the calibration and validation periods, with NSE values greater than 0.7 and IoA values exceeding 0.98 (Table 7, Figure 5).The retention curves of variants 1 and 2 exhibited a less steep increase, with n values less than or equal to 1.2 for soil layer 1, in contrast to n = 1.39 for variant 3 (Figure 4, Table 6).Variants 1 and 2 outperformed variant 3 in simulating h 30 during both the calibration and validation periods.Only variant 3 simulated θ 30 values comparable to the measured θ 30 values and displayed the best-fit measures for the calibration and validation periods.Variant 3 also produced the best results for h 60 because the WTD remained above 60 cm for a significant portion of the time series.Variants 1 and 2 provided good overall modeling results for the WTD, with NSE values exceeding 0.60 and IoA values greater than 0.98.This is a significant outcome since the WTD is a sensitive variable in the water balance of shallow groundwater sites and plays a role in various eco-hydrological processes (Bechtold et al., 2014).Small changes in water storage can lead to relatively large fluctuations in WTD at these sites, especially under near-saturated conditions (van Dam & Feddes, 2000).Therefore, inaccuracies in soil hydraulic models or their parameters, initial conditions, and layer boundaries can have a substantial impact on the accuracy of simulations.The simulated underestimation of the measured WTD both before and after significant rainfall events in June and August cannot be attributed to measurement errors because the WTD measurements were verified through measured soil water pressure heads h 60 and h 90 (Figure 3).Visually, the simulated WTD closely followed the storage change curve (Figure 5).This suggests that the processes within the vadose zone were not fully replicated by the model, leading to the observed discrepancies.
Variant 3 consistently underestimated the modeled inflow (R in ) throughout the observation period (Figure 5), indicating an overestimation of the water storage change in this variant.This overestimation was primarily caused by an overestimated available porosity in the model.At the beginning of the time series (Figure 3), the lysimeter was nearly saturated and then dried out during the summer months due to high ET a .The model's overestimated soil water storage capacity, which allowed for a larger soil water extraction, led to the underestimation of the simulated inflow at the lower boundary.Parameter θ r was considered less relevant (Wallor, Herrmann et al., 2018) because conditions in peatland rarely reach dry levels.However, the results presented here (Figure 5a) underscore the importance of the hydraulically active porosity, defined by θ s and θ r (Gnatowski et al., 2010), in more accurately simulating θ and flux time series values at shallow groundwater sites.
This study characterized the influence of LBC on simulated fluxes at shallow groundwater sites.A flux LBC (variants 1 and 2) combined with a flux upper BC provided fluxes at both ends of the soil profile.Consequently, the result is a perfect simulation of the storage change of the entire soil Vadose Zone Journal profile (Figure 5c) because there is no degree of freedom for the water budget components.However, this does not mean that the fluxes, soil water pressure heads, and volumetric water content were correct within the entire profile, especially in the unsaturated zone.For example, on June 6, 2019, variants 1 and 2 simulated an upward flux from the bottom to the surface in the entire soil profile (Figure 6) and an almost correct WTD (Figure 5b).The highest fluxes are at −40-cm depth at the bottom of the rooting zone.Simultaneously, the modeled volumetric water contents overestimate the measured values in the upper part of layer 1 (0.4 cm 3 cm −3 vs. 0.55 cm 3 cm −3 , Figure 6d-f).This shows that despite the seemingly correct modeling of the fluxes and the WT, other processes have to be considered, as discussed below.
In contrast, the variant utilizing the pressure head LBC (variant 3) encountered issues when simulating flows at the lower boundary that were not in accordance with measurements, especially when the pressure head LBC experienced significant changes (Figures 6 and 7).For instance, in Figure 6, variant 3 predicted a small downward flux of −0.006 cm h −1 in layer 3, resulting in an outflow at the lower boundary.However, measurements at this point indicated an inflow at the lower boundary in the opposite direction.In variant 3, the decreasing WTD served as an input to the model and led to the simulation of downward fluxes, essentially reversing the causality of the process.The underlying reason for the different behaviors of fluxes and WTD, depending on the selected LBC, lies in the distinct water storage characteristics of the model variants.The overestimation of available porosity in variant 3 caused additional drainage of the soil profile to reach the prescribed WTD change.Another example of a discrepancy between the simulated and measured fluxes in variant 3 was observed during a rainfall event on June 10, 2019 (Figure 7).While the measured WTD increased due to the rainfall event, variant 3 did not predict a downward flux reaching the capillary layer or the WT.Instead, the entire flux was stored in the upper centimeters of the first layer due to the overestimated storage capacity.Consequently, an inflow was simulated at the lower boundary to replicate the prescribed increase in the WTD.
Another important aspect to consider is the influence of root zone depth on vertical fluxes within the soil profile (Figure 6) during periods of significant transpiration (T a ) extraction, reaching up to 0.11 cm h −1 .For all variants, the simulated upward flux was most pronounced at the lower end of the simulated root zone, typically at a depth of approximately 40 cm.The rooting depth parameterization was based on field observations.From this depth to the surface, the simulated flux gradually decreased.Capillary flow from the saturated zone alone was insufficient to fulfill the entire root water extraction demand, necessitating the simulation of soilstorage water extraction in the vadose zone.Another notable factor impacting the simulated vertical fluxes and volumet-ric water content (θ) in this shallow groundwater site was the parameterization of root water uptake, particularly during periods with high T a .The rooting depth and root density play a crucial role in determining whether there is a continuous upward flux from the shallow WT and capillary fringe; the root model parameters affect the upward flux and the contribution of water extraction from the unsaturated zone (Grimaldi et al., 2015).Additionally, adjusting the value of the water stress index, also known as the root adaptability factor (Jarvis, 1989), could potentially enhance the simulations for variants 1 and 2 (Figure 5).
Discrepancies between measured and simulated data may also arise from potential issues with the measurement system.Notably, the lysimeter exhibited significant fluctuations in the WT during extreme rainfall events, which could have various explanations.For instance, the rapid increase in WT height (over 115 cm in 7 h) observed during the rainfall event on August 25, 2019, could be attributed to macropore flow in the lysimeter, a phenomenon not accurately captured by the model.However, this increase was so abrupt that it resembled the Lisse effect (Weeks, 2002).The Lisse effect involves the trapping of air between a shallow groundwater table and the wetting front resulting from a rainfall event, which can lead to a substantial rise in the WT in groundwater observation wells (Zhang et al., 2009).In this study, the WT measured in an observation well served as one of the inputs to the lysimeter control system (Dietrich et al., 2016), potentially causing the recorded increase in the WT in the lysimeter to be misleading.Another explanation for the pronounced WT rise could be short-term air entrapment within soil pores, resulting in hysteresis effects not accounted for in the modeling.Nevertheless, HYDRUS-1D can only simulate long-term effects of soil air during pore water drainage and imbibition when hysteresis is explicitly incorporated (Šimůnek et al., 2013).The error in simulated inflow (R in ) for variant 3 during the largest rainfall event also suggests the possibility of air entrapment effects.While a significant WT peak was enforced (as a LBC) in this variant, no corresponding substantial inflow was observed in the measurements (Figure 5a,b).It is important to consider this discrepancy when modeling vertical flows, especially in the context of reactive geochemical modeling applications, when employing a WT LBC.
A bimodal model (Durner, 1994) was implemented to consider the dual-porosity system observed in peatland soils (Dettmann et al., 2014;Kleimeier et al., 2014;Liu et al., 2017).The implementation of a bimodal soil water retention curve was not successful in achieving a larger improvement in variant 2 compared to the unimodal curve in variant 1.However, this contradicts the results of previous studies conducted on different types of peatland soils (Dettmann et al., 2014;Dimitrov et al., 2010;Rezanezhad et al., 2012;Weber et al., 2017).The present study suggests that the parameterization of SHPs based on a unimodal model, root zone depth, and selection of LBCs has a larger influence on the simulation of fluxes and volumetric water content (θ) in the vadose zone of a shallow groundwater site than the implementation of the bimodal model (Durner, 1994).Addressing these issues using methods such as global inverse simulation (Vereecken et al., 2016) and root zone parameterization would lead to even better simulation results, whereby a parameterized bimodal model could lead to greater improvements in critical events such as heavy rainfall.
The phenomena of air entrapment after heavy rainfall and overestimated porosity with the resulting underestimated inflow in variant 3 could only be observed here because of the available data on boundary fluxes and water storage changes in the lysimeter and their subsequent input into the model.In the literature, many model setups do not use the flux input values for all BCs.This means that the model setup is not a closed system in which there are no degrees of freedom for the components of the water balance.While they set a flux upper BC, they did not set a flux lower BC, but rather a variable pressure head lower BC (Akhtar et al., 2013;Forkutsa et al., 2009;Kelleners et al., 2005;Schwärzel, Šimůnek, van Genuchten et al., 2006).A lower BC variable flux is less common in the literature, as it is less important in sites with deep WT without a vertical upward flux from the WT into the vadose zone.Beyond that, a variable flux lower BC is more complex to measure and requires special equipment, such as a groundwater lysimeter (Dietrich et al., 2016).The issue with the more common setup and the lack of water balance and water storage change measurements is that the errors in water balance modeling may have moved from state variables towards water storage change.If the calibrated model contains an incorrect water-balance component, the simulation of fluxes throughout the soil column may be unsatisfactory.This was visible in variant 3 in this study (Figure 5a), where, despite good calibration results, the total inflow was underestimated, which affected the fluxes in the entire soil column.This error would not have been visible without the flux measurements at the lower boundary of the lysimeter.
The different BCs that can be measured and selected have advantages and disadvantages that affect model performance.The variable flux LBC, as in variants 1 and 2, defines the relative water storage change of the soil profile in the lysimeter and the fluxes in the saturated zone, whereas the soil water pressure head LBC, representing the WTD in shallow groundwater tables, as in variant 3, does not.However, with variants 1 and 2, there were errors in the WTD simulation, as the relative water storage change was known, but not the distribution of the water storage change in the soil profile.The pressure head LBC in shallow groundwater sites always defines the size and boundary of the vadose zone, which would improve θ and h simulation, whereas any errors in WTD simulation in variants 1 and 2 shift the boundary of the vadose zone and influence θ and h simulation negatively.In fact, variant 3 did perform best at simulating θ 30 and θ 60 , but not h 30 .It should be noted that calibrated model variant 3 achieved its results with an underestimated influx at the lower boundary, and thus, an incorrect water storage change.In this study, the hydraulically active available porosity of each soil layer was estimated from the inverse optimization of the SHPs.Supplementary available porosity measurements may improve this factor; however, in layered soils, there may still be larger inaccuracies in the overall soil water storage.However, the setup of variant 3, with a measured and defined variable WTD at the lower boundary, is a more practical alternative when lysimeter-based data are unavailable, such as in Forkutsa et al. (2009).
The measurements used in this study contained varying degrees of data uncertainty.We considered the water content measurements to be the most critical values because the profile probes only represent a limited measurement volume and can have large measurement errors, particularly in organic soils (Dietrich & Steidl, 2021).The errors were addressed by the soil-specific field calibration of the sensors used in this study, as described in Dietrich and Steidl (2021).Under variably saturated and water-unsaturated conditions with high suction, the measured tensiometer-based h values at water saturation agreed well with the measured WTD values.However, the tensiometer at 30-cm depth displayed measurement failures; therefore, the correlation between water content and soil water pressure head could not be recorded well for soil layer 1 under these conditions.The soil water pressure head and water content were measured in several different soil volumes.Owing to the spatial distances between sensors and the smallscale heterogeneity in the soil, the functional relationship between measurements deviates from expectations.
In this study, the estimated values of K s and SHPs, including the bimodal function, were inversely determined using local inverse optimization based on the presented procedure for the selection of initial parameters.However, local optimization can produce a nonunique solution, which is a local optimum derived from the initial parameters that may not be the globally optimal parameter set (Vereecken et al., 2016).Global inverse modeling can be applied to determine even more optimal retention and hydraulic conductivity curves (Hopmans et al., 2002;Iden & Durner, 2007;Kelleners et al., 2005).

CONCLUSIONS
The primary aim of this study was to characterize the soil hydraulic processes at a degraded peatland site with a shallow groundwater table.This characterization was essential as a precursor for quantifying solute transport and geochemical processes within such systems.The study specifically examined the influence of uni-and bimodal porosity models and different LBCs on quantitative aspects of the water balance, including water budget components, WTD, soil water pressure head (h), volumetric water content (θ), and vertical water fluxes.The study employed a numerical model calibrated using data obtained from a groundwater lysimeter.Key findings from this study underscore the critical role of selecting an appropriate LBC in conjunction with accurately characterizing the SHPs when modeling such complex systems.Notably, as follows: • Choosing measured fluxes at the lower boundary as the LBC offers the advantage of ensuring precise representation of water balance components, as it limits the degrees of freedom for these components.Conversely, utilizing measured water pressure head or WTD as the LBC can result in erroneous partitioning between inflow and outflow, thereby affecting the simulation of fluxes throughout the soil profile.

A C K N O W L E D G M E N T S
We thank our colleague, Mario Weipert, for supervising the lysimeter.We also thank the anonymous reviewers for their helpful comments.Special thanks go to Prof. Wolfgang Durner for the detailed review that helped us present our results more precisely and improve the manuscript.

C O N F L I C T O F I N T E R E S T S T A T E M E N T
The authors declare no conflicts of interest.

F I G U R E 3
Measured meteorological and lysimeter data for the period May to December.(a) cumulative potential evapotranspiration (ET p cum.) and precipitation (P cum.),(b) cumulative water flux at upper boundary (E a -P cum.) calculated from difference between actual evaporation E a and precipitation P, cumulative actual transpiration T a cum., and cumulative water flux at lower boundary (R in cum.) (positive values mean a flow into the soil monolith, negative values mean a flow out of the soil monolith), (c) water storage change ΔS, (d) water table depth (WTD), (e) soil water pressure head in 30-, 60-, and 90-cm depth h 30 , h 60 , h 90 , and (f) volumetric water contents in 30-, 40-, 50-, 60-, and 90-cm depth θ 30 , θ 40 , θ 50 , θ 60 , θ 90 .Please note the different scales of the y-axis for (d) WTD and (e) soil water pressure head.

F
Cumulated measured and simulated water fluxes at the lower boundary cum.R in (positive values mean a flow in the soil monolith, negative values mean a flow out of the soil monolith) (a), water table depth WTD (b), water storage change ΔS (c), soil water pressure heads in 30-and 60-cm depth h 30 (d) and h 60 (e), and volumetric water contents in 30-and 60-cm depth θ 30 (f) and θ 60 (g).The depiction of h 30 (d) has a different scale than that of h 60 (e) and the depiction of θ 30 (f) has a different scale than that of θ 60 (g).
Simulated Darcian velocity of vertical fluxes (a, b, and c) and simulated and measured volumetric water content θ (d, e, and f) in the upper 100 cm of the soil profile every 2 h on June 6, 2019, for each model variant.The selected date represents a high actual evapotranspiration (ET a ) and resulting upward flux during summer.Upward fluxes are defined positive, downward fluxes are defined negative.The diagrams are restricted to −100-cm depth for better visibility.Values of fluxes below this depth remain stable.E a , actual evaporation; P, precipitation; T a , actual transpiration.

Figure
Figure6d-f show that ET a causes a simulated θ decrease mainly in layer 2 for variants 1 and 2 and mainly in layer 1 for variant 3.Figure7a,c shows the behavior of water fluxes in the soil profile during and after a rainfall event.At 10:00 (Figure7a), water infiltrated the soil profile at a Darcian velocity of approximately 1.36 cm h −1 .This corresponded to the input data that were identical for all three variants.The flux decreases with increasing depth.The downward flux reached soil layer 2 in variants 1 and 2 within the first hour of the rainfall event.In variant 3, downward flux was observed only in the first 10 cm (Figure7a).Two hours later (12:00 a.m., Figure7b), variants 1 and 2 display a downward flux between −0.2 and −0.5 cm h −1 in the vadose zone up to 50-cm depth.In the capillary fringe at 55-cm depth, the values decrease down to −0.08 cm h −1 throughout layer 3, as determined by the R in LBC (Figure7b).Variant 3 displays smaller downward fluxes of −0.4 cm h −1 in layer 1 up to 25-cm depth.Beyond, variant 3 displays a small upward flux of 0.02 cm h −1 in layer 3 (Figure7b).In the next time step at 02:00 p.m., variants 1 and 2 display no relevant fluxes (Figure7c).In contrast, variant 3 displays an upward flux of 0.47 cm h −1 in soil layer 3 (Figure7f).In this time step, WTD increases in variant 3 because WTD is the LBC and an input for this model variant.

F
Simulated Darcian velocity of vertical fluxes (a, b, and c) and simulated and measured volumetric water content θ (d, e, and f) in the upper 100 cm of the soil profile every 2 h on the June 10, 2019, for each model variant.The selected date represents infiltration during a 5-h, 3-cm rainfall event.Upward fluxes are defined positive, downward fluxes are defined negative.The diagrams are restricted to −100-cm depth for better visibility, the fluxes remain stable below this depth.E a , actual evaporation; P, precipitation; T a , actual transpiration.

T A B L E 2
Overview of the Spreewald wetland lysimeter and weather monitoring station and the spatial and temporal resolution of measurements and data intervals for simulations.
c LAI was interpolated linearly between measurements.
T A B L E 5Abbreviationsa Determined using RETC.
• Even in models employing flux-based LBCs, incorrect unsaturated zone fluxes can occur if the soil hydraulic model is inadequate.For shallow groundwater sites, imposing a specified flux at the lower boundary can help maintain accurate fluxes within the saturated zone up to the vadose zone boundary.Notably, the model variant featuring a pressure head LBC showed superior performance compared to the flux-based LBC variant.Achieving more precise verification of model results necessitates the availability of more comprehensive, spatially resolved, and temporally frequent measurements across all layers of the soil monolith.In this study, the absence of data for a critical middle layer may have contributed to some level of uncertainty.