An integrated assessment of impacts to ecosystem services associated with prairie pothole wetland drainage quantifying wide-ranging losses

Our work focuses on the Prairies Ecozone of southwestern Canada, lying at the northern extent of the Great Plains of North America. The Prairies Ecozone spans southern portions of the provinces of Alberta, Saskatchewan, and Manitoba, and most of this area lies in the Prairie Pothole Region of North America, which is widely known for the myriad depressional wetlands left by recent (10 000 ybp) deglaciation. Topography across much of the region ranges from flat to hummocky, with the depressions playing an important hydrological and biogeochemical role; large areas of the region do not connect to higher order streams in most years (non-effective areas), and drainage networks are poorly established (Stichling and Blackwell 1957; Shook et al. 2015). This regular hydrological disconnection is important for trapping water to permit groundwater recharge, and solutes to limit nutrient loading to downstream water bodies. The regional climate is described as semi-arid to sub-humid, with conditions varying strongly by year as the region experiences dry and wet phases. Mean annual temperatures vary from 1 to 6 °C (Mekis and Vincent 2011) and annual precipitation ranges from 350–650 mm (Gray and Landine 1988), but is typically less than potential evapotranspiration (600–800 mm year^–1; Morton 1983).

Under these regional climate conditions, frozen soils limit infiltration during snowmelt, which makes it the key hydrologic period for producing upland runoff. Snow is redistributed across the landscape during the winter, where it accumulates in depressions (Pomeroy et al. 2007). Because of the excess snow and runoff available during spring melt, depressional wetlands fill, and surface runoff networks can be activated when wetland ponds reach their spill point (Leibowitz and Vining 2003; Shaw et al. 2012). Under these conditions, pond volumes peak in the spring and decline over the course of the warm season. Many of these depressional wetlands have water ponded above the surface for periods ranging from weeks to years (Stewart and Kantrud 1971). These ponds provide critical habitat for a wide variety of biota. This includes many bird species (Mantyka‐Pringle et al. 2019; Elliott et al. 2020), and the Prairie Pothole Region attracts a disproportionate number of North America's breeding ducks, producing up to 50% from only 10% of the land mass (Greenwood et al. 1995; Doherty et al. 2018).

Delineation of small (∼100 km²) watersheds across the Prairie Ecozone with the HydroSHEDS global dataset has identified more than 4000 distinct watersheds (Lehner et al. 2008). These have been classified into seven different watershed types using a hierarchical clustering analysis according to their biophysical characteristics, including wetland characteristics and land cover (Wolfe et al. 2019a). Subsequent investigation of the hydrology of each watershed class (He et al. 2023) revealed that each displays different hydrological behaviours and sensitivity to climate. In this study, we focus on one of these classes, Pothole Till, which is the largest class by area (Fig. 1), features the highest wetland pond density (Wolfe et al. 2019a), and continues to experience widespread wetland drainage. The extent of wetland ponds is temporally dynamic and influenced by both climate and anthropogenic activities (Johnson et al. 2005; Doherty et al. 2018). We used the Global Surface Water dataset (Pekel et al. 2016) to characterize baseline wetland pond extent, which identifies water bodies on the landscape according to the presence of water during the 1980–2010 period.

The foundation of our analysis is a catchment classification-based virtual modelling approach. Using this approach, we modelled the typical sensitivity of ecosystem services to drainage of wetlands in the Pothole Till watershed class under a characteristic climate for the region. This approach has been tested and described in detail for assessing hydrological sensitivity to climate and wetland drainage (Spence et al. 2022a, 2022b; He et al. 2023). The foundation of this approach is the Cold Regions Hydrological Modelling Platform (Pomeroy et al. 2007; Pomeroy et al. 2022), which was used to simulate the hydrological behaviour of a typical Pothole Till watershed. The Cold Regions Hydrological Modelling Platform is a process-based, spatially semi-distributed hydrological model, featuring key cold region and warm season hydrological processes of import in this region. In this analysis, we build on the wetland drainage–based investigation described by Spence et al. (2022b).

Briefly, the Cold Regions Hydrological Modelling Platform was configured to represent the typical 100 km² Pothole Till watershed. Hydrological response units are used in the model to represent the flow of water through the catena of land covers observed across the Pothole Till landscape (e.g., fallow fields, cultivated fields, grasslands, shrublands, woodlands, wetlands, and stream channels), and their fractional coverage was defined using the median of all watersheds in this class (e.g., Wolfe et al. 2019a). The virtual watershed model was set up to have (1) an effective area catena generating runoff directly to the stream, and (2) a non-effective area catena featuring a wetland complex composed of 46 wetlands where runoff from uplands is subject to depressional storage; routing order details for this structure are available in Spence et al. (2022b) and He et al. (2023). Forty-six wetlands were used, as this has been shown to be sufficient for representing wetland control of runoff transmission from Prairie watersheds (Pomeroy et al. 2014; Spence et al. 2022a), while also affording some parsimony. The effective and non-effective proportions of the watershed area and the wetland density (∼17 km^–2) were established according to the typical character of the Pothole Till watershed class (Wolfe et al. 2019b). A generalized pareto distribution with shape and scale parameters of 0.87 and 2227, respectively (e.g., Spence et al. 2022b), was used to characterize individual wetland pond sizes. Runoff behaviour is scaled within the Cold Regions Hydrological Modelling Platform from a smaller area (∼2.7 km⁻²) representing effective and non-effective catenas where appropriate wetland density can be represented with 46 wetlands, to the 100 km² Pothole Till watershed scale of interest.

Long-term simulations (46 years) of the virtual watershed were run with the Cold Regions Hydrological Modelling Platform using Adjusted and Homogenized Canadian Climate Data (Vincent et al. 2012) from Yorkton, Saskatchewan, which lies near the geographic centre of the Pothole Till class, as the forcing dataset. Because depressional water storage can have multi-year hydrological memory in this region, we excluded the first 5 years of the simulation, limiting our analysis to a 41-year period (1965–2005). The virtual watershed model was run for 11 different wetland drainage scenarios, from a baseline condition with no (0%) drainage to complete (100%) drainage, increasing by increments of 10%, with the wetland ponds drained by size, from smallest (low drainage scenarios) to largest (high drainage scenarios). This approach approximates drainage practices, which often target small, easy-to-drain wetland ponds first (Serran and Creed 2016). Further, Spence et al. (2022b) showed simulated hydrologic responses for the Pothole Till to be largely independent of the sequence of wetland drainage. These increments are nominal, as the wetland area drained was approximated to the nearest 10% according to the size of remaining wetland ponds available to drain in each subsequent scenario. In the model, water flows from wetlands when the depressional storage capacity is exceeded, mimicking the fill-and-spill behaviour characteristic of this landscape (Shaw et al. 2012), with spilling water flowing to downstream wetland depressions. When (former) wetland hydrologic response units are drained in one of the drainage scenarios, the depressional storage is set to zero and other parameters are changed to match those of cropland hydrologic response units, mimicking the land transformations occurring as part of agricultural wetland drainage (see also Spence et al. 2022b). Manning's n was similarly set to match that of the stream channel, on the assumption that wetland drainage leaves ditches between wetlands to rapidly convey water off the landscape.

Hydrological ecosystem service indicators included surface water storage and upland runoff. Storage of surface water in wetland ponds was quantified for each scenario by identifying the annual daily maximum water stored in each pond (where drained ponds hold no water) and summing this across all (remaining) ponds. We also investigated potential changes to runoff from the watershed associated with different wetland drainage levels, using annual return period flows calculated using a rank technique (Spence and Mengistu 2019), as flood frequencies in the Canadian Prairie do not follow a single frequency distribution (Zhang et al. 2020). We report runoff depths for return periods from 1 in 2 (mean annual flood) to a maximum of 1 in 41 years, which is the upper limit we can quantify using model simulations for the period of record.

We investigated total phosphorus export as a biogeochemical ecosystem service indicator. This was investigated using the runoff amounts described above, estimated for each year and scenario. This analysis used available total phosphorus concentration data from 14 streams draining Pothole Till watersheds in Manitoba and Saskatchewan (Cavaliere et al. unpublished data; Fasching et al. 2019); data were collected between 2011 and 2015, with total observations for each site ranging from 24 to 49. Discharge–concentration relationships have been reported for various locales outside of the Prairies (e.g., O'Connor et al. 2011; Bowes et al. 2015), but evidence of these relationships for streams in the region is limited (Casson et al. 2019; Cavaliere et al. unpublished data), and thus we assume chemostasis. Peak daily flows from these sites (∼800 000 m³ day^–1) were also comparable to but slightly higher than from our modelling (∼625 000 m³ day^–1), possibly because these sites were in some cases larger (25–500 km²) than the virtual watershed used in our analysis (100 km²). Total phosphorus concentrations in these streams had relatively low seasonal variation (median concentrations of 233 and 294 µg L^–1 for spring (March–May) and summer (June–September), respectively). Seasonal concentrations were used to calculate exports using modelled daily runoff during these periods. Winter (January–February) and fall (October–December) flows are extremely rare, so nutrient export during cold seasons is typically negligible. Nonetheless, we used spring and summer total phosphorus concentrations, respectively, and modelled daily flows to approximate the total phosphorus export associated with these small and infrequent events during winter and fall. Annual total phosphorus exports for all years and scenarios were calculated as the sum of daily exports.

We used pond riparian area as an indicator of ecosystem services provided by riparian habitat. To estimate the loss of pond riparian areas due to wetland drainage, we assume spherical ponds. As these pond areas can be quite dynamic, both within and among years, we used the 90th percentile of individual pond areas during the open-water season (April–September) across all years of each drainage scenario. Because the Cold Regions Hydrological Modelling Platform represents ponds as hydrologic response unit areas with a uniform depth, dynamic pond areas for each time step were calculated according to pond volume–area–depth relationships derived from sites in the Pothole Till watershed class (Hayashi and van der Kamp 2000) using CRHMr (Shook 2016; function: wetland_pond_area). Riparian buffer width was estimated from pond size, where width was set as the minimum of the square root of pond diameter, or 10 m (Fig. S1). Riparian area was then calculated as the difference in areas of the 90th percentile pond size plus the riparian buffer and the 90th percentile pond size.

Finally, the indicators of biodiversity ecosystem services were the numbers of breeding duck pairs and abundances of wetland-associated birds. These were calculated using statistical relationships with pond areas. Numbers of breeding pairs of five different dabbling duck species were calculated according to species-specific eqs. 1–5 from Bartzen et al. (2017):where DD, PA, and PD are conspecific duck density (6.56, 3.75, 9.79, 1.97, and 5.46, for eqs. 1–5, respectively), area of individual ponds, and pond density (10.5), respectively. The above equations were generated using 10 years of Waterfowl Breeding Population and Habitat Survey data for 1208 unique wetlands in the Canadian Prairie Pothole Region (see Bartzen et al. (2017) for details); this period spanned both historically dry and wet years on record. We used the mean pond size for each individual pond during 20 April to 10 June (which approximates the data collection periods) in each year and scenario. The number of dabbling ducks was summed across all ponds in the virtual Pothole Till watershed and reported here as population estimates on a per-area basis for the 100 km² watershed. Similarly, we were able to quantify the abundance of wetland-associated birds (which does not explicitly include ducks) using a statistical relationship derived from observational data for 300 wetland pond sites in Alberta (see Mantyka‐Pringle et al. (2019) for details on these data):where zPA represents a normalized wetland pond area (pond area z-score), calculated for each individual pond and year. This relationship was expected to be applicable to other areas of the PPR (e.g., Bidwell et al. 2014; Elliott et al. 2020); however, we note that the exact form of this relationship could vary regionally due to changes in avian species composition. As with ducks, results were calculated for each pond and year and summed to provide an estimate for the Pothole Till watershed. For all ponds and scenarios, we used the mean and standard deviation of pond sizes (late April–early June) for all years of the baseline (0%) wetland pond drainage scenario to calculate zPA from Cold Regions Hydrological Modelling Platform simulations of pond sizes. This approach was used for two reasons: (1) the pond size distribution changes across each drainage scenario (e.g., one or more pond sizes become zero at each successively higher level of drainage), making it inappropriate to use pond areas normalized according to scenario-specific results to assess change across scenarios; and (2) the data informing this relationship were collected for conditions that reasonably approximate the baseline scenario (e.g., baseline scenario pond sizes are very infrequently (∼1% of instances) larger than the upper limit of pond size (∼30 000 m²) for which bird abundance estimates are available from eqs. 1–6).

The study area spans parts of a region that suffers from numerous water security pressures (Gober and Wheater 2014), and our results highlight the added risk to water resources coincident with wetland drainage. The hydrography of the Canadian prairie region is complex and highly disconnected. Wetland drainage will fundamentally alter the hydrological behaviour. As wetlands are drained, the connected area increases (Spence et al. 2022b). With an increasing number of channels dug to remove water from the landscape, depressional water storage (e.g., in wetland ponds; Fig. 2) is lost. We illustrate here important changes in return period flows resulting from an increasingly connected landscape with wetland drainage, and increasing the water yield from small watersheds will contribute to downstream flooding, as reported extensively in the media (e.g., CBC News 2014; Saskatoon StarPhoenix 2015; Simes 2023). Flooding may occur locally at the scales we have modelled, but at larger scales, the cumulative impact of draining numerous wetland basins across multiple watersheds is also an important consideration and one that warrants further investigation. Drainage of wetlands across many watersheds of the scale we modelled here is expected to have additive and potentially synergistic effects. Our modelling framework, however, is not suited to predicting the flooded areas resulting from wetland drainage, both because our models are virtual in nature and because hydraulic models are needed to capture riverine processes that the Cold Regions Hydrological Modelling Platform is not designed for. That said, our results imply that removing 50% of the wetland area results in the 1:41 year flood occurring more than four times as often (Fig. 2b), with concomitant impacts on flooding costs.

Coincident with the increasing movement of surface water as wetlands are drained, is the increased likelihood that high nutrient export will occur when conditions are sufficiently wet. This is an important concern in this region with naturally nutrient-rich lakes and waterbodies (Quinlan et al. 2002), with already substantive issues of harmful algal blooms. Edge-of-field total phosphorus export is highly variable, but seems to centre around ∼0.25 kg ha⁻¹ for Prairie sites (Liu et al. 2021), which is higher than our watershed scale estimates herein, even under complete wetland loss (∼0.1 kg ha⁻¹). Estimates of nutrient runoff are sensitive to local runoff concentrations (our analysis uses median concentrations); hence, areas of higher (or lower) impact should be expected. Differences between edge-of-field and watershed-scale nutrient export estimates can also be attributed to more diverse land uses at the larger spatial scale considered in our analysis. The presence of grassland and shrubland areas where total phosphorus export may be lower than cropland, and other features, including wetlands with notable areal coverage at the watershed scale, can contribute to lower estimates. Ultimately, evidence of higher edge-of-field nutrient runoff estimates suggests our watershed scale estimates could be underestimated.

The magnitude of the increase in phosphorus export driven by wetland drainage, although uncertain, is clearly substantive. This increase is expected to have important ecological effects that can cascade through aquatic ecosystems. Algal blooms, while not a new phenomenon to the region, are increasingly topical and, in recent years, have disrupted drinking water treatment processes, leading to challenges in supplying major cities with adequate water supplies (Williams et al. 2019). Our results indicate that nutrient (phosphorus) export is expected to increase with higher levels of wetland drainage in all but the driest years (when conditions are too dry for the landscape to become hydrologically connected); if half of the historical complement of wetlands is lost, P export is predicted to approximately double. These estimates, however, do not account for changing fertilizer inputs and other changes associated with the drainage and cultivation of wetlands; hence, the real magnitude of change is likely to be greater than shown herein. Downstream mobilization of phosphorus associated with wetland drainage is a critically important consideration in a region where the agricultural industry has a strong focus on beneficial management practices. Importantly, many of these practices have limited efficacy in this cold, snowmelt-dominated region (Baulch et al. 2019). It remains unproven that any tools (beneficial management practices or otherwise), either alone or in concert, could be implemented in a way that can counteract what is anticipated to be a massive increase in nutrient delivery to lakes and other water bodies associated with the loss of wetlands from the landscape.

There is a direct link between wetland drainage and loss of habitat. Our predictions highlight the loss of riparian areas, dabbling ducks, and wetland-associated birds with progressive loss of areas of ponded water, showing their high sensitivity to drainage (Fig. 6). The analysis showed immediate and rapid loss of birds occurs as small wetland ponds are drained first, producing a pattern that begins to level off at moderate and high drainage levels when the wetland complex features only larger ponds, and very few of these organisms remain (Fig. 5). When considering the requirements of dabbling ducks and other wetland-associated birds, small ponds support higher bird densities (Bartzen et al. 2017), although wetland complexes composed of ponds of varying size and permanence are needed to meet the requirements of a diverse prairie waterbird community (e.g., Elliott et al. 2020). In our analysis, which removes wetlands in order of increasing size, the loss of water storage capacity in these small wetlands results in the remaining wetland ponds receiving more water. The hydrological behaviour of this system dictates that ponds that are present at higher drainage thresholds tend to be larger and deeper as more wetlands are drained, and consequently, they support lower densities of birds than prior to drainage, owing to strong and well-established relationships between pond size and bird abundances (Bidwell et al. 2014). At these higher drainage levels, where remaining wetlands are larger and more lake-like, dramatic changes in species composition are also expected (e.g., Elliott et al. 2020).

Unlike our analyses for nutrients and water, the mobility of biota can lead to interactive effects with adjacent watersheds. While our models allow for the movement of birds with changing pond area dynamics, (limited) relocation to habitat in undrained adjacent watersheds is possible but would be biased toward species tolerant of reduced habitat heterogeneity and subject to limits imposed via competition for available resources. Our models did not include these spatial scale effects, or those that could become evident as an increasing number of wetland basins are drained and converted to cropland, lengthening distances between wetland ponds. As fewer wetlands are retained in the landscape, wetland ponds become more isolated, and local complexes of wetlands needed to attract and support some bird species may become disfunctional (Fairbairn and Dinsmore 2001; Naugle et al. 2001; Blann et al. 2009). Wetland isolation might strongly affect species with relatively low mobility, such as amphibians (Lehtinen et al. 1999; Environment Canada 2012 [leopard frog]; Ruso et al. 2019). Importantly, these effects could accentuate biotic impacts and changes in community composition, especially at lower wetland retention levels.

Our focus here has been on dabbling ducks and wetland-associated birds, but it is important to acknowledge that riparian habitat loss at pond margins will have a profound impact on a wider range of biota. Wetland margins represent the only remaining natural upland cover in many intensively cropped landscapes, and can be important reservoirs for beneficial invertebrate populations, including predatory arthropods (Robinson et al. 2021) and pollinators (Vickruck et al. 2019). It is possible that margins lower contaminant exports from ponds via emerging aquatic insects (Kraus et al. 2021), as riparian buffers serve to reduce agrochemical runoff to pond water (Main et al. 2015). For agrochemical runoff reaching pothole ponds, these inundated areas also play an important role in enhancing the breakdown of pesticides via photolysis (Zeng and Arnold 2013). Aquatic insects can provide highly nutritious foods for terrestrial species such as bats and birds using these riparian areas (Hixson et al. 2015; Twining et al. 2016). Species that rely strongly on wetland ponds, such as amphibians, will be seriously impacted by wetland riparian habitat loss, altered wetland function and quality, and greater wetland isolation as drainage progresses. Even large ungulates often select for wetland habitats in prairie landscapes for shelter, food, and water (e.g., LaForge et al. 2016), so the removal of wooded riparian areas could concentrate animals in less remaining habitat and eventually produce density-related reductions in survival and reproduction (Bonenfant et al. 2009). Key resources associated with wetland ponds and adjacent riparian habitats may explain why some species can occur in intensively cropped agroecosystems (Janke et al. 2019; Nelson and Gillan 2020; Berzins et al. 2022). The potential for broad impacts on wetland riparian habitat and the associated benefits these habitats can provide to the agroecosystem, in concert with changes to water resources, suggest the need for a more comprehensive understanding of the economic pressures around wetland pond drainage.

The major impetus for draining pothole wetlands is private benefit to the landowner. Depending on their size and location, wetland ponds can reduce field trafficability, increase the time required to seed, fertilize, and harvest, and lead to higher input costs where equipment size is not ideally suited. Draining wetlands reduces both field obstruction costs and the opportunity cost of non-productive land (Cortus et al. 2011; Lawley 2014). The profitability of cultivating drained wetlands in the Canadian PPR, however, is highly variable for individual wetlands and, in some cases, leads to financial losses (Clare et al. 2021). In contrast, ecosystem services provided by pothole wetlands can be highly valuable (e.g., Dias and Belcher 2015; Pattison-Williams et al. 2018), and the impacts associated with drainage, some of which have been quantified here, can be costly to society. Cropland expansion in the United States is associated with increasingly marginal yields, and costs to wildlife are high (Lark et al. 2020). While there is inherent uncertainty around these costs and benefits (e.g., Hansen and Loesch 2017), a challenge also lies in there being a mismatch between those who bear the costs of wetland conservation and those who benefit (Dias and Belcher 2015). This can contribute to enduring conflict (Breen et al. 2018). While there are inherent uncertainties in quantifying and valuing ecosystem services, and these uncertainties can be an obstacle to assessing tradeoffs among land-use decisions (Johnson et al. 2012), there nonetheless remains a strong need for further efforts to both quantify and use economic approaches to value pothole wetlands. The Canadian prairie region features a range of approaches to wetland management, from no net loss (of wetland area) (e.g., Manitoba) to widespread permitting of wetland drainage (e.g., Saskatchewan). Given the examples herein of the loss of ecosystem services when both small and large wetland ponds are lost, these approaches can be expected to contribute to the loss of ecosystem services, but in different ways. In some instances, e.g., Saskatchewan's forthcoming agricultural water stewardship policy, economic analyses of wetland drainage consider only private benefits to landowners. Thus, there remains a strong need to comprehensively value the direct and indirect costs of wetland loss to allow for those costs to society to be pragmatically weighed against the benefits of wetland drainage that currently accrue for private landowners pursuing this practice. This approach is needed to better inform and design wetland conservation policy measures. The costs of habitat restoration and species recovery are generally far more costly and difficult than protecting existing habitat. These costs and others, including the costs of flooding and efforts to reduce nutrient loading, are rarely considered in long-term economic assessments of wetland loss resulting from agricultural expansion.

We focused herein on key biophysical ecosystem services, spanning provisioning (waterfowl), regulating (water storage, nutrient retention), and supporting (habitat) ecosystem services. These are among the biophysical services identified previously as being at risk of impact due to pothole wetland drainage (Baulch et al. 2021). While we believe that this analysis effectively examines the potential severity of the loss of ecosystem services, it stops short of quantifying other services, for example, supporting services such as carbon storage, floral resources, or amphibian habitat that have been investigated for U.S. regions of the Prairie Pothole Region (Mushet and Roth 2020). We have likewise not investigated the loss of cultural services here, including recreational experiences, and the significance of wetlands to Indigenous peoples (e.g., Clarkson et al. 2013; Department of the Environment 2016). In contrast, drained wetlands are typically replaced with cropland, which can provide other ecosystem services; however, these are widespread on this landscape and would not add to the diversity of services provided. It can be inferred that there would be at best a close to one-to-one replacement with land used to grow plants for human or animal consumption, and we have not included this potential for ecosystem service switching in our analysis.

It is worth noting that our approach is one designed to typify the behaviour of a watershed class. This is advantageous in that it allows us to regionalize behaviour on a very local scale (e.g., 100 km² is the scale of a few large farms). The tradeoff, however, is that this approach does not capture very well the potential range in loss of ecosystem services with wetland drainage for a particular watershed class. For example, we could expect that within the Pothole Till class, for (real) watersheds with an areal wetland coverage higher than that modelled here, we would see greater absolute loss of ecosystem services (e.g., water storage, total phosphorus export, bird abundance) with each additional level of drainage. The converse would be true of locations with a somewhat lower proportion of wetlands naturally present on the landscape than modelled here. Further, while we limit our analysis to a single representative climate station for the Pothole Till class, the analysis is deliberate in its use of the highest quality climate record available, which, at a length of more than 40 years, captures much of the most extreme variability observed to date. Given the inherent interannual variability in this region, which alternates between multi-year wet and dry phases, analysis of this length of record allows us to capture the dynamic behaviour of multiple wet and dry phases, including the most severe drought on record. We caution that, given the climatic variability within the region, this analysis should not be expected to illustrate the full range in watershed response to drainage for this class. The behaviour illustrated (e.g., Figs. 2 and 3) could be expected to shift somewhat should the analysis be repeated using a drier or wetter climate station lying closer to the periphery of the Pothole Till watershed class (e.g., to lower or high ranges, respectively). Nonetheless, given strong natural temporal variability from wet, well-connected conditions to dry conditions with limited runoff generation at local scales, the interannual variability shown here is likely to capture much of the spatial variability associated with the region. As such, these data should not be used to provide local scale assessments, but instead are expected to provide a robust representation of change at large spatial scales. We also rely on several assumptions in our analysis, most notably that ponds are round, to estimate riparian areas. While this will not be strictly true, it was a suitable simplification to approximate these relationships. We also acknowledge that our predictions are associated with the strength of the underlying hydrological modelling. He et al. (2023) explore this for the virtual-watershed hydrological modelling in some detail, illustrating that the modelling framework is robust. Given the importance of the hydrological behaviour of this landscape to the characterization of ecosystem services, we believe our approach, featuring robust hydrological representation, provides a strong underpinning to our work. Our focus on strong hydrological representation confers advantages over other off-the-shelf approaches that do not capture water storage and potential changes in water movement associated with the loss of wetlands and do not capture the dynamic interannual behaviour of wetlands.

In future work, our analysis could be extended to larger scale ecological models that can account for potential indirect effects. This could include, for example, the impact on ecosystems in other regions that could see fewer migrating birds as wetlands are drained. Given the large magnitude of impacts quantified through this analysis, extending this work to consider the economic value of the loss of ecosystem services should be pursued. This might include quantifying the potential for indirect economic impacts, e.g., tourism associated with duck hunting, as one example. While we have not quantified changes to carbon storage here, loss of carbon with wetland drainage is expected, given that cultivated wetland depressions have much lower carbon stores than uncultivated ones (Bedard-Haughn et al. 2006); however, the magnitude of impacts to greenhouse gas budgets for these systems has not yet been described with certainty (Baulch et al. 2021). Improved estimates of the impact on carbon storage and net climate forcing as pothole ponds are drained and converted to agriculture warrant a fuller investigation. With rising interest in maintaining natural carbon sequestration capacity, future assessments of the loss of ecosystem services using this framework should be expanded to include this aspect. We also acknowledge that our approach to nutrient export modelling is based on evidence of chemostatic phosphorus behaviour in the region and does not include explicit process representation. While available evidence suggests that water volume is the key control on phosphorus export, given evidence of higher edge-of-field total phosphorus export than reported herein at the watershed scale, further examination of how the role of wetland density, their connectivity, and water residence times in individual wetlands can affect nutrient retention in these systems should be considered in future work. Similarly, further work to understand and simulate drivers of localized variation in runoff chemistry will help refine our understanding of downstream impacts on aquatic ecosystems.

Finally, as a region undergoing strong climate change, our modelled response to wetland drainage for a period spanning the last half of the twentieth century may not reflect the consequences of wetland drainage under warmer and perhaps wetter or drier conditions. Wetland-driven evapotranspiration contributes to moisture recirculation and several degrees of cooling in wetland-rich areas (e.g., Zhang et al. 2022), with the potential that loss of these features can exacerbate hot and dry conditions and further raise water use demand (e.g., through irrigation). These interactions are important, as predictions of the change in bird abundance associated with pond conditions in different parts of the Prairie Pothole Region under future climate have suggested both decreases (e.g., Sorenson et al. 1998) and increases (McIntyre et al. 2019) are possible. Examining how wetland drainage can be expected to affect ecosystem services under future climate conditions in the region should be considered in future analyses. In this regard, we see a need to extend this analysis from the Pothole Till to other watershed classes that exhibit a range in wetland coverage (Wolfe et al. 2019a), hydrological behaviours, and sensitivity to climate (He et al. 2023).