Riparian Zone Nitrogen Management through the Development of the Riparian Ecosystem Management Model (REMM) in a Formerly Glaciated Watershed of the US Northeast

: The Riparian Ecosystem Management Model (REMM) was developed, calibrated and validated for both hydrologic and water quality data for eight riparian buffers located in a formerly glaciated watershed (upper Pawcatuck River Watershed, Rhode Island) of the US Northeast. The Annualized AGricultural Non-Point Source model (AnnAGNPS) was used to predict the runoff and sediment loading to the riparian buffer. Overall, results showed REMM simulated water table depths (WTDs) and groundwater NO 3 -N concentrations at the stream edge (Zone 1) in good agreement with measured values. The model evaluation statistics showed that, hydrologically REMM performed better for site 1, site 4, and site 8 among the eight buffers, whereas REMM simulated better groundwater NO 3 -N concentrations in the case of site 1, site 5, and site 7 when compared to the other ﬁve sites. The interquartile range of mean absolute error for WTDs was 3.5 cm for both the calibration and validation periods. In the case of NO 3 -N concentrations prediction, the interquartile range of the root mean square error was 0.25 mg/L and 0.69 mg/L for the calibration and validation periods, respectively, whereas the interquartile range of d for NO 3 -N concentrations was 0.20 and 0.48 for the calibration and validation period, respectively. Moreover, REMM estimation of % N-removal from Zone 3 to Zone 1 was 19.7%, and 19.8% of N against actual measured 19.1%, and 26.6% of N at site 7 and site 8, respectively. The sensitivity analyses showed that changes in the volumetric water content between ﬁeld capacity and saturation (soil porosity) were driving water table and denitriﬁcation.


Introduction
Riparian zones occur at the interface of terrestrial and aquatic components of the landscape. They regularly receive and process large amounts of excess nitrogen (N), draining out of agricultural fields towards open water bodies. They are often characterized as "filters" or "buffers" and are vital elements in watershed management schemes for water quality maintenance and stream ecosystem habitat protection [1][2][3].
Agriculture (cropland, pasture, managed forest) is an important component of many watersheds of the USA Northeast where N losses to major estuaries is of substantial concern. Decades of research on riparian zone hydrology and biogeochemistry has shown that riparian zones can serve as best management practices (BMPs) to mitigate the impact of agriculture (excessive leaching of nutrients, mostly N) on the quality of our waters [4][5][6]. Nevertheless, the buffering capacity of riparian zones (mostly for N) varies enormously due to the hydrogeomorphic setting such as topography, depth to water table, soil type, and surficial geology of the riparian zone [7][8][9][10][11][12]. Upland land use/land cover affects both the water quantity and quality of the water entering the riparian zone. Hydrogeomorphic model's capacity for water table depth and groundwater nitrate concentration simulation through calibration of the developed model by means of comparing model outputs with field data collected from eight buffer sites in RI, (iii) validation of model's output with the improved calibrated parameters, (iv) conduct a parameter sensitivity analysis. Ultimately, this approach will facilitate the use of this model in this region and improve its functionality with respect to nitrogen transformations and flux.

Site Description
Our study focused on eight riparian sites from the state of RI, USA. All the sites are located in upper Pawcatuck River Watershed, Washington County, set in the New England Hydrologic Region of southern RI (41 • 32 30 N, 71 • 35 W) ( Figure 1) and all were monitored for hydrology and water quality. The area of this watershed is about 258 km 2 . It consists mainly of forests (above 65% of the total watershed) and agricultural fields (about 32%). The soil parent materials in the watershed are comprised mostly of glacial till, glacial outwash, and organic and alluvial deposit [13]. Agricultural lands (mostly turf farms) are predominately located on loess soils over glacial outwash. Forested settings are usually on till. The elevation of the watershed ranges from 16 m (shoreline) to 144 m (gently rolling hills inland). It has a humid continental climate, with warm summers and cold winters. The temperature in the watershed ranges from a 30-year (1982- [41]; stream lines-shapefile was downloaded from open source Rhode Island Geographic Information System [42]; Pawcatuck River Watershed boundary was generated by a subset of TOPAZ, TOPAGNPS (the set of TOPAZ modules used for AGNPS)).
All sites were forested riparian wetlands dominated by red maple (Acer rubrum L.). Half the sites were located on first or second streams, one site was located along a pond, two sites were on a 4th order stream and one riparian site bordered an intermittent stream. The upland land use at seven of the sites was for commercial turf operations and only one site had forested uplands. Irrigation was routinely applied on a number of turf farms. Details on the soils, slope and buffer dimensions are provided in Tables 1 and 2.  [41]; stream lines-shapefile was downloaded from open source Rhode Island Geographic Information System [42]; Pawcatuck River Watershed boundary was generated by a subset of TOPAZ, TOPAGNPS (the set of TOPAZ modules used for AGNPS)).
All sites were forested riparian wetlands dominated by red maple (Acer rubrum L.). Half the sites were located on first or second streams, one site was located along a pond, two sites were on a 4th order stream and one riparian site bordered an intermittent stream. The upland land use at seven of the sites was for commercial turf operations and only one site had forested uplands. Irrigation was routinely applied on a number of turf farms. Details on the soils, slope and buffer dimensions are provided in Tables 1 and 2.

Description of The REMM Model
REMM is a field-scale process-based, two dimensional, daily time-step model that simulates interactions between hydrology, nutrient dynamics, sediment transport, and vegetation growth. REMM computes the loading of water, sediments, carbon, and nutrients coming from the upland into the riparian buffer. The model was designed such that water Agriculture 2021, 11, 743 5 of 24 and total N are transported from upland to field edge (Zone 3), field edge to mid-buffer (Zone 2), mid-buffer to stream edge (Zone 1), and ultimately from stream edge to open water body by means of surface runoff, seep flow, and subsurface flow [15,16].
REMM file version 0.1.1.46 (United States Department of Agriculture (USDA)-Agricultural Research Service (ARS)) was used for all simulations. A broad description of the REMM model is available in several publications, including [12,15,16,35]. Briefly, REMM is a computer simulation model of riparian forest buffer systems. The structure of REMM is consistent with the three zone riparian system as mentioned in [4]. REMM was originally field tested using a five-year hydrologic and nutrient dataset collected from an experimental riparian buffer site in Tifton, Georgia [12,35]. We used this model, and prepared the model set-up for each riparian site, parameterized, calibrated and validated for both hydrologic and nutrient simulation.
Within the model, the riparian system is considered to consist of three zones (parallel to a stream) between the field and the water body. However, zone 2 is not visibly distinguished at the riparian sites in this study, as they tend to move abruptly from zone 1 to zone 3 ( Figure 2). Each zone includes litter and three soil layers (through which the vertical and horizontal movement of water takes place) that terminate at the bottom of the plant root system, and a plant community that can include six plant types in two canopy levels. The riparian system characterized in REMM was originally designed to represent increasing levels of management away from the stream [15]. REMM is written in the C++ programming language.
Agriculture 2021, 11, x FOR PEER REVIEW 6 of 25 represent increasing levels of management away from the stream [15]. REMM is written in the C++ programming language. In the REMM module, movement of water and storage is defined by several processes, i.e., interception, evapotranspiration (ET), infiltration, vertical drainage, surface runoff, subsurface lateral flow, upward flux from the water table in response to ET, and seepage or exfiltration. These processes are simulated for Zone 3, Zone 2, and Zone 1. The water movement and storage between the zones is based on a combination of mass balance and rate-controlled approaches. Equation (1) presents the mass balance of water within each soil layer: where (mm) is the soil moisture on day t, (mm) is the soil moisture from the preceding day, (mm) is the addition of water as a result of infiltration in case of the upper soil layer, or drainage from upper soil layer for intermediate soil layers, (mm) is drainage out of the layer, (mm) is contribution because of lateral subsurface flow, (mm) is the outflow of water to lateral subsurface downslope flow, and ET (mm) is evapotranspiration [12]. REMM hydrologic outputs generated from the water balance simulation include daily surface and subsurface losses to the water body, evapotranspiration and deep seepage. Deep seepage is specified by the user for each zone as input only, such that the water that is lost through the deep seepage never comes back into REMM computations [16,35]. In the REMM module, movement of water and storage is defined by several processes, i.e., interception, evapotranspiration (ET), infiltration, vertical drainage, surface runoff, subsurface lateral flow, upward flux from the water table in response to ET, and seepage or exfiltration. These processes are simulated for Zone 3, Zone 2, and Zone 1. The water movement and storage between the zones is based on a combination of mass balance and rate-controlled approaches. Equation (1) presents the mass balance of water within each soil layer: where SM (t) (mm) is the soil moisture on day t, SM (t−1) (mm) is the soil moisture from the preceding day, Q V−in (mm) is the addition of water as a result of infiltration in case of the upper soil layer, or drainage from upper soil layer for intermediate soil layers, Q V−out (mm) is drainage out of the layer, Q H−in (mm) is contribution because of lateral subsurface flow, Q H−out (mm) is the outflow of water to lateral subsurface downslope flow, and ET (mm) is evapotranspiration [12]. REMM hydrologic outputs generated from the water balance simulation include daily surface and subsurface losses to the water body, evapotranspiration and deep seepage. Deep seepage is specified by the user for each zone as input only, such that the water that is lost through the deep seepage never comes back into REMM computations [16,35]. Reference [35] described the equation used in REMM model to simulate denitrification. Denitrification is calculated as the function of the interaction of factors representing the degree of anaerobiosis, temperature, nitrate-N, and available carbon: where K d is the rate of denitrification under optimal conditions (kgcm −1 ha −1 ), S d is the depth of the soil layer (cm), A f is the scaler factor representing the effect of anaerobiosis on denitrification (0-1), T denitrification,t is the scaler factor representing the effect of temperature on denitrification (0-1), N f is the scaler factor representing the effect of nitrate-N on denitrification (0-1), C f is the scaler factor representing the effect of available carbon on denitrification (0-1), and α is a coefficient determining the influence of nitrate on denitrification set at 0.19. According to [43], REMM has been developed as a hillslope-scale mechanistic model to predict how the width and composition of riparian habitats impact material loadings to streams. Particularly, when loadings of sediment and nutrients to the riparian zone are known, then REMM can be used to simulate the effect of riparian buffers on stream chemistry. Similar to other mechanistic models, the application of REMM is limited by its complexity and large needs for input data. REMM operates at the hillslope rather than catchment scale. It cannot forecast effects on instream ecological endpoints.

REMM Model Input Data
Parameters are input into the model through four basic files that include: (1) contributions of daily outputs from the field draining into the riparian system including surface runoff and associated eroded sediment, organic material and plant nutrients (*.FIN), (2) weather data (*.WEA), (3) vegetation (*.VEG) characteristics, and (4) soil physical and chemical parameters (*.BUF). The last two of these describe in detail the morphology of the buffer system being modeled. One of the key reliable inputs needed to simulate riparian system in the REMM model is the water and nutrients originating from the upland source area. These can be field measured or simulated using models such as AnnAGNPS (AGricultural Non-Point Source) and APEX (Agricultural Policy/Environmental eXtender). We had to rely on AnnAGNPS simulated upland data due to lack of field measured data (surface runoff and sediment data) at the edge of the field. We calibrated and validated the AnnAGNPS generated runoff against USGS measured streamflow at the watershed outlet and used as field input to REMM. These calibrated inputs will increase the reliability of REMM buffer simulations.

Upland Inputs (*.FIN File)
Upland inputs comprise the daily flux of upland water, associated sediment, sedimentborne chemicals, and dissolved chemicals entering the upper side of the buffer system during the period of simulation. REMM requires daily subsurface and surface flow data from the field or upland area. In general, this information is missing from most REMM studies, which can hinder the ability of REMM to properly predict riparian functions at the site scale [39]. To offset this issue, a field-scale hydrological model, Annualized AGricultural Non-Point Source (AnnAGNPS), was used to predict the runoff and sediment loading to the riparian buffer [44,45]. AnnAGNPS [46,47] is a daily time step, watershed scale, pollutant-loading, distributed model developed to simulate long-term runoff, sediment, nutrients, and pesticide transport from agricultural watersheds [48][49][50]. The AnnAGNPS model defines cells of various sizes; contaminants are routed from these cells into the associated reaches, and the model either deposits pollutants within the stream channel system or transports them out of the watershed [47]. TOPAZ (one of the modules of AnnAGNPS) is the TOpographic PArameteriZation program which generates cell and stream network information from the watershed digital elevation model (DEM). It also provides all of the topographic related information for AnnAGNPS. DEM was acquired from the United States Geological Survey (USGS)-The National Map Viewer (TNM Viewer version 2.0, USGS, Washington, DC, USA)7.5-min digital elevation models (DEMs)-with a 10-m horizontal, 7-m vertical resolution [51]. The soil map [52] and land use map [53] were also incorporated for watershed delineation. The simulated runoff and sediment loading from AnnAGNPS was calibrated by comparing with observed data (USGS gauge). The calibrated daily runoff and the sediment loading were used as input data in .FIN (Field Data File). The rationale behind the use of AnnAGNPS model rather than other hydrological models, such as SWAT, is that AnnAGNPS divides the entire watershed into a number of cells and generates cell wise upland input data file compatible to *.FIN file format in REMM model. For this study, AnnAGNPS was used with a cell size of 1.8-0.4 km 2 (interquartile range, IQR = Q 3 − Q 1, of cell size among a total of 185 cells in the watershed) to simulate input cells representing upland inputs to each of the eight riparian buffers within upper Pawcatuck River watershed. Details of the application of AnnAGNPS model on the study area can be found in our companion study [45]. Table 1 represents the upland characteristics for all the riparian sites. The sites are all in glacial outwash with sandy loam or fine sandy loam soil.

Weather Data (*.WEA File)
REMM requires daily rainfall, maximum and minimum air temperature, solar radiation, wind speed, and dew point temperature as its climate input files. All these data are obtained from the closest available national weather service cooperative observer program (NWS COOP) weather stations [54,55]. The data for three daily climate parameters, i.e., minimum air temperature, maximum air temperature, and precipitation, were acquired from PRISM website at 4 km spatial resolution (PRISM Climate Group, Oregon State University, http://prism.oregonstate.edu, accessed on 4 February 2004).

Vegetation Data (*.VEG File)
The vegetation (.VEG) data file contains plant specific information. Regional data sets are being developed that define plant characteristics typical of a region. The vegetation database created by REMM developers (USDA-ARS) separately for northeastern region was used during the simulation period. It is also accessible to other users. Maximum rooting depths (MRDs) influenced water uptake and plant transpiration which influenced ET and simulated WTDs. A maximum rooting depth of 200 cm was used for all zones for all plant species, as all riparian sites were forested, with hardwood red maple (Acer rubrum L.) the dominant species. This was the default value used by the model developers and also by the study of [33]. Therefore, whenever no local data were available, literature values were used. Initially we looked at the influence of MRDs on the simulated outputs, but found WTDs and nitrate concentration did change at a very low rate, so we did not include MRDs in our sensitivity part. The use of MRD as 200 cm in our study provided reasonable output for our area of study, so we did not change this parameter in our model. Ref. [15] also did not include MRDs as a parameter change input in their sensitivity analysis for streamflow, total N out, denitrification and N uptake. Besides, Ref. [56] stated from their REMM sensitivity study that REMM's comparatively low sensitivity to vegetation parameters supports the use of regional vegetation datasets that would make model implementation simple without compromising results. A specific leaf area of 0.0045 ha/kg C was used for all zones for the buffer in all sites. Most of the modification required for simulation was performed within the main data file (.BUF). Plant, litter and soil layer information are given in this data file. Soil characteristics for each of the three zones with three soil layers are entered in BUFFER DATA FILE (*. BUF). The characteristics include S10 Fraction, bubbling pressure, pore size distribution index, layer thickness, wilting point, field capacity, soil porosity, permeability (cm/hr), % sand, % silt, % clay, bulk density, pH, base saturation, etc. The data for most of these soil characteristics were obtained from REMM user manual based on the soil texture [57]. Also whenever there is no data available default value or value from published literature was used. The soil layers in the model are intended to correspond with horizons in the soil profile. According to [35], REMM keeps the soil physical properties constant during the simulation period. Soil parameters (saturated hydraulic conductivity, bulk density) which might change during conversion of cropland to permanent buffers may be changed at model user defined points in the simulation or can be set as intermediate point between field and mature buffer conditions. Ref. [15] stated that REMM is designed to be used at a hillslope scale to simulate the effects of buffer systems on edge-of-field loadings. REMM takes upland outputs supplied by the user and calculates loadings of water, nutrients, sediment, and carbon based on actual area of the zones of a buffer system. Similarly for plant and litter information, modifications were based on the guidelines described in the REMM User's Manual [57]. Whenever no local data were available, literature values were used. Table 2 shows the site characteristics of the modeled buffer for all riparian sites. All sites were located on hydric soils with alluvial or glacial outwash parent materials. Although the land uses varied among watersheds, all riparian sites were forested, with red maple (Acer rubrum L.), the dominant species [58]. More details about land use are available in [11,58].

Collection of Field Data and Other Essential Inputs for REMM
All original or field data (including site characteristics) was collected by a team led by the authors of this manuscript, with decades of published research on riparian zones in glaciated settings of the USA Northeast regions, for the development and calibration of the REMM model. Water table levels were recorded in water table wells in the riparian zone, biweekly during spring and fall when water table depths were expected to change most rapidly, and bi-monthly during summer and winter. A network of mini-piezometer nests was installed across the riparian zone from the upland to the stream. The minipiezometers allowed the collection of nitrate samples at discrete depths and the examination of groundwater denitrification in situ using the push-pull method [59]. The information regarding in-situ groundwater denitrification capacity measurement in the study sites 1, 2 [11] and sites 3,4 [58] were based on the studies done by [11,58]. For sites 1&2, the in-situ groundwater denitrification rates measured are within the range reported by previous studies [59]. For sites 3&4, in situ groundwater denitrification capacity measured in the shallow wells was not significantly different from that measured in deep wells; thus, shallow and deep wells were pooled (i.e., combined) for statistical analysis. A soil pit was dug and soil samples were taken from all soil horizons and analyzed for carbon content. Particle size distribution (percentage sand, silt, and clay) and soil carbon were determined for samples taken from the soil pits. The percentage sand, silt, and clay were used as inputs into REMM. The depth of the stream as 0.305 m has been used for all the sites [6].
Groundwater samples were analyzed for NO 3 -N using the SM 4500 NO 3 F automated cadmium reduction method on an Alpkem RFA 300 Rapid Flow Auto-analyzer (O.I. Analytical, Wilsonville, OR, USA) [11,58]. Soil samples were also examined for denitrification enzyme activity (DEA). The DEA is the potential denitrification measured under fully anaerobic conditions and excess NO 3 -N and available carbon. The denitrification rate constant (K d ) is based on the denitrification potential measurement [60]. This was the only user input in REMM for simulating denitrification in all zones and layers within the buffer. For sites (1 and 2) and sites (3 and 4), field data (depth to water table, groundwater NO 3 -N concentrations) were collected for a five-year period (1999-2003) and a two-year period (2004)(2005), respectively [11,58,61]. For sites (5, 6, 7 and 8), a simplified methodology (three-well approach) had been followed according to procedures described in [62] for the collection of field data (WTDs and NO 3 -N concentrations in the groundwater) for year 2018-2019.
The erosion factors, topsoil condition, Manning's n, soil structural characteristics, soil pH and temperature, permeability class (referring to the saturated hydraulic conductivity in the soil profile), surface condition, inter-rill roughness, soil and litter characteristics, including physical properties related to soils by soil texture (bulk density, porosity, field capacity, wilting point), pore size distribution index and bubbling pressure were obtained from the REMM user's manual based on soil texture of riparian sites [55]. Additional buffer topographic inputs consisting of buffer width, slope, and length and stream depth were obtained from a detailed topographic survey conducted in the field at the beginning of the study.

Model Assessment
The performance of the model was evaluated by comparing field collected/measured data and REMM modeled data for both WTD and groundwater NO 3 -N concentration.
The assessment of the model was accomplished for daily WTD and groundwater NO 3 -N concentration. Assessment of model performance for WTD and groundwater NO 3 -N concentration included both qualitative and quantitative methods. Qualitative methods included comparing graphs of measured and modeled data.
We used the mean absolute error (MAE) to statistically compare the simulated and measured WTDs by quantitatively assessing the goodness-of-fit between simulated and measured WTDs. Equation (3) was used to determine the MAE where W m is measured WTD (cm), W s is simulated WTD (cm), and n is the number of observations. The simulated groundwater NO 3 -N concentrations were compared against the field measured by using the root mean square error (RMSE) and MAE. The RMSE was computed as: where M i is the measured NO 3 -N (mg/L), S i is the simulated NO 3 -N (mg/L), and n is the number of observations. The Willmott's index of agreement (d) between simulated and measured data (WTDs and groundwater NO 3 -N concentration) was also calculated. The value of d is dimensionless and varies between 0 and 1, with an index of 1 corresponding to perfect agreement between simulated and measured data [63]. Equation (5) was used to compute Willmott's index of agreement: where d is the Willmott's index of agreement, M i is the measured data, S i is the simulated data, and M is the mean of the measured data.

REMM Model Calibration, Validation and Sensitivity Analysis
The Riparian Model was developed and evaluated first by testing the hydrologic component (measured WTDs) followed by the nutrient cycling component (measured groundwater NO 3 -N concentrations). The model was calibrated and validated for both hydrology and nutrient cycling in zone 1 and its three layers (similar to the procedures defined in [15,33]. The simulation period varies by site due to different field data collection periods. Simulation period ranged from 1999 to 2005 for sites 1 and 2 [11] and sites 3 and 4 [58] and from 2018 to 2019 for sites 5-8. The calibration and validation dates for each site is shown in tabulated form in Table 3. The model was manually calibrated by changing the values of input parameters one at a time. The range of input parameters were either defined by field/laboratory measurements or obtained from literature or the REMM user's manual.

•
We used field measured daily WTDs in order to calibrate and validate the hydrologic component of REMM.

•
Soil inputs of the upland area were first calibrated, but buffer parameters were kept constant. Soil parameters (soil porosity, field capacity, and wilting point) were then modified within recommended ranges consistent with the soil texture to reduce difference between simulated and measured WTDs.

•
We needed to adjust the soil layer thickness so that REMM generated buffer runoff and AnnAGNPS calibrated runoff, and simulated and measured WTDs were in close agreement.

•
For the improvement of REMM predictions of WTDs, saturated hydraulic conductivities were also adjusted. The REMM permeability class of 2 (saturated hydraulic conductivities ranging from 42-141 µm/s) was used for all the sites. Hydraulic conductivities significantly affected horizontal water movement between riparian zones and the vertical gravity drainage between soil layers [16].

•
The simulated WTDs were also sensitive to deep seepage from the bottom of the third layer (especially when the simulated water table was within layer 3) and were adjusted to improve model predictions of WTDs. Potential deep seep of 0.2 mm/day and 0.1 mm/day were used for all zones for site 5 and site 7, respectively. However, the other six sites had no potential deep seep in the model.

•
After the hydrologic calibration, the litter and soil carbon and nitrogen pools needed to be stabilized. Otherwise REMM might calculate irrational drop in soil organic carbon and associated high N mineralization. Followed by [64], several 35-year simulations (a period selected based on available local historical weather data) were performed by varying percentage of active, slow, and passive pools. Using the initial residue and humus pools, simulations were run and the carbon and nitrogen pools at the end of the period were then used as initial pool values for new simulations. The model was again rerun for another 35-year period which helped to stabilize the carbon and nitrogen pools. After stabilizing these pools, the denitrification rate constant (K d ) was modified to improve the goodness-of-fit between simulated and measured NO 3 -N concentrations in groundwater. The calibrated soil physical buffer inputs, and calibrated K d inputs are available in Supplemental Tables S1 and S2, respectively. • The calibrated model that achieved the best goodness of fit with observed conditions for both WTDs and groundwater NO 3 -N concentrations had previously been saved.
All the calibrated parameters were used without further changes to validate the model for the validation period. Model assessment guidelines defined in Section 2.4 were used to judge goodness of fit for WTDs and groundwater NO 3 -N concentrations in both calibration and validation phases. • Sensitivity analyses were performed to determine the effects of changing a number of key parameters associated with plant growth, nutrient cycling, surface runoff, and soil physical properties for REMM's hydrological and nutrient simulation in Zone 1.
We evaluated the sensitivity value of the most sensitive parameters (soil porosity, field capacity, wilting point) used for WTD estimation. In addition, we also evaluated the sensitivity value of the most sensitive parameters (soil porosity, field capacity, K d ) used for ground water nitrate concentration estimation. Each parameter was changed by +10% and −10% from the values used as the best estimates for each riparian site during calibration as described in [15]. Field capacity was always kept less than soil porosity during the change of these parameters. We utilized the integration of a local method into a global sensitivity method (the random one-factor-at-a-time) design proposed by [65]. This method consists of repetitions of a local method whereby the derivatives are calculated for each parameter by adding a small change to the parameter. The change in model outcome can then be measured by some lumped measure such as total mass export, sum of squares error between modeled and observed values or sum of absolute errors. The following equation has been used to perform sensitivity test for each parameter change-% change (in daily average WTD or Nitrate) = I − C I × 100% where I is the initial calibrated daily average WTD or Nitrate, and C is the changed daily average WTD or Nitrate after parameter change.
The sensitivity analysis and the calibration of REMM model were manually calibrated as in other studies [15,64,66,67].

Water Table Depths Calibration and Validation
Field measured and REMM simulated daily WTD (cm below surface) dropped from field edge Zone 3 to Zone 1. Figure 3 displays the field measured and REMM simulated daily mean WTD (cm below surface) coming from field edge (Zone 3) to stream edge (Zone 1) for site 5, site 6, site 7, and site 8. Due to lack of measured field edge (Zone 3) data, we could not compare the Zone 3 REMM simulated WTDs for site 1-4. Simulated daily Water Table Depths (WTDs) were compared with those measured in the field in Zone 1 (closest to the stream) of the riparian buffer for all the sites. Simulated and measured WTDs in Zone 1 for the calibration and validation periods are shown in Figure 4. Simulated and field measured WTDs in Zone 1 were compared using mean absolute error (MAE) and Willmott's index of agreement (d). In order to measure the variability of MAE and d among the eight riparian sites, the interquartile range (IQR) is shown as the difference between 75th and 25th percentiles, IQR = Q 3 − Q 1 (Table 4). Agriculture 2021, 11, x FOR PEER REVIEW 13 of 25      In general, among the riparian sites, the first Quartile (Q1), the second Quartile or median (Q2), and the third Quartile (Q3) of the MAE for WTDs in Zone 1 was 5.5 cm, 7.5 cm, and 9 cm for the calibration period, respectively. Likewise, for the validation period, the first Quartile (Q1), the second Quartile or median (Q2), and the third Quartile (Q3) of the MAE for WTDs in Zone 1 was 5 cm, 7 cm, and 8.5 cm, respectively. In case of Willmott's index of agreement (d) between daily measured and simulated WTDs, the Q1, the Q2, and the Q3 of the d for WTDs in Zone 1 was 0.12, 0.34, and 0.62 for the calibration period, respectively. Similarly, for the validation period, the Q1, the Q2, and the Q3 of the d for WTDs in Zone 1 was 0.33, 0.64, and 0.75, respectively. The value of d was within the acceptable limit between 0 and 1 for both the calibration and validation periods.

Groundwater NO3-N Concentrations Calibration and Validation
We noticed a decline in field measured and REMM simulated daily groundwater NO3-N concentrations in Zone 1 coming from the field edge of Zone 3. Figure 5 presents the field measured and REMM simulated daily mean NO3-N concentration (mg/L) coming from field edge (Zone 3) to stream edge (Zone 1) for site 5, site 6, site 7, and site 8. Hence, we could obtain an estimate of percentage N-removal from Zone 3 to Zone 1. In particular, for site 5, site 6, site 7, site 8, REMM showed 4.3% increase, no change, 19.7%  In general, among the riparian sites, the first Quartile (Q 1 ), the second Quartile or median (Q 2 ), and the third Quartile (Q 3 ) of the MAE for WTDs in Zone 1 was 5.5 cm, 7.5 cm, and 9 cm for the calibration period, respectively. Likewise, for the validation period, the first Quartile (Q 1 ), the second Quartile or median (Q 2 ), and the third Quartile (Q 3 ) of the MAE for WTDs in Zone 1 was 5 cm, 7 cm, and 8.5 cm, respectively. In case of Willmott's index of agreement (d) between daily measured and simulated WTDs, the Q 1 , the Q 2 , and the Q 3 of the d for WTDs in Zone 1 was 0.12, 0.34, and 0.62 for the calibration period, respectively. Similarly, for the validation period, the Q 1 , the Q 2 , and the Q 3 of the d for WTDs in Zone 1 was 0.33, 0.64, and 0.75, respectively. The value of d was within the acceptable limit between 0 and 1 for both the calibration and validation periods.

Groundwater NO 3 -N Concentrations Calibration and Validation
We noticed a decline in field measured and REMM simulated daily groundwater NO 3 -N concentrations in Zone 1 coming from the field edge of Zone 3. Figure 5 presents the field measured and REMM simulated daily mean NO 3 -N concentration (mg/L) coming from field edge (Zone 3) to stream edge (Zone 1) for site 5, site 6, site 7, and site 8. Hence, we could obtain an estimate of percentage N-removal from Zone 3 to Zone 1. In particular, for site 5, site 6, site 7, site 8, REMM showed 4.3% increase, no change, 19.7% removal, 19.8% removal of N against actual measured 46.7% removal, 28.9% removal, 19.1% removal, and 26.6% removal of N, respectively. Due to lack of measured field edge (Zone 3) nitrate data, we could not compare the Zone 3 REMM simulated nitrate concentration for site 1-4. Yet, it can be mentioned that the REMM simulated daily mean NO 3 -N coming from Zone 3 to Zone 1 were 8.4 mg/L to 0.4 mg/L, 25.7 mg/L to 1.1 mg/L, 8.8 mg/L to 3.9 mg/L, and 10.9 mg/L to 1.4 mg/L for site 1, site 2, site 3, and site 4, respectively. So the approximate simulated percentage N-removal rate becomes 95.2%, 95.7%, 55.7%, and 87.2% for site 1, site 2, site 3, and site 4, respectively. It is noted that the mean value is calculated only for the simulation dates for each site. Simulated and measured daily groundwater NO 3 -N concentrations in Zone 1 during the calibration and validation period are shown in Figure 6 and are compared statistically using the root mean square error (RMSE), the mean absolute error (MAE), and the Willmott's index of agreement (d). In order to measure the variability of RMSE, MAE, and d among the eight riparian sites, the interquartile range (IQR) is shown as the difference between 75th and 25th percentiles, IQR = Q 3 − Q 1 ( Table 5).  55.7%, and 87.2% for site 1, site 2, site 3, and site 4,, respectively. It is noted that the mean value is calculated only for the simulation dates for each site. Simulated and measured daily groundwater NO3-N concentrations in Zone 1 during the calibration and validation period are shown in Figure 6 and are compared statistically using the root mean square error (RMSE), the mean absolute error (MAE), and the Willmott's index of agreement (d).
In order to measure the variability of RMSE, MAE, and d among the eight riparian sites, the interquartile range (IQR) is shown as the difference between 75th and 25th percentiles, IQR = Q3 − Q1 (Table 5).    19.1% removal, and 26.6% removal of N, respectively. Due to lack of measured field edge (Zone 3) nitrate data, we could not compare the Zone 3 REMM simulated nitrate concentration for site 1-4. Yet, it can be mentioned that the REMM simulated daily mean NO3-N coming from Zone 3 to Zone 1 were 8.4 mg/L to 0.4 mg/L, 25.7 mg/L to 1.1 mg/L, 8.8mg/L to 3.9 mg/L, and 10.9 mg/L to 1.4 mg/L for site 1, site 2, site 3, and site 4, respectively. So the approximate simulated percentage N-removal rate becomes 95.2%, 95.7%, 55.7%, and 87.2% for site 1, site 2, site 3, and site 4,, respectively. It is noted that the mean value is calculated only for the simulation dates for each site. Simulated and measured daily groundwater NO3-N concentrations in Zone 1 during the calibration and validation period are shown in Figure 6 and are compared statistically using the root mean square error (RMSE), the mean absolute error (MAE), and the Willmott's index of agreement (d).
In order to measure the variability of RMSE, MAE, and d among the eight riparian sites, the interquartile range (IQR) is shown as the difference between 75th and 25th percentiles, IQR = Q3 − Q1 (Table 5).    Overall, for the calibration period, the Q1, the Q2, and the Q3 of the RMSE for groundwater NO3-N concentrations in Zone 1 were 0.50 mg/L, 0.63 mg/L, and 0.75 mg/L among the riparian sites, respectively, whereas a slightly higher value of RMSE was found for the validation period, including the Q1, the Q2, and the Q3 of the RMSE for NO3-N in Zone 1 as 0.38 mg/L, 0.65 mg/L, and 1.07 mg/L, respectively. For the calibration period, the Q1, the Q2, and the Q3 of the MAE for groundwater NO3-N concentrations in Zone 1 were 0.45 mg/L, 0.56 mg/L, and 0.68 mg/L among the riparian sites, respectively.  Overall, for the calibration period, the Q 1 , the Q 2 , and the Q 3 of the RMSE for groundwater NO 3 -N concentrations in Zone 1 were 0.50 mg/L, 0.63 mg/L, and 0.75 mg/L among the riparian sites, respectively, whereas a slightly higher value of RMSE was found for the validation period, including the Q 1 , the Q 2 , and the Q 3 of the RMSE for NO 3 -N in Zone 1 as 0.38 mg/L, 0.65 mg/L, and 1.07 mg/L, respectively. For the calibration period, the Q 1 , the Q 2 , and the Q 3 of the MAE for groundwater NO 3 -N concentrations in Zone 0.36 mg/L, 0.61 mg/L, and 0.99 mg/L among the riparian sites, respectively. In case of Willmott's index of agreement (d) between daily measured and simulated groundwater NO 3 -N concentrations, the Q 1 , the Q 2 , and the Q 3 of the d for NO 3 -N in Zone 1 were 0.40, 0.44, and 0.60 for calibration period, respectively. Similarly, for the validation period, the first Quartile (Q 1 ), the second Quartile or median (Q 2 ), and the third Quartile (Q 3 ) of the d for WTDs in Zone 1 was 0.17, 0.52, and 0.65, respectively. The value of d was within the acceptable limit between 0 and 1 for both calibration and validation periods.

Parameter Sensitivity Analysis for Water Table Depths and Groundwater NO 3 -N Concentrations Simulation
The sensitivity analysis indicated that simulated predictions of Water Table Depths and groundwater NO 3 -N concentrations could be very sensitive to selected soil physical properties and the denitrification rate constant, respectively. Change in soil porosity caused the greatest change in WTD in all the sites except Site 3 and Site 5 (where change in field capacity caused slightly higher percentage change in WTD than that caused by change in soil porosity; Figure 7). Several of the sites did not exhibit >±10% in response to a 10% change in the input parameters.

Parameter Sensitivity Analysis for Water Table Depths and Groundwater NO3-N Concentrations Simulation
The sensitivity analysis indicated that simulated predictions of Water Table Depths and groundwater NO3-N concentrations could be very sensitive to selected soil physical properties and the denitrification rate constant, respectively. Change in soil porosity caused the greatest change in WTD in all the sites except Site 3 and Site 5 (where change in field capacity caused slightly higher percentage change in WTD than that caused by change in soil porosity; Figure 7). Several of the sites did not exhibit >±10% in response to a 10% change in the input parameters. In Figure 8, the percentage change in daily average groundwater NO3-N concentrations from initial for calibration period is shown to be due to corresponding parameter change. Not only change in soil porosity and field capacity but also change in the denitrification rate constant, Kd greatly affected modeled denitrification and simulated and measured NO3-N concentrations. In Figure 8, the percentage change in daily average groundwater NO 3 -N concentrations from initial for calibration period is shown to be due to corresponding parameter change. Not only change in soil porosity and field capacity but also change in the denitrification rate constant, K d greatly affected modeled denitrification and simulated and measured NO 3 -N concentrations. In Figure 8, the percentage change in daily average groundwater NO3-N concentrations from initial for calibration period is shown to be due to corresponding parameter change. Not only change in soil porosity and field capacity but also change in the denitrification rate constant, Kd greatly affected modeled denitrification and simulated and measured NO3-N concentrations.

Discussion
The MAE between measured and simulated daily WTDs is comparable to the average absolute error of 14 to 36 cm achieved in two previous REMM modeling studies using three to five years of field data from a riparian site in Georgia [12] and two sites in North Carolina [33,34]. The Willmott's index of agreement (d) can be compared to the same statistic found in [33], as 0.72 to 0.92 (yearly scale). In more than 50% of the calibrations and validations, the value of d was greater than 0.5, indicating good agreement with measured data.
The simulated WTDs generally followed seasonal patterns (deeper during drier months, June to November and rising to the surface during wet months, December to May) of measured WTDs for all the sites. As for sites 1-4, REMM simulated WTDs were not under-predicted or over-predicted in a constant manner throughout the entire study period. Conversely and relatively, sites 5-8 showed a particular pattern. At site 5 and site 8, REMM over-predicted the WTDs than the measured from the late spring to summer while under-predicted from late summer to fall season. On the other hand, REMM over-predicted the summer WTDs than the measured but under-predicted the fall WTDs at site 6 and site 7. However, there were some discrepancies between simulated and measured WTDs during some time periods as a result of the spatial variability in rainfall between sites and the considered PRISM climate station, overestimation of ET by model etc. Summer field conditions were modified by irrigation of the uplands at sites 7 and 8. During field sampling, sites 5 and 6 were in unirrigated hay production (with deeper roots systems than turf). As a result, the daily water table depths were significantly deeper at sites 5 and 6 compared to other six sites.
The RMSE between measured and simulated daily groundwater NO3-N concentrations is comparable relatively to the 1.05 to 1.50 mg/L obtained in the field testing study using 5 years of data from a riparian site in North Carolina Coastal Plain [33]. The range of d is also close (0.34 to 0.68) to what was found by [33] from their study on the North Carolina Coastal Plain. The mean absolute error (MAE) for NO3-N concentrations in Zone 1 for calibration and validation periods was reasonably similar to the less than 1 mg/L of absolute error found by [15]. The difference between measured and simulated NO3-N concentrations during some time periods was most likely due to the low frequency of field data collection. This was caused by dry groundwater wells, especially during the summer months, and submergence of one of the groundwater wells in the river during heavy rainfall in November 2018. These constraints resulted in a lower Willmott's index of agreement in case of Site 3, Site 4, and Site 8. In addition, in the case of Site 1, the NO3-N concentration was very low (beyond the detection limit of 0.02 mg/L) for several times, which restricted the frequency of data collection.

Discussion
The MAE between measured and simulated daily WTDs is comparable to the average absolute error of 14 to 36 cm achieved in two previous REMM modeling studies using three to five years of field data from a riparian site in Georgia [12] and two sites in North Carolina [33,34]. The Willmott's index of agreement (d) can be compared to the same statistic found in [33], as 0.72 to 0.92 (yearly scale). In more than 50% of the calibrations and validations, the value of d was greater than 0.5, indicating good agreement with measured data.
The simulated WTDs generally followed seasonal patterns (deeper during drier months, June to November and rising to the surface during wet months, December to May) of measured WTDs for all the sites. As for sites 1-4, REMM simulated WTDs were not under-predicted or over-predicted in a constant manner throughout the entire study period. Conversely and relatively, sites 5-8 showed a particular pattern. At site 5 and site 8, REMM over-predicted the WTDs than the measured from the late spring to summer while under-predicted from late summer to fall season. On the other hand, REMM over-predicted the summer WTDs than the measured but under-predicted the fall WTDs at site 6 and site 7. However, there were some discrepancies between simulated and measured WTDs during some time periods as a result of the spatial variability in rainfall between sites and the considered PRISM climate station, overestimation of ET by model etc. Summer field conditions were modified by irrigation of the uplands at sites 7 and 8. During field sampling, sites 5 and 6 were in unirrigated hay production (with deeper roots systems than turf). As a result, the daily water table depths were significantly deeper at sites 5 and 6 compared to other six sites.
The RMSE between measured and simulated daily groundwater NO 3 -N concentrations is comparable relatively to the 1.05 to 1.50 mg/L obtained in the field testing study using 5 years of data from a riparian site in North Carolina Coastal Plain [33]. The range of d is also close (0.34 to 0.68) to what was found by [33] from their study on the North Carolina Coastal Plain. The mean absolute error (MAE) for NO 3 -N concentrations in Zone 1 for calibration and validation periods was reasonably similar to the less than 1 mg/L of absolute error found by [15]. The difference between measured and simulated NO 3 -N concentrations during some time periods was most likely due to the low frequency of field data collection. This was caused by dry groundwater wells, especially during the summer months, and submergence of one of the groundwater wells in the river during heavy rainfall in November 2018. These constraints resulted in a lower Willmott's index of agreement in case of Site 3, Site 4, and Site 8. In addition, in the case of Site 1, the NO 3 -N concentration was very low (beyond the detection limit of 0.02 mg/L) for several times, which restricted the frequency of data collection.
The sensitivity analyses showed that changes in the volumetric water content between field capacity (FC) and saturation (i.e., soil porosity) was driving the water table and denitrification dynamics. This observation is consistent with the relationships of [68] which show that percentage saturation (the percentage of water-filled pore space, as determined by water content and total porosity) is closely related to denitrification. The lower the difference between FC and soil porosity, the more the changes in water table response to precipitation events [69] and, the larger this difference, the less responsive water tables will be to infiltration-and these two parameters also influence percentage of water-filled pore space. With a high FC and low porosity, the percentage of water-filled pore space at FC might be high enough to regularly generate denitrification.
In general, the results from the sensitivity analyses demonstrated that the percentage change in WTDs (−36% to 25%) was comparatively less than that of groundwater nitrate concentration (−60% to 60%) when responding to the input parameter change (Figures 7  and 8). The percentage change in WTDs was least in site 1 and site 6 (−7% to 7%) while site 2, site 7, and site 8 had the greatest percentage change in WTDs (−36% to 25%). In terms of NO 3 -N concentration, percentage change was the least in site 5 and site 7 (−4.5% to 3.5%) whereas site 2 and site 4 had the most percentage change (−60% to 60%).

Conclusions
In this study, we successfully calibrated and validated the REMM model by coupling upland inputs from a distributed model (AnnAGNPS) with field-measured hydrologic and N data from multiple buffer sites located in a formerly glaciated watershed of Rhode Island. Both the hydrologic and nutrient estimation of REMM showed that it captured well the daily measured WTDs and groundwater NO 3 -N concentrations in Zone 1 and in Zone 3 for the study periods. The sensitivity analyses demonstrated that changes in the volumetric water content between field capacity and saturation (soil porosity) was directing water table and denitrification dynamics.
This modeling study indicated the suitability of REMM to simulate the basic hydrologic and nutrient cycling processes happening in real-world buffers, particularly in glacial geomorphic settings (where REMM applicability has not yet been tested). The use of distributed model AnnAGNPS provided better estimates of upland inputs. The calibrated parameters and model outputs of this study establishes the base for site-specific parameters required to evaluate management and design of riparian buffers effectively in other sites of a similar setting. The site specific and design parameters are soil porosity, field capacity, wilting point, denitrification rate constant, riparian width, vegetation type, etc. The riparian zone hydrologic and nutrient quantification through REMM also contributes to keeping a check on the rates of change occurring in the essential ecological processes (water cycle, biogeochemical or nutrient cycling, the flow of energy, etc.) in ecosystems.