Soil oxygen dynamics: a key mediator of tile drainage impacts on coupled hydrological, biogeochemical, and crop systems

Tile drainage removes excess water and is an essential, widely adopted management practice to enhance crop productivity in the US Midwest and throughout the world. Tile drainage has been shown to significantly change hydrological and biogeochemical cycles by lowering the water table and reducing the residence time of soil water, although examining the complex interactions and feedbacks in an integrated hydrology–biogeochemistry–crop system remains elusive. Oxygen dynamics are critical to unraveling these interactions and have been ignored or oversimplified in existing models. Understanding these impacts is essential, particularly so because tile drainage has been highlighted as an adaptation under projected wetter springs and drier summers in the changing climate in the US Midwest. We used the ecosys model that uniquely incorporates first-principle soil oxygen dynamics and crop oxygen uptake mechanisms to quantify the impacts of tile drainage on hydrological and biogeochemical cycles and crop growth in corn–soybean rotation fields. The model was validated with data from a multi-treatment, multi-year experiment in Washington, IA. The relative root mean square error (rRMSE) for the corn and soybean yield in validation is 5.66 % and 12.57 %, respectively. The Pearson coefficient (r) of the monthly tile flow during the growing season is 0.78. Plant oxygen stress turns out as an emergent property of the equilibrium between the soil oxygen supply and biological demand. The impact of tile drainage on the system is achieved through a series of coupled feedback mechanisms. The model results show that tile drainage reduces the soil water content and enhances soil oxygenation. It additionally increases the subsurface discharge and elevates inorganic nitrogen leaching, with seasonal variations influenced by climate and crop phenology. The improved aerobic condition alleviates crop oxygen stress during wet springs, thereby promoting crop root growth during the early growth stage. The development of greater root density, in turn, mitigates water stress during dry summers, leading to an overall increase in the crop yield by ∼6 %. These functions indicate the potential of tile drainage in bolstering crop resilience to climate change and the use of this modeling tool for large-scale assessments of tile drainage. The model reveals the underlying causal mechanisms that drive the agroecosystem response to drainage on the coupled hydrology, biogeochemistry, and crop system dynamics.

1 Introduction

Agricultural subsurface drainage, commonly referred to as tile drainage, is one of the most important agricultural management practices that enables timely planting and enhances crop productivity in the US Midwest (Moore, 2016; Shen et al., 2013; Skaggs et al., 1994). Over 80 % of the tile-drained fields in the US are concentrated in six US Midwestern states, covering one-third of the region's cropland (NASS-USDA, 2017; Valayamkunnath et al., 2020; Zulauf and Brown, 2019). Notably, nearly half of the fields in the 3I states (i.e., Iowa, Indiana, and Illinois) are tile-drained, and the adoption rate of tile drainage continues to grow (NASS-USDA, 2017; Valayamkunnath et al., 2020; Zulauf and Brown, 2019). Tile drainage improves drainage conditions by removing excessive water and lowering the water table (Kalita et al., 2007), which benefits seed germination and crop growth (Ashraf, 2012; Nóia Júnior et al., 2023). Tile drainage also helps reduce the risk of delays in crop planting by enabling the timely operation of farm machinery during wet spring months, and thereby extending the crop growing period (Kucharik, 2008; Shirzaei et al., 2021). Tile drainage will become even more important under climate change, as the US Midwest is expected to experience wetter springs with more frequent and intense late-spring storms (Lesk et al., 2016; Li et al., 2019; Lobell et al., 2014; Schmidhuber and Tubiello, 2007; Seneviratne et al., 2022; Zhou et al., 2022). Thus, understanding the tile drainage impacts and managing the hydrological condition over the Midwestern agroecosystem are critically needed.

Extensive studies have explored the impact of tile drainage from different perspectives, such as hydrology, soil biogeochemistry, and crop growth. Hydrologically, tile drainage induces changes in both water storage and water fluxes (Blann et al., 2009; Boland-Brien et al., 2014; Hanrahan et al., 2020). Tile drainage has been shown to lower the water table, reduce the soil water content, and increase the temporal soil water storage capacity, which might enhance percolation, reduce surface runoff, and mitigate flooding at the field scale (Blann et al., 2009; Rahman et al., 2014; Skaggs et al., 1994; Yimer et al., 2023). These local-scale changes additionally alter watershed hydrology (Hansen et al., 2013; Miller and Lyon, 2021; Woo and Kumar, 2019). The impacts of tile drainage on hydrology are complicated by the interacting environmental conditions and management practices, soil properties, and antecedent soil moisture (Blann et al., 2009; Cain et al., 2022; Stops et al., 2022; Thomas et al., 2016; Wiskow and van der Ploeg, 2003). For instance, tile drainage may either increase baseflow or result in a flashier hydrograph, depending on the specific meteorological and physical characteristics (Adelsperger et al., 2023; Miller and Lyon, 2021; Schilling et al., 2012; Schilling and Helmers, 2008; Thomas et al., 2016; Valayamkunnath et al., 2022). The impact of tile drainage on evapotranspiration (ET) is generally more associated with land conversion (Ma et al., 2023; Wiskow and van der Ploeg, 2003), and the ET in tile-drained fields may be similar to that in undrained fields in the same crop systems (Khand et al., 2017; Yang et al., 2017). Tile drainage degrades stream water quality by increasing both field nitrogen and phosphate leaching (Castellano et al., 2019; David et al., 2010; Grenon et al., 2021; Ma et al., 2023; Ren et al., 2022; Sims et al., 1998). Further, tile drainage fosters a more aerobic soil condition, which would largely alter soil microbe activities, i.e., mineralization and immobilization (Brown et al., 2017; Jacinthe et al., 2015; Kumar et al., 2014). Notably, the impacts of tile drainage on both hydrology and biogeochemistry exhibit seasonal variation (Lam et al., 2016; Macrae et al., 2007; Ma et al., 2023; Williams et al., 2015). Despite the substantial attention on tile drainage, the impacts of tile drainage on hydrology, biogeochemistry, and crop growth are often studied separately within disciplinary boundaries, preventing an integrated understanding of their complex interactions.

Figure 1 Hypothesis for how tile drainage impacts hydrology, soil biogeochemistry, and crop growth. Tile drainage directly changes the field hydrology, and then changes the soil biogeochemistry and crop growth. Among those, soil oxygen is a key mediator that links hydrology, soil biogeochemistry, and crop growth in the coupled system.

Process-based models are promising for their capabilities of integrating physical, chemical, and biological processes, thus providing a cost-efficient and time-efficient means to advance scientific understanding and offer decision/policy support compared to field experiments (Jones et al., 2017). Tile drainage modeling has a long history and has recently gained attention for its integration into comprehensive hydrology and land surface models (Bailey et al., 2022; De Schepper and Therrien, 2017; Hansen et al., 2013; Li et al., 2010; Muma et al., 2017; Rumph Frederiksen and Molina-Navarro, 2021; Smith et al., 2020; Valayamkunnath et al., 2022). However, many of these models are specialized for particular processes, and, in some cases, they either omit or oversimplify other critical processes. Hydrologic models primarily focus on hydrological responses but do not represent soil biogeochemistry and crop growth, such as the National Water Model, the Soil Water Assessment Tool (SWAT), and the Tsinghua Hydrological Model (THREW) (Arnold et al., 2012; Li et al., 2010; Valayamkunnath et al., 2022). Similarly, reactive transport models, such as the PFLOTRAN and Advanced Terrestrial Simulator (ATS), simulate water and nutrient transport and biogeochemical transformation but lack the capability of representing crop growth and agricultural management activities (Hammond et al., 2014; Li et al., 2017, 2021). Crop growth models have recently been used to illuminate the impacts of waterlogging, but they often lack a comprehensive representation of the interconnections between hydrology, crop dynamics, and soil biogeochemistry. For instance, many such models, such as the DRAINMOD model, the Soil-Water-Atmosphere-Plant (SWAP) model, the Agricultural Production Systems sIMulator (APSIM), the decision support system for agrotechnology transfer (DSSAT) model, and the Environmental Policy Integrated Climate (EPIC) model, do not adequately represent root respiration, a key process influencing root development, maintenance, and nutrient uptake. Instead, they rely on soil water content as a proxy for oxygen stress, potentially neglecting important nuances in the interconnections between these critical processes (Ebrahimi-Mollabashi et al., 2019; Liu et al., 2020; Pasley et al., 2020; Skaggs et al., 2012; Sharpley, 1990).

Representations of ecophysiological and biochemical mechanisms under excessive water are critical in using process-based models to understand the interlinkage and interaction between hydrology, biogeochemistry, and crop growth. The limited soil oxygen diffusion under water-saturated conditions has been recognized as a key factor altering crop growth and soil biogeochemistry (Elzenga and van Veen, 2010; Pan et al., 2020; Rubol et al., 2013). Restricted soil oxygen under excessive water will suppress root respiration, decrease root activity, and even lead to root senescence, which further affects crop yield (Pan et al., 2020). Soil oxygen availability regulates soil biogeochemistry and composition by altering the redox potential, influencing microbial processes, nutrient cycling, and the mobility of elements (Elzenga and van Veen, 2010; Rubol et al., 2013). Under saturated conditions, anaerobic microbial respiration is favored to consume electron acceptors such as nitrate, sulfate, or iron instead of oxygen for respiration, producing greenhouse gases, and changing the availability of essential nutrients for crop uptake. However, those processes are not well represented in the aforementioned process-based models.

Here we aim to use a process-based model, ecosys, that incorporates all the essential physical mechanisms, especially oxygen-related dynamics, to understand the role of tile drainage in the integrated hydrology–biogeochemistry–crop agroecosystem dynamics by addressing the following questions: (1) how does tile drainage alter the agroecosystem hydrology, biogeochemistry, and crop growth? (2) More importantly, how are those impacts on the three aspects interrelated? (3) How do seasonal precipitation patterns influence tile drainage and agricultural production? We hypothesize that tile drainage positively influences crop growth by altering in-field hydrology and soil biogeochemistry, with soil oxygen acting as the critical mediator that links these three components. Specifically, tile drainage improves soil aeration by reducing the soil water content. The resulting increase in oxygen availability stimulates root development and microbial activity, subsequently enhancing overall crop growth. In turn, this greater biological activity creates a feedback loop: enhanced microbial activity accelerates organic matter mineralization, and the root system alters the hydrology through increased water and nutrient uptake (Fig. 1). We further hypothesize that tile drainage could bolster agricultural production and potentially serve as an efficient adaptation strategy in the context of climate change. We first validate the ecosys model using data from a multi-year field experiment at a research and demonstration farm in Washington, Iowa. In Sect. 2, we provide an overview, and some key processes related to tile drainage and soil oxygen simulations in the ecosys model and introduce the data used in this study and hypothetical numerical experiment designs. Sect. 3 presents the model calibration and simulation results. In Sect. 4, we specifically discuss and answer the above-mentioned questions and draw conclusions in Sect. 5.

2 Model and data

2.1 Ecosys model

2.1.1 Overview of ecosys model

The ecosys model is an agroecosystem model with essential mechanistic representations of hydrology, soil biogeochemistry, and crop growth in the soil–vegetation–atmosphere continuum at the hourly time step (Figs. 2 and S1–S3 in the Supplement) (Grant, 2001; Grant et al., 2017). Ecosys can be configured to run in 1D, 2D, or 3D. For typical field-scale applications, it operates as a highly detailed, multi-layered 1D model. And 2D and 3D configurations have been employed to investigate topography control on ecosystem processes (Grant et al., 2017). It has shown promising performance in simulating water fluxes (e.g., ET), biogeochemistry (e.g., soil carbon storage, greenhouse gas emission), and crop growth (e.g., gross primary productivity and crop yield) in different cropping systems (Grant, 1993, 1998; Li et al., 2022; Mezbahuddin et al., 2016; Qin et al., 2021; Yang et al., 2022; Zhang et al., 2021; Zhou et al., 2021). In this study, we configured the model to study the impact of tile drainage to run in a 1D setting, which also ensures a greater numerical stability.

Figure 2 Examples of the interplays between hydrology, soil biogeochemistry, and crop growth in the ecosys model.

Ecosys simulates the movement of water through the soil–crop–atmosphere continuum, with the representation of crop interception of precipitation, irrigation, soil and residue evaporation, crop transpiration, infiltration, surface runoff, subsurface discharge, and snow (Fig. S1). All the water fluxes in both the soil and crop are driven by water potential and are tightly coupled with energy cycles (Grant et al., 1999; Mezbahuddin et al., 2016).

Soil biogeochemistry in ecosys is simulated by tracking the flow of carbon (C), nitrogen (N), and phosphorus (P) among various organic states within the soil (Fig. S2). The model represents organic matter into six organic states (i.e., solid organic matter, soluble organic matter, sorbed organic matter, acetate, microbial communities, and microbial residues). Each state is further divided into components with varying vulnerability to hydrolysis by microbial populations. Microbes are the agents that control the C, N, and P transformation. Microbial activity is simulated based on the energetics of oxidation-reduction reactions, driving processes such as decomposition, nitrification, denitrification, and methanogenesis. Meanwhile, the energy generated in those processes and the nutrients will also be used for microbe maintenance and growth. Microbes also undergo decomposition (Grant et al., 1993a, b).

Ecosys simulates crop growth by representing the crop as a collection of individual branches and organs. The growth of branches and organs is driven by the balance between carbon fixation through photosynthesis and carbon losses through respiration and senescence. Carbon fixation happens in the leaves via the Farquhar model, and the fixed carbon is then mobilized to other branches and organs (Grant, 1994). Water and nutrient uptake (i.e., of P and N) is simulated with a hierarchical root system (Grant, 1993, 1998) as affected by temperature, nutrient availability, and soil oxygen concentration. Similarly, the nutrients absorbed by the root (i.e., N and P) are remobilized to other branches and organs for crop growth. For example, the N mobilized to the leaves determines the specific activities and surficial concentrations of leaf rubisco and chlorophyll, further affecting the CO₂ fixation rate.

The three components (hydrology, soil biogeochemistry, and crop growth) are tightly interconnected within ecosys. For example, carbon assimilation in crop growth is tightly coupled to canopy transpiration, as stomatal conductance, affected by canopy turgor potential, determines both the transpiration rate and photosynthesis rate (Grant and Pattey, 1999). Root water uptake is the driver for root nutrient uptake in the dispersivity–diffusivity processes; see details in Sect. S3. The hydrological cycle in the model is also tightly linked to soil biogeochemistry. Soil water movement drives the movement of soil nutrients, determining nutrient vertical distribution and loss. In addition, the movement of water also drives the movement of soil gas, e.g., soil oxygen, and subsequently changes both root respiration and microbe activities; see details in Sects. S4–S8. Further, microbial activities control the release of nutrients from organic matter, influencing the availability of nitrogen and phosphorus for crop uptake, which dynamically links soil biogeochemistry with crop growth.

We used ecosys to evaluate the impact of tile drainage on field hydrology, soil biogeochemistry, and crop growth. We provide details about soil oxygen–related processes and tile drainage processes in the following section as soil oxygen is a critical component linking hydrology, biogeochemistry, and crop growth. More detailed processes of the various components, like ecosystem–atmosphere energy exchange, canopy carbon fixation, etc., of the ecosys model can be found on GitHub (https://github.com/jinyun1tang/ECOSYS/blob/master/ecosys_documentation.pdf, last access: January 2024; Grant, 2023).

2.1.2 Soil oxygen dynamics

Oxygen is represented in two phases in ecosys: the gaseous oxygen in the air-filled porosity and the dissolved oxygen in the soil water. The vertical transport of both gaseous and dissolved oxygen, the transfer between dissolved and gaseous oxygen, and the oxygen consumption by both the crop roots and soil microbes are explicitly represented in the model. The volatilization–dissolution transfer between dissolved and gaseous oxygen is driven by the oxygen difference in the two phases, and is determined by the diffusive transfer coefficient, and the air–water interfacial area.

$$\begin{array}{}\text{(1)}& Q_{\mathrm{d},{\mathrm{O}}_{\mathrm{2}}}=a_{\mathrm{g}}{D}_{\mathrm{d},{\mathrm{O}}_{\mathrm{2}}}\left(S_{{\mathrm{O}}_{\mathrm{2}}}^{\prime }{f}_{\mathrm{t},\mathrm{s},{\mathrm{O}}_{\mathrm{2}}}\left(C_{\mathrm{g},{\mathrm{O}}_{\mathrm{2}}}-C_{\mathrm{s},{\mathrm{O}}_{\mathrm{2}}}\right)\right)\end{array}$$

where $Q_{\mathrm{d},{\mathrm{O}}_{\mathrm{2}}}$ is the volatilization–dissolution of O₂ between the solute and gaseous phases [g m⁻² h]; $D_{\mathrm{d},{\mathrm{O}}_{\mathrm{2}}}$ is the O₂ volatilization–dissolution transfer coefficient [m² h⁻¹]; a_g is the air–water interfacial area [m² m⁻²]; $S_{{\mathrm{O}}_{\mathrm{2}}}^{\prime }$ is the Ostwald solubility coefficient of O₂ at 30 °C [–]; $f_{\mathrm{t},\mathrm{s},{\mathrm{O}}_{\mathrm{2}}}$ is the temperature dependence function of $S_{{\mathrm{O}}_{\mathrm{2}}}^{\prime }$ [–]; $C_{\mathrm{g},{\mathrm{O}}_{\mathrm{2}}}$ is the gaseous concentration of O₂ in the soil [g m⁻³]; $C_{\mathrm{s},{\mathrm{O}}_{\mathrm{2}}}$ is the corresponding solute concentration in the soil [g m⁻³].

The vertical transport of dissolved and gaseous oxygen in the soil is calculated from the convective–dispersive equation,

$$\begin{array}{}\text{(2)}& {\displaystyle }& {\displaystyle }{Q}_{\mathrm{g},{\mathrm{O}}_{\mathrm{2}}}=-Q_{\mathrm{w}}{C}_{\mathrm{g},{\mathrm{O}}_{\mathrm{2}}}+D_{\mathrm{g},{\mathrm{O}}_{\mathrm{2}}}{\displaystyle \frac{\partial C_{\mathrm{g},{\mathrm{O}}_{\mathrm{2}}}}{\partial L}}\text{(3)}& {\displaystyle }& {\displaystyle }{Q}_{\mathrm{s},{\mathrm{O}}_{\mathrm{2}}}=Q_{\mathrm{w}}{C}_{\mathrm{s}}+D_{\mathrm{s},{\mathrm{O}}_{\mathrm{2}}}{\displaystyle \frac{\partial C_{\mathrm{s},{\mathrm{O}}_{\mathrm{2}}}}{\partial L}}\text{(4)}& {\displaystyle }& {\displaystyle }{D}_{\mathrm{g},{\mathrm{O}}_{\mathrm{2}}}={\displaystyle \frac{D_{\mathrm{g},{\mathrm{O}}_{\mathrm{2}}}^{\prime }{f}_{\mathrm{t},\mathrm{g}}{\mathit{\theta }}_{\mathrm{g}}^{\mathrm{2}}}{{\mathit{\theta }}_{\mathrm{p}}^{\mathrm{0.67}}}}\text{(5)}& {\displaystyle }& {\displaystyle }{D}_{\mathrm{s},{\mathrm{O}}_{\mathrm{2}}}=D_q\left|Q_{\mathrm{w}}\right|+D_{\mathrm{a}}^{\prime }{f}_{\mathrm{t},\mathrm{s}}{\mathit{\theta }}_{\mathrm{w}}\mathit{\tau }\end{array}$$

where $Q_{\mathrm{g},{\mathrm{O}}_{\mathrm{2}}}$ is the gaseous flux of O₂ in the soil [g m⁻² h]; Q_w is the water flow rate in the soil [m³ m⁻² h⁻¹]; $C_{\mathrm{g},{\mathrm{O}}_{\mathrm{2}}}$ is the gaseous concentration of O₂ in the soil [g m⁻³]; $D_{\mathrm{g},{\mathrm{O}}_{\mathrm{2}}}$ is the gaseous diffusivity of O₂ in the soil [m² h⁻¹], determined by its gaseous diffusivity at 0 °C ($D_{\mathrm{g},{\mathrm{O}}_{\mathrm{2}}}^{\prime }$) [m² h⁻¹], temperature dependence function for gaseous diffusivity (f_t,g) [–], air-filled porosity (θ_g) [m³ m⁻³], and soil porosity (θ_p) [m³ m⁻³]; $\frac{\partial C_{\mathrm{g},{\mathrm{O}}_{\mathrm{2}}}}{\partial L}$ is the concentration gradient of gaseous O₂ in the soil [g m⁻³ m⁻¹]; $Q_{\mathrm{s},{\mathrm{O}}_{\mathrm{2}}}$ is the solute flux of O₂ in the soil [g m⁻² h]; $C_{\mathrm{s},{\mathrm{O}}_{\mathrm{2}}}$ is the solute concentration of O₂ in the soil [g m⁻³]; $D_{\mathrm{s},{\mathrm{O}}_{\mathrm{2}}}$ is the solute diffusivity of O₂ in the soil [m² h⁻¹], determined by the dispersivity in the soil (D_q) [m], its solute diffusivity at 0 °C ($D_{\mathrm{a}}^{\prime }$) [g m⁻² h], Q_w, temperature dependence function for solute diffusivity (f_t,s) [–], soil water-filled porosity (θ_w) [m³ m⁻³], and soil tortuosity (τ) [–].

Soil oxygen will be used by crops, mycorrhizal, and microbes for their maintenance and growth. See Sect. 2.1.3 for a detailed description of the crop oxygen uptake. See the Supplement and the online document on GitHub (https://github.com/jinyun1tang/ECOSYS/blob/master/ecosys_documentation.pdf) for a detailed description of the microbial growth and oxygen uptake.

2.1.3 Root respiration and crop oxygen demand

Roots play a critical role in crop growth by acquiring necessary resources, including water and nutrients (i.e., nitrogen, phosphorus, etc.), from the soil for crop development and stabilizing the crop body structure (Hodge et al., 2009). Understanding the interactions between the root system and soil is essential to quantify the impacts of different environmental factors on crop growth (Jin et al., 2020). Ecosys explicitly simulates the root system with a representation of vertical primary axes and horizontal secondary axes (details in Grant, 1993, 1998). In the model, root growth and maintenance are driven by root respiration, and the rate of root respiration at maximum turgor in each soil layer is controlled by the available carbon storage, soil moisture, temperature, oxygen availability, and nutrient status,

$$\begin{array}{}\text{(6)}& R_{\mathrm{T}}=Q_{\mathrm{R}}{C}_{\mathrm{R}}{f}_{\mathrm{t},\mathrm{R}}{f}_{\mathrm{o},\mathrm{R}}{f}_{\mathit{\lambda },\mathrm{R}}\end{array}$$

where R_T is the root respiration for maintenance and growth [gC m⁻² h⁻¹]; Q_R is the specific respiration of CH₂O [g g⁻¹ h⁻¹]; C_R is the nonstructural CH₂O in the root [gC m⁻²]; f_t,R is the temperature function for respiration [–]; f_o,R is the oxygen function for respiration, represented as the ratio of O₂ uptake to O₂ demand [–], and will be described in detail below; f_λ,R is the nutrient status function for respiration [–]. The actual respiration rate is further adjusted by the root turgor and soil strength (Grant, 1993, 1998). Nutrient uptake (NO${}_{\mathrm{3}}^{-}$, NH${}_{\mathrm{4}}^{+}$, PO${}_{\mathrm{4}}^{\mathrm{3}-}$) also respires CH₂O,

$$\begin{array}{}\text{(7)}& R_{\mathrm{U}}=\mathit{\alpha }\mathrm{\Sigma }{U}_{\mathit{\lambda }}\end{array}$$

where R_U is the respiration for nutrient uptake [gC m⁻² h⁻¹]; α is the specific respiration rate for nutrient uptake [–]; U_λ is the uptake rate of nutrient Z (NO${}_{\mathrm{3}}^{-}$, NH${}_{\mathrm{4}}^{+}$, PO${}_{\mathrm{4}}^{\mathrm{3}-}$) [gN m⁻² h⁻¹ or gP m⁻² h⁻¹]. The total root respiration is, then, the total respiration for root maintenance, root growth, and root nutrient uptake (R_T+R_U). The crop oxygen demand ($U_{\mathrm{O}}^{\prime }$) is defined as the oxygen uptake rate without soil oxygen limits,

$$\begin{array}{}\text{(8)}& U_{\mathrm{O}}^{\prime }={\displaystyle \frac{\frac{R_{\mathrm{T}}}{f_{\mathrm{o},\mathrm{R}}}+R_{\mathrm{U}}}{R_{\mathrm{Q}}}}\end{array}$$

where R_Q is the respiratory quotient [gC(gO₂)⁻¹].

Figure 3 Representation of (a) oxygen dynamics and root oxygen uptake, and (b) subsurface tile flow in the ecosys model. r_R is the root radius [m]; r_W is the thickness of the water film [m]; r_P is the radius of the root porous core [m]. D₁: water table depth in the field [m]; D₂: tile depth [m]; 2d₁: tile spacing [m]; θ_k: soil water content in the kth soil layer [m³ m⁻³]; θ_s,k: saturated soil water content in the kth soil layer [m³ m⁻³]. In this schematic, the field water table is in the ith soil layer, and tile drainage is installed in the mth soil layer. The vertical dots represent the intermediate soil layers not shown. The diagram only represents the case that the water table in the field is higher than the tile drainage.

The oxygen uptake rate in ecosys is controlled by both the soil oxygen supply (dissolved oxygen transport rates to the root surfaces) and the ability of roots to take up oxygen (active uptake rates at the root surface where respiration is modeled). The conceptualization of crop roots is depicted in Fig. 3a, with a porous core in the middle, surrounded by an aqueous zone where respiration happens, then encased in a water film. Gaseous and dissolved oxygen transport in both the root porous core and the soil contributes to root respiration. The movement of oxygen is assumed to be radial, and the rate of oxygen moving from the soil water to the root surface and the rate of oxygen moving from the aqueous zone of the root porous core to the root surface are obtained from Eqs. (9) and (10), respectively.

$$\begin{array}{}\text{(9)}& {\displaystyle }{U}_{\mathrm{O},\mathrm{s}}=U_{\mathrm{w}}{C}_{{\mathrm{O}}_{\mathrm{2}},\mathrm{s}}+\mathrm{2}\mathit{\pi }{D}_{{\mathrm{O}}_{\mathrm{2}}}L{\displaystyle \frac{\left(C_{{\mathrm{O}}_{\mathrm{2}},\mathrm{s}}-C_{{\mathrm{O}}_{\mathrm{2}},\mathrm{R}}\right)}{\left(\frac{r_{\mathrm{R}}+r_{\mathrm{W}}}{r_{\mathrm{R}}}\right)}}\text{(10)}& {\displaystyle }{U}_{\mathrm{O},\mathrm{P}}=\mathrm{2}\mathit{\pi }{D}_{{\mathrm{O}}_{\mathrm{2}}}L{\displaystyle \frac{C_{{\mathrm{O}}_{\mathrm{2}},\mathrm{R}}-C_{{\mathrm{O}}_{\mathrm{2}},\mathrm{P}}}{l\left(\frac{r_{\mathrm{R}}}{r_{\mathrm{P}}}\right)}}\end{array}$$

where U_O,s is the rate of oxygen uptake by the root from the soil [g m² h⁻¹]; U_O,P is the rate of oxygen uptake by the root from the root porous core [g m² h⁻¹]; U_w is the root water uptake rate [m³ m⁻² h⁻¹], determined by the soil and root water potential and root resistances (Grant, 1998); $D_{{\mathrm{O}}_{\mathrm{2}}}$ is the dispersivity–diffusivity of dissolved oxygen [m² h⁻¹] (Bresler, 1973); L is the root length [m m⁻²]; $C_{{\mathrm{O}}_{\mathrm{2}},\mathrm{s}}$ is the dissolved oxygen concentration in the soil [g m⁻³]; $C_{{\mathrm{O}}_{\mathrm{2}},\mathrm{R}}$ is the oxygen concentration at the respiration site [g m⁻³]; $C_{{\mathrm{O}}_{\mathrm{2}},\mathrm{P}}$ is the dissolved oxygen concentration in the root porous core [g m⁻³]; r_R is the root radius [m]; r_W is the thickness of the water film [m] (Kemper and Rollins, 1966); r_P is the radius of the root porous core [m]. The active oxygen uptake rate by the roots is modeled in the Michaelis–Menten format,

$$\begin{array}{}\text{(11)}& U_{\mathrm{O}}={\displaystyle \frac{U_{\mathrm{O}}^{\prime }{C}_{{\mathrm{O}}_{\mathrm{2}},\mathrm{R}}}{C_{{\mathrm{O}}_{\mathrm{2}},\mathrm{R}}+K_{\mathrm{O}}}}\end{array}$$

where U_O is the root oxygen uptake rate [g m² h⁻¹], K_O is the Michaelis–Menten constant for the root oxygen uptake [g m⁻³]. U_O is solved iteratively from Eqs. (6)–(8), with $U_{\mathrm{O}}=U_{\mathrm{O},\mathrm{s}}+U_{\mathrm{O},\mathrm{P}}$. All dissolved oxygen concentrations are driven by the oxygen transport in gaseous phases, and by dissolution from gaseous to aqueous phases in the soil and roots, which will be affected by the soil drainage conditions. Details of the oxygen transport and dissolution (i.e., solute and gaseous) in the soil and root could be found in Grant (1993). Then, the oxygen stress indicator (f_o,R) in Eq. (6) is defined as the ratio between the U_O and $U_{\mathrm{O}}^{\prime }$,

$$\begin{array}{}\text{(12)}& f_{\mathrm{o},\mathrm{R}}={\displaystyle \frac{U_{\mathrm{O}}}{U_{\mathrm{O}}^{\prime }}}.\end{array}$$

2.1.4 Water balance

The water balance in the simulated field is given by

$$\begin{array}{}\text{(13)}& \mathrm{\Delta }{S}_i=P_i-Q_i-R_i-{\mathrm{ET}}_i.\end{array}$$

All terms are defined as fluxes per unit area [m³ m⁻²] at the time step i. Specifically, ΔS_i is the change in the soil water storage; P_i is the precipitation; R_i is the surface runoff, referring to the water leaving the system above the soil surface, and it is estimated with Manning's equation (Sect. S1 in the Supplement); Q_i represents the total water exchange across the subsurface boundary and tile drainage. A positive value indicates outflow (subsurface discharge) and a negative value indicates inflow (subsurface recharge). Although the model operates on an hourly time step, in the following analysis, the subsurface exchange was aggregated and attributed to discharge and recharge on a daily basis. The lateral water exchange is controlled by a lateral subsurface boundary condition, with a specific external water table depth and a specific lateral distance over which the lateral subsurface water flow occurs (Fig. S4). The external water table depth represents the water table depth of the surrounding environment; for instance, the water table depth at the field boundary, channel, or nearby lakes.

2.1.5 Tile drainage

The vertical soil water flow is governed by the Richards' equation,

$$\begin{array}{}\text{(14)}& {\displaystyle \frac{\partial {\mathit{\theta }}_{\mathrm{w}}}{\partial t}}={\displaystyle \frac{\partial }{\partial z}}\left(K\left({\mathit{\theta }}_{\mathrm{w}}\right)\left({\displaystyle \frac{\partial h}{\partial z}}+\mathrm{1}\right)\right)-S\end{array}$$

where θ_w is the soil water content [m³ m⁻³]; K(θ_w) is the soil hydraulic conductance at θ_w [m h⁻¹]. S is the soil water sink term [m³ m⁻³ h⁻¹], including the crop and mycorrhizal water uptake, lateral water fluxes to/from the external water table, and lateral discharge to the tile drainage. To solve the equation, the soil column is discretized into several user-specified layers (Fig. 3b).

Water fluxes to the tile pipes are estimated laterally with Darcy's law in saturated soil layers,

$$\begin{array}{}\text{(15)}& v={\displaystyle \frac{K\mathrm{\Delta }D}{d}}\end{array}$$

where v is the flow velocity [m h⁻¹]; K is the saturated hydraulic conductance [m h⁻¹]; ΔD is the pressure drop [m] over a distance d [m]. The pressure drop is defined as the difference between the internal water table depth and the tile depth (D₂–D₁), and the distance is then defined as half of the tile space (d₁). Similar to the Dupuit–Forchheimer assumption, we assume water flow horizontally. Thus, the tile flow is only calculated in the soil layers above and containing tile drains. There is no tile flow if the water table in the field is below the tile pipes. The water table in the field is in the lowest unsaturated soil layer below which all soil layers are saturated. Specifically, the water table in the field is estimated with

$$\begin{array}{}\text{(16)}& D_{\mathrm{1}}=\left(d_{\mathrm{s}}-{\displaystyle \frac{L_i{\mathit{\theta }}_i}{{\mathit{\theta }}_{i,\mathrm{s}}}}\right)\end{array}$$

where D₁ is the water table depth [m]; d_s is the depth to the top of the uppermost saturated soil layer [m]; L_i is the thickness of the lowest unsaturated soil layer [m]; θ_i is the volumetric soil water content of the lowest unsaturated soil layer [m³ m⁻³]; θ_i,s is the saturated volumetric soil water content of the lowest unsaturated soil layer [m³ m⁻³].

Figure 4 Tile fractions over the US Midwest region, and the location and layout of the selected experiment field. The yellow boxes represent sub-fields without tile drainage and the green boxes represent sub-fields with conventional drainage.

2.2 Model setups

2.2.1 Field data

Data from an experimental field site (Fig. 4) in the Iowa State University Southeast Research and Demonstration Farm in Washington County (41.20, −91.49°), were used for the model setup and validation (Chighladze et al., 2021). The major soil types in this site are the Tanitor and Kalona soils (silty clay loam soil) (Fig. S5). The study site consists of four tile drainage treatments: conventional drainage, shallow drainage, controlled drainage, and no drainage. Each of these treatments had two replicates with corn–soybean rotations. Border tiles without monitoring were installed to reduce the interaction between adjacent plots. Only the conventional drainage and no drainage plots were used in this study. Tile drainage was installed in 2006 and the tile flow, crop yield, and daily water table depth were monitored from 2007 to 2017. The tile diameter, tile depth, and tile spacing between tiles drains are 25.4, 122, and 18.3 m, respectively. Management practices, such as tillage and fertilizer application, were documented and used as model inputs. On-site daily precipitation was monitored from 2007 to 2017. The precipitation data in 2007 were removed due to quality issues (Fig. S6). All these data can be accessed on a website of Iowa State University (https://datateam.agron.iastate.edu/td/, last access: May 2023; Chighladze et al., 2021).

2.2.2 Model calibration, validation, and experiment design

Soil properties, weather, management practice data, and tile drainage settings are required to drive the ecosys model. The North American Land Data Assimilation System (NLDAS-2) dataset was used as the major meteorology driver, including temperature, solar radiation, humidity, and wind speed (Xia et al., 2012). Daily precipitation data from on-site observations for the years 2008 to 2017 were substituted for the NLDAS-2 dataset to better capture the local rainfall pattern. As hourly data were unavailable, daily precipitation totals were uniformly distributed over two separate hours on each rainy day. This approach was chosen because common downscaling techniques, such as AWE-GEN (Ivanov et al., 2007), require sub-daily precipitation statistics that were not available for our site. Furthermore, this simplification is justified as our analysis focuses on the system's response at monthly and annual time scales. The soil information was obtained from the Gridded Soil Survey Geographic Database (gSSURGO) dataset, and the detailed soil information used in the study can be found in Fig. S5 and Table S2 (Soil Survey Staff, 2023). The settings for the tile drainage in the tile field are shown in Table 1. We referred to one observation well in the undrained field to set the subsurface boundary condition. Specifically, the external water table depth was set as the mean value of the observed water table depths (1.00 m), and the distance to the external water table was set as 50 m, which is around half of the length of the experimental field. The tile drainage depth was set as the lowest point of the tile drainage, which was the depth of the tile drainage plus the radius of the tile drainage, 1.35 m. The tile drainage spacing was 18.3 m.

Table 1 Tile drainage parameter settings.

The ecosys simulation was started in  1990, with the initial 17 years (1990–2006) as the initialization period to stabilize the model, followed by an 11-year analysis period (2007–2017). Due to the limited length of the observational record, the model was calibrated against the undrained field and validated against the tile-drained field. The configuration of the model relied on established crop parameters in the US Midwest conditions (Li et al., 2022). Key parameters sensitive to local conditions were calibrated using observed crop yields from the undrained field, such as the maturity group and maximum rate of carboxylation (VCMX). The parameter ranges and final values were shown in Table S1 in the Supplement. The Pearson coefficient (r), percent error (PE), root mean square error (RMSE), and relative root mean square error (rRMSE) were used to assess the model's performance, comparing the simulated yield with the observed yield. The tile flow simulation was assessed on a monthly basis by comparing the simulated and observed values through r, RMSE, and the Nash–Sutcliffe model efficiency coefficient (NSE),

$$\begin{array}{}\text{(17)}& {\displaystyle }& {\displaystyle }r={\displaystyle \frac{\mathrm{\Sigma }\left(x_i-\stackrel{\mathrm{‾}}{x}\right)\left(y_i-\stackrel{\mathrm{‾}}{y}\right)}{\sqrt{\mathrm{\Sigma }{\left(x_i-\stackrel{\mathrm{‾}}{x}\right)}^{\mathrm{2}}\mathrm{\Sigma }{\left(y_i-\stackrel{\mathrm{‾}}{y}\right)}^{\mathrm{2}}}}}\text{(18)}& {\displaystyle }& {\displaystyle }\mathrm{PE}={\displaystyle \frac{\stackrel{\mathrm{‾}}{x}-\stackrel{\mathrm{‾}}{y}}{\stackrel{\mathrm{‾}}{x}}}\text{(19)}& {\displaystyle }& {\displaystyle }\mathrm{RMSE}=\sqrt{{\displaystyle \frac{{\mathrm{\Sigma }}_{i=\mathrm{1}}^N\left(x_i-{\left(y_i\right)}^{\mathrm{2}}\right)}{N}}}\text{(20)}& {\displaystyle }& {\displaystyle }\mathrm{rRMSE}=\sqrt{{\displaystyle \frac{\frac{{\mathrm{\Sigma }}_{i=\mathrm{1}}^N{\left(x_i-y_i\right)}^{\mathrm{2}}}{N}}{{\mathrm{\Sigma }}_{i=\mathrm{1}}^N{\left(x_i\right)}^{\mathrm{2}}}}}\text{(21)}& {\displaystyle }& {\displaystyle }\mathrm{NSE}=\mathrm{1}-{\displaystyle \frac{{\mathrm{\Sigma }}_{i=\mathrm{1}}^N{\left(x_i-y_i\right)}^{\mathrm{2}}}{{\mathrm{\Sigma }}_{t=\mathrm{1}}^T{\left(x_i-\stackrel{\mathrm{‾}}{x}\right)}^{\mathrm{2}}}}\end{array}$$

where x_i is the observation, y_i is the simulation, and $\stackrel{\mathrm{‾}}{x}$ and $\stackrel{\mathrm{‾}}{y}$ are the mean value of observation and simulation, respectively. N is the number of observations.

To investigate the effects of tile drainage in wetter conditions, we ran the calibrated model to simulate hydrological processes, biogeochemical dynamics, and crop growth across a spectrum of precipitation scenarios. Specifically, we manually adjusted the precipitation inputs with a scale factor, λ, to evaluate wetter conditions (Fig. S7),

$$\begin{array}{}\text{(22)}& {\displaystyle }& {\displaystyle }\mathit{P}=\left(p_{\mathrm{1}},p_{\mathrm{2}},\mathrm{\dots },p_n\right)\text{(23)}& {\displaystyle }& {\displaystyle }{\mathit{P}}_{\mathit{\lambda }}=\left(\mathit{\lambda }{p}_{\mathrm{1}}\mathit{\lambda }{p}_{\mathrm{2}},\mathrm{\dots },\mathit{\lambda }{p}_n\right)\end{array}$$

where P is the original precipitation time series that drives the model and p_i is the precipitation at ith time step [mm]. P_λ is the precipitation time series with the scale factor $\mathit{\lambda }\in \mathit{\{}\mathrm{0.9},\mathrm{1.0},\mathrm{1.1},\mathrm{1.2},\mathrm{1.3}\mathit{\}}$ corresponding to five hypothetical scenarios. Given that the temperature, which dominates the vapor pressure deficit, is more influential on the crop yield in drought over the US Midwest (Lobell et al., 2014), this hypothetical experiment is mainly targeted at understanding the system's response to flooding with more precipitation.

3 Results

3.1 Model validation

The ecosys model's crop yield simulations were evaluated against field observations (Fig. 5 and Table 2). Overall, the Pearson coefficient r for the yield simulation was over 0.95 for both calibration and validation. Specifically, the PE for corn in calibration and validation was −0.50 % and −3.10 %, respectively. The PE for soybean in calibration and validation was −8.75 % and −1.57 %, respectively. The rRMSE for corn and soybean in validation was 5.66 % and 12.57 %, respectively. Both the observations and simulation showed the benefit of tile drainage for crop yield. For corn, tile drainage increased the yield by 0.830 t ha⁻¹ (6.97 %) and 0.515 t ha⁻¹ (4.20 %) in the model simulation and observations, respectively. For soybean, the tile drainage increased the yield by 0.193 t ha⁻¹ (5.37 %) and 0.498 t ha⁻¹ (13.66 %) in the model simulation and observations, respectively (Table 3). Besides, the model successfully captured the seasonal pattern of more tile flow in the late spring and early summer (Fig. 6). The observations suggested that there was no tile flow in January and February, and we hypothesized that this was due to the low temperature that disabled the measurement device. Thus, we only validated the tile flow in the growing season (April to October). The r, RMSE, and NSE for the monthly tile flow simulation in the growing season were 0.785, 28.42 mm per month, and 0.43, respectively. These results indicate that the ecosys model is promising for estimating both crop yield and tile drainage flow.

Figure 5 Validation for the crop yield and tile flow. Comparison of ecosys-simulated and ground-measured (a) maize (15 % moisture) and soybean (13 % moisture) grain yield, and (b) monthly tile flow in the growing season (April to October). (c) Time series of the simulated and observed yield. Corn years (light gold) and soybean years (light green) are indicated with different background colors.

Table 2 Goodness-of-fit statistics for crop growth for calibration and validation. PE: percent error; rRMSE: relative root mean square error; RMSE: root mean square error.

Table 3 Observed and simulated crop yield under the tile and undrained conditions.

Figure 6 Ecosys-simulated and observed tile flow. The boxes represent 25 %–75 % of the tile flow for the simulated period (2007–2017), and the triangles represent the mean tile flows. The triangles represent the multi-year mean tile flow in a certain month.

3.2 The impacts of tile drainage on hydrology

We first evaluated the water fluxes and partitioning with the ecosys model under the tile and undrained conditions at the study site. The annual mean precipitation was 881 mm from 2007 to 2017, and most of the precipitation (66.4 %) occurred from April to August (Figs. 7 and S8). The simulations indicated that the tile drainage primarily altered the subsurface water exchange (Fig. 7). Compared to the undrained condition, the tile drainage increased the annual mean subsurface discharge from 216 to 276 mm and increased the subsurface recharge from 34 to 104 mm. In contrast, the tile drainage had a limited impact on other fluxes. The annual mean ET was similar under the tile (659 mm) and undrained (655 mm) conditions, as was the surface runoff (48 and 46 mm, respectively). As a direct consequence of the increased subsurface drainage, the simulations showed that the tile drainage reduced the soil water content (Fig. S11a).

Figure 7 Ecosys-simulated annual water balance under the tile and undrained conditions. Overall, tile drainage increases both subsurface discharge (water coming out of the field) and subsurface recharge (water going into the field), and the ET and surface runoff are similar under the tile and undrained conditions. The imbalance between influxes and outfluxes is subject to storage change.

The results indicated that the effects of the tile drainage followed a seasonal pattern (Fig. 6). For the study site, the tile drainage actively removed excess water in the spring, leading to an increase in the net subsurface discharge from 120 to 137 mm, corresponding to high precipitation and low ET during those months. Less water was drained by the tile drainage in the summer due to high crop water consumption despite high precipitation (Figs. 6, 8, and S10). On average, the net subsurface discharge was −53 mm in the summer under tile conditions, indicating a significant recharge from the surrounding soils to the tile-drained field. Tile drainage increased the net subsurface discharge from 17 to 28 mm and from 48 to 60 mm in the autumn and winter, respectively, when compared to the undrained condition. Tile drainage slightly raised the ET in the summer months from 226 to 238 mm, and slightly reduced the ET in other months (Fig. 8).

Figure 8 Ecosys-simulated water fluxes under the tile and undrained conditions from 2007 to 2017. (a) Boxplot of the quarterly net subsurface discharge (subsurface discharge – subsurface recharge); (b) boxplot of the quarterly net subsurface discharge difference between the tile and undrained conditions; (c) boxplot of the quarterly surface runoff; (d) boxplot of the quarterly surface runoff difference between the tile and undrained conditions; (e) boxplot of the quarterly ET; (f) boxplot of the quarterly ET difference between the tile and undrained conditions. The upper and lower parts of the boxplots indicate the 25 % and 75 % quantiles, and the boxes indicate the interquartile variation. The triangles indicate the mean values. Delta is the difference between the tile and undrained conditions. The corresponding monthly results are shown in Fig. S10.

3.3 The impacts of tile drainage on the soil biogeochemistry and subsequent crop growth

In addition to altering water fluxes, the model showed that tile drainage significantly changed the fate and transport of inorganic nitrogen (IN). At the study site, the annual mean total IN loss was 2.72 and 1.89 g Nm⁻² for the tile and undrained conditions, respectively (Fig. 9). Tile drainage primarily increased the subsurface IN loss, with values rising from 1.89 to 2.45 g Nm⁻², while the surface IN loss remained relatively constant, with no significant differences noted between the tile and undrained conditions, around 0.27 g Nm⁻². These values were within the range of riverine nitrogen yield in the central US Midwest reported by David et al. (2010). Over 85 % of IN left the system through the subsurface in both the tile and undrained conditions. Most IN loss happened in the spring, coinciding with fertilizer application and the peak precipitation period. The impacts of tile drainage on IN loss also exhibited a seasonality. Our model revealed that subsurface IN loss had the most substantial increase in the spring with an increase of 0.35 g Nm⁻², while the increase of subsurface IN loss in the summer was only 0.08 g Nm⁻². Furthermore, the model results showed that nitrogen loss increased with the total precipitation in the hypothetical experiments (Fig. S12).

Figure 9 Ecosys-simulated inorganic nitrogen (IN) lost under the tile and undrained conditions from 2007 to 2017. (a) Boxplot of the quarterly and annual subsurface IN discharge; (b) boxplot of the difference of the quarterly and annual subsurface IN loss between the tile and undrained conditions; (c) boxplot of the quarterly and annual surface IN loss; and (d) boxplot of difference of the quarterly and annual surface IN loss between the tile and undrained conditions.

The model suggested that tile drainage increased the soil oxygen concentration (Fig. S11), further affecting the soil biogeochemistry and crop growth. Fig. 10a and b suggest that tile drainage promoted soil microbe activity and accelerated soil organic nitrogen mineralization and soil heterotrophic respiration (Brown et al., 2017; Castellano et al., 2019). Our model showed that the annual mean IN generated in mineralization–immobilization processes increased by 0.36 g Nm⁻², from 9.39 to 9.75 g Nm⁻², with tile drainage. The elevated oxygen concentration also promoted crop root growth and nutrient uptake, as indicated by the increase of mean annual root respiration from 182 to 193 g Cm⁻² (Fig. 10c).

Figure 10 Ecosys-simulated soil root and soil microbe activities from 2007 to 2017. Boxplot of ecosys-simulated (a) annual soil net mineralization (mineralization – immobilization) under the tile and undrained conditions and their difference; (b) annual root respiration under the tile and undrained conditions and their difference; and (c) annual heterotrophic respiration under the tile and undrained conditions and their difference. The upper and lower parts of the boxplots indicate the 25 % and 75 % quantile, and the boxes indicate the interquartile variation. The triangles indicate the mean values. Delta is the difference between the tile and undrained conditions.

To understand the role of tile drainage on crop growth under excessive precipitation, we specifically looked into a typical wet year, 2009. This year was characterized by high annual precipitation, and June in particular received 197 mm of rain, surpassing the multi-year average of 158 mm. Both the simulation and the observation showed that the tile drainage helped to increase the soybean yield during 2009 (Fig. S11). The profile of the soil water content (SWC), O₂, and root density in the soil column on 30 June are presented in Fig. 11. Under the undrained conditions, the soil profile was nearly saturated below a 0.4 m depth. In contrast, the tile-drained profile was not saturated until below 1.1 m, which demonstrated that the tile drainage effectively mitigated excess water accumulation in the soil. Consequently, the soil oxygen concentration was higher under the tile condition (Fig. 11b), which provided aerated conditions for crop root growth. This improved the oxygen availability, contributing to denser and deeper root development (Fig. 11c). Fig. 11e shows the time series of the O₂ stress indicator, defined in Sect. 2.1.3, which showed that the tile drainage helped to reduce the O₂ stress. We observed that the tile promoted crop transpiration, especially during the early growth season (Fig. 11f and g). Similar impacts of tile drainage on the crop yield and soil conditions were observed for corn during another wet year, 2010 (Fig. S14).

Figure 11 Ecosys-simulated soil profile and time series of water stress, oxygen stress, and ET in typical wet and dry soybean years. The profile of (a) soil water content, (b) soil O₂ concentration, and (c) root density profiles in the soil column on 30 June 2009. Time series of (d) minimum canopy water potential, (e) crop actual O₂ uptake rate/O₂ demand (potential O₂ uptake rate under a non-limiting O₂ condition), (f) ET, (g) crop transpiration in the 2009 growing season (a typical wet year for soybeans). The profile of (h) soil water content, (i) soil O₂ concentration, and (j) root density profiles in the soil column on 15 August 2013. Time series of (k) minimum canopy water potential, (l) crop actual O₂ uptake rate/O₂ demand (potential O₂ uptake rate under non-limiting O₂ condition), (m) ET, (n) crop transpiration in the 2013 growing season (a typical year with a wet spring and dry summer for soybean). The vertical red lines indicate the days to get soil profiles. The x axis of (c) and (j) is in log-scale; see Fig. S29 for the plots showing root density in the linear scale.

The results also showed the potential of tile drainage to enhance crop resilience to drought in summer with a wet spring. The annual precipitation in 2013 was 874 mm, slightly below the annual mean value of 881 mm. In 2013, a severe summer drought occurred, with 71 mm of precipitation recorded during the summer months (Fig. S8), and the drought led to a yield drop (Fig. S13). The precipitation in May reached 230.61 mm, which saturated the soil in the early stage of crop growth. Both the observation and model showed an increase in the soybean yield under the tile drainage condition, although our model underestimated the yield (Fig. 5 and Sect. S9). Our results suggested that the tile drainage reduced the soil water content while increasing the soil oxygen concentration and promoting root growth (Fig. 11h–j). Fig. 11k shows that the tile drainage increased the minimum canopy water potential modeled in the summer, indicating that soybeans suffered less water stress under the tile conditions. Our model also suggested that the crop transpired more water during the summer drought under the tile condition, even with less soil water content (Fig. 11h–n), which evidenced that the more developed root system helped crops access the soil water in deeper soil (Fan et al., 2017; Schenk and Jackson, 2005; Steudle, 2001).

3.4 The impact of tile drainage on crop growth under different precipitation amounts

The validated model was then used to assess the impact of tile drainage under various precipitation amounts in hypothetical numerical experiments. Our simulation revealed that the mean crop yield over the assessment period decreased for both corn and soybeans, as precipitation levels increased. The yield reductions with the increase of precipitation were more pronounced under the undrained conditions, and the yield difference between the tile and the undrained conditions became increasingly substantial with rising precipitation levels (Fig. 12). Those findings indicated that tile drainage provided a yield benefit, and this benefit became even more pronounced in conditions of higher precipitation in the sites that already had relatively abundant precipitation. Tile drainage increased the resilience of crops to precipitation increase, indicated by higher crop yield variation with the change of precipitation (Fig. S18).

Figure 12 Ecosys-simulated crop yield in the hypothetical numerical experiment under different precipitations. (a) Multi-year mean soybean yield under the tile and undrained conditions. (b) The soybean yield benefit (yield/grain carbon difference between the tile and undrained conditions). (c) Multi-year mean corn yield under the tile and undrained conditions. (d) The corn yield benefit (yield/grain carbon difference between the tile and undrained conditions). The x axis, rainfall amount ratio, is the scale factor in the hypothetical numerical experiment. For example, 1.3 represents that the precipitation amount at each time step is 1.3 times of the original precipitation.

To examine the impacts of increased precipitation and tile drainage on crop growth, we examined the responses of biogeochemistry and crop growth to varying precipitation levels under the tile and undrained conditions in a typical wet year for soybeans (2009). Fig. 13a shows that the mean soil O₂ content in the top 1 m soil during June decreased with rising precipitation, under both the tile and undrained conditions ($r=-\mathrm{0.943}$ and −0.977 for the tile and undrained conditions, respectively). Higher oxygen concentration further led to lower crop oxygen stress, indicated by elevated values of the O₂ stress indicator under conditions of high soil O₂ (Fig. 13b; r=0.931 and 0.995 for the tile and undrained conditions, respectively). Crops suffering from less oxygen stress tended to develop denser root systems (Fig. 13c; r=0.982 and 0.984 for the tile and undrained conditions, respectively). Our results revealed that the grain carbon reduced as precipitation increased (Fig. 13d). Furthermore, the soil O₂ concentration, crop O₂ stress indicator, root density, and grain carbon under the tile conditions were consistently higher compared with those under the undrained condition (Fig. 13), which indicates that tile drainage would benefit crop growth by elevating the soil O₂ content and then reducing the crop oxygen stress. Besides, the steeper slopes in Fig. 13 under the undrained conditions suggested that the crop system without tile drainage exhibited higher sensitivity to changes in precipitation. This implies that tile drainage could bolster the system's resilience to precipitation variability, and the benefits of tile drainage become more pronounced with more precipitation. Similar results were also shown in the typical wet year for corn (Fig. S19).

Figure 13 Ecosys-simulated responses of the biogeochemistry and crop growth to the precipitation amounts in a typical wet year for soybeans (2009) under the tile and undrained conditions. The relationships between (a) soil O₂ concentration in the top 1 m soil column and precipitation in June; (b) crop O₂ stress indicator and soil O₂ concentration in the top 1 m soil column in June; (c) root density (0.6–0.8 m soil layer) and crop actual O₂ uptake rate/potential O₂ uptake rate under a non-limiting O₂ condition in June; (d) grain carbon and root density (0.6–0.8 m soil layer) in June; and (e) grain carbon and crop O₂ stress indicator in June.

4 Discussion

4.1 On the necessary processes to simulate the tile drainage impacts and the effectiveness of the ecosys model

A central challenge in quantifying the full impact of tile drainage is capturing the complex interactions between hydrology, soil biogeochemistry, and crop growth. In this study, we utilized the ecosys model for its first-principles approach to simulating the coupled feedbacks between the soil water, aeration, biogeochemistry, and crop growth (Grant, 2001), which offers a unique lens to investigate the underlying causal mechanisms that drive the agroecosystem response to drainage. To evaluate the benefits of improved soil aeration, this study relies on the model's mechanistic simulation of the complete soil oxygen cycle, including its transport, phase changes, and biological consumption (Sect. 2.1.2). This provides a direct biophysical connection between a lowered water table and the oxygen available for root function. In ecosys, oxygen stress is an emergent property of the balance between supply and biological demand, while many crop models either ignore this stress or approximate it by applying empirical functions tied to soil water status (Ebrahimi-Mollabashi et al., 2019; Feddes et al., 1978, 2001; Šimůnek and Hopmans, 2009; Vrettas and Fung, 2017). The crop's response to drainage also emerges from a series of simulated feedbacks in the model. For instance, soil oxygen availability regulates root respiration, which provides the energy for root maintenance and growth, and root nutrient uptake (Sect. 2.1.3). Then, the impact on the root system further cascades to the whole crop through processes, such as the coupled plant hydraulics and photosynthesis model over the soil–crop–atmosphere continuum, nitrogen remobilization in the crop root, shoot, and leaf, etc. (Sect. 2.1.1). The effectiveness of this approach is demonstrated in our results, which show how improved soil oxygen under tile drainage alleviates crop oxygen stress, particularly during wet springs, leading to denser root systems (Figs. 11 and 13). By mechanistically linking the hydrological modification to these detailed physiological responses, this modeling approach allows us to dissect the causal chains driving crop yield benefits, providing a more transferable and systematic process-level understanding (Warren et al., 2015).

The application of the model, despite recent advances in modeling capability and process realism, is limited by the availability of observation data. Ecosys relies on hourly weather input to drive the water and energy cycles. However, accurate hourly weather data are not always available. Here we downscaled the daily in-situ precipitation data to provide hourly inputs by assuming even precipitation within two hours in a day. However, precipitation intensity is a key factor that determines the runoff generation mechanisms (i.e., infiltration-excess runoff and saturation-excess runoff) (Horton, 1933; Nanda and Safeeq, 2023; Tromp-van Meerveld and McDonnell, 2006). Our results under different precipitation amounts also show that the proportion of water leaving the system through surface runoff increases with the increase of precipitation amount/intensity (Fig. S20). The simple downscaling method (Sect. 2.2.2) inadvertently reduces the occurrence of intense precipitation events while increasing the frequency of smaller, milder precipitation events, which might underestimate the surface runoff and overestimate the subsurface discharge. Further, to fully leverage the capability of the ecosys and improve its accuracy, a wealth of observational data are necessary for both model calibration and validation. For instance, we did not find work that directly contrasts the root system and soil oxygen status in the tile and undrained conditions, potentially due to the difficulty in measuring the quantities. We suggest that future field and greenhouse experiments prioritize the systematic collections of data on various variables such as water fluxes, solute nutrient fluxes, greenhouse gas emissions, root development, above-ground crop biomass, and more. These datasets would serve a dual purpose: facilitating model validation and performance assessment while deepening our understanding of the underlying physical processes. This improved understanding can then be leveraged to refine model mechanisms and parameterization (Liu et al., 2020; Nóia Júnior et al., 2023; Warren et al., 2015). Expanding the availability of such data would be invaluable in advancing our modeling efforts and increasing their applicability to real-world scenarios. In this case, observation, experiments, and measurement are integrated together, which aligns with the DOE well-proposed model–data experimentation (ModEx) framework (Hoffman et al., 2017).

4.2 The impact of tile drainage on hydrology

Our results demonstrate that tile drainage has a pronounced impact on field hydrological cycles, influencing both water storage and water fluxes (Fig. 14). Overall, tile drainage functions as an efficient conduit for expediting subsurface water drainage (Gramlich et al., 2018; Miller and Lyon, 2021; Pluer et al., 2020; Schilling et al., 2012), directly contributing to an increase in the subsurface discharge (Fig. 7). Concurrently, tile drainage leads to a reduction in the soil water content (Fig. S11). Besides, we found that precipitation alone cannot sustain both the tile flow and field ET, and recharge from an external source is required to close the system water balance at the study site (Fig. S21). Our model indicates that tile drainage increases the subsurface recharge (Fig. 7), and tile drainage has a limited impact on the surface runoff at our study sites (Fig. 7). Though the model predicted less soil water content under the tile conditions, the topsoil layer might be more affected by the soil atmospheric interaction and the soil water content there is similar (Fig. S13). The infiltration rate is largely affected by the soil water content in the top soil layer, which might account for the similar surface runoff in the tile and undrained conditions. Our simplified 1D column representation of the field limits the surface water holding capacity, which might in turn diminish the perceived role of tile drainage in managing the surface water. We also acknowledge that we currently do not consider the macropores that directly connect the tile drainage with the surface soil, which might underestimate the effects of tile drainage on the surface runoff (Askar et al., 2020; Williams et al., 2023). We found that tile drainage does not significantly change the annual ET at the study site, similar to some previous studies (Khand et al., 2017; Yang et al., 2017).

Figure 14 Schematic of the impact of tile drainage on the hydrology, soil biogeochemistry, and crops in the US Midwest agroecosystems with sufficient precipitation in the spring.

The impact of tile drainage on the hydrological cycle exhibits a seasonality, coinciding with the seasonality of climate and crop phenology. Tile drainage actively functions from May to June, corresponding to high precipitation in this period, and tile drainage significantly increases the subsurface discharge (Figs. 6 and 8). In the summer, crops actively draw water from the soil, which reduces the soil water content. High ET, coupled with reduced soil water storage under the tile conditions, results in an increase in the subsurface recharge (Fig. 8). Tile drainage slightly increases the ET during the peak growing seasons, which is balanced by the decrease in the early growing season. In the early growing season, the soil evaporation might be reduced due to the reduction of soil water content under the tile conditions (Yang et al., 2017). In the summer, crop transpiration dominates the total ET (Paul-Limoges et al., 2022; Song et al., 2018), and higher crop productivity under the tile conditions (Fig. 11) would increase the ET in the peak growing season (Beer et al., 2009; Guerrieri et al., 2016; Yang et al., 2017). The hydrological and crop phenological seasonality has informed the design of the controlled drainage and subirrigation system. Specifically, excessive water was removed in wet springs, and retained in dry summers and the non-growing season in controlled drainage (Helmers et al., 2022). With an irrigation system, the removed water can be cycled and irrigate the field crop in summers with less precipitation (Singh et al., 2022).

4.3 The impact of tile drainage on biogeochemistry

The impacts of tile drainage on hydrology further lead to downstream ramifications for the soil biogeochemistry (Fig. 14). Tile drainage has long been recognized as a major contributor to nitrate exporting from agricultural landscapes (David et al., 1997, 2010). The results here similarly suggest that tile drainage increases the subsurface IN loss by 28.5 %, accompanied by a 29.6 % increase of the subsurface discharge in water partitioning at the study site (Figs. 7 and 9). Furthermore, the impact of tile drainage on nitrogen loss exhibits a seasonal variation, mirroring the seasonality in water partitioning, with a more pronounced increase in the early growing season, corresponding to high tile flow (Figs. 6 and 9) (Ma et al., 2023; Williams et al., 2015). The reduction of the soil water content under the tile conditions also leads to an increase in the soil oxygen concentration, which subsequently promotes soil microbe activities (Linn and Doran, 1984), indicated by higher heterotrophic respiration (Fig. 10b). The elevated soil oxygen content also hastens both mineralization and immobilization (Castellano et al., 2019; Randall and Mulla, 2001), ultimately increasing the net mineralization (Fig. 10a). The increased net mineralization compensates for the decreased soil IN through subsurface loss.

4.4 The impact of tile drainage on crop growth

The impacts of tile drainage through hydrology and biogeochemistry on crop growth and yield are intricate and multifaceted. Tile drainage reduces the soil water content, which may limit crop water availability while reducing crop oxygen stress (Fig. 11). The increased IN loss by tile drainage reduces the soil IN content, while the increased net mineralization, contradictorily, increases the soil IN content (Figs. 9 and 10). They, together, might also alter the temporal variation of the soil IN content (Castellano et al., 2019; Drinkwater and Snapp, 2007). Crops and soil microbe communities can also adapt to environmental changes (Fan et al., 2017; Waldrop and Firestone, 2006). As a result, the intricate interplay between tile drainage, hydrology, biogeochemistry, and crop responses can collectively exert a significant influence on crop growth. Also, the intricate interactions are likely to change with environmental variations and make it challenging to gain a full understanding of their impacts on crop growth.

Here we use the root system as a proxy to understand the impact of the tile drainage on crop growth, as crop roots are the key mediator between soil hydrological and biogeochemical changes and crop growth (Hodge, 2004; Hodge et al., 2009; Jochen Schenk, 2005; Wang et al., 2021). In wet springs, tile drainage alleviates crop oxygen stress by reducing the soil water content and elevating the soil oxygen concentration (Figs. 11 and 13), which guarantees the early growth of both the root system and the above-ground part of the crop under excessive precipitation. Deeper and denser root systems are observed under the tile drainage conditions (Fig. 11), which further benefits the crop water and nutrient uptake (Ebrahimi-Mollabashi et al., 2019). Further, we observed that the developed root system also helped to reduce water stress in the dry summer, potentially due to accessibility to water in the deep soil with a developed root system (Fig. 11f–j). The lower oxygen stress and developed root system together benefit crop growth (Fig. 14). The findings are consistent with previous field and greenhouse studies showing shallower root development under waterlogged conditions and high water tables (Ebrahimi-Mollabashi et al., 2019; Follett et al., 1974; Ordóñez et al., 2018; Ren et al., 2016). This modeling study also provides valuable insights into these complex and difficult-to-observe processes and provides a reference for future field investigations.

In summary, our study reveals that tile drainage significantly enhances the production of both corn and soybeans at the study site that has abundant precipitation in the spring. Firstly, tile drainage proves beneficial for crop yield by directly mitigating crop oxygen stress during wet years (Figs. 11 and S14). Secondly, the crop with better developed roots under the tile conditions would have a better ability to absorb soil water and thus reduce crop water stress and benefit the crop yield (Fig. 11). Lastly, our hypothetical numerical experiments indicate that the yield benefit of tile drainage amplifies with increasing precipitation across various amounts (Fig. 12).

4.5 The implications of tile drainage for climate change adaptation

Our results at the study site indicate that tile drainage might be a valuable adaptation strategy to enhance agricultural production under climate change. Our results at the study site have demonstrated that tile drainage has the potential to increase crop yield under excessive precipitation conditions (Figs. 11, 12, and S14), and it might play a more critical role in sustaining high crop yields in the future, especially given the projected increase in spring precipitation and the likelihood of more intense precipitation events in the US Midwest (Seneviratne et al., 2022; Zhou et al., 2022). Furthermore, our results suggest that tile drainage enhances yield stability under different precipitation amounts (Fig. S18), which implies that tile drainage might also help mitigate the risks associated with variable weather conditions, especially the excessive precipitation conditions. Our results also reveal that tile drainage has the potential to sustain a high crop yield under a projected increase in the summer drought due to the better developed root systems under the tile drainage conditions (Fig. 11) (Zhou et al., 2022). While tile drainage benefits crop yield in a severe drought (2013), the crop still faces high water stress, resulting in a relatively low yield (Figs. 12 and S13). We advocate other approaches to sustain high crop yield with increasing summer drought risks, for instance, subirrigation systems and control drainage systems (Singh et al., 2022; Youssef et al., 2023).

However, tile drainage also poses threats to downstream water quality under climate change. Tile drainage increases the nitrogen loss to freshwater systems and, ultimately, coastal regions, degrading downstream and coastal water quality (David et al., 2010; Ma et al., 2023). Under climate change, the increased spring precipitation may flush more nitrogen through the tile drainage (Fig. S12), further burdening impaired water systems (Jiang et al., 2020; Sinha et al., 2017; Wang et al., 2015). Recent initiatives are focused on mitigating nitrate loading in tile-drained systems through within-field management practices (e.g., improved fertilizer management and cover crops) and edge-of-field practices (e.g., controlled drainage, saturated buffers, woodchip bioreactors, and subirrigation systems) (Mitchell et al., 2023; USDA NRCS, 2017, 2023). While many studies suggest promising outcomes of these conservation practices in terms of reducing nitrogen loss and enhancing other ecosystem services, debates persist on their effectiveness for controlling the nutrient loss under different environments and socioeconomic feasibility for a broad adoption (Frankenberger et al., 2023; Mitchell et al., 2023). For instance, the controlled drainage, involving a water control structure at the tile drainage system outlet, holds water in the field when drainage is unnecessary, which may help reduce the nitrogen loss and potentially provide yield benefits under dry conditions (Delbecq et al., 2012; Ghane et al., 2012; Singh and Nelson, 2021; Youssef et al., 2023). Controlled drainage could potentially enhance yield with more available water in such cases. Nevertheless, an existing study shows that controlled drainage might have negative impacts on yield during wet seasons (Youssef et al., 2023). Moreover, the efficacy of controlled drainage in reducing nitrogen loads remains highly uncertain (Mitchell et al., 2023; Ross et al., 2016; Shedekar et al., 2021). While controlled drainage directly reduces nitrogen loading by retaining water, uncertainties arise as the retained water and nitrogen may exit the system through other pathways, such as surface/subsurface runoff, adjacent tile systems, or deep percolation (Lavaire et al., 2017; Ross et al., 2016; Shedekar et al., 2021). Furthermore, higher financial costs for the control structure installation also prevent its adoption by farmers and landowners. Similar issues of high uncertainty and additional financial costs are faced by other practices aiming to reduce nitrogen loading. Consequently, more research is needed to comprehensively understand the impacts of conservation practices on agricultural productivity, nutrient loss reduction, and other ecosystem services as well as their tradeoffs and balance in the context of climate change mitigation and adaptation in tile-drained agricultural ecosystems.

5 Conclusions

In this study, we used a process-based model to evaluate the impact of tile drainage on the coupled hydrology, biogeochemistry, and crop system. Our results reveal that soil oxygen dynamics are the key emergent property and central mediator that link the physical alterations in hydrology to the profound responses in biogeochemical cycles and crop growth in tile-drained systems. Specifically, crop oxygen stress in the root system emerges as the competition between the crop oxygen demand and the soil oxygen supply, governed by oxygen transport, phase changes, and biological consumption mechanisms in the model. Tile drainage influences the system firstly by altering the water fluxes and storage, which subsequently changes the soil aeration and further regulates the crop oxygen stress. Its impacts on the root zone further cascade to the whole crop system through the tightly coupled soil–crop–atmosphere processes. Meanwhile, the biological oxygen consumption also drives the change in the soil biogeochemistry processes and further affects crop growth. The model performance is validated with field crop yield data and tile flow observation. We systematically compared model simulations under both the tile-drained and undrained conditions to quantitatively evaluate the influence of tile drainage on the hydrology, biogeochemistry, and crop growth of the agroecosystem. Further, through a series of numerical experiments, we revealed the pivotal role of tile drainage in the face of climate change, considering various precipitation scenarios.

• The impact of tile drainage on hydrology: We found that tile drainage firstly modifies the hydrological cycles, influencing both the water storage and water fluxes. At the study site, our results reveal that tile drainage reduces the soil water content, and increases the annual subsurface discharge and subsurface recharge, while it does not significantly change the surface runoff and ET. Those impacts on hydrology exhibit a seasonality, controlled by the seasonality of the climate and crop phenology. Specifically, tile drainage mainly increases the subsurface discharge when there is high precipitation or low ET and increases the subsurface recharge when the crops actively transpire water from the soil.

• The impact of tile drainage on soil biogeochemistry: The changes in hydrology further propagate through the agroecosystem, instigating ramifications within the biogeochemical cycles. Specifically, tile drainage increases the subsurface nitrogen leaching with the increase of subsurface discharge. Tile drainage also elevates the soil oxygen content, as fewer soil pores are occupied by water. The elevated soil oxygen content further increases the soil net mineralization.

• The impact of tile drainage on crop growth and its implications under climate change: Those changes in the hydrology and biogeochemistry substantially benefit crop growth under both wet springs and dry summers. High soil oxygen concentration under the tile-drained conditions provides an aeration condition that mitigates crop oxygen stress, promoting robust root development and overall crop growth in wet springs. The developed root system also enhances crop resilience to summer drought. We also found that the yield benefit of tile drainage increases with the increase of precipitation and higher crop resilience to precipitation variation under the tile drainage conditions.

In conclusion, this study reveals the underlying causal mechanisms that drive the agroecosystem response to drainage on hydrology, biogeochemistry, crop growth, and their interconnections with oxygen dynamics as the key mediator. Our study also highlights tile drainage as a promising and adaptable climate change adaptation and mitigation strategy with the potential to enhance agricultural resilience in the US Midwest agroecosystems.