Development of an Integrated Modelling System for Evaluating Water Quantity and Quality Effects of Individual Wetlands in an Agricultural Watershed

A GIS-based fully-distributed model, IMWEBs-Wetland (Integrated Modelling for Watershed Evaluation of BMPs—Wetland), is developed to simulate hydrologic processes of site-specific wetlands in an agricultural watershed. This model, powered by the open-source GIS Whitebox Geospatial Analysis Tools (GAT) and advanced database technologies, allows users to simulate and assess water quantity and quality effects of individual wetlands at site and watershed scales. A case study of the modelling system is conducted in a subbasin of the Broughton’s Creek Watershed in southern Manitoba of Canada. Modelling results show that the model is capable of simulating wetland processes in a complex watershed with various land management practices. The IMWEBs-Wetland model is unique in simulating the water quantity and quality effects of individual wetlands, which can be used to examine location-specific targeting of wetland retention and restoration at a watershed scale.


Introduction
Wetland retention and/or restoration is an important best management practice (BMP) in agricultural watersheds as it provides critical hydrologic functions including flood attenuation, groundwater recharge, and contaminant filtering.Various hydrologic models have been applied to study wetland processes and evaluate their impacts on water quantity and quality at a watershed scale.However, most of these models have a semi-distributed structure, which lump all wetlands in a subbasin into a single functionally equivalent wetland for parameterization [1].Due to the complex and dynamic inter-connections of wetlands within a watershed, it is challenging to characterize the spatial heterogeneity of wetlands with a semi-distributed model.In recent years, several efforts have been made to address this challenge in watershed hydrologic modelling.These model enhancements have contributed to an improved understanding of wetland functions and services at a watershed scale.
Wang et al. [2] proposed a hydrologic equivalent wetland (HEW) concept in the Soil and Water Assessment Tool (SWAT) to characterize multiple wetlands within a subbasin.This approach lumps multiple wetlands within a subbasin into one wetland based on the aggregation of their geometric characteristics including wetland surface area, volume, and drainage area.The HEW function is achieved through calibrating wetland parameters including the fraction of the subbasin area that drains into wetlands, the volume of water stored in the wetlands when filled to their normal water level, the volume of water stored in the wetlands when filled to their maximum water level, the longest tributary channel length in the subbasin, and the Manning's n value for tributary and main channels.The HEW approach has the advantage of characterizing the non-linear functional relations between runoff and wetlands [3] and is reasonable for simulating hydrological processes of prairie pothole wetlands [4,5].Different from the HEW approach, Mekonnen et al. [6] proposed a probability distribution routine in the SWAT model to represent multiple wetlands within a subbasin.This approach applied a probability density function to characterize the spatial heterogeneity of landscape depression storages in a watershed and incorporated the seasonal variation of soil erodibility to account for the change of erosion rate during soil freeze and thaw.An application of the approach showed an improved simulation of sediment export in a Canadian prairie watershed [7].
Nasab et al. [8] developed an alternative depression characterization approach in the SWAT model to characterize multiple wetlands within a subbasin.This approach utilizes topographic characteristics and distribution of depressions to establish hierarchical relationships of depressions.This approach generates detailed hydro-topographic characteristics of depressions to parameterize SWAT pothole features by merging lower level depressions into higher level depressions.An application of the approach to a watershed with numerous potholes in North Dakota showed an improved model performance in simulating stream flows at the watershed outlet.Another approach in enhancing SWAT characterization of wetlands is to couple the SWAT model with a United States Environmental Protection Agency (USEPA) field scale model called System for Urban Stormwater Treatment and Analysis Integration (SUSTAIN) that characterizes drained water from upland flow [9,10].In addition to wetland modules developed in SWAT, improvements have also been made in the semi-distributed Soil and Water Integrated Model (SWIM) [11,12] and the distributed hydrological model HYDROTEL [13,14] through incorporating wetland flow, nutrient, and groundwater related dynamics to characterize wetland processes.
While the semi-distributed watershed models have had various enhancements and applications, there still exists the limitation on spatially explicit characterization of wetlands and their connectivity.In the SWAT model, areas with the same land use, soil and slope class in a subbasin are grouped into hydrologic response units (HRUs) to characterize their spatial heterogeneity.However, HRUs are not spatially connected within a subbasin.Evenson et al. [15,16] developed an alternative approach to characterize geographically isolated wetlands (GIWs) by redefining SWAT HRUs to conform with the mapped GIWs and their drainage boundaries.New model input files were constructed to direct the simulation of GIW fill-spill hydrology and upland flows to GIWs.The enhanced SWAT was applied to a North Carolina watershed and a North Dakota watershed to examine wetland effects on stream flow, baseflow and peak flows with a satisfactory model performance at the watershed outlet.Modelling application to a coastal plain watershed in North Carolina showed that increased extent of isolated and riparian wetlands had significant effects on decreasing seasonal and annual flows [17].This approach of wetland representation in the SWAT model contributes to a more realistic simulation of wetland effects.However, the lumping at the subbasin level due to inherent semi-distributed model structure prevents the examination of wetland effects at specific locations.There is a need to develop a fully distributed watershed model to explicitly characterize individual wetlands and simulate their effects on flow and water quality at a watershed scale.
In the period of 2004-2013, Agriculture and Agri-Food Canada's (AAFC) Watershed Evaluation of BMPs (WEBs) program conducted BMP assessments in nine experimental watersheds across Canada [18].In the WEBs program, the Guelph Watershed Evaluation Group developed a cell-based integrated modelling system for watershed evaluation of BMPs (IMWEBs) to conduct BMP assessments at site, field, farm, and watershed scales [19].The foundation of IMWEBs was the Water and Energy Transfer between Soil, Plant and Atmosphere (WetSpa) model, a fully distributed hydrologic model for flood prediction and watershed management [20].Several other well-known agricultural watershed models such as SWAT [21], Agricultural Policy Environmental Extender (APEX) [22] and Riparian Ecosystem Management Model (REMM) [23] were also referenced in the IMWEBs development.The IMWEBs model has an object-oriented and modular-based structure to perform Water 2018, 10, 774 3 of 18 dynamic watershed modelling supported by five databases including geospatial, hydro-climate, BMP, parameter, and model output.Hydrologic processes simulated in the model include climate, snowmelt, water balance, plant growth, soil erosion, nutrient cycle, and channel routing, while general agricultural BMPs including crop management, fertilizer management, and tillage management are incorporated in the model simulation.The IMWEBs model was firstly applied to the 75 km 2 South Tobacco Creek watershed in southern Manitoba of Canada and satisfactory modelling results were obtained [19].
The objective of this paper is to present the development of a wetland module as one of the BMP components in the IMWEBs model (IMWEBs-Wetland).A case study in a small agricultural watershed is also presented to demonstrate its applicability, performance, and usefulness in characterizing wetlands in an agricultural watershed.With a cell-based structure, the IMWEBs-Wetland model is specifically designed for simulating and evaluating water quantity and water quality effects of individual wetlands at a watershed scale.This makes the model more convenient and straightforward for supporting site-specific wetland retention and restoration.The modelling system also has novelties of using open-source GIS and databases and a modular structure to characterize different hydrologic processes.An object-oriented computer interface is developed to facilitate watershed delineation, model setup, model parameterization, scenario design, and output analysis.The model was applied to a 15.7-km 2 watershed in southern Manitoba of Canada to evaluate wetland effects.In addition to simulating the hydrologic processes operating in each wetland, the IMWEBs-Wetland model can also produce spatial distributions of various hydrologic variables at user-defined spatial and temporal scales.This makes the model an effective tool for spatial watershed management, particularly the assessment of site-specific wetland retention and restoration scenarios.Discussions of the model performance and future development perspectives are provided at the end of this paper.

System Design
In addition to other existing modules in the IMWEBs modular library, IMWEBs-Wetland is designed to simulate and assess wetland hydrologic processes at both individual and watershed scales.Processes simulated in the wetland module include water balance, sediment balance, and nutrient balance.In this study, wetlands that are isolated from mainstreams and dominated by fill and spill processes [1] are simulated in the IMWEBs-Wetland model, while riparian wetlands along mainstreams are simulated with channel processes in the model.
In the IMWEBs-Wetland model, a watershed is delineated into subbasins, and one subbasin contains one isolated wetland at the subbasin outlet.These isolated wetlands are spatially connected in the model through a DEM delineated flow path once they are filled.Other subbasin outlets can be defined at major tributary confluences, monitoring stations, watershed outlet, and user defined locations assisted by the modelling interface.In addition to the general IMWEBs inputs including climate, DEM, soil, land use, and land management, wetland inventory data with associated wetland parameters are required for model setup.Wetland parameters such as surface area, volume, and drainage area are estimated based on available DEM and wetland inventory data with the IMWEBs-Wetland interface prior to model simulation.Outputs from the IMWEBs-Wetland model include time series of flow, sediment, and nutrient concentration at any user defined wetland or reach locations, and spatial distribution of wetland or watershed hydrologic variables at user-defined spatial and temporal scales determined during model setup.A diagram for a delineated wetland subbasin characterized in the IMWEBs-Wetland model is shown in Figure 1.
To follow the modelling structure of the IMWEBs model, cells within a wetland polygon are simulated the same as upland cells and are not grouped into one unit.Processes of runoff and pollutant generation for cells within the wetland polygon are simulated using IMWEBs algorithms based on the land cover and DEM data.However, soil parameters, particularly the wetland bottom hydraulic conductivity are modified during IMWEBs-Wetland model setup.As such, the wetland module is on the top of the existing IMWEBs model with inputs from its drainage areas, and outputs the same as those of the reservoirs.Wetlands are connected through surface water and groundwater in the model.Surface water flows along the pathway derived from the DEM once the wetland is filled or above its normal storage for drained-altered wetlands.Groundwater is simulated separately from the wetland module at a subbasin scale using a non-linear reservoir method based on its contributing area.Descriptions on groundwater simulation in the model are provided in Sections 2.3 and 2.5.To follow the modelling structure of the IMWEBs model, cells within a wetland polygon are simulated the same as upland cells and are not grouped into one unit.Processes of runoff and pollutant generation for cells within the wetland polygon are simulated using IMWEBs algorithms based on the land cover and DEM data.However, soil parameters, particularly the wetland bottom hydraulic conductivity are modified during IMWEBs-Wetland model setup.As such, the wetland module is on the top of the existing IMWEBs model with inputs from its drainage areas, and outputs the same as those of the reservoirs.Wetlands are connected through surface water and groundwater in the model.Surface water flows along the pathway derived from the DEM once the wetland is filled or above its normal storage for drained-altered wetlands.Groundwater is simulated separately from the wetland module at a subbasin scale using a non-linear reservoir method based on its contributing area.Descriptions on groundwater simulation in the model are provided in Sections 2.3 and 2.5.

Wetland Classification and Characterization
Five types of wetlands including drained-altered, drained-consolidated, drained-lost, undrained-altered, and undrained-intact are classified and simulated in the IMWEBs-Wetland model based on the Ducks Unlimited Canada (DUC) wetland inventory data (Figure 2).The drainedaltered wetland has an outlet drain, while drained-consolidated wetland has an inlet drain, both with riparian and aquatic vegetation present.The drained-lost wetland has an outlet drain without riparian and aquatic vegetation present.The undrained-altered wetland does not have drains with riparian and aquatic vegetation absent or disturbed, while undrained-intact wetland does not have drains with riparian and aquatic vegetation present.The drained-lost and drained-altered wetlands are similar in that they are all drained but one is altered for cultivation and the other is vegetated with native plants.This multi-temporal classification approach defines wetlands that are in either drained or undrained state.The wetland classification provides a basis for identifying existing wetlands for retention and drained-lost wetlands for restoration during IMWEBs-Wetland scenario assessment.

Wetland Classification and Characterization
Five types of wetlands including drained-altered, drained-consolidated, drained-lost, undrained-altered, and undrained-intact are classified and simulated in the IMWEBs-Wetland model based on the Ducks Unlimited Canada (DUC) wetland inventory data (Figure 2).The drained-altered wetland has an outlet drain, while drained-consolidated wetland has an inlet drain, both with riparian and aquatic vegetation present.The drained-lost wetland has an outlet drain without riparian and aquatic vegetation present.The undrained-altered wetland does not have drains with riparian and aquatic vegetation absent or disturbed, while undrained-intact wetland does not have drains with riparian and aquatic vegetation present.The drained-lost and drained-altered wetlands are similar in that they are all drained but one is altered for cultivation and the other is vegetated with native plants.This multi-temporal classification approach defines wetlands that are in either drained or undrained state.The wetland classification provides a basis for identifying existing wetlands for retention and drained-lost wetlands for restoration during IMWEBs-Wetland scenario assessment.
The drained-consolidated and undrained-intact wetlands have similar geometric features for which no outlet drains exist, and flow out of the wetlands follows the fill-and-spill process.For these two types of wetlands, the normal wetland surface area and normal storage are assumed equal to the maximum wetland surface area and maximum wetland storage.The following volume-area regression equations are used to estimate wetland storage based on the wetland surface area [24]: V = 2.85A 1.22 when A ≤ 70 ha (1) where A is the wetland surface area (ha) and V is the corresponding full supply volume (10 3 m 3 ).The constant in Equations ( 1) and ( 2) can be adjusted for specific wetlands if observation data are available.Discharge out of the wetland is calculated after the wetland is filled, and all excess water discharges into the downstream reach during the time step.
Water 2018, 10, 774 5 of 18 The drained-altered wetland is different from above in that the wetland is drained and altered but still maintain water in the wetland with a reduced holding capacity.Flow out of the wetland occurs when water level in the wetland is higher than outlet drains.To characterize this type of wetland, the maximum storage is calculated using Equations ( 1) and ( 2) based on the wetland inventory data, while the normal wetland surface area is assumed to be 1/3 of the maximum surface area if no field survey data are available, and the corresponding normal storage is calculated using Equations ( 1) and (2).Discharge out of the wetland when water volume is between normal storage and maximum storage is calculated using equation ( 3) from the SWAT model [25]: where V flowin and V flowout are the water volume entering and flowing out of the wetland over the time step (m 3 ), V 1 is the end storage of previous time step (m 3 ), V normal is the wetland normal storage (m 3 ), and C is a wetland discharge coefficient with a default value of 10.No outflow occurs if sum of inflow and existing storage is less than normal wetland storage.When calculated wetland water volume is over maximum storage, all excess water discharges during the time step.The wetland discharge coefficient can be adjusted for specific wetlands when field observation data are available.
Water 2018, 10, x FOR PEER REVIEW 5 of 19 The drained-consolidated and undrained-intact wetlands have similar geometric features for which no outlet drains exist, and flow out of the wetlands follows the fill-and-spill process.For these two types of wetlands, the normal wetland surface area and normal storage are assumed equal to the maximum wetland surface area and maximum wetland storage.The following volume-area regression equations are used to estimate wetland storage based on the wetland surface area [24]: where A is the wetland surface area (ha) and V is the corresponding full supply volume (10 3 m 3 ).The constant in Equations ( 1) and ( 2) can be adjusted for specific wetlands if observation data are available.Discharge out of the wetland is calculated after the wetland is filled, and all excess water discharges into the downstream reach during the time step.The drained-altered wetland is different from above in that the wetland is drained and altered but still maintain water in the wetland with a reduced holding capacity.Flow out of the wetland occurs when water level in the wetland is higher than outlet drains.To characterize this type of wetland, the maximum storage is calculated using Equations ( 1) and ( 2) based on the wetland inventory data, while the normal wetland surface area is assumed to be 1/3 of the maximum surface area if no field survey data are available, and the corresponding normal storage is calculated using Equations ( 1) and (2).Discharge out of the wetland when water volume is between normal storage and maximum storage is calculated using equation ( 3) from the SWAT model [25]: where Vflowin and Vflowout are the water volume entering and flowing out of the wetland over the time step (m 3 ), V1 is the end storage of previous time step (m 3 ), Vnormal is the wetland normal storage (m 3 ), and C is a wetland discharge coefficient with a default value of 10.No outflow occurs if sum of inflow and existing storage is less than normal wetland storage.When calculated wetland water volume is

Water Balance
The water balance for a wetland is: where V 2 is the water volume in the wetland at the end of the time step (m 3 ), V pcp is the precipitation volume on the wetland over the time step (m 3 ), V evap is the evaporation/evapotranspiration from the wetland surface over the time step (m 3 ), and V seep is the water volume lost from the wetland by seepage (m 3 ).To avoid duplicate calculation, V pcp , V evap , and V seep are calculated in the IMWEBs-Wetland model and summed for cells within the wetland polygon.For each wetland cell, based on land cover and wetland inventory data, V evap is estimated with a potential evaporation/evapotranspiration rate in which evaporation is calculated for wetlands with an open water surface, and evapotranspiration is calculated for wetlands with a vegetation cover including forest.V seep is calculated based on wetland bottom conductivity and water availability calculated in the wetland module.Drained-consolidated wetland is different from drained-altered or drained-lost wetlands in that water is drained through surface inlet drains connected to a tile underlying the wetland rather than at the wetland outlet.To simulate this process in the model, flow at the subbasin outlet is calculated Water 2018, 10, 774 6 of 18 by summing surface runoff from non-wetland cells within the subbasin and bypasses the wetland into the downstream reach.Water balance in the consolidated wetland is maintained with inputs of precipitation in the wetland area and lateral flow from upland fields, while outflow from the wetland occurs when the wetland is filled, and evapotranspiration and seepage are calculated the same as those for drained-altered wetlands.Output from the wetland including flow, sediment and nutrient loading replaces the original reach output of the subbasin and joins the routing of downstream reaches.

Sediment and Nutrient Balance
A mass balance approach in the SWAT model [25] is used for wetland sediment and nutrient routing.The governing equation is: For sediment balance, M 2 and M 1 are the amount of sediment in the wetland at the end and at the beginning of the time step (ton), M flowin is the amount of sediment added to the wetland with inflow (ton), M flowout is the amount of sediment transported out of the wetland with outflow (ton), and M stl is the amount of sediment removed from the water body by settling (ton).M flowin is obtained from IMWEBs-Wetland reach output of the subbasin.As the wetland is located at the outlet of a subbasin, all cells within the subbasin including wetland area contribute runoff, sediment and nutrient yields at the subbasin outlet.Calculation of sediment deposition is based on actual sediment concentration and equilibrium sediment concentration of the wetland using the SWAT approach [25].Sed flowout is the suspended sediment out of the wetland calculated by outflow multiplied by sediment concentration.The calculated sediment loading replaces the reach sediment output and joins the routing of downstream reaches.This method assumes a uniform distribution of deposited sediment on the wetland bottom without accounting for its spatial distributions.
The SWAT wetland water quality algorithms are used to simulate nutrient balance for each wetland.The mass balance for nitrogen and phosphorous is similar as the sediment mass balance, where M 2 and M 1 in Equation ( 5) are the amount of nutrients in the wetland at the end and at the beginning of the time step (kg), M flowin is the amount of nutrients added to the wetland with inflow (kg), M flowout is the amount of nutrients transported out of the wetland with outflow (kg), and M stl is the amount of nutrients removed from the water body by settling and seepage (kg).Detailed descriptions of the method can be found in [25].

Groundwater
In the IMWEBs model, baseflow is calculated at subbasin scale using a non-linear reservoir method.This approach is applicable for a relatively large subbasin.However, when delineated subbasin is very small, the groundwater flow may not contribute to the subbasin outlet but at a location in downstream reaches.This pattern causes challenge for groundwater characterization in the IMWEBs-Wetland model, because individual wetlands are located at subbasin outlets, and their drainage areas may be very small for a prairie watershed with numerous pothole wetlands.As a result, groundwater may bypass the outlet from underground and does not join the flow at the subbasin outlet.
To solve this problem, a threshold drainage area (ha) is incorporated in the IMWEBs-Wetland model for groundwater simulation at the reach outlet.This threshold area can be estimated based on field observations at a site where groundwater flow is initiated during snowmelt period or after heavy storms.It can be also determined through model calibration if flow data are available at monitoring stations with small contribution areas within the watershed.If the reach contribution area is less than the threshold value, no groundwater outflow exists at the reach outlet, and the calculated groundwater flow for this subbasin is accumulated to the next downstream reach.Once the reach drainage area is greater than the threshold value, the accumulated groundwater flow returns to the channel.The groundwater flow calculated for this subbasin and the groundwater flow accumulated from upstream subbasins are summed and added to the inflow of the river reach for channel routing.These assumptions are workable for upstream subbasins based on reach drainage areas calculated during model setup.For small subbasins in middle and downstream areas adjacent to the main channel, if the subbasin area is less than the threshold value, the calculation follows the same way as upstream subbasins, and the calculated groundwater flow is added to the main channel outlet.

Framework and Database Management
A computer interface was developed to assist IMWEBs-Wetland input data preparation, watershed delineation, model setup, parameterization, and result visualization.The interface, powered by the Whitebox Geospatial Analysis Tools (GAT) and SQLite database technologies, facilitates the user to simulate water quantity and quality effects of individual wetlands at site, field, farm and watershed scales.Whitebox GAT is an open-source desktop GIS and remote sensing software package for general applications of geospatial analysis and data visualization and is intended to provide a platform for advanced geospatial data analysis with applications in both environmental research and the geomatics industry [26].The purpose of the interface is to conduct pre-and post-processing for the IMWEBs-Wetland model.These include: (a) to delineate watershed subbasins accounting for each wetland; (b) to compute model parameters for each wetland; (c) to display wetland and watershed results; and (d) to manage all associated wetland and watershed input, output, and parameter databases.
The IMWEBs-Wetland model has five databases and one modular library (Figure 3).The geospatial database is a collection of geometric entities and attributes of the watershed such as DEM, soil, land use, streams, boundaries, climate and hydrologic stations, wetlands, reservoirs, and hydraulic structures.

Drainage Delineation
To simulate the hydrologic processes of individual wetlands, their contribution areas and the drainage network that links each wetland need to be properly defined to characterize interconnections of wetlands during model simulation.Traditionally, watershed drainage delineation requires a filling algorithm to remove depressions and flat areas in the watershed so that a continuous stream network can be generated using a GIS.This approach eliminates the depression

Drainage Delineation
To simulate the hydrologic processes of individual wetlands, their contribution areas and the drainage network that links each wetland need to be properly defined to characterize interconnections of wetlands during model simulation.Traditionally, watershed drainage delineation requires a filling algorithm to remove depressions and flat areas in the watershed so that a continuous stream network can be generated using a GIS.This approach eliminates the depression information, e.g., wetland, in a grid DEM, which is not appropriate for wetland drainage delineation.
In the IMWEBs-Wetland model interface, a customized watershed delineation method was developed to solve this problem.Two types of input data, wetland boundary and wetland outlet (optional), are incorporated in the delineation process.If wetland outlet location data are not available, the interface can automatically generate the outlet information using the priority-flood operation approach [27] based on the DEM.For a large wetland with surface area greater than a user defined threshold value, multiple outlets are allowed for an individual wetland during the delineation process.The wetland boundary and outlet vector layers are firstly rasterized into grids containing wetland identifier values.Combining with the DEM, a D8 flow-direction raster is generated and is used to calculate flow accumulation and generate a drainage network.
Figure 4 illustrates how this updated flow direction raster is created based on the DEM, wetland boundary, and wetland outlet information.For non-boundary cells (white cells in Figure 4), flow directions are determined by examining grid elevation values using the D8 method [28], i.e., flow direction is pointed to the steepest downward neighbor cell based on elevation.By overlaying wetland boundary with the DEM, wetland boundary cells (light gray cells in Figure 4) are detected.For these wetland boundary cells, flow directions are enforced to circle the wetland within these boundary cells until a wetland outlet cell (dark gray cell in Figure 4) is reached.The subdivided wetlands within the original large wetland are treated as individual wetlands with their own outlets and flow pathways to the downstream reaches in the IMWEBs-Wetland model.Because the model generates spatial distribution results for each wetland, outputs of these subdivided wetlands are summarized automatically back after model simulation using the weighted average approach based on their surface areas within the original large wetland.

Wetland Parameterization
Parameters for each wetland are estimated based on wetland inventory data and are prepared after watershed delineation.Table 1 summarizes the wetland parameter name, units and their descriptions.Wetland ID, type, and maximum surface area are read from the wetland inventory attribute table.The parameter of operation year is used in the IMWEBs-Wetland model for simulating wetland loss and restoration scenarios.Parameters of subbasin ID and contribution area are obtained from the watershed delineation results.Normal water volume and maximum water volume are

Wetland Parameterization
Parameters for each wetland are estimated based on wetland inventory data and are prepared after watershed delineation.Table 1 summarizes the wetland parameter name, units and their descriptions.Wetland ID, type, and maximum surface area are read from the wetland inventory attribute table.The parameter of operation year is used in the IMWEBs-Wetland model for simulating wetland loss and restoration scenarios.Parameters of subbasin ID and contribution area are obtained from the watershed delineation results.Normal water volume and maximum water volume are calculated using Equations ( 1) and (2), while normal water volume corresponds to the normal storage for drained-altered wetlands over which spill flow would occur.Other parameters are read from the default wetland parameter table in the BMP database.All these parameter values can be adjusted during model calibration or re-edited if field observation data are available.

Study Area
A case study of IMWEBs-Wetland modelling was conducted in a 15.7 km 2 small watershed, which is a subbasin of the Broughton's Creek watershed in southern Manitoba of Canada (Figure 5).The Broughton's Creek flows southeasterly into the Little Saskatchewan River, joining the Assiniboine River, and eventually entering Lake Winnipeg.The watershed has an average slope of 1.40% with a range from 0.0% to 10.8% based on a 30 m DEM and is dominated by the moderately well drained Newdale soils formed in calcareous, loamy glacial till of limestone, granite and shale origin.Agriculture is the major land use in the watershed (75.3%) with dominant crop types of spring wheat and canola.Other land use types include grassland (3.4%), wetland (20.4%), road (0.7%), and deciduous forest (0.2%).The study watershed has hundreds of undrained depressions ranging from potholes to large sloughs.Compared with potholes, sloughs are relatively large, shallow, and typically elongated northwest to southeast.In addition, there are also several small lakes present.Limited by available detailed information, this study did not differentiate potholes, sloughs, and lakes, and modeled them as wetlands.Based on the DUC wetland inventory data, a total of 492 wetlands were identified in the study watershed including 293 existing (Figure 6) and 199 drained-lost.Among the 293 existing wetlands, 85 were identified as drained-altered, 11 drained-consolidated, 133 undrained-altered, 46 undrained-intact, and 18 riparian wetlands.The total wetland area including drained-lost is 550 ha, which means 35.0% of the watershed was covered by wetland before agricultural development.Detailed descriptions of their surface areas and volumes are provided in Table 2.
The study area has a semi-arid climate, with a pronounced seasonal variation in temperature and precipitation.Based on the 1981-2010 climate data recorded at Strathclair station, located about 10 km northwest of the study watershed (Figure 5), the average yearly temperature was 1.5 • C, with the highest monthly temperature of 17.7 • C in July and the lowest monthly temperature of −17.1 • C in January.Average annual precipitation was 475 mm, of which 118 mm (25%) was snowfall, lasting from November to the following April.The average annual daily discharge at the watershed outlet was 0.065 m 3 /s, ranging from 0.00 to 1.57 m 3 /s based on the observed data collected at the EC9 station (Figure 5) from 2009 to 2013.The average annual runoff was 140 mm, with an average runoff coefficient of 0.29.More than 80% of this runoff and all annual peak discharges were observed in spring (late April to early May) over the monitoring period because of snowmelt.Baseflow was a small portion of the total runoff (<10%) and provided little contribution to the total flow and sediment transport.
The study watershed has suffered flooding and nutrient loading problems during spring snowmelt and heavy summer storms due to activities of wetland drainage, road construction, and land clearing for agricultural production.Therefore, retention of existing wetlands and restoration of drained-lost wetlands are important for addressing these environmental problems.The entire Broughton's Creek watershed had been selected to study the effects of wetland loss and restoration on stream water quality using the SWAT model [4,5].However, due to the SWAT's semi-distributed model structure, the assessment of wetland loss and restoration effects was conducted at a subbasin scale but not for individual wetlands.

Study Area
A case study of IMWEBs-Wetland modelling was conducted in a 15.7 km 2 small watershed, which is a subbasin of the Broughton's Creek watershed in southern Manitoba of Canada (Figure 5).The Broughton's Creek flows southeasterly into the Little Saskatchewan River, joining the Assiniboine River, and eventually entering Lake Winnipeg.The watershed has an average slope of 1.40% with a range from 0.0% to 10.8% based on a 30 m DEM and is dominated by the moderately well drained Newdale soils formed in calcareous, loamy glacial till of limestone, granite and shale origin.Agriculture is the major land use in the watershed (75.3%) with dominant crop types of spring wheat and canola.Other land use types include grassland (3.4%), wetland (20.4%), road (0.7%), and deciduous forest (0.2%).The study watershed has hundreds of undrained depressions ranging from potholes to large sloughs.Compared with potholes, sloughs are relatively large, shallow, and typically elongated northwest to southeast.In addition, there are also several small lakes present.Limited by available detailed information, this study did not differentiate potholes, sloughs, and lakes, and modeled them as wetlands.Based on the DUC wetland inventory data, a total of 492 wetlands were identified in the study watershed including 293 existing (Figure 6) and 199 drainedlost.Among the 293 existing wetlands, 85 were identified as drained-altered, 11 drainedconsolidated, 133 undrained-altered, 46 undrained-intact, and 18 riparian wetlands.The total wetland area including drained-lost is 550 ha, which means 35.0% of the watershed was covered by wetland before agricultural development.Detailed descriptions of their surface areas and volumes are provided in Table 2.The study area has a semi-arid climate, with a pronounced seasonal variation in temperature and precipitation.Based on the 1981-2010 climate data recorded at Strathclair station, located about 10 km northwest of the study watershed (Figure 5), the average yearly temperature was 1.5 °C, with the highest monthly temperature of 17.7 °C in July and the lowest monthly temperature of −17.1 °C in January.Average annual precipitation was 475 mm, of which 118 mm (25%) was snowfall, lasting from November to the following April.The average annual daily discharge at the watershed outlet was 0.065 m 3 /s, ranging from 0.00 to 1.57 m 3 /s based on the observed data collected at the EC9 station (Figure 5) from 2009 to 2013.The average annual runoff was 140 mm, with an average runoff coefficient of 0.29.More than 80% of this runoff and all annual peak discharges were observed in spring (late April to early May) over the monitoring period because of snowmelt.Baseflow was a small portion of the total runoff (<10%) and provided little contribution to the total flow and sediment transport.
The study watershed has suffered flooding and nutrient loading problems during spring snowmelt and heavy summer storms due to activities of wetland drainage, road construction, and land clearing for agricultural production.Therefore, retention of existing wetlands and restoration of drained-lost wetlands are important for addressing these environmental problems.The entire Broughton's Creek watershed had been selected to study the effects of wetland loss and restoration on stream water quality using the SWAT model [4,5].However, due to the SWAT's semi-distributed

Model Setup and Calibration
The IMWEBs-Wetland model for the study watershed was setup based on the geospatial data of DEM, land use, soil, and wetland inventory obtained from the DUC.A total of 515 subbasins were delineated, of which 511 have wetlands at the subbasin outlets and four are in the mainstreams.The original land use layer has 15 land use classes ranging from agricultural cropland to roads and trails.To convert this land use layer into IMWEBs-Wetland format, a lookup table was created linking original land use categories to the IMWEBs-Wetland land use code.Accordingly, a user-defined soil parameter database was developed based on available soil attribute data.Crop management is one of the key factors in controlling runoff, sediment and nutrient yields from an agricultural watershed.Nitrogen (N) and phosphorus (P) from agricultural land are major sources to the receiving wetlands and streams.Because no detailed crop management data were available, a two-year representative crop rotation of spring wheat and canola was assumed for the study watershed (Table 3).Crop management parameters including seeding and harvest date, fertilizer application rate and date, and tillage type and date were referenced from available literature values [29,30].
The calibration period for the study watershed was from 2009 to 2013 at a daily scale, whereas the period from 2000 to 2008 was used for warming up.This period was selected for calibration because observed flow and water quality data were available at the subbasin outlet EC9 station (Figure 5).Precipitation and temperature data over the simulation period were obtained from the Strathclair and Rivers stations, while wind speed and wind direction data were obtained from the Brandon-A station (Figure 5).Data of solar radiation and relative moisture for the study area were downloaded from NASA's online database [31].The solar radiation data was further validated using an image from Energy, Mines, and Resources Canada [32].A manual calibration was conducted for those parameters deemed most sensitive based on a parameter sensitivity analysis of the IMWEBs-Wetland model.These include runoff and water balance parameters (interception capacity for different land covers, evapotranspiration correction factor, interflow scaling factor, field capacity and porosity of top soil layer, baseflow constant and exponent, potential surface runoff coefficient, snowmelt threshold temperature and degree-day coefficient, frozen soil moisture and temperature), soil erosion and sediment transport parameters (soil erodibility and practice factor in the universal soil loss equation, stream flow peak rate adjustment factor, critical velocity for channel erosion, constant and exponent for calculating the maximum amount of sediment that can be transported in a reach segment), and nutrient yield parameters (initial soil NO 3 and soluble P concentration, organic N and organic P enrichment ratio, phosphorous soil partitioning and percolation coefficient, and nitrate percolation coefficient).Other parameters were set to their default values and were not adjusted during the process of model calibration.Because the purpose of this paper is to introduce the IMWEB-Wetland model development and its modelling ability, detailed descriptions of model calibration for the study watershed are not given here.
For calibration of flow, sediment, and nutrient parameters, model performance was evaluated graphically and statistically based on model bias (BIAS), Nash-Sutcliffe coefficient (NSC), root mean square error (RMSE), and correlation coefficient (CORR) at the EC9 station.BIAS can be expressed as the relative mean difference between observed and predicted results reflecting the ability of reproducing water, sediment and nutrient balance.NSC describes how well the predictions are produced by the model that is commonly used for model evaluation [33].The calibration objective for flow was to maximize NSC and CORR while simultaneously attempting to reduce BIAS.Calibration of sediment, total nitrogen (TN) and total phosphorous (TP) were conducted for their loadings on sampling days.Observed sediment loading was calculated by multiplying observed sediment concentration by observed flow of the day.Observed TN and TP loadings were calculated by multiplying sampled TN and TP concentrations by observed flow of the day.The model was calibrated firstly for stream flow, then sediment, and finally TN and TP.A summary of model performance at EC9 station for the 2009-2013 calibration period is provided in Table 4, and a graphical comparison between observed and simulated flow at the EC9 station is shown in Figure 7. Overall, the IMWEBs-Wetland simulated stream flows, sediment and nutrient loadings at the outlet station matched the observed data reasonably well based on the statistical assessment results and graphical comparisons.
Under existing condition over the period of 2009-2013, the model predicted an average runoff of 127 mm/year with a runoff coefficient of 0.24, sediment loading 0.04 t/ha, TN loading 2.56 kg/ha, and TP loading 0.47 kg/ha at the watershed outlet (Table 5).Figure 8 shows the simulated spatial distribution of TP yield in the study watershed for the year 2010.Clearly, higher TP losses were from crop lands over the watershed, in which the highest TP losses were in areas with steep slopes.A spatial distribution of wetland TP concentration for the year 2010 is given in Figure 9.The simulated TP concentration was highly variable among the existing wetlands ranging from 0.0 to 0.15 mg/L depending upon their geometric characteristics, the size of upstream contribution areas and their land management practices.However, because no wetland monitoring data were available in the study watershed, the model calibration was conducted at the watershed outlet but not at specific wetland sites.This would cause uncertainties for the wetland modelling results due to input data, model structure, and model parameter estimation.produced by the model that is commonly used for model evaluation [33].The calibration objective for flow was to maximize NSC and CORR while simultaneously attempting to reduce BIAS.Calibration of sediment, total nitrogen (TN) and total phosphorous (TP) were conducted for their loadings on sampling days.Observed sediment loading was calculated by multiplying observed sediment concentration by observed flow of the day.Observed TN and TP loadings were calculated by multiplying sampled TN and TP concentrations by observed flow of the day.The model was calibrated firstly for stream flow, then sediment, and finally TN and TP.A summary of model performance at EC9 station for the 2009-2013 calibration period is provided in Table 4, and a graphical comparison between observed and simulated flow at the EC9 station is shown in Figure 7.
Overall, the IMWEBs-Wetland simulated stream flows, sediment and nutrient loadings at the outlet station matched the observed data reasonably well based on the statistical assessment results and graphical comparisons.Under existing condition over the period of 2009-2013, the model predicted an average runoff of 127 mm/year with a runoff coefficient of 0.24, sediment loading 0.04 t/ha, TN loading 2.56 kg/ha, and TP loading 0.47 kg/ha at the watershed outlet (Table 5).Figure 8 shows the simulated spatial distribution of TP yield in the study watershed for the year 2010.Clearly, higher TP losses were from crop lands over the watershed, in which the highest TP losses were in areas with steep slopes.A spatial distribution of wetland TP concentration for the year 2010 is given in Figure 9.The simulated TP concentration was highly variable among the existing wetlands ranging from 0.0 to 0.15 mg/L depending upon their geometric characteristics, the size of upstream contribution areas and their land management practices.However, because no wetland monitoring data were available in the study watershed, the model calibration was conducted at the watershed outlet but not at specific wetland sites.This would cause uncertainties for the wetland modelling results due to input data, model structure, and model parameter estimation.

Scenario Development and Assessment
To demonstrate the ability of the IMWEBs-Wetland model for assessing the effects of wetlands on runoff, sediment and nutrient yields, five scenarios were constructed as listed in Table 5. Scenario I is the baseline scenario with all existing wetlands and existing land management practices.Scenario II is an extreme scenario assuming all existing wetlands are drained and lost for cultivation.Scenario III is another extreme scenario assuming all drained-lost wetlands are restored.Scenario IV and V

Scenario Development and Assessment
To demonstrate the ability of the IMWEBs-Wetland model for assessing the effects of wetlands on runoff, sediment and nutrient yields, five scenarios were constructed as listed in Table 5. Scenario I is the baseline scenario with all existing wetlands and existing land management practices.Scenario II is an extreme scenario assuming all existing wetlands are drained and lost for cultivation.Scenario III is another extreme scenario assuming all drained-lost wetlands are restored.Scenario IV and V are two spatial targeting scenarios for TP reduction to identify the top 10 most effective wetlands from the existing wetlands for retention and to identify the top 10 most effective ones from the drained-lost wetlands for restoration in the study watershed.The selection of these wetlands was performed by calculating the TP reduction efficiency with the equation: where TP_E is the TP reduction efficiency of the wetland (kg/year/ha), TP_R is the annual average wetland TP reduction (kg/year) calculated by subtracting outflow TP from inflow TP of the wetland modelling outputs, and Wet_A is the wetland surface area (ha).The selected 10 most effective existing wetlands for retention and 10 most effective drained-lost wetlands for restoration are shown in Figure 10.The assessment of wetland loss/restoration impacts on runoff, sediment and nutrient loading for the study watershed was performed based on the 2009-2013 climate and land management data and by changing the wetland operation year in the wetland BMP database.Modelling results were compared to the baseline scenario and are provided in Tables 5 and 6, respectively.Scenario II (loss of all existing wetlands) would increase total runoff at the watershed outlet by 21.3%, sediment by 35.4%, TN by 23.1%, and TP by 20.6% respectively.Scenario III (restoration of all drained-lost wetlands) would decrease total runoff at the watershed outlet by 18.9%, sediment by 29.8%, TN by 17.7%, and TP by 14.5%, respectively.The assessment of wetland loss/restoration impacts on runoff, sediment and nutrient loading for the study watershed was performed based on the 2009-2013 climate and land management data and by changing the wetland operation year in the wetland BMP database.Modelling results were compared to the baseline scenario and are provided in Tables 5 and 6, respectively.Scenario II (loss of all existing wetlands) would increase total runoff at the watershed outlet by 21.3%, sediment by 35.4%, TN by 23.1%, and TP by 20.6% respectively.Scenario III (restoration of all drained-lost wetlands) would decrease total runoff at the watershed outlet by 18.9%, sediment by 29.8%, TN by 17.7%, and TP by 14.5%, respectively.Scenario IV (loss of 10 most effective wetlands) would increase runoff by 3.94%, sediment by 6.07%, TN by 3.98%, and TP by 3.23%, respectively.Scenario V (restoration of 10 most effective wetlands) would decrease runoff by 7.87%, sediment by 13.0%, TN by 8.21%, and TP by 6.33% respectively.The full wetland loss scenario (II) and full restoration scenario (III) have average TN and TP increase rates of 2.91 and 0.48 kg/ha and average TN and TP reduction rates of 3.94 and 0.60 kg/ha, respectively.The relatively small reduction rates of TN and TP after full restoration (Scenario III) are associated with different N and P forms, for which part of dissolved N and P percolated from the wetlands would return to mainstreams with groundwater flow.As a result, the reduction rates for runoff, TN, and TP are on the same order and smaller than the sediment reduction rate based on model simulation (Table 6).In contrast, Scenario IV (loss of top 10 most effective wetlands for retention) and Scenario V (restoration of the 10 most effective wetlands) have average TN and TP increase rates of 8.0 and 1.2 kg/ha and average TN and TP reduction rates of 6.47 and 0.92 kg/ha, respectively.Considering the hundreds of existing and drained-lost wetlands in the study watershed, spatial targeting of those wetlands with higher TP reduction rates for retention and restoration has the potential of improving the effectiveness of wetland conservation programs.

Discussion and Conclusions
A cell-based modular modelling system, IMWEBs-Wetland, is developed for simulating and assessing the water quantity and water quality effects of individual wetlands at a watershed scale.The model is supported by one modular library and five databases (geospatial, hydro-climate, BMP, parameter, and output), which are managed by a Whitebox GAT based user interface.The model simulates processes of climate, flow, sediment, and water quality by incorporating land management practices.Compared to other watershed wetland models, the IMWEBs-Wetland model has distinguishing features of: (1) setting up the model based on project objective, watershed characteristics, data availability, and outputs with interest; (2) simulating hydrologic processes with more spatial details and providing both time series and spatial distribution outputs at user-defined spatial and temporal scales; (3) integrating with economic and ecologic models more easily for cost-effective assessment of wetland loss and restoration scenarios due to its cell based structure; and (4) interfacing with an open-source Whitebox GIS and SQLite database.However, some limitations also exist for the IMWEBs-Wetland model.For example, the model needs more detailed BMP distribution and operation data for model setup, and site-specific observation data for model calibration.Because the model runs for each grid cell, it would take a long computational time and a large memory space for a watershed with a small cell size.
A case study of IMWEBs-Wetland modelling was conducted in a small watershed in southern Manitoba of Canada with hundreds of existing and drained-lost wetlands.The model was setup based on existing wetlands and land management conditions.Calibration results demonstrated that the model performed well for flow, sediment, and water quality simulation at the watershed outlet.A simulated TP spatial distribution showed that TP concentrations were highly variable among existing wetlands in the study watershed depending wetland and its contribution area characteristics.However, because the model was not calibrated at wetland sites, this may cause uncertainties in the modelling results.Four wetland loss and restoration scenarios were constructed and evaluated with the calibrated IMWEBs-Wetland model.Compared to Scenario II (loss of all existing wetlands) and Scenario III (restoration of all drained-lost wetlands), the spatial targeting scenario IV (retention of the top 10 existing wetlands with the highest TP reduction efficiency) and scenario V (restoration of the top 10 drained-lost wetlands with the highest TP reduction efficiency) are more effective in reducing sediment, TN, and TP yields at the watershed outlet for wetland retention and restoration.As indicated by Fossey et al. [13] and Evenson et al. [15], field monitoring data are critical in validating distributed models and applying modelling results for watershed management.Because no monitoring data were available to validate the model at representative wetland sites in this study, aggregated calibration was conducted at the watershed outlet and uncertainties of modelling results would exist for individual wetlands.Therefore, site-specific observation data for enhanced model calibration are essential for improving the reliability of modelling results.
With a cell-based structure, modular system, open-source GIS and database, and the use of advanced modelling techniques, the IMWEBs-Wetland model is capable of modelling and assessing wetland and other BMP scenarios for spatial watershed management and decision-making.
The modelling results have the potential to improve the effectiveness of wetland conservation programs.

Water 2018 , 19 Figure 1 .
Figure 1.An open-water wetland and its drainage area delineated in the IMWEBs-Wetland model.

Figure 1 .
Figure 1.An open-water wetland and its drainage area delineated in the IMWEBs-Wetland model.

Figure 4 .
Figure 4.A schematic of wetland delineation with multiple outlets.

Figure 5 .
Figure 5. Location and land use of the study watershed.Figure 5. Location and land use of the study watershed.

Figure 5 .
Figure 5. Location and land use of the study watershed.Figure 5. Location and land use of the study watershed.

Figure 6 .
Figure 6.Existing wetlands in the study watershed.

Figure 6 .
Figure 6.Existing wetlands in the study watershed.

Figure 7 .
Figure 7. Observed and simulated flow at the EC9 monitoring station.

Figure 7 . 19 Figure 8 .
Figure 7. Observed and simulated flow at the EC9 monitoring station.Water 2018, 10, x FOR PEER REVIEW 14 of 19

Figure 8 .
Figure 8. Simulated spatial distribution of TP yield in the study watershed for the year 2010.

Figure 8 .
Simulated spatial distribution of TP yield in the study watershed for the year 2010.

Figure 9 .
Figure 9. Simulated spatial distribution of wetland TP concentration for the year 2010.

Figure 9 .
Figure 9. Simulated spatial distribution of wetland TP concentration for the year 2010.

Figure 10 .
Figure 10.Simulated top 10 wetlands for retention and top 10 wetlands for restoration.
Future research and development perspectives of the model include: (1) improvement of the model by developing new algorithms for specific hydrologic process modules; (2) improvement of the model by adding more BMP modules; (3) improvement of the model to perform an effective sensitivity

Table 1 .
Wetland parameters for each wetland.

Table 2 .
A summary of wetland type, number, area, and volume in the study watershed.

Table 2 .
A summary of wetland type, number, area, and volume in the study watershed.

Table 3 .
Crop management practices in the study watershed.

Table 4 .
Model performance at EC9 station for the 2009-2013 calibration period.

Table 4 .
Model performance at EC9 station for the 2009-2013 calibration period.

Table 5 .
Wetland loss and restoration scenarios for the study watershed.

Table 6 .
Runoff, sediment, TN and TP changes at the watershed outlet for different scenarios.