Reach-Scale Model of Aquatic Vegetation Quantifies N Fate in a Bedrock-Controlled Karst Agroecosystem Stream

: In-stream fate of nutrients in karst agroecosystems remains poorly understood. The significance of these streams is recognized given spring/surface water confluences have been identified as hotspots for biogeochemical transformations. In slow-moving streams high in dissolved inorganic nutrients, benthic and floating aquatic macrophytes are recognized to proliferate and drastically impact nutrient fate; however, models that quantify coupled interactions between these pools are limited. We present a reach-scale modeling framework of nitrogen dynamics in bedrock-controlled streams that accounts for coupled interactions between hydrology, hydraulics, and biotic dynamics and is validated using a multi-year, biweekly dataset. A fluvial N budget with uncertainty was developed to quantify transformation dynamics for the dissolved inorganic nitrogen (DIN) pool using a GLUE-like modeling framework, and scenario analyses were run to test for model function over variable environmental conditions. Results from a 10,000 run uncertainty analysis yielded 195 acceptable parameter sets for the calibration period (2000–2002), 47 of which were acceptable for the validation period (2003) (Nash-Sutcliffe Efficiency (NSE) > 0.65; percent bias (PBIAS) < ±15), with significantly different posterior parameter spaces for parameters including denitrification coefficients and duckweed growth factors. The posterior solution space yielded model runs with differing biomass controls on DIN, including both algae and duckweed, but suggested duckweed denitrifies at a rate that would place the bedrock agroecosystem stream on the high-end of rates reported in the literature, contradicting the existing paradigm about bedrock streams. We discuss broader implications for watershed-scale water quality modeling and implementation strategies of management practices for karst agroecosystems, particularly with respect to stream restoration.


Introduction
Agricultural watersheds with karst terrain are highly vulnerable to nutrient leaching, and have been a significant source of dissolved inorganic nitrogen (DIN) loads to regional waterbodies such as the Gulf of Mexico and Chesapeake Bay [1][2][3][4][5]. Small tributaries streams in these landscapes are important given they have high capacity to remove nitrate efficiently because of their high ratios of streambed area to water volume and can account for much of the stream length within a drainage network [6][7][8]. Much work in recent years has focused on the impact of hyporheic connectivity to impact in-stream N removal (e.g., [9,10]). Surface headwater streams in karst agroecosystems are often bedrock controlled and lack prominent hyporheic interaction, which would tend to suggest limited potential for DIN attenuation and transformation [11]. However, surface waters immediately downstream of surface-groundwater interfaces (e.g., springs) in karst landscapes have been identified as hotspots for biogeochemical transformations [12]. In-stream aquatic vegetation such as benthic algae and floating aquatic macrophytes can flourish in these environments and influence transient storage and transformation of nutrients, yet their impact on the stream N budget remains poorly constrained [11,13,14].
Previous research on N removal in low-gradient agricultural stream channels has emphasized algal uptake and denitrification in sediment and detrital organic matter as governing mechanisms of DIN removal. Open-canopied agricultural streams provide favorable conditions for autochthonous algal production, which has been found to drive gross primary production and ecosystem respiration in many streams [15]. Algal uptake reflects a transient sink for N and may ultimately be regenerated to the stream channel or transported downstream as detrital organic matter or scoured biomass [16,17]. Scour and downstream export of algae has been highlighted as a significant component of the fluvial N budget in low-gradient agricultural streams [18]. Regarding permanent removal via denitrification, several studies have found benthic sediments and organic matter to support significantly higher rates of denitrification than plant material [19][20][21]. This may partially stem from decomposition of the particulate organic carbon in detrital organic matter that enhances denitrification through reduction of dissolved oxygen concentrations and creation of an anaerobic habitat [22]. Nevertheless, denitrification does occur in algal mats, and algae enhances denitrification rates because of the supply of a labile source of organic C to fuel heterotrophic bacteria [16,23]. Despite findings of high denitrification rates in agricultural reaches for benthic sediments, and to a lesser extent benthic algae, many studies note that in-stream denitrification did not substantially reduce total nitrate export from watersheds [21,[24][25][26].
Although limited research has been conducted on the influence of duckweed on fluvial N budgets in surface streams, research on duckweed-algae interactions in other landscapes suggest high N removal rates that can exceed algae and detrital N removal. Much of the nutrient removal research regarding duckweed has focused on wastewater ponds, stormwater detention basins, and constructed wetlands [27][28][29][30][31][32][33]. Results from these landscapes have shown that duckweed grows rapidly in N-rich environments and is highly efficient at removing N over long-periods of time, with active life-spans of mats exceeding 25 days [34]. As a result of rapid growth rates, duckweed is often harvested to optimize nutrient removal [27]. Denitrification rates in duckweed have also been found to be high, and are impacted by biomass to water volume ratios, velocity, and nutrient enrichment and have been found to be on the same order of magnitude as biomass uptake [28,31,35,36]. In studies where both algae and duckweed are present, duckweed has been found to have a stronger effect on permanent N removal rates [28,30,31,36]. Given duckweed can grow rapidly, it may deplete the available N pool in the water body which can induce a nutrient limitation stress to algae [37][38][39]. Further, several studies have found that the surface cover of duckweed prevents sunlight from penetrating into the underlying water column, inhibiting algal growth [30,31,36,40]. This can result in temporal variability of biological controls that are impacted by environmental conditions of temperature and light availability [31]. A need exists to quantify the N budget for streams impacted by both duckweed and algal dynamics.
Reach-averaged numerical modeling provides an economically feasible alternative to continuous monitoring for quantifying stream N budgets; however, existing models have limitations that restrict their ability to simulate N dynamics in karst agroecosystem streams. In reviewing existing numerical models with regard to their capability of simulating aquatic vegetation pools and the associated impacts on the stream N cycle, we find three primary limitations. First, parsimonious models of duckweed-algae interactions and their influence on the N cycle have been developed and applied primarily for wastewater ponds and treatment wetlands [30,31,[41][42][43]. Existing reach-scale stream models that incorporate both benthic algae and floating aquatic macrophytes, such as WASP and AQUATOX, are complex and require extensive parameterization [44,45]. Second, existing models simulating duckweed biomass assume similar first-order growth kinetics to algae [44,45]. Laboratory studies have shown that duckweed requires a variable-order equation derived from a second-order function to capture the rapid growth dynamics and physical surface area limitations on duckweed growth [41]. Third, denitrification is often assumed to vary as a function of detrital organic matter or sediment carbon content, or it is lumped as a single parameter to account for benthic and water column denitrification [30,[44][45][46]. However, field-scale and laboratory studies have found denitrification rates are often directly proportional to the biomass pools and more explicit consideration of denitrification as a function of biomass pools is needed [21,22,28,31].
Improved modeling capabilities will lead to an improved understanding of N removal dynamics and can provide insight into how dynamic environmental conditions and spatial variability of processes may impact N dynamics in the future. Springs are recognized to have high spatiotemporal variability in nutrient concentrations, which stem from soil variability and upland management practices [3,11]. Likewise, spring-surface water interfaces are hotspots for biogeochemical processes, and downstream reaches will experience gradients in temperature and water chemistry that could alter N removal dynamics [12]. Environmental drivers have also been found to have a significant impact on stream N removal dynamics. Flow variability and timing of stormflows has been found to significantly influence biomass residence times in the region [16]. Changing precipitation patterns are strongly expected to impact N loading dynamics from uplands in the future [47], and changing temperatures may alter in-stream vegetation dynamics [16]. A need exists to evaluate the behavior of vegetation pools and N removal dynamics as a function of environmental drivers and spatial heterogeneity of processes, which can ultimately result in improved N management strategies.
The overarching objective of this work was to investigate and quantify impacts of aquatic vegetation and associated biota on DIN dynamics in karst agroecosystem streams with bedrock control through development and application of a numerical model. Specific objectives included: (1) develop a parsimonious model that couples algae, duckweed, and organic matter recycling that explicitly links biomass pools to denitrification dynamics in order to quantify stream N dynamics in bedrock controlled agroecosystems, (2) evaluate the model in a karst bedrock agroecosystem stream and quantify relative roles of uptake and denitrification removal processes from different pools, and (3) determine the sensitivity of the N removal processes under varying environmental conditions to improve understanding of impacts of climate, landcover change, and spatial heterogeneity on N removal in these landscapes.

Model Formulation
To meet the objectives of this study, a reach-scale numerical model was developed that simulates the influence of vegetative pools on the fluvial N cycle for low-gradient bedrock stream channels. The sections below detail the model framework used to reflect the aforementioned conceptual model, and the equations and numerical methods used to quantify dynamics.

Framework for N Cycling in Bedrock Streams
In an effort to capture relevant processes, while also considering model parsimony, this model simulates the fluvial N cycle in vegetated bedrock streams using four state-variables that capture dynamics in the biomass pools and their associated impact on the total DIN pool ( Figure 1). The four state variables are (1) DIN: total dissolved inorganic nitrogen concentration, (2) Alive: living algal biomass, (3) DWlive: living duckweed biomass, and (4) OMDet: detrital organic matter biomass. For each state variable, a governing mass-balance is included that considers presumed important processes. Microbial biomass is not explicitly considered as a state variable; however, N transformation dynamics mediated by microbes are explicitly accounted for by considering biochemical fluxes and associated rate limitations, e.g., biomass saturation, light, nitrate availability, and temperature [16,30,44,48].

Model Equations
The numerical model simulates the aforementioned model framework using one-dimensional mass-balance equations that consider advection into and out of the stream reach and biochemical processes that impact the fate of the specified state variables. The model uses a finite differencing numerical scheme to solve the governing mass-balance equations, which is a common approach for similar reach-scale nutrient models [16,50,51]. Spatial discretization is handled through simple routing between reaches based on user-supplied hydrologic time-series at reach boundaries. The following sections define the spatially (j) and temporally (i) discretized mass-balance equations for DIN, Alive, DWlive, and OMDet (gN).

Dissolved Inorganic Nitrogen Mass-Balance
Dissolved inorganic nitrogen, DIN (gN), was simulated as a function of advective inputs and outputs, assimilative uptake (U) and regeneration (R) by biomass pools, and permanent removal via denitrification (Den) as follows.
where, Qin is the inflow at the upstream boundary of the reach (m 3 d −1 ), is the DIN concentration at the upstream boundary (gN m −3 ), Qout is the outflow at the downstream boundary of the reach (m 3 d −1 ), is the DIN concentration at the outlet of the reach and is assumed equal to the DIN concentration in the reach during the previous timestep (gN m −3 ), R is the regeneration of N from a biomass pool by endogenous respiration, decomposition or hydrolysis (gN d −1 ), m is the biomass pool (Alive, DWlive, or OMDet), U is the DIN uptake and assimilation rate and is assumed to equal zero for the detrital pool (gN d −1 ), Den is the denitrification rate associated with each biomass pool (gN d −1 ), and ∆ is the model timestep (d). Flowrates at the upstream and downstream boundaries of all reaches (Qin and Qout, respectively) are supplied to the model. Concentrations of DIN into a reach is supplied at the upper most reach input boundary and calculated for subsequent downstream reaches. Average concentrations of DIN in a reach at a given timestep are estimated as follows.
where, V is the volume of water in the reach and is calculated at each timestep using a water massbalance.
Denitrification was simulated as a function of biomass pool size, temperature and nutrient limitation [30,52]. Denitrification flux (gN d −1 ) associated with microbes living on a given biomass pool (Denm) was simulated as follows.
where, is the denitrification rate constant for each biomass pool, m (gN gNm −1 day −1 ), is the Arrhenius constant, is the water temperature (°C), and is the half saturation constant of available N for microbial denitrification (gN m −3 ).

Algal Mass-Balance Model
Algal N dynamics are simulated using a modified formulation from Rutherford et al. [53]. Briefly, Rutherford et al. [53]'s model was developed to simulate the epilithic algal C biomass at the reach-scale for nutrient-rich agroecosystem streams, similar to the landscape in this study. While the original model assumed negligible nutrient limitation impacts on growth rate, we modify the growth term to consider N limitations and then apply C:N ratios to quantify N biomass dynamics. We also incorporate a slough and death term that considers the explicit mass flux to the detrital organic matter pool. The mass of algae, Alive (gN) is estimated at each spatial and temporal step using the following mass-balance.
where, is the algal colonization rate (gN d −1 ), and DA is the death/sloughing rate of algae that is added to the detrital pool (gN d −1 ).
Scour of algae occurred during high flow events, and it was assumed to be completely lost from the system given that algae is relatively neutrally buoyant and would not be expected to settle out of flows that are high enough to cause scour [16]. This is accounted for using a piecewise function based on a critical discharge threshold, QcA (m 3 d −1 ). When flowrate in a reach exceeded the critical threshold, algal biomass was reset to a near-zero seed value to ensure recolonization [17].
The assimilative uptake rate of algal N ( ) is estimated using a modified version of the Rutherford et al. [53] growth model considering light, temperature, population saturation, and nutrient limitations as follows: where, is the maximum uptake rate constant (gC m −2 d −1 ), is the nitrogen to carbon ratio of algae (gN/gC), is the light intensity limitation coefficient (dimensionless), is the temperature limitation coefficient (dimensionless), is the population saturation limitation coefficient (dimensionless), is the half saturation constant of available DIN (gN m −3 ), and is the surface area of the streambed for the specified reach.
Light limitations are estimated based on photosynthetically available radiation [53] and consider the impact of floating aquatic macrophytes to attenuate light.
where, IA is the photosynthetically available radiation incident on the surface of the algal mat (µmol m −2 s −1 ) and is the saturating radiation constant for algae (µmol m −2 s −1 ). The available radiation at the algal mat is estimated as the available radiation at the water surface ( ) minus a linear attenuation term that accounts for the relative amount of duckweed biomass.
where, I is the photosynthetically available radiation incident on the water surface (µmol m −2 s −1 ), DWlive is N duckweed biomass (gN), and DL is the duckweed mat density limit (gN m −2 ). Temperature limitations are assumed to follow an asymmetrical Gaussian distribution when temperature of the water deviates from optimum temperature for epilithic algal growth [53].
where, is the water temperature (°C), is the optimum water temperature for epilithic algae (°C), is the minimum temperature for epithilic algae (°C), is the maximum temperature for epithilic algae (°C). Assuming As population increases, light availability to the entire algal mat decreases as deeper cells become shaded [53]. To account for this, population consequences are accounted for as follows.
where, is the saturation density-dependence coefficient (gC m −2 ). The algal colonization rate refers to the rate of colonization by algal cells from upstream reaches and is calculated as follows.
where, is the algal colonization rate (gC m −2 d −1 ). The death and sloughing term for algae (DA) refers to the transfer of live algal biomass to the detrital organic matter pool and is simulated using first-order kinetics [49]. The term collectively accounts for losses due to grazers, sloughing, and death of algal biomass.
where is the periphyton death and sloughing rate (day −1 ) from Chapra et al. [49]. Regeneration of algal biomass to the DIN pool occurs through endogenous respiration of algal biomass and is simulated using an analogous approach to Rutherford et al. [53] as follows.
where is respiration rate constant (day −1 ) and is the reference temperature (°C) at which is measured, and is the temperature coefficient for algal respiration (analogous to the Arrhenius constant).

Duckweed Mass-Balance Model
Simulation of duckweed N mass (DWlive) was based on processes observed, and models developed in, wetland and wastewater pond environments [30,41,42]. The model considers uptake, endogenous respiration, advective transport, and mortality of duckweed in the mass balance as follows.
where, is the death/mortality rate of duckweed (gN d −1 ). Scouring due to high flows was accounted for using a piecewise function based on a critical discharge threshold, QcDW (m 3 d −1 ). When flowrate in a reach exceeded the critical threshold, duckweed biomass was reset to a near-zero seed value to ensure recolonization ( ) [17,44] in the same way as algal biomass. The DIN assimilation rate of duckweed ( ) follows variable-order kinetics [41] and is estimated as follows.
where, is the duckweed mat density limit (gN m −2 ), and is intrinsic growth rate of duckweed (day −1 ) which varies as a function of temperature, light intensity, and available DIN as follows.
where, is the maximum growth rate (day −1 ), is the temperature coefficient (Arrhenius constant) for duckweed growth and death, is the reference temperature for duckweed (℃), is the saturating radiation constant for duckweed (µmol m −2 s −1 ), and = half saturation constant of available DIN for duckweed (gN m −3 ).
Death rate of duckweed follows first order kinetics [30] and reflects the detrital OM component for duckweed.
where, is the mortality rate of duckweed (day −1 ) and follows a piecewise function in order to account for severe environmental conditions in which = 0.05, Regeneration of DIN to the stream from duckweed occurs through endogenous respiration and is simulated using similar first-order kinetics to algae [53] due to the influence of temperature on respiration [54].
where, is duckweed respiration rate (day −1 ), and PKrespDW is the temperature coefficient for duckweed respiration (analogous to the Arrhenius constant).

Detrital OM Mass-Balance Model
Detrital organic matter (OMDet) receives inputs from the live organic matter pool and is balanced by regeneration of inorganic N to the DIN pool, and scour and subsequent downstream transport.
where, ROM is the regeneration of detrital organic matter to the DIN pool (gN d −1 ). Regeneration of detrital organic matter to the stream channel occurs through microbial decomposition and hydrolysis and hence is modeled as first-order reaction as follows.
where = hydrolysis and decomposition rate (day −1 ) from Chapra et al. [49]. Scouring of detrital organic matter due to high flows was accounted for using a piecewise function based on a critical discharge threshold, QcOM (m 3 d −1 ). When flowrate in a reach exceeded the critical threshold, the detrital organic matter was reset to a near-zero value to preserve model continuity [17] in the same way as algal and duckweed biomass.

Study Site and Materials
Camden Creek is a spring fed, bedrock-controlled stream located within a karst agroecosystem watershed in the Inner-Bluegrass Region of Central Kentucky. The Camden Creek watershed (drainage area of 10.69 km 2 ; Figure 2) is characterized by broad, shallow sinkholes, low relief, broad valleys and ridges, sparse rock outcrops, and thick, fertile, limestone and shale residual soils over phosphatic Ordovician limestone [55]. Camden Creek and surface tributaries in the watershed are shallow, emanate from springs, flow over limestone bedrock, and are generally unshaded through grazed pasture with some riparian vegetation [56], with low streambed sediment storage on exposed bedrock [11]. As a result of spatial variability in land use, nutrient concentrations varied across spring inputs to the stream channel. In general, nutrient concentrations at springs were high. As an example, measured nitrate-N concentrations at spring site Sp1 ranged from 3.8-13.6 mgN L −1 , with an average of 6.36 mgN L −1 . Stream sites (ST8, ST7, and ST4) had nitrate levels below 1 mg L −1 and sometimes below detection limits during the summer months in multiple years, contrasting the high nitrate concentrations found at the spring sites. This highlights the importance of in-stream N removal in the stream and the potential for N-limiting conditions during summer. Regarding dissolved reactive phosphorus (DRP), Ford et al. [11] found average slow flow concentrations of 0.233 mg L −1 , nearly an order of magnitude higher than eutrophic thresholds of 0.02-0.03 mgP L −1 for freshwater algal proliferation [39,57], suggesting DRP is likely not a rate-limiting nutrient in this system. A roughly 1 km stretch of stream was selected for the model application and was discretized into two stream reaches based on tributary inputs. Reach 1 (~450 m) and 2 (~570 m) refer to the reach of Camden Creek between the junction of ST8 and Sp1 through ST4 ( Figure 2). This model domain was selected to evaluate the aforementioned model because (1) nutrient concentrations are high and benthic sediment storage is low, (2) duckweed, algae, and detrital biomass are all well recognized to proliferate in this section of the channel based on long-term monitoring conducted by the authors, (3) availability of long-term flow and nutrient datasets for model evaluation, and (4) nitrate levels decrease longitudinally downstream, reflecting a significant impact by stream vegetation [11]. Flow data was collected at ST4, ST7 and ST8 across different periods of the long-term monitoring effort. Flowrate data for ST4 was available periodically from November 1994 through December 2004 and continuously from 2000-2003, ST7 from May to December 1999, and ST8 periodically from September 1997 to June 1998 and April to December 1999. Flow depth data was collected using pressure measurements from an ISCO 4220 flow meter with pressure transducer at the weirs for each sampling site. The weirs at all three locations were 90° v-notch, 10 feet wide, with flowrate (m 3 s −1 ) estimated using a piecewise function to account for flow within and overtopping the v-notch weir.
Water quality data, including nitrate (NO3 − ) and total ammoniacal-N (TAN), were collected at the specified stream (ST) and spring (Sp) sites throughout the watershed ( Figure 2). Of importance to this modeling effort, data was collected at ST8, Sp1, ST7, ST6, and ST4. Grab sampling began in October 1996 and was conducted through June 2007, with unpreserved samples for nitrate collected in 250 mL amber glass bottles and preserved samples for ammoniacal-N collected in 250 mL clear glass bottles, with all samples placed on ice immediately after sampling and delivered to Kentucky Geological Survey (KGS) laboratories within 6 h of collection [11]. Nitrate was analyzed on a Dionex Ion Chromatograph within 48 h of sample collection and ammonia-N was analyzed colorimetrically using a UV Vis spectrometer by Varian within 28 days of sample collection [11].

Model Inputs and Parameterization
Hydrologic and hydraulic inputs into the model included flowrates at upstream and downstream boundaries of each reach and channel geometry. Regarding geometry, the channel cross-section was assumed rectangular, which was determined to be appropriate based on our site visits. Average channel width and length were estimated using aerial imagery from the site for each reach. Regarding flowrates, there were three primary hydrologic inputs into the modeled stream reaches, ST8 and SP1 at the upstream of Reach 1 and ST6 at the upstream boundary of Reach 2. A four-year span of continuous flowrate data collected at ST4 (January 2000-December 2003) was utilized as the downstream flowrate boundary condition at the outlet of Reach 2. In absence of continuous flow data at the three hydrologic inputs, comparisons of flowrate data from overlapping data collection efforts at ST8, ST7, and ST4 from April 1999-December 1999 were used to establish relationships between flowrates at upstream locations using ST4 as a response variable. To partition flow inputs to Reach 2 from ST7 and ST6, the relationship between the flow ratio (RST7 = QST7/QST4) and flowrate at ST4 was used. The flow ratios from 0.037 m 3 s −1 to 0.34 m 3 s −1 showed no discernable trends, and an average of 76% was used for ST7, with the remaining 24% coming from ST6 for those flow conditions. Partitioning between ST8 and SP1, flow ratios (RST8 = QST8/QST7) showed a linear increase (R 2 = 0.72) with flowrate for QST7 > 0.017 m 3 s −1 . This linear regression model was used to calculate flow ratios for flows above the threshold. For low flow conditions, two parameters ( R and R ) were included given the high uncertainty in flow ratios. The model timestep was set to a 30-min interval based on average flow travel times in the two reaches.
Environmental inputs to the model include water temperature and photosynthetically available radiation (PAR) which were obtained using a mixture of atmospheric data from a nearby gauging station and correlation with field measurements. The water temperature input was estimated continuously over the four-year period based on long-term air temperature measurements at the research location, as well as a linear regression model of air temperature vs. water temperature using high-resolution measurements from 29 August 2018 through 1 January 2019 at the watershed outlet (ST1) and air temperature data from a NOAA weather station co-located at the study site. Simple linear regression suggests the NOAA air temperature explains much of the variability in the water temperature dataset (Water Temp = 0.4174 (Air Temp) + 9.9629; R 2 = 0.83). Water temperature inputs within the modeled stream reach were corroborated by comparing the model inputs with biweeklymonthly water temperature measurements at ST8, Sp1, and ST4 from 2000-2003. The PAR was estimated from solar radiation data collected at the Blue Grass Airport (roughly 10 miles to the west of the research site) and managed by the National Solar Radiation Database (NSRDB). The total solar radiation from NSRDB was converted from W m −2 to PAR in µmol m −2 s −1 using a conversion factor of 2.02 [58].
Boundary conditions for DIN concentrations were considered using flow weighted averages of hydrologic inputs at upstream nodes in the model and considering seasonality, flow, and annual variability of concentration dynamics at the flow sources. The boundary DIN concentrations were compared with flowrate at ST4. The DIN concentration for each of the input sites (ST8 and Sp1) was plotted against flowrate at ST4 at that time for each of the seasons outlined above over the four-year span 2000-2003. These seasonal concentration-flowrate relationships were used to calculate the DIN concentration at each time based on the flowrate for each season in each year to better constrain the boundary conditions for year-to-year variability. To avoid over-prediction of concentrations at high flows when only low flow conditions are present in the measured dataset, the equations were capped at the highest DIN concentration for that particular season and year. While concentrations are recognized to dilute at peak flow conditions this would provide a conservative (low-end) estimate of N removal percentages for the system. The inputs from ST6 were left as seasonal averages due to the lack of flowrate data for 2004-2005 when nutrient concentrations were collected.
The parameters used for model calibration were obtained from a mixture of published studies on algal and floating vegetation models, published data on growth/denitrification rates from stream and pond systems, and physical system constraints (shown in Table 1). Since parameter distributions were unknown, all prior parameter ranges were assumed to have uniform distributions, which is typical of other uncertainty analyses of stream water quality models [51].
The algal and detrital biomass input parameters were obtained from previous modeling and monitoring studies in nutrient rich stream channels [21,22,49,51,53]. Ranges for algae growth and respiration parameters (Pmax, IkA, Tmin, Topt, Tmax, Psat, PrespA, TrefA, PKrespA, Pcol) were obtained from Rutherford et al. [53] and have been successfully applied in a nearby agricultural stream [50]. Parameters for death and decomposition/hydrolysis of algal and detrital organic matter (kd, kh) were derived from Chapra et al. [49]. The half-saturation constant of available DIN for algae (khsA) was assumed to vary from a minimum value of zero to a maximum value reflecting a conservative half saturation constant for nitrate [30], given ammonium was rarely present in the system. This range encompasses other algal half saturation constants for N used by parsimonious modeling studies that simplify DIN to a single pool [49]. The ranges for algal and detrital organic matter denitrification parameters (kDenA, kDenOM) were intended to be as conservative as possible, using zero as the minimum and maximum denitrification rates for benthic plant material and benthic organic matter reported in the literature [21,22]. The half saturation constant of nitrate for denitrification (in all pools), was obtained from studies in the wastewater, activated sludge, and bioreactor communities [59][60][61] and compares favorably with half-saturation values reviewed in Arango et al. [22]. The maximum critical stream flowrate (QcA, QcOM) values used for advective transport of benthic algae and dead organic matter were based on Kazama and Watanabe [17], which suggest complete removal of all biomass or matter when the flow reaches or exceeds the top five percent of annual flows. The minimum bound was set to zero to be as conservative as possible. The duckweed input parameters (DL, rmaxDW, θDW, km, TrefDW) were obtained from duckweed nutrient removal studies conducted in laboratory, wetlands, and wastewater ponds [28,30,41,42]. The respiration rate, temperature coefficient, and minimum biomass (PrespDW, PKrespDW, DWmin) were assumed analogous to algal biomass, with the respiration temperature coefficient range consistent with the range for the Arrhenius constant. The half-saturation constant of available DIN for duckweed uptake (khsDW) was based on results Lasfar et al. [42]. The duckweed denitrification parameter (kDenDW) range was intended to be as conservative as possible, using zero as the minimum and the maximum obtained from the maximum reported denitrification rate for agricultural streams [62] normalized by maximum N content of duckweed (gN g dry −1 ) from Körner and Vermaat [28]. The formulation for the maximum critical stream flowrate value for duckweed (QcDW) was also based on Kazama and Watanabe [17], as duckweed floats on the surface and would be washed downstream when the flow exceeds a certain threshold. The maximum bound for the removal of the duckweed biomass was again assumed to be the top five percent of annual flows, analogous to algae and organic matter. The minimum bound was set to 6000 cubic meters per day (m 3 d −1 ), roughly 2.5 cubic feet per second (ft 3 s −1 ), which is the flowrate through the 90-degree weirs at ST8, ST7, and ST4 when the water level is at the top of the weir. When the water level is below the top of the weir, the flow is backed up and the surface generally remains calm enough for duckweed to accrue in the reach.

Model Evaluation Procedures
For the model evaluation procedures, a GLUE-like uncertainty analysis was performed using a Monte Carlo simulation with randomized parameter inputs for each run to compare measured and modeled results and generate a posterior solution space of acceptable parameter ranges [16,51,63]. Based on iterative model improvements and Monte Carlo simulations with each iteration, 10,000 randomized runs presented the best posterior solution space relative to simulation run time and robustness of the prior solution space. We performed two primary tasks to evaluate the model performance for our case study. First, we performed a sensitivity analysis to evaluate how different components of the model reflected downstream data based on model inputs. Next, we performed a robust model calibration and uncertainty analysis of the posterior parameter and solution spaces using well-accepted model evaluation statistics and performance criteria.
Sensitivity analysis was performed for the case-study application to identify potential impacts of predominant N transformation pools including algae, duckweed, and its associated detrital biomass as well as microbial denitrification associated with each pool. Sensitivity was performed by running a series of scenario analysis on our un-calibrated model to identify potentially important components of the numerical model structure. Scenarios included: an inert conduit with no biotic activity or benthic detritus (1), algal biomass and detritus without denitrification (2) and with denitrification (3), duckweed biomass and detritus without denitrification (4) and with denitrification (5), and both algal and duckweed biomass and detritus without denitrification (6) and with denitrification (7). The scenarios listed above will highlight the differences in the vegetation pools' impact on the stream DIN concentration at the outlet of the reach (ST4) and how that compares with measured data. The exclusion and inclusion of denitrification for each of the specific vegetation pools is intended to highlight the extent of denitrification impact from each pool.
Model calibration and validation was performed for the four-year case-study using wellaccepted water quality modeling statistics. The model response variable used for calibration was CDIN (mg/L) at the reach outlet (site ST4). We compared the biweekly DIN concentrations from 2000 through 2002 (3 years) for calibration, and biweekly DIN concentrations from 2003 (1 year) for validation. The performance criteria used for evaluation of the randomized runs and creation of the posterior solution space were Nash-Sutcliffe Efficiency (NSE) and percent bias (PBIAS). Considering their inclusion in the hydrologic and water quality model calibration guidelines outlined in Moriasi et al. [64], and their recommended uses, these two performance measures require the model to have strong goodness of fit without bias (over or under-estimation). Given the robustness of our boundary condition parameterization, we required statistical criterion to meet 'very good' thresholds for model evaluation statistics [64]. For NSE this required a numerical value of 0.65 or greater for both calibration and validation periods. For PBIAS this required a value of ±15%. For parameterizations that resulted in acceptable model statistics for both calibration and validation, the parameter sets and solutions were accepted into a 'posterior space'.
The posterior parameter space was compared with the prior parameter space. To identify if the model parameters were sensitive in calibration, we used statistical tests to evaluate if statistical differences existed between the prior and posterior parameter spaces. Due to the non-parametric nature of the parameter spaces, the Mann-Whitney U, or Wilcoxon rank sum test in SigmaPlot was used assuming a 5% significance level. We test for equal variance in SigmaPlot [65] using the Levene Median test along with the rank sum test to check for the difference in variability between the prior and posterior solution spaces. These were also verified in Matlab using the Brown-Forsythe test, assuming 5% significance, which is an adaptation of Levene's test that uses the median of the values as opposed to the mean and can provide better performance on heavily skewed distributions [66]. The posterior solution space was used to quantify uncertainty in fluxes and dynamics for in-stream N dynamics.

Model Evaluation Analysis
Results of the sensitivity analysis highlight that inclusion of both vegetation pools and their associated denitrification terms are important in order to capture the temporal variability observed in the four-year dataset (Figure 3). In scenario 1, the minimum-maximum range for DIN concentrations (CDIN) constrains the datapoints for winter months reasonably well with datapoints typically falling between the minimum and median lines. Deviation between measured and modeled values occur in the spring through the fall and are particularly high in the summer with minimum values over-predicting the measured data by as much as 1.66 mgN L −1 (11/8/2000) in the fall. These findings highlight the assumption that an inert conduit is inappropriate for the specified stream reach. Sensitivity analysis of the algae, its detrital organic matter, and associated microbial denitrification show improved capabilities to predict concentrations, particularly in spring and summer, but some limitations in fall (Scenarios 2-3). Minimum CDIN in Scenario 2 were able to bound much of the measured concentrations of DIN during the spring and summer season except for some periods in the fall (see fall of 2000), although median values still vastly over-predicted the data in the growing season. Comparing results of Scenario 3 and Scenario 2 illustrate the effects of denitrification associated with algae and its detritus are subtle and mainly only impact the magnitude of diel fluctuations in CDIN (see the maximum lines during the growing season). Sensitivity analysis for duckweed, its detrital organic matter, and associated denitrification also show the ability to capture the low CDIN observed during the growing season with some improved predictive capabilities over algae, especially in the fall (Scenarios 4-5). Minimum CDIN values in Scenario 4 was able to capture much of the variability in the dataset, and when including denitrification (Scenario 5) median concentrations were also able to capture much of the variability, suggesting that CDIN is highly sensitive to denitrification in the duckweed biomass pool. Comparing Scenario 5 and Scenario 3, we found that duckweed showed improved capacity to capture fall DIN concentrations compared with algae scenarios. Results for both duckweed and algae (Scenario 6) and their associated denitrification (Scenarios 7) more closely reflect the duckweed scenarios, although diel fluctuations were altered because algae was able to offset some of the N regenerated by duckweed (see difference in maximum lines from Scenario 4 and 6). Cumulatively these results suggest that both algae and duckweed may describe dynamics reasonably well, but the CDIN response variable is more sensitive to denitrification associated with duckweed than it is with algae and detrital organic matter.   Figure 4, suggesting highest rates of biotic activity most accurately predict CDIN in this system. However, the acceptable model runs showed the potential for strong diel oscillations during these periods, particularly in dry years (e.g., 2000), suggesting the relationship between DIN removal through denitrification and DIN regeneration through respiration and decomposition may be of particular importance in dry periods of the summer and fall. Of the 23 variable parameters in the uncertainty analysis, 9 were found to have statistically significant differences between the prior and posterior solution spaces by either the Mann-Whitney Rank Sum Test, Brown-Forsythe Test, or both. The posterior parameter spaces found to be significant by the Mann-Whitney Rank Sum Test were kh (p = 0.016), PrespDW (p < 0.001), rmaxDW (p < 0.001), kDenDW (p < 0.001), C-to-N Ratio (p = 0.015), QcDW (p < 0.001), RST7low (p < 0.001), and RST8low (p < 0.001). The posterior parameter spaces found to be significant by the Brown-Forsythe Test were PrespDW (p < 0.001), DL (p = 0.027), rmaxDW (p < 0.001), kDenDW (p < 0.001), QcDW (p < 0.001), RST7low (p < 0.001), and RST8low (p < 0.001).
Histograms for each of these statistically significant parameters are shown in Figure 5. , highlighting the potential for denitrification in the duckweed mats to exert a strong influence on N dynamics in the system. The carbon to nitrogen ratio trends towards the lower end, with a median of 8.99 gC gN −1 , higher than the value of 5.56 gC gN −1 used in Chapra et al. [49]. The decomposition/hydrolysis rate of organic matter, kh, also trends towards the lower end with a median of 0.04 d −1 , indicative of the sensitivity of regenerated DIN on the overall stream DIN concentration. The posterior spaces of both RST7low and RST8low only gave acceptable model solutions for the upper 80th and 50th percentiles of the prior parameter spaces, respectively, emphasizing the impact of main-stem DIN concentrations inputs on overall N dynamics in the study reach. The critical discharge threshold for duckweed, QcDW, skewed towards the upper half of the parameter range with a median of 28,370 m 3 d −1 , indicating the potential for duckweed mats to survive higher flow conditions.   The modeled DIN concentrations, biomass pools and associated denitrification rates are presented from a representative model run (45) in the posterior solution space (Figure 6-right column). Regarding CDIN we found that the calibrated model solution most reflects the minimum values for Scenario 5 (duckweed dominated system with denitrification) reflecting the aforementioned findings of insensitive algal model parameters in calibration. This result is further supported by the modeled biomass values, which highlight differences in the prevalence of the different vegetation pools. Although duckweed biomass dominates total biomass and typically peaks in the summer, algae may be present throughout the winter and begins growing earlier in the spring and continues to grow later in the fall, providing a longer growing season for N uptake and periphyton based denitrification, particularly in the spring of 2000 and 2001 and the fall of 2003. The organic matter pool can be seen to eclipse algal biomass as duckweed increases, as it is made up of detrital material from both pools, and it continues to increase until removed by scour. The modeled denitrification rates for each vegetative pool further highlights the effect of the dominant time periods for each of these vegetative pools, as the associated denitrification is shown to occur concurrently. Duckweed associated denitrification dominates overall denitrification by two orders of magnitude over algae and organic matter.

Nitrogen Budget Results
A comprehensive budget of modeled N transformations is provided (Table 2) input during this summer season was considerably lower than the other summer periods, making the total median removal rate 0.13 kgN d −1 , which was considerably smaller than in previous summers. This lower input DIN decreased overall removal loads through nutrient stress on biota. The lowest median percent removal in summer occurred during 2003 (8.08 percent in the summer). As evidenced in Figure 6, this summer had frequent high-flow events which washed out biomass, limiting residence time for denitrification to occur. These findings suggest moderate flow conditions promote optimal conditions for DIN removal during warm months to minimize environmental stressors on the biomass pools. Comprehensive budget results for biomass pools displayed the importance of uptake for both algal and duckweed pools, but showed denitrification was primarily associated with the duckweed pool (see Supplementary Materials Tables S1 and S2). The DIN assimilation through vegetation uptake was more variable than denitrification, showing dominance by algae or duckweed at different points in time and under different seasonal conditions. Duckweed growth in the winter periods of each year is negligible while algae dominates, contributing all of the DIN removal during those periods (median rates of 0.34 to 0.60 mgN m −2 h −1 ), reflecting the difference in optimum temperatures for the two pools. The spring periods show a similar trend, with median duckweed rates beginning to make an impact (6.04 × 10 −2 to 3.27 mgN m −2 h −1 ), but algal biomass still dominating the overall assimilation (median range of 1.81 to 3.13 mgN m −2 h −1 ). Duckweed was the prominent uptake pool in summer with median rates of 1.16 to 5.49 mgN m −2 h −1 compared to 0.57 to 1.88 mgN m −2 h −1 for algal assimilation, likely indicating the potential for duckweed populations to grow rapidly and shade the benthos during the warmer summer months. In the fall, median algae and duckweed uptake rates were comparable (0.24 to 1.24 mgN m −2 h −1 and 5.55 × 10 −3 to 1.63 mgN m −2 h −1 , respectively). When considering the influence of each potential denitrification pool on the total denitrification removal, duckweed dominates the total rate (median range of 3.58 × 10 −6 to 8.86 mgN m −2 h −1 ), with rates for algae (1.38 × 10 −9 to 0.22 mgN m −2 h −1 ) and organic matter detritus (1.76 × 10 −10 to 0.80 mgN m −2 h −1 ) often several orders of magnitude smaller than duckweed.

Implications for Stream N Modeling
The parsimonious modeling framework from this study was able to provide acceptable modeling statistics for simulation of DIN dynamics in the low-gradient bedrock agroecosystem stream and highlighted the importance of modeling components that are unique to this study. Results of the model calibration and validation highlight in-stream vegetation and microbial N removal processes were able to capture dynamics in spring-fall that were unable to be explained by boundary conditions alone (Figures 3 and 4). Duckweed had the most sensitive model parameters in calibration, with 5 of 12 parameters having a statistically significant difference between their prior and posterior parameter space. This included both the limit density and intrinsic growth rate, as well as the associated denitrification rate of duckweed, which was based on the variable-order growth model of duckweed and the biomass-specific denitrification terms that have not been explicitly considered in existing stream N models [16,44,45,48]. The results of the uncertainty analysis confirm the importance of duckweed and its associated denitrification for stream N removal, as the median duckweed biomass and associated denitrification was often orders of magnitude higher than algae and detrital organic matter. These findings underscore the importance of duckweed biomass and its substratespecific denitrification rates for simulating N removal in other stream ecosystems where floating aquatic macrophytes occur.
The importance of accurate representation of growth kinetics and environmental stressors of biomass pools is further highlighted by the results our sensitivity and uncertainty analysis whereby duckweed was observed to outcompete algae for nutrients and light under favorable environmental conditions, resulting in higher rates of duckweed production and, subsequently, denitrification. The results of the sensitivity analysis in this study showed that both algae and duckweed could bound observed gradients in DIN measurements. Nevertheless, when considering both pools during model calibration, duckweed often outcompeted algae for nutrients and sunlight as evidenced by 37 of the 47 acceptable models runs having higher duckweed biomass than algae. The use of a variable-order growth rate for duckweed biomass enabled a faster proliferation of duckweed that outcompetes algae for DIN, which then exerts a control on light availability. Cumulatively, our results suggest that rapid uptake of duckweed may decrease algal biomass by creating N and light limiting conditions, particularly during summer and fall. Our visual inspection of the stream channel throughout the year qualitatively support these findings, as we often found thick duckweed mats at the water surface during summer that was underlain by senescing or detrital algal biomass. Further, numerous studies in the wetland and wastewater communities have highlighted that duckweed grows rapidly in comparison to algae and that in wastewater ponds and wetlands duckweed has the capability to reduce the abundance of algae due to shading and nutrient limitations during certain periods of the year [29][30][31][36][37][38][39][40].

Stream N Dynamics in Low-Gradient Bedrock Agroecosystem Streams
Broader comparison of model results with other agricultural streams suggest bedrock karst agroecosystem headwater streams can be hotspots for permanent N removal, which contrast existing perceptions. Cumulatively, average denitrification rates for the bedrock stream was in the upper 50th percentile of rates reported for agricultural streams in Mulholland et al. [62], indicating these systems are hotspots for permanent N removal. These findings are counterintuitive given the results of Argerich et al. [67], which showed metabolic activity in bedrock sections of a low-order stream in Oregon were significantly lower than an adjacent alluvial reach. The differences likely are reflective of the steep gradients, canopy cover, and low-disturbance conditions in Argerich et al. [67], which create unfavorable conditions for vegetation proliferation. Our results emphasize the importance of considering duckweed in fluvial N budgets of low-gradient disturbed headwater streams.
Permanent N removal via denitrification in the bedrock agroecosystem stream contrast dynamics observed in wastewater ponds as well as sediment dominated streambeds in the region. Our results showed that denitrification accounted for an average of 46 percent of total N removal in the studied stream reach which was higher than rates reported is wastewater ponds ranging from 10-40% of total N removal [28,31,36]. This finding likely reflects harvesting operations that are commonly performed in wastewater ponds that promote higher removal through biomass uptake as well as the higher velocities of the stream system which allow nutrients to advect to anoxic microsites where denitrification can occur. The budget results highlight the dominance of duckweed on overall denitrification rates in the spring, summer, and fall, often close to 100 percent of total denitrification, but found minimal contributions in winter. This was likely influenced by the favorable environmental conditions during the warm periods, particularly in the summer, creating longer residence times and increased organic carbon retention for more efficient denitrification [26,[68][69][70]. The lack of importance of denitrification in winter contrasts results from a nearby third-order stream with extensive fine sediment deposits and presence of surficial fine-grained laminae [16] which found more than 70% of N removal to be associated with denitrification in winter. Given that the study of Ford et al. [16] was a higher-order stream, our results suggest that there may be variable in-stream control points along the fluvial continuum on N removal dynamics throughout the year in karst agroecosystems.
Results of the study provide insight regarding the coupled impacts of environmental drivers and nutrient gradients on N removal potential of bedrock agroecosystem streams. The optimal DIN removal percentage (52.6%) and removal rate (1.21 kgN d −1 ) occurred during summer of 2001 when temperature and solar radiation were maxima, and flow conditions were moderate, i.e., higher average flow and DIN loading than 2000 or 2002 but less than the wet summer of 2003. Under low flow periods of 2000 and 2001, assimilative and dissimilative DIN removal processes occur at faster rate than DIN inputs, creating rate-limiting nutrient conditions. Conversely, the high flows in 2003 resulted in continuous flushing of duckweed biomass, limiting the amount of denitrification that can occur. As a result of the temporal variability in environmental drivers, we found N removal can vary by an order of magnitude on a year-to-year basis. These findings have implications for landscape variability and behavior of N removal dynamics in the future given land-use change may alter flowrates and nitrate loadings delivered to stream channels, and climatic variables such as temperature and precipitation may also alter the hydrologic and nitrate loading dynamics [47,[71][72][73][74].
We further explored how changes in nitrate loading and temperature may impact the N removal using a scenario analysis of the calibrated model (Figures 7 and 8). Regarding DIN input concentrations, increasing constant DIN input concentrations results in rapid initial increases in biomass and denitrification for all three biomass pools, but each begins to plateau as the concentrations rise, starting with algae (due to nutrient and light limitation by duckweed). Both organic matter and duckweed begin to level off at higher concentrations (due to the impact of duckweed mortality on the detrital pool), indicating that although nutrients are available for uptake, the biomass begins to reach a saturation point for the available stream area (limit density, DL). The denitrification rate in duckweed shows a more sustained increase at the higher DIN inputs due to the available concentrations remaining above the half-saturation constant for denitrification, however, the denitrification rate is still tied to the biomass pool and would be limited by the saturated biomass were concentrations to continue to increase. Regarding temperature, adjusting the average temperature over a ±5 °C range had contrasting impacts on algal and duckweed biomass due to the difference in optimal temperatures for growth. While algal biomass is low under colder and warmer temperatures, duckweed biomass increases under warmer average temperatures, which in turn increases duckweed denitrification. The relatively low biomass values for each pool compared to those under higher DIN loadings in Figure 7, however, suggest the influence of temperature is less important to overall biomass than DIN availability, meaning that the influence of the two conditions should be taken together when considering future scenarios (i.e., higher DIN loads and temperatures together could inhibit algae and substantially increase duckweed growth up to the limiting biomass density). These findings indicate that under potential changes in land-use and climatic regimes, small headwater streams with higher nutrient loading, temperatures, and flowrates could see significant increases in vegetation growth capable of partially offsetting increased inputs. However, population saturation conditions of biomass and flow threshold exceedance for scour and washout may result in reductions in nutrient removal potential under high DIN loading regardless of temperature, which may enhance downstream eutrophication of receiving waterbodies. These findings thus have implications for management strategies which are further discussed in the next section.

Broader Implications
The findings of this study provide insight into the potential spatial and temporal variability of N removal processes in bedrock-controlled streambeds of karst agroecosystems. Briggs and Hare [12] note spring-surface water interfaces as potential hotspots for biogeochemical reactions due to buffered temperature, high nutrients (particularly DIN in agriculturally impacted systems), and low dissolved oxygen. The parameterization of the modeling structure presented here, the results of duckweed literature suggesting growth is optimal at higher temperatures than algae [36,42], and the scenario analysis results that higher temperatures will promote high production of duckweed (and subsequent denitrification) indicate that spring-surface water interfaces may not be optimal locations for duckweed proliferation. Although the incoming nitrate loading may be considerably higher than downstream locations, the lower temperatures after mixing could hinder duckweed growth. To increase denitrification in these locations, enhanced residence times at the interface may be necessary to enable surface radiation to warm the water to more favorable temperatures for duckweed growth. Such improvements may be achieved through implementation of treatment wetlands or stream restoration at these interfaces.
In this regard, stream restoration in disturbed landscapes promote favorable conditions for duckweed biomass. Lorenz et al. [75] found a prevalence of Lemna minor, or common duckweed, in restored reaches in Germany following restoration even in systems that did not have detectable levels prior to restoration. Griffiths et al. [70] notes that while implemented restoration strategies, namely increased floodplain connectivity, in Midwestern agricultural streams promote bank stability and decrease erosion [76][77][78], they may also result in increased residence times, allowing for denitrification of excess nitrate. Increased residence time may also result in increased regeneration of ammonium through decomposition of detrital organic matter, which could induce more coupled nitrification/denitrification in the duckweed mats, resulting in higher permanent removal of nitrate [31]. With the potential of duckweed to become a prominent feature in restored reaches, and the high denitrification rates associated with duckweed shown in this study, the potential exists for duckweed to exert strong controls on N dynamics in restored reaches during low-flow periods, and will be an important area of investigation in future work.
The model developed for this study provides a validated tool that may be used to help inform sustainable management strategies in restored and natural streams with abundant duckweed biomass. Harvesting of duckweed biomass is common in wastewater and wetland treatment systems [27,[29][30][31]34,36,40], and has indicated that periodic, planned duckweed harvest could improve the overall N uptake by reducing the biomass periodically to allow for rapid regrowth of duckweed mats. This would also permanently remove the assimilated N from the stream, reducing the loss of living organic matter to receiving waterbodies during storm events that cause catastrophic scour. Harvested duckweed also has the potential to be a feed supplement, as highlighted in Körner et al. [29] and Cheng and Stomp [79]. Although the harvest of duckweed would likely be more difficult in a stream than wastewater ponds designed for surface skimming, the restoration process for impacted streams often includes widening the channel and connecting the stream with its floodplain, which increases the surface area and, potentially, accessibility for harvest. Modeling results may help to inform optimum timing of harvesting. For instance, in our study, we found DIN removal percentages were lowest in late fall through early spring (0.04-4.75 percent). Harvesting may therefore have the greatest impact in late spring, when the N supply is abundant, and other conditions are non-rate-limiting (e.g., temperature and light availability). The numerical model used in this study could provide a tool for site-specific harvest scheduling plans based on anticipated environmental conditions.

Model Limitations and Future Work
Notwithstanding the important findings of this study to inform vegetation impacts on N cycling in bedrock-controlled karst agroecosystem streams, we observed broad ranges in our uncertainty analysis which reflects the infrequent measurements used for model evaluation purposes and suggests a need for improved databases for model evaluation. In particular, we were unable to properly constrain diel fluctuations in CDIN, which can be seen to be substantial in both the sensitivity analysis ( Figure 3) and the calibrated and validation model output (Figure 4), showing fluctuations as high as roughly 1 mg L −1 . The use of high-frequency water quality sensors and high-resolution data, though, has been seen to provide important insights into diel fluctuations of nitrate (see [80] and references within) as well as provide estimates of primary productivity and gross primary production using coupled dissolved oxygen (DO) and nitrate diel variability [81,82]. Given the large diel fluctuations in these modeling results, the inclusion of high-resolution data and methodologies could better constrain these CDIN oscillations, which in turn could provide improved estimates of biotic uptake and removal and overall DIN budget improvements.
Further, while the modeling structure presented here offers a parsimonious representation of DIN, which is applicable to our system due to low levels of total ammoniacal nitrogen (TAN) yearround, explicit modeling of TAN and nitrate could be an important consideration elsewhere. While both TAN and nitrate are biologically available, aquatic vegetation has been shown to prefer uptake of ammonium when both species are abundant [83]. In a duckweed pond N transformation model, Peng et al. [30] considers ammoniacal-N and nitrate separately, and has different half-saturation coefficients for the uptake of each individual N pool. Surface streams with high organic runoff from pastures or point source contributions of wastewater would benefit from explicit consideration of each N form, which could aid in more appropriate estimates of biomass growth, overall N removal, and more realistic loading estimates. This could be of particular importance when determining management strategies for consistent inputs or estimating the effects of accidental overflow or wastewater system failure.

Conclusions
The model developed in this study offers a parsimonious approach for inclusion of algae and floating aquatic macrophytes into stream water quality modeling and provides an alternative approach to the estimation of denitrification rates in streams with negligible benthic sediments (i.e., streams with bedrock control). The results indicate the potential of this model to capture the controls on N dynamics in karst agroecosystem streams, offers insights into the seasonality of competing controls on DIN concentrations, provides inferences into changing N dynamics under differing hydrologic or climatic regimes, as well as highlights the possibility for duckweed to control overall denitrification rates. The modeling results also suggest these streams have the capacity to perform denitrification on the same order of magnitude as other agricultural streams, often considered to have higher denitrification potential due to their extensive benthic sediments and hyporheic exchange, highlighting their potential importance to watershed and regional N fate and transport studies.
Supplementary Materials: The following are available online at www.mdpi.com/2073-4441/12/9/2458/s1, Table  S1: Seasonal average DIN removal rates, DIN vegetation uptake rates, denitrification rates, and regeneration rates for median (minimum-maximum) DIN values across the 47 posterior solutions for the calibrated model, Table S2: Seasonal average DIN assimilation and DIN denitrification rates in each biotic pool for median