Modeling the Effects of Artificial Drainage on Agriculture‐Dominated Watersheds Using a Fully Distributed Integrated Hydrology Model

In agriculture‐dominated watersheds where natural drainage is poor, agricultural ditches (narrow engineered channels) and tile drains (perforated pipes) are widely employed to enhance surface and subsurface drainage, respectively. Despite their relatively small scale, these features exert substantial control over the hydro‐biogeochemical function of watersheds and their effects need to be represented in the models. We introduce a novel strategy to incorporate the effects of artificial agricultural drainage into a fully distributed basin‐scale integrated surface‐subsurface hydrology models. In our approach, narrow agriculture ditches for surface drainage are resolved efficiently using ditch‐aligned computational meshes that are hydrologically conditioned to ensure connectivity in the stream/ditch network. For tile drainage in the subsurface, we use the physically based Hooghoudt's drainage equation as a subgrid model and route the water drained through tiles to the nearest ditch. Without site‐specific calibration, this model reproduced observed streamflow in the Portage River Watershed (>1,000 km2) as recorded by a USGS gauge with good accuracy (normalized KGE = 0.81) and outperformed a calibrated SWAT model (normalized KGE = 0.68). Numerical experiments confirm that artificial drainage reduces surface inundations and effectively controls the water table. At the watershed scale, artificial drainage increases baseflow but has little effect on watershed discharges above the 90th percentile. The strong physical underpinnings and reduced need for calibration allow us to study the impacts of artificial drainage on distributed hydrological response in terms of fluxes and states and provide a platform for investigating watershed‐scale nutrient transport.


Introduction
There is an ever-growing concern about the environmental impacts of agriculture on natural water bodies (Metcalfe et al., 2017;Rasul, 2016).Excess inputs of nitrogen (N) and phosphorous (P) from fertilizer application, artificial drainage, and livestock grazing in agricultural watersheds threaten critical drinking water sources and aquatic ecosystems through eutrophication and hypoxia (Daniel et al., 1998;Lintern et al., 2018Lintern et al., , 2020;;Selman et al., 2008).Artificial drainage via tile drains and agricultural ditches contributes significantly to the mobilization and delivery of N and P to receiving waters by creating a network of fast-flowing hydrologic pathways (Blann et al., 2009;Sun et al., 2021).To control nutrient inputs and enhance water quality, it is essential to understand the intricacies of hydrological and transport pathways (King et al., 2015;Macrae et al., 2019;Rallapalli et al., 2022).The success of strategies aimed at improving water quality relies on understanding local hydrological conditions and complex spatial and temporal patterns of non-point sources (Macrae et al., 2019).This necessitates the use of models that accurately depict hydrological connectivity, including tiles and ditches, as well as nutrient transport processes at varying spatial scales (Kay et al., 2009;Lintern et al., 2018;Rallapalli et al., 2022).In this study, our primary objective is to develop a basin-scale, fully distributed hydrology model capable of resolving heterogeneous impacts of tiles and ditches on the hydrological processes with the potential to provide valuable insights into water quality management.
Agricultural watersheds present unique challenges to fully distributed models because of ubiquitous small-scale yet hydrologically dominant drainage features.Tile drains and agricultural ditches (henceforth referred to as tiles and ditches for brevity) are used extensively to manage agricultural water and improve crop productivity in regions with shallow water tables and poor natural drainage (flat topography and low soil permeability) like the fertile glaciated Midwestern United States ("Midwest").Tiles are perforated pipes laid down typically 0.5-2 m below the surface that allow the egress of the excess water from the root zone (Schilfgaarde, 1999).Through tiles, water leaving the topsoil drains into open ditches and streams.Ditches are narrow channels constructed around the farm fields to channel water from the tiles and surface into nearby streams.Natural stream systems are often straightened and enlarged to drain the incoming water from tiles and ditches (Blann et al., 2009).Interactions between tiles and ditches are known to amplify the impact above their individual contributions (Blann et al., 2009;Schwab et al., 1963;Skaggs et al., 1994).Tile drains shift a major component of flow from surface runoff to the subsurface, which tends to reduce sediment and phosphorous export to ditches and streams as sediment and associated phosphorous concentrations are lower in the tile-drained water compared to surface runoff.However, tile-drained waters are also more likely to bypass buffer strips and other features installed to trap nutrients in the surface runoff (Blann et al., 2009;Haygarth et al., 1999;Sims et al., 1998).Conversely, tile-drained water typically carries higher loads of soluble contaminants like nitrates and dissolved phosphorous (Blann et al., 2009; accordance with the DOE Public Access Plan: (http://energy.gov/downloads/doepublic-access-plan) Water Resources Research 10.1029/2023WR035993 King et al., 2014).These interactions and competing effects highlight the need to represent the effects of both tile drains and surface ditches in models.
Modeling tiles in 3D fully distributed hydrological models is challenging, specifically because the exact locations of tile pipes are generally unknown (De Schepper et al., 2015;Hansen et al., 2013;Schilling et al., 2019).Recently, thermal imaging has been shown to be a promising tool for mapping tiles (Allred et al., 2018;Woo et al., 2019), however, these experiments have not been scaled to basin scales.Even in the cases when the information on tile locations may be available, modeling individual tile pipes requires high refinement of mesh around tiles for accurate simulations, which is not practical in basin scale models due to the required mesh size and computation burden (De Schepper et al., 2015, 2017).Several approaches have been adopted in previous studies to represent tile effects implicitly at the catchment scales.A widely used approach to represent the effect of tile drains in fully distributed models is to use internal sinks or seepage nodes at the tile elevation.Details of that approach vary, depending on how the seepage nodes are implemented, whether they correspond to actual tile locations or are distributed across the entire drained areas, and how the tile efflux is handled.CATHY (Muma et al., 2017) andFLUSH (Nousiainen et al., 2015;Turunen et al., 2013;Warsta et al., 2013) use constant pressure head internal boundary condition at tile nodes to compute the tile drain flow and simply remove the drained water without placing it back in the surface flow system.MIKE-SHE (Hansen et al., 2013) uses a head-dependent seepage condition and computes tile-drain flow based on the water table above the drain level and a userdefined characteristic time constant, with the option of routing the drain efflux to the surface flow system.HGS applications range from field-scale studies with explicit representation of flow in tile piles using 1D linear elements (De Schepper et al., 2015;Thomas et al., 2016) to catchment-scale studies with the seepage node approach (Boico et al., 2022;De Schepper et al., 2017).A more approximate approach is to deploy a highpermeability layer, also known as an equivalent medium, to mimic tile effects through a lateral preferential flow (Carlier et al., 2007;De Schepper et al., 2015;Rozemeijer et al., 2010;Thomas et al., 2016).The highpermeability layer approach has the disadvantage of mixing tile flow and non-tile flow, which limits its utility in applications that involve solute transport, where it is important to separate transport paths that go through tiles from those that do not.In a recent modeling study of an agricultural catchment in Denmark, Boico et al. (2022) found that multiple conceptual models performed similarly at reproducing observed streamflow.However, nutrient transport is more sensitive to partitioning between surface and subsurface flow paths, which highlights the importance of realistically representing the hydrological processes that govern partitioning among competing flow paths.
Previous efforts studying the effects of drainage in spatially resolved models have largely neglected agricultural ditches despite their strong hydrologic controls and interactions with the tile system.Explicitly resolving narrow agricultural ditches is challenging because of the computational cost of the required mesh resolution.Standard triangulated irregular networks (TIN) or structured raster-based meshes would require small mesh elements, making watershed-scale models computationally infeasible.Ditches are often not fully resolved in the DEMs on which ISSHMs rely to route the overland flow.In addition, DEMs often show other artifacts like false discontinuities in ditches and streams where they pass under highway and railroad embankments.Rallapalli et al. (2022) highlighted the importance of accurate simulation of hydrological connectivity for planning effective BMPs.A few studies attempted to overcome these challenges through effective parameterizations recreating effects of surface drainage in catchment-scale models, like using an effective Strickler coefficient in pseudo-ditch approach (Dunn & Mackay, 1996), or internal sinks (Nousiainen et al., 2015;Warsta et al., 2013).Such approaches undermine the physical basis of the models as effective parameters rely on calibrations and may not reflect how conveyance through ditches depends on dynamic local conditions.Therefore, there is a need for a new approach that can resolve the heterogeneous effects of ditches explicitly as narrow land-surface depressions in the watershed-scale models.
Here we explore high-resolution ISSHMs as tools for modeling agriculture-impacted river basins.Our objective is to extend ISSHM models of agriculture-impacted regions from the scale of small catchments to river basin scales while retaining sufficient spatial resolution to represent the landscape.We address the challenge of representing managed drainage at scale through a meshing strategy that combines ditch-and stream-aligned mixedelement meshes with hydrologic conditioning of the ditch/stream network.That surface drainage representation is combined with a nonlocal, implicit representation of drainage through subsurface tile drains.The implicit representation is similar to subgrid drainage models used previously in semi-distributed models (Guo et al., 2018;Valayamkunnath et al., 2022).Here we adapt the subgrid drainage representation to fully distributed models, add routing of the efflux to the surface flow system, and implement the model in the parallel Advanced Terrestrial Simulator (ATS) software.Using river discharge in the heavily managed Portage River Basin, Ohio as an example, we evaluate the performance of ATS with and without our new representations of managed drainage and compare it to the calibrated semi-distributed model SWAT (Neitsch et al., 2011).Thereafter, we perform highresolution numerical experiments to delineate the effects of surface and subsurface drainage on key hydrological fluxes and states at river basin scales (>1,000 km 2 ).

Field Site
As a test case, we selected a watershed in the Portage River basin (see Figure 1), which is about 1,024 km 2 containing 12 HUC12 units located in the midwestern US in an area suitable for agriculture due to highly productive soil (Mehaffey et al., 2012) and abundant rainfall (Williams & King, 2020).For convenience, we will call it Portage Watershed.A USGS gauge station at the outlet of the watershed near Woodville (OH) provides discharge data for model evaluation.Like most of the watersheds in this region, the Portage Watershed is dominated by agriculture supported by nutrient-rich glaciated soil.Landscapes in this region are poorly drained due to flat topography and low soil permeability, severely limiting the water egress from the surface and subsurface, respectively.Networks of agricultural ditches and tile drains are ubiquitous in the region to drain the landscape efficiently and make agriculture feasible.This fast-flow system dominates watershed hydrology for these agricultural landscapes and their effects must be represented if ISSHMs are to accurately reflect watershed hydrology.For the 2010-2019 period, the mean annual precipitation (Daymet, Thornton et al., 2020), evapotranspiration (SSEBop, Senay, 2018), and discharge (USGS gage) in the Portage watershed are 1035.08,410.74, and 441.90 mm, respectively.

Integrated Surface-Subsurface Hydrology Simulator
We use the Advanced Terrestrial Simulator (ATS) (Coon et al., 2019) for fully distributed integrated surface/ subsurface hydrology simulations.ATS is an open-source code offering a flexible multiphysics framework to solve complex ecosystem-based, integrated, distributed hydrology (Coon et al., 2019(Coon et al., , 2020)).In this study, we solve for water fluxes and storage changes in a 3D domain, including surface, subsurface, canopy, and snowpack.Subsurface flow is represented by the 3D Richards equation, and overland flow is modeled with the diffusive wave equation.Surface and subsurface flow are coupled through flux and pressure continuity conditions (Coon  , 2020;Kollet et al., 2017).Additionally, snow and canopy water dynamics at each surface grid cell are represented by simple ordinary differential equation-based bucket models.The aboveground water fluxes are solved using a land surface process model that includes snow interception, accumulation, sublimation, and melt for the snowpack domain, interception and transpiration for the canopy domain, and bare ground evaporation.These models for aboveground fluxes are based on a multichannel variant of the Priestley-Taylor model for potential evaporation and single-layer radiation balance.For more details and the mathematical formulation of these evapotranspiration models, see Bhanja et al. (2023).ATS builds upon the Arcos multiphysics framework (Coon et al., 2016) to allow the rapid implementation of new processes through composable "evaluators" representing fluxes, sources, or other terms in the water conservation equations.To solve these equations, ATS uses the computational infrastructure of Amanzi that includes an unstructured-mesh library, data structures, discretization schemes, and solvers (Moulton et al., 2012) along with its third-party libraries like Trilinos for parallel data structures and linear solvers (Heroux et al., 2003).Bhanja et al. (2023) recently evaluated this ATS model configuration for integrated hydrology simulations over various natural watersheds in the eastern and central USA.They found a good match of ATS-simulated streamflow and evapotranspiration with the gauge observation and MODIS-derived evapotranspiration, respectively, without any catchment-specific calibration.In this study, we take a similar approach and use only public data products and ATS defaults to the extent possible with the goal of an ISSHM for agricultural watersheds that minimizes or eliminates the need for calibration so we can perform numerical experiments to understand how watershed system will respond to changing land use and climatic conditions.

Stream/Ditch-Aligned Mixed-Element Mesh Resolving Agricultural Ditches
Agricultural ditches are narrow channels constructed around farm fields to collect water from tile-drain outflows and surface runoff.To resolve them at the watershed scale, we adopt a stream/ditch-aligned mixed-element meshing strategy, as shown in Figure 3a.Under this approach, we place elongated quadrilateral elements along the mapped streams/ditches.The short edges of the quadrilaterals are determined by the width of the stream/ditch, and the length of the long sides are determined by a user-defined discretization scale.To avoid small triangles at stream junctions, the first element downstream of a stream junction is modified to become a pentagon.The first elements in headwater streams are pinched out to become an elongated triangle.This strategy avoids exposing short edges to the non-stream regions so that once the stream network is meshed, a TIN mesh coincident with the stream mesh can be constructed for the rest of the surface without generating excessively small triangles.The TIN mesh is then typically coarsening away from the streams.We demonstrate (Rathore et al., 2024) the effectiveness of this approach in resolving the lowest parts of the river valley (river corridor) with an accuracy comparable to locally highly refined TIN mesh but at a fraction of the computation cost and improved flexibility to include additional river-specific hydrologic details.In the context of agricultural watersheds in this study, the approach of mixed-element mesh allows us to resolve narrow ditches efficiently, enforce hydrologic connectivity in stream/ditch networks, and route tiles water into ditches.
Hydrological connectivity and water conveyance in our stream/ditch-aligned mesh are challenged by multiple factors specific to agricultural watersheds.Some examples include: -Agricultural ditches are mapped in the National Hydrography Database (NHD Plus High Resolution (HR) Flowlines), but their mapped locations are often slightly displaced from their actual locations (Figure 2a) and inconsistent with high-resolution DEMs, leading to incorrect elevations in mesh elements in the stream/ditch.-Discontinuities in ditch elevations can occur due to the stream/ditch passing through culverts and underpasses beneath road/rail embankments (Figure 2b) or due to resolution mismatch between DEM and narrow ditches (Figure 2c) -Close proximity of the stream/ditch to detention ponds or wetlands can create spurious connections between ditches and those features (Figure 2d) To ensure hydrological connectivity, a multistep conditioning method for the stream/ditch mesh is adopted.Ditch elements are placed on NHD flowlines, assigned an elevation by smoothing the DEM in the neighborhood of the element, then displaced downward by a user-defined amount that is consistent with local agriculture practices.Hydrologic conditioning is then used to remove obstructions to flow in the stream/ditch network that can arise from the data artifacts mentioned above.Finally, the bank integrity of the ditches in proximity to ponds and pits is checked and restored by upward displacement if needed.This conditioning ensures water is routed in the stream/ ditch network and is consistent with regular maintenance practice involving channelization, dredging, and cleaning of these drainage networks to ensure uninterrupted conveyance of water.See Appendix A for more details on conditioning.
Figure 3 shows the simulation-ready mixed-element mesh for the Portage watershed.As can be seen in Figures 3c  and 3e, the narrow ditches are effectively resolved by quad elements.Clearly a very large number of small triangles would be needed to resolve these ditches, leading to large computational burden.The ditch cells delineated as quad elements are also instrumental in implicit tile modeling, as described in the following subsection.To delineate the tile and ditch effects, we performed numerical experiments with two additional model configurations: (a) without tile drains (including only surface ditches), (b) without tile drains and ditches (natural stream network only, no drainage) (Figure 3b).In the no-drainage case, the reaches from the NHDFlowline layer of the NHD(HR) data set that had "None" as the GNIS_Name attribute were assumed to be ditches and were not resolved using quad-meshing (Figure 3b); the resulting disconnect pits were removed during the pit-filling step.This semi-automatic workflow for generating mixed-element mesh is implemented and released in Watershed Workflow v1.5.

Subgrid Tile-Drainage Model
We adopt a subgrid tile-drain modeling approach that assigns a distributed sink to the catchment associated with each reach in the stream network.Each catchment-based sink comprises multiple cells in the subsurface layer, placed at the depth of 1 m from the surface uniformly throughout the domain (Figure 4).Tile drain fluxes, as described below, are established from each catchment-based sink to the associated stream or ditch that will receive the drain effluent.The catchment delineation for each reach is taken from NHDPlus(HR) data sets which are determined by the DEM (Figure 4b), and, therefore, may or may not be the actual connection between the land surface and reach.Whether the catchment is tile-drained or not is determined based on the data set by Valayamkunnath et al. (2020), which, on a 30 m resolution raster, maps if the pixel is likely tile-drained or not, as pixel values of 1 and 0, respectively.From this data set, if the averaged pixel values within the catchment polygon are greater than 0.25, we marked that catchment as tile drained.
The water flux associated with the tile drainage from the subsurface domain can be calculated from Hooghoudt's drainage model (Ritzema, 1994;van Schilfgaarde et al., 1957).The assumptions for the derivation of this model include homogeneous soil, horizontal flow (vertical equipotential lines), uniform drain spacing, steady-state flow, and impermeable layer underlying the drains.For equally spaced drains when the water table elevation midway between drains is above the drain elevation but below the ground surface elevation, where q tile [m/s] is the steady state volumetric discharge per unit area, K b and K a are the hydraulic conductivities below and above the drain, respectively, d is the equivalent distance from the drain to the bottom impermeable boundary, h m is the water table elevation relative to drain elevation midway between drains, and L is the drain spacing.
where r 0 is the radius of the tile drain, γ = 2πD L is a nondimensional parameter relating the actual distance to the bottom impermeable boundary D to the drain spacing L. The F(γ) term in the denominator is expressed as a series (van der Molen & Wesseling, 1991) that converges rapidly for γ > 0.5, which is met in most applications, and an approximate closed form expression for γ ≤ 0.5: The intent here is to implement q tile as a nonlinear drain without resolving individual tiles.At this coarser resolution, we have access to the average water table elevation in each grid cell but not the subgrid-scale quantity h m .
Assuming the water table elevation between drains varies parabolically from 0 at the drain elevation to h m midway between drains, h m is 3/2 the average height of the water table above the drains h ave and Equation 1becomes Rewriting in terms of the average water pressure (p) within a cell instead of head (h ave ) to enable implementation in ATS, we have Where and permeability k a and k b are permeability above and below the tiles, μ is the viscosity of water [Pa s], and p is pressure for a particular cell in the tile subsurface layer.
The total sink Q tile [mol/s] from a given catchment is then obtained by integrating over the tile layer in the catchment where η is the molar density of water [mol/m 3 ], and H is the Heaviside function.
Water collected from all the tile cells in the subsurface domain of a catchment is instantaneously introduced in the corresponding stream/ditch reach cells as a uniformly distributed source.This approach allows tile drainage to be represented at the watershed scale without modeling individual tile pipes.Similar subgrid models based on Hooghoudt's drainage equation have been implemented in semi-distributed models (Valayamkunnath et al., 2022), but as far as we are aware, this is the first implementation within a fully distributed model.
Implementation in a fully distributed model requires aggregating the drainage across each catchment, then distributing that aggregated flow into the surface quad cells of the associated ditch or stream reach.While conceptually straightforward, efficient implementation in a parallel simulation based on domain decomposition requires some care.In general, domain decomposition will not be consistent with the mapped drain catchments.Drain cells and receiving ditch cells for a given catchment/ditch pair may reside on different CPUs.To avoid the bookkeeping overhead, each core keeps information about all the accumulated tile sinks.The source-sink Water Resources Research 10.1029/2023WR035993 mapping between the surface quad-cells of a reach or ditch and subsurface tile cells in the corresponding catchment is provided through a file containing preprocessed global cell IDs.In the case of the Portage Watershed, we have a total of 492 catchment-reach pairs.This source-sink mapping is implemented as an Arcos evaluator and provided to the Richards (subsurface) and overland flow equations to be aggregated with all other sources and sinks of water.

Model Input Data Sets
Various publicly available data sets hosted by different federal and state government agencies were used to define spatially distributed model inputs and evaluate model performance.Globally defined model parameters like parameters appearing in the evapotranspiration models were set to the ATS default values as in Bhanja et al. (2023), which were chosen based on literature values and some limited model sensitivity exploration undertaken prior to the Bhanja et al. ( 2023) study.See Supporting Information S1 for tables listing data sources and values of key model parameters.Manning's coefficients for ditches/streams and the rest of the domain (predominantly croplands) were set to 0.02 and 0.025, respectively.These stream/ditch-specific details could be conveniently assigned to the mesh because of our use of a stream-aligned mixed-element mesh.From the Daymet repository (Thornton et al., 2020), daily maximum and minimum temperature, precipitation, incoming solar radiation, vapor pressure, snow water equivalent, and day length data at 1 × 1 km 2 spatial resolution were used.Our study includes an additional data set by Valayamkunnath et al. (2020) to determine if the catchment is likely tile drained or not.These data sets were retrieved, processed, and mapped onto computational mesh by Watershed Workflow (Coon & Shuai, 2022).Figure 5 shows a few example data sets mapped onto the computational mesh for the Portage Watershed for spatially resolved integrated hydrology simulations.Tile spacing and tile radius in Hooghoudt's drainage model were set to 15 m and 5 cm, respectively (Miller & Lyon, 2021).These tile-drain parameter values fall within the mid-range of typical configurations in Ohio watersheds (Miller & Lyon, 2021).Additionally, sensitivity analysis also indicated that tile spacing of 15 m effectively constrained water table depth below the tile depth in the Portage Watershed.The depth to the bottom impermeable layer and permeability values in Hooghoudt's model were assumed to be uniform and were computed by taking the average over the domain.

SWAT Model for Portage Watershed
For comparison purposes, we also performed SWAT simulations for the Portage Watershed.SWAT is a semidistributed model that relies on a watershed-specific calibration (Neitsch et al., 2011).Here, we used the version of SWAT that does not solve explicitly solve groundwater flow but instead deploys linear storage functions for shallow and depth groundwater flow paths.For our case, we created 1529 hydrologic response units (HRUs) using slope class, land class and soil class with a 10% threshold value.The SWAT model was driven by Daymet meteorological forcings.Usually, the tile is installed at the cropland with poor drainage soil and flat slope (Daggupati et al., 2015); we assumed that HRUs meeting three conditions have tile drains: (a) a slope less than 1.4%, (b) land use class as soybean or corn, and (c) soil hydrologic group as C or D. Out of 1529 HRUs, 673 HRUs were identified to have tile drainage.We selected the CN-based surface runoff option and the Penman-Monteith evapotranspiration, which is downregulated based on moisture availability.
The lateral flow component is computed using a kinematic storage model and shallow groundwater uses the available storage and drainage coefficient.The deep groundwater was assumed to bypass the watershed.The flow in the stream is routed using a variable storage model and tile drain flow is based on soil moisture, depth to the water table, and time lag between drainable water and tile pipes.Appendix B provides details of the SWAT model setup for the Portage Watershed, including a comparison of HRU maps with and without tile drainage criteria.
For the period of WY 2010-2014, we calibrated 24 SWAT parameters (see SI for the list of calibrated parameters) using a dynamically dimensioned search (DDS) (Mudunuru et al., 2022;Tolson & Shoemaker, 2007), a globally optimized scheme, by comparing the modeled and measured streamflow.Further details on SWAT calibration methodology can be found in SI.This study used a total of 5000 model runs with 10 different random seeds and 500 model runs per random seed.Among 5000 model parameter sets, we chose the behavioral parameter sets that show higher streamflow accuracy; for example, the Nash-Sutcliffe efficiency (NSE) and log NSE (NSE value of the log-transformed flow) are greater than 0.55, and mKGE is greater than 0.6, and percent bias (PBIAS) is less than 10%.The selected top 7 parameter sets were used to evaluate the modeled streamflow in the validation periods of WY 2015-2020.

Model Evaluation
We consider results for water years (WY) 2011-2019 for model evaluation and comparison.To get the initial conditions right for the transient simulation, a three-step spin-up process consisting of steady-state, cyclic steadystate, and transient spin-up steps, is performed.Steady-state spin-up involves 20,000 days of simulation using constant precipitation (60% of the mean precipitation).In cyclic steady-state spin-up, a typical year of forcing, computed as a day-by-day average of the full record, is repeated for 15 years.The transient spin-up step involves 6 years of simulation using real data (2005-2010), which is not included in the model evaluation.
The models were evaluated by comparing daily continuous streamflow data from 2011 to 2019 from the USGS archive (https://waterdata.usgs.gov)using Kling-Gupta efficiency (KGE), a multi-criteria model performance metric (Gupta et al., 2009;Kling et al., 2012) developed for hydrology applications.Knoben et al. (2019) noted that KGE for the case when the mean of the observation is used as a predictor, a standard threshold benchmark for model performance, becomes 0.41, which makes its interpretation less intuitive.We use a modified KGE (mKGE, see Equation 7) as in Bhanja et al. (2023), where an mKGE score of 0 indicates that the model has the same explanatory power as the mean of the observations, and a score of 1 indicates a perfect match, analogous to Nash-Sutcliffe Efficiency.
To further evaluate our ATS model, we also compared our results to the top 5 cm soil moisture levels from the SMAP-HydroBlocks (SMAP-HB) product (Vergopolan et al., 2020(Vergopolan et al., , 2021)).SMAP-HB integrates NASA's SMAP-L3E data (9 km resolution) (Chan et al., 2018) with advanced land surface modeling, radiative transfer modeling, machine learning, and ground-based observations to achieve hyper-resolution.SMAP-HB is a model result that assimilates satellite observations and comes with its own uncertainties.Thus, the comparison is best Water Resources Research 10.1029/2023WR035993 thought of as an intermodel comparison not as a comparison to observations.The comparison is further complicated by the fact that SMAP-HB and ATS results are both available only as snapshots that are not coincident in time.For that reason, we focus on qualitative comparisons of spatial patterns and quantitative comparisons of annual mean values of the distributions in time.
Intermodel comparisons seek to establish robustness of a model result (e.g., Baumberger et al., 2017), an important line of evidence to build support for a process-based model.

Model Evaluation
Model spinups and production runs for parallel ATS simulation took 292 hr on 748 cores on NERSC-Perlmutter HPC.SWAT used 56 cores with two threads per core for 5007 runs (5000 calibration runs and 7 validation runs).
With 5 min per SWAT run, the total run time of the simulation campaign is 417.25 hr.Based on mKGE and NSE values (Table 1), we find that uncalibrated ATS outperformed the calibrated SWAT model for the Portage Watershed.Figures 6a-6d compare ATS simulated streamflow and observed at the USGA Woodville gage.Despite underestimating some peaks (particularly in the spring season), ATS performed well at capturing most peaks and their timings (Figure 6a).Base flow is also generally reproduced well (Figures 6b and 6d).Based on the performance rating by Moriasi et al. (2007), ATS's performance can be classified as good and SWAT's as unsatisfactory.
The ATS model particularly benefited from our detailed effort at ensuring hydrological connectivity across the managed landscape of the Portage Watershed.Resolving agricultural ditches and tile drains that exert strong controls on the hydrology of the Portage watershed, characterized by flat topography and low permeability soils, allowed water conveyance as expected.A snapshot of the ponded depth field in Figure 7 shows that ditches drain the peripheral area and route the water effectively.
In Figure 8, we present a comparison of the top 5 cm soil moisture levels from the ATS and SMAP-HydroBlocks (SMAP-HB) products.To enable a direct comparison, we projected the ATS results onto a 30 m × 30 m SMAP-HB grid for the year 2019.Developed areas identified by the NLCD data set (urban areas and roads) were masked because a representation of impermeable surfaces was not included in the ATS model.The median spatial difference (ATS minus SMAP-HB) was observed to be 0.04 (Figure 8a), likely due to sampling discrepancies and differences in modeling assumptions.The SMAP-HB soil moisture data were adjusted to remove the bias to facilitate visual comparisons of the spatial patterns.The annual mean comparison (Figures 8b and 8c) reveals a good resemblance in fine-scale spatial patterns of soil moisture attributable to local topography, soil characteristics, and land use.In particular, both maps have sparse northeast-trending bands of drier soils north of latitude

Effects of Tile Drainage on Subsurface Water
To isolate the tile effects, we reran the model with the tile model turned off.Note that in the no-tiles case, we still have ditches represented through quad-elements in the mesh.Figures 9a and 9b show the 10-year average water table depth fields, demonstrating reduced water table depths in scenarios where tile drainage is incorporated into the model.Locations in the tile-drained scenarios where water table depths are higher correspond to urban regions where tile drainage has not been implemented.This model instance does not include a representation of impermeable surfaces of urban areas; had impermeable surfaces been included it is likely the water table in the small fraction of the domain with urban areas would have been lower.Figure 9c illustrates the efficacy of the tile drainage model by comparison distributions of spatially averaged water table depths between scenarios with and without tile implementation.These comparisons reveal the model's capacity to mimic the role of tile drains in regulating water table levels.

Effects Tiles and Ditches on Surface Water
Tiles and ditches significantly influence surface hydrology, particularly in reducing surface inundations.Figures 10a-10c show the 10-year mean surface ponded depth fields across three scenarios: combined surface and subsurface drainage, exclusive surface drainage, and no artificial drainage.Implementation of ditches markedly decreased inundation by facilitating the rapid conveyance of surface water to streams, while the addition of tiles augments this reduction through enhanced infiltration.Water Resources Research

10.1029/2023WR035993
To quantify these effects, we use exceedance probability plots of inundation fraction shown in Figure 10d.The probability of less than 0.5% of the total area being inundated remains approximately the same in all three model configurations (approximately 0.93).In other words, during very dry conditions, the presence of tiles and ditches does not significantly impact the inundation behavior of the system.
However, the introduction of artificial drainage (both tile drains and ditches) substantially reduces the probability of larger areas being inundated.For instance, the probability of 5% or more of the watershed area getting inundated decreases from 50.5% in the absence of artificial drainage to 23.9% in the case of surface drainage only and further to 3.6% when both surface and subsurface drainage were implemented.

Effects Tiles and Ditches on the Flow Extrema at the Outlet
We compare how drainage impacts low flows (flows lower than the 10th percentile), and high flows (flows greater than the 90th percentile) in Figure 11.Tiles and ditches increase baseflow while shifts in high flows (flows greater than the 90th percentile) appear to be minor.High flows particularly peak flows in response to storm events need further investigation as other factors may play a role in how tile drains shift the peak flow response (Boland-Brien et al., 2014;Sloan et al., 2016Sloan et al., , 2017)).Low flows (flows lower than the 10th percentile) increased by about 60% in the median by artificial drainage, with both ditches and tile drains contributing (see Figure 11).It is important to note that high-flow response is dependent on other factors as well (Sloan et al., 2016(Sloan et al., , 2017)), including the spacing of the tiles, as investigated below.

Sensitivity to the Tile Spacing
The Hooghoudt drainage model enables the representation of various physical parameters of the tile-drain network in our models.We investigate tile drain effects for a range of tile spacing values.These effects arise  from intricate interactions between the tile drains and the dynamics of surface-subsurface flow.The implementation of the Hooghoudt drainage model in our fully distributed model simplifies and facilitates such investigations.Widening the tile spacing leads to a reduction in the potential tile-drainage flux, which, in turn, causes an increase in surface inundations, elevated surface water content, decreased topsoil water content, and a general increase in infiltration.For extremely large tile spacing values, all variables asymptotically approach the scenarios with no tiles, which aligns with expectations (see Figure 12).
However, the streamflow response presents an intriguing trend.When the tile spacing is increased, the recession curve initially lengthens until a specific threshold spacing value is reached.Beyond this point, the recession curve begins to shorten again.Another interesting observation from Figure 12 is the sensitivity of event response to tile spacing during winter but not in the summer.

Relationship to Previous Modeling Approaches
Our subsurface drainage representation uses nonlinear drain nodes placed at the drain depth, which has similarities to the seepage node approach that has been used before to represent tile drains in fully distributed  Water Resources Research 10.1029/2023WR035993 hydrology models.However, our approach is unique for fully distributed models in that we use the Hooghoudt drainage model to represent flow to the drain nodes.This approach is physically based and takes as input soil permeability and tile properties (spacing, depth, radius, and distance above an impermeable bottom boundary) instead of a calibrated drainage conductance term.Moreover, the Hooghoudt model splits the flow to the drain into two terms, one representing flow to the drain from below and one from above.The latter term is nonlinear in the pressure head at the drain elevation and is neglected in previous models that used seepage nodes with a linear reservoir approximation.Although not addressed here, the Hooghoudt drainage model has been extended to include entry resistance at the drain (Oosterbaan & Ritzema, 1992).That enhancement and additional algorithms to include finite drain capacity and active drainage water management would further improve the physical realism of the model and range of applicability.
Another important feature of this drainage model is that we collect tile efflux from the contributing area of a surface drainage ditch and then redistribute that flux into that ditch.This routing of tile efflux through the surface system has been done in an approximate way with MIKE-SHE (Hansen et al., 2013) but was typically ignored in most previous studies.Explicit routing through the surface flow system becomes important in modeling events in larger watersheds where the propagation time for a flow pulse through the stream network cannot be ignored.Moreover, routing will become important as we move from flow-only simulations to reactive transport simulations, where it is important to represent processing in the stream corridor.

Role of Artificial Drainage
It is apparent from the literature and our results that the prediction of drainage effects on the hydrological response of a watershed is challenging, given the strong dependencies on local conditions (Blann et al., 2009;King et al., 2014;Macrae et al., 2019;Miller & Lyon, 2021;Robinson & Rycroft, 1999;Skaggs et al., 1994;Sloan et al., 2017;Thomas et al., 2016).Using a spatially resolved watershed-scale model, we explore underlying factors that control the shifts in the hydrologic variables due to drainage.The effect of drainage depends on the scale of the drainage area, local conditions, and event size; thus, they are neither uniform in space nor steady in time (Blann et al., 2009).
Tile drains have a notable impact on both subsurface and surface water dynamics.In the subsurface, they effectively control the water table, keeping it below the topsoil level, and promote enhanced infiltration leading to a reduction in surface inundations and overland flow.Furthermore, our findings reveal that the recession curve is highly responsive to changes in tile spacing in the wet season.Widening the tile spacing initially results in a lengthening of the recession curve.However, as the spacing becomes large enough to weaken the effects of the tiles, saturation excess overland flow begins to occur, leading to a shortening of the recession curve.
The unexpected trend of recession curve shifts with tile spacing in Figure 12 introduces an interesting aspect to our overall findings.One possible explanation for this phenomenon is as follows: In the two extreme cases, that is, closely placed tiles (high potential tile drainage flux) and no tiles, precipitated water quickly reaches streams or ditches albeit through different pathways-saturation excess overland flow in the former and tile-drained flow in the latter.As a result, the streamflow response to precipitation events is primarily controlled by the stream/ditch network.However, in the intermediate cases, the presence of tile drains reduces the saturation of the topsoil, which hinders the initiation of saturation-excess flow, yet the infiltrated water takes longer to reach the stream due to the lower tile drainage flux.The patterns observed in infiltration, surface inundations, and surface water content in Figure 12 provide supporting evidence for this hypothesis.
This disparity in the sensitivity of event response in the wet and dry seasons can be attributed to water table behavior in these two periods.In the summer, the water table remains below the location of the tiles, regardless of whether tile drains were present or not.As a result, the tile drains do not significantly influence the water table during this period and the event response shows minimal sensitivity to tile spacing during the summer months.Conversely, in winter, as the water table rises closer to the tile locations, the efficiency of tile drainage and its effect on the event response become more pronounced.
The effects of ditches are more prominent on surface water, allowing excess water on the surface to quickly leave through ditches and significantly reducing the chances of surface inundations.Because of the explicit representation of the drainage ditches in our model, we resolve the spatial patterns and quantify the reducing effect of ditches on surface inundations (Section 3.3, Figure 10d).As seen in Figures 10a-10c, ditches convey the excess surface water effectively through channels, while in the absence of ditches, pools of slow-moving or stagnated water form.This reduced surface storage capacity and availability of faster flow routes yields higher peak flows.Ditches also indirectly help reduce soil saturation by allowing quick egress of surface water that otherwise would have to infiltrate.Additionally, water in the topsoil periphery of the ditches can seep into ditches.This reduction in the ponded water and soil saturation due to combined surface and subsurface drainage not only reduces the risk of crop damage from excess water but also helps with drier conditions for harvest and farm machinery usage and provides farmers with a broader choice of crop varieties (Blann et al., 2009;Spaling & Smit, 1995).

Process-Based Spatially Resolved Models at Watershed Scales
Due to challenges in resolving ditches and tile drains, the applications of spatially resolved ISSHMs have been limited to farm and small catchment (Carlier et al., 2007;De Schepper et al., 2015;Hansen et al., 2013;Rozemeijer et al., 2010;Thomas et al., 2016).However, the impacts of tile drains are known to vary with the scale of the watershed and event size (Schilling et al., 2019;Sloan et al., 2016), and the insights gained from field-scale studies might not be transferable to larger scales.We provide a comprehensive modeling strategy to scale out to a full watershed scale (>1,000 km 2 ) while resolving tile-drains and agricultural ditches effects at scale.
Fully distributed ISSHMs used in most previous farm or small-catchment scale studies are calibrated for specific sites.Due to computational cost, calibration may not be viable for a full watershed/basin-scale ISSHM.More importantly, reduced reliance on calibration improves confidence when simulating the effect of changing environmental conditions like different agricultural practices, hydrologic intensification, or shifting storm patterns under changing climates.Hence, our model emphasizes resolving dominant hydrological processes like artificial drainage, heterogeneous soil properties, and meteorological forcings, to name a few, as accurately as the available data allow to generate reliable results with reduced need for calibration.Without site-specific calibration, our fully distributed ATS ISSHM for the Portage watershed performed remarkably well compared to the USGS gauge data (Figure 6, normalized KGE = 0.81).Moreover, the ATS model outperformed the calibrated semi-distributed SWAT model (normalized KGE = 0.68).In addition, the ATS simulations produce shallow soil moisture patterns that are quite consistent with SMAP-derived patterns.
The process of building confidence in distributed integrated hydrologic models is an iterative one that relies on multiple lines of evidence.Good performance against integrated metrics like streamflow is one line of evidence, which is strengthened in our case by the fact that site-specific calibration was not required.Consistency with the SMAP-HB product, an independent model that assimilates satellite observations, demonstrates model robustness, another important line of evidence.The fact that ISSHMs are built from component models that have each been extensively studied by the scientific community provides additional support.Additional iterations of this confidence-building activity might compare to groundwater levels, evapotranspiration fluxes, and multiple stream gauges in areas where those are available.

Conclusions
Our strategy for representing the effects of managed agriculture drainage in fully distributed integrated surfacesubsurface hydrology models (ISHHMs) combines ditch-aligned computational meshes, hydrologic conditioning to ensure connectivity in the stream/ditch network, a physically based subgrid model for drainage through tile drains, and an algorithm for redistributing drainage effluent into the surface flow system.The strategy is implemented in Watershed Workflow and the parallel Amanzi-ATS code.Using the new capabilities, we demonstrated high-resolution ISSHM simulations of agricultural watersheds with artificial drainage at scales >1,000 km 2 .Below are key takeaways from our study: • Without any site-specific calibration, our model reproduced streamflow in the Portage River Basin as recorded by USGS gage (normalized KGE = 0.81) and outperformed a calibrated SWAT model.These results lend further support to recent suggestions (Bhanja et al., 2023) that high-resolution ISSHMs models backed by high-quality data products may reduce the need for site-specific calibration.• Our novel strategy of ditch-aligned computational meshes that are hydrologically conditioned to ensure connectivity in the stream/ditch network is an effective and tractable way to represent surface drainage ditches explicitly.This preprocessing workflow is implemented in Watershed Workflow, which builds input for Amanzi-ATS, but the strategy applies to any ISSHM that supports arbitrary polyhedral cells.
Water Resources Research 10.1029/2023WR035993 • Hooghoudt's drainage equation applied at the subgrid level with tile-drained water routed back to ditches is an effective way to represent tile drainage at scale in ISSHMs.The approach is effectively a generalization of previous approaches based on distributed seepage nodes or internal boundary conditions but is more physically based with more direct links to the physical characteristics of the drainage system.Moreover, Hooghoudt's equation has been extended to represent entry resistance and finite drain capacity, which would be a useful extension of the capability presented here.• The strong physical underpinnings and reduced need for calibration allow us to study the impacts of artificial drainage on distributed watershed hydrological response in terms of fluxes and states.
With their spatially explicit representations of land surface conditions, hydrologic states, andwater fluxes, ISSHMs have a unique role as tools for understanding the effects of agricultural management practices.An accurate representation of the surface and subsurface flow systems is a necessary first step.The good results obtained here for flow in a highly managed system covering more than 1,000 km 2 suggest that high-resolution ISSHMs are promising tools for modeling agriculture regions at river-basin scales while still retaining sufficient spatial resolution to represent nutrient transport processes.Future work will leverage the surface/subsurface multicomponent reactive transport (Molins et al., 2022) and multiscale river corridor capabilities in Amanzi-ATS (Jan et al., 2021;Painter, 2018) and further evaluate the potential of contemporary ISSHMs as platforms for a new generation of agricultural watershed models.
This study does not account for the potential reverse flow from ditches into tile drainage systems, which could become problematic should water conveyance in streams and ditches be obstructed by debris or due to backwater effects originating downstream.Given that the study area within the Portage River basin is approximately 30 km from the nearest waterbody, that is, Lake Erie, it has been presumed that lake-induced backwater effects do not significantly influence our domain, and that the drainage channels are routinely maintained to ensure the unimpeded flow of water.Moreover, the framework of this study allows for the integration of drainage management practices, which could be considered and adapted in subsequent applications as required.We assume that the tile drainage is located where soybean/corn is planted, the slope is less than 1.4, and the soil hydrologic group is C or D.

Figure 1 .
Figure 1.The study site is in the Portage River basin (black boundary) with a USGS gage station at Woodville, OH (site ID 04195500) at the domain outlet shown by red pin drop and internal HUC 12 boundaries shown in the red.

Figure 2 .
Figure 2. NHD-flowline, Google Earth images (on 18 March 2023), and NED 3m DEM demonstrating challenges to hydrological connectivity for the stream mesh in an agricultural watershed such as Portage.

Figure 3 .
Figure 3. (a) Conceptual figure showing different types of mesh elements along the stream/ditch network; (b) NHDPlusHR Flowlines for streams and ditches (c) 3D domain with the elevation map from NED 3m DEM; (d) simulation-ready mixed-element mesh for Portage watershed generated from semi-automated meshing workflow showing quad elements and junction; (e) 3D mesh view showing how narrow ditches are resolved using narrow quad elements.

Figure 4 .
Figure 4. (a) Conceptual figure showing areal sink in the tile layer of the subsurface mesh represents effects of the tile-drain system implicitly (left) and areal sinks in the subsurface mesh corresponding to catchments mapped to corresponding reaches in the surface mesh using catchment ids (right); (b) catchments for each reach used for source-sink mapping from reach cells to individual tile cells in the corresponding catchment; zoomed in windows of showing mesh-cells in different catchments with colors of their corresponding mapped reaches.

Figure 5 .
Figure 5. Data sets for soil porosity, soil permeability, soil thickness, tile drainage, depth to bedrock, and land cover mapped onto the computational mesh for Portage Watershed.Data sources and resolutions are provided in the Supporting Information S1.

Figure 6 .
Figure 6.Comparison of ATS-and SWAT-simulated streamflow with the USGS gage data from Woodville Station: (a) time series of flow and precipitation; (b) time series of flow on a logarithmic scale; (c) exceedance probabilities of flow (logarithmic); (d) scatter plot with 1:1 line in red.

Figure 7 .
Figure 7. Snapshots of surface ponded depth field demonstrating that ditches resolved by quadrilateral elements drain peripheral area and channel water effectively.

Figure 8 .
Figure 8. Comparative Analysis of the 2019 Annual Mean Soil Moisture Distribution between SMAP-HB and ATS.Developed areas are excluded (indicated by white spaces) and the SMAP-HB is adjusted for the median difference between SMAP-HB and ATS to focus on the spatial patterns.

Figure 9 .
Figure 9. Effects of Tile Drainage on Subsurface Water Dynamics.Panels (a, b) compare the mean water table depth fields for the 2010-2019 period, with (a) showing the scenario with tile drainage and (b) without.Panel (c) presents the probability distribution of the spatial average watertable depth, demonstrating that the water table is consistently maintained below the tile layer depth when tile drainage is implemented.

Figure 10 .
Figure 10.Effects of tiles and ditches on surface water; (a) comparison of 10-year mean surface ponded depth fields under different drainage scenarios; (b) comparison of the exceedance probability distributions under different drainage scenarios.

Figure 11 .
Figure 11.Comparison of low and high flows from simulations for three drainage scenarios and from USGS gage-boxplots of flow smaller than the 10th percentile, higher than the 90th percentile.

Figure 12 .
Figure 12.Sensitivity of streamflow, surface water content, topsoil water content, average infiltration rates and inundation fractions to tile spacing parameter in Hooghoudt drainage model implemented in ATS.

Figure A1 .
Figure A1.Three example reaches demonstrating hydrological conditioning of the quad mesh: first row: NHD flowline and 3m DEM, with red the red arrow showing the flow direction; second row: elevation profile of reach before and after hydrological conditioning with distance measured from the downstream end of the reach.

Figure B2 .
Figure B2.The location of the tile drainage in the SWAT hydrological response unit (HRU) maps: (a) actual HRU map; (b) HRU with tile drainage.81% of the area is assumed to have tile drainage.The selection of HRU with tile drainage is based on soil hydrologic group, slope, and land cover.We assume that the tile drainage is located where soybean/corn is planted, the slope is less than 1.4, and the soil hydrologic group is C or D.