1 Towards a Dynamic Effective Drainage Area Map for the Canadian Prairie: Sensitivity of Contributing Area to Wetland Storage Capacity Stephanie Bacsua and Christopher Spenceb* aStephanie Bacsu Consulting, Saskatoon, Saskatchewan, Canada; bEnvironment and Climate Change Canada, 11 Innovation Blvd., Saskatoon, Saskatchewan, 5 Canada, S7N 3H5, chris.spence@ec.gc.ca Accepted to the Canadian Water Resources Journal November 2023 2 TOWARDS A DYNAMIC EFFECTIVE DRAINAGE AREA MAP FOR THE CANADIAN PRAIRIE: SENSITIVITY OF CONTRIBUTING AREA TO WETLAND 10 STORAGE CAPACITY Abstract Wetlands that occupy topographic depressions are a defining feature of the Canadian Prairie. These features control hydrological connectivity as they contain high storage capacity relative to 15 typically available precipitation and runoff. Altering wetland distribution changes the frequency with which areas can become hydrologically connected to the catchment outlet. There are several methods that have proven successful in estimating contributing area response to changes in wetland storage, but none have been applied widely across the Canadian Prairie. The objective of this study was 1) to determine if the rate at which contributing area changes with loss of wetland 20 storage capacity is related to wetland distribution, and if so, 2) to map sensitivity across the region. To do so, a GIS analysis was employed in which iterative measurements were made of contributing area expansion with simulated wetland storage capacity loss. Results show that those catchments with more small wetlands have more sensitive contributing areas to storage capacity loss. Extrapolation of this relationship across the Canadian Prairie shows that areas in 25 western Manitoba and southeastern Saskatchewan are among the most sensitive. These results provide a better understanding of how contributing area may change with loss of wetland storage capacity. This may inform implementation of beneficial management practices, transportation network design, and management of streamflow and nutrient transport to downstream lakes across the region. 30 Resumé Les milieux humides occupant une dépression topographique sont trés frèquents dans les des Prairies canadiennes. Ils contrôlent la connectivité hydrologique, du fait qu’ils possèdent une 3 capacité de stockage élevée relative aux précipitations et aux ruissellements typiquement disponibles. Toute altération de la distribution des milieux humides change la fréquence à 35 laquelle les différentes zones peuvent devenir connectées hydrologiquement à la sortie du bassin versant. Plusieurs méthodes ont donné de bons résultats pour estimer la réponse des zones contributives aux changements du stockage des milieux humides, mais aucune n’a été appliquée à large échelle dans les Prairies canadiennes. L’objectif de cette étude était 1) de déterminer si le taux auquel une zone contributive change avec la perte de capacité de stockage des milieux 40 humides est lié à leur distribution, et si c’est le cas 2) de cartographier leur vulnérabilité dans la région. Pour ce faire, une analyse par système d'information géographique des mesures itératives de l’expansion des aires contributives ont été faites en utilisant une simulation de la perte de capacité de stockage des milieux humides. Les résultats montrent que les bassins versants possédant davantage de milieux humides de petites dimensions ont davantage de zones 45 contributives vulnérables à la perte de capacité de stockage. Une extrapolation l’ensemble des Prairies montre que les zones de l’ouest du Manitoba et du sud-est de la Saskatchewan font partie des plus vulnérables. Ces résultats permettent de mieux comprendre comment une zone contributive peut changer en fonction d’une perte de capacité de stockage des milieux humides. Ceci peut éclairer la mise en œuvre de meilleures pratiques de gestion la conception routière 50 ainsi que la gestion et la gestion de l’écoulement et du transport des nutriments jusqu’aux lacs situés en aval des milieux humides. Keywords: Prairie, contributing area, storage capacity, streamflow, wetlands 55 4 Introduction The Canadian Prairie is a distinct landscape in which wetlands are a prominent feature. They are unevenly distributed throughout the region both spatially and temporally. In addition, they vary 60 in size, shape, and pond permanence (Ameli and Creed, 2017; van der Kamp et al., 2016). Wetlands are typically identified based on the presence of wetland vegetation and soils, created by the presence of standing water and/or saturated conditions. It is not necessary for all wetlands to contain ponded water and ponded periods are dynamic (Pennock et al., 2014; van der Kamp et al., 2016; Van Meter and Basu, 2015). Antecedent moisture conditions, surrounding topography, 65 and precipitation inputs determine the surface water storage state in wetlands (Shaw et al., 2013; Hayashi et al. 2016). The storage state of the wetland complex controls the extent of the streamflow network and the total area that can provide runoff to the catchment outlet (Stichling and Blackwell, 1957; Golden et al. 2016). During periods of drought, surface storage deficits are common among individual wetlands, and this prevents runoff from proceeding downslope. 70 During periods of deluge storage deficits can be overwhelmed by water inputs, increasing hydrological connectivity among neighbouring wetlands. This expands the stream network and contributing area. Connectivity among wetlands and areas contributing runoff, therefore, are variable in timing and duration (Amando et al., 2018; Ameli and Creed, 2017; Hayashi et al. 2016; Shaw, Pietroniro and Martz, 2013; Shaw et al., 2012). The portion of water contributed 75 downstream is controlled by such fill-and-spill or fill-and-merge mechanisms, which are threshold mediated responses (Ehsanzadeh et al., 2012; Haque and Badiou, 2018; Leibowitz, Mushet and Newton, 2016; Spence, 2007). ’Gate-keeping’ wetlands have large storage capacity relative to runoff production capability of the upslope landscape (Baulch et al., 2021; Phillips, Spence and Pomeroy, 2011) and disproportionately control downstream streamflow response and 80 hydrological connectivity. These processes that control variability in hydrological connectivity 5 and contributing area makes the Canadian Prairie one of the most difficult landscapes in which to predict hydrologic behaviour (Baulch et al., 2021; Spence et al., 2019). Estimating contributing area in particular is difficult. It is a function of both static parameters, 85 such as topography, and dynamic variables, such as antecedent wetness and precipitation. Wilson et al. (2019) calculated median contributing areas for several Manitoba catchments assuming that there is a threshold stream density associated with the median contributing area. However, the relationship between storage (as manifested in stream density) and contributing area is both non- linear and hysteretic (Shook et al., 2013). The non-linearity can be accounted for using methods 90 such as the Simple Pothole terraIn anaLysis aLgorithm (SPILL) (Shaw, Pietroniro and Martz, 2013). The success of SPILL is subject to the resolution and accuracy of the elevation data. The low relief of the Canadian Prairie means that elevation data often need to be of fine resolution (i.e., < 1 m), with high vertical accuracy, to represent the gentle rolling topography of the region (Mengistu and Spence, 2016; Montgomery et al., 2021; Shaw et al., 2013). Unfortunately, much 95 of the region still lacks such data (Mekonnen et al., 2014). The fully distributed Wetland DEM Ponding Model (WDPM) allows for explicit representation of the storage-area hysteresis, which is parameterized in the Pothole Cascade Model (PCM) and this makes these two among the more robust approaches for estimating contributing area (Shook and Pomeroy, 2011). 100 The industry standard for representing the usual state of the contributing area on the Canadian Prairie is the effective drainage area map drafted by Agriculture Canada in the 1970s (Agriculture Canada, 1983; Godwin and Martin, 1975; Martin, 2001). The effective drainage area is a term uniquely used in the Canadian Prairie to define the portion of the gross drainage 6 area (i.e., the maximum area of a catchment delineated by the topographic divide between 105 neighbouring catchments) that contributes water to the basin outlet during a median year (Agriculture Canada, 1983; Ehsanzadeh, van der Kamp and Spence, 2016; Godwin and Martin, 1975; Martin, 2001). Hydrological and biogeochemical studies in the region have relied on the effective drainage area maps but it is generally understood that these maps have become out of date because of changes in climate that have altered the median annual flow, as well as the 110 removal of wetlands for agricultural expansion that has increased runoff efficiency and expanded the effective drainage area (Wilson et al., 2019). A variety of approaches, such as filling, ditching, draining, tiling and consolidating, have been used to remove 40-71% of wetland area in some parts of the Canadian Prairie and comparable areas in the United States (Amando et al., 2018; Baulch et al., 2021; Conly and van der Kamp, 2001; Haque and Badiou, 2018; Krapu, 115 Kumar and Borsuk, 2018; McCauley et al., 2015; Spence et al., 2022; van Meter and Basu, 2015). It is a challenging landscape for researchers, policy makers, and water managers to determine how wetland drainage has influenced catchment hydrological connectivity and streamflow response because, as noted above, neither the climate nor the landscape has remained stationary, and wetland inventories are sporadic in space and time (Badiou, Page and Boychuk, 120 2018; 2019). With the preferential loss of smaller wetlands (McCauley et al., 2015) and changes in the distance between neighbouring wetlands and river networks (van Meter and Basu, 2015) should come changes in how often areas contribute to streamflow. The methods described above to estimate contributing area all bring value and enhanced 125 understanding of the behaviour of contributing area in this region. However, there has been no wide application of the methods described above to evaluate how storage capacity loss dictates 7 contributing area across the Canadian Prairie. There is evidence that the distribution of storage in wetlands in prairie landscapes influences contributing area (Spence, 2007) but quantifying this relationship has only been done with theoretical drainage networks (Shook et al., 2021). A map 130 of sensitivity of contributing area to storage state across the region would be helpful in updating the effective drainage area maps and help expand these maps to floods of different return periods. This could identify possible flood and nutrient hotspots which could inform water management policy, identify best locations for agriculture BMP (beneficial management practice) adoption, and support infrastructure design. Therefore, the objectives of this study 135 were to 1) determine the rate at which contributing area changes as wetland storage capacity is satisfied or removed from the landscape among a sample of representative prairie basins and 2) determine if there is a relationship between this rate and antecedent wetland distribution, and if so, map sensitivity of contributing area to wetland storage capacity loss across the region. 140 Study Catchments and Methods Catchment selection and flood frequency analysis Catchment traits needed for this study included the availability of streamflow data from 1961- 1990. This period was used as it was expected to reflect wetland distribution state when the late 1970’s era non-contributing area mapping was performed (as described in the next section). 145 Streamflow needed to be unregulated. Gross drainage area was to be less than 1000 km² to minimize the chance of large withdrawals via reservoir evaporation and human consumption (i.e., agricultural, domestic and industrial use). The wetland distribution among the selected catchments needed to represent the range of wetland distributions documented across the Prairie ecozone. This range was derived from results of the catchment classification of Wolfe et al. 150 8 (2019). Wolfe et al. (2019) reported two metrics that represent wetland distribution which were adopted here; the scale (β) and shape (ξ) parameters of the Pareto distribution of wetland area (Shook et al., 2013). The scale is an index of dispersion, similar to the standard deviation. The shape is how quickly the distribution declines from the peak frequency, and a smaller value reflects tendency towards smaller wetlands in the distribution. 155 Twelve catchments gauged by the Water Survey of Canada in southern Manitoba and southeastern Saskatchewan met these requirements. From these twelve, sensitivity to the loss of wetland storage capacity was evaluated for eight (Figure 1; Table 1). The four that were removed extended beyond the geographic boundaries of the Wolfe et al. (2019) classification, 160 and so were lacking wetland distribution data, but were used to develop the regional flood curve. Daily data for these stations were obtained from the Water Survey of Canada (WSC, 2023). Zhang, Stadnyk and Burn (2020) determined that there may not be a best statistical distribution to apply to flood frequency analyses on the Canadian Prairie. Therefore, flood frequency curves were generated for each catchment to estimate the median flood (Q2) assuming both Gumbel 165 Extreme Value and Log Pearson Type III statistical distributions. The distribution that produced the best fit (i.e., highest regression coefficient) regional median annual flood - effective area curve was applied in subsequent analyses. Regression was performed using the lm function in R (Wilkinson and Rogers, 1973). The slope of this curve was used as the rate at which the Q2 may increase with effective area with the loss of wetland storage capacity. 170 9 Figure 1: Location of the eight natural non-regulated gauged basins selected for the study. The inset map is of the Canadian Prairie Provinces; Alberta (AB), Saskatchewan (SK) and Manitoba (MB). The yellow box in southern Manitoba in the inset map denotes the extent of the main map. 175 Table 1: Selected watersheds for analysis. Water Survey of Canada (WSC) hydrometric data were obtained from WSC (2023). Station Number Station Name Gross Drainage Area (km²) Original Effective Drainage Area (km²) Median annual flood (m3/s) Scale () Shape () 05LL007 Pine Creek 620.2 470.5 8.2 1765 0.9 05LL013 Boggy Creek 470.4 201.5 4.0 2643 0.6 05NG012 Elgin 500.1 221.8 7.8 1625 0.7 05MG003 Gopher 301.3 152.8 3.5 1358 0.8 05NF008 Graham 740.7 177.9 0.4 1892 0.6 05MH006 Little Souris 644.1 262.1 12.0 1374 0.9 05NG010 Oak Creek 982.7 439.6 8.0 1970 1.4 05OA008 Pembina 337.9 269.7 7.9 2005 1.1 10 Rate of effective drainage area expansion Two specific contributing area states were mapped across much of the Canadian Prairie by 180 Agriculture Canada (1983). The first state is the gross drainage area, which is what most hydrologists would define as catchment area, i.e., the plane area defined by surface topographic divides (Agriculture Canada, 1983; Martin, 2001). Gross drainage area was manually mapped on National Topographic System 1:50000 scale topographic maps by hand drawing drainage divides perpendicular to contour lines that defined the height of land. In many years much of this area 185 may not contribute streamflow to the catchment outlet because upslope runoff is intercepted and retained in depressions. The second state is the effective drainage area, which was originally mapped with planimeters on the 1:50,000 scale maps. The original hard copy maps have since been digitized (AAFC, 2023). Effective drainage area was estimated from the gross drainage area by excluding those areas draining into depressions not expected to permit the transmission 190 of runoff during a year of median wetness. Aerial photographs, interviews with landowners, and field inspections were used to corroborate topographic information (Godwin and Martin, 1975). In this study, a GIS analysis was used to estimate contributing area expansion beyond the original Agriculture Canada defined effective drainage area (Ae0) as surface storage capacity is 195 lost. There are two predominant ways in which surface storage capacity in a prairie pothole depression is overcome so that contributing area can expand. The first is satisfaction of the storage capacity with water from local snowmelt in the depression, upslope runoff, direct precipitation and groundwater discharge. The second is removal of the storage capacity via drainage (e.g., ditching, tiling or infilling). Wetlands were deleted from the dataset iteratively to 200 simulate either. Wetlands located on the periphery of the antecedent effective drainage area were 11 removed first. One stream order at a time, the area of all wetlands that drained towards that stream reach was summed (Aw) and the effective drainage area was subsequently expanded (Ae) until the gross drainage area was reached. This effective drainage area expansion can be interpreted as being due to the satisfaction of storage capacity during wet conditions, or the 205 removal of storage capacity via drainage. Stream order is a function of the scale of the underlying drainage network that is selected. The drainage network used to identify stream order was generated using ArcMap version 10.7.1 Spatial Analyst toolbox from ~30 m resolution elevation data provided by the Shuttle Radar Topography Mission DEM (SRTM). The number of stream orders present in the catchment dictated the number of iterations necessary to expand 210 the contributing area to the gross drainage area. This varied from four to seven iterations, with the example of the Little Souris illustrated in Figure 2. 12 Figure 2: The expansion of effective drainage area in the Little Souris watershed (05MH006). Each panel illustrates the effective drainage area at each iteration of wetland storage capacity 215 removal exercise. The black line denotes the gross drainage area (shown for reference in Panel A). The grey lines delineate catchments from the Wolfe et al (2019) classification. The effective drainage area is in yellow, and wetlands are shown in dark blue. The embedded text in each panel details the contributing area and the fraction of wetland area remaining after each iteration. 220 13 Wetland area was measured using the Global Surface Water dataset (Pekel et al., 2016) which indicates where water was detected at any time using Landsat satellites from 1984-2015, at a resolution of 30 m. This dataset more correctly identifies depressions that have periodic or permanent inundated area, rather than wetlands, which are more properly defined by soils and vegetation. However, the dataset is indicative of the location, size and shape of depressions 225 where water pools in this landscape. This type of product was selected for two reasons. First, it provides a measurement of the number and extent of depressions in which water can collect so that it is prevented from reaching the catchment outlet. This is an indication of one of the static properties (i.e., topography) that controls contributing area, and in turn, annual daily maximum streamflow. Second, it provides a standardized data product across the entire ecozone. It 230 represents the most contiguous wetland open source dataset that covers the Canadian Prairie. While area was used as a metric of wetland size, available volume is what contains any water that would otherwise be transmitted downstream. Methods to determine volume exist (Minke, Westbrook and van der Kamp, 2010; Hayashi and van der Kamp, 2000), however, these would have consistently adjusted the wetland metric such that the same relationships would be found, 235 but with different parameters. As such, for this exercise, wetland area was chosen as the focus. Regional sensitivity to wetland storage capacity removal The scale (β) and shape (ξ) parameters of the Pareto distribution of wetland area were calculated for gross and effective drainage areas for the eight catchments. These data were extracted from 240 those calculated by Wolfe et al. (2019) for every 100 km2 catchment across the Prairie ecozone using data from the Global Surface Water dataset (Pekel et al., 2016). As three forms often found when relating streamflow and topographic parameters, least squares linear, power and 14 logarithmic relationships between  and  and the change in both contributing area and median annual streamflow were derived with the R package lm (Wilkinson and Rogers, 1973). These 245 were used to determine the best shape of the relationship between antecedent wetland distribution and the sensitivity of contributing area and streamflow to wetland storage capacity loss. The relationship with the best regression coefficient was used to extrapolate sensitivity across the catchments included in the Wolfe et al. (2019) classification to measure regional sensitivity of contributing area and streamflow to wetland storage deficit loss. While there may 250 be other landscape characteristics that control contributing area, the focus here was on determining if wetland distribution was significant. Results Rate of effective drainage expansion 255 The rate at which contributing area increased with loss of wetland storage capacity varied among study catchments (Table 2; Figure 3). The highest rate of increase was in the Gopher catchment in which effective drainage area increased 30.5 km2 for every km2 of wetland area with lost storage capacity. This contrasts with the Pembina River, with a rate of increase an order of magnitude smaller at 3.8 km2 / km2. Similar slow rates of effective drainage area expansion in 260 the Oak and Pembina were not a function of the original effective drainage area, as the original Oak effective drainage area (439.6 km2) approached twice the size of the Pembina (269.7 km2). The Oak and Pine were among the catchments (Figure 3) with the largest original effective drainage areas, but they increased in size at quite different rates, 7.4 km2/ km2 and 23.7 km2/ km2, respectively. The original Graham effective drainage area was among the smallest of the 265 eight catchments, but the Graham expanded at a rate similar to that of the Oak, which had the 15 second largest original effective drainage area. The rate of effective drainage area expansion was not a function of the antecedent effective drainage area. Table 2: Antecedent (Ae0) and final (Aex) effective drainage areas, antecedent (Aw0) and final (Awx) wetland areas, and the rate of change in effective drainage area and median annual flood with 270 loss of surface water storage capacity. Station name Ae0 (km²) Aex (km²) Aw0 (km²) Awx (km²) Ae/Aw (km2/km2) Q2/Aw (m3/s/ha) Pine 470.5 508.1 4.7 3.4 23.7 0.6 Boggy 201.5 336.7 10.3 2.1 15.9 0.3 Elgin 221.8 384.2 9.8 0.3 17.2 0.3 Gopher 152.8 217.5 3.1 0.8 30.5 0.6 Graham 177.9 693.9 12.4 2.0 47.9 1.0 Little Souris 262.1 473.7 22.3 1.4 8.8 0.2 Oak Creek 439.6 746.7 62.1 16.2 7.4 0.1 Pembina 269.7 296.5 11.1 4.2 3.8 0.1 Figure 3: Scattergram to illustrate the increase in effective drainage as wetland storage capacity is lost for the eight study catchments. 275 16 A power relationship exhibited the best fit between rate of change in effective drainage area with loss of wetland storage capacity (Table 3). All three relationships performed comparably. The logarithmic relationship (Table 3; Figure 4) was selected for the extrapolation exercise: ln((∆𝐴𝑒)/(∆𝐴𝑤 ))=-1.83∙ln()+ln(2.34) (1) 280 where Ae is change in effective drainage area and Aw is the change in wetland area with storage capacity loss. The standard error associated with Eq. (1) implies there is 95% certainty that estimates of the rate of change in effective drainage area are ±2 km2/km2. The reason Eq. 1 was selected over the power function is because when the latter was applied across the Prairie ecozone, the extreme non-linear nature of the relationship produced unrealistically high estimates 285 of the rate of effective drainage area expansion in catchments with low  values (< 0.4). While all eight catchments had skewed values of  closer to zero (Figure 5), the form of Eq. (1) implies it is those catchments with smaller values of  and higher numbers of small wetlands that are more sensitive to changes in effective contributing area with loss of storage capacity on the landscape. Those catchments with more large wetlands will have slower rates of change in 290 effective drainage area and than those with more small wetlands. Table 3: Comparison of relationships between wetland distribution and the rate of change in effective drainage area and loss of surface water storage capacity. The terms “m” and “b” are the multiplier and intercept, respectively in the linear and logarithmic equations and the multipler and exponent in the power function. 295 Relationship form r2 p m b Linear 0.48 0.06 -0.81 1.14 Power 0.64 0.04 0.23 -2.214 Logarithmic 0.54 0.04 -2.16 2.4 17 Figure 4: The relationship between the natural log of the Pareto wetland area distribution shape parameter (ln()) and the rate at which effective drainage area increases with the area of wetland storage capacity loss (ln(Ae/Aw)). The black dots denote each of the study basins. The blue 300 line is the line of best fit, and the gray area its 95% confidence limits. Figure 5: Kernel densities of wetland area for each of the gauged basins selected for the study. 18 Rate of increase in median annual flood The regional median flood frequency curve shown in Figure 6 can be represented by: 305 𝑄2 = 0.02 ∙ 𝐴𝑒 + 1.8293 (2) The regional relationship between effective area (Ae; km2) and the median annual flood (Q2; m3/s) is significant (p=0.015) with an r2 of 0.46. Confidence intervals evident in Figure 6 show the relationship exhibits less uncertainty when effective drainage areas are between 150 and 300 km2. Applying Eq. 2 to estimate the change in median annual flood with expansion in 310 contributing area suggests a very consistent increase of ~2 m3/s for every km2 of additional contributing area. The rate at which the contributing area increases as wetland surface storage capacity is lost varies widely among the eight catchments (Figure 3; Table 2). The Graham catchment has the most sensitive median annual flood with an increase of 1 m3/s for every km2 of wetland storage capacity that is removed or satisfied. This rate is an order of magnitude higher 315 than that in the Pembina catchment. This range demonstrates the high variability in sensitivity that is controlled by wetland distribution. Figure 6: Regional median annual flood frequency curve derived from the selected gauged catchments. The black dots denote basins. The blue line is the line of best fit, and the gray area 320 its 95% confidence limits. 19 Regional sensitivity to wetland storage capacity loss Combining Eqs. 1 and 2 can permit the calculation of Q2 as wetland storage capacity is satisfied or removed. Catchments in the Wolfe et al. (2019) classification all had gross drainage areas of 325 approximately 100 km2. There is 95% certainty that the best fit line between effective drainage area and median annual flood ranges between 1.25 and 7 m3/s at an area of ~100 km2. The uppermost range of uncertainty in Eq.1 equates to approximately ±7.5 km2/km2. Because of the uncertainty in the regional streamflow:area curve (Eq.2) in addition to that associated with predicted rate of change in effective drainage area with change in wetland area (Eq. 1) it was 330 considered inappropriate to extrapolate sensitivity of mean annual streamflow to wetland storage capacity loss. Only the effect on contributing area was mapped. Figure 7 illustrates the results of the extrapolation exercise applying Eq. 1 to each of the 4000 catchments evaluated in the Wolfe et al. (2019) analysis. The most sensitive areas are in the eastern portions of the Prairie, particularly portions of southwest Manitoba and east central Saskatchewan. In these areas 335 storage capacity loss can expand contributing area more than 20 km2 / km2 of wetland area. This does not mean other areas are not sensitive, but that sensitivity is not as intense as in these two areas. 20 Figure 7: Map illustrating the rate of change in effective drainage area with loss of wetland 340 storage capacity across the Prairie ecozone. The three provinces are labeled; Alberta (AB), Saskatchewan (SK) and Manitoba (MB). Discussion There are a few key assumptions made during the study that have implications for the 345 interpretation of the results. A limitation of the study is there is the inexact overlap among the periods used to collect the Global Surface Water dataset, the hydrometric data, and the original effective drainage area mapping. Another limitation is the 30 m resolution of the Global Surface Water dataset. Dictated by the relatively coarse resolution of Landsat imagery, the dataset does not include the smallest wetlands known to be present in the Canadian Prairie. Because every 350 iteration of effective drainage area (Ae) is influenced by this lack of small wetlands, and that all 21 the small wetlands would have been removed following the methodology, the overall estimate of Ae/Aw should be the same with or without the small wetlands. The direction of any bias in the estimate during a specific iteration of wetland storage capacity loss will depend on the spatial distribution of all wetlands through the catchment. 355 The study methodology used in this work does not properly account for the influence of individual wetland spatial extent on controlling incremental changes in contributing area. The original effective drainage area mapping exercise identified those depressions that acted as gatekeepers (Phillips, Spence and Pomeroy, 2011) during upland runoff conditions less than the 360 median annual flood. Here, it was assumed the entire wetland complex within each incremental stream order exhibited gatekeeping behaviour. The estimated incremental rate of change, Ae/Aw, would be an underestimate if a gatekeeper were in a specific iteration of wetland removal. Conversely, if several wetlands were removed, but none had gatekeeping abilities, Ae/Aw, would be an overestimate. It is difficult to estimate a minimum wetland size for a 365 gatekeeper, as this depends on the size relative to the upland drainage area, and its relative location is also important. As the iterations of wetland removal included the entire range of catchment areas from the effective to the gross drainage area (Figure 2), it may be that the overall rate Ae/Aw reported in Table 2 remains robust. It is likely the relationships in Figure 3 are less linear than implied there. 370 The rates of change in effective drainage area presented in Table 2 do not address the known hysteretic relationships between contributing area with antecedent wetness conditions on the Canadian Prairie (Mengistu and Spence 2016; Shook and Pomeroy, 2011). Larger basins, such as 22 those selected here, may be more amenable to assuming negligible hysteresis (Shook, Papalexiou 375 and Pomeroy, 2021). Future work could build on what has been reported here and include the role of hysteresis in estimating both expansion and contraction of contributing area in this landscape to attain a truly dynamic effective drainage area map for the Canadian Prairie. This research examined the relationship between wetland landscape metrics and streamflow 380 response to determine those types of catchments that are more responsive to wetland storage capacity loss. The shape parameter, , from Wolfe et al. (2019) proved indicative of catchment sensitivity. The results demonstrate that a catchment with a greater number of small wetlands will demonstrate an increased sensitivity to wetland storage capacity loss. This is counterintuitive to the results summarized in Shaw et al. (2012) who indicated that few larger 385 wetlands are more effective than many smaller wetlands for changing contributing area. However, while the storage capacity of an individual small wetland is less than a similar large wetland, many more small wetlands on the landscape results in many more possible connections. It may be catchments that have enhanced structural connectivity because of the number of higher possible hydrological connections (Ambroise, 2004) will have more sensitive effective drainage 390 areas. How functional connectivity manifests at any given time, however, will be dictated by rates of evapotranspiration, groundwater recharge and upslope runoff. The idea that those catchments with more small wetlands are more sensitive is supported by Shook et al. (2013). They applied the Wetland Digital Ponding Model (WDPM) to DEMs of 395 three catchments across the Prairie (Vermillion, St. Denis and Smith Creek) to calculate the drainage area of each wetland. All three sites exhibited power law relationships between 23 wetland and upland area with exponents less than one. These values of less than one imply drainage area increases more slowly than wetland area. This means larger wetlands tend to have contributing areas that are smaller relative to the depression. Small depressions have relatively 400 larger contributing areas. Shook et al. (2013) found the pattern holds for Prairie lakes across Saskatchewan, Manitoba and Alberta. This type of scaling relationship is common for fractal objects, of which Prairie lakes and wetlands are (Mandelbrot, 1982; Minke et al., 2010). The implication is that satisfying the storage capacity of small wetlands increases catchment contributing area faster than satisfying the storage capacity of large wetlands, and those 405 catchments with more small wetlands have contributing areas that will more readily expand with the loss of wetland storage capacity. This study starts to fill a knowledge gap for quantifying how effective drainage area may change with land cover across the Canadian Prairie. The Agriculture and Agri-food Canada effective 410 drainage area map has remained useful as the best available resource to interpret how dynamic contributing areas influence hydrologic response. However, it is known to be outdated in many places because of anthropogenic change (i.e., expansion of drainage) and climatic influences (i.e., trends in rainfall frequency and magnitude). Demonstrating the sensitivity of an area to wetland storage capacity loss may be the first step towards updated effective drainage area and 415 dynamic contributing area mapping. Contributing area and streamflow will respond to landscape alteration and hydrological condition variation across Canadian Prairie. Place matters, as some catchments exhibit faster rates of change with loss of surface storage capacity. Ali and English (2019) also show that there are 420 24 some areas of the Prairie ecozone that preferentially connect hydrologically. The impacts of wetland drainage (as summarized in Baulch et al., 2021) may be more significant in some areas than others (Figure 7). Even without land cover changes, a change in climate will alter the median annual flood. Expected changes in temperature and precipitation across the region (Li et al., 2019) will be felt differently as some areas are more sensitive than others. Policies, programs 425 and practices should be implemented cognizant of this. For instance, it has been estimated that to meet nutrient targets for Lake Winnipeg (Manitoba Agriculture and Resource Development, 2020) streamflow to the lake will need to be reduced by 10% (Zhang and Rao, 2012). Those regions identified as most sensitive to changes in effective drainage area in Figure 7 may be places to focus efforts. These may be catchments where wetlands should be retained on the 430 landscape, as they have a disproportionate societal benefit for flood mitigation and nutrient transport and a focus for encouraging implementation of wetland retention policies and opportunities for wetland retention. Wetlands in catchments that are less sensitive to wetland loss still provide value. The information in Figure 7 may provide information on where to prioritize efforts and as such could be a valuable resource for researchers, policy makers, and 435 catchment managers to weigh decision making. Conclusions This research examined the rate of contributing area change with wetland storage capacity loss in eight natural non-regulated gauged basins in southern Manitoba and southeastern Saskatchewan. 440 There was wide variation among the eight catchments in sensitivity of contributing area to storage capacity loss. This permitted a relationship to be found between this rate of change and antecedent wetland distribution. The shape parameter of the Pareto distribution, known to 25 represent wetland area distribution, was found to be a good predictor of this sensitivity. Extrapolating across the region using this predictor revealed spatial variability in sensitivity 445 across the region. Swaths across southern Manitoba and east central Saskatchewan are most sensitive to changes in contributing area with surface storage capacity loss. These regions are perhaps areas where efforts to mitigate regional floods and nutrient transport with wetland retention or restoration will be most effective. The research and data products produced here represent only the first step towards a dynamic contributing area map for the Canadian Prairie. 450 The methods could be improved by developing ways to account for surface water storage – contributing area hysteresis, improving access to publicly available spatially widespread high resolution elevation data, and developing wetland inventories. It is hoped that in the interim, identifying probable rates at which effective drainage areas can change with storage capacity loss will improve the ability of Prairie communities to make water and land management decisions. 455 Acknowledgments The authors wish to thank Kim Hodge at Agriculture and Agri-Food Canada for his leadership developing the portion of the Eastern Prairies Living Labs program of Agriculture and Agri- Food Canada that funded this research. He also provided some very helpful comments on an 460 earlier version of this manuscript. The authors are grateful to Environment and Climate Change Canada for supporting this research and acknowledge the role of the Global Water Futures Program Prairie Water in developing the Prairie catchment classification system, which was a crucial resource for this research. Thank you to two anonymous reviewers who helped improve the manuscript. 465 26 Data Availability Statement Data are available from the authors upon request. References 470 AAFC, 2023. Gross and effective drainage area boundaries of the AAFC Watersheds project – 2013. https://open.canada.ca/data/en/dataset/063ee9b6-b3f2-45ab-9bed-d330880064d5, date last accessed November 16, 2023. Agriculture Canada. 1983. The determination of gross and effective drainage areas in the Prairie provinces, Hydrology Report #104. Prairie Farm Rehabilitation Administration, 475 Engineering Branch, Regina, 23 pp. Ali, G. and C. English (2019). "Phytoplankton blooms in Lake Winnipeg linked to selective water-gatekeeper connectivity." Scientific Reports 9(1): 8395. Ali, G., Wilson, H., Elliott, J., Penner, A., Haque, A., Ross, C., and Rabie, M. 2017. “Phosphorus export dynamics and hydrobiogeochemical controls across gradients of scale, topography 480 and human impact.” Hydrological Processes 31(18): 3130-3145. Doi:10.1002/hyp.11258. Ambroise, B. 2004. “Variable ‘active’ versus ‘contributing’ areas or periods: a necessary distinction”. Hydrological Processes 18:1149–1155. Amando, A. A., Politano, M., Schilling, K., and Weber, L. 2018. “Investigating hydrologic connectivity of a drained Prairie Pothole Region wetland complex using a fully 485 integrated, physically-based model.” Wetlands 38(2): 33-245. Doi:10.1007/s13157-016- 0800-5. Ameli, A. A., and Creed, I. F. 2017. “Quantifying hydrologic connectivity of wetlands to surface water systems”. Hydrology and Earth System Sciences, 21(3): 1791-1808. Doi:10.5194/hess-21-1791-2017. 490 Badiou, P, B. Page and L. Boychuk. 2018. Souris River Watershed Wetland Inventory and Change Detection: Estimating the effects of wetland distribution and loss on water quality and quantity in a large prairie watershed. A report prepared for the USEPA – Grant Number: CD-96805501-3. Institute for Wetland and Waterfowl Research, Ducks Unlimited Canada, 171 pp. 495 Badiou, P, B. Page and L. Boychuk. 2019. Nutrient Export from Camrose Creek: A Moderately Drained Watershed in Alberta. Institute for Wetland and Waterfowl Research, Ducks Unlimited Canada, 30 pp. Baulch, H., Whitfield, C., Wolfe, J., Basu, N., Bedard-Haughn, A., Belcher, K., Clark, R., et al. 2021. “Synthesis of science: findings on Canadian Prairie wetland drainage.” Canadian 500 Water Resources Journal 46(4): 229-241. Doi:10.1080/07011784.2021.1973911. https://open.canada.ca/data/en/dataset/063ee9b6-b3f2-45ab-9bed-d330880064d5 27 Brunet, N. N., and Westbrook, C. J. 2012. “Wetland drainage in the Canadian prairies: Nutrient, salt and bacteria characteristics.” Agriculture, ecosystems and environment, 146(1): 1-12. Doi:10.1016/j.agee.2011.09.010. Conly, F. M., and Van der Kamp, G. 2001. “Monitoring the hydrology of Canadian prairie 505 wetlands to detect the effects of climate change and land use changes.” Environmental Monitoring and Assessment, 67(1-2): 195-215. Doi:10.1023/A:1006486607040. Ehsanzadeh, E., van der Kamp, G., and Spence, C. 2016. “On the changes in long-term streamflow regimes in the North American Prairies.” Hydrological Sciences Journal, 61(1): 64-78. Doi:10.1080/02626667.2014.967249. 510 Ehsanzadeh, E., Spence, C., van der Kamp, G., and McConkey, B. 2012. “On the behaviour of dynamic contributing areas and flood frequency curves in North American Prairie watersheds.” Journal of Hydrology, 414: 364-373. Doi:10.1016/j.jhydrol.2011.11.007. Godwin, R. B., and F. R. J. Martin. 1975. Calculation of gross and effective drainage areas for the Prairie Provinces. In Canadian hydrology symposium – 1975 Proceedings, 11–14 515 August 1975, Winnipeg, Manitoba, 219–223. Associate Committee on Hydrology, National Research Council of Canada, Ottawa. Golden, H.E., Sander, H. A., Lane, C.R., Zhao, C., Price, K., D’amico, E. and Christensen, J.R. 2016. “Relative effects of geographically isolated wetlands on streamflow: a watershed- scale analysis.” Ecohydrology, 9, 21-38. 520 Haque, A., Ali, G., and Badiou, P. 2018. “Hydrological dynamics of prairie pothole wetlands: Dominant processes and landscape controls under contrasted conditions.” Hydrological Processes, 32(15): 2405-2422. Doi:10.1002/hyp.13173. Hayashi, M. and van der Kamp, G. 2000. Simple equations to represent the volume-area-depth relations of shallow wetlands in small topographic depressions. Journal of Hydrology, 525 237: 74-85. Hayashi, M., van der Kamp, G., and Rosenberry, D. O. 2016. “Hydrology of prairie wetlands: understanding the integrated surface-water and groundwater processes.” Wetlands, 36(2): 237-254. Doi:10.1007/s13157-016-0797-9. Krapu, C., Kumar, M., and Borsuk, M. 2018. “Identifying wetland consolidation using remote 530 sensing in the North Dakota Prairie Pothole Region.” Water Resources Research, 54(10): 7478-7494. Doi:10.1029/2018WR023338. Leibowitz, S. G., Mushet, D. M., and Newton, W. E. 2016. “Intermittent surface water connectivity: fill and spill vs. fill and merge dynamics.” Wetlands, 36(2): 323-342. Doi:10.1007/s13157-016-0830-z. 535 Lim Y., Li, Z., Zhang, Z., Chen, L., Kurkute, S., Scaff, L. and Pan, X. 2019. “High resolution regional climate modeling and projection over western Canada using a weather research 28 forecasting model with a pseudo-global warming approach”. Hydrology and Earth System Sciences 23: 4635-4659. Mandelbrot B. 1982. The Fractal Geometry of Nature.W. H. Freeman and Company: New York, 540 New York, USA. Manitoba Agriculture and Resource Development. 2020. Lake Winnipeg: Nutrients and Loads Status Report 2020. Winnipeg, 12 pp. Martin, F. 2001. The determination of gross and effective drainage areas in the Prairie provinces, Addendum #8, Prairie Farm Rehabilitation Administration, Regina, 109 pp. 545 McCauley, L. A., Anteau, M. J., van der Burg, M. P., and Wiltermuth, M. T. 2015. “Land use and wetland drainage affect water levels and dynamics of remaining wetlands.” Ecosphere, 6(6): 1-22. Doi:10.1890/ES14-00494.1. McKenna, O. P., D. M. Mushet, D. O. Rosenberry and J. W. LaBaugh 2017. “Evidence for a climate-induced ecohydrological state shift in wetland ecosystems of the southern Prairie 550 Pothole Region.” Climatic Change 145(3): 273-287. Mekonnen, M. A., Wheater, H. S., Ireson, A. M., Spence, C., Davison, B., and Pietroniro, A. 2014. “Towards an improved land surface scheme for prairie landscapes.” Journal of Hydrology, 511: 105-116. Doi:10.1016/j.jhydrol.2014.01.020. Mengistu, S. G., and Spence, C. 2016. “Testing the ability of a semidistributed hydrological 555 model to simulate contributing area.” Water Resources Research, 52(6): 4399-4415. Doi:10.1002/2016WR018760. Minke, A.G., Westbrook, C..J, and van der Kamp, G. 2010. “Simplified Volume-Area-Depth Method for Estimating Water Storage of Prairie Potholes”. Wetlands 30: 541–551. Doi:10.1007/s13157-010-0044-8 560 Montgomery, J., Mahoney, C., Brisco, B., Boychuk, L., Cobbaert, D., and Hopkinson, C. 2021. “Remote Sensing of Wetlands in the Prairie Pothole Region of North America.” Remote Sensing, 13(19), 3878. https://doi.org/10.3390/rs13193878. Pekel, J. F., Cottam, A., Gorelick, N., and Belward, A. S. 2016. “High-resolution mapping of global surface water and its long-term changes.” Nature, 540(7633): 418-422. 565 Doi:10.1038/nature20584. Pennock, D., Bedard-Haughn, A., Kiss, J., and van der Kamp, G. 2014. “Application of hydropedology to predictive mapping of wetland soils in the Canadian Prairie Pothole Region.” Geoderma, 235: 199-211. Doi:10.1016/j.geoderma.2014.07.008. Phillips, R.W., C. Spence and J.W. Pomeroy. 2011. “Connectivity and Runoff Dynamics in 570 Heterogeneous Basins.” Hydrological Processes, 25: 3061-3075. Doi:10.1002/hyp.8123. https://doi/ 29 Shaw, D. A., Vanderkamp, G., Conly, F. M., Pietroniro, A., and Martz, L. 2012. “The fill–spill hydrology of prairie wetland complexes during drought and deluge.” Hydrological Processes, 26(20): 3147-3156. Doi:10.1002/hyp.8390. Shaw, D. A., Pietroniro, A., and Martz, L. W. 2013. “Topographic analysis for the prairie 575 pothole region of Western Canada.” Hydrological Processes, 27(22): 3105-3114. Doi:10.1002/hyp.9409. Shook, K. R. and Pomeroy, J. W. 2011. “Memory effects of depressional storage in Northern Prairie hydrology.” Hydrological Processes, 25(25): 3890-3898. Doi:10.1002/hyp.8381. Shook, K., Pomeroy, J. W., Spence, C. and Boychuk, L. 2013. “Storage dynamics simulations in 580 prairie wetland hydrology models: evaluation and parameterization.” Hydrological Processes, 27(13): 1875-1889. Doi:10.1002/hyp.9867. Shook, K., Papalexiou, S. and Pomeroy, J.W. 2021 “Quantifying the effects of Prairie depressional storage complexes on drainage basin connectivity”. Journal of Hydrology 593: 125846. 585 Spence, C. 2007. “On the relation between dynamic storage and runoff: A discussion on thresholds, efficiency, and function.” Water Resources Research, 43(12). Doi:10.1029/2006WR005645. Spence, C., Wolfe, J.D., Whitfield, C.J., Baulch, H.M., Basu, N.B., Bedard-Haughn, A.K., Belcher, K.W., et al. 2019. “Prairie water: a global water futures project to enhance the 590 resilience of prairie communities through sustainable water management.” Canadian Water Resources Journal, 44(2): 115-126. Doi:10.1080/07011784.2018.1527256. Spence, C., He, Z., Shook, K. R., Mekonnen, B. A., Pomeroy, J. W., Whitfield, C. J. and Wolfe, J. D. 2022. Assessing hydrological sensitivity of grassland basins in the Canadian Prairies to climate using a basin classification–based virtual modelling approach. Hydrology and 595 Earth System Sciences, 26: 1801-1819. Doi:10.5194/hess-26-1801-2022. Stichling, W., and S. R. Blackwell. 1957. Drainage area as a hydrologic factor on the glaciated Canadian Prairies. IUGG Proceedings, Volume 111, Toronto, Ontario, 12 pp. van der Kamp, G., Hayashi, M., Bedard-Haughn, A., and Pennock, D. 2016. “Prairie pothole wetlands–suggestions for practical and objective definitions and 600 terminology.” Wetlands, 36(2): 229-235. Doi:10.1007/s13157-016-0809-9. Van Meter, K. J. and Basu, N. B. 2015. “Signatures of human impact: size distributions and spatial organization of wetlands in the Prairie Pothole landscape.” Ecological Applications, 25(2): 451-465. Doi:10.1890/14-0662.1. Wilkinson, G. N. and Rogers, C. E. 1973. “Symbolic descriptions of factorial models for analysis 605 of variance”. Applied Statistics, 22, 392–399. Doi: 10.2307/2346786. Wilson, H. F., Casson, N. J., Glenn, A. J., Badiou, P. and Boychuk, L. 2019. “Landscape controls on nutrient export during snowmelt and an extreme rainfall runoff event in 30 northern agricultural watersheds.” Journal of environmental quality, 48(4): 841-849. Doi:10.2134/jeq2018.07.0278. 610 Wolfe, J. D., Shook, K. R., Spence, C., and Whitfield, C. J. 2019. “A watershed classification approach that looks beyond hydrology: application to a semi-arid, agricultural region in Canada.” Hydrology and Earth System Sciences, 23(9): 3945-3967. Doi:10.5194/hess-23- 3945-2019. WSC, 2023. https://wateroffice.ec.gc.ca/index_e.html; last date accessed November 16, 2023. 615 Zhang, W. and Y.R Rao. 2012. Application of a eutrophication model for assessing water quality in Lake Winnipeg. Journal of Great Lakes Research, 38, 158-173. Zhang, Z., Stadnyk, T. A., and Burn, D. H. 2020. “Identification of a preferred statistical distribution for at-site flood frequency analysis in Canada.” Canadian Water Resources Journal, 45(1): 43-58. Doi:10.1080/07011784.2019.1691942. 620 https://wateroffice.ec.gc.ca/index_e.html