Inventory and Connectivity Assessment of Wetlands in Northern Landscapes with a Depression-Based DEM Method

Wetlands, including peatlands, supply crucial ecosystem services such as water purification, carbon sequestration and regulation of hydrological and biogeochemical cycles. Peatlands are especially important as carbon sinks and stores because of the incomplete decomposition of vegetation within the peat. Good knowledge of individual wetlands exists locally, but information on how different wetland systems interact with their surroundings is lacking. In this study, the ability to use a depression-based digital elevation model (DEM) method to inventory wetlands in northern landscapes and assess their hydrological connectivity was investigated. The method consisted of three steps: (1) identification and mapping of wetlands, (2) identification of threshold values of minimum wetland size and depth, and (3) delineation of a defined coherent area of multiple wetlands with hydrological connectivity, called wetlandscape. The results showed that 64% of identified wetlands corresponded with an existing wetland map in the study area, but only 10% of the wetlands in the existing map were identified, with the F1 score being 17%. Therefore, the methodology cannot independently map wetlands and future research should be conducted in which additional data sources and mapping techniques are integrated. However, wetland connectivity could be mapped with the depression-based DEM methodology by utilising information on upstream and downstream wetland depressions, catchment boundaries and drainage flow paths. Knowledge about wetland connectivity is crucial for understanding how physical, biological and chemical materials are transported and distributed in the landscape, and thus also for resilience, management and protection of wetlandscapes.


Introduction
Wetlands play an important role in the world's freshwater system [1]. As a result of changes in land and water use driven by humans and climate change, wetlands have declined globally [2,3]. According to Ramsar [4], 64% of the world's wetlands disappeared in the 20th century, with 40% of them occurring in the last 50 years. Wetlands are crucial in supplying ecosystem services such as the regulation of However, knowledge of how different wetland systems interact with their surroundings is lacking. Although previous studies have used technology to assess preferential water pathways in wetlands, such as dye tracer [39], computed tomography [40] and nuclear magnetic resonance [41], these methodologies are used at a small scale to understand the relationship between hydrological processes and soil properties. Few studies have focused on flow connectivity at a larger scale, based on field instrumentation, but still covering a relatively limited extent of the landscape [42]. Model frameworks have been developed recently to explore landscape permeability or dispersal characteristics within wetlands, but focusing on landscape connectivity for wildlife [43]. Existing models do not have the capacity to consider geometry or quantify subsurface connectivity in a wetlandscape perspective [32]. Comprehensive inventories of wetland areas would be useful in assessing the effectiveness of different regulatory and management activities, knowledge that could guide policymakers and stakeholders [22].
The aim of this study was to investigate the capability of a depression-based DEM method for mapping northern landscape wetlands and, in particular, to increase understanding of wetland dynamics through mapping hydrological connectivity among wetlands. Specific objectives were to: (1) identify wetland areas using a mapping technique based on DEM; (2) evaluate threshold values of size and depth to identify wetlands; and (3) assess the hydrological connectivity of wetlands between catchments. The method developed was based on the premise that surface depressions often give rise to temporary, seasonal or permanent wetlands, since runoff from upslope areas is retained at downslope sites within the landscape [44][45][46].

Study Site
The study site was a 214 km 2 area near the village of Kaamanen in the municipality of Inari in northern Finland (69 • 8.435 N, 27 • 16.189 E) ( Figure 1). The region is located in the northern boreal vegetation zone and in the subarctic climate zone (Dfc zone in the Köppen climate classification system), which is characterised by a cold climate without dry seasons and cold summers [47]. The study site is adjacent to Lake Inari, the third largest lake in Finland ( Figure 1) and just south of the Sammuttijänkä-Vaijoenjänkä mires, a Ramsar site. Based on data from the closest nearby weather station (Inari Ivalo airport; 68 • 36.65 N, 27 • 24.833 E) for the period 1981-2010, mean annual temperature at the site is −0.4 • C and mean annual precipitation is 472 mm [48]. There is no permafrost in the area, and the topography is flat, characterised by flat peatlands and some forested hills and esker ridges, with elevation varying between 138 and 242 m a.s.l. Managed pine (Pinus sylvestris) forests dominate the area, but birch (Betula pubescens) and mixed pine-birch forests are also present [49]. Lakes and ponds cover 8% of the study area, while approximately 36% is covered by peatlands, including both tree-covered and open peatlands.

Mapping Wetlands
The depression-based DEM method investigated for mapping wetlands and their hydrological connectivity included the following three steps: (1) identification and mapping of wetlands; (2) identification of threshold values; and (3) delineation of wetlandscapes using hydrological connectivity.

Data Acquisition and Pre-Processing
The data used consisted of a DEM, an aerial image and a topographic database managed by the National Land Survey (NLS) of Finland, which are freely available at https://tiedostopalvelu. maanmittauslaitos.fi/tp/kartta?lang=en. NLS produced the DEM from LIDAR data and manually inspected it using aerial image interpretation. The data had spatial resolution of 2 m × 2 m and vertical accuracy of 0.15-0.3 m, and were acquired on 12 July 2016. The NLS DEM was further pre-processed with a 3 × 3 cell moving window mean filter to remove artefacts which, unlike actual surface depressions, often have a more irregular shape and size and are shallower, and originate from construction operators when producing the DEM [50]. A four-band (blue, green, red, near infra-red) 0.5 m spatial resolution aerial image of the study area was acquired on 27 June 2016. The image was orthorectified by NLS and had positional accuracy of 0.5-2.0 m. The topographic database obtained included land cover information about lakes, streams and wetlands in a vector format and 1:10,000 scale. These terrain attributes were used as validation and comparison data specifically to identify threshold values, such as wetland minimum depth and size, to delineate wetlands.
Water 2020, 12, x FOR PEER REVIEW  4 of 22 are shallower, and originate from construction operators when producing the DEM [50]. A four-band (blue, green, red, near infra-red) 0.5 m spatial resolution aerial image of the study area was acquired on 27 June 2016. The image was orthorectified by NLS and had positional accuracy of 0.5-2.0 m. The topographic database obtained included land cover information about lakes, streams and wetlands in a vector format and 1:10,000 scale. These terrain attributes were used as validation and comparison data specifically to identify threshold values, such as wetland minimum depth and size, to delineate wetlands.

Identification and Mapping of Wetlands
The procedure of delineating potential wetlands using depression-based DEM was done using an addition to the existing ArcGIS toolbox called Wetland Hydrology Analyst, developed by Wu and Lane [51], which is freely accessible for download at https://figshare.com/articles/software/Wetland_Hydrology_Analyst_Toolbox/8866025. The different steps followed in the modelling work in the toolbox are illustrated as a workflow in Figure 2, and described below.

Identification and Mapping of Wetlands
The procedure of delineating potential wetlands using depression-based DEM was done using an addition to the existing ArcGIS toolbox called Wetland Hydrology Analyst, developed by Wu and Lane [51], which is freely accessible for download at https://figshare.com/articles/software/Wetland_ Hydrology_Analyst_Toolbox/8866025. The different steps followed in the modelling work in the toolbox are illustrated as a workflow in Figure 2, and described below.

Delineation of Wetlands
In the first modelling step, a depressionless DEM was generated by identifying and filling the depressions in the DEM. This task was performed using the priority-flood algorithm developed by Wang and Liu [52], where filling proceeds to the level where water would start to overflow from the perimeter. The geometrical boundary of the wetland was thereby defined by the "pond" created by flooding the depression. The filled depressions were then used as a template to subtract the depressions that led the defined threshold values, including minimum size and depth. This is a computationally more efficient procedure than searching for, identifying and selecting depressions from the original DEM [52]. In a second step, an elevation difference grid was obtained by subtracting the original DEM from the depressionless DEM, where each cell value represented the depression depth. In a third step, the difference grid was converted to a binary image where all cells with value greater than 0 (depressions) were set to 1 and all cells with 0 (non-depressions) were set to remain 0. In step four, the binary image was used as a template to subtract the original smoothed DEM data for the depression areas. In a fifth step, a region group algorithm was applied to the DEM subset, which is a method of connecting the continuous set of depression cells together, forming regions and naming them [53]. In step six, the zonal statistics of depressions were calculated from the input raster [54]. Once the grouping and statistics for each depression had been obtained, in step seven the potential wetlands were identified from the DEM subset using threshold values. Since depressions show many shapes and sizes in the landscape, defining the threshold values for minimum depression size and associated depth enabled depressions that are potential wetlands to be located. A separate vector file was obtained containing polygons with depressions that fulfilled the criteria, in other words potential wetlands.

Delineation of Wetlands
In the first modelling step, a depressionless DEM was generated by identifying and filling the depressions in the DEM. This task was performed using the priority-flood algorithm developed by

Delineation of Wetland Catchments
Once the wetlands were defined, their corresponding potential catchment was delineated. This was done by first calculating the flow direction based on the DEM, where the flow direction from a cell was defined by its steepest downslope neighbour [55], i.e., the D8 algorithm by O'Callaghan and Mark [56]. All depressions that did not fulfil the threshold values were assumed to be filled, to prevent water from being trapped in these incorrect wetlands. The flow direction was used to obtain the flow accumulation, determined by the number of cells that flowed into each cell downstream [57], which in turn was used to delineate flow paths in the following step. Cells with high values are likely to represent channels where water will accumulate, whereas cells with low values are higher up and represent landforms such as ridges [57]. The wetland catchment was obtained by the surrounding cells with a flow direction towards the wetland.

Delineation of Drainage Flow Paths
The potential flow path is the path that water takes when a wetland is filled and water starts to spill out from the wetland as runoff. To obtain the flow path, the depressionless DEM was used together with the vector file of identified wetlands and a filling parameter, which was rainfall intensity. The rainfall intensity parameter only affected the time taken for each wetland to be filled in the model, the so-called ponding time, and not the actual flow path. The ponding time indicated the order in which a wetland was filled in relation to adjacent wetlands, which was used to determine the flow path. The rainfall intensity was assumed to be equally distributed over the study site, which meant that the value of rainfall intensity did not have an impact on the delineation of wetlands. The least-cost-path search algorithm was used to obtain the flow path that water takes from a wetland when it spills out [58]. The vector file was used to locate the centroid of each depression, which was used as the starting point for the least-cost-path search algorithm. The distance from the wetland centroid to its spilling point was later removed. The ponding time was calculated from the ratio between potential water storage of the wetland (maximum ponding volume before water spills out, determined through Equation (1), and catchment wetland area and rainfall intensity, determined according to Equation (2).
where V [m 3 ] is potential water storage, n [-] is number of cells that fall within the depression, C [m] is elevation of the cell from which the water spills out, Z i [m] is the elevation of the grid cell and R [m] is the spatial resolution. where

Identification of Threshold Values
In order to identify potential wetlands using the model, the most representative thresholds of minimum wetland size and depth were (i) examined by first studying the output results from several model runs performed with varying threshold values randomly attributed, and (ii) compared with the NLS Finland map of wetlands in the study area. After initial parameter testing, in which the minimum wetland size varied from 50 m 2 to 10,000 m 2 , the model was run with two values (minimum wetland sizes of 100 and 1000 m 2 ) that showed correspondence with the previous mapping of wetlands by NLS Finland. Two minimum sizes were chosen so that comparison could be made between analysis that either included (minimum wetland size 100 m 2 ) or excluded (minimum wetland size 1000 m 2 ) small wetlands. Six values of minimum depth (0.1, 0.2, 0.3, 0.4, 0.5 and 0.6 m) were used to study the depth that showed the highest similarity with the existing mapping of wetlands by NLS Water 2020, 12, 3355 7 of 21 Finland, resulting in 12 outputs. The coverage of potential wetlands obtained from each model run was compared against the NLS Finland map of wetlands and lakes in the study area [59]. Since| identification of wetlands basically consists of finding depressions in the landscape, it is likely that wetlands will appear where lakes are located. However, the LIDAR technique used to obtain the DEM cannot penetrate water, which means that the water surface of lakes and ponds appeared as a solid surface when the model searched for depressions. Therefore, the wetlands that might coincide with lakes represented the volume that would be added if the lake were filled to the brim, before water poured out to a downstream depression. Based on the definition of wetland used in this study [25], the area that coincided with lakes was not included in wetland mapping.
A comparison was made using Jaccard similarity index [60], a measurement of similarity, of how well the wetlands identified by the model coincided with those in existing maps created by NLS Finland. Jaccard similarity index is obtained by dividing the size of the intersection between two sets by the combined size of the same sets (Equation (3)), as illustrated in Figure 3.
where A is one set (here NLS Finland topographic database) and B another set (here wetlands identified by the model).
In order to identify potential wetlands using the model, the most representative thresholds of minimum wetland size and depth were (i) examined by first studying the output results from several model runs performed with varying threshold values randomly attributed, and (ii) compared with the NLS Finland map of wetlands in the study area. After initial parameter testing, in which the minimum wetland size varied from 50 m 2 to 10,000 m 2 , the model was run with two values (minimum wetland sizes of 100 and 1000 m 2 ) that showed correspondence with the previous mapping of wetlands by NLS Finland. Two minimum sizes were chosen so that comparison could be made between analysis that either included (minimum wetland size 100 m 2 ) or excluded (minimum wetland size 1000 m 2 ) small wetlands. Six values of minimum depth (0.1, 0.2, 0.3, 0.4, 0.5 and 0.6 m) were used to study the depth that showed the highest similarity with the existing mapping of wetlands by NLS Finland, resulting in 12 outputs. The coverage of potential wetlands obtained from each model run was compared against the NLS Finland map of wetlands and lakes in the study area [59]. Since identification of wetlands basically consists of finding depressions in the landscape, it is likely that wetlands will appear where lakes are located. However, the LIDAR technique used to obtain the DEM cannot penetrate water, which means that the water surface of lakes and ponds appeared as a solid surface when the model searched for depressions. Therefore, the wetlands that might coincide with lakes represented the volume that would be added if the lake were filled to the brim, before water poured out to a downstream depression. Based on the definition of wetland used in this study [25], the area that coincided with lakes was not included in wetland mapping.
A comparison was made using Jaccard similarity index [60], a measurement of similarity, of how well the wetlands identified by the model coincided with those in existing maps created by NLS Finland. Jaccard similarity index is obtained by dividing the size of the intersection between two sets by the combined size of the same sets (Equation (3)), as illustrated in Figure 3.
where A is one set (here NLS Finland topographic database) and B another set (here wetlands identified by the model).
To further validate how accurately the modelled wetlands corresponded to wetlands mapped by NLS Finland, a confusion matrix was used. In addition, user accuracy (precision), producer accuracy (recall) and F1 score were calculated, using the following equations: To further validate how accurately the modelled wetlands corresponded to wetlands mapped by NLS Finland, a confusion matrix was used. In addition, user accuracy (precision), producer accuracy (recall) and F1 score were calculated, using the following equations: where TP is true positives, FP false positives and FN false negatives obtained from the confusion matrix.

Delineation of Wetlandscape Using Hydrological Connectivity
Surface depressions cause spatially varying and independently localised hydrological systems by trapping water, thereby breaking the connectivity in a landscape [61]. Since surface depressions can vary in size, they fill at different rates, affecting the hydrological connectivity by regulating the flow. Information on the filling-merging-spilling processes within a wetland is vital for characterising their hydrological connectivity [46]. When a wetland fills (step 1 in Figure 4a) with water and reaches its maximum storage, it can either merge (step 2 in Figure 4a) with an adjacent wetland or allow water to spill (step 3 in Figure 4a) to a downstream wetland (Figure 4b) [26].

Identified Potential Wetlands
The number of wetlands detected by the model varied depending on the threshold values. When using small minimum wetland size (100 m 2 ) and the shallowest minimum wetland depth (0.1 m), 27,000 wetlands were detected, covering 16% of the study area ( Figure 5). When the threshold values were more constrained for both minimum wetland size (1000 m 2 ) and depth (0.6 m), the number of wetlands detected decreased to around 700, covering around 4.5% of the study area ( Figure 5). The number of detected wetlands did not vary greatly when changing the minimum depth from 0.3 to 0.6 m. Above a minimum wetland depth of 0.3 m, the proximity between the graphs representing the proportion and number of identified wetlands started to increase. In the methodology developed, the hydrological connectivity obtained by the delineated drainage flow paths represents the maximum water capacity of the system (i.e., maximum change over time), since the model fills the wetlands to the brim to produce drainage flow paths. The flow paths from each identified wetland together make up an interconnected set of drainage routes, where upstream wetlands drain to downstream wetlands. In this way, information on adjacent wetlands that are connected, and how, is obtained.
To validate the delineated drainage flow paths, stream waters mapped by NLS Finland were used for comparison. Since information about the spatial extent of each wetland within the catchment and the specific drainage flow path was acquired after running the model, the hydrologically coupled system and the total hydrological catchment were obtained. For each wetland identified, information on the wetland located downstream was used to track how the water flowed in the wetlandscape. Some wetlands did not have a direct drainage flow path connection to a downstream wetland, but information on the wetland located downstream was still obtained. The wetlandscape was determined by selecting the connected wetlands (located just upstream and releasing water to downstream wetlands, or located just downstream and receiving water from upstream wetlands), with their individual catchments (Figure 4). This information was obtained in the resulting file after executing the tool. The wetlandscape obtained was thereby a type of catchment or sub-catchment where all wetland catchments that had interconnections resulted in a watershed with a common drainage outlet. The delineation of the wetlandscape was not automated and involved manual tracing of the connections in the data.

Identified Potential Wetlands
The number of wetlands detected by the model varied depending on the threshold values. When using small minimum wetland size (100 m 2 ) and the shallowest minimum wetland depth (0.1 m), 27,000 wetlands were detected, covering 16% of the study area ( Figure 5). When the threshold values were more constrained for both minimum wetland size (1000 m 2 ) and depth (0.6 m), the number of wetlands detected decreased to around 700, covering around 4.5% of the study area ( Figure 5). The number of detected wetlands did not vary greatly when changing the minimum depth from 0.3 to 0.6 m. Above a minimum wetland depth of 0.3 m, the proximity between the graphs representing the proportion and number of identified wetlands started to increase.

Similarity Analysis and Selection of Model Parameters
Wetlands identified by the model differed in spatial extent from the land cover map of wetlands from NLS Finland. Wetlands mapped by NLS Finland occupied 36% of the study area (Table 1) and coincided to 15% with the identified wetlands considering a minimum size of 100 m 2 and depth of 0.1 m, and to around 11% with wetlands considering a minimum size of 1000 m 2 and depth of 0.1 m ( Figure 6). The lakes mapped by NLS Finland covered 8% of the study area (Table 1), and coincided to a greater extent with the identified wetlands, with almost 31% similarity when using threshold values of 100 m 2 minimum size and 0.1 m depth, and 37% when using 1000 m 2 minimum size and 0.1 m depth ( Figure 6). When the minimum depth was ≤0.3 m, the similarity with lakes was higher for the minimum area of 1000 m 2 , while at minimum depths >0.3 m, the similarity with lakes was equally high for both minimum area thresholds. In general, similarities with both wetlands and lakes mapped by NLS Finland decreased when the minimum depth increased, except for minimum threshold values up to 0.2 m of minimum depth for wetlands and 0.3 m for lakes by NLS Finland.

Similarity Analysis and Selection of Model Parameters
Wetlands identified by the model differed in spatial extent from the land cover map of wetlands from NLS Finland. Wetlands mapped by NLS Finland occupied 36% of the study area (Table 1) and coincided to 15% with the identified wetlands considering a minimum size of 100 m 2 and depth of 0.1 m, and to around 11% with wetlands considering a minimum size of 1000 m 2 and depth of 0.1 m ( Figure 6). The lakes mapped by NLS Finland covered 8% of the study area (Table 1), and coincided to a greater extent with the identified wetlands, with almost 31% similarity when using threshold values of 100 m 2 minimum size and 0.1 m depth, and 37% when using 1000 m 2 minimum size and 0.1 m depth ( Figure 6). When the minimum depth was ≤0.3 m, the similarity with lakes was higher for the minimum area of 1000 m 2 , while at minimum depths >0.3 m, the similarity with lakes was equally high for both minimum area thresholds. In general, similarities with both wetlands and lakes mapped by NLS Finland decreased when the minimum depth increased, except for minimum threshold values up to 0.2 m of minimum depth for wetlands and 0.3 m for lakes by NLS Finland. Table 1. Total area of water surface land cover in terms of spatial coverage (km 2 ) and percentage of study area. The wetlands identified by the model were obtained using thresholds of 1000 m 2 minimum size and 0.2 m minimum depth.

Similarity Analysis and Selection of Model Parameters
Wetlands identified by the model differed in spatial extent from the land cover map of wetlands from NLS Finland. Wetlands mapped by NLS Finland occupied 36% of the study area (Table 1) and coincided to 15% with the identified wetlands considering a minimum size of 100 m 2 and depth of 0.1 m, and to around 11% with wetlands considering a minimum size of 1000 m 2 and depth of 0.1 m ( Figure 6). The lakes mapped by NLS Finland covered 8% of the study area (Table 1), and coincided to a greater extent with the identified wetlands, with almost 31% similarity when using threshold values of 100 m 2 minimum size and 0.1 m depth, and 37% when using 1000 m 2 minimum size and 0.1 m depth (Figure 6). When the minimum depth was ≤0.3 m, the similarity with lakes was higher for the minimum area of 1000 m 2 , while at minimum depths >0.3 m, the similarity with lakes was equally high for both minimum area thresholds. In general, similarities with both wetlands and lakes mapped by NLS Finland decreased when the minimum depth increased, except for minimum threshold values up to 0.2 m of minimum depth for wetlands and 0.3 m for lakes by NLS Finland. Considering the vertical accuracy of the DEM (0.15-0.3 m) and the shallowness of a wetland with 0.1 m depth, it is not reliable to use this value as a minimum depth to detect potential wetlands. However, the vertical accuracy refers to absolute values, which indicates that within a small area, the relative vertical accuracy can be higher. Minimum depths ≥ 0.4 m did not affect the number of identified wetlands or significantly increase the similarity with the wetlands or lakes mapped by NLS Finland. A minimum depth of 0.2 m and 0.3 m showed relatively high similarity with both the NLS Finland wetlands and the lakes, although with slightly lower similarity for lakes. The difference in similarity when using 1000 m 2 or 100 m 2 was small, but a minimum wetland size of 1000 m 2 showed higher similarity with the NLS Finland lakes and wetlands at a depth of 0.2 m. Based on these results, the thresholds of 0.2 m minimum depth and 1000 m 2 minimum size were selected for further analyses of wetlands.

Mapping of Identified Wetlands
Overall, 2122 wetlands were identified in the study area by the model, covering an area of 23 km 2 . Areas where wetlands, lakes and rivers overlapped (11 km 2 ) were not considered, resulting in total coverage of 12 km 2 of wetlands ( Table 2). Of these wetlands, 8 km 2 corresponded to NLS Finland wetlands (64% user accuracy), while 69 km 2 of the wetlands were missed (10% producer accuracy) (Table 2), with the F1 score for the wetland class being 17%. Figure 7 shows the cover of identified wetlands and their overlap with the wetlands and lakes mapped by NLS Finland.   Figure 8 presents selected parts of the study area to illustrate the delineated potential wetlands on a smaller scale, with more details. On one hand, the methodology developed identified some potential wetland areas not mapped by NLS Finland. In Figure 8a, the terrain around lakes/ponds identified as wetland by the model, but not classified as wetland by NLS Finland, was partly peatland and partly upland forest, according to visual interpretation of the aerial image. In addition, four lakes/ponds observed in the aerial image to be inundated area (Figure 8a) were not identified by the NLS Finland mapping (as a lake or wetland), but they were correctly identified by the model.  Figure 8 presents selected parts of the study area to illustrate the delineated potential wetlands on a smaller scale, with more details. On one hand, the methodology developed identified some potential wetland areas not mapped by NLS Finland. In Figure 8a, the terrain around lakes/ponds identified as wetland by the model, but not classified as wetland by NLS Finland, was partly peatland and partly upland forest, according to visual interpretation of the aerial image. In addition, four lakes/ponds observed in the aerial image to be inundated area (Figure 8a) were not identified by the NLS Finland mapping (as a lake or wetland), but they were correctly identified by the model. Figure 8b shows three small ponds that can give rise to wetlands, which were identified by the methodology, but left out of the NLS Finland map. Figure 8c shows an example of how the model sometimes covered larger areas than other mapping methods. In this example, the small patches of wetlands (NLS Finland) located at the bottom of the identified wetlands might be larger, according to visual interpretation of the aerial image. For example, the small pond located just above the wetland patches was not identified by NLS Finland, indicating that its mapping methodology missed some wetland areas. On the other hand, in some cases, the model failed to identify wetlands. Figure 8d,e show potential wetlands identified by the model that seem to be upland forest in visual interpretation of the aerial image. In addition, NLS Finland wetlands were not identified in the area. Figure 8f shows a situation where the identified wetlands coincided with wetlands by NLS Finland, but missed relatively large parts of the surrounding area, whereas the wetlands mapped by NLS Finland covered larger coherent areas of wetland. Figure 8g illustrates how the methodology sometimes had difficulty in mapping the wetland area, by identifying several patches spread within the study area. The wetlands by NLS Finland, which are not included in the images to enable the terrain below to be seen, provided coherent cover of the whole area. Figure 9 shows the results obtained from the connectivity analysis, illustrated for some selected areas (shown in Figure 9a). Figure 9b,c show the connectivity among wetlands in selected parts of the study area, where the modelled drainage flow paths from the identified wetland are delineated together with stream waters and wetlands mapped by NLS Finland. Figure 9d shows an example of a wetlandscape identified in the study area, illustrating how water in the landscape flows from upstream wetlands towards a larger identified wetland area, which surrounds a lake, through other identified wetlands and lakes located along its way. Water discharges at the outlet in the southernmost part of the lake. The yellow/orange area in Figure 9d is the spatial extent of the wetlandscape, which is made up of watersheds from each identified wetland forming a sub-catchment in the study area. The areas surrounding the delineated wetlandscape do not have a direct connection to the wetlands within the wetlandscape or to the outlet. In general, according to the flow path results, almost 60% of the identified wetlands in the study area were shown to receive water from no other wetland, comprising isolated wetlands. Around 23% received water from one other upstream wetland, while 16% of the identified wetlands received water from two to eight other upstream wetlands. The remaining 1% of identified wetlands received water from nine to 42 other upstream wetlands. The delineated flow paths extended over 410 km, of which 250 km were on NLS Finland wetlands, representing almost 61% of the delineated flow paths. Approximately 11% of the stream water mapped by NLS Finland was similar to the delineated flow paths.  Figure 9 shows the results obtained from the connectivity analysis, illustrated for some selected areas (shown in Figure 9a). Figure 9b,c show the connectivity among wetlands in selected parts of water from no other wetland, comprising isolated wetlands. Around 23% received water from one other upstream wetland, while 16% of the identified wetlands received water from two to eight other upstream wetlands. The remaining 1% of identified wetlands received water from nine to 42 other upstream wetlands. The delineated flow paths extended over 410 km, of which 250 km were on NLS Finland wetlands, representing almost 61% of the delineated flow paths. Approximately 11% of the stream water mapped by NLS Finland was similar to the delineated flow paths.

Model Performance in Identification of Wetlands
The methodology developed allowed the partial identification of wetlands and their interconnectivity. Of the wetlands identified using model threshold values of 1000 m 2 for minimum size and 0.2 m for minimum depth, 64% were also mapped as wetlands by NLS Finland. According to visual interpretation of the aerial image, the remaining 36% of identified wetlands comprised actual wetland areas not mapped by NLS Finland and areas incorrectly classified as wetlands by the model. However, only 10% of the wetlands by NLS were identified, and the F1 score was only 17%, indicating that the methodology did not cover the full extent of wetlands in the area. Nevertheless, earlier studies have shown that topographic maps such as the used NLS topographic database have uncertainties and have misclassified wetland areas [37]; therefore, the accuracy assessment results have uncertainties and should be treated critically.
Due to low producer accuracy in particular, the methodology cannot be used alone to map peatland and wetland areas in northern landscapes. One reason for the low coverage and low producer accuracy obtained was likely to be the assumptions made when modelling the complex hydrology and characteristics of wetlands. The methodology only uses DEM and two parameters, minimum wetland size and depth, while other factors that can indicate wetland area are not considered. The wetlands in the study area mostly consist of aapa mires [62], i.e., hydrologically continuous areas alternating with flarks (hollows) and strings (intervening ridges) [63] (Figure 8g). Aapa mires are not necessarily located in depressions and some can even be located on slopes [64], maintained by good access to water. These types of wetlands are difficult to map by assuming that wetlands are located in depressions and by using only two parameters (size and depth). Therefore, these assumptions might affect the model's accuracy in delineation of wetlands. Additionally, identification of wetlands in this study was based on the premise that surface depressions often give rise to temporary, seasonal or permanent wetlands, since runoff from upslope areas accumulates at downslope sites within the landscape. Thus, for study areas flatter than the site in this study, the method might have more difficulties in identifying wetlands. Furthermore, the study site is known to be a wetland-rich area, which increased the probability of actually identifying wetlands.
The approach used in this study for delineating potential wetlands using DEM was first developed by Wu and Lane [51] to detect actual wetland landscapes in the Pipestem river sub-basin in the prairie pothole region (PPR) of North Dakota. It is therefore not developed specifically to identify potential wetlands in a flat northern boreal landscape with a subarctic climate. However, there are some similarities between PPR and northern boreal landscapes, since PPR is characterised by small (median size 1600 m 2 ) and shallow (generally less than 1 m) wetlands [51], while the wetlands identified in this study had a median size of 2200 m 2 and a mean depth of 0.36 m.
The methodology can be used as a complement to other mapping methods, such as field work, aerial image interpretation, other DEM-based methods [36,37] or more automated remote sensing-based land cover classifications [34,35]. Moreover, since many of the existing topographic maps, such as the NLS Finland topographic database used as a source of reference data in this study, miss out wet areas, DEM-based methods can be used to improve the existing maps. In a study in Sweden, Lidberg et al. [37] used fieldwork data as training and validation data, and classified wetlands with the help of existing maps and topographic indices, including topographic wetness and depth to water. The inclusion of topographic indices increased the accuracy of mapping compared with existing maps, with overall accuracy increasing from 74% to 84% and producer accuracy increasing from 36% to 75%. High classification accuracies can also be obtained with different remote sensing datasets, in particular if different types of data are included in the classification task. In a study in northern Sweden, Karlson et al. [38] obtained classification accuracy of 91% when delineating peatlands with Sentinel-1 synthetic aperture radar and LIDAR DEM data. In some other studies, moderate to high producers' and total classification accuracies (>60%) have been obtained in wetland delineation and boreal wetland land cover classifications with different remote sensing datasets, such as optical, LIDAR and synthetic aperture radar data [36,65,66]. Accuracies can be much lower, even <20%, in wetland identification when using only a single data source [66], which shows that multiple data sources should be used in wetland delineation, including DEM, optical and synthetic aperture radar data. Most of the remote sensing-based wetland mappings require training data that can be collected through field work or an existing map of wetlands. However, some DEM-based approaches, such as topographic indices and the method developed in this study, do not necessarily require training data, which is definitely a strength of DEM-based methods. The method may be improved by combining it with other methods, such as the topographic wetness (besides using additional data). Future studies should test how the depression-based DEM method used in this study compares with other DEM and remote sensing-based methods in different landscapes, and how well it can complement other mapping and datasets. It is also suggested to combine DEM and other data to improve the method of mapping wetlands.
While the threshold values set may limit the identification of wetlands, the inclusion of additional factors that can indicate wetland, such as soil type, soil moisture, vegetation type and groundwater exchange, can relatively easily increase the potential of the methodology. In future research, the focus should be on developing the methodology in that regard. An advantage of the methodology is the use of DEM data, which are usually available. A number of countries collect high-resolution DEM, while more global high-resolution datasets such as the ArcticDEM [67] are also freely accessible to different practitioners, and thus can be used to identify wetlands. The use of only two parameters (minimum wetland size and depth) makes the method simple and, according to the results, enables it to identify wetlands in areas where other methods do not, specifically areas around lakes and ponds. Field work could be targeted at these areas to validate the findings. Despite the regional importance of Arctic wetlands, the datasets currently available (e.g., those reported by Arino et al. [68]; Chen et al. [69]; Kobayashi et al. [70]; Lehner and Döll [71]; Schroeder et al. [72] do not have enough accuracy and full coverage of the Arctic region. Due to substantial differences between wetlands, finding a standard method to outline all wetlands can be challenging, although there are several indicators available to define location of potential wetlands. Such indicators include different type of vegetation, hydrology, soil type, proximity to water, and landform [73] that should be delineated from multiple different data sources. Several databases, including national, regional, and global, to some extent can define wetland areas, but most available databases don't have sufficient resolution to accurately identify future vulnerable wetland areas. Moreover, none of the high-resolution datasets e.g., for soil type and soil moisture covers the entire Arctic region, indicating a need to consider several datasets.
The methodology developed in this study can also be used to increase understanding of wetland dynamics. For example, results could be regenerated using different DEMs collected at different times, such as multiple times during a year or from multiple years, and comparing these. Such analysis could show how the landscape and thereby the connectivity between wetlands changes, e.g., which flow paths are transient and which are more permanent. Furthermore, in flat landscapes such as the study area used in this research, catchment boundaries are not necessarily stable over time and water from a particular wetland might flow to different catchments, depending on the temporal dynamics.

Hydrological Connectivity of Wetlands and Wetland Management
Regarding the results from the connectivity analysis, delineation of the wetlandscape was relatively straightforward once the wetlands were identified and their catchments and drainage flow paths were delineated. Furthermore, the delineated flow paths indicated in many cases flow directions from the identified wetlands towards lakes or stream waters. These flow paths were only partly similar (11%) to streams mapped by NLS Finland, but were largely (61%) located in wetlands mapped by NLS Finland. This suggests that water flows and also spreads along the mapped flow paths. Hence, even if the identified wetlands do not cover all wetlands, the identified wetlands together with their delineated drainage flow paths define a more coherent spread of wetlands. However, it has to be considered that the delineation of the wetlandscape is carried out with the help of topographical data only and is based on the identified wetlands. The methodology has managed to delineate a wetlandscape in practice from a conceptual idea suggested by Thorslund et al. [38]. Further investigations regarding the mapped wetlandscapes and their specific characteristics were not conducted in this study. In future studies, the connectivity analysis could be extended by connecting the findings to data obtained with techniques such as dye tracers [39]. Moreover, the connectivity analysis could include groundwater fluctuations by monitoring the groundwater levels with pressure transducers [42].
Modelling hydrological connectivity is important, since it acts as the mechanism by which physical, biological and chemical materials are transported and distributed in the landscape [26,74].
These elements are crucial for the ecosystems that exist in and around wetlands, contributing more than 20% of the total value of ecosystem services globally [75]. These ecosystem services have the ability e.g., to sequester carbon, protect water quality, regulate flooding, support biodiversity and, not least, manage climate change. Understanding connectivity is a way to understand interactions between freshwater systems and prepare to meet current and future challenges, as has been especially emphasised in recent studies [2,38,[76][77][78]. The ability of the method to delineate wetlandscapes based on connectivity allows the total hydrological catchment of multiple wetlands to be analysed. The method could thus be used for tracking flows of water through a wetlandscape and identifying those that are more vulnerable and need more protection, e.g., flows that are more or less sensitive to changes as a result of climate change, land use or water use. Application of the method to compare DEMs acquired at different times could also be used to easily assess wetland connectivity over time, which will vary with climate change. In this way, the method can be used for remote monitoring of climate impacts on wetland connectivity. Another area of application where the method could be used is the poorly studied land-to-water (lateral) transport of carbon, which affects ecosystems and the biochemistry of freshwater [79]. Wetlands receive large amounts of carbon from the terrestrial environment and play a major role in transporting, transforming, storing and filtering water, sediments and solutes [80,81]. This functioning will be affected by climate-induced changes in the hydrological cycle [82,83], thereby altering the carbon cycle. The ability to map the connectivity between wetlands could be useful in research to quantify future transport of carbon and the carbon balance of watersheds.

Conclusions
In this research, DEM data were used to inventory wetlands in northern landscapes, and a depression-based method was used to apply the conceptual wetlandscape concept to wetland connectivity analysis. However, the wetlands identified in the study area did not cover the full extent of wetlands based on existing wetland maps. Thus the method cannot be used alone to map peatland and wetland areas in northern landscapes, and should be complemented with other datasets and mapping methods, such as field work and remote sensing-based mappings. Importantly, the delineated wetlandscapes achieved with the depression-based DEM method improved understanding of the hydrologically coupled system of several wetlands and their connectivity. Future research is required to improve the method, which could benefit from including more data, such as soil moisture and vegetation, and from validating the findings through field studies given the uncertainties associated with current NLS Finland map of wetlands. In the present application, the method was used to map wetlands and their hydrological connectivity, which is fundamental and valuable information for the management and protection of northern wetland landscapes.
Author Contributions: Z.K. initiated the research and planned it together with E.S. and A.R. E.S. conducted the data analyses under supervision by Z.K. and A.R. E.S. wrote the first draft of the manuscript and finalised it together with Z.K., C.S.S.F. and A.R. All authors have read and agreed to the published version of the manuscript.
Funding: This research was funded by Arctic Avenue and ASIAQ mobility grants from Stockholm University.