Testing the EPIC Richards submodel for simulating soil water dynamics under different bottom boundary conditions

Most biogeochemical models simulate water dynamics using the tipping bucket approach, which has been often found to be too simplistic to represent vadose zone dynamics adequately under shallow groundwater conditions. Recently, a solution to the Richards equation using the Mualem–van Genuchten model (Rich‐vGM) has been added into the EPIC (Environmental Policy Integrated Climate) model to address this shortfall. Its performance was tested using lysimeters operating under free drainage (FD) and at a shallow water table (60‐ [WT60] and 120‐cm depth [WT120]). Model accuracy was also compared with the upgraded tipping bucket‐based method implemented into EPIC (the variable saturation hydraulic conductivity method [VSHC]). Soil water content (SWC) data were split into calibration and validation subsets. Model evaluation also included annual evapotranspiration (ET), percolation (PRK), and upward water movements to assess underlying soil water balance factors. The submodels provided accurate and similar results upon comparison with SWC measures under FD (Nash–Sutcliffe coefficient [NSE] = 0.26 and 0.61 using VSHC and Rich‐vGM, respectively). The Rich‐vGM model accurately reproduced observed SWC and ET (e.g., NSE = 0.70 and percentage bias [PBIAS] = −3.7% for WT120, respectively) although it slightly overestimated PRK (PBIAS = 47.8%, on average). Instead, VSHC proved unable to correctly simulate shallow groundwater conditions (e.g., NSE = −1.85 for WT60 SWC). Under shallow groundwater conditions, the Rich‐vGM method is recommended, despite the additional data required and the need to define the bottom boundary conditions according to water table fluctuations. In conclusion, the Richards solver introduced and tested in EPIC improved the model's ability to represent complex biophysical and biogeochemical processes in terrestrial ecosystems associated with the hydrological balance.


INTRODUCTION
Since their inception in the 1960s, process-based models have been among the best tools to advance scientific understanding and support policy decisions . Such models rely on a myriad of equations to simulate agrisystem changes in response to external drivers (weather and management practices), and emerging interactions among other system components (e.g., soil and landscape properties). As expected, the most accurate models result from integrating daily or hourly updates into the calculations characterizing the flows and transformations of water, C, and N. Conversely, to reduce code complexity, input parameter complexity, and execution time, simplified algorithms are used to supply approximate representations of the biogeochemical processes (Zhao et al., 2019).
Simpler models rely primarily on the "capacity or tipping bucket" method, a widely used modeling approach to describe how water moves through the soil profile. The method uses the concept of field capacity (FC), meaning each soil layer has a specific water storage capacity. When soil water content (SWC) exceeds FC, or when macropore-held water approximates air-entry potential, gravitational drainage occurs. This is the approach used by some of the most utilized biogeochemical models, (e.g., DeNitrification-DeComposition [DNDC]; Li et al., 1992), although exceptions embedding more physically based methods also exist, such as the Richards equation in Cropsyst (Dhakar et al., 2018) and the Agricultural Production Systems sIMulator (APSIM; Huth et al., 2012). Moreover models that are coupled with Richards-based hydrologic models were also extensively developed (Kanda et al., 2018;Siad et al., 2019) .
Although current models simulate the process of drainage quite well, modeling the upward movement of water has been more challenging (Emerman, 1996;Tenreiro et al., 2020). Among the exceptions, irrigation scheduling models such as the Agricultural Water Productivity Model for Shallow Groundwater (AWPM-SG; Gao et al., 2017) have been implemented with simulation of water upward movement. Compared with the tipping bucket method, Richards-based approaches are physically based and allow flow to be driven either upward or downward depending on water potential gradients (Zha et al., 2019).
In biogeochemical models, the tipping bucket approach generally well simulates water dynamics at low SWC and deep water table levels (Liu et al., 2011;Wang et al., 2011). In contrast, inadequate simulations can occur under high SWC conditions because of the overestimated rapid drainage and the often-missing upward flux representations (Diekkrüger et al., 1995;Li et al., 2006). Given its direct link to hydraulic conductivity also, texture may affect model performance, which is among the main drivers for describing the drainage

Core Ideas
• EPIC Richards solver was tested in a 9-yr experiment under different groundwater conditions. • The new Richards algorithms were compared with the tipping-bucket based method. • Both methods predicted similar soil water content under free drainage. • Water content and ET were better predicted by Richards under shallow water table. • Under specific conditions, percolation was overestimated for both approaches.
response time. Consequently, the generally effective approach is notably inaccurate when simulating SWC redistribution in water tables shallower than 5 m, since modeled SWC outputs rarely exceed FC, whereas observed SWCs do (Kröbel et al., 2010). Furthermore, neglecting upward fluxes under shallow groundwater (GW) conditions may introduce errors into the description of crop evapotranspiration (ET) (Gao et al., 2017) or N transport within the soil profile (Morari et al., 2012). Indeed, ET is positively affected by water table depth, since water from deeper layers serves as an additional supply for ET (Luo & Sophocleous, 2010). This was also noticed by Schwaerzel and Bohl (2003), who experimentally observed water upward movements from deeper layers once the top water layer had evaporated. Regarding N fluxes, Brilli et al. (2017) reported SWC simulation errors to be the main cause of general discrepancies concerning C and N emissions in many of the studies they reviewed. Moreover, while using the DNDC model, Uzoma et al. (2015) found N 2 O emissions were underestimated due to the inability of the tipping bucket method to simulate SWC above the level of FC. Finally, in a three-model comparison by Zhang et al. (2015), some large spikes in N 2 O were highly underestimated when using the tipping bucket method. The Environmental Policy Integrated Climate model (EPIC) (Williams et al., 1984) was originally developed by the USDA to assess the impact of soil erosion on agricultural productivity. The model has evolved over the years by incorporating new functions and by extending its application to pesticide fate simulation and microbial denitrification (Gassman et al., 2004;Izaurralde et al., 2017). Recently, Doro et al. (2018) found that the EPIC tipping bucket method misrepresented percolation (PRK) dynamics within the soil profile. These results indicated that denitrification was also mischaracterized, which became evident when SWC rarely exceeded FC despite rapid water drainage. Furthermore, under such conditions, an anaerobic environment could not be established, limiting the development of a realistic denitrification submodel.
To remedy water dynamics simulation weaknesses in the EPIC model, two new submodules have recently been added to its code (Doro et al., 2018). The first submodule refines the original master PRK component, in which saturated hydraulic conductivity is nonlinearly reduced, according to SWC. The second submodule is a numerical solution of Richards equation following the approach presented by Ross (2003)  . Although the latter is more advanced from a physical property perspective, the former is an attractive alternative to more data-intense, complex simulation models. Indeed, simpler models minimize data collection efforts (Pachepsky et al., 2006), such those requested for long-term or extensive upscaling simulations.
Agrienvironmental measures are practices subsidized by the European Union to foster environmentally friendly farming techniques. Recently, EPIC has been identified as a suitable model to be coupled with a GIS platform to quantify the performances of agrienvironmental measures under the current Rural Development Programme (2014)(2015)(2016)(2017)(2018)(2019)(2020) in the Veneto Region of Northeast Italy. However, a considerable portion of the Veneto low plain is occupied by shallow GW (ARPAV, 2014) that fluctuates between 0.8 and 1.0 m in the winter and 2.0 and 3.0 m in the summer.
In 2011, an experimental station consisting of an array of lysimeters has been established to investigate the effect of water table depth on the main biogeochemical cycles. Due to their confined environment and the measuring instrumentation they may be equipped with, lysimeters are often considered excellent tools for studying water flow and chemical transport, and their data are often used to calibrate and test numerical models. Biogeochemical models capable of representing shallow water table conditions are necessary to describe the complex crop-soil interactions in the vadose zone of the Veneto low plain (Dal Ferro et al., 2016;Longo et al., 2021). Considering our literature review, our hypothesis is that the usage of the numerical solution of Richards equation is a prerequisite for adequate simulations of shallow water table conditions with EPIC, despite the higher requirements in terms of input data. The specific aims were (a) to reproduce a complex 9-yr-period lysimeter experiment that involved several treatments and factors combined along the period; (b) to calibrate and validate model water dynamics in terms of water content and evaluate the model accuracy for simulating water balance components and plant biomass production; and (c) to compare the Richardsbased submodel with the upgraded tipping bucket-based submodel under shallow water tables and free drainage (FD) conditions.

Experimental data description
The data used to test EPIC were collected from three experiments conducted between June 2011 and December 2019 at the experimental farm of the University of Padova in northeastern Italy (45˚19′ N, 11˚31′ E, 8 m asl).
The study site was originally set up in 1984 and consisted of 20 1-m × 1-m × 1.5-m (length × width × depth) nonweighable drainable lysimeters, made of reinforced concrete (Giardini et al., 1988). The lysimeters were buried into the soil to avoid artifacts caused by temperature gradients. To prevent border effects, crops were also grown outside lysimeters, creating a 2-m buffer.
The principle of communicating vessels was used to control water table conditions (Supplemental Figure S1). The bottom of each lysimeter is funnel-shaped and connected to an equally tall external column (150-cm height) with a 3-cm-diam. hose to avoid the effect of capillary action, fitted with a valve to regulate both the water table level and leached discharge.
The column was equipped with a water level sensor, to avoid the water level to deviating more than ±10 cm from the set reference level. Accordingly, well water was manually added by mean of the external columns to compensate upward flux, when water level was at least 10 cm lower than the reference level. On the contrary, PRK water was collected through the bottom valve when, following rainfall or irrigation, the water level raised >10 cm compared with the reference level. Both terms were measured with a 0.01-L precision on a daily basis in spring-summer and twice a week in autumn-winter.
Each lysimeter was filled in 1984 with soil excavated from the adjacent experimental farm using a method that preserved the original soil horizons. To facilitate water drainage and prevent soil washout, a 15-cm layer of 30-to-50-mm-diam. gravel was used to cover the bottom of each lysimeter. The lysimeter soil is a Fluvi-Calcaric Cambisol (CMcf; WRB, 2015; Table 1) and is representative of ∼50% of the lower Venetian plain. Soil particle size distribution (sand, silt, and clay content) was determined using a particle size analyzer (Mastersizer 2000, Malvern Panalytical;Bittelli et al., 2019); pH and electrical conductivity were measured by an electrode in soil suspensions with a soil/water ratio of 1:2.5 (Kabała et al., 2016) and 1:2 (w/v) (ISO 11265; ISO 1994) respectively; total N was analyzed with the Kjeldahl method (ISO 11261;ISO 1995); soil organic C content was measured with an elemental analyzer (VARIO MACRO, Elementar Analysensysteme); total and active carbonates were measured with the Dietrich-Fruehling calcimeter (ISO 10693;ISO 1995) and ammonium oxalate titration method (Jeffery & Hutchinson, 1981), respectively; available P was measured with the T A B L E 1 Soil physical and chemical properties in the 0-to-60and 60-to-135-cm layers at the study site Olsen method (ISO 11263;ISO 1994); and exchangeable K was measured by the extraction with ammonium acetate (Schollenberger & Simon, 1945). Field capacity (−33 kPa) and permanent wilting point (WP, −1,500 kPa) were derived from the soil water retention curve, which was estimated through the wind's evaporation method and dewpoint potentiometer. In addition, soil texture, bulk density and soil organic C were available for each lysimeter in the following layers: 0-30, 30-55, 55-75, 75-95, and 95-135 cm.

Treatments and management techniques
The experiment tested several combinations across many years. Given the considerable portion of the Veneto Region occupied by shallow water table (ARPAV, 2014), three GW conditions were tested and remained fixed variables throughout all the trials (2011-2019): FD conditions, a shallow water table set at 120-cm depth (WT120), and one at 60-cm depth (WT60). The variables that differed in the trials included different fertilization levels, N organic fertilizer types, and different agricultural practices (Table 2). For the 2011-2014 years, we tested three GW conditions against two N input levels (250 vs. 378 kg N ha −1 yr −1 ). In 2015, the GW conditions were combined with two organic fertilizer types (solid digestate vs. slurry). In both experiments (2011-2014 and 2015), the main crop was maize (Zea mays L.), which was grown with continuous cropping practices and according to conventional tillage operations. Soil was left bare during winter. Finally, in 2018 and 2019, we considered the combination of the three GW conditions and two management systems (conservation vs. conventional) ( Table 2). Conservation agriculture protocol prescribed no-till, cover cropping practices, and crop rotations. Conventional agriculture consisted of simulated harrowing operations performed at a 25-cm depth through a spring spade crop before seeding, residue incorporation, and bare soil between the main crops. Main crops were the same for both treatments and included maize in 2018 and grain At the end of each growing season, crop grain and residues were collected and weighted. Crop grain and residue samples (i.e., 400 g) were dried at 65˚C in a forced draft oven for 72 h for dry weight determination.
In all cases, the trials followed a randomized block design with two replicates, such that 12 lysimeters were used (three bottom boundary conditions [BBC] × two treatment levels × two replicates). Surface water input was regulated throughout the 2011-2015 period by a plastic roof that automatically closed to cover the lysimeter site during rain events. The cover not only prevented natural uncontrolled rainfall into the lysimeters, but it also protected the crop from extreme weather, such as hailstorms. In order to prevent the heating up of the air above the lysimeter station, this automatically reopened just after rain events ended. During this period, water inputs were provided by a series of simulated rainfall events applying amounts of water through irrigation, which were previously precisely weighted. Summer simulated-rainfall events were applied according to the average crop water needs measured in FD treatments; in the other seasons, water input was randomly distributed over the typical rainy months. In both cases, manual water application was uniform among lysimeters. From 2016, the plastic roof was closed only for extreme weather events, and surface water inputs came from both irrigation and rainfall events, with rainfall measured by the site meteorological station through an ECRN-100 rain gauge (METER Group) connected to a EM50 datalogger (METER Group) allowing a resolution of 0.2 mm. The same station measured air temperature, relative humidity, solar irradiance, and wind speed and direction at a height of 2 m from the ground.
During each week from 2011 to 2013, an active time domain reflectance (TDR) sensor (Moisture Point MP-917, ESI Environmental Sensors) connected to previously installed vertical waveguides measured the SWC across three different soil layers (0-15, 15-30, and 60-90 cm). Starting in 2013, every lysimeter was equipped with an automated monitoring system that allowed the continuous measurement of SWC.
The system was composed of CS635 TDR probes (Campbell Scientific) installed at 15-, 30-, and 60-cm depths. At the bottom of each WT lysimeter, a pressure gauge (TL-700, TECNO.El) was added for registering the water level. The TDR and pressure gauge probes were connected to a CR-10X datalogger (Campbell Scientific) through a series of multiplexers. Measurements were taken every 30 min. Actual ET (mm) was estimated on an annual basis as the residual term of the water balance with: where IN is the surface water input (rainfall and irrigation, mm), PRK is the percolation (mm), QIN is the upflux from the water table (mm), and Wt 1 and Wt 2 are the soil water storage in the root zone (mm) at the beginning and end of the considered period. All the previous terms were experimentally observed. Free drainage lysimeters did not consider QIN term. Water balance losses such as runoff (Q) and subsurface flow (SSF) were not considered in the lysimeter experiment since they are null in this kind of setup.

Description of the simulation experiment
Information about file structure and content needed to perform EPIC runs with a FORTRAN-based EXE file can be obtained from http://epicapex.tamu.edu/manuals-andpublications/. The EPIC version used in this study (v1102) included the two main submodels for modeling soil water dynamics (Doro et al., 2018).
The most physically based method (a Richards solver using the Mualem-van Genuchten model [Rich-vGM]) is a solution of the Richards equation (Richards, 1931): where SWC is the soil water content (cm 3 cm −3 ), t is time (s), z is vertical coordinate (cm), K is hydraulic conductivity (cm s −1 ) and ψ is hydraulic capillary head (cm). The embedded algorithm is based on the soil-water flow model developed by Ross (2003). Unsaturated flux was described by the matric flux potential φ (Campbell, 1985;Gardner, 1958), given by: where ℎ is pressure head (cm), ℎ e is pressure head at air entry, λ and η define the shapes of the water retention and conductivity curves, S is degree of effective saturation, and ϕ e = s ℎ e ∕(1 − λη). This is also known as the Kirchhoff transform and is frequently used in both analytical and numerical work because it linearizes the matric potential term in Darcy's Law for vertical downward flow: where q is water flux (m s −1 ) and z is depth, which allows simple calculations of flux to be accurate over larger distances as opposed to using h. The implementation drives vertical water flows through soil layers on a subdaily variable time step after daily simulated infiltration and evaporation from the surface soil layer and root water extraction from root-penetrated soil layers. Horizontal inflows and outflows are assumed to balance out in the EPIC model due to its independent modeling framework.
Initially, the Rich-vGM module presented by Doro et al. (2018) included only free gravitational drainage as lower boundary conditions. In our study, the approach was expanded to include both FD and variable pressure head conditions. One-dimensional flow with FD was assumed for the FD treatment. Variable head conditions were assumed for the WT treatments, with the bottom hydraulic head set daily with pressure gauge values. For characterizing soil water retention and unsaturated hydraulic conductivity functions, the model parameterization follows the modified van Genuchten-Mualem model (Schaap & van Genuchten, 2006): where SWC r and SWC s are the residual and saturated water contents, respectively (cm 3 cm −3 ), α (cm −1 ) is related to the inverse of the air-entry pressure, ℎ (cm) is the water tension, L (-) is an empirical pore-connectivity parameter, n (-) is a measure of pore-size distribution, and SAT is the saturated hydraulic conductivity (cm s −1 ). Soil hydraulic characteristics can be provided as inputs or can be estimated by pedo-transfer functions included in the code. In this work, the former option was chosen because soil water retention parameters and hydraulic conductivity at the site were previously estimated by Dal Ferro and Morari (2018) at the depth of 0-20, 20-45, and 45-145 cm.
In addition to Rich-vGM, another submodel was recently implemented: the variable saturation hydraulic conductivity method (VSHC). The submodule, which is a refinement of the original PRK method, assumes a FD condition as BBC and does not estimate the permeability of individual soil layers (Doro et al., 2018). Water movement above the level of FC is limited based on soil water status using the following approach: where K eff and K SAT are effective and saturated hydraulic conductivities (mm h −1 ), ST is the actual SWC, PO is soil porosity, and ε is a soil characteristic exponent which can range from 1 to 5 according to model documentation. In addition, when water exceeds FC, upward movements are modeled as follows: where UF is upward flow (mm) from layer l j to layer l j − 1 , and −1 are water tensions (mm of water) in layers l j and l j−1 , and β 1 and β 2 are two coefficients calculated as a function of WP and FC. Similarly, is estimated as a function of WP, FC, and actual SWC (Doro et al., 2018) Coefficients inside Equations 8 and 9 cannot be modified by the user. Independently by the selection of the PRK submodel, daily ET is internally calculated by the model (Williams, 1990). Firstly, the potential ET is calculated according to the equation selected. This is used by the model to calculate the soil evaporation, adjusting the value to consider the rainfall intercepted by the soil cover and, when necessary, snow cover. Concurrently, plant transpiration is calculated as plant water uptake and is affected by the crop growth and the crop available water. Finally, actual ET is calculated as the sum of soil evaporation and plant transpiration. Here, potential ET (PET) was calculated with the Hargreaves method (Hargreaves & Samani, 1985). In EPIC, the original coefficient and exponent have been increased from 0.5 and 0.0023 to 0.6 and 0.0032, respectively, in order to improve measurement representations. Calculations did not differ among treatments and submodels.
Soil layer depths were discretized into the model to fit to the depth of SWC probes (i.e., 15, 30, and 60 cm). A total number of 10 soil layers were simulated, divided at the cumulative depth of 12.5, 17.5, 27.5, 32.5, 45, 57.5, 62.5, 75, 95, and 135 cm. With the exception of K SAT and water retention curve parameters, soil layers and characteristics were kept fixed among the two submodules (i.e., 12 soil profiles corresponding to the 12 lysimeters). Bottom drainage layer (135-150 cm) was not considered.
In EPIC, crop phenological growth is based on daily heat unit accumulation (Williams et al., 1989). Further details about the EPIC growth submodel (Williams et al., 1989) can be found in the supporting materials. Root distribution within the soil profile is not directly affected by the degree of water saturation in the growth submodel. Consequently, potential root depth of each crop (RDMX) was adjusted to be within the unsaturated soil profile (0-60 cm) for WT60 simulations. No adjustment was necessary for WT120 and FD, since maximum root growth depth was already within 120 cm and could not reach the saturated profile.
In addition to RDMX and the potential heat units, few other crop parameters affecting plant development were modified. The critical aeration factor can range from 0.0 to 1.0, where the former maximizes the stress and the latter nulls it. In our study, it was set to 1.0 for both maize and sorghum. Moreover, their harvest index and the biomass energy ratio were set to 0.60 and 45.00 kg ha −1 MJ −1 (Causarano et al., 2008;Giardini et al., 1998). Similarly, also the maximum potential leaf area index was modified and set to 6.5.

Calibration and evaluation of model performance
The calibration procedure for the Rich-vGM method focused on the water retention curve parameters (α, n, SWC r , SWC s , K SAT ), even though they were previously estimated by inverse modeling . In the optimization process, layers were grouped according to the previously estimated parameter value layers (i.e., three: 0-20, 20-45, and 45-135 cm). Parameter ranges were 0.003-0.030 cm −1 for α, 1.05-1.45 for n, 0.00-0.15 m 3 m −3 for SWC r , 0.40-0.50 m 3 m −3 for SWC s , and 0.0-100.0 mm h −1 for K SAT . Application of the VSHC method required calibration of ε in Equation 7 without distinction among the different GW conditions. The term ε is a general parameter not linked to a specific soil layer, which can vary from 1 to 5 according to model documentation (https://epicapex.tamu.edu/epic/). Optimization of the VSHC method was restricted to ε as a result WP, FC, and PO calibration producing values of no hydraulic significance for silt-loam soil. For instance, to increase the simulated SWC to match the observed values under WT60, the algorithm increased WP to 0.25 m 3 m −3 , against an observed 0.13 m 3 m −3 , FC to 0.45, against an observed 0.31 m 3 m −3 , and PO to 0.55 m 3 m −3 , against an observed 0.43 m 3 m −3 .
EPIC is often used to simulate open field conditions (Gassman et al., 2004), where water balance components are generally less available than SWC (QIN and PRK in particular). Therefore, we used 2013-2017 SWC data for model calibration and 2018-2019 SWC data for validation. In addition, validation was also performed on the basis of annual ET, PRK, QIN, and crop yield and residues of the whole simulation period (2011-2019).
The multiobjective genetic parameter estimation algorithm Non-dominated Sorting Genetic Algorithm II (NSGA-II; Deb et al., 2002) was used to find the best parameter set. The optimization procedures aimed to maximize the Nash-Sutcliffe coefficient of efficiency (NSE; Nash & Sutcliffe, 1970) of SWC by layers, so that the procedure simultaneously optimized three NSE values (Supplemental Table S1). The parameter set that maximized the average NSE was then selected. The NSE coefficient is expressed as follows: where Obs is the ith observation of the variable being evaluated, Sim is the ith simulated value of that variable, Obs is the mean of observed values, and n is the total number of observations. Model outcomes produce good approximations of field data when the index is >0, with NSE = 1 representing the best fit.
The NSE and the percent bias (PBIAS) were calculated to serve as metrics of model performance (Tonitto et al., 2010). The PBIAS measures the average tendency of the simulated data to be larger or smaller than the respective observed value and is expressed as: The optimal value of PBIAS is 0, with low magnitude values indicating an accurate model simulation. Positive values indicate underestimation bias, whereas negative values indicate overestimation bias.

Experimental results
The SWC in FD was highly variable during the experimental period. Values were as high as 0.41 m 3 m −3 and as low as 0.06 m 3 m −3 to produce an overall average of 0.28 m 3 m −3 (Supplemental Figure S2). Despite the wide range, water dynamics were seasonally similar for 2013 to 2015, years when maize was grown with continuous cropping practices.
The combination of bare soil and favorable weather conditions (low temperatures, high precipitation) led to high SWC values during the winters. Conversely, summer SWC values were low and remained so until autumn. Similar trends were observed in 2018 and 2019 when only a slight alteration was noted due to the presence of a rye winter cover crop. Under shallow water table conditions, prolonged periods of saturation or near-saturation were observed even at relatively low depths (30 cm). Notably, SWC values in both WT120 and WT60 during the summer dipped only slightly and never fell below 0.25 m 3 m −3 (Supplemental Figure S2). Only in the topsoil during the warmer months did SWC drop rapidly (to 0.12 m 3 m −3 in WT120 and 0.23 m 3 m −3 in WT60), driven by high ET and insufficient upward GW fluxes. On the contrary, SWC under WT120 showed low dynamics, highlighting the high contributions given by capillary rise even considering the plant water uptake. Water surface input provided by natural rainfall and artificial rainfall (through irrigation) were constant among treatments and averaged 868.5 mm yr −1 , ranging from 639 mm (2015) and 1,219 mm (2017). Water balances were influenced by ET (Supplemental Table S2), for which observed mean annual values ranked lowest in FD (814.1 mm yr −1 ), followed by WT120 and WT60, which were relatively similar (1,061.4 and 1,120.3 mm yr −1 , respectively), averaging 1,007.6 mm yr −1 in total. Overall, annual ET values ranged from 608.3 mm yr −1 in FD (2015) to 1,326.9 mm yr −1 in WT60 (2012). Across the T A B L E 3 Average yield and residues observed (Obs) and simulated (Sim) using Richards equation (Rich-vGM)  simulation period, average PRK was similar between WT60 and WT120 (218.8 and 213.2 mm yr −1 ) but differed considerably in FD (37.8 mm yr −1 ). On average, PRK equaled 165.5 mm yr −1 despite being highly variable among lysimeters (standard deviation was ±90.6 mm yr −1 ) (Supplemental Table S2). Across the simulation period, QIN values varied between 287.2 and 669.8 mm yr −1 under shallow water table conditions, with similar annual averages for WT60 and WT120 (473.7 and 413.8 mm yr −1 ). Therefore, the annual average contribution to ET from the upward movement of water under WT120 and WT60 equaled 45%. Input provided to FD as QIN were null. On average, 313.4 mm yr −1 was input into the system as QIN. On the other hand, losses due to subsurface flow (SSF) and runoff (Q) both equaled zero due to the confined environment of lysimeters. Finally, average annual changes in SWC storage were slightly greater under FD (16.6 mm yr −1 ) than WT120 and WT60 (7.7 and 3.1 mm yr −1 , respectively). Groundwater conditions affected also plant biomass production. Average yield dry matter was indeed greater under WT120 and WT60 (11.0 and 9.9 Mg ha −1 , respectively), whereas lower values were observed in FD (8.0 Mg ha −1 ) (Table 3). Similarly, residue dry matter was affected by the presence and depth of the water table, with the optimal conditions observed again under WT120 (16.8 Mg ha −1 ) and FD and WT60 averaging 14.4 and 15.0 Mg ha −1 , respectively (Table 3).

Modeling results: Calibration
In estimating water retention curve parameters for the Rich-vGM method, the NSGA-II algorithm generally increased α, particularly for the top soil layer of WT120 and WT60, which varied from 0.010 to 0.021 and from 0.013 to 0.027 (Supplemental Table S3). The SWC r value was also increased, espe-  Table S3). The ϵ value resulting from NSGA-II optimization in the VSHC method equaled the maximum of its range (ε = 5) recommended for equation physical reliability. The Rich-vGM method satisfactorily described water dynamics for most of the period (Figure 1), showing particular improved capacity to represent elevated SWC conditions under FD and WT120. At depths of 15 and 30 cm, simulated SWC followed the WT60 > WT120 > FD trend, albeit with a slight overprediction at low values in the top soil layer in WT60 and WT120 (Figure 1). The submodel correctly reproduced FD dynamics with PBIAS in the calibration period equal to 0.7%. The satisfactory results were also confirmed by NSE coefficient (0.61). Good performances were also demonstrated under shallow water table conditions. Here, average simulated SWC values were extremely close to the observed ones (Table 4), ranging from 0.35 to 0.39 m 3 m −3 (simulated) and 0.36 to 0.39 m 3 m −3 (observed). The model failed to reproduce only a few events at a depth of 15 cm. For instance, EPIC underpredicted SWC during summer 2014 due to an overprediction of plant biomass, and in turn, ET (Supplemental Figure S2). The opposite occurred in summer 2015, when simulated values correctly captured the observed trend, but translated them at higher values. Overall, the Rich-vGM submodel performed well at values close to saturation, especially at the deepest layers. This was demonstrated by PBIAS, which never exceeded 1%.
The VSHC submodule predicted SWC better in FD than in WT120 and WT60, as indicated by NSE values equal Vadose Zone Journal F I G U R E 1 Daily soil water contents (SWCs, m 3 m −3 ) observed and simulated using the Richards solution (Rich-vGM) approach in lysimeters under free drainage (FD), and with water tables set at 120-and 60-cm depths (WT120 and WT60) to 0.26 (Table 4), against values equal to −1.12 and −3.06 in WT120, and −1.85 and −2.33 in WT60. Indeed, simulated dynamics were comparable with those observed in FD (PBIAS = −4.4%), especially at depths of 15 and 30 cm ( Figure 2). Occasionally on bare soil, SWC was overpredicted, such as during 2016 and during the 2013-2015 winters. In WT120, SWC conditions were comparable only in the top soil layer, and a large underrepresentation was observed at greater depth. Overall, measured SWC was 0.35 m 3 m −3 , on average, while simulated average 0.28 m 3 m −3 . Consequently, the WT120 PBIAS was the worst among all the simulations. Not only PBIAS metric was bad, but also NSE (−1.12) (Table 4). Similarly, a large underestimation was noted in WT60 despite its better PBIAS compared to those in WT120 (−12.3%). Predictions improved at a depth of 30 cm, for which winter dynamics mirrored observed dynamics (Figure 2). Average SWCs in the calibration period equaled 0.34 vs. measured values of 0.39 m 3 m −3 . The NSE was negative and equal to −1.85 (Table 4).

Rich-vGM VSHC
To conclude, the lack of constant head as BBC under VSHC submodule often led to excessive SWC underestimation, phenomena that was less common using Rich-vGM Vadose Zone Journal F I G U R E 2 Daily soil water contents (SWCs, m 3 m −3 ) observed and simulated using the variable saturation hydraulic conductivity (VSHC) approach in lysimeters under free drainage (FD), and with water tables set at 120-and 60-cm depths (WT120 and WT60) (Supplemental Figure S2). Indeed, the VSHC-SWC tended to be higher than did the Rich-vGM SWC in the topsoil of FD, whereas the contrary was observed for WT60 and WT120 (Supplemental Figure S2). As for deeper layers, the Rich-vGM SWC was always higher than the VSHC-SWC across all treatments.

Modeling results: Validation
Generally, performances were similar for calibration and validation periods for both submodules. Good agreements between observed and simulated water dynamics were produced using the Rich-vGM method, despite SWC values were overpredicted for summer 2018 in FD (at all depth) and topsoil WT120 and WT60 (Supplemental Figure S2). Consequently, metrics were particularly good under shallow water table, where NSE equaled 0.68 for both WT60 and WT120 and percentage error did not overcome ±1%. Slightly worse performances, despite still satisfactory, were produced in FD (Figure 1), which NSE was 0.35 and values were overestimated on average of 8.5% (0.27 m 3 m −3 observed SWC vs. 0.29 m 3 m −3 simulated).
Despite VSHC led to satisfactory results under FD in the calibration period (Figure 2), performances worsened in the validation period (Table 4) for SWC, with a NSE coefficient equal to −0.38. Indeed, the submodule highly overpredicted topsoil SWC from June 2018 to April 2019. On the other hand, the contrary was observed at 30-and 60-cm depth, where EPIC underpredicted SWC in both summer seasons considered. Similar to calibration period, SWC conditions in WT120 were comparable only at 15 cm, whereas a large underprediction was observed at the deeper layers where SWC often remained stable at WP. Therefore, observed average was 0.36 m 3 m −3 and simulated 0.26 m 3 m −3 , which resulted in a low PBIAS (−28.0%). Similar phenomena were observed for WT60, with NSE = −2.33 and PBIAS = −16.9%. Under these conditions, some overprediction events were again simulated as a consequence of upward flux overrepresentations.
Lysimeter PET estimated with Hargreaves method, which did not vary among submodels, ranged from 1,382 to 1,545 mm yr −1 . The VSHC approach generally underestimated ET losses under all treatments (Figure 3), when simulations averaged 710.7 mm yr −1 vs. observations equal to 1,007.6 mm yr −1 , corresponding to PBIAS = −31.8%. Considering GW conditions, closer results for ET were observed in FD (694.6 mm yr −1 vs. observed 814.1 mm yr −1 ), whereas greater discrepancies were produced under WT conditions. On average, simulations appeared more realistic using the Rich-vGM approach as simulated ET losses rose to 998.7 mm yr −1 (Figure 3). Consequently, all the statistical indices improved as compared with those using the VSHC approach (Supplemental Table S2).
In terms of PRK, good agreement was observed in WT60 and WT120 using VSHC, despite its high overprediction in FD (320.9 mm yr −1 simulated vs. 37.8 mm yr −1 observed) as a consequence of the unpredicted plant biomass Vadose Zone Journal F I G U R E 3 Water fluxes (ET, PRK, and QIN; mm yr −1 ) observed vs. simulated with the variable saturation hydraulic conductivity (VSHC) and the Richards solution (Rich-vGM) approaches. ET: evapotranspiration, PRK: percolation, QIN: upflux from the water-table, FD: free drainage, WT120: water-able set at 120-cm depth, WT60: water table set at 60-cm depth ( Figure 3). This fouled its results overall, making it miss the average trend (PBIAS = 41.4%) and lead to a negative NSE (−1.14). The PRK was also overpredicted using Rich-vGM method, although to a lesser degree. Even though performances were overall satisfactory (NSE = 0.30), PRK was also overpredicted using Rich-vGM (PBIAS = 47.8%), in particular for WT60 (218.8 mm yr −1 observed vs. 351.8 mm yr −1 simulated).
Variable saturation hydraulic conductivity showed no ability to reproduce QIN, with high underpredictions (450.3 mm yr −1 observed vs. 70.4 mm yr −1 simulated) (Figure 3; Supplemental Table S2). The general PBIAS value dropped to −89.3%, which caused the QIN/ET ratio to be considerably underestimated (10%). On the other hand, the Rich-vGM submodel represented QIN (average = 537.7 mm yr −1 ) better and delivered a slightly general improvement in all indicators (Supplemental Table S2), despite a general overprediction. In particular, QIN resulted highly overpredicted in WT60, which average trend was 57.6% greater than observed one. Consequently, the QIN/ET ratio averaged 54%.
The EPIC water balance also considered SSF and Q losses, items that were null in the lysimeter experiments due to the input management and the lysimeter installation, which restricted movement to only vertical water movement. As such, much effort was taken to lower simulated values as much as possible by minimizing lateral hydraulic conductivity and runoff curve numbers in the input files. The results were satisfactory using the VSHC method, as Q and SSF were negligible (0 and 0.01 mm yr −1 ). It was more difficult to lower Q (43.6 mm yr −1 ) using the Rich-vGM submodel.
Finally, the accurate SWC simulated with Rich-vGM submodel proved the method ability at reproducing yield dry matter. Indeed, NSE averaged 0.36 (Table 3), with better values in WT60 and WT120, despite the model slightly missing the remarkably high yields experimentally observed, resulting in a trend underestimation of −12.6 and −5.0%, respectively. Conversely, NSE was negative for FD (−0.06), but the trend was captured by the model resulting in average 8.1 Mg ha −1 against the 8.0 Mg ha −1 observed (PBIAS = 1.3%). Model performance was poorer in terms of residues. Indeed, the extremely high values that were observed were not reproduced by EPIC, which produced a PBIAS coefficient equal to −51.7%, on average (15.3 Mg ha −1 observed vs. 7.4 Mg ha −1 simulated) (Table 3). Simulations were overall similar among GW conditions, ranging from 6.9 (FD) to 7.8 Mg ha −1 (WT120), and NSE was always negative.
Concerning VSHC, performance worsened for both yield and residues, for which NSE never reached positive values. Indeed, the SWC underestimation led to low yield, which was notably similar among GW, ranging from 6.4 Mg ha −1 in FD and WT60 to 6.7 Mg ha −1 in WT120, which resulted in a very low PBIAS (−33.7%, on average; Table 3). The PBIAS was even lower considering the model ability to reproduce residues, and the lack of constant head as BBC made the coefficient drop to −57.8% while NSE equaled −2.55. Notably, residues values were similar to yield ones, varying from 6.4 (FD) to 6.8 Mg yr −1 (WT60), sensibly lower than observed ones, the range of which were 14.4−16.8 Mg yr −1 .

DISCUSSION
In our study, we demonstrated the ability of the Richards solution integrated in EPIC to realistically represent relationships between GW level and root zone SWC. Aside from different model comparisons, tipping bucketand Richards-based approaches have been rarely compared within a single model. Eitzinger et al. (2004) Faria and Bowen (2003), after introducing the steady-state flow Darcy-Buckingham equation into DSSAT, concluded that model performance was significantly enhanced over that delivered with the tipping bucket approach. Furthermore, direct comparison of DNDC and Daisy by Kröbel et al. (2010) reported that the former, which uses a tipping bucket model, was unsuitable for accurate soil water simulation. In the case of the latter, which uses Richards equation, they found it quite suitable. Recently, Shelia et al. (2018) reported that integrating the HYDRUS-1D model into DSSAT provided better SWCs compared with the tipping bucket method, despite the appearance of submodel discrepancies from crop variables. In the recent comparison of the performances of 12 models made by Groh et al. (2020), models that used Richards-based methods achieved better results in terms of SWC and net water flux at the lysimeter bottom than models that did not integrate them. Moreover, this study further validated our methodology by emphasizing the importance of lysimeters for providing the necessary soil-related data for testing soil-vegetationatmospheric models.
Although Rich-vGM showed better ability to simulate SWC than VSHC, a few issues still need attention. For instance, the more complex code behind the former method did not allow us to completely null Q and SSF. Furthermore, QIN was greatly overestimated under WT60, which in turn affected PRK, leading to a great overrepresentation of the latter. The relationship between the two terms is also substantiated by the results observed using VSHC under WT60 and WT120. Indeed, although SWCs were not accurately reproduced and a rapid drainage was simulated, PRK was adequately captured by the model. This suggests that the underestimation of the QIN, and therefore the decrease of the overall water input, might be significantly relevant for the process. If adequate simulations of PRK are relevant for accurate simulations of N and P leaching, precise SWC predictions are essential not only for plant production studies but also for greenhouse gas emissions and other important biogeochemical processes (Brilli et al., 2017).
A substantial difference in the two submodels was made by including different lower boundary conditions. Indeed, only the Rich-vGM submodel included both FD and hydraulic head conditions, whereas VSHC included just the former. This explains also why FD performances were the only satisfactory for both approaches, whereas VSHC failed to reproduce SWC under shallow GW.
The VSHC submodel calibration initially included the PO, FC, and WP coefficients; however, the automatic algorithm maximized them to values without hydraulic significance. In fact, by the maximization of these terms, the model increased the SWC in the attempt to describe the SWC profile associated with shallow water conditions. The submodel inaccuracy to simulate the high SWC was also substantiated in the final automated calibration, when the term ε reached its upper limit. With the increase of this term, the procedure maximized the hydraulic conductivity reduction as a function of SWC with the aim of decreasing water PRK. As stated by the submodule developers, VSHC performance may improve if the exponent ε would be differentiate by soil layers and based on their physical characteristics (Doro et al., 2018). Furthermore, performance could also benefit from specific calibrations of coefficients inside Equations 8 and 9, but currently this is not possible for users.
In this work, model performance was improved by supplying Rich-vGM code with BBC set as hydraulic head pressures measured daily by transducers. Rapid drainage was simulated when FD was set as BBC (Supplemental Figure S2). It is not common to drive Richards-based models with measured variable hydraulic head data at such a frequency, especially when EPIC is upscaled to a regional or territorial tier. Nonetheless, EPIC could be further improved by adding a subroutine to simulate daily water table fluctuations, such as that implemented in EPIC-WT by Sabbagh et al. (1991). This modified version of the EPIC model was developed with the specific aim of reproducing drainage system behavior. Water table depth calculations in EPIC-WT differ from those in EPIC, as fluctuations are not only a function of precipitation and GW storage, but also of hourly soil water leaving the system (i.e., Q, SSF, and PRK). Integration of upward flux data exiting the system would surely enhance calculation of the water table depth in the model. Such an enhancement could greatly increase Rich-vGM method relevance, since outputs would not be dependent from the provision of BBC values, boosting the significance for large-scale modeling.
Inside the Rich-vGM submodel, water retention curve parameters can also be estimated by inner pedotransfer equations, such as those of Saxton and Rawls (2006) and Weynants et al. (2009). Similarly, models that interpolate water retention curves using basic soil properties such as FC, WP, and SWC s also exist. The SWIM3 model embedded into APSIM makes use of these characteristics for fitting curves through the usage of monotonic cubic Hermite splines (Huth et al., 2012). In our study, the ability of the Richards method to outperform the VSHC method was also due to the modification of the water retention curve parameters. This necessary adjustment was substantiated by a study in northern Italy by Bonfante et al. (2010), in which a comparison of the SWAP, MACRO, and CropSyst models was conducted. They found that CropSyst needed adjusted parameters, even though the implemented solution was based on the Campbell soil-waterretention model (Campbell, 1985). Similarly, Žydelis et al. (2018) found that the AgroC model required calibration of K SAT values, since it did not satisfactorily reproduce the in situ SWC using laboratory values. This might also be related to the differences between field and laboratory water retention curves. For instance, Herbrich and Gerke (2017) found that field drying retention curves were generally steeper than those estimated in the laboratory, which was claimed by the authors to be related to the smaller SWCs at which drying branches start at field conditions. Experimental actual ET was calculated from the lysimeter water balance. Despite being a simplistic method, Denager et al. (2020) found results from this approach to be comparable with more sophisticated methods like eddy covariance. However, ET calculations can be affected by rain gauge-related errors, since they have difficulty measuring non-rainfall water such as dew and hoar frost (Hoffmann et al., 2016). Furthermore, some authors highlighted that non-rainfall water such as dew and hoar frost can be relevant and substantial in the water balance (Groh et al., 2018). In our lysimeters, shallow GW conditions increased experimental ET by 34% as compared with FD. This result confirmed those of Allmaras et al. (1975), who reported that upward flux contributed 40-60% of ET. Camarotto et al. (2018) also estimated contributions ranging between 21 and 60% for the Veneto Region. Finally, Gao et al. (2017) simulated the QIN/ET ratio of different water table depths in a semiarid area by coupling the EPIC crop growth model to WIPE (Watershed Irrigation Potential Estimation), a model to estimate soil moisture in the root and vadose zones. Study results indicated a contribution to ET by the water table as high as 41% when the latter was at a 100-cm depth. Moreover, positive values were measured down to 200-cm depth when no irrigation was supplied. In their modeling work, Kollet and Maxwell (2008) found that, if in the 1-to-5-m depth range, the GW dynamics influenced the components of energy balance, which coincided with the root zone defined in the model. In another modeling study, same authors found that the "critical zone" for observing accentuated correlations between water table depth and land-surface energy was instead in the 2-to-5-m range depth . On the contrary, at deeper water table depths, GW was dislinked from the land surface, which was in dynamic equilibrium with atmospheric forcing. Moreover, the concept of a "critical zone" was substantiated by Groh et al. (2016), who reported that ET and PRK were highly sensitive to water table depth only if included in a specific depth range that varies depending on soil texture, as confirmed by Soylu et al. (2011).
This could also partially explain why observed ET did not considerably differ between WT60 and WT120. Indeed, occurrence of steady state upward flux can be assumed at the lysimeters under a shallow water table. In these conditions, the ET rate is mainly dictated by the external atmospheric drivers compared with the water-transmitting properties (Skaggs, 1978;Hillel, 1998).
Under a shallow water table, annual ET reached values extremely and unusually high (up to 1,500 mm yr −1 ) for the Veneto Region, where the average Hargreaves ET0 is around 800 mm yr −1 (Veneto Region Environmental Protection Agency, ARPAV). The uncommon ET experimentally observed may be attributed to factors related to lysimeterspecific artifacts, such as border and oasis effect. Other authors have reported unexpectedly high ET from lysimeters (Oberholzer et al., 2017) despite the relatively large size (1-m 2 surface area) and that the surroundings were covered by lawn. In our study, the same crop grown in the lysimeters was also grown outside them to prevent border and oasis effects. Indeed, such effects typically occur if plant and atmospheric conditions of the lysimeters are different from those of the surrounding (Wegehenkel & Gerke, 2013).
In general, modeled ET appears to be notably sensitive to water submodels, as suggested by Doro et al. (2018). Indeed, as a direct function of SWC in the root zone, marked differences were observed in a comparison of the two approaches, especially in WT60 and WT120, where the Richards average ET was >50% greater than that with VSHC. Considering the latter method, the misrepresentation of the QIN component was one of the main reasons for the underprediction of ET values. The satisfactory ET predictions using Rich-vGM may also suggest that constraining root depth depending on the water table depth may be an effective way to simulate crop water uptake. Using DNDC, Smith et al. (2019) attributed poor SWC predictions to the lack of an algorithm for root distribution over the full soil profile, as well as an inability of the model to simulate a heterogeneous soil profile and water table. In the recently revised DNDC version (Smith et al., 2020), root density functions were included and shown to be particularly effective at improving predictions of soil water dynamics. In this same version of the model, the tipping bucket method was upgraded to estimate a SWC-dependent effective unsaturated conductivity, analogous to the VSHC approach in EPIC. Wang et al. (2014) combined EPIC and CHAIN_2D to improve simulation of dynamic root growth, root water uptake, and crop yield under furrow irrigation. In this work, the enhancements were aimed at horizontal dynamics and vertical root growth was not constrained in accordance with SWC. In our study, aeration stress factor was maximized to 1.0 for both maize and grain sorghum in order to null the effect of waterlogging on root growth. In addition to EPIC, other models that integrate ways for reducing total biomass production, potential transpiration, or photosynthesis in concomitance of shallow GW or waterlogging events have been developed (Shaw et al., 2013;Liu et al., 2020). On the other hand, enhanced root growth modules able to account for excessive moisture effects are quite rare (Ebrahimi-Mollabashi et al., 2019).
Given the widely usage of EPIC for crop simulation studies (Perez-Quezada et al., 2003;Balkovič et al., 2013;Wang et al., 2015), great importance should be placed on the ability of the model to reproduce with accuracy crop yield and residues. In our study, the VSHC submodel highly underestimated yield values as a consequence of underpredicted SWC and, therefore, water stress conditions. Specific and deep calibration of crop parameters could probably improve submodule performance, but this was beyond the goal of the study. On the other hand, crop residue production was unrepresented also using Rich-vGM, as confirmed by the similar and low NSE value. This can be attributed to the EPIC growth model, which is more simplistic with respect to models specifically designed for modeling plant biomass and phenology, such as DSSAT (Yakoub et al., 2017).
Although the coupling of hydraulic models to agricultural system models has been well reported (Kanda et al., 2018;Siad et al., 2019), it can be speculated that using Richards solutions as an internal option is more efficient and better aligned with common applications of agricultural systems model users. For instance, the incorporation of the EPIC crop growth module into the HYDRUS-1D model (Wang et al., 2015) omitted the full capabilities of EPIC in simulating biogeochemical processes (e.g., soil water erosion, soil organic C dynamics). Additionally, coupling models may lead to further conceptualization and computational problems, which increases the uncertainty levels at modeling (Antle et al., 2001).

CONCLUSIONS
Simulation models capable of considering shallow water table conditions are necessary to describe the complex crop-soil interactions that occur in the vadose zone under such conditions. Recent enhancements to the hydrological EPIC components, specifically the tipping bucket-and Richards-based submodels, were investigated in a set of lysimeter experiments under different shallow water table conditions. Results suggested that both submodels are reliable for predicting processes that affect the Venetian soil profiles under deep GW conditions, where FD conditions prevail. In this case, the tipping bucket-based model may offer a preferable combination of simplicity and computational efficiency compared with the Richards solution. Under shallow GW conditions, more important differences were revealed. In these instances, only the Richards method was capable of reproducing accurate SWC above the water table, despite upward water movement from the water table being overestimated. Under these circumstances, the Richards method is strongly encouraged, even after recognizing the additional data required to run the model and the need to define the bottom boundaries according to water table fluctuations. Furthermore, the method represents a unique functionality added to a widely used model that is not available via existing model couplings.
Coupling a modified EPIC with a GIS platform will create a powerful tool for quantification of agrienvironmental performance measures and identification of sustainable agricultural systems, even in peculiar environments such as that of the low plain of Veneto Region.

A C K N O W L E D G M E N T S
This research was co-funded by the Rural Development Programme for Veneto 2014-2020.