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

Riparian zones are often used as best management practices due to their ability to remove nitrate (NO3−) 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 (NO3− and phosphate (PO43−), concentration and removal in subsurface flow, total phosphorus (TP) removal in overland flow, nitrous oxide (N2O), methane (CH4), and carbon dioxide (CO2) 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, PO43− and CO2) 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 NO3− concentration, −7.83% and 20.82 for NO3− removal, 6.64% and 0.35 for PO43− concentration, 2.55% and 0.17 for TP removal, 40.33% and 0.23 for N2O, 72.68% and 0.18 for CH4, and −34.98% and 0.91 for CO2. 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.


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 3 − ) loading on water quality [1,2].
Studies have noted the potential for riparian zones to have NO 3 − 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 3 − ) 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 2 O), carbon dioxide (CO 2 ), and methane (CH 4 ) 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 2 O and CH 4 in riparian zones [17][18][19]. Within riparian zones susceptible to flooding, an inverse relationship between soil moisture and N 2 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 2 sinks [22], while others are a source of CO 2 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 3 − removal and that changes in vegetation cover have limited impact on NO 3 − removal [27]. In general, REMM predicted riparian zone water table depths and NO 3 − 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 3 − concentration at the field edge and NO 3 − removal in riparian subsurface flow, but also phosphate (PO 4 3− ) concentration, GHG emissions (N 2 O, CO 2 , CH 4 ), 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.

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). 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). 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). 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 ‡  4 = methane emissions at the soil-atmosphere interface; CO 2 = carbon dioxide emissions at the soil-atmosphere interface; N removal = nitrate removal in subsurface flow; NO 3− = nitrate concentration at the field edge; N 2 O = nitrous oxide emissions at the soil-atmosphere interface, PO 4 3− = phosphate concentration at the field edge; PO 4 3− removal = phosphate removal in subsurface flow; TP = total phosphorus removal in overland flow; water table = water table depth below ground surface. only considered in subsurface flow since research suggests subsurface flow is the primary transport mechanism for both PO 4 3− and NO 3 − 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 2 ) 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 2 ( adj R 2 ) values to reduce over-fitting the model as would occur if simple R 2 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.

Model Validation
The predictive ability of the RZ-TRADEOFF model was validated using additional independent data from the published literature (See Tables 3 and 4 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/cdoweb/) 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. 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     Peat (mg C m −2 d −1 ) n/a n/a n/a 0 †: CH 4 = methane emissions at the soil-atmosphere interface; CO 2 = carbon dioxide emissions at the soil-atmosphere interface; NO 3− removal = percent nitrate removal in subsurface flow; NO 3− = nitrate concentration at the field edge; N 2 O = nitrous oxide emissions at the soil-atmosphere interface, PO 4 3− = phosphate concentration at the field edge; PO 4 3− 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.
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].

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 2 equivalents, groundwater flux, NO 3 − flux, and the mass of NO 3 − 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 −1 ) per meter of stream length. Ks is the saturated soil hydraulic conductivity (m day −1 ) 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

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). Four of the models, CO 2 , PO 4 3− , and NO 3 − 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.  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.
The water table depth model developed for RZ-TRADEOFF is a linear model including two temporal and four hydrogeomorphic variables ( Adj R 2 = 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 ( concentration at the field edge, and PO 4 3− 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 4 3− concentration at the field edge model had an Adj R 2 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 4 3− removal in subsurface flow, which means that the variables used in our database cannot explain PO 4 3− 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 2 O, CH 4 , and CO 2 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 4 emissions in non-peat sites while emissions in peat sites were represented by a constant value ( Table 5). The CO 2 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).

Model Validation
Data compiled from a total of 33 independent published studies (i.e., not used for model development; Tables 3 and 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). The PBIAS calculated from the validation observations ranged from −34.98% for the CO 2 model to 72.68% for the CH 4 ( Table 6). The MAPE ranged from 6.75% for TP removal to 100.10% for N 2 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 3 − removal (PBIAS −7.83%; MAPE 18.73%). The TP removal and CH 4 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 2 O and CH 4 models at 0.23 and 0.18 respectively, while high for CO 2 at 0.91 (Table 6). 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. 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%).

To What Extent Does RZ-TRADEOFF Advance Our Ability to Predict N Related Riparian Functions?
The RZ-TRADEOFF models for NO 3 − concentration at the field edge and percent NO 3 − 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 3 − 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  [32]. The RZ-TRADEOFF RMSE value for NO 3 − concentration at the field edge was consistent with the midpoint of this range at 3.7 mg L −1 , 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 3 − 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 3 − removal and NO 3 − 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).

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 4  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 4 3− 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 4 3− concentration model however, is more accurate than precise, and may be more useful when applied at the watershed scale.

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 2 O, CH 4 , and CO 2 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 2 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 4 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 4 fluxes (4 orders of magnitudes). A constant is used for CH 4 fluxes in peat sites because the values in the development dataset were much higher (>600 mg m −2 d −1 ) than that of the non-peat sites (−0.5-26.67 mg m −2 d −1 ). Therefore, including the peat sites in the Michaelis-Menten model would reduce its sensitivity for the other sites. Some studies report mean peatland CH 4 emission values lower than 50 mg m −2 d −1 [96,97], while others are between 100-200 mg m −2 d −1 [98,99].  [100]. This range and variability of fluxes is not captured by the development dataset, therefore the constant value reported here for CH 4 emissions may not present a wholly accurate representation of CH 4 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 4 peat model, the GHG components of the RZ-TRADEOFF model for N 2 O, CO 2 , and CH 4 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 2 O and CH 4 models, indicating that these two models may provide useful estimates at the watershed scale nonetheless.
Of note is the positive correlation between N 2 O emissions and confining layer depth in our model (see Table 5). Several studies have noted that groundwater fluctuations and flooding regimes influence N 2 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 2 O emissions.

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 3 − 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 3 − 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.

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 2 O, PO 4 3− concentration, and CO 2 models have low (<0.3) adj R 2 values, with the best goodness-of-fit being an adj R 2 of 0.65 for the water table model (Table 5). While these results are an improvement over previous modeling efforts, the low adj R 2 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.

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.