A multi-Criteria Wetland Suitability Index for Restoration across Ontario’s Mixedwood Plains

: Signiﬁcant wetland loss (~72%; 1.4 million hectares) in the Province of Ontario, Canada, has resulted in damage to important ecosystem services that mitigate the e ﬀ ects of global change. In response, major agencies have set goals to halt this loss and work to restore wetlands to varying degrees of function and area. To aid those agencies, this study was guided by four research questions: ( i ) Which physical and ecological landscape criteria represent high suitability for wetland reconstruction? ( ii ) Of common wetland suitability metrics, which are most important? ( iii ) Can a multi-criteria wetland suitability index (WSI) e ﬀ ectively locate high and low wetland suitability across the Ontario Mixedwood Plains Ecozone? ( iv ) How do best sites from the WSI compare and contrast to both inventories of presettlement wetlands and current existing wetlands? The WSI was created based on seven criteria, normalized from 0 (low suitability) to 10 (high suitability), and illustrated through a weighted composite raster. Using an Analytical Hierarchy Process (AHP) and importance determined from a scoping review of relevant literature, soil drainage had the greatest meaning and weight within the WSI (48.2%). The Getis-Ord Gi* index charted statistically signiﬁcant “hot spots” and “cold spots” of wetland suitability. Last, the overlay analysis revealed greater similarity between high suitability sites and presettlement wetlands supporting the severity of historic wetland cannibalization. In sum, this transferable modeling approach to regional wetland restoration provides a prioritization tool for improving ecological connectivity, services, and resilience.


Introduction
With growing intensity and speed, the Earth is undergoing a massive environmental transformation best described as global change [1,2]. Attributed to anthropogenic factors such as resource extraction, growing infrastructure footprint, and land conversion, global change can largely be summarized into two overarching categories: land-cover change and climate change [3][4][5][6][7]. The combination of climate change, land-cover change, and in turn, habitat loss, is largely understood to be a major driver in decreasing biodiversity worldwide [8][9][10]. Consequently, biodiversity is recognized as an Wetland loss in Ontario is the result of many small planning decisions in favor of agriculture, or industrial and residential development over these natural landscape features [19,31,32,34]. A larger-scale perspective of Ontario's wetland reconstruction potential could help to facilitate planning for wetland projects provincewide. Around the world, wetland indices have been used to assess wetland loss and understand where wetland restoration may be most valuable and viable; however, such an index does not yet exist for Southern Ontario. This research aims to create a multi-criteria index that prioritizes areas best suited for wetland reconstruction in Ontario's Mixedwood Plains. Therefore, to aid major agencies in their macroscale wetland restoration efforts, this study was guided by four research questions: (i) Which physical and ecological landscape criteria represent high suitability for wetland reconstruction? (ii) Of the common wetland suitability metrics, which are most important? (iii) Can a multi-criteria wetland suitability index (WSI) effectively locate high and low wetland suitability across the Ontario Mixedwood Plains Ecozone? (iv) How do best sites from the WSI compare and contrast to both inventories of presettlement wetlands and current existing wetlands? An overarching goal of this paper was to deliver landscape planners, ecological restoration scientists, regional planners, and environmental managers and decision-makers an applied example for systematically evaluating wetland suitability across space in temperate climate zones.

Study Area
The Mixedwood Plains Ecozone within the Province of Ontario defines the study area boundary for this research ( Figure 1). For clarification, in 1976 the Canada Committee on Ecological Land Classification (CCELC) defined and delineated differing ecological (biophysical) areas into a multi-scale land management system based on abiotic and biotic characteristics [49]. This classification serves as a systematic resource for sustainable planning across Canada; ecozones are the largest of the environmental management classifications. This study area excludes the portion of the Ecozone crossing into Quebec and the Manitoulin Island region of Ontario, yet it is still viewed as one of the most valuable regions when considering its environmental and economic importance. Collectively, the entire Mixedwood Plains Ecozone makes up about 1.2% of Canada's total landmass [34]; however this study's area equates to roughly 0.8% of Canada's governing territory. Specifically, 81,134 km 2 (68%) of the entire 118,870 km 2 Mixedwood Plains Ecozone [34] is incorporated into this study. The Ecozone is bound between three Great Lakes: Erie, Huron, and Ontario. More precisely, it is located in Southern Ontario below the Hudson Bay Lowlands and Ontario Shield ecozones and spans across the southernmost regions of the province, including the St. Lawrence River Valley and the postglacial floodplains of the three aforementioned Laurentian Great Lakes [34,50]. The Ecozone is characterized by a temperate climate of warm summers and cool winters. The mostly glacial-shaped terrain is flat to gently rolling and includes an extensive network of waterways, surface and groundwater features, moraines, and the Niagara Escarpment [34,51]. Wetlands in the Ontario portion of the Mixedwood Plains are considered critically important landscape features, with six recognized as wetlands of international significance by the Ramsar Convention: Long Point National Wildlife Area, St. Clair National Wildlife Area, Point Pelee National Park, Mer Bleue Conservation Area, Matchedash Bay Provincial Wildlife Area, and Minesing Swamp [19]. Last, with wetlands playing a dominate role, the ecosystem goods and services provided by the Mixedwood Plains Ecozone have been conservatively estimated at $84 billion (CAD) per year [34].
The Mixedwood Plains is the most populous and industrious ecozone in Canada; it is considered the second most fragmented ecozone in Canada [52]. Thirteen major urban areas are found within this ecozone, including the first, second, and fourth most populated Canadian cities of Toronto, Montreal, and Ottawa, respectively. Between 1971 and 2006, the human population of the Mixedwood Plains increased from approximately 11 million people to more than 16 million people; the rapid population growth was projected to increase by 30% by 2031 [53]. As a result, the Mixedwood Plains has seen the largest loss in wetlands, which in turn impacts the largest populations of Ontarians. Despite suffering a significant loss in wetland area (roughly~72%) since presettlement, Ontario remains home to 6% of the Earth's total wetland area. For comparison, the conterminous United States has experienced an estimated loss of 53% of its presettlement wetlands, contributing to the large-scale erasure of wetlands and restoration needs in North America [54]. Specific to the Mixedwood Plains Ecozone, the wetland land cover classification was estimated at 25% of the total area presettlement and reduced to 7% by 2002 [32]. Wetlands have been expensed in favor of other types of land use due to poor management initiatives and small development decisions by various governmental bodies in order to accommodate the 92% of Ontarians (35% of Canadians) and support industrial, residential, and agricultural needs over time [55]. This creates a necessity for more careful planning from large-scale perspectives in order to understand landscape composition and configuration, prevent habitat fragmentation, and restore larger wetland patches [44]. This phenomenon, in combination with the OMNRF goals to protect and restore Ontario's wetlands, make for an ideal study area for creating a multi-criteria wetland suitability index (WSI).  The Mixedwood Plains is the most populous and industrious ecozone in Canada; it is considered the second most fragmented ecozone in Canada [52]. Thirteen major urban areas are found within this ecozone, including the first, second, and fourth most populated Canadian cities of Toronto, Montreal, and Ottawa, respectively. Between 1971 and 2006, the human population of the Mixedwood Plains increased from approximately 11 million people to more than 16 million people; the rapid population growth was projected to increase by 30% by 2031 [53]. As a result, the Mixedwood Plains has seen the largest loss in wetlands, which in turn impacts the largest populations of Ontarians. Despite suffering a significant loss in wetland area (roughly ~72%) since presettlement, Ontario remains home to 6% of the Earth's total wetland area. For comparison, the conterminous United States has experienced an estimated loss of 53% of its presettlement wetlands, contributing to the large-scale erasure of wetlands and restoration needs in North America [54]. Specific to the Mixedwood Plains Ecozone, the wetland land cover classification was estimated at 25% of the total

Data Preparation
Criteria chosen for the wetland suitability index (WSI) included biophysical characteristics already present that would best support a wetland, land cover representing areas of need for wetland ecosystem services, and constraint variables. Biophysical variables present in similar wetland indices were selected, including soil drainage, groundwater level, slope, and proximity to streams [18,20,24,31,33,45,46]. In previous studies, ecosystem services were considered by incorporating areas of need such as Sustainability 2020, 12, 9953 5 of 21 agricultural land and frequently flooded areas to make best use of a wetland's flood mitigation and water filtration services [25,27,46]. Last, major waterbodies and built-up areas were included as obstacles to restoration in the index. In an effort to best model present conditions of the Mixedwood Plains landscape, the most current, complete, and highest possible resolution data were selected (Table 1). A raster data type was chosen as the medium for calculating and illustrating the WSI. Raster data allow for pixel-by-pixel comparison within a landscape, allowing for finer resolution results when compared to vector boundaries. Each variable was resampled to 50 m, because it was the lowest resolution of the input data ( Table 2). This resolution was applied to all of the input data layers and provided the most detail without oversimplifying higher resolution data. For policy and planning purposes, using the finest resolution creates more spatially precise results that can later be aggregated and applied to different management scales [25]. All input layers were projected to the same georeferencing system to avoid any mismatch issues. In ESRI's [56] ArcMap 10.7, each variable was processed using Model Builder to ensure uniform parameters such as projection, snap raster, resolution, and pixel alignment. Within that geographic information systems (GIS) software, variables underwent differing geoprocessing to ensure they represented the intended goals (see Table 2).

Suitablity Criteria and Normalization
The WSI was created based on seven criteria, normalized from 0 (low suitability) to 10 (high suitability), and illustrated through a weighted composite raster ( Figure 2). A 0-10 normalization allowed for the rescaling of data of different ranges, some of which could not be divided into equal numbers of classes due to the qualitative nature of certain data sources. The 0-10 works well, as it is more likely to be universally recognized as a scale of quality. Additionally, it provides a larger and more gradual range of potential site suitability that makes best use of the fine raster resolution. Raster analysis allows for each grid cell to have an assigned value for each criterion. This simplistic method works for quantifying the cell's relative suitability for a wetland restoration/reconstruction project [44]. In the forthcoming section, weighting was determined by reviewing relevant literature on the applicable biophysical properties for wetland reconstruction [18,25,[44][45][46]. The seven criteria variables-agriculture, built-up, soil drainage, groundwater, major water bodies, slope, and stream proximity-are further justified and explained in the next paragraph.  (i) Agriculture: Due to the wetland loss to agricultural land use in Ontario specifically [31,32], it was included as a variable in the model with a ranking of 10 for agriculture and 0 for all other land use. Agricultural land offers a blank slate free of infrastructure or other barriers for wetland restoration, increasing the likelihood that it will serve as a suitable site [27,46]. (ii) Built-up: Roads, impervious surfaces, and infrastructure, which made up the built-up variable, were ranked into two classes, representing built-up, with a score of 0, or not, with a score of10. If the landscape is already built with impervious associated surfaces it has low probability of being reconstructed as wetlands. (i) Agriculture: Due to the wetland loss to agricultural land use in Ontario specifically [31,32], it was included as a variable in the model with a ranking of 10 for agriculture and 0 for all other land use. Agricultural land offers a blank slate free of infrastructure or other barriers for wetland Sustainability 2020, 12, 9953 7 of 21 restoration, increasing the likelihood that it will serve as a suitable site [27,46]. (ii) Built-up: Roads, impervious surfaces, and infrastructure, which made up the built-up variable, were ranked into two classes, representing built-up, with a score of 0, or not, with a score of10. If the landscape is already built with impervious associated surfaces it has low probability of being reconstructed as wetlands. (iii) Soil drainage: Many studies cited soil drainage to be prioritized in suitability models, as wetlands often consist of poorly drained soils. Soil type was ranked by drainage type, from rapidly drained (worst) with a score of 2.5 to very poorly drained (best) with a score of 10. In addition, frequently inundated soils were included among the "very poorly drained" soil class to prioritize flood mitigation to sites that need it most. (iv) Groundwater: Sivakumar and Ghosh [18] prioritized shallower groundwater in their overlay analysis wetland mitigation plan. For the purpose of this study, a groundwater surface was interpolated using the Provincial Groundwater Monitoring Network (PGMN), consisting of 367 points [57]. Kriging was used to interpolate a continuous raster surface from the 367 PGMN groundwater depth values using ESRI's [56] ArcGIS 10.7. In addition to the interpolated surface, kriging produces an output variance layer, which unlike other interpolation methods, provides the user with the confidence values to statistically validate the predicted pixel values (Appendix A  Table A1). Kriging spatial interpolation models with a mean prediction error (MPE) greater than and less than 0 represent under-and overestimation of kriged values in comparison to the actual measured values [58,59]. Underestimation of variability in the kriged surfaces is identified when the standardized root-mean square prediction error (SRMSPE) is >1, while a SRMSPE <1 represents an overestimation of variability [59,60]. The resulting groundwater surface was then reclassified into three quantile classes, ranking shallow with a score of 10 and a deep water table with a score of 0. (v) Slope: Wetlands require flat or gently sloping topography in order for water to collect [33,[44][45][46]. Slope ranged from 0-41.54 degrees with a mean of only 1.12 degrees. Quantile classes were used to divide slope into three rank classes, with a score of 10 for a low/flat degree of slope to a high/steep degree of slope with a score of 0. (vi) Stream proximity: Darwiche-Criado [45] and Huang et al. [24] used 500 m as the threshold for proximity to streams in their wetland site selection processes. As freshwater streams are abundant in the Mixedwood Plains Ecozone, and most pixels were within 500 m of a stream network, 250 m and 500 m stream buffers were applied to create two classes of closeness, scoring 10 and 5, respectively. (vii) Major waterbodies: The waterbodies variable required only two classes (acceptable and not acceptable). This variable acted as a constraint in the model (locations where reconstructed wetlands could not be placed), unlike the six other variables of varying suitability [44,61]. The waterbodies variable was assigned a value of "NoData" for major waterbodies, while the rest of the pixels were assigned a 0. This ensures that permanent waterbodies are not included in the final WSI surface.

Literature Analysis for Index Weighting
There remains no agreed method for weighting indicators (e.g., wetland suitability criteria) or sub-indices when creating a multi-metric composite index (e.g., wetland suitability index; WSI). Therefore, a heuristic method for weighting the six non-constraint wetland suitability criteria was created through a scoping review of relevant literature. Scoping reviews of relevant literature lend way to informed decision-making based on collective knowledge [62]. This review will be used to determine the relative importance, and therefore the weighting, of each variable in the index. By considering the relative importance of each suitability variable through popularity in studies, many professional opinions were considered at once, effectively reducing bias and crowdsourcing consensus. To begin the scoping review, a query was designed to capture articles in ProQuest's Science and Technology database relevant to this study to ensure a wide search without duplicates [63]. The query ("Wetland reconstruction" OR "wetland restoration" OR "wetland suitability" OR "wetland creation" OR "wetland site suitability") AND ("multi-criteria" OR "index" OR "indicator" OR "GIS") was then further narrowed down by including only peer-reviewed articles, in English, from the past 10 years. The search included the entire contents of each item and was not limited geographically (as this resulted in too few items). Next, content analysis on the selected articles was carried out to find the Sustainability 2020, 12, 9953 8 of 21 number of relevant articles that considered each variable [64]. To ensure consistency between searches, search terms were determined using one word that best described the variable.
The relative importance of each wetland suitability criterion was decided by the number of articles where it appeared in the relevant literature. Using the Analytical Hierarchy Process (AHP), the relevant importance of each variable can be used to determine each variable's weight in the model and therefore its influence on the result [61,65]. Similar studies on creating raster surfaces for environmental planning scenarios have used Saaty's [65] AHP to find the importance of criteria in relation to each other, and to the project's goal (in this case wetland reconstruction) [44,61,66]. As the variable major waterbodies acted as a constraint in the model, it only represented the inclusion (rank = 0) or exclusion (rank = NoData) and therefore was not weighted with the other six variables of varied rankings. Every remaining relationship was then classified on Saaty's [65] nine-point reciprocal scale (Table 3). Based on the following differences in their relative counts, relationships between criteria were classified into equal groupings: 1-equal (+/−0), 3-moderate (+/−75-150), 5-strong (+/−225-300), 7-very strong (+/−375-450), or 9-extreme (+/−>525) in their relative importance to the goal of wetland reconstruction (Table 3). Variables were then assigned the corresponding value (1,3,5,7,9) if they were of more importance in the relationship, or the reciprocal value (1, 1/3, 1/5, 1/7, 1/9) if they were of less importance. Intermediate values (2, 4, 6 . . . etc.) were assigned if differences in count fell between the above categories. The resulting pairwise comparison matrix represents each variable's relative importance to wetland suitability based on the literature content analysis (Table 3). Following the creation of the pairwise comparison matrix, each value must be normalized by the sum of its corresponding column, creating the normalized pairwise metrics (Appendix A Table A2). The criteria weights were then calculated by summing each row and dividing by the number of variables (see Table 3). For the resulting weights to be considered acceptable, the consistency within the weights was checked using the following standards [65]: The consistency ratio (CR) describes the level of inconsistency in the model. A consistency ratio of <0.1 is considered the acceptable threshold. For this process, each value in the pairwise comparison matrix is multiplied by its criteria weight. When solved, the resulting consistency matrix can be used to calculate the weighted sum value (or sum of each row; Appendix A Table A3). The weighted sum value is divided by the corresponding criteria weight, and the result is summed and averaged, which = λmax. To calculate the consistency index (CI), If n = 6 and λmax = 6.4783 The consistency ratio (CR) is CI divided by the random index (RI)-determined by n: Sustainability 2020, 12, 9953 9 of 21 Therefore, the CR is 0.0771 and the proportion of inconsistency is considered acceptable. The final ranking of each criterion represents its suitability to wetland landscapes while weighting represents each criterion's relative importance to wetland suitability based on prevalence in relevant literature.

Wetland Suitablity Index Validation
High suitability sites, also referenced as "best sites," from the WSI were determined by selecting those locations one standard deviation above the mean WSI score from the weighted composite raster. Comparison of WSI best sites with historic and existing wetland databases helps to understand WSI predictions of most suitable sites for wetland reconstruction. Utilizing the location of historic wetlands helps to focus efforts on wetland reconstruction by confirming the existence of suitable wetland criteria where present-day wetlands do not exist [67]. To account for the extensive wetland loss in Southern Ontario, best sites were compared with both the Ducks Unlimited Canada (DUC) historic wetland layer and the Land Information Ontario (LIO) existing wetlands layer. To remain consistent with the 50 m cell size of the rest of the study, wetlands under 1275 m 2 (or half of the 2500 m 2 cell size) were deleted from each vector file. This ensured that wetlands that were too small to be of at least half of one 50 m pixel were not transformed into much larger wetlands when rasterized. Next, binary rasters were created for both sets: (i) best sites and historic, and (ii) best sites and existing. Binary similarity coefficients can be used to show similarity in the distribution of features between two areas in ecological studies [68][69][70].
For the purpose of this study, the Jaccard index was used to measure the shared presence of wetlands in the WSI best sites raster with existing and historic wetlands. The Jaccard index [71] is one of the most widely used indices for quantifying similarity of binary data in ecology and biogeography; albeit, it does not consider negative matches [72]. In ESRI's [56] ArcMap 10.7, the raster calculator was used to find the intersection (& operator) and union (| operator) between WSI best sites and both existing and historic wetlands. The Jaccard index was calculated by dividing the intersection by the union in each case. The Jaccard's index shows complete dissimilarity where Jaccard = 0, and complete similarity where Jaccard = 1.

Wetland Suitablity "Hot Spot" Analysis
The first law of geography asserts that near things in space are more similar (spatially autocorrelated) than farther apart things [73]. To help answer this study's guiding research questions, a two-step spatial analysis was conducted to test the global level of spatial autocorrelation and illustrate statistically significant spatial clusters of wetland suitability using Getis-Ord General G and Getis-Ord Gi*, respectively. Getis-Ord General G is considered a global test of spatial autocorrelation; Getis-Ord Gi* is recognized as a local index of spatial association (LISA) test and is often used as a "hot spot" analysis [74]. Both global and local Getis-Ord tests compute z-scores and p-values, specifying whether features are non-randomly distributed (global) and where those feature's patterns are located (local) [75]. For statistically significant positive z-scores, the larger the z-score, the greater the clustering of high values (hot spot); for negative z-scores, the smaller the z-score, the greater the clustering of low values (cold spot; [56]).
Census subdivisions [76] and quaternary watershed boundaries [77] were chosen as aggregation boundaries for the spatial univariate tests, as they had a similar number of features for results comparison (261 and 296, respectively) without oversimplifying changes in wetland suitability over the study area. From an environmental management perspective, the use of both political and natural boundaries provides municipal and regional planning decision-makers context, which remained consistent with the goals of this study. The Zonal Statistics to Table tool was used to aggregate and calculate the mean WSI value for each polygon, which could then be attribute joined to its corresponding polygon. For both the global General G and local Gi* tests, distance thresholds were set to 32.2 km and 33.0 km for census subdivisions and watersheds, respectively; the Calculate Distance Band, within the Neighbor Count tool, revealed the kernel search distance when eight neighbors were reached for all features. The local Gi* statistic allowed each feature's z-score to be illustrated, and was used to display geographic "hot spots" and "cold spots" of mean WSI across the study area census subdivisions and watersheds. The aforementioned spatial analysis procedures and Getis-Ord non-stationarity tests were conducted using ArcMap 10.7 Spatial Statistics tools [56].

Index-Ranking and Weighting
The literature review drew 1253 results of relevant documents using the search query ("Wetland reconstruction" OR "wetland restoration" OR "wetland suitability" OR "wetland creation" OR "wetland site suitability") AND ("multi-criteria" OR "index" OR "indicator" OR "GIS"). These were limited to English, peer-reviewed full texts from the last 10 years. All searches were within the text (not limited to the title only) and duplicates were removed. The content analysis for each variable drew in the following counts: (agriculture) = 754, (urban) = 653, (soil) = 964, (groundwater) = 447, (slope) = 358, and (stream*) = 608. After undergoing the Analytical Hierarchy Process (AHP), the three criteria with the highest influence (weighting) on the index were soil drainage (48.2%), agriculture (21.8%), and built-up (12.3%). The remaining three criteria accounted for the outstanding 18% of the model (Table 4). The wetland suitability index (WSI) was created based on the aforementioned justified seven criteria, normalized from 0 (low suitability) to 10 (high suitability), and illustrated through a weighted composite raster. Therefore, using the symbols from Table 4, the WSI equation is as follows:

Wetland Restoration Suitablity
The WSI was created using seven variables from a variety of sources that represent geographical needs for wetland reconstruction, reflect natural wetland form, and make best use of ecosystem services ( Figure 3). These data were combined to create this index for which suitable sites for wetland construction could be determined. After summing the aforementioned criteria, the resulting WSI raster surface had a suitability range of (worst) 0.000-9.999 (best) with a mean suitability of 6.40 and a standard deviation of 1.96. Raster cells with WSI index values one standard deviation above the mean (>8.35) were chosen as "best sites." These made up 14.529% of pixels in the final wetland suitability surface. If each cell is 0.25 ha at a 50 m resolution, then best sites have a combined area of 993,299.25 ha. The resulting surface accentuated high wetland suitability in the southernmost region of Ontario, the Niagara Peninsula, and the northeast part of the study area south of Ottawa. Areas of very low suitability can largely be attributed to the presence of urban areas, and not applicable/rapidly drained soils.

Validation with Historic and Existing Wetlands
Most suitable sites (best sites) established from the WSI were compared and contrasted with historic and existing wetlands (Figure 4). The isthmus found between Lake St. Clair and Lake Erie in the southernmost region of the study area has undergone significant wetland loss since presettlement; a majority of suitable sites indicated by WSI intersected with Ducks Unlimited Canada's historic wetlands here ( Figure 4A). The northeast region of the study area, especially lands south of Ottawa, also had a high degree of overlap between WSI best sites and the historic wetlands data set. That said, greater corroboration between WSI and once-present wetlands was found throughout the Mixedwood Plains Ecozone. When contrasting WSI with existing wetlands, the intersect analysis rendered sparse but consistent results throughout the study area ( Figure 4B). A greater number of WSI best sites intersected with existing wetlands in the central portion of the study area. When contrasting WSI best sites with historic and existing wetland data sets, the Jaccard coefficient was used to evaluate spatial data overlap. The Jaccard coefficient scores between 1 and 0, with 1 being completely similar or overlapped and 0 having no similarity or overlap. Validation results between WSI best sites and historic wetlands resulted in Jaccard = 0.231, indicating some similarity. Validation results between WSI best sites and present existing wetlands resulted in Jaccard

Validation with Historic and Existing Wetlands
Most suitable sites (best sites) established from the WSI were compared and contrasted with historic and existing wetlands (Figure 4). The isthmus found between Lake St. Clair and Lake Erie in the southernmost region of the study area has undergone significant wetland loss since presettlement; a majority of suitable sites indicated by WSI intersected with Ducks Unlimited Canada's historic wetlands here ( Figure 4A). The northeast region of the study area, especially lands south of Ottawa, also had a high degree of overlap between WSI best sites and the historic wetlands data set. That said, greater corroboration between WSI and once-present wetlands was found throughout the Mixedwood Plains Ecozone. When contrasting WSI with existing wetlands, the intersect analysis rendered sparse but consistent results throughout the study area ( Figure 4B). A greater number of WSI best sites intersected with existing wetlands in the central portion of the study area. When contrasting WSI best sites with historic and existing wetland data sets, the Jaccard coefficient was used to evaluate spatial data overlap. The Jaccard coefficient scores between 1 and 0, with 1 being completely similar or overlapped and 0 having no similarity or overlap. Validation results between WSI best sites and historic wetlands resulted in Jaccard = 0.231, indicating some similarity. Validation results between WSI best sites and present existing wetlands resulted in Jaccard = 0.028, indicating greater dissimilarity.
While the similarity between existing wetlands and best sites is lower, it is important to note that many of these landscapes have undergone land-cover change and their presettlement wetlands have been repurposed into other land uses. In sum, less than 30% of presettlement wetlands remain present-day for existing overlap validation, which is particularly evident within the isthmus of the southernmost region of the study area.
Sustainability 2020, 12, x FOR PEER REVIEW 13 of 22 = 0.028, indicating greater dissimilarity. While the similarity between existing wetlands and best sites is lower, it is important to note that many of these landscapes have undergone land-cover change and their presettlement wetlands have been repurposed into other land uses. In sum, less than 30% of presettlement wetlands remain present-day for existing overlap validation, which is particularly evident within the isthmus of the southernmost region of the study area.

Hot Spot Analysis
The global Getis-Ord General G analysis revealed that the patterns of mean WSI, aggregated to both the census subdivisions (n = 261) and quaternary watersheds (n = 296), were not significantly different than a random spatial distribution. Regarding the governmental subdivisions, the observed Getis-Ord General G score was 0.031 and the z-score was 1.137. Concerning the watersheds, the observed Getis-Ord General G score was 0.029 and the z-score was 0.133. That said, the local Getis Ord Gi* statistic identified spatial clusters of low and high wetland suitability for both the census subdivision and quaternary watershed boundaries ( Figure 5).

Hot Spot Analysis
The global Getis-Ord General G analysis revealed that the patterns of mean WSI, aggregated to both the census subdivisions (n = 261) and quaternary watersheds (n = 296), were not significantly different than a random spatial distribution. Regarding the governmental subdivisions, the observed Getis-Ord General G score was 0.031 and the z-score was 1.137. Concerning the watersheds, the observed Getis-Ord General G score was 0.029 and the z-score was 0.133. That said, the local Getis Ord Gi* statistic identified spatial clusters of low and high wetland suitability for both the census subdivision and quaternary watershed boundaries ( Figure 5).
Of the 261 census subdivision features, 60 were of significant spatial clustering with less than 1% chance of being random, and 35 were spatially clustered to some degree of confidence ( Figure 5A). Four of five 99% confidence hot spots were present in the southernmost region's Middlesex and Lambton census divisions, signifying high wetland suitability in those political management boundaries on the west side of the study area. Hot spots reporting lower degrees of confidence were found surrounding the municipalities of Ottawa, Stratford, and Waterloo, signifying marginal wetland suitability around those cities. Cold spots made up 25 of the significant census subdivisions. All four clustered subdivision cold spots of high (99%) confidence were within the Toronto census metropolitan area (CMA), surrounded by five clustered subdivisions at 95 and 90% confidence level. Multiple negative clustered subdivisions were also found surrounding the City of Cambridge and within the Peterborough census division.
Of 296 watershed boundaries, 125 were spatially clustered to a degree of confidence, including 18 watersheds at 99% confidence, 30 at 95% confidence, and 16 at 90% ( Figure 5B). Notably, four of the 18 hot spot watersheds at the 99% confidence level encompassed the Ausable River watershed, signifying high wetland suitability in those environmental management boundaries on the west side of the study area. Other high-confidence hot spots surrounded the Sydenham River, the Thames River outlet, and the Teeswater River outlet. Cold spots accounted for 61 of the 125 significantly clustered watersheds. Of the 33 negatively clustered watersheds, 13 were grouped at 99% confidence along the northwest shore of Lake Ontario around Toronto, including the Credit River, Don River, and Humber River watersheds. Other negative or "cold spot" clustered watersheds (with less than 1% chance of being random) were found in the Georgian Bay peninsula in the northwest of the study area. Notable cold spot watersheds, at the 95% confidence level, were also found in the Kawartha Lakes, Gull River, and Crowe Lake watersheds in the north-central part of the study area.
intersection of a Land Information Ontario (LIO) existing wetlands and best sites designated by the WSI.

Hot Spot Analysis
The global Getis-Ord General G analysis revealed that the patterns of mean WSI, aggregated to both the census subdivisions (n = 261) and quaternary watersheds (n = 296), were not significantly different than a random spatial distribution. Regarding the governmental subdivisions, the observed Getis-Ord General G score was 0.031 and the z-score was 1.137. Concerning the watersheds, the observed Getis-Ord General G score was 0.029 and the z-score was 0.133. That said, the local Getis Ord Gi* statistic identified spatial clusters of low and high wetland suitability for both the census subdivision and quaternary watershed boundaries ( Figure 5).

Wetland Suitablity Criteria and Their Importance
The first research question asked: Which physical and ecological landscape criteria represent high suitability for wetland reconstruction? Cited articles for similar studies often named shallow groundwater, proximity to streams, soil drainage, and low slope as ideal properties for wetland reconstruction. Criteria were also considered based on previous wetland studies in the Mixedwood Plains. The Detailed Soil Survey (DSS) from the Canada Land Inventory (CLI) contained updated data from previous CLI soil surveys. These data are an updated version of that used by Snell [31] to create the first wetland inventory in Ontario. Additionally, Duck's Unlimited Canada [32] utilized the CLI soil survey (combined with quaternary geology) to identify poorly and very poorly drained soils as indicators of historic wetlands. By using a current soil survey and ranking poorly and very poorly drained soils highest, consistency was maintained with these previously conducted studies. Last, to incorporate ecosystem services, two factors were considered areas of need: agriculture, and frequently inundated soils. In Ontario, where many wetlands have been eliminated due to intensive agriculture, the landscape is not only suitable for reintroduction of wetlands, but needs them [27]. Agriculture characterized by flood risk, runoff, soil erosion, and low biodiversity was included in the criteria in contrast to a wetland's ability to mitigate floods, introduce and maintain biodiversity, slow erosion, and store and filter water [78,79]. Frequently inundated soils were included among the highest-ranking soil class to prioritize flood mitigation to sites that need it most. The ranking of each criterion represents its suitability to a wetland project with a range from 0 (worst) to 10 (best).
Existing wetlands were left out of the wetland suitability index (WSI) criteria for several reasons. Few WSI best sites intersected with existing wetlands, as revealed by the Jaccard index; greater similarity was found between best sites and historic wetlands, corroborating the severity of historic wetland cannibalization. The mass eradication of approximately 72% of Ontario's wetlands to dredging, filling, and alterations to hydrology meant that existing wetlands could not be considered representative of Ontario's suitable land for wetland sites. Similarly, the remaining existing wetlands are not necessarily emblematic of the most ideal quality in Ontario, but more likely land that was left in its original condition, as it was not desirable for another land use. In sum, if the WSI were to incorporate the minimally existing wetlands into its criteria then it would be biased away from landscapes more suitable for reconstruction. These important findings were evident from the Jaccard spatial overlay analysis found within this study. Therefore, due to the extreme eradication of wetlands globally and their sparse spatial patterning across study areas, caution should be used when incorporating existing wetland inventories into future suitability analyses.
The second research question asked: Of common wetland suitability metrics, which are most important? In this paper, a novel technique for advancing wetland suitability and prioritization was created by analyzing relevant connected literature via scoping review. This served as an objective approach for understanding each criterion's relative importance to wetland suitability by considering the collective contents of many research articles at once. Using Saaty's [65] Analytical Hierarchy Process (AHP) for determining multi-criteria weights, the relative importance of each variable from the literature review was transformed into criteria weights for the WSI. The criteria "soil drainage" and "agriculture" drew the most results, with weightings of 48.2% and 21.8%, respectively. The criteria search term "slope" had the fewest results and was the least important criterion, weighted at only 2.8%. In retrospect, it was considered that this literature review method for creating sub-index weights may be capturing research popularity rather than accepted scientific need. A fertile area for future research exists with using literature reviews most effectively for index weighting. It is hoped through the example of this research that other environmental management index initiatives will consider adopting this mixed-method technique to suit their needs (i.e., restoration of the Prairie Pothole Region).

Wetland Suitability Patterns and Characteristics
The third research question asked: Can a multi-criteria WSI effectively locate high and low wetland suitability across the Ontario Mixedwood Plains Ecozone? Getis-Ord Gi* charted statistically significant "hot spots" and "cold spots" of wetland suitability. Best sites were clustered around the City of Ottawa, within the isthmus found between Lake St. Clair and Lake Erie, and on the Niagara Peninsula. These regions are characterized homogenously with low slope, high amounts of agriculture land, and poor/very poor soil drainage. Medium suitability was scattered throughout central and northern Ontario and reflects the spatial patterns of agriculture, shallow groundwater depth, higher slopes, and well-drained soils. Roads and urban areas had the lowest values. Getis Ord Gi* analysis indentified statistically signifcant local clusters of high and low suitabilty. High confidence (99%) cold spots occur in Toronto due to its dense urban landscape. Statistically significant hot spots made up 13.4% of census subdivisions and 21.62% of quaternary watersheds. Hot spots of varying (90-99%) confidence were present for both census subdivision and quaternary watershed bondaries surrounding the City of Ottawa and in the southernmost portion of the Mixedwood Plains Ecozone moving north. These are largely agricultural areas with phyiscal characteristics that could support wetlands such as poorly drained soils and shallow groundwater.
Last, the fourth research question asked: How do best sites from the WSI compare and contrast to both inventories of presettlement wetlands and current existing wetlands? Ducks Unlimited Canada presettlement wetlands and Land Information Ontario existing wetlands were used in an overlay analysis for validating WSI best sites. Due to the high ranking and weighting (10 at 21.8%) of agriculture lands in the index, best sites were largely located on agricultural landscapes compared to any other land cover. This resulted in the index identifying few existing wetlands as they fell into the "Not" agriculture class ranked at 0. This, in contrast with the much higher similarity between WSI best sites and historic wetlands, testifies to the mass erasure of wetlands in Ontario specifically due to agriculture. This research has demonstrated that many of the areas in which wetlands have been removed from the landscape may still be suitable for sustaining a wetland habitat. Areas with suitable sites and few existing wetlands should be prioritized in wetland restoration and reconstruction planning. Suitable sites bordering existing wetlands could represent opportunities to increase patch sizes and bolster connectivity between existing habitats. Areas where best sites intersect with presettlement wetlands can be viewed as opportunities to restore historic wetlands, as these sites still consider present-day obstacles such as infrastructure and were formerly successful wetlands.

Recommendations and Future Considerations
Canada's Federal Policy on Wetland Conservation of 1991 recognizes the socioeconomic value of wetlands in the form of tourism, sustaining water resources, hunting, fishing, and forestry [80]. The federal policy aims to protect these functions by managing national parks, promoting research, and encouraging the conservation of wetlands in decision-making processes [81]. In 2014, the National Wetland Conservation Fund (NWCF) was created as part of the Government of Canada's National Conservation Plan (NCP); however, the $50 million (CAD) allocation was redistributed to other causes in 2019 when the plan was cut short [82]. That said, the United States Fish and Wildlife Service, through the North American Wetlands Conservation Act (P.L. 101-233) of 1989 [83], continues to generate funding for the protection of wetland habitats for migratory birds in the United States, Mexico, and Canada. At the provincial level in Ontario, a convoluted system remains for administering and implementing policies and legislation governing wetlands. Since the establishment of the first provincial wetland policy in 1984, the management of these keystone landscape features is now done through a variety of policies that include over 20 unique legislation items governed and executed by five provincial ministries, two federal departments, a provincial agency (Niagara Escarpment Commission), 36 conservation authorities, and 444 municipalities [19]. Despite a strong promise by governmental and non-governmental organizations, it remains financially and logistically impossible to restore all wetlands to their presettlement landscape conditions. Therefore, effective wetland restoration, management, and planning requires uncomplicated near-, mid-, and long-term policies that simplify management processes and prioritize limited financial resources.
Despite the noteworthy research contributions from this study, there are limitations that inform future research directions. The limitations largely involve methodological issues with data from various sources and of various quality. While the Mixedwood Plains Ecozone includes Manitoulin Island and surrounding smaller islands, these had to be removed from the study due to a lack of data. However, the island was not included in other studies such as Ducks Unlimited Canada's Southern Ontario Wetland Analysis and thus had no loss of continuity in comparing the final suitability surface with preceding studies. Additionally, some data sources were over a decade old, of varying resolutions, or were compilations of old and new data sets. These temporal and resolution incongruities may have compromised the results, but best practices were always practiced. These data include the Agricultural Resource Inventory (ARI), as well as the Detailed Soil Survey (DSS) from the Canada Land Inventory (CLI). To account for lower resolution datasets, finer resolution data had to be simplified to the determined resolution of 50m. This resulted in some data loss, yet increased confidence in WSI outcomes. Last, the creation of a groundwater raster surface involved kriging interpolation, which can provide only a best estimate of groundwater depth. With the growing quantity of higher resolution geospatial data, future research can be improved by including more criteria or using more current data at a finer scale. For instance, these could include landscape ecology metrics that capture additional ecosystem goods and services such as habitat physical connectedness, biodiversity, and recreation [84,85]. Progressing this research to the next level using the aforementioned ideas will help to meet the Ontario Ministry of Natural Resources and Forestry goals, which call for a reinstatement of wetland area and function to landscapes where loss has been the greatest.

Conclusions
Wetlands provide diverse habitats and ecosystem goods and services, which are increasingly valuable under the current conditions of global change. The overwhelming loss of wetland area in Southern Ontario calls for macroscale conservation planning initiatives. In response, this paper was created to deliver landscape planners, ecological restoration scientists, regional planners, and environmental managers and decision-makers an applied example for systematically evaluating wetland suitability across space in temperate climate zones. The use of geographic information systems (GIS) for decision-making has allowed for the creation of replicable processes, which can be altered to suit the needs of this project and its stakeholders. In addition, considering multiple wetland suitability criteria allowed for cost-effective, environmentally conscious planning. Through the incorporation of ecosystem services into these criteria, sites that represent the best possible reconstruction and restoration outcomes were identified. In other regions around the world, wetland indices have been implemented to identify suitable wetland reconstruction sites, however, such an index did not previously exist for Ontario. Therefore, this study provides the first macroscale study of wetland suitability across the Ontario Mixedwood Plains Ecozone by investigating and answering its four research questions.
Restoring and maintaining wetlands are critical for building resilience to global change, and knowledge, information, and new methods are required to prioritize environmental management efforts. Therefore, to aid major Canadian agencies in their macroscale wetland restoration efforts in Southern Ontario, this study had the applied goal of determining where the most suitable wetland restoration sites were spatially distributed in the Mixedwood Plains and how these best sites compare to historic and present models of wetlands in the area. The generation of the Wetland Suitability Index (WSI) was based on seven criteria, normalized from 0 (low suitability) to 10 (high suitability), and illustrated through a weighted composite raster. Similarities between historic wetland extent and the best sites generated by the WSI illustrated Ontario's potential to recover from historic wetland loss. Specifically, the overlay analysis revealed greater similarity between high suitability sites and presettlement wetlands, supporting the magnitude and severity of wetland cannibalization. The WSI framework provides a useful tool for understanding the spatial distribution of wetland suitability on a large scale for wetland restoration in Ontario's Mixedwood Plains. More specifically, a novel mixed-methods modeling approach to regional wetland restoration provides a prioritization tool for enhancing macroscale ecological connectivity, services, and resilience. Further, returning wetlands to the Mixedwood Plains Ecozone landscape of Canada would help to combat the negative impacts of global change on biogeochemical systems by reintroducing natural water filtration, stabilizing temperatures, increasing biodiversity, improving ecological services, decreasing runoff, and reducing flood risk. This research has provided a roadmap for informing policies and planning that meet this goal. In sum, studies like this one help to mend humanity's bond with its life-supporting biogeochemical systems. Funding: This study was funded by Ryerson University. The University had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.