A coupled surface water storage and subsurface water dynamics model in SWAT for characterizing hydroperiod of geographically isolated wetlands

Advances in Water Resources 131 (2019) 103380

Contents lists available at ScienceDirect

Advances in Water Resources journal homepage: www.elsevier.com/locate/advwatres

A coupled surface water storage and subsurface water dynamics model in SWAT for characterizing hydroperiod of geographically isolated wetlands Junyu Qi a, Xuesong Zhang a,b,∗, Sangchul Lee c,d, Glenn E. Moglen d, Ali M. Sadeghi d, Gregory W. McCarty d a

Earth System Science Interdisciplinary Center, University of Maryland, College Park, MD 20740, USA Joint Global Change Research Institute, Pacific Northwest National Laboratory and University of Maryland, College Park, MD 20740, USA c Department of Environmental Science & Technology, University of Maryland, College Park, MD 20742, USA d USDA-ARS Hydrology and Remote Sensing Laboratory, Beltsville, MD 20705-2350, USA b

a r t i c l e Keywords: Richards equation Hydroperiod Inundation Soil moisture SWAT

i n f o

a b s t r a c t Wetlands play an important role in watershed hydrology and biogeochemistry. A geographically isolated wetland (GIW) module for Soil and Water Assessment Tool (SWAT) was developed to couple surface water storage and subsurface water dynamics to characterize hydroperiod. The new GIW module includes the following features: (1) a flexible geometric formula to characterize wetland surface water area, volume, and depth; (2) a revised algorithm to account for evaporation from both water and soil surfaces in wetlands; (3) application of the Richards equation to couple surface water storage and subsurface water dynamics; (4) use of Darcy’s law with an effective hydraulic conductivity parameter to simulate groundwater discharge. We tested the GIW module using observed daily water level data from four wetlands at two sites (including restored and natural wetlands with and without a low-permeability soil layer) in the Delmarva Peninsula, USA. The results show that the wetland module reasonably reproduced observed water levels for both restored and natural wetlands with and without a lowpermeability soil layer. The module was also able to reasonably simulate saturated and unsaturated portions of the soil corresponding to wet and dry periods. The ability of the GIW module to describe inundation conditions for wetlands holds promise to enhance the understanding and quantification of hydrological and biogeochemical roles of GIWs in a watershed context.

1. Introduction Wetlands, as a unique hydrologic feature of the landscape, generally form in flat areas, depressions, or shallow slopes in a watershed with water lying near the land surface (Keddy, 2010; Zedler and Kercher, 2005). Wetlands share many characteristics of aquatic and terrestrial environments due to their special transitional role between aquatic and terrestrial ecosystems (Keddy, 2010). As a result, wetland environments create unique hydrological and biogeochemical conditions that distinguish them from aquatic and terrestrial settings (Keddy, 2010). Wetlands support a variety of ecosystem functions including flood water storage, reducing delivery of sediment and nutrients to rivers, supporting biogeochemical cycling, and providing habitat for diverse plants and animals (Smith et al., 1995). For example, wetlands serve as a moderator of landscape water flow through storing flood water and reducing flow velocities (Dodds et al., 2008; Hefting et al., 2013; Ullah et al., 2005). In addition, wetlands are capable of trapping sediment and nutrients delivered from the surrounding upland areas due to their shallow

∗

depths and low slopes. Many studies have demonstrated that wetlands improved water quality by accumulating nutrients, trapping sediment, and transforming a variety of chemicals (Mitsch and Ewel, 1979; Jordan et al., 2003; Kuenzler, 1989; Faulkner and Richardson, 1989; Johnston, 1991). Wetlands also influence the global carbon (C) cycle and therefore climate change and global warming (Bridgham et al., 1995; Bridgham et al., 2006). Inundation conditions of wetland ecosystems limit the availability of oxygen to soil microbes and, as a result, organic matter decomposition proceeds slowly (Schlesinger, 2005). This condition further causes wetlands to accumulate organic matter in the substrate, such as the organic-rich soil of peatlands (Chanton et al., 1995). Wetlands have accumulated approximately one third of the global terrestrial soil C due to slow rates of anaerobic decomposition (Bridgham et al., 2006). Inundation conditions also result in the production of methane (CH4 ) by methanogenic microbes. It is estimated that wetlands are responsible for 15–40% of current global CH4 emissions (Solomon, 2007).

Corresponding author. E-mail addresses: [email protected], [email protected] (X. Zhang).

https://doi.org/10.1016/j.advwatres.2019.103380 Received 13 March 2019; Received in revised form 12 June 2019; Accepted 10 July 2019 Available online 10 July 2019 0309-1708/© 2019 Elsevier Ltd. All rights reserved.

J. Qi, X. Zhang and S. Lee et al.

The hydrological and biogeochemical functions of wetlands are determined by various factors such as climate, hydrology, underlaying substrate, as well as their position and area in the landscape (Keough et al., 1999; Zimmer et al., 2000). A geographically isolated wetland (GIW) is a freshwater wetland which is surrounded by uplands at the local scale (Tiner, 2003). About 8.3 million GIWs cover 6.5 million hectares and comprise approximately 16% of freshwater wetlands within the conterminous United States (Lane and D’Amico, 2016). GIWs generally lack permanent surface water connectivity to other aquatic systems, and as a result, they may have minimal influence on downstream flows because they are poorly connected with local river systems (Dahl, 2014). However, GIWs may have strong interactions with underlying groundwater systems depending on the hydraulic properties of the underlying substrate (McLaughlin et al., 2014). GIWs may be connected by intermittent surficial waters or drained via tiles and ditching for agricultural purposes or urban development (Lane and D’Amico, 2016; Rains et al., 2016). As a result, GIWs have been lost, to a large extent, due to human disturbance especially since the industrial age (Dahl et al., 1991; Rubec, 1994). Owing to the hydrological and ecological importance of wetlands, particularly under the influence of intensified of agricultural activities and climate change, quantitative tools are needed for assessing hydrological and biogeochemical effects of GIWs (Evenson et al., 2016). GIWs can be assessed through in-situ monitoring. However, instrumenting much larger watersheds having many GIWs is neither practical nor cost-effective. Hydrological and biogeochemical modeling is therefore an alternative approach for wetland functions assessment. Wetlands are directly or indirectly being included in many watershed scale models including the Soil and Water Assessment Tool (SWAT) model (Arnold et al., 1998). The original wetland module within the SWAT model can reproduce the streamflow hydrology of watersheds having a considerable portion of wetlands (Wu and Johnston, 2008; Wang et al., 2008; Martinez-Martinez et al., 2014; Javaheri and BabbarSebens, 2014; Babbar-Sebens et al., 2013). To simulate the hydraulic interactions between riparian wetlands, rivers, and aquifers, a previous study modified the SWAT model with an enhanced riparian and depressional wetland module (Rahman et al., 2016). Remotely sensed data were also employed in an enhanced SWAT model to simulate the spatial dynamics of inundation extent of riparian wetlands at the catchment scale (Lee et al., 2017). The SWAT model was also modified to examine the hydrological characteristics of GIWs (Evenson et al., 2016; Lee et al., 2018). However, none of these studies explicitly investigated the hydroperiod of wetlands with respect to the interaction between subsurface soil water and surface water storage, which is particularly important for modeling hydrological and biogeochemical processes in wetlands. For example, the water table and its fluctuation within wetlands is the primary factor driving soil organic C (SOC) decomposition, plant C fixation, CH4 production and consumption, and other biogeochemical processes (Zhang et al., 2002). An accurate description of hydroperiods regarding soil water and surface water interaction benefits water quantity and quality assessment from site to watershed scales. In the present study, we propose a GIW module with enhanced capabilities describing the interaction between surface water storage and subsurface water dynamics so that the hydroperiod can be explicitly modeled. Specifically, the objective of this study is to develop and test a GIW module in the SWAT model using monitored water level data from restored and natural wetlands with and without low-permeability soil layers in the coastal plain of the Chesapeake Bay. 2. Materials and methods 2.1. Soil and Water Assessment Tool The SWAT model is a continuous, physically-based, and semidistributed watershed scale model which has been successfully used to simulate water quantity and quality in a wide range of watersheds

Advances in Water Resources 131 (2019) 103380

across the world (Qi et al., 2017a, 2017b; Srinivasan et al., 2010; Zhang et al., 2008; Abbaspour et al., 2017). This model has been widely used to evaluate water quality and assess effectiveness of BMPs as affected by landuse and climate change (Li et al., 2014; Qi et al., 2018a; Zhang et al., 2017; Qi et al., 2019a, 2019b; Liang et al., 2019). It is currently one of the most widely used hydrological models for water resource assessment and watershed management. Major SWAT model components include weather, hydrology, soil temperature, sedimentation, nutrients, pesticides, pathogens, plant growth, and land management (Neitsch et al., 2011a). Model inputs include the physical characteristics of a watershed defined by local meteorology, soil type, and topography. Either simulated or recorded weather data (e.g., precipitation, temperature, and solar radiation) may be used by SWAT (Qi et al., 2019c). The model analyzes small or large watersheds by discretizing them into subbasins, which are then further subdivided into hydrological response units (HRUs) with homogeneous land use and soil properties and slope. The model calculates the water balance (i.e., surface and subsurface runoff, percolation and base flow, and evapotranspiration (ET) and transmission losses), crop growth, nutrient cycling, and pesticide movement at the HRU scale. Water flow, sediment, and nutrient loadings from each HRU in a subbasin are summed and the resulting loadings are then routed through channels, ponds, and reservoirs to the watershed outlet. Water balance within the soil is a function of precipitation, surface runoff, lateral flow, evapotranspiration, and percolation (Neitsch et al., 2011a). Surface runoff is calculated using the Soil Conservation Service (SCS) Curve Number (CN) method (Neitsch et al., 2011a). Water infiltrated into the soil profile is simulated using a storage routing technique (Narasimhan et al., 2005). Downward flow in the soil occurs when a soil layer exceeds its field capacity and the layer below is not saturated. Soil water can evaporate, be taken up and transpired by plants, and flow into water bodies through lateral flow, or percolate into groundwater. The water entering the vadose zone flows into shallow aquifer becoming unconfined groundwater. Subsequent flow may occur as baseflow to the main stream, as confined groundwater flow in a deep aquifer, or may evaporate if there is insufficient soil water to support ET demand (Neitsch et al., 2011a). Three ET methods are included in the SWAT model, i.e., Penman-Monteith method (Monteith, 1965), Priestley-Taylor method (Priestley and Taylor, 1972), and Hargreaves method (Hargreaves and Samani, 1985). These methods vary in the amount of required inputs and conditions they can be used. The Priestley-Taylor method is suitable for conditions where surface areas are wet, and as a result, we employed the Priestley-Taylor method for GIW modeling. 2.2. GIW module The GIW module was developed at the HRU scale because of finer scale resolution better aligned with our objective of simulating surface water storage and soil water dynamics for GIWs. We incorporated a robust, generalized, and flexible wetland geometric formula to characterize wetland surface water area, volume, and depth. We developed a wetland evaporation algorithm to account for evaporation from both water and soil surfaces of GIWs through variable water / soil surface areas. We improved the calculation of seepage from the wetland bottom by redefining the soil hydraulic conductivity. Most importantly, we used the Richards equation to characterize soil moisture movement in the soil profile. These modifications dramatically improved soil water simulation and facilitated the coupling of surface water storage and subsurface water dynamics (Qi et al., 2018b). Finally, we also modified the groundwater algorithms corresponding to two different types of GIWs, i.e., GIWs with a low-permeability soil layer and GIW without a lowpermeability soil layer. Our efforts made for a better representation of groundwater flow by applying Darcy’s law with an effective conductivity. The newly developed GIW module simulates conditions at an hourly time step and is illustrated in detail as follows.

J. Qi, X. Zhang and S. Lee et al.

Advances in Water Resources 131 (2019) 103380

2.2.3. Surface runoff Surface runoff is generated when surface water storage SWS (mm) is greater than the maximum surface water storage SWSmx (mm) of a GIW, calculated as, 𝑅𝑜𝑢𝑡 = 𝑆 𝑊 𝑆 − 𝑆 𝑊 𝑆𝑚𝑥

(6)

where SWSmx is determined by, 𝑆 𝑊 𝑆𝑚𝑥 = 1000 ⋅

𝑆 𝑉𝑚𝑥 𝐺𝐼𝑊 𝐴

(7)

where SVmx (m3 ) is the maximum surface water volume, and is calculated based on Eqs. (1) and (2) as, 𝑆 𝑉𝑚𝑥 =

Fig. 1. Conceptualized GIW shape with plan view (upper circles) and profile view. Black lines and letters indicate GIW maximum surface area GIWA, radius rmx , and depth Dmx . (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

2.2.1. GIW spatial structure A GIW is conceptualized as having a circular boundary as shown in Fig. 1. We incorporated a generalized and flexible geometric formula to determine the relationships between surface water volume SV (m3 ), surface area SA (m2 ), and depth D (i.e., water depth at the center of GIWs; m) for GIWs. This formula was developed by Hayashi and van der Kamp (2000) who tested it for a range of GIWs with non-unique shapes. The mathematical form of the model is shown as, )2∕𝑝 ( 𝐷 𝑆𝐴 = 𝑏 ⋅ (1) 𝐷0 ( ) 𝑏 𝐷1+2∕𝑃 𝑆𝑉 = ⋅ (2) 2∕𝑃 1 + 2∕𝑝 𝐷0 where b and p are scale and shape factors, and D0 is unit depth (i.e., 1 m). The scale factor b can be calculated using the user specified maximum values of GIW water surface area GIWA (m2 ) and depth Dmx (m), given by, 𝑏=

𝜋 ⋅ 𝑟2𝑚𝑥 (𝐷𝑚𝑥 ∕𝐷0 )2∕𝑝

(3)

where rmx is the maximum radius for GIWs (Fig. 1) and it is calculated by, 𝐺𝐼𝑊 𝐴 1∕2 (4) ) 𝜋 The shape factor, p, is used to calibrate wetland shape. For example, increasing the value of p represents a more cylindrically-shaped wetland. Eqs. (1) and (2) can be used to calculate GIW surface water area SA (m2 ) and depth D (m) given a value of GIW surface water volume SV (m3 ), which is converted from surface water storage SWS (mm) simulated by the GIW module at the HRU scale. 𝑟𝑚𝑥 = (

2.2.2. Surface water storage The HRU water balance for surface water storage SWS (mm) is calculated by, Δ𝑆𝑊 𝑆 = 𝑃 − 𝐸 − 𝑅𝑜𝑢𝑡 − 𝑆

(5)

where ΔSWS is the change in surface water storage (mm), P is precipitation (including rainfall and snowfall; mm), E is evaporation (mm), Rout is surface runoff (mm), and S is seepage (mm) into underlying soils.

𝑏 1+2∕𝑝 2∕𝑝 ⋅ 𝐷𝑚𝑥 ∕𝐷0 (1 + 2∕𝑝)

(8)

2.2.4. Evaporation The GIW module takes evaporation from both water surface and soil surface into account, calculated by, ( ( ) ) 𝑆𝐴 𝑆𝐴 𝐸 = 𝐸𝑤 ⋅ + 𝐸𝑠 ⋅ 1 − (9) 𝐺𝐼𝑊 𝐴 𝐺𝐼𝑊 𝐴 where Ew and Es are the evaporation from the water and soil surfaces (mm), respectively. They are determined as a function of plant leaf area index, LAI, and the maximum leaf area index, LAImax , as follows, ( ) 𝐿𝐴𝐼 𝐸𝑤 = 𝑃 𝐸𝑇 ⋅ 1 − (10) 𝐿𝐴𝐼𝑚𝑎𝑥 ( ) 𝐿𝐴𝐼 𝐸𝑠 = 𝐸𝑠𝑜 ⋅ 1 − (11) 𝐿𝐴𝐼𝑚𝑎𝑥 where PET and Es o are the potential evapotranspiration (mm) and soil evaporation (mm), respectively. Both are determined by the PriestleyTaylor method within the SWAT model (Neitsch et al., 2011a). 2.2.5. Seepage Initial seepage, Sini , (mm) from the bottom of a GIW is calculated by, ( ) 𝑘 ⋅ 𝑆𝐴 𝑆𝑖𝑛𝑖 = min 𝑆 𝑊 𝑆 , 1 𝐺𝐼𝑊 𝐴

(12)

where k1 is the hydraulic conductivity (mm) for the first soil layer. Maximum soil pore space, Scap , (mm) is determined as, 𝑆𝑐𝑎𝑝 =

𝑘𝑘 ∑

𝑃 𝑜𝑟_𝑒𝑓 𝑓𝑖

(13)

1

where Por_eff is the effective porosity (air-filled space of soil layers; mm), i indicates soil layers from top (i = 1) to bottom (i = last soil layer number), and kk is determined by the first soil layer whose saturated hydraulic conductivity is less than Sini or the first saturated soil layer when searching from the top soil layer to the bottom soil layer. The final seepage S (mm) is then calculated by, ( ) 𝑆 = 𝑚𝑖𝑛 𝑆𝑖𝑛𝑖 , 𝑆𝑐𝑎𝑝 + 𝑘𝑚𝑖𝑛 (14) where kmin (mm h−1 ) is the minimum saturated hydraulic conductivity of soil layers. 2.2.6. Soil water simulation Soil water is simulated by the modified Richards equation (Zeng and Decker, 2009), i.e., [ ( ( ) )] 𝜕 ℎ − ℎ𝑒 𝜕𝜃 𝜕 = 𝑘⋅ −𝑄 (15) 𝜕𝑡 𝜕𝑧 𝜕z where 𝜃 is the volumetric soil water content (mm3 mm−3 ), t is time (h), z is the depth below soil surface (mm; positive downwards), k is the hydraulic conductivity (mm h−1 ), h is the soil matric potential (mm), Q is a soil water sink term (mm3 mm−3 h−1 ), and he is the equilibrium soil matric potential (mm). Detailed information regarding the development and integration of the Richards equation-based soil moisture module in

J. Qi, X. Zhang and S. Lee et al.

Advances in Water Resources 131 (2019) 103380

Fig. 2. Water table height wt relative to the main stream (reference elevation) and the distance from the wetland to the main channel L in a subbasin for wet and dry periods.

the SWAT model can be found in (Qi et al., 2018b). In this study, we further modified the Richards equation-based soil moisture module for GIW simulation. The modifications are: (1) The upper boundary flux for the Richards equation is the seepage Eqs. (12)–((14)) from GIWs when surface water exists. When there is no surface water storage, an evaporation flux occurs based on Eq. (9) (i.e., Es ). (2) Plant transpiration is the only sink term in the Richards equationbased soil moisture module for wetlands with a low-permeability layer, while for wetlands without a low-permeability layer, an extra sink term is added based on the predetermined groundwater flow (see Eq. (16)). (3) Saturated hydraulic conductivity is calculated based on the method used by the Community Land Surface (CLM) model (Oleson et al., 2010). This is because the measured data from the SWAT model database i.e., USDA Natural Resources Conservation Service (NRCS) Soil Survey Geographic Database (SSURGO) may not be accurate for soils in GIWs. (4) The lower boundary condition is modified to simulate two types of GIWs based on soil impermeability characteristics. When there is a low-permeability soil layer, flux is limited via calibrating the saturated hydraulic conductivity for the low-permeability layer Wet_K. The layer with highest content of clay is defined as the lowpermeability layer. When there is no low-permeability soil layer (hydrologically connected soil water and groundwater system), the shallow aquifer was included in the model to retain hydraulic interaction between soil water and the unconfined groundwater. Specifically, a numerical node below the soil bottom was added in the solution of the Richards equation to calculate recharge to the shallow aquifer or the upward movement of groundwater due to soil evaporation and uptake by deep-rooted plants (Oleson et al., 2010). The hydraulic properties for the additional node were adapted from those of the last soil layer (excluding organic material properties) to estimate hydraulic conductivity and soil matric potential (Qi et al., 2018b). The lower boundary condition is set at the interface layer of shallow and deep aquifers with a low flux which needs to be calibrated via saturated hydraulic conductivity, Wet_K. 2.2.7. Groundwater flow Since there are low-permeability layers, a relatively weak hydraulic connection exists between soil water and the shallow aquifer. The recharge to the shallow aquifer is determined by the limited flux at the low-permeability layer in the Richards equation-based soil moisture module. In this case, the groundwater discharge Gw (i.e., baseflow in mm h−1 ;) is calculated using the original algorithm in the SWAT model (Neitsch et al., 2011b). To simulate groundwater discharge, Gw, (mm h−1 ) for GIWs without a low-permeability layer, Darcy’s law was used, given by, 𝐺𝑤 = 𝑘𝑒𝑓 𝑓 ⋅ 𝐴𝑐 ⋅

𝑤𝑡 𝐿

(16)

Where Ac is the cross-section area (m2 ) of water flow between the wetland and the main channel in a subbasin, wt is the water table height

(m) relative to the reference elevation (Fig. 2), L (m) is the distance from the wetland to the main channel in a subbasin (Fig. 2), and keff is the effective hydraulic conductivity (mm h−1 ) which needs to be calibrated. Since we did not know the exact cross section area between the wetland and the main channel in a subbasin, we assumed that the Ac = GIWA, similar to the concept using in (Rahman et al., 2016). By this assumption, the water balance and groundwater dynamics can be conveniently simulated at the HRU scale (with all water fluxes in unit of mm h−1 and volume in unit of mm). Note that this assumption has impacts on the value of keff and we will discuss this point in Section 4.2. The reference elevation is defined at the main stream Elestr (m). The water table height above the bottom of a GIW can be determined by, 𝑤𝑡 = 𝐸 𝑙𝑒𝑏𝑜𝑡 − 𝐸 𝑙𝑒𝑠𝑡𝑟 + 𝐷

(17)

where, Elebot is the invert elevation of GIWs (lowest land surface elevation on the wetland bottom; m). When D = 0, the water table height is at the bottom of a GIW. The water table height below the bottom of a GIW can be determined by, 𝑤𝑡 = 𝐸 𝑙𝑒𝑏𝑜𝑡 − 𝐸 𝑙𝑒𝑠𝑡𝑟 − 𝐷𝑠𝑜𝑙

(18)

where Dsol (m) is the water level depth (positive) when the water table is in the soil (typically during dry periods). Elestr , Elebot , and L can be derived from the digital elevation model (DEM; see more information in Section 2.4). Note that in Eq. (16), the model can simulate recharge from adjacent stream to wetlands when wt is negative. Because our observations do not include conditions where wetland water table is below stream elevation, this function remains to be validated in future studies. The water stored in the unconfined aquifer has a prescribed maximum value (5000 mm). When the water table is below the soil bottom, the change in water storage in the unconfined aquifer is updated as the difference between recharge, Gr, (mm h−1 ; determined from the Richards equation), and discharge, Gw, (mm h−1 ). The change in water table is then calculated by, 𝐺 − 𝐺𝑤 Δ𝑤𝑡 = 𝑟 , (19) 𝑆𝑦 where Sy is the specific yield for the unconfined aquifer (m3 m−3 ) (Masterson et al., 2016). 2.2.8. GIW module parameterization The GIW module has several physically-based calibration parameters. The shape factor p Eqs. (1) and ((2)) is used to calibrate the GIW shape, which is important for simulation of water level at the center of a GIW. The saturated hydraulic conductivity, Wet_K, is used to calibrate limited water flux for the low-permeability layer for both types of GIWs. The effective hydraulic conductivity, keff , is used to calibrate groundwater discharge for GIWs without a low-permeability layer (Eq. (16)). The GIW module also needs several measured input parameters including maximum wetland water surface area, GIWA, and maximum depth at the center of the wetland, Dmx . In addition, the reference elevation at the main stream, Elestr , the elevation of wetland bottom, Elebot , and the distance from the wetland to the main channel, L, as well as the specific yield, Sy , for the unconfined aquifer are needed for GIWs without a low-permeability layer.

J. Qi, X. Zhang and S. Lee et al.

Advances in Water Resources 131 (2019) 103380

Fig. 3. The location of Site #1 and # 2 in the Chesapeake Bay Watershed as well as images of natural and restored wetlands in the two sites.

2.3. Study sites and data collection The present study involves a watershed defined by a U.S. Geological Survey (USGS) gauge station on the Choptank River near Greensboro (USGS#01491000), Maryland, which hereafter will be referred to as the Greensboro Watershed (GW; 290 km2 ; Fig. 3). It is located on the headwaters of the Choptank River Watershed (CRW) in the coastal plain of the Chesapeake Bay in Maryland (Fig. 3). Due to flat topography, closeness to the groundwater table, and high annual precipitation, wetlands are abundant in this region. Those wetlands are reported to provide various ecosystem services, water quality improvement (Jordan et al., 2003), water quantity control (Lee et al., 2018), and provision of natural habitat (Russell and Beauchamp, 2017). In the GW, major land uses are forest (48.3%) and agriculture (36.1%). A large portion of soils (74.5%) in the GW are poorly-drained (hydrological soil groups C and D) in which 67.2% of croplands are located (Ator and Denver, 2015). The SWAT model requires detailed information on the climate, soils, and land use for the study site. Precipitation and temperature data were collected from National Climatic Data Center (NCDC) at the Greensboro station (US1MDCL0009) within the GW (Fig. 3). Other climatic variables, such as solar radiation, relative humidity, and wind speed, were generated by the weather generator (WXGEN) embedded within the SWAT model due to data unavailability. A soil map was derived based on the SSURGO database. Topography was delineated by a 10 m Light Detection and Ranging (LiDAR)-based DEM. The land use map and the scheduling of crop rotations were generated using 2008–2012 data from the USDA-National Agriculture Statistics Service (NASS) Cropland Data Layer (CDL) (Lee et al., 2016). In 2016, a well was installed at the wetland invert to monitor the water level for restored and natural wetlands at Sites #1 and #2 as shown in Fig. 3 (Lee et al., 2019). The well was installed with a depth of 0.9 m below the wetland invert elevation. The water level, when it is below the wetland invert elevation, is defined as the depth of the interface

between unsaturated and saturated soil layers (these will be negative values as the wetland invert elevation is the reference datum). According to hydrogeologic data, Site #2 is characterized by soils with large saturated hydraulic conductivity, while Site #1 includes multiple lowpermeability layers (Table 1) (Lee et al., 2019). Daily water level data were observed from 2016 to 2017 at Site #2 and for 2016 at Site #1 (Table 2). Hourly precipitation data measured by a tipping-bucket rain gauge (Campbell Scientific TE525-L10) at each wetland were used for investigating the variation of water level over time (Lee et al., 2019). Daily precipitation was calculated as the sum of hourly precipitation for each day and used to simulate water level dynamics. Restored wetlands are those that have been converted from active agricultural use and have had a more natural wetland hydroperiod restored through a variety of methods, including plugging ditches with earthen mounds. Plant communities in restored sites resembled natural wetlands based on the percentage of species that were native and hydrophytic. However, natural sites were forested, while restored wetlands were primarily herbaceous. Restored wetlands tend to be more affected by anthropogenic disturbance compared to natural sites. They can develop diverse native wetland plant communities within a decade, but they remain different from natural wetlands (Yepsen et al., 2014).

2.4. Model evaluation Site # 2 is located within HRU # 170016 in the GW (Fig. 3). The original HRU is characterized as: FRSD (forest landuse), MD011HoB1 (soil type), and 0–2 (slope class) within the SWAT model. For the natural GIW, we changed FRSD to WETF (Wetlands-Forested), while for the restored GIW we used RNGE (Range-Grasses). The same HRU was used for GIWs at Site # 1, but with a different soil type (Table 1). The GIW module was used to simulate daily water level dynamics in the HRU. The water level simulated is the water depth (positive when SWS is

J. Qi, X. Zhang and S. Lee et al.

Advances in Water Resources 131 (2019) 103380

Table 1 Physical properties for soils at Site # 1 and 2 based on USDA-NRCS. Site

Site #1

Site #2

Depth (mm)

Org. Matter (%)

Sand (%)

Silt (%)

Clay (%)

Bulk Density (g cm−3 )

Satu. Conductivity (mm h−1 )

203 305 610 1397 1575 2032 100 360 760 1000 1500 2030

2.50 3.25 0.25 0.25 0.25 0.25 0.70 0.15 0.23 0.23 0.23 0.12

17 13 5 6 24 47 80 82 65 65 65 88

71 73 62 64 52 35 14 14 19 19 19 7

12 14 33 29 24 18 6 4 16 16 16 5

1.45 1.45 1.65 1.65 1.65 1.70 1.55 1.40 1.55 1.55 1.55 1.60

32 32 3 3 3 256 331 331 100 100 100 1520

Table 2 GIW input parameters and water level data availability for restored and natural GIWs at two sites. Site

Wetland

Dmx (m)

GIWA (m2 )

Elestr (m)

Elebot (m)

L (m)

Sy (m3 m−3 )

Data Availability

Site #1

Natural Restored Natural Restored

0.37 0.4 0.65 0.43

1552 2536 3570 3000

– – 5.7 5.7

– – 15.2 15.5

– – 1842 1842

– – 0.2 0.2

2016

Site #2

2016–2017

greater than zero and negative when the soil is unsaturated) referenced at the wetland invert elevation. The input parameters required, i.e., maximum wetland surface water area, GIWA, and depth, Dmx , as well as reference elevation at the main stream, Elestr , the elevation of GIW bottom (the wetland invert elevation), Elebot , and the distance from the wetland to the main channel, L, as well as the specific yield, Sy , for the unconfined aquifer are shown in Table 2. We choose the elevation at the outlet of the subbasin #17 (containing HRU # 170016) to represent the reference elevation of the main stream (Fig. 4), Elestr . Thereafter, the distance from the wetland to the outlet of the subbasin, L, was determined using ArcGIS for Site # 2 (Fig. 4). Elestr and Elebot were determined for the restored and natural GIWs at Site # 2 via ArcGIS using DEM data. GIWA and Dmx were estimated based on our observations at the two sites. The value for the specific yield, Sy , for the unconfined aquifer was adopted from (Masterson et al., 2016). The GIW module was initialized with saturated soils and maximum surface water storage because the model started simulation on the first day of the year when wetlands were filled with full volume of water at the study area. We simulated conditions from 2014 to 2016/2017 with 2014–2015 (two years) serving as the warm up (spin-up) period for natural and restored wetlands at the two sites. The weather data from NCDC were used for spin-up period. The GIW module was calibrated using the Sequential Uncertainty Fitting algorithm version 2 (SUFI-2) procedure in the SWAT Calibration and Uncertainty Program (SWATCUP) (Abbaspour et al., 2007). Specifically, the GIW module was executed 1000 times for each iteration. We found that no distinct improvement in objective function value (i.e., the Nash-Sutcliffe coefficient (NS) (Nash and Sutcliffe, 1970)) after 2–3 iterations. Calibration was conducted using data from 2016 for both the natural and restored wetlands at the two sites. Model validation was conducted only for natural and restored wetlands at Site # 2 using data from 2017 (Table 2). We do not have additional data for model validation at Site # 1 (Table 2). Performance of the GIW module was assessed using three coefficients of accuracy, i.e., Bias, NS, coefficient of determination (R2 ), given as: 𝐵𝑖𝑎𝑠 = 𝑃𝑎𝑣𝑔 − 𝑂𝑎𝑣𝑔 , ∑𝑛 (

𝑁𝑆 = 1 − ∑ 𝑖=1( 𝑛 𝑖=1

𝑂𝑖 − 𝑃𝑖

(20) )2

𝑂𝑖 − 𝑂𝑎𝑣𝑔

)2 ,

(21)

Fig. 4. The location of Site # 2 and associated subbasin # 17 and the outlet in the Greensboro Watershed. The elevation at the outlet of the subbasin #17 is used to represent the reference elevation of the main stream, Elestr ; the distance from Site # 2 to the outlet of the subbasin # 17, L, was also determined using ArcGIS (Table 2).

J. Qi, X. Zhang and S. Lee et al.

Advances in Water Resources 131 (2019) 103380 2

3.2. Model calibration and validation

) ( ) ⎛ ⎞ ∑𝑛 ( ⎜ ⎟ 𝑖=1 𝑂𝑖 − 𝑂𝑎𝑣𝑔 ⋅ 𝑃𝑖 − 𝑃𝑎𝑣𝑔 𝑅2 = ⎜ [ , ( )2 ∑𝑛 ( )2 ]0.5 ⎟ ∑ 𝑛 ⎜ ⎟ − 𝑂 − 𝑃 𝑂 𝑃 𝑖 𝑎𝑣𝑔 𝑖 𝑎𝑣𝑔 𝑖=1 𝑖=1 ⎝ ⎠

(22)

where Oi and Pi are the observed and predicted values, respectively, and Oavg and Pavg are the average of the observed and predicted values, respectively. 2.5. Sensitivity analysis Sensitivity analyses were conducted on the four calibration parameters. Parameter sensitivities were determined using the following multiple regression equation, based on results running SWAT-CUP multiple times, given as: ∑𝑚 g=α+ 𝛽𝑖 ⋅ 𝑏 𝑖 (23) 𝑖=1

where g is the objective function value, 𝛼 and 𝛽 i are regression coefficients, bi is the calibration parameters, and m is the number of parameters considered. NS was used as the objective function value and student t-tests were used to identify the statistical significance of each parameter. This sensitivity analysis method estimates the change in the objective function resulting from changes in each parameter while all other parameters are changing. The procedure does not provide an absolute measure of the sensitivity but rather a relative sensitivity. This method has the advantage of being quite fast compared to similar procedures. The p-value was used to identify the relative sensitivity of parameters and a p-value less than 0.05 was taken to represent a relatively sensitive parameter in the present study.

Model performance evaluation statistics for simulating daily water level during calibration and validation are shown in Table 5. For Site # 1, the GIW module performed well for water level simulation for both the natural and restored wetlands based on the three statistics during calibration. Values of R2 and NS were all greater than 0.68 and Bias values were less than 11 mm (which is relatively small compared to the average annual precipitation about 1200 mm). The GIW module performed slightly better for the restored wetland than for the natural wetland with higher values of NS and lower Bias (Table 5). For Site # 2, the GIW module performed much better for both the natural and restored wetlands during calibration (values of R2 and NS were all greater than 0.83) compared with the model performance at Site # 1 based on R2 and NS (Table 5). The value of Bias was −22 mm (approximately 2% of average annual precipitation) for the restored wetland indicating an underestimation of water level during calibration. During validation, values of R2 and NS were greater than 0.76 suggesting that the GIW module performed well for both the natural and restored wetlands (Table 5). The module overestimated water level with a Bias value of 30 mm (approximately 3% of average annual precipitation) in the restored wetland. Observed vs. simulated daily water level during calibration and validation for natural and restored wetlands at both sites are shown in Fig. 5. In general, the simulated water level curves matched observed curves well for both the natural and restored wetlands at the two sites during both calibration and validation. This is consistent with the results presented in Table 5 (as indicated by large R2 and NS values).

3. Results 4. Discussion

3.1. Parameter calibration and sensitivity analysis Calibrated parameter values for the GIW module at the two sites are shown in Table 3. In general, the calibrated shape factor, p, was within a narrow range i.e., from 1.09 to 1.75 for all four wetlands. Calibrated values for saturated hydraulic conductivity, Wet_K, in the low-permeability layer for Site # 1 and the interface layer between the shallow and deep aquifers for Site # 2 were very close to each other (ranging from 0.030 to 0.036 mm h−1 ). The values of effective hydraulic conductivity, keff , were rather close for the natural and restored wetlands (0.875 and 3.3 mm h−1 , respectively) at Site # 2 (Table 3). Sensitivity analysis results are shown in Table 4. The p-values for shape factor p was greater than 0.05 while those for Wet_K and keff were both less than 0.05, indicating that hydraulic properties in the low-permeability layer and in the lateral direction (from wetlands to streams) are important parameters influencing water level dynamics in both natural and restored wetlands. Table 3 Calibrated GIW module parameters for the four wetlands.

Site # 1 Site # 2

Wetland

Shape Factor p

Wet_K (mm h−1 )

keff (mm h−1 )

Natural Restored Natural Restored

1.75 1.48 1.09 1.67

0.035 0.035 0.030 0.036

– – 0.875 3.300

Table 4 The p-value for parameter sensitivity four the four wetlands.

Site # 1 Site # 2

Wetland

Shape factor p

Wet_K

keff

Natural Restored Natural Restored

0.581 0.461 0.088 0.428

<0.05 <0.05 <0.05 <0.05

– – <0.05 <0.05

4.1. Overall model performance evaluation The usage of R2 and NS for model performance evaluation places more emphasis on larger errors while smaller errors tend to be neglected (error is the difference between the simulated and the observed variable at each time step) (Krause et al., 2005). Thus, R2 and NS are more sensitive to high flow (i.e., peak flow) rather than low flow (i.e., baseflow). This deficiency may be undesirable for studies focusing in detecting low flow, however, it becomes an advantage for detecting the largest values on the wetland water level fluctuation curve. Here, the largest departures from zero on the water level curves are created mainly during dry periods when the surface water dries up and the top soil layers become unsaturated. During these periods, large negative values on the water level curve occurs (for example, see water levels in Sep.–Oct. of 2016 in Fig. 5a–d). The GIW module developed here was intended to accurately describe these sensitive periods when surface water disappears and subsurface water declines. In fact, these periods reflect key processes influencing biogeochemical reactions in soils. For example, accurate simulations of the beginning of surface water disappearance and surface water accumulation after dry periods, as well as unsaturated soil depth during dry periods, contribute to the accurate prediction of CH4 generation, transportation, and emission from wetlands. Based on model performance evaluation as indicated by R2 and NS during calibration and validation at the two sites, we found that the GIW model can produce well the wetland water level curves for both types of wetlands, with and without a low-permeability soil layer. Evaluation statistics indicated that the GIW module had large Bias for the restored wetlands at Site # 2 (Table 5). This is probably because the hydraulic conductivity of the restored wetland is affected by the restoration process including excavations. Soils within restored wetlands have higher soil bulk density than natural wetlands due to increased soil resistance resulting from construction (Campbell et al., 2002; Bantilan-Smith et al., 2009). This anthropogenic impact on the restored wetlands might not be fully reflected in the soil data used in this study.

J. Qi, X. Zhang and S. Lee et al.

Advances in Water Resources 131 (2019) 103380

Table 5 GIW module performance evaluation during calibration and validation. Calibration

Site # 1 Site # 2

2

Wetland

R

Natural Restored Natural Restored

0.71 0.72 0.90 0.88

Validation

NS

Bias (mm)

Mean Sim. (mm)

Mean Obs. (mm)

R2

NS

Bias (mm)

Mean Sim. (mm)

Mean Obs. (mm)

0.68 0.72 0.90 0.83

11 6 5 −22

207 226 481 171

196 220 476 193

– – 0.84 0.84

– – 0.80 0.76

– – −10 30

– – 583 52

– – 593 22

Fig. 5. Observed vs. simulated daily water level for wetlands at Site # 1 and 2: (a) natural wetland at Site # 1 during calibration; (b) restored wetland at Site # 1 during calibration; (c) natural wetland at Site # 2 during calibration and validation; (d) restored wetland at Site # 2 during calibration and validation.

Fig. 6. shows the simulated evapotranspiration and seepage by the GIW module for the four wetlands. Daily seepage varied between 0.0 mm day−1 (when surface water disappears) and ca. 1.6 mm day−1 (during winter) for natural wetland at Site #1 and restored wetlands at both Site #1 and 2 Fig. 6a, b, and d, respectively). Meanwhile, the maximum seepage value for the natural wetland at Site # 2 was ca. 0.7 1.6 mm day−1 which was lower than those for the other three wetlands. The results are consistent with calibrated Wet_K values in Table 3 where the natural wetland at Site # 2 had the lowest value while other three wetlands had higher values reflecting that the seepage is not only controlled by surface area (daily seepage varied consistently with water level; Figs. 5 vs. 6), but also determined by the minimum hydraulic con-

ductivity of fluid medium (Eqs. (12)–((14)). The two main forces that drive surface water dynamics are rainfall and evaporation. Rainfall is evenly distributed throughout the year in the study area with relatively more dry days in late summer and early fall, which caused surface water to disappear during this time (Fig. 5) instead of during the peak of evapotranspiration (Fig. 6). 4.2. Parameter estimation The most sensitive parameter is the saturated hydraulic conductivity, Wet_K, for the low-permeability soil layer (in the wetlands at Site # 1) and the interface layer between the shallow aquifer and the deep

J. Qi, X. Zhang and S. Lee et al.

Advances in Water Resources 131 (2019) 103380

Fig. 6. Simulated daily evapotranspiration (ET) and seepage by the GIW module for the four wetlands from during calibration and validation.

aquifer (in the wetlands at Site # 2). The calibrated values for Wet_K for the low-permeability soil layer at Site # 1 are about two orders of magnitude lower than reported values by SSURGO database (about 0.03 vs. 3 mm h−1 ; Tables 1 and 3). Our calibrated values are within the range of those reported in peatland soil studies (Holden and Burt, 2003; Rycroft et al., 1975; Ingram, 1967). It has been reported that peatland soils share similar hydraulic properties with wetland soils (Holden et al., 2004). These results indicate that SSURGO data may not be accurate for wetlands. Calibrated values of saturated hydraulic conductivity at the interface layer between the shallow and the deep aquifer at Site # 2 are of similar magnitude to reported values of low-permeability aquifers (Bear, 2012), and represent reasonable estimates. The GIW module uses a simple groundwater algorithm to simulate groundwater discharge for wetlands without the low-permeability soil layer (e.g., Site # 2). The effective hydraulic conductivity, keff , needs to be calibrated if no measurement is available. The calibrated effective hydraulic conductivity ranged from 0.875 to 3.3 mm h−1 for natural and restored wetlands in Site # 2. These calibrated values are less than the range of commonly measured saturated hydraulic conductivity of soils (Klute and Dirksen, 1986) but greater than those measured from peatlands (Holden and Burt, 2003). The heterogeneity of soil properties and geological features along the hydraulic gradient from the wetland to the outlet of a subbasin may explain the somewhat small calibrated keff values determined in the present study. In addition, we assumed that the cross-section area (Ac ) of groundwater flow between the wetland and the main channel was equal to the maximum surface area GIWA, which would impact on the calibrated value of keff according to Eq. (16). Given all those uncertainties, keff is used as a fitting parameter in the present study. 4.3. Model performance at site # 1 vs. Site # 2 Model performance at Site # 2 was better than that at Site # 1 during calibration (Table 5). A major reason is that the weather data (i.e., solar radiation, air temperature, and relative humidity) used to calculate PET

with the Priestley-Taylor method were derived from the weather station within the GW for both sites. The natural and restored wetlands at Site # 1 are located outside the GW so that larger uncertainty in weather data can be expected. This can also explain the discrepancies between simulated and observed water level after Oct. 2016 (Fig. 5a and b). This result demonstrates the importance of accurate weather input data for wetland modeling. The ability to reproduce wet or dry periods (seasonal variation of water level) is mostly driven by accurate modeling of precipitation and evapotranspiration. 4.4. Impacts of spin-up periods on simulation “Spin-up” refers to the practice of allowing the model to start simulations with weather observations for a period of time before the modeling period of interest commences. This allows state variables within the model to approximate actual conditions by the time the desired modeling period begins. We investigated how different spin-up periods would affect GIW module performance. In doing so, we used the natural wetland at Site # 1 as an example and assessed the performance of the GIW module following the same procedure in Section 2.4 under six spin-up scenarios, i.e., zero, one, two, three, four, and five spin-up years. Model performance evaluation results are shown in Table 6. Calibrated GIW

Table 6 GIW module performance evaluation for different spin-up scenarios. Spin-up Years

R2

NS

Mean Sim. (mm)

Bias. (mm)

0 1 2 3 4 5

0.65 0.84 0.71 0.69 0.76 0.73

0.56 0.82 0.68 0.66 0.72 0.70

229 204 207 205 205 205

33 8 11 9 9 9

Note: mean observed water level is 196 mm.

J. Qi, X. Zhang and S. Lee et al.

Advances in Water Resources 131 (2019) 103380

Table 7 Calibrated GIW module parameters for different scenarios. Spin-up Years

Shape Factor p

Wet_K (mm h−1 )

0 1 2 3 4 5

1.89 1.88 1.75 1.77 1.89 1.86

0.044 0.038 0.035 0.036 0.035 0.034

module parameters for different scenarios are shown in Table 7. Values of R2 ranged from 0.65 to 0.84, NS ranged from 0.56 to 0.82, and Bias ranged from 8 to 33 mm; Shape factor p ranged from 1.75 to 1.89 and Wet_K ranged from 0.034 to 0.044 mm h−1 for the six spin-up scenarios. Except for the no spin-up (i.e., zero spin-up time) and one spin-up year scenarios, the ranges of the three coefficients were narrow and there are consistent patterns in model performance and calibrated parameters between different scenarios (two to five spin-up years). Therefore, a two-year spin-up period is sufficient for the GIW module to stabilize with respect to water level simulation. Similar results were observed for other three wetlands, which are not shown here for the purpose of saving space. 4.5. Limitations and potential use of the GIW module Fig. 7 shows the conceptualized wet and dry conditions for surface water storage, soil moisture, and a shallow aquifer as simulated by the GIW module for two types of GIWs. For the GIW without a low-permeability soil layer, the soil and aquifer below the water ta-

ble is always saturated during wet or dry periods (Fig. 7a and b). For the GIW with a low-permeability soil layer, the soil portion above the low-permeability layer is variably saturated during wet or dry periods (Fig. 7c and d). Due to the existence of a low-permeability layer, the soil portion below the low-permeability layer and above the groundwater table is unsaturated. When extremely dry periods occur, the developed GIW module can extract water from the shallow aquifer to address ET demand for both types of GIW. Saturated and unsaturated portions of the soil in wetlands has profound impacts on various biogeochemical processes especially during the transitional periods from wet to dry conditions or from dry to wet conditions. The newly developed GIW module can provide reasonable estimation of saturated and unsaturated portions of the soil for these sensitive periods. The ability to accurately describe hydroperiods in GIWs is important for biogeochemical modeling. For example, wetlands tend to produce more CH4 from saturated soils which are ideal anaerobic environments for fermentation and methanogenesis. Therefore, accurate estimation of CH4 emissions requires well-simulated inundation conditions for wetlands. The GIW module developed in the present study has great potential for enhancing simulation of biogeochemical cycles not only at the site scale, but also at the watershed scale. The GIW module does not consider lateral groundwater inflow (controlled by regional groundwater fluxes) which is one of the key drivers that may affect the fluctuation of water levels in the coastal plain region (Lee et al., 2019). Thus, considering lateral groundwater inflow may improve simulation of wetland hydrologic processes and understanding of wetland behaviors. The newly developed GIW module does not consider lateral groundwater inflows because the module is developed at the HRU scale where no hydraulic connectivity to adjacent HRUs is allowed in the current version of the SWAT model. The two sites used in the present study may have little external influence from regional groundwater fluxes so that the GIW module performed well. Further studies are needed to test the GIW module in wetlands under the influence of regional groundwater fluxes. In addition, the GIW module is not applicable to simulate riparian wetlands which have strong hydraulic interactions with adjacent streams and other water bodies (Rahman et al., 2016). Note that, there are existing models that employ physically based algorithms (i.e. Richard’s equation) to solve coupled surface and subsurface water dynamics. Examples include Parflow (Parallel Flow) (Maxwell et al., 2015), Cathy (CATchment HYdrology) (Camporese et al., 2010), and The HydroGeoSphere (HGS) model (Therrien et al., 2010). For example, the HGS model has been applied to conduct three-dimensional modeling to depict non-forested wetland dynamics (Ameli and Creed, 2017; Frei et al., 2010; Liu et al., 2016). Usually, those three-dimensional models require high resolution elevation data to characterize structure of wetland and flow paths, and are demanding in terms of computational resources. The GIW module is designed to efficiently perform forested GIW modeling for which high quality elevation data is often not available due to vegetation cover. The intended use of the GIW module presented here is to lower the computational cost and data quality demand for large scale watershed modeling. In particular, we tested the module within the SWAT model framework (which has been widely used in wetland modeling (Evenson et al., 2016; Wu and Johnston, 2008; Wang et al., 2008; Rahman et al., 2016; Lee et al., 2018)) and demonstrated its improvement to the existing SWAT wetland module. However, given the importance of GIW and other type of wetlands in hydrological and biogeochemical cycles, we acknowledge the need of comparison of existing models with respect to different types of wetlands, and providing assessment of different models for different conditions in the future. 5. Conclusions

Fig. 7. Conceptualized wet and dry conditions for surface water storage, soil moisture, and a shallow aquifer for the two types of GIW.

The focus of this study was to test a new GIW module within SWAT model. The novelty of the newly developed GIW module is that it

J. Qi, X. Zhang and S. Lee et al.

coupled surface water storage and subsurface water dynamics for varying inundation conditions. The GIW module was tested using observed daily water level data from two sites on the Delmarva Peninsula, i.e., Site # 1 and # 2. Each site includes a restored and a natural GIW. The two GIWs at Site #1 have a low-permeability soil layer, while the two GIWs at Site # 2 have no low-permeability soil layer. Most of the input parameters of the GIW module can be determined based on readily available data, except that we need to calibrate three parameters. Sensitivity analyses indicated the saturated hydraulic conductivity, Wet_K, for the lowpermeability soil layer (in wetlands at Site # 1) and the interface layer between the shallow aquifer and the deep aquifer (in wetlands at Site # 2) was the most sensitive parameter. The calibrated values of Wet_K were about two order of magnitude lower than reported (about 0.03 vs. 3 mm h−1 ) in the USDA Soil Survey Geographic Database (SSURGO) suggesting that SSURGO data may not accurate for wetlands. Another sensitive parameter for the wetlands without a low-permeability soil layer was the effective hydraulic conductivity keff . We recognized that there was a large uncertainty in calibrated keff which was preferably used as a fitting parameter in the present study. In addition, after examining the impact of different spin-up periods on module performance, we found that twoyear spin-up period was enough for the GIW module to make reliable simulation of water levels. The simulated daily water level and observed data were compared and evaluated using three statistics, i.e., Bias, the Nash-Sutcliffe coefficient (NS), and the coefficient of determination (R2 ). We used monitored water level data in 2016 to calibrate the GIW module at two sites, and data in 2017 from Site # 2 were used to validate the module. Model calibration and validation results show that, in general, the newly developed GIW module reproduced well hydroperiods for both restored and natural wetlands at the two study sites. The module also performed well for wetlands with and without a low-permeability soil layer. The GIW module could reasonably capture the beginning of surface water disappearance and surface water accumulation after dry periods; and saturated and unsaturated portions of the soil corresponding to wet and dry periods were well described by the GIW module. The ability to accurately describe inundation conditions for wetlands is important for biogeochemical modeling. Therefore, the newly developed GIW module has a strong potential to enhance simulation of biogeochemical cycles at both the site and watershed scale. As the SWAT model has been used widely in understanding hydrologic and biogeochemical significance of wetlands in the watershed context, the new GIW holds potential to benefit future studies in this aspect. We also note that the GIW module may not be applicable to GIWs strongly influenced by larger-scale groundwater fluxes and would likely not be suitable to simulate riparian wetlands. Acknowledgments We greatly appreciate the valuable comments from three anonymous reviewers, which helped us to improve the quality of the paper. The funding support for this project was provided by NASA (NNX17AE66G and NNH13ZDA001N), and USDA (2017-67003-26485), and NSF INFEWS (1639327). Funding was also provided in part by the USDA Natural Resources Conservation Service - Conservation Effects Assessment Project (NRCS-CEAP). Disclaimer The U.S. Department of Agriculture is an equal opportunity provider and employer. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government. Supplementary materials Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.advwatres.2019.103380.

Advances in Water Resources 131 (2019) 103380

References Abbaspour, K., Vejdani, M., Haghighat, S., 2007. SWAT-CUP calibration and uncertainty programs for SWAT. MODSIM 2007 International Congress on Modelling and Simulation, Modelling and Simulation Society of. Abbaspour, K.C., Vaghefi, S.A., Srinivasan, R., 2017. A Guideline for Successful Calibration and Uncertainty Analysis for Soil and Water Assessment: A Review of Papers from the 2016 International SWAT Conference. Multidisciplinary Digital Publishing Institute. Ameli, A.A., Creed, I.F., 2017. Quantifying hydrologic connectivity of wetlands to surface water systems. Hydrol. Earth Syst. Sci. 21. Arnold, J.G., Srinivasan, R., Muttiah, R.S., Williams, J.R., 1998. Large area hydrologic modeling and assessment part I: model development. J. Am. Water Resou. Assoc. 34, 73–89. Ator, S.W., Denver, J.M., 2015. Understanding the Nutrients in the Chesapeake Bay Watershed and Implications for Management and Restoration: The Eastern Shore. US Department of the Interior, US Geological Survey. Babbar-Sebens, M., Barr, R.C., Tedesco, L.P., Anderson, M., 2013. Spatial identification and optimization of upland wetlands in agricultural watersheds. Ecol. Eng. 52, 130–142. Bantilan-Smith, M., Bruland, G.L., MacKenzie, R.A., Henry, A.R., Ryder, C.R., 2009. A comparison of the vegetation and soils of natural, restored, and created coastal lowland wetlands in Hawai ‘i. Wetlands 29, 1023–1035. Bear, J., 2012. Hydraulics of Groundwater. Courier Corporation. Bridgham, S.D., Johnston, C.A., Pastor, J., Updegraff, K., 1995. Potential feedbacks of northern wetlands on climate change. Bioscience 45, 262–274. Bridgham, S.D., Patrick Megonigal, J., Keller, J.K., Bliss, N.B., Trettin, C., 2006. The carbon balance of North American wetlands. Wetlands 26, 889–916. Campbell, D.A., Cole, C.A., Brooks, R.P., 2002. A comparison of created and natural wetlands in Pennsylvania, USA. Wetlands Ecol. Manage. 10, 41–49. Camporese, M., Paniconi, C., Putti, M., Orlandini, S., 2010. Surface-subsurface flow modeling with path-based runoff routing, boundary condition-based coupling, and assimilation of multisource observation data. Water Resour. Res. 46. Chanton, J.P., Bauer, J.E., Glaser, P.A., Siegel, D.I., Kelley, C.A., Tyler, S.C., et al., 1995. Radiocarbon evidence for the substrates supporting methane formation within northern Minnesota peatlands. Geochim. Cosmochim. Acta 59, 3663–3668. Dahl, T.E., Johnson, C.E., Frayer, W., 1991. Wetlands, Status and Trends in the Conterminous United States Mid-1970′s to mid-1980′s. US Fish and Wildlife Service. Dahl, T.E., 2014. Status and Trends of Prairie Wetlands in the United States 1997 to 2009. US Fish and Wildlife Service. Dodds, W.K., Wilson, K.C., Rehmeier, R.L., Knight, G.L., Wiggam, S., Falke, J.A., et al., 2008. Comparing ecosystem goods and services provided by restored and native lands. AIBS Bull. 58, 837–845. Evenson, G.R., Golden, H.E., Lane, C.R., D’amico, E., 2016. An improved representation of geographically isolated wetlands in a watershed-scale hydrologic model. Hydrol. Process. 30, 4168–4184. Faulkner, S.P., Richardson, C.J., 1989. Physical and chemical characteristics of freshwater wetland soils. In: Hammer, D.A. (Ed.), Constructed Wetlands for Wastewater Treatment: Municipal, Industrial, and Agricultural. Lewis Publishers, Chelsea, pp. 41–72. Frei, S., Lischeid, G., Fleckenstein, J., 2010. Effects of micro-topography on surface–subsurface exchange and runoff generation in a virtual riparian wetland—A modeling study. Adv. Water Resour. 33, 1388–1401. Hargreaves, G.H., Samani, Z.A., 1985. Reference crop evapotranspiration from temperature. Appl. Eng. Agric. 1, 96–99. Hefting, M.M., van den Heuvel, R.N., Verhoeven, J.T., 2013. Wetlands in agricultural landscapes for nitrogen attenuation and biodiversity enhancement: opportunities and limitations. Ecol. Eng. 56, 5–13. Holden, J., Burt, T., 2003. Hydraulic conductivity in upland blanket peat: measurement and variability. Hydrol. Process. 17, 1227–1237. Holden, J., Chapman, P., Labadz, J., 2004. Artificial drainage of peatlands: hydrological and hydrochemical process and wetland restoration. Prog. Phys. Geogr. 28, 95–123. Ingram, H., 1967. Problems of hydrology and plant distribution in mires. J. Ecol. 711–724. Javaheri, A., Babbar-Sebens, M., 2014. On comparison of peak flow reductions, flood inundation maps, and velocity maps in evaluating effects of restored wetlands on channel flooding. Ecol. Eng. 73, 132–145. Johnston, C.A., 1991. Sediment and nutrient retention by freshwater wetlands: effects on surface water quality. Crit. Rev. Environ. Sci. Technol. 21, 491–565. Jordan, T.E., Whigham, D.F., Hofmockel, K.H., Pittek, M.A., 2003. Nutrient and sediment removal by a restored wetland receiving agricultural runoff. J. Environ. Qual. 32, 1534–1547. Keddy, P.A., 2010. Wetland Ecology: Principles and Conservation. Cambridge University Press. Keough, J.R., Thompson, T.A., Guntenspergen, G.R., Wilcox, D.A., 1999. Hydrogeomorphic factors and ecosystem responses in coastal wetlands of the Great Lakes. Wetlands 19, 821–834. Klute, A., Dirksen, C., 1986. Hydraulic conductivity and diffusivity: laboratory methods. In: Methods of Soil Analysis: Part 1—Physical and Mineralogical Methods. Soil Science Society of America, American Society of Agronomy, pp. 687–734. Krause, P., Boyle, D., Bäse, F., 2005. Comparison of different efficiency criteria for hydrological model assessment. Adv. Geosci. 5, 89–97. Kuenzler, E.J., 1989. Value of forested wetlands as filters for sediments and nutrients. In: Proc Symposium: The Forested Wetlands of the Southern United States, pp. 85–96. Lane, C.R., D’Amico, E., 2016. Identification of putative geographically isolated wetlands of the conterminous United States. J. Am. Water Resour. Assoc. 52, 705–722. Lee, S., Yeo, I.-Y., Sadeghi, A.M., McCarty, G.W., Hively, W.D., Lang, M.W., 2016. Impacts of watershed characteristics and crop rotations on winter cover crop nitrate-nitrogen

J. Qi, X. Zhang and S. Lee et al. uptake capacity within agricultural watersheds in the Chesapeake Bay region. PLoS ONE 11, e0157637. Lee, S., Yeo, I.-Y., Lang, M., McCarty, G., Sadeghi, A., Sharifi, A., et al., 2017. Improving the catchment scale wetland modeling using remotely sensed data. Environ. Model. Softw. in press. https://doi.org/10.1016/j.envsoft.2017.11.001. Lee, S., Yeo, I.-Y., Lang, M., Sadeghi, A., McCarty, G., Moglen, G., et al., 2018. Assessing the cumulative impacts of geographically isolated wetlands on watershed hydrology using the SWAT model coupled with improved wetland modules. J. Environ. Manage. 223, 37–48. Lee, S., McCarty, G.W., Moglen, G.E., Lang, M.W., Sadeghi, A.M., Green, T.R., et al., 2019. Effects of subsurface soil characteristics on wetland–groundwater interaction in the coastal plain of the Chesapeake Bay watershed. Hydrol. Process. 33, 305–315. Li, Q., Qi, J., Xing, Z., Li, S., Jiang, Y., Danielescu, S., et al., 2014. An approach for assessing impact of land use and biophysical conditions across landscape on recharge rate and nitrogen loading of groundwater. Agric. Ecosyst. Environ. 196, 114–124. Liang, K., Qi, J., Liu, E., Jiang, Y., Li, S., Meng, F.-R., 2019. Estimated potential impacts of soil and water conservation terraces on potato yields under different climate conditions. J. Soil Water Conserv. 74, 225–234. Liu, G., Schwartz, F.W., Wright, C.K., McIntyre, N.E., 2016. Characterizing the climate– driven collapses and expansions of wetland habitats with a fully integrated surface— subsurface hydrologic model. Wetlands 36, 287–297. Martinez-Martinez, E., Nejadhashemi, A.P., Woznicki, S.A., Love, B.J., 2014. Modeling the hydrological significance of wetland restoration scenarios. J. Environ. Manage. 133, 121–134. Masterson, J.P., Pope, J.P., Fienen, M.N., Monti Jr, J., Nardi, M.R., Finkelstein, J.S., 2016. Documentation of a Groundwater Flow Model Developed to Assess Groundwater Availability in the Northern Atlantic Coastal Plain Aquifer System from Long Island. US Geological Survey, New York, to North Carolina. Maxwell, R., Condon, L., Kollet, S., 2015. A high-resolution simulation of groundwater and surface water over most of the continental US with the integrated hydrologic model ParFlow v3. Geosci. Model Dev. 8, 923. McLaughlin, D.L., Kaplan, D.A., Cohen, M.J., 2014. A significant nexus: geographically isolated wetlands influence landscape hydrology. Water Resour. Res. 50, 7153–7166. Mitsch, W.J., Ewel, K.C., 1979. Comparative Biomass and Growth of Cypress in Florida Wetlands. American Midland Naturalist, pp. 417–426. Monteith, J.L., 1965. Evaporation and environment. Symp. Soc. Exp. Biol. 4. Narasimhan, B., Srinivasan, R., Arnold, J., Di Luzio, M., 2005. Estimation of long-term soil moisture using a distributed parameter hydrologic model and verification using remotely sensed data. Trans. ASAE 48, 1101–1113. Nash, J., Sutcliffe, J., 1970. River flow forecasting through conceptual models part I—a discussion of principles. J. Hydrol. 10, 282–290. Neitsch, S.L., Arnold, J.G., Kiniry, J.R., Williams, J.R., 2011. Soil and Water Assessment Tool Theoretical Documentation Version 2009. Texas Water Resources Institute. Neitsch, S.L., Williams, J.R., Arnold, J.G., Kiniry, J.R., 2011. Soil and Water Assessment Tool Theoretical Documentation Version 2009. Grassland, Soil and Water Research Service, Temple, Texas, USA. Oleson K.W., D.M. Lawrence, B. Gordon, M.G. Flanner, E. Kluzek, J. Peter, et al. Technical description of version 4.0 of the Community Land Model (CLM). (2010). Priestley, C., Taylor, R., 1972. On the assessment of surface heat flux and evaporation using large-scale parameters. Mon. Weather Rev. 100, 81–92. Qi, J., Li, S., Jamieson, R., Hebb, D., Xing, Z., Meng, F.-R., 2017. Modifying SWAT with an energy balance module to simulate snowmelt for maritime regions. Environ. Model. Softw. 93, 146–160. Qi, J., Li, S., Yang, Q., Xing, Z., Meng, F.-R., 2017. SWAT setup with longterm detailed landuse and management records and modification for a microwatershed influenced by freeze-thaw cycles. Water Resour. Manage. 31, 3953–3974. https://doi.org/10.1007/s11269-017-1718-2. Qi, J., Li, S., Bourque, C.P., Xing, Z., Fan-Rui, M., 2018. Developing a decision support tool for assessing land use change and BMPs in ungauged watersheds based on decision rules provided by SWAT simulation. Hydrol. Earth Syst. Sci. 22, 3789–3806. Qi, J., Zhang, X., McCarty, G.W., Sadeghi, A.M., Cosh, M.H., Zeng, X., et al., 2018. Assessing the performance of a physically-based soil moisture module integrated within the Soil and Water Assessment Tool. Environ. Model. Softw. 109, 329–341.

Advances in Water Resources 131 (2019) 103380 Qi, J., Zhang, X., Wang, Q., 2019. Improving hydrological simulation in the Upper Mississippi River Basin through enhanced freeze-thaw cycle representation. J. Hydrol. 571, 605–618. Qi, J., Zhang, X., Cosh, M.H., 2019. Modeling soil temperature in a temperate region: a comparison between empirical and physically based methods in SWAT. Ecol. Eng. 129, 134–143. Qi, J., Wang, Q., Zhang, X., 2019. On the use of NLDAS2 weather data for hydrologic modeling in the Upper Mississippi River Basin. Water 11, 960. Rahman, M.M., Thompson, J.R., Flower, R.J., 2016. An enhanced SWAT wetland module to quantify hydraulic interactions between riparian depressional wetlands, rivers and aquifers. Environ. Model. Softw. 84, 263–289. Rains, M., Leibowitz, S., Cohen, M., Creed, I., Golden, H., Jawitz, J., et al., 2016. Geographically isolated wetlands are part of the hydrological landscape. Hydrol. Process. 30, 153–160. Rubec, C., 1994. Canada’s federal policy on wetland conservation: a global model. In: Mitsch, W.J. (Ed.), Global Wetlands: Old World and New. Elsevier, pp. 909–917. Russell, K.N., Beauchamp, V.B., 2017. Plant species diversity in restored and created Delmarva bay wetlands. Wetlands 37, 1119–1133. Rycroft, D., Williams, D., Ingram, H., 1975. The transmission of water through peat: I. Review. J. Ecol. 535–556. Schlesinger, W.H., 2005. Biogeochemistry. Elsevier. Smith, R.D., Ammann, A., Bartoldus, C., Brinson, M.M., 1995. An Approach for Assessing Wetland Functions Using Hydrogeomorphic Classification, Reference Wetlands, and Functional Indices. Army Engineer Waterways Experiment Station, Vicksburg, MS. Solomon, S., 2007. Climate Change 2007-The Physical Science Basis: Working Group I Contribution to the Fourth Assessment Report of the IPCC.. Cambridge University Press. Srinivasan, R., Zhang, X., Arnold, J., 2010. SWAT ungauged: hydrological budget and crop yield predictions in the Upper Mississippi River Basin. Trans. ASABE 53, 1533–1546. Therrien, R., McLaren, R., Sudicky, E., Panday, S., 2010. HydroGeoSphere: A Three-Dimensional Numerical Model Describing Fully-Integrated Subsurface and Surface Flow and Solute Transport. Groundwater Simulations Group. University of Waterloo, Waterloo, ON. Tiner, R.W., 2003. Geographically isolated wetlands of the United States. Wetlands 23, 494–516. Ullah, S., Breitenbeck, G., Faulkner, S.P., 2005. Denitrification and N 2 O emission from forested and cultivated alluvial clay soil. Biogeochemistry 73, 499–513. Wang, X., Yang, W., Melesse, A.M., 2008. Using hydrologic equivalent wetland concept within SWAT to estimate streamflow in watersheds with numerous wetlands. Trans. ASABE 51, 55–72. Wu, K., Johnston, C.A., 2008. Hydrologic comparison between a forested and a wetland/lake dominated watershed using SWAT. Hydrol. Process. 22, 1431–1442. Yepsen, M., Baldwin, A.H., Whigham, D.F., McFarland, E., LaForgia, M., Lang, M., 2014. Agricultural wetland restorations on the USA Atlantic coastal plain achieve diverse native wetland plant communities but differ from natural wetlands. Agric. Ecosyst. Environ. 197, 11–20. Zedler, J.B., Kercher, S., 2005. Wetland resources: status, trends, ecosystem services, and restorability. Annu. Rev. Environ. Resour. 30, 39–74. Zeng, X., Decker, M., 2009. Improving the numerical solution of soil moisture–based Richards equation for land models with a deep or shallow water table. J. Hydrometeorol. 10, 308–319. Zhang, Y., Li, C., Trettin, C.C., Li, H., Sun, G., 2002. An integrated model of soil, hydrology, and vegetation for carbon dynamics in wetland ecosystems. Glob. Biogeochem. Cycles 16. Zhang, X., Srinivasan, R., Debele, B., Hao, F., 2008. Runoff simulation of the headwaters of the Yellow River using the SWAT model with three snowmelt algorithms. J. Am. Water Resour. Assoc. 44, 48–61. Zhang, C., Li, S., Qi, J., Xing, Z., Meng, F., 2017. Assessing impacts of riparian buffer zones on sediment and nutrient loadings into streams at watershed scale using an integrated REMM-SWAT model. Hydrol. Process. 31, 916–924. Zimmer, K.D., Hanson, M.A., Butler, M.G., 2000. Factors influencing invertebrate communities in prairie wetlands: a multivariate approach. Can. J. Fish. Aquat. Sci. 57, 76–85.