RZ-TRADEOFF: A New Model to Estimate Riparian Water and Air Quality Functions

Abstract

Riparian zones are often used as best management practices due to their ability to remove nitrate (NO₃⁻) from subsurface flow. Research suggests that beyond local biogeochemical controls, the impact of riparian zones on nitrogen removal and other functions, such as phosphorus dynamics and greenhouse gas emissions, largely depends on land-use/land-cover, hydrogeomorphology, and weather. In this study, we therefore present RZ-TRADEOFF, a novel and easily applicable model that connects multiple riparian functions and characteristics (NO₃⁻ and phosphate (PO₄³⁻), concentration and removal in subsurface flow, total phosphorus (TP) removal in overland flow, nitrous oxide (N₂O), methane (CH₄), and carbon dioxide (CO₂) emissions, water table) to landscape hydrogeomorphic characteristics, weather, and land-cover/land-use. RZ-TRADEOFF was developed with data from past studies and digital databases, and validated with data collected from the literature. Three functions (water table, PO₄³⁻ and CO₂) were observed to be significantly influenced by climate/weather, while the others were primarily influenced by hydrogeomorphology and land use. The percent bias and normalized root mean square error respectively were −3.35% and 0.28 for water table, 16.00% and 0.34 for NO₃⁻ concentration, −7.83% and 20.82 for NO₃⁻ removal, 6.64% and 0.35 for PO₄³⁻ concentration, 2.55% and 0.17 for TP removal, 40.33% and 0.23 for N₂O, 72.68% and 0.18 for CH₄, and −34.98% and 0.91 for CO₂. From a management standpoint, RZ-TRADEOFF significantly advances our ability to predict multiple water and air quality riparian functions using easily accessible data over large areas of the landscape due to its scalability.

1. Introduction

Riparian zones (RZ) play a critical role in regulating the movement of nutrients from terrestrial to freshwater environments, and are often considered for use as best management practices (BMPs) to mitigate the impacts of agricultural and urban nitrate (NO₃⁻) loading on water quality [1,2]. Studies have noted the potential for riparian zones to have NO₃⁻ removal efficiencies of 90% or greater [3,4,5]. However, N removal in riparian zones can vary widely depending on hydrology and geomorphology [6]. In riparian buffers across a European climate gradient, the hydraulic gradient and nitrate (NO₃⁻) load were observed to have a greater impact on N removal than differences in weather or vegetation [7]. Similarly, results from a GIS modeling study showed that while vegetation is important, nutrient export is more strongly influenced by hydrology [8]. Therefore, to fully understand the management implications of riparian zones, the influence of hydrogeological attributes and their interactions with weather and land cover must be recognized and considered [9].

While research regarding the use of riparian zones as BMPs has largely focused on N removal, hydrogeomorphology can also influence other aspects of riparian zone biogeochemical functions. Phosphorus (P) loading contributes to water quality degradation, yet research on P retention in riparian zones, both in subsurface and overland flow, is far more limited [10]. In some cases, riparian zones have been noted to retain up to 50% of P input as soluble reactive phosphorus (SRP) in subsurface flow [2], while in others they act as a source of SRP [10]. For P in overland flow, total phosphorus (TP) retention rates have been linked to buffer width [11,12]. Slope has also been seen to influence TP removal, with low gradients that decrease flow retaining as much as 63% of overland flow phosphorus, while steep slopes show minimal or even negative phosphorus retention [13].

Riparian zones are also a small, yet important, natural source of greenhouse gases (GHGs) nitrous oxide (N₂O), carbon dioxide (CO₂), and methane (CH₄) to the atmosphere [14,15,16]. Soil water content and flooding regimes, as determined by hydrogeomorphology and temperature, have been identified as dominant drivers of the production of both N₂O and CH₄ in riparian zones [17,18,19]. Within riparian zones susceptible to flooding, an inverse relationship between soil moisture and N₂O emissions has been observed [20,21]. Meanwhile, manipulation of flooding regimes in constructed riparian wetlands showed that pulsing water regimes and a low static water table led to significantly lower emissions than in those that were continuously flooded [22]. Further, some riparian wetlands and peatlands function as long-term CO₂ sinks [22], while others are a source of CO₂ produced by root respiration, microbial activity, and the oxidation of organic matter [23].

Over the past few decades, several different modeling approaches described below have been used to evaluate the efficacy of riparian zones for water quality management [8,9,12,24,25,26,27,28,29,30]. In most cases, existing models for assessing riparian zone functions are however too simplistic (i.e., large error) or too complex (i.e., too hard to parameterize) for practical management-related use. The majority also focus solely on N removal without considering the impact of riparian zones on P or GHG emissions. These modeling approaches include conceptual models used to develop landscape-scale metrics linking N removal and sink function in riparian zones to permeable soil depth, soil texture, and slope [9,26,31]. Another approach is GIS modeling, used for instance by Baker et al. (2001) to predict subsurface riparian hydrology, which was then linked to nutrient export [8]. In another study, Dosskey and Qiu (2010) determined through GIS models that the optimal location for riparian zones can differ depending on whether placement based on topographic or soil indices is prioritized [30]. Process-based models, such as the Riparian Ecosystem Management Model (REMM), offer an alternative to conceptual or GIS based models [25]. REMM modeling suggests wider riparian buffers are more effective at NO₃⁻ removal and that changes in vegetation cover have limited impact on NO₃⁻ removal [27]. In general, REMM predicted riparian zone water table depths and NO₃⁻ concentrations are reported to be in good agreement with measured values [32], with one case showing a 25% difference between REMM predicted and measured surface flow nitrogen loads [28]. However, REMM has limitations, as it is unbounded and requires a large amount of site-specific data. Finally, a limited number of studies have also used statistical modeling to identify trends and quantify relationships between different riparian zone characteristics and functions. In one statistical meta-analysis, linear and non-linear regression models were developed to examine the relationship of N removal with buffer zone width and land cover type [29]. Results from the models showed that wider buffers are consistently more effective at removing N. Similarly, results from a statistical model developed by Zhang et al. (2010) showed that buffer width accounts for 44% and 35% of the variation in removal efficiency for N and TP in overland flow, respectively [12]. Results also indicate no effect of vegetation cover type [29] and a favorable slope condition of approximately 10% for N removal in overland flow [12].

Based on our current understanding of riparian functions and this short overview of existing riparian modeling approaches, we argue that the next generation of riparian models should be generalizable, easy to use, and address multiple contaminants. To this end, the objectives of this study were to: (1) Develop a set of comprehensive and generalizable statistical models (combined in one model in RZ-TRADEOFF) that can be used to predict not only NO₃⁻ concentration at the field edge and NO₃⁻ removal in riparian subsurface flow, but also phosphate (PO₄³⁻) concentration, GHG emissions (N₂O, CO₂, CH₄), water table depth, and P removal (both in subsurface and overland flow) based on easily accessible hydrogeomorphic, land use/land cover, and weather/climate data; (2) Assess model accuracy through validation with additional data not used for model development, and (3) Create a user-friendly model interface that provides the model equations and allows the user to enter attributes and receive model predictions for a riparian zone of interest. RZ-TRADEOFF was developed for dominant riparian zone types primarily found in the formerly glaciated areas of the North American Northeast and Midwest, and was designed to only require digitally available data (see methods for details) as input data for the purpose of generability and ease of use. It is expected that the model presented here (RZ-TRADEOFF) will help land use managers evaluate water and/or air quality trade-offs in riparian zones and make informed decisions regarding the use of riparian zones for water quality management.

2. Materials and Methods

2.1. Model Development

Individual models included in RZ-TRADEOFF were developed using a database compiled from 27 studies published in peer reviewed journals and one conference presentation paper (

Table 1. Summary of studies used to obtain data for the development of individual models for each riparian function included in the riparian zone (RZ)-TRADEOFF water and air quality model. One site per study used unless indicated otherwise.

Paper	Location †	Type of Data ‡
Abu-Zreig et al. (2003) [11]	Elora, ON	TP
Addy et al. (1999) [33]	Kingston, RI and Charleston, RI (2 sites)	NO3−, water table
Addy et al. (2002) [34]	Kingston, RI and Charleston, RI (2 sites)	NO3−, water table
Daniels and Gilliam (1996) [2]	NC Piedmont (2 sites)	TP
Duval et al. (2007) [35]	Southern ON (2 sites)	NO3−, N removal, water table
Goel et al. (2007) [36]	university of Guelph, ON	TP
Gold et al. (1998) [37]	Kingston, RI	NO3−, water table
Gomez et al. (2016) [38]	Adirondack Region, NY	N2O, CO2, CH4
Jacinthe et al. (2012) [20]	Central IN (3 sites)	N2O, water table
Jacinthe et al. (2015) [39]	Central IN (2 sites)	CO2, CH4
Jacinthe et al. (2017) [40]	Central IN (3 sites)	NO3−, N2O, CH4, CO2, water table
Kellogg et al. (2005) [41]	Pawcatuck River Watershed, RI (4 sites)	NO3−, water table
Kellogg et al. (2008) [42]	Pawcatuck River Watershed, RI (4 sites)	NO3−, water table
Lee et al. (2000) [43]	Bear Creek, IA	TP
Lee et al. (2003) [44]	Bear Creek, IA	TP
Liu et al. (2014) [10]	Indianapolis, IN (2 sites)	NO3−, N removal, PO43−, PO43− removal, water table
Magette et al. (1989) [45]	Chesapeake Bay, MD	TP
Mankin et al. (2007) [46]	Greary County, KS	TP
Nelson et al. (1995) [47]	Kingston, RI	N removal, water table
Simmons et al. (1992) [48]	Kingston, RI (3 sites)	N removal, water table
Vidon and Hill (2004) [49]	Southern ON (8 sites)	NO3−, N2O, N removal, water table
Vidon and Hill (2004) [26]	Southern ON (8 sites)	NO3−, N removal, water table
Vidon and Hill (2004) [50]	Southern ON (8 sites)	water table
Vidon and Serchan (2016) [51]	Adirondack Region, NY	NO3−, N2O, CO2, CH4
Vidon and Serchan (2016) [52]	Adirondack Region, NY	NO3−, PO43−, N2O, CO2, CH4,
Vidon et al. (2014) [53]	Indianapolis, IN	NO3−, N removal, PO43−, PO43− removal, N2O, CO2, CH4, water table
Vidon et al. (2016) [54]	Freeville, NY (3 sites)	NO3−, N removal, PO43−, PO43− removal, N2O, CO2, CH4, water table
Watson et al. (2010) [55]	Central and Southern RI (6 sites)	NO3−, water table

†: IA, Iowa; IN, Indiana; KS, Kansas; MD, Maryland; NC, North Carolina; NY, New York; ON, Ontario; RI, Rhode Island. ‡: CH₄ = methane emissions at the soil–atmosphere interface; CO₂ = carbon dioxide emissions at the soil–atmosphere interface; N removal = nitrate removal in subsurface flow; NO³⁻ = nitrate concentration at the field edge; N₂O = nitrous oxide emissions at the soil–atmosphere interface, PO₄³⁻ = phosphate concentration at the field edge; PO₄³⁻ removal = phosphate removal in subsurface flow; TP = total phosphorus removal in overland flow; water table = water table depth below ground surface.

). The database was populated primarily by original data collected and published by the co-authors of this study. We believe this was the most efficient way to build the database because most other sources, including published studies, generally only provide aggregated values and/or limited hydrogeomorphic information about study sites [26]. Altogether, the original dataset spanned 88 data years and 30 sample sites representing the dominant riparian zone types found in Southern Ontario (Canada) and the US Northeast and Midwest regions (Table 1,

Table 2. Summary of data used for the development of individual models in RZ-TRADEOFF.

Variable †	Min. Value	Max. Value	Range	# of Observations (# of Sites)
WT depth (cm)	2.8	186.5	183.5	307 (16)
NO3− (mg L−1)	0.005	20.7	20.7	276 (23)
NO3− Removal %	−100	100		203 (14)
PO43− (mg L−1)	0.005	0.47	0.47	91 (6)
TP Removal %	22	99		20 (8)
PO43− Removal %	−100	92		88 (6)
N2O (mg N m−2 d−1)	0.01	7.01	7.00	169 (8)
CO2 (g C m−2 d−1)	−0.51	11.17	11.68	170 (8)
CH4 (mg C m−2 d−1)	−6.89	2612.27	2619.16	171 (8)

†: WT = water table; NO³⁻ = nitrate concentration at the field edge; NO₃⁻ removal % = nitrate removal in subsurface flow (%); PO₄³⁻ = phosphate concentration at the field edge; TP = total phosphorus removal in overland flow; PO₄³⁻ removal % = phosphate removal in subsurface flow (%); N₂O = nitrous oxide emissions at the soil–atmosphere interface; CO₂ = carbon dioxide emissions at the soil–atmosphere interface; CH₄ = methane emissions at the soil–atmosphere interface.

). The data published by the co-authors of this study were supplemented with additional data retrieved from six different peer reviewed studies and a conference presentation paper reporting values for total phosphorus (TP) removal in overland flow. These studies were conducted in a similar geographic region as the original data, with sites in Iowa, Kansas, Ontario, and Maryland, except for one site located in North Carolina (Table 1).

The database contained measurements for 6 riparian biogeochemical functions (CO₂, CH₄, and N₂O fluxes at the soil–atmosphere interface, NO₃⁻ and PO₄³⁻ percent removal in subsurface flow, and TP removal in overland flow), water table depth, and NO₃⁻ and PO₄³⁻ concentrations at the field edge, which represent input from the upland area into the riparian buffer. Nitrate and PO₄³⁻ removal were only considered in subsurface flow since research suggests subsurface flow is the primary transport mechanism for both PO₄³⁻ and NO₃⁻ in riparian zones [6,26,56,57,58]. The 14 attribute characteristics in the database were buffer width, confining layer at each site and immediate upland area, dominant soil type at each site and immediate upland area, surficial geology (i.e., glacial till, outwash, alluvium), seeps, topography (topography type at the site, overall slope, and edge of field slope), land cover at each site and in the immediate upland area, 30-year mean normal temperature, and 30-year mean normal precipitation. The database also included 13 weather variables (total 1, 7, 14, 30 day antecedent precipitation, mean 7, 14, 30 day antecedent precipitation, maximum and minimum 1 day temperature, mean 1, 7, 14, and 30 day antecedent temperature) calculated by RZ-TRADEOFF from daily temperature and precipitation values entered by the user. Data collected by the co-authors were entered in both summary form (averaged yearly values) and original form (daily values), while the additional TP data were only recorded in summary form as provided by the literature. The weather data were collected from the nearest National Weather Service cooperative observer program station, the National Climate Data Center (NCDC) online database, or the Government of Canada Historic Climate database. Categorical attribute variables (soil type, land cover, topography type, and geology) were coded with numeric values for the purpose of using them as dummy categorical variables in the model. The soil type variables were coded in order of increasing hydraulic conductivity, with 1 for peat and 24 for coarse sand gravel. Upland land cover was coded from 1 for forested to 5 for fertilized agricultural land to reflect potential impacts on water quality with respect to nitrogen, and the site land cover was coded from 1 for herbaceous to 20 for softwood forests to reflect the typical vegetation gradient observed in riparian zones. The topography variable was coded from 1 for flat to 3 for depressed topography. If not directly known by the user, the data described above can be obtained for use in the model from web sources including the United States Department of Agriculture Natural Resources Conservation Service (USDA-NRCS) web soil survey (https://websoilsurvey.sc.egov.usda.gov) for soil/surficial geology, the NCDC database (https://www.ncdc.noaa.gov/cdo-web/) for weather data, various digital elevation models (DEM) (e.g., National Digital Elevation dataset or Lidar derived DEMs where available) for topography, the National Land Use/Land Cover (LULC) dataset for LULC information, or satellite photography for individual sites when needed (e.g., Google Earth).

Individual models were developed from data subsets for each function using the PAST statistical software package [59]. Statistical methods such as principal component analysis (PCA) and correlation analysis were used to reduce data dimensionality and identify influential attributes and weather variables. Specifically, PCA was used to determine the combination of variables that explained dataset variability (and implicitly those that don’t), while correlation analyses were used to identify and analyze the strength of linear relationships between riparian zone attributes and functions, and eliminate independent variables that were highly correlated (see below) from the list of candidate variables for model development. Linear correlations with a p-value of less than 0.05 were considered significant. To avoid any instances of multicollinearity, correlations between attributes with a correlation coefficient (R²) of 0.8 or higher were considered strongly correlated and were not used together in the models. Together, PCA and correlation analysis therefore allowed us to reduce the number of candidate variables for model development and start model development with only variables that have the potential to affect variability in the dataset without being auto-correlated. Subsequently, linear regression models were built using the influential attributes from the PCA and correlation analysis, with backwards elimination used to identify significant predictors (p < 0.05) and estimate coefficients. The final models were selected and assessed based on goodness of fit represented by adjusted R² (_adjR²) values to reduce over-fitting the model as would occur if simple R² values were used [60].

Bivariate models were also considered in cases when initial data analysis (PCA, correlation) showed a strong relationship between only one attribute or weather variable, and a riparian zone function or characteristics. Depending on the type of relationship observed, bivariate linear, Michaelis–Menten, exponential, and logarithmic models were tested with the best fit determined by likelihood using Akaike’s Information Criterion (AIC) values. If no significant hydrogeomorphic, land use/land cover, or weather predictors were identified, the mean was calculated from the development data and reported as a constant value.

2.2. Model Validation

The predictive ability of the RZ-TRADEOFF model was validated using additional independent data from the published literature (See

Table 3. Summary of studies used to obtain data for the validation of the individual models included in the RZ-TRADEOFF water and air quality model.

Source	Location †	Type of Data ‡
Archibald (2010) [61]	Hartford, NY	PO43−, water table
Clausen et al. (2000) [62]	Muddy Brook, CT	N% Removal
Daniels and Gilliam (1996) [2]	NC Piedmont	N Removal
DeSimone et al. (2010) [63]	Beverly Swamp, ON	N2O
Dillaha et al. (1988) [13]	Blacksburg, VA	TP
Gerwig et al. (2001) [64]	Dupin County, NC	PO43−
Gleason et al. (2009) [65]	Stutman County, ND	CO2
Groffman et al. (2002) [66]	Gwynns Falls watershed, MD	water table
Groffman et al. (2009) [67]	Gwynns Falls watershed, MD	N2O
Hefting and Bobbnik (2003) [17]	Eastern Netherlands	water table
Hernandez and Mitch (2006) [68]	Columbus, OH	N2O
Jaynes and Isenhart (2014) [69]	Bear Creek, Story County, IA	water table
Jordan et al. (1993) [4]	Southeast Creek, MD	NO3−, PO43−, N Removal
Kim et al. (2009) [70]	Bear Creek, Story County, IA	NO3−, N2O
Kim et al. (2010) [71]	Bear Creek, Story County, IA	CH4
Klocker et al. (2009) [72]	Mineback Run, MD	N Removal
Kronvang et al. (2005) [73]	Southern Norway	TP
Lee et al. (2003) [44]	Bear Creek, Story County, IA	N Removal
Martin et al. (1999) [5]	Carol Creek, ON	N Removal
Morse et al. (2012) [16]	Albermarle Peninsula, NC	N2O, CO2, CH4, water table
Peterjohn and Correll (1984) [3]	Rhode River (near Annapolis, MD)	N% Removal
Schilling et al. (2004) [74]	Walnut Creek, Jasper County, IA	water table
Schilling et al. (2006) [75]	Walnut Creek, Jasper County, IA	NO3−
Schmitt et al. (1999) [76]	Mead, NE	TP
Schoonover and Willard (2003) [77]	Cypress Creek Watershed, IL	NO3−, N Removal
Soosaar et al. (2010) [15]	Tartu County, Estonia	N2O, CO2, CH4
Syverson and Borch. (2005) [78]	South Eastern Norway	TP
Teiter and Mander (2005) [79]	Tartu County, Estonia	CO2
Tomer et al. (2007) [80]	Traynor, IA	PO43−, water table
Tufekcioglu et al. (1999) [81]	Bear Creek, Story County, IA	CO2
Ullah and Moore (2011) [82]	Montreal, QC	N2O, CH4
Woltemade and Woodward (2008) [83]	Shippensburg, PA	NO3−, N% Removal
Young and Briggs (2008) [84]	Onondaga County, NY	PO43−, water table

†: CT = Connecticut; IA = Iowa; IL = Illinois; MD = Maryland; NC = North Carolina; ND = North Dakota; NE = Nebraska; NY = New York; OH = Ohio; ON = Ontario; PA = Pennsylvania; QC = Quebec; VA = Virginia; VT = Vermont. ‡: CH₄ = methane emissions at the soil–atmosphere interface; CO₂ = carbon dioxide emissions at the soil–atmosphere interface; N removal = percent nitrate removal in subsurface flow; NO³⁻ = nitrate concentration at the field edge; N₂O = nitrous oxide emissions at the soil–atmosphere interface, PO₄³⁻ = phosphate concentration at the field edge; PO₄³⁻ removal = phosphate removal in subsurface flow; TP = total phosphorus removal in overland flow; water table = water table depth below ground surface.

and

Table 4. Summary of data used to validate the individual models in RZ-TRADEOFF.

Variable †	Summary of validation data ‡
	Min. Value	Max. Value	Range	# of Observations (# of Sites)
Water table depth (cm)	30	349	319	22 (5)
NO3− (mg L−1)	2	13	11.00	25 (5)
NO3− Removal %	40	99	59.00	11 (9)
PO43− (mg L−1)	0.005	1.6	1.60	28 (5)
TP Removal %	55	90	35.00	9 (4)
PO43− Removal %	n/a	n/a	n/a	0
N2O (mg N m−2 d−1)	0.01	4.80	4.79	23 (5)
CO2 (g C m−2 d−1)	0.48	19.20	18.72	35 (4)
CH4 No Peat (mg C m−2 d−1)	−0.5	26.67	27.17	24 (4)
CH4 Peat (mg C m−2 d−1)	n/a	n/a	n/a	0

†: CH₄ = methane emissions at the soil–atmosphere interface; CO₂ = carbon dioxide emissions at the soil–atmosphere interface; NO³⁻ removal = percent nitrate removal in subsurface flow; NO³⁻ = nitrate concentration at the field edge; N₂O = nitrous oxide emissions at the soil–atmosphere interface, PO₄³⁻ = phosphate concentration at the field edge; PO₄³⁻ removal = phosphate removal in subsurface flow; TP = total phosphorus removal in overland flow; water table = water table depth below ground surface. ‡: n/a and zero values indicate that although model estimates were compared with known literature values in the discussion section, no independent validation sites directly applicable to the RZ-TRADEOFF model were found in the literature.

for a complete list of validation data). Remaining consistent with the development data, most of the validation sites were located in the Midwest and Northeast regions of the United States. However, since the availability of data for some riparian zone functions reported in a format that was appropriate for the models was limited, the geographic region was expanded to include a few sites in Northeastern Europe and one site in North Carolina. We argue that the European sites were appropriate to use for validation as they were located in previously glaciated settings and currently have a similar climate and geography than that of the area of interest in North America. If all the characteristics necessary for a model were not reported in the study, the information was retrieved from the USDA-NRCS Web Soil Survey for attribute data (https://websoilsurvey.sc.egov.usda.gov), and either the NCDC (https://www.ncdc.noaa.gov/cdo-web/) or the Government of Canada Historic Climate database (http://climate.weather.gc.ca/) for climate/weather data using the weather stations closest to the location of the study site.

Percent bias (PBIAS) between the model output and the actual values of the validation data was calculated and used as a measure of model accuracy [85,86,87]. Root mean square error (RMSE), which accounts for the range of actual values was used as an additional accuracy measure [88]. The normalized RMSE was also calculated by dividing the RMSE by the range to allow for a comparison of error between models [89]. Mean absolute percent error (MAPE) was also calculated and used to represent model precision. Nash–Sutcliffe coefficients were not used because they require multiple observations from one site, whereas the studies used for validation commonly reported aggregated data [86]. To further increase accuracy, the models were calibrated based on PBIAS. The average PBIAS of a randomly selected two-thirds of the validation observations for a model was calculated and rounded to the nearest 10%. Model outputs were then increased (or decreased if negative) by the calculated percentage for the entire validation dataset. If the average PBIAS was improved by the calibration, the percent change was kept as a calibration factor for the final model [90].

2.3. RZ-TRADEOFF User Interface Development

A user-friendly interface for RZ-TRADEOFF was created in Excel [91] for the purpose of allowing users to easily apply RZ-TRADEOFF. Model information, the RZ-TRADEOFF model, and the user interface can be accessed at https://dataverse.harvard.edu/dataverse/esf. The RZ-TRADEOFF Excel interface allows the user to input attribute and weather variables and calculates the predicted values using the individual model equations. In addition, RZ-TRADEOFF provides supplementary calculations based on the model output including minimum and maximum model value for uncertainty analysis by the user, GHG output converted to CO₂ equivalents, groundwater flux, NO₃⁻ flux, and the mass of NO₃⁻ removed. Specifically, RZ-TRADEOFF is programmed to calculate the groundwater flux between the field edge and the stream using the one-dimensional form of Darcy’s law: q = Ks(δh/lδ)Z per unit length of a stream bank. In this equation, q is the water flux (L day⁻¹) per meter of stream length. Ks is the saturated soil hydraulic conductivity (m day⁻¹) and is entered by the user either based on field measurement (if available) or based on soil texture class [92]. In addition, δh/δl is the topographic gradient as a proxy for water table gradient [93], and Z is the effective flow depth (confining layer depth minus water table depth) in meters. The flux is expressed over 1 meter of stream length and a unit conversion factor. The NO₃⁻-N flux entering the riparian zone is then calculated as the groundwater flux multiplied by the model predicted NO₃⁻ concentration at the field edge. The mass of NO₃⁻ -N removed across the riparian zone (N sink) is calculated as the NO₃⁻ flux entering the riparian zone multiplied by the model predicted percent NO₃⁻ removal in subsurface flow. The N sink values calculated by the model are assigned a rank of small (<0.1 g N d⁻¹ m⁻¹ stream length), medium (0.1–1 g N d⁻¹ m⁻¹ stream length), or large (>1 g N d⁻¹ m⁻¹ stream length) based on known variability in N sink values from the literature [10,26]. RZ-TRADEOFF also calculates all GHG flux values in CO₂ equivalents using the Intergovernmental Panel on Climate Change (IPCC) methodology and CO₂, N₂O, and CH₄ global warming potentials of 1, 298, and 25, respectively [94].

3. Results

3.1. RZ-TRADEOFF Model Structure and Equations

Hydrogeomorphic attributes frequently observed to have a significant influence on the modeled variables were soil type, slope, land cover type, and depth to confining layer (

Table 5. Description of individual models for each riparian function represented in the RZ-TRADEOFF water and air quality model.

Variable †	Model Type	Model Equation ‡	AdjR2 *
Water table depth (cm)	Multivariate Linear	143.24 + (−5.476 × LC Site) + (5.731 × CL Site) + (−10.267 × Topo Site) + (1.593 × Mean30dt) + (−2.762 × Mean14dp) + (−3.369 × Soil Site)	0.65
NO3− (mg L−1)	Multivariate Linear	(−3.464) + (−0.005 × Julian Date) + (0.236 × CL Site) + (1.160 × Edge Slope) + (0.187 × Soil Upland) + (0.522 × LC Upland)	0.38
NO3− Removal %	Multivariate Linear	1.3 × ((−330) + (0.896 × RZ width) + (−5.379 × CL site) + (2.804 × Edge Slope) + (14.306 × Soil Site) + (−11.205 × Soil Upland) + (59.934 × LC Upland))	0.45
PO43− (mg L−1)	Multivariate Linear	0.053 + (0.0047 × CL Site) + (−0.0025 × Mean 30dt)	0.29
TP Removal %	Logistic	1.03 × (90.056/(1 + 3.7939 × exp(−0.482 × RZ Width)))	n/a
PO43− Removal %	Constant Value	Constant 1.5%	0
N2O (mg N m−2 d−1)	Univariate Linear	0.388 + (0.024 × CL-Upland)	0.20
CO2 (g C m−2 d−1)	Multivariate Linear	1.8 × [(−0.844) + (0.180 × Nor Ann T) + (0.112 × Mean30dt)]	0.29
CH4 non-peat (mg C m−2 d−1)	Michaelis-Menten	(−0.240 × Soil Site)/(−5.045 + Soil Site)	n/a
CH4 peat (mg C m−2 d−1)	Constant Value	Constant 621.5	0

*: Significant at the 0.05 probability level; n/a = not applicable; multivariate linear refers to linear models with multiple independent predictors and one dependent (response) variable. †: CH₄ = methane emissions at the soil–atmosphere interface; CO₂ = carbon dioxide emissions at the soil–atmosphere interface; NO³⁻ removal = percent nitrate removal in subsurface flow; NO³⁻ = nitrate concentration at the field edge; N₂O = nitrous oxide emissions at the soil–atmosphere interface, PO₄³⁻ = phosphate concentration at the field edge; PO₄³⁻ removal = phosphate removal in subsurface flow; TP = total phosphorus removal in overland flow; water table = water table depth below ground surface. ‡: n/a and zero values indicate that although model estimates were compared with known literature values in the discussion section, no independent validation sites directly applicable to the RZ-TRADEOFF model were found in the literature. ‡: CL site = depth to the confining layer (CL) below ground in the riparian zone (m); CL upland = depth to the confining layer below ground in the immediate upland area; edge slope = percent slope at the field edge; LC site = land cover in the riparian zone; LC upland = land cover type in the immediate upland area; mean14dp = mean 14-day antecedent precipitation (mm); mean30td = mean 30-day antecedent temperature (°C); Nor Ann T = normal annual temperature in the location of the riparian zone (°C); RZ width = width of the riparian zone (m); soil site = dominant soil type in the riparian zone; soil upland = dominant soil type in the immediate upland area; topo site = topography type at the riparian zone. All model variables are further defined in the Excel based RZ-TRADEOFF model.

). Four of the models, CO₂, PO₄³⁻, and NO₃⁻ concentration at the field edge, and water table depth were influenced by temporal variables (Julian day, temperature, or precipitation), and thus will produce output values that vary by date. For the other riparian zone functions, temporal variables were not observed to be significant and are predicted with models using only static geomorphic or land use/land cover characteristics. Therefore, the model outputs are not dependent on time of year/date.

The water table depth model developed for RZ-TRADEOFF is a linear model including two temporal and four hydrogeomorphic variables (_AdjR² = 0.65). The output will vary by date, since two temporal variables (mean 30-day antecedent temperature and mean 14-day antecedent precipitation) were significant in the model. The land use/land cover and hydrogeomorphic characteristics in the water table depth model were land cover, soil type, confining layer depth at the site, and topography type (Table 5).

The N related components of RZ-TRADEOFF included linear models for NO₃⁻ removal in subsurface flow and NO₃⁻ concentration at the field edge as well as calculations for NO₃⁻ flux and mass of N removed (N sink strength). The six predictor-variables in the NO₃⁻ removal model in the subsurface accounted for 43.5% of the variability in the model (Table 5). This model was largely influenced by land cover in the upland, soil type, and slope. The four significant land use/land cover and hydrogeomorphic variables in the NO₃⁻ concentration model (_AdjR² = 0.38) were confining layer depth, field edge slope, soil type at the site, and land cover in the upland (Table 5). While the NO₃⁻ concentration model did not include any weather variables, Julian day was observed to have a significant effect, thus indicating that the NO₃⁻ concentration will vary depending on the time of year.

Three separate P related models were developed for TP removal in overland flow, PO₄³⁻ concentration at the field edge, and PO₄³⁻ removal in subsurface flow. TP removal was best estimated by a logistic model, where the removal percentage increases with increasing riparian width (Table 5). The PO₄³⁻ concentration at the field edge model had an _AdjR² value of 0.29 and was driven by confining layer depth at the site and mean 30-day antecedent temperature (Table 5). No significant model (p < 0.05) could be developed from the database for PO₄³⁻ removal in subsurface flow, which means that the variables used in our database cannot explain PO₄³⁻ removal and that the mean of the observed values is therefore the best predictor based on the data available (Table 5).

Greenhouse gas models were developed for N₂O, CH₄, and CO₂ emissions at the soil atmosphere interface in riparian zones. Nitrogen oxide emissions were best predicted by a bivariate linear model, as the confining layer depth in the upland was observed to be the sole variable with significant influence (Table 5). Methane was best predicted by two separate models for peat sites and non-peat sites. A Michaelis–Menten model driven by soil type at the site was used for CH₄ emissions in non-peat sites while emissions in peat sites were represented by a constant value (Table 5). The CO₂ emission model is a linear model driven by two weather/climate variables, with emissions increasing with both higher normal annual temperature and mean 30-day antecedent temperature (Table 5).

3.2. Model Validation

Data compiled from a total of 33 independent published studies (i.e., not used for model development; Table 3 and Table 4) were used to validate each of the models included in RZ-TRADEOFF. Four different error values (PBIAS, mean absolute percent error (MAPE), RMSE, and normalized RMSE) were calculated to assess and compare the models (

Table 6. Error values of the individual models included in RZ-TRADEOFF.

Variable †	Model Error ‡§
	PBIAS	MAPE	RMSE	Normalized RMSE ¶
Water table depth (cm)	−3.35	39.56	90.30	0.28
NO3− (mg L−1)	16.00	53.18	3.70	0.34
NO3− Removal %	−7.83	18.73	20.82	0.35
PO43− (mg L−1)	6.46	93.14	0.70	0.44
TP Removal %	2.55	6.79	6.00	0.17
PO43− Removal %	n/a	n/a	n/a	n/a
N2O (mg N m−2 d−1)	40.33	100.10	1.10	0.23
CO2 (g C m−2 d−1)	−34.98	75.48	17.0	0.91
CH4 No Peat (mg C m−2 d−1)	72.68	294.4	4.80	0.18
CH4 Peat (mg C m−2 d−1)	n/a	n/a	n/a	n/a

†: See Table 5 for acronym definition in the “variable” column. ‡: PBIAS = percent bias; MAPE = mean absolute percent error; RMSE = root mean square error. §: n/a = see definition in Table 5. ¶: normalized RMSE = root mean square normalized by the range; n/a = not applicable.

). The PBIAS calculated from the validation observations ranged from −34.98% for the CO₂ model to 72.68% for the CH₄ (Table 6). The MAPE ranged from 6.75% for TP removal to 100.10% for N₂O (Table 6). The models with the lowest PBIAS and MAPE values overall, and therefore the most accurate and precise, were TP removal (PBIAS 2.54%; MAPE 6.75%) and NO₃⁻ removal (PBIAS −7.83%; MAPE 18.73%). The TP removal and CH₄ models had the lowest relative error as represented by the normalized RMSE, with values of 0.17 and 0.18. All three GHG models had high PBIAS and MAPE values compared to the other models (Table 6). However, RMSE values were low for both N₂O and CH₄ models at 0.23 and 0.18 respectively, while high for CO₂ at 0.91 (Table 6).

4. Discussion

4.1. To What Extent Does RZ-TRADEOFF Contribute to Improving Riparian Zone Water Table and Groundwater Flow Modeling?

RZ-TRADEOFF includes the first empirical model for estimating water table depth in riparian zones. Additionally, RZ-TRADEOFF uses the water table depth estimates and Darcy’s law to calculate the groundwater flux in the riparian zone in a field-to-stream direction. By only using landscape scale variables, the RZ-TRADEOFF water table model presents an alternative complement to REMM, which can simulate water table depth, surface runoff, and hydrologic budgets, but requires finer-scale variables such as infiltration and soil moisture. In a study by Inamdar et al. (1999) [95], the mean absolute residual error between REMM simulated water table depths and actual values was 0.27 m and 0.14 m for two separate sites. The mean absolute residual error for the RZ-TRADEOFF water table model is somewhat less favorable at 0.64 m (data not shown in this format in Table 6). However, this represents the error of the model when applied to several sites, not one individual site as in the Inamdar et al. study [95]. These numbers are therefore not directly comparable. Furthermore, the low PBIAS of −3.35% indicates that the RZ-TRADEOFF predicted values are in good agreement with the observed values. Ultimately, because the water table model in RZ-TRADEOFF has a good model performance (i.e., low PBIAS = −3.35%), RZ-TRADEOFF is likely to be useful to estimate water table depth values for multiple dates and sites across broad areas, even though the precision on a site-per-site basis is only moderately good (i.e., MAPE = 39.56%).

4.2. To What Extent Does RZ-TRADEOFF Advance Our Ability to Predict N Related Riparian Functions?

The RZ-TRADEOFF models for NO₃⁻ concentration at the field edge and percent NO₃⁻ removal in the subsurface and the calculations for the N sink function using Darcy’s law together improve on previous modeling efforts with regards to riparian N related functions. For instance, the meta-analysis conducted by Mayer et al. (2007) presented several non-linear bivariate regression models that assess the individual influence of width, vegetation type, and flow path on NO₃⁻ removal [29]. RZ-TRADEOFF builds on these models by incorporating buffer width and land cover type, as well as additional variables such as soil type and slope, into one linear model for predicting NO₃⁻ removal in subsurface flow. This linear NO₃⁻ removal model has an _AdjR² value of 0.45, improving on the _AdjR² values of the Mayer et al. (2007) models ranging from 0.00–0.39. RZ-TRADEOFF also compares favorably and offers some management advantages to the data-intensive REMM model [29]. For NO₃⁻ concentration, Tilak et al. (2014) found that the RMSE for REMM predicted values ranged from 1.04–5.92 mg L⁻¹ [32]. The RZ-TRADEOFF RMSE value for NO₃⁻ concentration at the field edge was consistent with the midpoint of this range at 3.7 mg L⁻¹, but RZ-TRADEOFF is less data intensive than REMM. Inamdar et al. (1999) found that REMM simulations of surface N flow load range from one standard deviation to a full order of magnitude difference from the actual values [95]. Additional studies have conducted sensitivity analyses [27] as well as tested the ability of REMM to simulate total Kjeldahl nitrogen load during storm events [28]. However, there are no studies that we know of that assess the accuracy and precision of REMM in simulating subsurface NO₃⁻ removal. Although comparison data are scarce, comparison with existing models suggests that the N components of the RZ-TRADEOFF model are an improvement over existing model offerings because of better model performance and/or by generating a broader range of N related outputs. In addition, the accuracy and the RMSE of the model are good for NO₃⁻ removal and NO₃⁻ concentration at the field edge, suggesting again that one of the primary benefits of RZ-TRADEOFF for estimating N related functions might be for watershed management for many sites (easy-to-use, good accuracy).

4.3. How Does RZ-TRADEOFF Improve Our Ability to Estimate P Dynamics in Riparian Zones?

RZ-TRADEOFF is also the first empirical model with the ability to predict PO₄³⁻ concentration at the field edge, TP removal in overland flow, and PO₄³⁻ removal in subsurface flow. As there are no models for PO₄³⁻ concentration or PO₄³⁻ removal, there is no basis to compare the results of RZ-TRADEOFF with respect to these functions. Overall, however, the TP model is both highly accurate and precise, with PBIAS and MAPE <10% and a low normalized RMSE value. Based on the PBIAS and MAPE, the PO₄³⁻ concentration model has high accuracy but low precision. Indeed, subsurface PO₄³⁻ removal is reported as constant mean value. This does not imply that the PO₄³⁻ removal efficiency of all riparian zones is identical, but instead that the mean is a better estimate than a model prediction. This is consistent with literature findings that show a wide variability in riparian zone PO₄³⁻ removal capacity. In Liu et al. (2014), for instance, one riparian site was observed to be a source of SRP with increased outputs ranging from 1100% to 44% above field edge values, while another had removal efficiency between 32–63% [10]. Other cases show riparian zones increasing SRP as much as 3650% [4] and removals of up to 50% [95]. With respect to TP removal in overland flow, Zhang et al. (2010) modeled the relationship between buffer width and TP removal and predicted a maximum removal efficiency of 89.5% within the buffer width range of 20–25 m [12]. RZ-TRADEOFF predicts a TP removal of 93% for riparian zones that are 20–25 m wide. Both RZ-TRADEOFF and the model developed by Zhang et al. (2010) show the importance of buffer width in predicting the removal efficiency of TP in overland flow and suggest that relatively short buffers (i.e., 20–25 m) are expected to have removal efficiencies near 90% [12]. Overall, the TP removal in the overland flow model has very good model performance, suggesting it is highly effective at providing estimates at both the single-site and watershed scale. The PO₄³⁻ concentration model however, is more accurate than precise, and may be more useful when applied at the watershed scale.

4.4. What is the Extent of RZ-TRADEOFF’s Ability to Estimate the Air Quality Impacts of Riparian Zones?

RZ-TRADEOFF provides a novel method of estimating N₂O, CH₄, and CO₂ fluxes from riparian zones to the atmosphere. Due to unique hydrologic and biogeochemical conditions, riparian zones can be potential hotspots of soil-to-atmosphere emissions of natural GHGs [39]. Because of the prominence of climate change as an environmental issue, the potential increase in natural GHG emissions are an important factor to account for when considering the use of riparian buffers for water quality management. When used together with the N and P models, predictions from the RZ-TRADEOFF GHG models allow for the assessment of trade-offs with regards to the air quality and water quality impacts of riparian zones. All three GHG models are less accurate and less precise than those for other variables in RZ-TRADEOFF. The CO₂ model has the most favorable PBIAS and MAPE, however the normalized RMSE is quite high in comparison to the range, indicating a relatively high residual error. Conversely, while the CH₄ Michaelis–Menten model for non-peat sites has the highest PBIAS and MAPE, the normalized RMSE is low, suggesting a small residual error due to the wide range of variability of CH₄ fluxes (4 orders of magnitudes). A constant is used for CH₄ fluxes in peat sites because the values in the development dataset were much higher (>600 mg m⁻² d⁻¹) than that of the non-peat sites (−0.5–26.67 mg m⁻² d⁻¹). Therefore, including the peat sites in the Michaelis–Menten model would reduce its sensitivity for the other sites. Some studies report mean peatland CH₄ emission values lower than 50 mg m⁻² d⁻¹ [96,97], while others are between 100–200 mg m⁻² d⁻¹ [98,99]. A synthesis by Turetsky et al. (2014) meanwhile reported mean fluxes across several peatlands ranging from 93.0–96.3 mg m⁻² d⁻¹ with maximum values from 1328.2–5722.3 mg m⁻² d⁻¹ [100]. This range and variability of fluxes is not captured by the development dataset, therefore the constant value reported here for CH₄ emissions may not present a wholly accurate representation of CH₄ emission in peat dominated riparian zones. Though the model performance of the RZ-TRADEOFF GHG models are highly varied and likely unsatisfactory in the case of the CH₄ peat model, the GHG components of the RZ-TRADEOFF model for N₂O, CO₂, and CH₄ non-peat sites nevertheless provide an important advancement in our ability to predict the air quality impacts of riparian zones. Further, the residual error and accuracy based on the RMSE is low for the N₂O and CH₄ models, indicating that these two models may provide useful estimates at the watershed scale nonetheless.

Of note is the positive correlation between N₂O emissions and confining layer depth in our model (see Table 5). Several studies have noted that groundwater fluctuations and flooding regimes influence N₂O emissions strongly in riparian zones [20,101]. Considering the known influence of upland confining layer depth on riparian water table dynamics [51,102], it is therefore logical that from a statistical standpoint, upland confining layer depth is important at controlling riparian N₂O emissions.

4.5. What Are the Implications and Applications of RZ-TRADEOFF?

The RZ-TRADEOFF model uses easy to obtain land use/land cover, hydrogeomorphic characteristics, and weather data to predict multiple key riparian zone functions and characteristics that can be used to assess trade-offs both within and between riparian zones. RZ-TRADEOFF builds off the wide body of literature documenting the importance of land use, riparian hydrogeomorphology, and weather at regulating riparian functions by quantifying the effects of these variables on key riparian functions. RZ-TRADEOFF also expands on previous models, such as those developed by Mayer et al. (2007) to incorporate models for P, GHGs, and water table depth in addition to NO₃⁻ removal [29]. The comprehensiveness of RZ-TRADEOFF advances our modeling ability by simultaneously estimating the efficacy of a specific riparian zone at, for instance, removing agricultural N runoff as well as any potential negative impacts such as GHG emissions. The ease with which the model can be applied from using digitally available data also allows for the application of the model to only one site but to multiple sites at the watershed scale.

Further, model output from RZ-TRADEOFF for theoretical “design” buffers can also be compared to determine which characteristics would be the most appropriate based on site characteristics and individual management goals so management guidelines for specific riparian buffer types can be developed [103]. For instance, in areas where implementing buffers is an issue due to their large footprint, RZ-TRADEOFF can identify the minimum width required to meet a desired NO₃⁻ removal efficiency. In regions such as the Chesapeake Bay area where buffers need to be implemented at a large scale, the model can be used to select the optimal locations for placement. If limiting air quality impacts is a management priority, the model output can assist in determining the characteristics and size required to keep GHG below a set limit while also maintaining any desired water quality benefits. In this way, applying RZ-TRADEOFF to identify potential sites and characteristics of interest can help prioritize where to conduct further analysis, and ultimately reduce the need for extensive field measurements and data collection.

RZ-TRADEOFF is also novel in its user accessibility. The model interface is an Excel spreadsheet in which site characteristics and weather data can be directly entered. As opposed to process-based models, RZ-TRADEOFF does not rely on hard to obtain parameters such as streamflow and sediment load data nor does it require extensive calibration. Furthermore, all the characteristics required to run the model are digitally available for free online through organizations such as the USDA-NRCS and the National Weather Service and can be readily accessed if the user does not have this information in their own records (see details on how to access these data in the RZ-TRADEOFF model accessible at https://dataverse.harvard.edu/dataverse/esf). The Excel spreadsheet format also allows the model results to be easily copied and exported. This provides the ability for the model output to be further analyzed with other statistical software to meet the needs of each user.

4.6. What Are the Limitations of RZ-TRADEOFF?

The primary limitation of RZ-TRADEOFF in its current state is the high uncertainty of the individual models. The N₂O, PO₄³⁻ concentration, and CO₂ models have low (<0.3) _adjR² values, with the best goodness-of-fit being an _adjR² of 0.65 for the water table model (Table 5). While these results are an improvement over previous modeling efforts, the low _adjR² values for many of the models indicate that there is still substantial variability not accounted for by RZ-TRADEOFF. This suggests there may be high uncertainty associated with model output obtained when applying RZ-TRADEOFF. Additional uncertainty may exist because of the input dataset used to develop the RZ-TRADEOFF models containing data for many variables from a relatively limited number of sites. Uncertainty that could result from the large number of variables in the development dataset was however addressed during the initial phase of the analysis by conducting PCA and correlation analysis to reduce the dimensionality of the dataset prior to model development [104,105]. Another limitation of RZ-TRADEOFF is that the models were developed from data collected from sites in the Northeast and Midwest and is therefore only applicable to riparian sites located in these regions. Further, the Northeast and Midwest are in themselves large geographic areas, and it is certainly possible that there are riparian sites with conditions and characteristics not captured by the sites included development dataset. This inherent natural variability in riparian buffers may not be accounted for in RZ-TRADEOFF, resulting in further uncertainty.

Due to these limitations, RZ-TRADEOFF may not be appropriately applicable to all riparian sites, or for addressing all management related concerns. Though we maintain that RZ-TRADEOFF can be a highly useful tool for predicting riparian functions, we also stress that users be aware of the limitations associated with RZ-TRADEOFF. The uncertainties associated with the models should be thoroughly considered when applying RZ-TRADEOFF and using the model predictions to inform management decisions.

5. Conclusions

An important step in the future will be to engage and collaborate with colleagues to test and optimize individual models within RZ-TRADEOFF. Improving RZ-TRADEOFF beyond its current state is key in further enhancing our ability to predict riparian zone functions considering the remaining uncertainty in the model. The format of RZ-TRADEOFF and the development database allows it to be easily updated with new data, thus additional riparian site information can be added to the database and individual models reworked in the future as necessary. For instance, if more efficient models are developed in future studies, the spreadsheet equations in RZ-TRADEOFF can be updated to a version that reflects the changes. While there are some limitations and uncertainty in the RZ-TRADEOFF models, the results of this analysis are important in highlighting the influence of hydrogeomorphology relative to weather for most riparian functions.