Monitoring and Modeling the Effect of Agricultural Drainage and Recent Channel Incision on Adjacent Groundwater-Dependent Ecosystems

Abstract

Channel incision isolates flood plains, disrupts sediment transport, and degrades riparian ecology. Reactivation and periodicity of incision may affect the water table and hydrological conditions far beyond the stream margin. Long-term incision and its recent acceleration along Iron Springs Creek, North Dakota, USA, has affected adjacent ecosystems. An agricultural surface drain empties directly into the original spring-fed source of the creek, which triggered channel erosion both up- and downstream. Historical maps, recent LiDAR, and field surveying were used to characterize incision since ditch excavation in 1911. Although the soils are sandy, small hydrological gradients impede natural drainage in the surrounding stabilized dunes. Incision resulting from expanded drainage and increased precipitation has been as much as 5 m. Numerical models of lateral groundwater profiles corroborated with field measurements show that the nearby water table responds quickly, becoming deeper and less variable. With 1 m of recent incision, model evapotranspiration rates are decreased 50% to 15% from the channel margin to 1 km, respectively, and the hydropattern disrupted >1 km. Species diversity is reduced and floristic quality is 25% less near the drain. A near-channel solution to erosion—fencing out cattle—failed to mitigate the problem because a broader watershed approach was necessary.

1. Introduction

Stream channel incision occurs naturally and under conditions of human-induced land cover change, including excavated agricultural drains that essentially constitute “controlled” incision. Channel deepening can occur as a result of many factors, including climate change [1,2], lowered stream base level [3], flashier stream flow (often following human development and greater runoff) [4], channelization [5], and changes in sediment load [5,6,7]. Incision triggered by land cover changes can occur rapidly, with unintended consequences, which include the channel being cut off from its flood plain, stream enlargement [8], severe erosion and deposition [9], and the destruction of riparian zones because of reduced flood water and lowering of the water table [10]. This last process can be profoundly deleterious to groundwater-dependent ecosystems [11]. In areas where soils and sediments are characterized by naturally large permeability, the effect of incision on groundwater drawdown can quickly propagate large distances from the stream channel. The deleterious effects may go unnoticed or attributed simply to other unrelated processes because of landscape heterogeneity and a lack of monitoring and quantitative observation.

In southeastern North Dakota, fine sandy Late Pleistocene sediments underlie the watershed of Iron Springs Creek (Figure 1). Abundant springs and seeps discharge groundwater along the margins of the Sheyenne River valley, where the river during the last 10,000 years has eroded through the fine sands and cut down into underlying clay till. As the name suggests, Iron Springs Creek was fed by a discrete spring before the area was settled and agricultural development commenced (Figure 2). The creek’s topographic watershed expanded greatly in 1911 when a jurisdictional ditch, Richland County Drain 10, was excavated and connected to the creek’s spring headwaters. The hydrological imbalance that ensued led to deep erosion, a permanently lower water table, and subsequent loss of Iron Springs. Within the last two decades, increasing precipitation in the region coupled with expansion of the drainage conveyed into Iron Springs Creek has greatly augmented flow in the original stream and uncontrolled headward migration of incision. Erosion along the channel has resulted in the continued and expanded loss of groundwater-dependent ecosystems in the vicinity of the drain and deeply incised stream.

1.1. Objectives

This report describes the case study of Iron Springs Creek and Drain 10 that shows how unintended consequences of drain excavation, stream incision, and headward erosion can strongly alter groundwater-dependent ecosystems that lie laterally outward from even a small and ephemeral waterway and channel. Objectives of this study include: (1) show how groundwater monitoring and modeling can portray hydrological conditions that have been disrupted by drainage and channel incision, (2) describe how these disturbances alter plant communities at distances far beyond those generally considered to be the channel’s riparian zone, and (3) outline methods to forecast the extent to which the seasonally shallow water table and accompanying evapotranspiration may be altered within the watershed if headward erosion remains unchecked.

1.2. Description of the Study Site

The Sheyenne grassland of southeastern North Dakota (Figure 1) covers roughly 55,000 ha and overlies Late Pleistocene delta and lake underflow sediments that formed within glacial Lake Agassiz [12]. Elevations range from about 290 m in the lower Sheyenne River valley to 350 m in the center of the grassland (Figure S1). As much as 30 m of unconfined, fine deltaic sand underlie the grassland [13], with an approximate 20 m thickness in the vicinity of Iron Springs Creek and Drain 10. Most of the groundwater recharge that occurs in the grassland area flows from the center of the area outward, discharging to evapotranspiration and to springs that create small streams which empty into the Sheyenne River on the north and Elk Creek to the south. The Sheyenne grassland, which covers only about 2% of the state of North Dakota, hosts more than 80% of the state’s native plant species [14]. This unusually large diversity reflects the wide array of ecosystems and contrasting niches that occur in the grassland.

Groundwater dynamics in the region is characteristic of a mid-continental, cold temperate climate, where snowfall accumulates during the winter and melts rapidly in the spring, driving a rapid and large rise of the water table. Much of the annual groundwater recharge occurs during the spring, especially when rain can infiltrate before the large rate of summer evapotranspiration develops. The water table will rise rapidly and then shallow groundwater will be released slowly to evapotranspiration during the summer. Although large summer storms may raise the water table, it is usually short lived with much of the excess water lost quickly to evapotranspiration [15].

Of particular importance are the small groundwater-fed streams and wetlands, which lie near the margins of the Sheyenne grassland and host regionally rare and unusual fish [16] and riparian ice-age plant relicts [14]. These groundwater-dependent ecosystems have been damaged by alteration of local cover, caused by excessive livestock grazing and diversion of excess surface water, as exemplified in this report.

2. Materials and Methods

2.1. Drainage History and Characterizing the Incision

LiDAR data collected in 2008 by the International Water Institute [17] and distributed in one-meter-resolution gridded format by the North Dakota State Water Commission [18] provide a benchmark for topographic analysis prior to and after LiDAR acquisition. To characterize changes since the acquisition of LiDAR in 2008, channel profiles were surveyed in the field during October 2017. To estimate incision, the 2008 channel of Iron Springs Creek and Drain 10, which flows into it, was delineated using LiDAR data (Figures S1 and S2). The gridded digital elevation model DEM panels were mosaicked, digitally filled to eliminate mapping artifact sinks along the channel [19], and processed to find the flow direction of cells and flow accumulation. To identify incipient channels within the DEM, a threshold watershed consisting of 20,000 cells (2 ha) was used. Raster cells representing Iron Springs Creek and the Drain 10 were then converted to a simpler line vector which was superimposed on scanned maps for further analysis.

To obtain an approximate historical record of surface drain and stream incision, historical topographic maps available from the U.S. Geological Survey were scanned at high resolution (Figure 2) and georeferenced using Public Land Survey lines mapped and digitized in a UTM zone 14N NAD 1983 projection. Survey line intersections on the most recent digital topographic map and the scanned historical map image were used as match points. Historical topographic maps at a scale of 1:125,000 and a topographic contour interval of 6.1 m (20 feet) for Wyndmere and Casselton are available for the entirety of Iron Springs Creek and Drain 10. Topographic surveys for the Wyndmere (1904) and Casselton (1894–95) quadrangles were conducted prior to the excavation of Drain 10 in 1911. More recent 7.5 minute, 1:24,000 scale topographic maps (McLeod and Coburn quadrangles), based on aerial imagery taken in 1952 and 1958 together with 1960–61 plane-table surveys, were used for acquisition of later data. To estimate the channel elevation on the historical topographic maps, the point closest on the current channel to the upstream inflection of the contour on map was assumed to be the elevation of the channel at the time the topography was surveyed.

To measure recent incision, profiles across Iron Springs Creek and Drain 10 were surveyed in 2017 using a Topcon laser level and receiver positioned at the end of a 7.62 m (25 foot) stadia rod. The level was set on an elevated point above the profile location and the relative elevations of least 12 points, each separated laterally by approximately 5 m, were used to establish elevations across the channel. GPS coordinates at each point were recorded using a Trimble Juno model 3B with Trimble’s GPScorrect extension installed. In addition, 2008 LiDAR grid elevations were interpolated for the corresponding x-y coordinates along each profile. This provided a correlated set of elevations for 2008 (LiDAR) and 2017 (field survey).

To estimate the amount of incision that has occurred since the acquisition of LiDAR data, little or no change in elevation was assumed to have occurred beyond the confines of the drain and stream during the 10 years following the survey. Correlating LiDAR points away from the drain is how the two sets of data were correlated and referenced to the same vertical datum (mean sea level). Nevertheless, interpolation of the scattered raw LiDAR points into a digital elevation model grid and 10 years of ground disturbance by grazing cattle, some small differences in elevation are expected. To minimize this uncertainty, survey profiles were extended to a distance at least the width of channel beyond the margins of the channel to provide elevation control for areas assumed to be unaffected. Based on the portions of the surveyed profiles beyond the outer margins of the channel, the maximum error in vertical incision is likely to be ±20 cm, while the error in lateral incision is less certain. The error estimate was determined by comparing the LiDAR and surveyed profiles varied beyond the banks of the drain, where they have most likely remained unchanged.

2.2. Selecting Blocks and Transects

The heterogeneity of the dune topography and its patchy relief of 2 to 7 m (Figure S3) indicates that the depth from the ground surface to the water table varies locally both temporally and spatially. Swales are natural depressions that lie between the dunes and host jurisdictional wetlands in most areas unaffected by excavated drains or tile. Visual appearance, such as droughty swales lacking a typical diversity of wet meadow plant species, suggested strongly that groundwater-dependent ecosystems in the Sheyenne grassland had been perturbed by nearby agricultural drains.

The approach to monitoring focused on collecting data to characterize the small-scale changes in hydrology and vegetation related to individual dune–swale systems, including areas affected by drain drawdown, in the context of the broader hydrology of the landscape. Selecting a single long and straight profile would confound the local detail required to characterize groundwater and vegetation in the dune–swale landscape. Therefore, a way was needed to monitor and define the variability of the water table level in the heterogeneous dunes. To do this, a GIS approach was used that depends only on topographic characteristics to identify dune crests and swale bottoms within the local landscape. Dune topography and hydrology were characterized within two quadrants, both up-gradient, with one in the disturbed zone close to the drain and the other beyond obvious effects from the drain (gray rectangles, Figure 1).

Initially, the topographic position index (TPI) [20,21] within the study’s two quadrants (Figure 1) was calculated using a digital elevation model (DEM) with one-meter raster cell dimensions, obtained from the 2008 aerial LiDAR data. The unitless TPI measures how the mean elevation varies within a specified neighborhood adjacent to each DEM raster cell. The smallest scores indicate a low area on the landscape with most immediate neighboring cells lying at a higher elevation. Similarly, the largest TPI scores correspond to cells with the majority of neighbors lying at a lower elevation. Following trials using different sized radial neighborhoods, a 50 m radius was applied for this application. Distinct dune tops and swale bottoms observed in the field and on topographic maps were readily identified and quantified using these parameters, suggesting that this selection of TPI scale and geometry were appropriate for the dune–swale landscape within the quadrants.

The resulting TPI values were divided into ten quantile classes, which were used to discriminate between dune crests (cells scoring within highest quantile, 90% of all the scores were smaller) and swale bottoms (cells scoring within the lowest two quantiles, 80% of all the scores were larger). Raster cells representing the dune crests tended to be more isolated and distinct than the interconnected and irregular swales. Therefore, dune crests were more easily captured statistically when compared to the swales. Only the one upper quantile was needed to categorize dune crests whereas the two lowest quantiles were necessary to delineate swales. Centroids for groups of adjacent swale cells and groups of adjacent dune crest cells were used to define a specific point that represented the discrete swale and dune crest locations. Potential monitoring transects throughout the quadrants were then established by connecting a line from each dune-crest centroid to the nearest swale centroid. The actual transects used in the analysis were selected randomly from the full set of several hundred dune crests identified within the quadrants (Figure 3).

2.3. Monitoring the Water Table

For hydrological monitoring, three transects were selected in the more distant reference quadrant, which was less affected by drainage, and four transects selected from within the disturbed quadrant adjacent to Drain 10. Monitoring points and corresponding vegetation characterization plots along the transects (Figure 3) were established in the field by locating the swale bottom, dune crest, the midpoint, and two other points along that continuum: low transition and upper transition. Identification of the transition points was based on a visual assessment of the topographic position and the distribution of upland versus wetland vegetation communities. The profiles ranged from 10 to 50 m in length, depending on the distance between the swale and dune summit.

Water-level piezometers and monitoring wells were installed by hand in late 2012, with the former constructed using a 2.6 cm diameter pipe attached to a 0.3 m screen at the base. Monitoring wells were constructed with 10.2-cm, schedule-20 PVC pipe and 1.5 m screens. The top of the screen for all was set just below the dry season level of the water table. Riser tops were surveyed to the nearest 0.003 m using a laser-level-line transect that extended across the two quadrants. Discrete water-level measurements in the wells and piezometers were made from May through October 2013. Continuous measurements were collected in selected wells during July–October 2013 using unvented Solinst level loggers and corrected for changes in barometric pressure.

2.4. Modeling Water Levels and Water Budget

The objective for creating the numerical model was to: (1) characterize the short- and long-term drawdown related to recent incision of the drain and (2) estimate the seasonal water-level variability that occurs within the area affected by drawdown near the drain and at more distal points. Because there is little or no surface water flow across the sandy soils of the Sheyenne region, the variability within the natural water budget arises from the interaction of infiltration, evapotranspiration, and recharge. Evapotranspiration is a function of water-table depth; its interaction with infiltration and recharge creates non-linearity in the water budget. Because of this coupled non-linearity, simple analytical solutions cannot adequately describe the dynamics of the system and a more robust modeling approach was necessary [22]. MODFLOW-2005 (version 1.11.00, U.S. Geological Survey, Reston, Virginia, USA), a full-featured finite-difference groundwater modeling code [23], was used to explore the effects of incision on original water budget conditions.

As required in MODFLOW’s evapotranspiration package with segments (ETS), a maximum rate of evapotranspiration occurs when the water table lies at the evapotranspiration surface (approximately the ground surface in this application) and diminishes piecewise linearly to zero as the depth of the water table falls to the extinction depth [24]. Parameters for the variation of evapotranspiration with depth below the evapotranspiration surface were obtained from the data compiled by Shah et al. [25]. Because the pattern of evapotranspiration is much more complex under natural conditions (e.g., [26]), this provides only an approximation and does not consider moisture dynamics above the water table.

Other characteristics of the MODFLOW models used for the site-specific application, such as the boundary conditions, stress periods and time steps, spatial discretization, and matrix solving methods are presented in the Results section. The basic input files used for the simulation are provided in Table S2.

2.5. Analysis of Plant Assemblages

Dune-swale transects were established as described previously (Section 2.2) up to a distance of 2000 m from Drain 10 (Figure 3). All plant species within small quadrats along the seven selected transects were cataloged and categorized (Table S3). It is likely that the water table varies spatially relative to the drain and temporally both intra- and inter-annually; thus, each transect included plots from four local topographic zones (dune crest, swale bottom, and two intermediate positions) along each transect. Seven transects of four or five quadrats each were sampled. The following protocol was used for characterizing the vegetation in 2012 and 2013:

• GPS was used to navigate to a dune crest at which point a transect was extended from the summit to the closest swale point previously identified in GIS.
• Five vegetation plots were examined along each transect:
  • crest point of the dune
  • swale point
  • midpoint between these points
  • an additional position referred to as the dune–swale transition was chosen to capture the vegetation continuum between the swale and dune slope, and/or
  • a low transition was selected as a point representing vegetation characterizing the margin of the lowest and wettest point on the dune–swale continuum.
• At each of the four or five sampling points on each transect, the following procedure was carried out:
  • all species present within a 1 × 1 m (one square meter) quadrat were recorded, including both vascular plants and bryophytes. This step was completed when the plants were actively growing and robust, after emergence and before senescence.
  • the local zone category was recorded (from Step 2 above).
  • the GPS coordinates of the point were checked again and re-recorded to the nearest meter.
• Data for the 2012 and 2013 field seasons were combined and analyzed.

Plant species for the quadrats in each of the seven transects were combined and used to describe qualitative statistics for the community, based on regional species information [27,28]. These included the number of species, ratio of native to non-native species, coefficient of conservatism (C), adjusted floristic quality index (aFQI), and a wetness index equal to the species count ratio of the obligate plus facultative-wet species to upland plus facultative upland species. No attempt was made to quantify the abundance of any one or more of the species along transects.

The coefficients of conservatism (C) assigned to individual species are based on tolerance to site degradation and need for natural remnant habitats. The C values range from 0 to 10, with the most conservative species (C values > 7) usually associated with each other in locations where the community and species evolved together under long-term, stable conditions. In contrast, the least conservative species (C values < 3) are adapted to strong anthropogenic or natural degradation that would eliminate both high- and mid-conservative species. Non-native species are given a coefficient of 0 [29].

The C values for individual species that occur at each transect were used to calculate the aFQI metric:

$$aFQI = 100 \left(\frac{{\overline{C}}_N}{10}\right) \left(\frac{\sqrt{n_N}}{\sqrt{n_T}}\right)$$

where ${\overline{C}}_N$ is the mean C for native species, $n_N$ is the number of native species, and $n_T$ is the total number of species. Because the coefficient of conservatism for non-native species is 0, the adjusted floristic quality index includes the contribution of non-native species and provides a better balance when comparing sites that have been disturbed [30].

3. Results

3.1. Drainage History and Characterizing the Incision

The location of the LiDAR-derived channel corresponded to within 175 m on the georeferenced 1894–95 Casselton map, except for the lowest 500 m reach which showed Iron Springs Creek flowing along a more western course immediately above its confluence with the Sheyenne River (Figure 2). It is unknown if the stream at that time flowed in this course and later moved into its present more northerly course that lies to the east. In addition, the 1894–95 survey shows no measurable difference between the current and past elevation of the Sheyenne River, although minor incision is likely to have occurred during the last century because of dam construction upstream and local river channelization. The linear longitudinal profile of Iron Springs Creek at that time suggests that the channel may have accommodated a steep drop near its outlet (Figure 4).

In response to needs for agricultural drainage, Drain 10 was let for contracts in 1907 [31] and excavated in 1911 from the original location of Iron Springs southward about 10 km. Although no record of the elevation of the natural spring exists, both the early 1894–1904 surveys and 1960 topographic maps suggest that its position was at 319.5 m above mean sea level and about 4950 m upstream along the natural channel from the Sheyenne River (Figure 4).

The 2008 LiDAR survey provides the only continuous longitudinal profile for Iron Springs Creek - Drain 10 (Figure 4). Its present strongly convex profile likely does not reflect natural, long-term steady-state conditions [32]. The 2008 profile also shows clearly a channel knick-point about 6 km upstream from the Sheyenne River (Road 1212 crossing). A rock and sheet pile weir, constructed in 2001, was present at this location in 2008 (Figure 5). A combination of snowmelt, rainfall, expanded drainage, and baseflow associated with unusually high groundwater levels throughout the watershed in early 2011 destroyed this structure. Rapid incision upstream to the concrete and rock crossing within the channel bed at Road 1211 ensued (Figure 4). Profiles surveyed during 2017 (Figure 6) show up to 1.6 m of incision directly below the Road 1211 crossing.

On the 1960 topographic data, the large gradient near the confluence continues to appear, although less than one-half the drop compared to 1897. Furthermore, the LiDAR-delineated channel corresponds much more closely to the 1960 map, having no more than a 32 and 2.5 m lateral and vertical separation, respectively. The reason for the lateral difference revealed on the earlier maps—whether a consequence of the natural lateral movement of the channel or error in the map surveys—cannot be easily determined. Regardless, the longitudinal profile is more revealing for characterizing the changes in Iron Springs Creek that have occurred since its connection to Drain 10 in 1911.

The only instance where the channel appears to lie at a higher elevation at an earlier date occurs near the Road 1212 crossing in 1960 when compared to 2008. In this case, it is likely that sediments were deposited upstream from the weir and sheet pile that were installed in 2001 (Figure 5). The most stable reach of the longitudinal profile lies near 3.5 km above the confluence with the Sheyenne River, where incision has been minimal since at least 1960 (Figure 4).

In summary, the pattern of the longitudinal profile, which appears to have been nearly linear prior to excavation of Drain 10 and its conveyance into Iron Springs Creek, has now become more convex, with the clear development of a knick-point that has migrated 4.5 km, from roughly a point 3.5 km above the Sheyenne River in 1960 to about 8 km in 2017. Below the hardened crossing of Road 1212, as much as 4 m of channel incision have occurred during the 60-year period and approximately 5 m since 1911 (Figure 4), with accelerated incision between the original spring and Road 1211 since 2008 (Figure 6).

3.2. Monitoring and Modeling the Effect on the Water Table

In June 2013, the water table level varied by 3 m between its lowest elevation near Drain 10 and highest elevation in the reference block that lies about 2 km southwest from the channel (Figure 3). This difference in level was reduced to about 1 m by September 2013, but maintained the same pattern. A change in the June hydraulic gradient is obvious at a distance of about 800 m from the drain, beyond which the water table gradient becomes more gentle (5 × 10⁻⁴) compared to that near the drain (3 × 10⁻³).

Groundwater flow was modeled to characterize the effect that incision has had on the water table and associated ecosystems between Road 1212 and Road 1211, west of Drain 10 (Figure 1), where approximately 1 m of recent incision has occurred. Drain 10 accumulates seepage as it intersects groundwater, which is conveyed to Iron Springs Creek, and thus draws down the water table along its margins. The 20 m thickness of generally homogeneous fine sandy soils facilitates a prompt temporal and spatial adjustment of the water table to the deepening channel. The natural water table in the Sheyenne grassland tends to be shallow [13,33,34] and often lies above the ground in the swales seasonally, which was observed in Spring 2013. Water infiltrates and percolates quickly within the conductive sand sediments that underlie the dune–swale landscape, as revealed by near-continuous logging of water levels under dunes adjacent to swales (Figure 7). Note that it is only brief occasions following rainfall when the water table is separated vertically by more than a few cm between adjacent swales and dunes.

Input conditions and properties for the numerical model of groundwater flow above Drain 10 are summarized in

Table 1. Temporal and spatial discretization, aquifer properties, and solution conditions used in the MODFLOW model.

Temporal Discretization
total time	20 years (the time required to reach essentially steady-state conditions)
stress periods	60—3 per year (repeats with a 30, 150, and 185 day period each year)
time steps	3, 8, and 9 steps for the 3 stress periods per year, time-step multiplier = 1.2
Spatial Discretization (Figure 8)
layer	fixed—one unconfined layer with the top set at 30 m
columns	fixed—26 (increasing width from the drain boundary)
column width	2, 3, 4.5, 7, 10.5, 13, 18, 27, 35, 50, 70, 100, 150, 200, 250, 360, 500, 700 ….
	1000, 1250, 1500, 1750, and 4 × 2000 m
row	fixed—1, length set at 1 m
drain	column 1, row 1 (2 m wide column), elevation = 20 (initial) 19 (post-disturbance), conductance = 2, 7, and 24.5 m3/day
Aquifer Properties
hydraulic conductivity	variable—2, 7, and 24.5 m/day
specific yield	variable—0.06, 0.12, and 0.24
recharge rate	fixed—0.008 m/day (applied during 30 day stress period)
evapotranspiration
rate	fixed—0.0013 m/day (applied during the 30 and 150 day stress periods)
surface	fixed—20 m (at the ground surface)
extinction depth	fixed—18 m (2 m below the surface)
exp decay parameters	yo = −0.015, d″ = 0.32 m, b = 2.6 m−1
yearly water deficit	fixed—1.0 (potential evapotranspiration to recharge ratio)
Solution
initial head	20 m (original condition: non-incised drain), 19 m (incised drain)
matrix solver	conjugate-gradient method (using default MODFLOW2005 variables)

and Figure 8. The model was constructed by incorporating the following assumptions:

• Recharge, flow, and discharge at points along the groundwater flow line from areas of recharge to discharge can be modeled satisfactorily by using a quasi-two-dimensional profile (Figure 8).
• A no-flow groundwater divide (zero-flux Neumann condition) occurs near the center of the regional groundwater mound that underlies the Sheyenne grassland. A narrow (2 m) drain comprises the boundary on the opposite end of the profile. A simple lowering of the drain elevation replicates the effect of drain excavation and incision.
• Surface water flow does not occur except to accommodate discharge at the drain. The drain does not recharge the groundwater model system, although monitoring results show there may be some infiltration and bank storage during peak flow events.
• The aeolian and deltaic sediments are assumed homogeneous and isotropic, and sufficiently coarse to result in a thin tension-saturated zone and a steep decrease of soil moisture above the water table. Therefore, modeling of variably saturated flow would not significantly affect the results.
• Local differences in the water table level between dunes and swales are ephemeral, brief, and do not change the broader patterns of groundwater flow (Figure 7).
• Recharge is modeled as a single annual pulse, 30 days long. This assumption is based on the observation that frost-related exfiltration occurs during the winter months, which has a similar effect on the water table as does evapotranspiration except that water is stored in the shallow vadose zone. Melting and infiltration of this stored water coupled with snowmelt and early spring rainfall creates a brief but strong recharge event [15].
• Evapotranspiration is distributed equally throughout a 150-day period during the year, which replicates the typical natural freeze-free growing season from April through October.
• The extinction depth, defined as the water table depth below which ET ceases to occur, was established at 2 m, with decay defined by piecewise linear increments [25] and the evapotranspiration surface set at ground level.

3.3. Spatial and Temporal Discretization

A single unconfined layer with variably spaced columns was used for the model (Figure 8), with a MODFLOW drain cell on the east (2 m wide) Finite-difference columns are progressively narrower toward the surface drain, following the general convention that cells be no more than 1.5 times the width of adjacent columns (Table 1). The width of the profile extended 16,000 m, which is roughly the distance from the incised drain west to the groundwater divide within the center of the Sheyenne grassland (Figure 1). Variable aquifer properties, hydraulic conductivity and specific yield, along with recharge that depended on the distribution of evapotranspiration, were adjusted to calibrate the model and test sensitivity relative to the observed field measurements. Calibration and sensitivity within a range of model properties was based on the water levels observed in the two monitoring quadrants adjacent to Drain 10, while maintaining the parsimonious assumptions noted previously.

The simulation was performed using two modes, one configured for the elevation of the original drain (20 m above the unconfined base) and then simulated again at 19 m, which represented the post-incision condition. Models were generated for 20 years, with each year consisting of three stress periods, with 3, 8, and 9 times steps, respectively (Table 1). Output from both the initial and post-incision modes showed that a new near-steady-state condition developed after about 5 years. Initial hydraulic heads from the original drain configuration were used as input for the post-incision simulation.

3.4. Flow System Properties

Calibration of the model was done by matching the general pattern of water levels from the center of the aquifer (within the center of the Sheyenne grassland) to the spring and seep discharge areas along the margins of the aquifer and the Sheyenne River. Data from the North Dakota State Water Commission [18] for 40 wells in the aquifer were selected from a set of 83 wells that include water-level information (Table S4). The 43 wells eliminated from the set had either incomplete data, data that ended before 2010, spatial duplicity, a screen completion at depths below the unconfined aquifer, or most commonly, were located in areas that were beyond the proximity to Drain 10.

To characterize the variability, the groundwater levels for the wells were analyzed and compared to the Palmer hydrological drought index (PHDI) (Figure 9). The secondary scale on the right side shows the variability, which is based on the water levels in 40 shallow groundwater wells, screened in the unconfined Sheyenne aquifer and recorded monthly (April–November) by the North Dakota State Water Commission (Table S4). Because of multiple wells and water levels lying at various elevations, the records for each observation well were divided by the maximum water-level range tabulated for the well during the period 1994–2018. These normalized results were combined for each measurement date, which provided the data to find the trend of the median and first quartile range. Results show that the water levels track closely with the PHDI (Figure 9). Plotted relative to measurements made on 9 May 2016, when the PHDI close to zero, shallow groundwater levels remained within the interquartile range, even for this period of strong PHDI variability, except for a few months during the exceedingly wet period in late 2010.

The highest and lowest water levels for the period 2010–18 were tabulated for the wells and showed median difference of 1.4 m, with an interquartile range of 0.3 m. The highest and lowest water levels were almost exclusively limited to May 2010 and August 2017, respectively. These dates correlate well with the pattern revealed by the PHDI for the southeastern North Dakota climate division at those times [35] (Figure 9). Although the PHDI was lower during 2012 and 2013, groundwater levels generally remained higher than those recorded during late summer 2017.

Model outputs indicate significant changes in both water levels and the spatial distribution of evapotranspiration rates (Figure 10 and Figure 11). Although the periodicity and amplitude of the wet (spring) and dry (fall) water-table levels remain very similar following incision, the effect of the relatively permeable sediments leads to new lower level for both periods at distances greater than 3 km from the drain (Figure 10). Most of the transient lowering of the seasonal water table occurs within 5 years after channel incision (about 90%), but reached a new quasi-steady-state equilibrium by 20 years. Water-level measurements during 2013 show a close correlation to the model levels, although some bank storage is evident at distances up to at least 600 m (Figure 10).

The model suggests that a significant change in the evapotranspiration rates occurs in the profile following incision (Figure 11). The deeper water table, during both the wet and dry periods, results in rates that are lowered by 50% during the dry part of the growing season (Figure 11). Disruption to evapotranspiration rates occurs to distances as great as 2 km, with at least a 15% change (0.0002/0.0013) up to 1000 m (Figure 11).

3.5. Plant Communities

Plant assemblages from seven dune crest to swale bottom profiles were characterized in 2012 and again in 2013. As described in the methods, these profiles had either four or five quadrats where plant species were identified and cataloged. Results show that the number of species recorded ranged from 46 (profile 19:1900 m from the drain) to 25 (profile 38:240 m from the drain), with a generally increasing ratio of native to non-native at greater distances from the drain (

Table 2. Summary of floristic data collected along the transect profiles (Figure 3). The X and Y coordinates denote the UTM zone 14N NAD83 coordinates, FQI is the floristic quality index, %C is the percentage of species along the profile that have a conservatism value ≥ 7, and % wtld (wetland) species column gives the obligate plus facultative-wet species as a percent of the total number of species.

Profile ID	X (m)	Y (m)	Distance to Drain 10 (m)	Adjusted FQI	%C Value ≥7	Total Species	Native/Non-Native	% Wtld Species
18	631447	5147380	2033	44.7	22	36	5.00	33
19	631565	5147437	1900	44.2	20	49	5.57	31
23	631769	5147553	1685	44.3	21	44	4.50	28
28	632135	5147570	1322	45.1	21	34	7.50	27
37	633094	5147729	356	39.0	11	28	4.60	41
38	633266	5147654	240	39.0	20	25	2.57	48
40	633460	5147572	67	34.2	7	30	1.73	14

). The adjusted floristic quality index (adj FQI) ranged from 45.1 (profile 28:1322 m from the drain) to 34.2 (profile 40:67 m from the drain), matching a pattern similar to the percent coefficients of conservatism greater than 7 (Table 2). The presence of robust invasive species, Euphorbia esula (leafy spurge) and Poa pratensis (perennial bluegrass), along with two regionally common sedge species (Carex inops subspecies heliophilia and C. pellita), were noted at all the profiles. Interestingly, the profile site closest to Drain 10 (profile 40:67 m from the channel), hosted the largest number of species unique to that single profile (Table 2), although all species but one (Packera paupercula) were common non-native species of low conservation value, suggesting that the channel may serve as a conduit and pathway for the introduction of exotic and potentially invasive species.

As suggested by the model results, the permanent lower water table near the incised ditch has put obligate and facultative-wet species out of reach of the water table. Surprisingly, no clear pattern was observed for the dominance of wetland and facultative-wet species within the distribution of the profiles. More time may be needed for the disappearance of these species. Alternatively, the wetland species may be able to adapt to non-phreatic, drier conditions, especially for the more common non-native species that have apparently been established within 1000 m of the ditch (Table 2).

4. Discussion

The results show that uncontrolled incision on Drain 10 has led to not only the obvious erosion and destruction of the natural seeps at the head of Iron Springs Creek, but to a lower water table and loss of wetlands. In the fine sandy sediments of the Sheyenne grassland, the drainage of shallow groundwater extends several km away from the incised drain, but the obvious effect of the drainage is perhaps obscured by the naturally large amplitude and interannual periodicity of water levels in the unconfined Sheyenne aquifer (Figure 12). As shown in other research in different environments, this fluctuation can result in disrupted and altered riparian ecosystems [36,37]. The results indicate that obvious hydrological changes to the integrity of groundwater-dependent ecosystems, which are unrelated to any natural riparian system, extend at least 2 km from the drain.

As noted in the results, even if the change in the periodicity of the water table remains largely unaltered, the pattern of evapotranspiration can change significantly depending on the seasonal depth of the water table (Figure 11). For example, in the case of Drain 10, evapotranspiration will be forced to adjust with a deepening water table. Even with as little as 1 m of incision, the affects will propagate to a distance more than 2 km from the channel. Modeling of species dynamics by Loheide and Booth [11] revealed that channel incision resulted in the loss of obligate species and replacement by facultative and facultative upland species. This change in the predominance of wetland species was not observed (Table 2) along Drain 10. Nonetheless, rapid adaption of the plant assemblages to varying moisture regime has been shown to occur in the Sheyenne grassland [38], and these data suggest that the disturbance provides a niche for the loss of diversity and the expansion of weedy species (Table 2), especially close to the channel.

The need for optimizing soil moisture conditions for agriculture [39] has led to the development and application of a variety of analytical tools for predicting the effect of drainage systems (e.g., [40,41,42]). Although these are generally easy to implement and used regularly by agencies that oversee and advise agricultural activities [43], none provide for more complicated conditions related to the effects of soil heterogeneity and the non-linear interaction of recharge, evapotranspiration, and a dynamic water table. Nonetheless, the results for Ditch 10 and Iron Springs Creek indicate the significance that these more complex conditions engage on natural plant communities, even in the unusually homogeneous Sheyenne grassland landscape (Figure 12). Therefore, a numerical approach to predicting the effects of drainage may be necessary for routine analysis of drainage projects, especially if groundwater-dependent ecosystems may be affected. For fine-grained soils, consideration of vadose processes effects may also be necessary [11,44].

The results of monitoring and modeling emphasize the importance of recognizing that local problems, such as channel incision and unintended groundwater drainage, often relate to a watershed problem that requires an integrated approach to an underlying large-scale problem. For example, in 2005 a group of representatives from federal and state agencies, local water and soil conservation board, a non-profit conservation organization, and land owners formed the “Sheyenne Valley Coordinated Resource Management Work Group”, which developed and partly accomplished a plan for mitigating the damage along the Iron Springs Creek and Drain 10 [45]. Although the proposal and plan dutifully recognized that changes in the upstream watershed were aggravating incision and erosion downstream, the funding and effort of the work went toward constructing a rock and sheet pile weir on Road 1212 (Figure 5) to reduce erosion rates and installing extensive riparian fencing to prevent cattle from damaging the stream and drain channel banks. The weir washed out in 2011 (Figure 5). Although controlling the access to the stream by cattle with fencing has improved forage conditions and mitigated direct bank erosion, it has had little to no effect on incision.

The problem recognized and not addressed was the large expansion of the Drain 10 and Iron Springs Creek watershed caused by new ditch excavation and drainage of excess water in the southwestern part of the watershed in 2010–11 (Figure 1). In addition, increased erosion and incision are likely as agricultural tile drainage expands in the Drain 10 watershed, which is currently developing rapidly in the Sheyenne River and upper Red River basins (e.g., [46,47]). Tile drainage contributes to initially sediment-free water and flashier downstream peak flow, both of which lead to accelerated erosion [6,48,49]. These are large, watershed-wide changes that require a broad management solution [50] focused on preventing the conveyance of excess water into Iron Springs Creek.

Hopkins and Running [38] suggest that the Sheyenne grassland plant communities are very well adapted to drought and deluge, which is further supported by the characterization of incision along Drain 10 and Iron Springs Creek. The plant community will change in response to changing hydrology, albeit the changes may be subtle. The observation that plant communities can adapt readily to the changes corroborates with other studies showing species communities correlate with changes in the depth to groundwater, but that vegetation cover is only partially coupled to groundwater depth [51].

5. Conclusions

Increased flow and flashiness of discharge have led to the incision of Iron Springs Creek in southeastern North Dakota, which demonstrates unintended consequences that occur following disruption of a natural hydrological system. The contributing factor was the excavation of Drain 10 in 1911 and conveyance of excess water into the original groundwater-fed stream. The lack of channel data preclude any detailed analysis of the incision process, although U.S. Geological Survey topographic maps from 1897–1907 and 1960, together with more recent LiDAR data gathered in 2008 and a field survey in 2017, suggest that up to 5 m of incision has occurred since excavation. Excessive runoff in 2011 washed out a rock weir that was temporarily stopping headward advance of incision. Recent expansion of drainage in the Drain 10 - Iron Springs Creek watershed has produced larger, flashier flows, which creates more potential for incision and further expansion of the watershed.

Landscape features along the margins of the drain channel locally vary greatly, with wetland swales as close as a few meters to xeric dune crests. Rather than attempt to monitor along a continuous profile, the analytical GIS approach used in this study was to identify high and low points of the dune tops and swale bottoms, respectively. The nearest adjacent high and low points within a landscape block were then used to define transects to characterize plant assemblages and measure water table elevations. This method provided a means to gather a robust set of field data that integrated the landscape heterogeneity while minimizing the effort expended to collect the data.

The spatial and temporal dynamics of water table were determined both through monitoring selected profile sites and by numerical modeling using boundary conditions derived from regional hydrogeology. Calibration of the model was based on water levels and aquifer properties within the Sheyenne grassland region recorded by the North Dakota State Water Commission. Hydrologic variability and contrast is mostly vertical rather than lateral, with large synchronic rise and fall of the water table throughout the area, although agricultural drains and incision (exemplified along Drain 10) dampen the effect significantly within at least 2 km from the channel.

Results of this study couple the loss of plant assemblage diversity to groundwater drainage caused by the Drain 10 and its recent incision. Although diversity is less near the drain, no clear differences were observed in the relative number of obligate and facultative-wet species within the disturbed and reference quadrants. Another feature influencing ecology in a subtle way is the large change in evapotranspiration rate. Although the relative seasonal rise and fall of the water table is generally maintained in areas beyond the channel, the water table is deeper below the local surface at all times. Therefore, drain excavation and channel incision causes a more stable water table depth that lies consistently deeper, putting roots out of phreatic contact and reducing evapotranspiration. For the soil and hydrological conditions near Drain 10, this process is likely to be most significant (greater than a 15% change) up to 1 km up-gradient from the channel, which diminishes to a negligible amount at a distance of 2 km.

Lastly, resource managers must be able to recognize the source and scale of the problems associated with drainage and incision. Although the local land management community met and acted on the problem of erosion along Iron Springs Creek, only mediocre mitigation of the problem resulted from their efforts. Unfortunately, the committee identified a problem that was local (e.g., cattle access to the creek), which was mistakenly identified as the primary source of the incision. The problem of altered watershed hydrology was recognized but overlooked when significant resources went toward fencing out the cattle.