Introduction

Wide application of industrially produced nitrogen fertilizer for agriculture has contributed pre-eminently to the doubling of nitrogen fluxes over the last 60 years (Galloway and others 2008). The quantification and management of nitrogen leaching from agricultural land remain as the main challenges because of diffuse transport pathways to water bodies (Alexander and others 2000). Indeed, the doubling of reactive nitrogen has tremendous direct effects on surface water eutrophication in lakes and streams as well as in estuarine (for example, Brittany) and marine (for example, North Sea and Baltic) environments. These environmental impacts led Rockstrom and others (2009) to consider this nitrogen increase as one of the two most pressing global environmental concerns together with the loss of biodiversity (Seitzinger 1988; Vitousek and others 1997; Seitzinger and others 2006; Diaz and Rosenberg 2008; Galloway and others 2008).

Estimations of nitrogen removal capacity based on budget balances have reported that about 30 to 40% of N input in river catchments is lost by denitrification (Seitzinger and others 2006): the microbially facilitated reduction of nitrate and other nitrogen oxides to dinitrogen. This microbial process uses nitrate as an electron acceptor during the oxidation of the organic matter (the electron donor). Denitrification is a well-known pathway and its proximal primary environmental drivers (that is, anoxia, presence of bioavailable organic carbon, and in situ nitrate concentrations) are well defined. Yet, direct measurement of in situ denitrification is very difficult due to the high spatio-temporal variability of its environmental drivers (Groffman and others 2006).

It has been known for several decades that denitrification is the main process removing nitrate in riparian zones, which can buffer upslope input of nitrate and limit aquatic ecosystem pollution (see reviews by Burt and others 2010; Ranalli and Macalady 2010). These particular landscape features possess all the characteristics to potentially host denitrification activity, that is, (i) anoxia in soil during stream flood events or groundwater rise, (ii) the presence of high organic matter concentration in soils generated by very productive riparian vegetation, and (iii) the potential nitrate input from surface and subsurface flow from adjacent agricultural lands (see review by Décamps and others 2004). Many case studies have been conducted during the last 30 years that confirmed the potential nitrogen buffering capacity of riparian zones and have reported removal rates of up to 90% (for example, Vidon and Hill 2004a). But these studies have also underlined substantial spatial variability caused by local geomorphic heterogeneity (Sabater and others 2003), local hydrogeological conditions (Vidon and Hill 2004a, b; Duval and Hill 2006), and temporal variability caused by stochastic groundwater table fluctuations and flood events (Burt and others 2002) or stream stage fluctuations. The large local variability in nitrogen buffering capacity makes robust extrapolation of site-specific nitrogen retention to the landscape and/or drainage basin level, a major research challenge. Similar challenges are faced in understanding the role of the riparian zone for other substances, including DOC (Grabs and others 2012). Model-based estimations of the potential riparian nitrate removal capacity under optimal conditions for denitrification to operate reveal that riparian zones could not contribute more than 15–20% of the removal of total nitrogen fluxes in agricultural catchments (Groffman and others 2006; Montreuil and Mérot 2006; Seitzinger and others 2006). Indeed, riparian zones are not the only landscape features in which nitrogen is removed very efficiently by denitrification. For instance, hyporheic zones (Hill and others 1998; Duval and Hill 2006; Krause and others 2009, 2013; Trauth and others 2014), ditches, potholes, upland slope breaks (Clément and others 2003), hedgerows (Viaud and others 2001), are other landscape structures that can promote microbial denitrification, at least seasonally. Knowledge about the extent of their contribution to overall nitrate turnover is critical for management. However, the identification of these hotspots and their role in a landscape context (Vidon and others 2010) requires further research to evaluate the overall removal capacity of these landscape features mentioned above at the drainage basin scale.

Existing Approaches

Mass Balance Method

It consists of comparing nitrogen input within the catchment with nitrogen output at the outlet from the drainage basin. Since Omernik’s studies (1976), several nitrogen input-output studies in large drainage basins, that is, greater than 500 km2, have shown a positive relationship between the percentage of agricultural land and fluxes of nitrogen at the outlet. The European Water Framework Directive (EU 2000), which stipulates monthly measurement of nitrate concentration at least in drainage basins of 100 km2 or larger promotes such an approach. The interpretation of such relationships in these mass balance approaches is difficult because it is based on a black-box approach, which does not provide much information on the actual removal capacity of the specific landscape features within the wider drainage basin. In fact, denitrification is simply estimated as the difference between nitrogen input and output (Groffman and others 2006; Seitzinger and others 2006) and does not provide a suitable spatio-temporal framework.

Proxies of Denitrification

These proxies include potential denitrification activity and potential wetted area, measured by biogeochemical and geographic methods, respectively, to quantify the spatio-temporal variability of denitrification. Potential denitrification activity (for a review of methods see Smith and Tiedje 1979; Burt and others 2010) provides a useful indication of the denitrifying enzyme activity present at the sampling time. Potential denitrification can reflect immediate environmental changes affecting denitrification activity such as soil moisture and aeration, for instance. It provides also an evaluation of the potential physiological denitrification capacity of the existing denitrifying community; as such it represents a more robust picture of current environmental conditions. Yet, in the past, several studies have over-interpreted the results of this method by confounding potential and actual denitrification activity. Indeed, potential denitrification can provide information on the potential ability of the existing community but not on the real rates because they also depend on the nitrate availability/input to the site which is primarily influenced by the local hydrogeological context (Sabater and others 2003).

At the drainage basin scale, potentially wetted areas or wetlands, where anoxic conditions and organic carbon accumulation prevail, have been often used as proxy for denitrification. Several types of models based on topographic indices, that is, the upslope contributing area per unit contour and the local slope angle (Beven and Kirkby 1979; Moussa 2009) have been used to evaluate the areal extent of potentially wetted areas within a drainage basin (Mérot and others 2003). This modeling approach is very useful for quantifying the potential for wetlands to develop within catchments (Poggio and Soille 2011). Although this approach provides a qualitative indicator for potential denitrification zones, it does not support any quantitative assumptions of nitrogen removal because it does not provide information about nitrate input or residence time. Indeed, coincidence of flow paths transporting nitrate into reactive anoxic zones are a prerequisite for denitrification to occur (McClain and others 2003). Several attempts to use a hydrogeomorphic unit approach based on Brinson and others’ work (1984) did not really improve our capacity to assess denitrification rates because of the lack of the hydrological component necessary to quantify nitrate input to the anoxic sites. Moreover, residence times of water in soils or sediments need to be quantified to move from potential denitrification zones to actual ones. Indeed, several in situ studies have demonstrated that increasing the time water resides within sediment or soil increased nitrogen retention and removal (Pinay and others 2002, 2007; Zarnetske and others 2011). It may seem surprising in a world where the nitrogen fluxes have doubled over the last 60 years, but in most cases, the key limiting factors for denitrification are (i) the nitrate input to the active sites and (ii) the time a nitrate molecule is being exposed to denitrifying conditions, that is, its exposure time (Oldham and others 2013). Hence, the local hydrogeological context is the key denitrification driver (Sabater and others 2003). Therefore, the move from potential denitrification to real denitrification activity evaluation at the drainage basin level requires the quantification of nitrate input into these potential denitrification sites (McClain and others 2003). At the landscape scale, the riparian buffer ratio as defined by McGlynn and Seibert (2003) in the hydrological landscape analysis, or the length of contact between dry and wet areas, including upslope/wetland lengths, should be a good proxy for nitrate input to a potential denitrification site. Indeed, many studies on the role of riparian zones in mitigating diffuse nitrate pollution revealed that when removal occurred, it was within the first few meters of the riparian zone seen from the upslope side (Clément and others 2003; Sabater and others 2003). The advantage of such a proxy is that it can include other types of wet/dry interfaces and can be easily quantified at large scales using remote sensing imagery.

Hydrologic Landscape Analysis

Hillslope hydrology is predominantly controlled by topography in catchments with shallow soil and poorly permeable bedrock, which are typical for some regions. McGlynn and Seibert (2003) examined the variability in, and controls on, hillslope inputs to stream networks and the potential for riparian zones to regulate hillslope inputs and thereby both quantitatively and qualitatively buffer, or modify, stream responses to hillslope hydrology. The ratio of riparian zone storage to hillslope inputs was the most important plot-scale measure of the buffering capacity of the riparian zone. Clearly, catchments will differ in the degree to which riparian zones buffer the delivery of water from hillslopes to streams, thereby affecting the amount, timing, and quality of hillslope water inputs expressed in streamflow (McGlynn and others 1999; McGlynn and McDonnell 2003). Hydrologic dynamics within and connections between landscape units partially control the sources, flowpaths, amount, and age of water exiting the catchment through each segment of the riparian zone. Each landscape unit type has characteristic hydrologic and geomorphic attributes that can be assessed through field investigations and topographic analysis of emergent patterns and connections between landscape units (landscape organization). For instance, Meybeck and Moatar (2012) propose a typology of catchment nutrient transfer responses based on relationships between the log of nutrient concentration and log of discharge, as well as flux estimations and their uncertainties (Moatar and others 2013). Hydrological landscape analysis (HLA) techniques enable such upscaling as demonstrated in a diverse variety of catchments in published research. HLA can, for instance, provide a quantification of the topographic control on water age or residence time distributions. Often it is expected that mean water age in runoff increases with catchment area but this could not be confirmed by either McGlynn and others (2004, 2005). McGlynn and others (2004) found the median subcatchment area to be correlated with mean residence time, whereas McGuire and others (2005) found the median flowpath length toward the stream network divided by its gradient to be best correlated with residence times. Duncan and others (2013) used HLA, in combination with local field measurements of N2 gas production, as a proxy for denitrification to upscale denitrification to a whole forested catchment. They found that denitrification hotspots in topographic hollows caused by riparian microtopography, as shown earlier by Frei and others (2012), has a significant influence on catchment-scale denitrification.

Model-Based Upscaling of Local Nitrate Removal Capacities

Deterministic, spatially detailed catchment models (for example, SWAT (Arnold and others 1998), SWIM (Krysanova and others 1998), INCA (Whitehead and others 1998), HBV-N (Arheimer and Brandt 1998) are applied at large-scale drainage basins (up to several 100,000 km2) to quantify nitrogen uptake and removal in a spatially distributed way (Sahu and Gu 2009; Lam and others 2010; Ficklin and others 2013; Poudel and others 2013). Yet, the spatial and temporal heterogeneities of the denitrification process in particular in riparian corridors, characterized by the spatially and temporally dynamic exchange between groundwater and surface water question their large-scale assessment reliability. As these models are designed to assess water quality at the scale of entire catchments, their underlying concepts are based on significant spatial averaging of catchment properties (for example, at the subcatchment scale) and processes (that is, they are semi-distributed) and simplified descriptions of the connected groundwater system (for example, INCA, SWAT). Furthermore, system dynamics are averaged in time by typically using relatively long time steps. Similarly, the spatial discretization of mesoscale catchment models such as SWAT or SWIM is organized in terms of hydrological response units (HRU’s) that are assumed to represent homogeneous landscape units of similar hydrological behavior. As model discretization and delineation of response units at large scales are considered to be static, the behavior at interface zones such as riparian corridors and wetlands (Hattermann and others 2006) are usually not implemented in a dynamic way. Due to these limitations such models are usually applied for predicting average system behavior (for example, monthly or annual loads over longer time periods) and frequently assess the impacts of climate and land-use changes, but are limited in their capacity to represent dynamic feedback functions at interfaces such as riparian zones, which would be needed, for example, to evaluate local measures to improve water quality. Process-rate variability leads to lack of precision in such models. Moreover, bias in denitrification estimation often arises when there is co-variation between denitrification activity in the riparian zone and nitrogen inputs from upslope soils to these riparian zones. This bias leads to serious inaccuracies.

The Existing Frameworks

Three main frameworks have been proposed to date to formalize the understanding of the forces driving denitrification and the upscaling of nitrogen removal capacity to the landscape and drainage basin organization levels (Figures 1, 2, 3). These three frameworks, that is, the classification of denitrification areas according to whether nitrate transfer is via diffusion or advection (Seitzinger and others 2006), the “biogeochemical hot spot” concept (McClain and others 2003), and the application of the Damkhöler ratio to evaluate riparian zone nitrogen removal efficiency (Ocampo and others 2006), provide complementary elements to evaluate nitrogen removal capacity in drainage basins.

Figure 1
figure 1

Schematic representation of the existing frameworks related to denitrification appraisal. A Diffusion of nitrate into the denitrification hotspot (for example, pond or lake sediment) or B advection (for example, riparian or hyporheic zone).

Full size image
Figure 2
figure 2

Schematic representation of the existing frameworks related to denitrification appraisal. The hotspot concept with A input of nitrate into the denitrification hotspot (for example, riparian zone), and B convergence of nitrate and organic carbon into the denitrification hotspot (for example, hyporheic zone). Adapted from McClain and others (2003).

Full size image
Figure 3
figure 3

The Damkhöler number as a function of nitrate buffer zone capacity of potential denitrification hotspots. Adapted from Ocampo and others (2006).

Full size image

Diffusion Versus Advection

In a review paper, Seitzinger and others (2006) classified sites along a continuum ranging from terrestrial to freshwater and marine environment as a function of the nitrate delivery mode to the denitrification zone. The classification of sites as a function of the delivery of nitrate, that is, by diffusion or advection, can be related to particular types of landscape/waterscape features (Figure 1). For instance, lake and pond sediments can be classified as zones where nitrate is delivered to the denitrification areas mostly by diffusion, whereas riparian or hyporheic zones are mostly characterized by a nitrate delivery by advection. The diffusion/advection framework touches upon the rate of delivery of nitrate, and as such, adds a temporal component. Indeed, diffusion rate is often slower than denitrification rate measured in lakes, pond, or stream sediments (Sébilo and others 2003). In such cases, the rate of nitrate delivery is the limiting factor for denitrification. On the contrary, nitrate advection can present a wide range of delivery rates, depending on factors driving the local water flow rate, that is, slope, matrix structure and texture, hydraulic conductivity. Therefore, in the advection case, either nitrate delivery or rate of activity of the denitrifiers can be the limiting factor for denitrification (Sébilo and others 2003).

Apparently, the time that is available for denitrification is completely different for the two cases. It therefore seems appropriate to define characteristic time scales at which either the advective or the diffusive system is being exposed to conditions favorable for denitrifying activities (Oldham and others 2013). Appropriate conditions for denitrification imply anoxic conditions.

Under advective delivery mode the exposure time to anoxic conditions, \( \tau_{{E , {\text{anoxic}}}}^{\text{adv}} \) is defined as

$$ \tau_{{E , {\text{anoxic}}}}^{\text{adv}} = \frac{{L_{\text{anoxic}} }}{v}, $$

\( \tau_{{E , {\text{anoxic}}}}^{\text{adv}} \) is a mean value and there is a probability distribution attached to it, where L anoxic is the length scale over which transport processes are operating under anoxic conditions (m), v is the mean ground water flow velocity (for example, in m d−1).

Alternatively, the exposure time scale to anoxia under diffusive conditions is given as

$$ \tau_{{E , {\text{anoxic}}}}^{\text{diff}} = \frac{{(L_{\text{diff}} )^{2} }}{{D_{\text{eff}} }}, $$

D eff is the effective diffusion coefficient [m2 s−1] and L diff [m] is the length scale of the system across which diffusion occurs.

Putting this concept into a more quantitative framework, diffusion versus advection can be addressed in terms of the Peclet number under the assumption of identical advection and diffusion scales (L diff = L anoxic) (Oldham and others 2013). The non-dimensional Peclet number, Pe, provides the balance between advection and diffusion under hydrological connectivity, as

$$ Pe = \frac{{\tau_{{E , {\text{anoxic}}}}^{\text{diff}} }}{{\tau_{{E , {\text{anoxic}}}}^{\text{adv}} }} = \frac{{v \cdot L_{\text{anoxic}} }}{{D_{\text{eff}} }}. $$

Typically, Pe numbers ≫ 1 are indicative of relatively slower diffusion than advection and thus advective transport, while Pe numbers ≪ 1 reflect relatively slower advection to diffusion and thus diffusive transport. Critical zones in terms of transport are systems where Pe is on the order of 1. Here diffusive mixing becomes relevant. Note that Pe balances with scale, that is, the smaller the scale the more the systems tend to become diffusion controlled.

At small scales, diffusion is a very efficient mixing process that tends to homogenize concentrations. With classical values of diffusion coefficients for nitrate in water of 2 10−9 m2/s, the typical diffusion times across pores of 100 μm to 1 mm are of the order of 5 s to 8 min. Advection transports solutes over long distances more quickly delivering them to potential turnover sites where diffusion can mix them with other reactants.

The characteristic mixing scale is much enlarged by hydrodynamic dispersion processes (Bear 1973). Hydrodynamic dispersion comes from the pore-scale and interpore-scale velocity fluctuations and is generally modeled as an effective diffusion process under Fick’s law. The equivalent dispersion coefficient D is proportional to the velocity v and to the dispersivity α. Taking the Peclet number with respect to dispersion rather than to diffusion leads to Pe = L anoxic/α. The Peclet number can also be interpreted as the ratio of the characteristic advection scale to the dispersivity (characteristic dispersion scale). As dispersivity ranges from centimeters to hundreds of meters depending on the geological material and of the scale itself (Gelhar and others 1992), the critical scale over which mixing equilibrates with transport is much larger than the pore scale and can reach meters to decameters.

Biogeochemical Hotspot

The second framework is the biogeochemical hotspot concept (Figure 2). It stipulates that biogeochemical hotspots are “areas (or patches) that show disproportionately high reaction rates relative to the surrounding area (or matrix)” (McClain and others 2003). These hotspots occur at the convergence of hydrological flowpaths carrying complementary reactants, or where a flowpath carries one reactant into a substrate containing the other reactant. The convergence of reactants can be found, for instance, in stream hyporheic zones where nitrate transported by groundwater upwelling merges with organic carbon provided by surface water down-welling into the hyporheic zone (Krause and others 2009, 2013; Trauth and others 2014). Nitrate being transported by surface or subsurface flow into the riparian zone represents the second case, that is, the nitrate reactant being carried into a zone where the other conditions necessary for denitrification (availability of organic carbon and anoxia) are present. Current hyporheic zone modeling approaches focus on a hydrologically controlled spatial extent of the hyporheic zone, but with few exceptions (Zarnetske 2011; Bardini and others 2012; Trauth and others 2014), do not take into account that electron donors and acceptors might be delivered via different pathways. Hence, from a kinetic point of view, biogeochemical hotspots are zones where turnover rates are high and can be characterized by short characteristic reaction time scales τ reaction. Characteristic reaction time scales can be derived from any kinetic rate law as the time that is required to establish a certain turnover rate. The most common approach is to use the inverse of the first-order rate intrinsic constant k [for example, 1/d], that is, τ reaction = 1/k.

The rate constant is the intrinsic constant valid at the molecular scale or the scale at which microbes are operating. This rate constant is scale invariant and should not be confounded with effective rate constants derived from fitting, for example, of the advection dispersion equation with a reaction term calculated from slurry experiments.

Damköhler Number

The third framework is related to the evaluation of the nitrate removal efficiency of riparian zones using the Damköhler number (Figure 3). The dimensionless Damköhler number is the ratio between the rate of transport (rate of nitrate input to the site) and the rate of reaction (denitrification in the site). This can be defined as the ratio between the exposure time to anoxic conditions and the reaction time scale of nitrate:

$$ Da_{{{\text{NO}}_{ 3}^{ - } }} = \frac{{\tau_{{E,{\text{anoxic}}}} }}{{\tau_{\text{reaction}} }}. $$

Thus, it is a measure of the competition between transport and reaction processes and it actually combines the two frameworks discussed above. Large Damköhler ratios

(\( Da_{{{\text{NO}}_{ 3}^{ - } }} \gg 1 \)) indicate reaction times much smaller than exposure times to anoxic conditions. The reaction has ample time to occur within the allocated exposure time and nitrate is efficiently removed. On the contrary, when the characteristic reaction time is much larger than the characteristic exposure time (\( Da_{{{\text{NO}}_{ 3}^{ - } }} \ll 1 \)), the reaction does not have time to occur and nitrate is not removed.

Ocampo and others (2006) used the Damköhler number to quantify the relative importance of transport versus reaction in the attenuation of nitrate concentrations within riparian zones. Their study was based on published results from several sites all over the world. They predicted an increase in nitrate removal rate \( R_{{\text{NO}}_{ 3}}\) from riparian zones with increasing Damköhler numbers. We have evaluated their data and obtained an exponential relationship between these two variables.

$$ R_{{{\text{NO}}_{ 3} }} [\% ] = 75 \cdot (1 - \exp ( - 0.4Da)) - 25. $$

Marzadri and others (2012), Zarnetske and others (2012a, b), and Briggs and others (2013) used a similar approach based on a Damköhler number for oxygen consumption to delineate zones of net denitrification from zones of net nitrification in the hyporheic zone.

These three frameworks are very complementary in their quest to appraise the overall denitrification activity or nitrate removal in landscapes and drainage basins; but none of them fully reached that goal, because of the daunting challenge of upscaling heterogeneous processes both in space and time (Groffman and others 2006). Yet, identification of potential biogeochemical hotspots, modes of local delivery, that is, diffusion versus advection and transport rate versus activity rate provide a sound basis for building up a new approach.

Gu and others (2007) have highlighted the interdependence of Pe and Da because inherent to Da numbers is also the flow rate. Nitrate removal rates are at a minimum for both low Pe and low Da, whereas they are at their maximum rates for Pe values greater than 10 and Da greater than 25 in their example. The effect of flow rate has an opposite effect on both numbers. For a given reaction time scale, removal rates are decreasing if the flow rates are too high (exposure time should be too short, Da number decreasing) but on the other hand, Pe numbers are increasing. Hence a combination of low turnover rate, low Pe (diffusive conditions), and low Da can only be achieved if the characteristic reaction time is low.

Toward an Integration of Frameworks

We propose to combine the concepts of residence time distribution (RTD) and Da/Pe relationships, which are strongly interdependent. Pore to Darcy scale hydrodynamic dispersion progressively transforms the Dirac distribution into more complex distributions, which span a large range of residence times, as described for example, by an exponential, power law or inverse Gaussian distribution (Table 1) (Bear 1973). Figure 4 presents the differences between the well-peaked Dirac distribution (blue) and the more extended exponential and inverse Gaussian distributions (red and black).

Table 1 Probability Distribution Functions Derived for the Damköhler Number Directly Derived from the Residence Time Distributions with Parameters µ for the Mean Damköhler Number and Pe for the Peclet Number if Relevant as Illustrated in Figure 5
Full size table
Figure 4
figure 4

A Illustration of three theoretical Damköhler density distributions p(Da), and B cumulative distributions (P(X > Da), with the same mean equal to 1. The Peclet number corresponding to the inverse Gaussian is Pe = 0.5.

Full size image

This requires tools for determining RTD in anoxic zones and also the spatial and temporal variability of Da and Pe. Both Da and Pe are functions of RTDs, that is, hydrological functions (Da via the exposure time and Pe via the ratio between advective and diffusive transport). But Da is also a function of the characteristic reaction time, that is, of biogeochemical functions. For kinetically controlled reactivity, it can be compared even more straightforwardly to the characteristic reaction time, or even to some distributed first-order reaction times if needed, by convolution of exposure time and reaction time distributions. For this, it is critical to define the intrinsic rate constants, their dependence on electron donors (that is, DOC), and their spatial variation (various species).

In hot-spots environments where bioavailable organic matter exists in excess, for example, peatlands, wetlands, and relict C deposits (Gurwick and others 2008), denitrification can be approximated to be a pseudo zero-order reaction with respect to organic matter (Reddy and Patrick 1984), which implies that

$$ \frac{{d{\kern 1pt} \,c({\text{NO}}_{3}^{ - } )}}{dt} = - k^{\prime}c({\text{NO}}_{3}^{ - } ), $$

where k′ combines the 2nd order rate constant for denitrification k den and the concentration of dissolved organic matter. Reddy and others (1982) determined mean values for k′ of 0.00147 ± 25% d−1 (µg C mL−1)−1 for mineral soils and 0.00155 ± 65% d−1 (µg C mL−1)−1 for organic soils.

Actually, good data are still lacking. The experiments by Reddy and others (1984) were performed at high initial nitrate concentration (200 mg/L) and we need to be careful to insert these k den values in the equation above. k den is certainly a good approximation and it is to our knowledge the best that we have. So what is currently missing is a reliable study to establish the rate constants for denitrification and possibly also Michaelis–Menten parameters, both under low and high nitrate concentrations as well as DOC concentration and quality. It may well be that with complex soil organic matter it is the hydrolysis step to generate low molecular weight organic compounds required for denitrification that is rate limiting.

DOC exhibits a wide range of concentrations in various systems. For instance Knorr (2013) determined an average DOC concentration of 25 mg L−1 in a minerotrophic fen in the Fichtelgebirge (Southern Germany), while average DOC values of 2 mg L−1 were determined in the hyporheic zone of a mountain stream in Western Oregon (Zarnetske and others 2011), which is probably at the lower end of DOC values in shallow subsurface environments. Hence, k′ values may vary between 0.002 and 0.04 d−1 and even higher values at higher DOC concentrations. Therefore, τ reaction (=1/k′) ranges between 25 days (or even shorter) and 500 days. This analysis illustrates that it is mainly the DOC concentration that is driving the reaction time scale rather than the rate constant, given the constant value of k den cited above.

If we assume DOC concentrations to be rather constant for a certain system (for example, the hyporheic zone discussed above) and also fluctuations of DOC at least under base-flow conditions not to be high (Strohmeier and others 2013), then, τ reaction can be assumed to be a roughly constant parameter for that system. Hence, the Damköhler number in this system (and the subsequent nitrate removal capacity) will increase when the transport time, that is, residence time, increases. If we allow for DOC to vary, but can predict this variability as has been successfully done in many environments (for example, Winterdahl and others 2011, 2014); then, there will be more to the temporal dimension than just the change in flow rates.

The broad variability of the nitrate removal capacity also comes from the large variability of the hydrologic and geologic conditions. The RTD in the different hydrologic compartments (unsaturated zone, aquifer, hyporheic, and riparian zones) are not well-peaked Dirac-like distributions but broader distributions that span over large time ranges (Figure 4). Flow heterogeneities from the pore scale to the watershed scale enhance dispersive and mixing processes as larger scales are considered. Additional local-scale trapping, diverting mechanisms in slow diffusive zones, chemical sorption, and exchanges between advective and diffusive zones critically enhance the solute spreading in all the spatial directions and the variability of the travel times (Berkowitz and others 2000; Carrera and others 1998; Haggerty and Gorelick 1995).

Such local-scale transport processes may be the key in effective nitrate removal capacities as advection dominated zones may coexist with diffusion dominated zones, differing critically in their residence times and potentially also in their chemical conditions (more anoxic conditions in the diffusive highly denitrifying zones than in the main advection paths). Practically, these mechanisms decouple flow from exposure times in anoxic zones. Both flows and exposure times may remain simultaneously high, while in more advective conditions they are inversely correlated. Nitrate delivery may remain high and nitrate can have long exposure times because of the abovementioned local diffusion processes. Initial synthetic tests indeed indicate that the hydrogeological complexity conditions effective denitrification rates (Green and others 2010, 2014).

Even if multiple scales, multiple domains, and multiple processes intervene in nitrate degradation, the exposure time in anoxic zones distribution remains a good integrative indicator of the transport-induced constraints on reactivity (de Dreuzy and others 2013). For kinetically controlled reactivity, it can be even more straightforwardly compared to the characteristic reaction time, or even to some distributed reaction times if needed, by simple convolution (Böhlke 2002). The ratio of the exposure time to the reaction time (Damköhler number) becomes also broadly distributed and can be further considered to estimate potential denitrification. When assuming a unique denitrification rate, the probability density function of the Damköhler number directly follows from the probability density function of the RTD described by the Dirac, exponential, and inverse Gaussian distributions given in Table 1 and illustrated in Figure 4A for the same mean Damköhler value of 1. It should be noted that the shape of the distribution may have a determinant effect on the denitrification capacity as expressed by the cumulative distribution of the Damköhler number (Figure 4B) even for distributions having been calibrated on the same data (Green and others 2014). This approach requires the a priori choice of a RTD shape but it does not require large amounts of data. Residence time distributions can be calibrated based on data, or derived from flow models and further used in a predictive manner. Although limited so far (Eberts and others 2012; Green and others 2014), modeling studies show that beyond some consistent choice of the distribution, the most limiting factor is the availability and the accuracy of tracer data. These modeling studies are important as the RTD is never directly accessible from field data but always inferred with some assumptions. Numerical simulations on realistic aquifer structures offer a sound alternative to test beforehand the sensitivity of denitrification to different geological structures and hydrological conditions.

Priority Research Questions

We have already underlined that up scaling from local assessments from each individual landscape feature is dubious because of the high local spatio-temporal variability of denitrification. But how can plot-scale measurements of the hydrologic processes that link hillslopes, riparian areas, and streams be scaled up to entire catchments? How can we transfer small catchment scale understanding to larger portions of the landscape or other catchments? How can we quantify the effects of different spatial patterns (for example, the distribution of riparian wetlands) within catchments? These questions highlight our current inability to transfer understanding of hydrological processes studied at the plot or hillslope and reach scale to small catchments (Figure 5). These questions beg integrated, multi-scale approaches that combine “landscape level” topographic analysis, process-based field investigations, and catchment-scale integration to identify the factors controlling hydrologic connectivity between source areas with different hydrochemical properties and the flow paths that link source areas to streams. The approach consists of a combination of hydrological landscape analysis to derive the RTD and an assessment of the reaction time scales in different landscape elements (for example, Jencso and others 2010). To these ends several tools are required.

Figure 5
figure 5

Schematic view of the small drainage basin and hypothetical log–log distribution (P) of nitrate time residence (years) in A the insaturated zone where percolation/retention patterns occur with a long tail distribution corresponding to biological retention and poorly hydrologically-connected porosity; B the shallow aquifer with dispersion being a function of the complexity of the landscape structure; and C in the deep aquifer where dilution occurs.

Full size image

In small drainage basins, that is, less than 10 km2, at any given point in time, the relationship between nitrogen fluxes at the outlet of the catchment and the percentage of land-use, a non-spatially explicit variable, does not hold any more (Burt and Pinay 2005; Strayer and others 2003). We hypothesize that large spatial variability of nutrient fluxes in small drainage basins (Bishop and others 2008) reveal the subtle changes in land cover and land uses as well as the importance of their spatial arrangements. Moreover this absence of relationship between nitrogen fluxes and the percentage of land-use reveals also how denitrification efficiency of the spatial architecture changes over time, both in the short term with runoff rates, and seasonally with biological activity, as well as in response to antecedent conditions. Indeed, it is not the total percentage of a given land cover which matters at these small scales but the local arrangement and connectivity of the different land covers and their interfaces, including riparian zones, within the landscape matrix at a given point in time. Therefore, this approach needs to be undertaken in small drainage basins, that is, stream order ≤3 which roughly corresponds to catchments between 1 and 10 km2. This is an important issue because this approach proposes as an underlying statement that landscape structure and land cover arrangement are the key parameters to evaluate the overall nitrate removal capacity at the drainage basin level (Laudon and others 2011).

Although detailed models coupling locally accurate physical and chemical processes can be set up on intensively studied local sites, such models are usually not tractable at small stream order drainage basin scales. We therefore propose an alternative approach that is based on lower data input but a higher degree of generalization and transdisciplinary synthesis from well-studied sites. Nitrate removal capacity is not expected to be highly dynamic; therefore, it will be amenable to estimation based on key proxies that determine the distribution of exposure time in anoxic zones and identify relevant nitrate removal drivers such as the ratio of unsaturated to saturated zones, the extension of the riparian zones, the average rate of quick runoff to slower infiltration and transfer in aquifers, the mixing capacity and chemical inventory (for example, pyrite and organic matter content) of the aquifers, and biological descriptors of nitrate removal efficiency. All these proxies should be identifiable from widely available field data and insights gained from detailed process studies. Such a proxy-based assessment of nitrate removal capacity is consistent with the increasing availability of spatially and temporally resolved data on relevant topographical, hydrological, and geochemical data. We argue that this approach of embedding the expertise and insight gained from more in-depth analysis of a broad diversity of sites should be more informative than purely correlative and data-mining analyses. Inter-comparison, benchmarking and in-depth studies conducted on a wide variety of sites are pre-requisites for this approach. For an improved synthesis it will also be advantageous to fully capitalize on the growing mathematical and numerical capacities for scenario simulations and modeling experiments.

Several proxies and tracers of processes that can be used to calibrate, validate, and develop this new approach are presented below. We believe that their use, in combination with hydrological landscape analysis as a means of predicting the denitrification rate at different flow rates through the biogeochemical environments that comprise the catchment, will allow improved models of how water flows through the landscape to continuously change the activity and connectivity of denitrification hotspots (Figure 5).

Tools and Methods Available

Remote Sensing and Embarked Imagery to Evaluate Landscape Features

At the landscape scale, a wide range of remote sensing data can be used to reveal the local arrangement of the different land covers and their interfaces, including riparian zones, within the landscape matrix (Rogan and Chen 2004; Goetz 2006). Indeed, satellite and aerial imagery can contribute to determine sites that are suitable for denitrification in given landscape arrangements in the drainage basins. The scale of perception of the drainage network influences the flowpaths and land/water interface length (Gold and others 2001), especially when analyzing them through data provided by remotely sensed imagery. Instruments such as airborne lasers light detection and ranging (LiDAR) provide innovative contributions to the detection and mapping of the drainage network as well as the wet/dry interface length at a fine scale. LiDAR data have been used to characterize the ground microtopography with centimeter accuracy, including ditches and streams, even under tree cover (James and others 2007) or herbaceous vegetation (Hopkinson and others 2005).

The accuracy of Digital Terrain Models (DTMs) generated with LiDAR data, which are most often constructed in a given spatial resolution according to the size of the study area, depends on the LIDAR point density. The drainage network can be automatically detected using form criteria derived from LiDAR data, either indirectly from the DTM (Pirotti and Tarolli 2010) or directly from the point cloud (Bailly and others 2008). Most studies conducted on the mapping of the drainage network with such data are focused on the identification of network elements, including small ditches and the historical drainage network (Werbrouck and others 2011) without characterizing them. However, information such as the width or depth of streams and ditches that are needed for biogeochemical functions like denitrification can also be determined with a high level of precision (Rapinel and others 2013). Bathymetry of the drainage network can be reconstructed from the slopes of the emerged banks, but this type of model requires many land surveys (Merwade and others 2008). Hence, the wet/dry interface length can be accurately mapped from DTMs generated from LiDAR data and integrated in models based on topographic indices to evaluate the surface of potentially wetted areas within a drainage basin. The spatio-temporal quantification of this proxy can be achieved using other active imaging sensors. Observations in near-surface aquifers may however be limited. Simple models should then be used to propose assumptions for their structure. It is especially the case for the weathered zone of crystalline basements. The free-surface aquifer is controlled by chemical weathering reactions and by percolation processes. For example, in the case of non-limiting chemical weathering, Rempe and Dietrich (2014) propose a “bottom-up” model of the weathered profile, which can be considered as a proxy of small drainage basin shallow aquifers. At equilibrium between erosion and uplift, the weathered zone is bounded above by the sediment export capacity from the hillslope to the river and below by the persistent saturation of the slowly drained fresh bedrock.

Process-Based Models of Hydrologically Mediated Turnover

Small-scale process-based flow and transport models are now a common tool to evaluate hydrologic dynamics and biogeochemical turnover at compartment interfaces in the landscape. Especially for the groundwater-surface water interface such models have been widely used to quantify groundwater-surface water exchange dynamics (for example, Frei and others 2012) and hyporheic flows in 2D (for example, Cardenas and Wilson 2007a, b) and 3D (for example, Tonina and Buffington 2007; Trauth and others 2013) as well as associated biogeochemical turnover (for example, Bardini and others 2012; Trauth and others 2014). These models are typically used as learning tools to derive a general mechanistic understanding of hydrologic dynamics and the resulting biogeochemical process patterns. For example, they can elucidate the dominant hydrologic controls for the development of biogeochemical hotspots in the riparian zone (Frei and others 2012) or biogeochemical zonations in meander bends (Boano and others 2010; Gomez and others 2012) and in the hyporheic zone (Marzadri and others 2011, 2012; Bardini and others 2012; Trauth and others 2014). Model results can replicate observed patterns and link them to the controlling processes and mechanisms. A more integral assessment of biogeochemical turnover may be achievable based on RTDs, which can be obtained from the hydrologic models (for example, Cardenas and others 2008). In a field study of an instream gravel bar (Zarnetske and others (2011), 2012a, b) demonstrated that nitrogen turnover in the hyporheic zone was clearly correlated with median residence time of the hyporheic water and that denitrification capacity of the gravel bar could be described in terms of Damköhler numbers. There is some evidence that such relationships may also hold for small catchments as long as they are structurally not too complex (van der Velde and others 2010, 2012) and predictable (Bishop and others 2011; Winterdahl and others 2014). A further complication is the transient nature of catchment-scale RTDs. However, van der Velde and others (2012) could show that for differently structured hill slopes, and presumably also for small catchments, transient RTDs can be collapsed into a unique, time invariant probability distribution function for the outflow of water of a specific age. If this function can be parameterized based on readily obtainable catchment characteristics, a process-based, integral description of catchment-scale denitrification might be possible without the need for a complex process-based model.

However, due to the inherent complexities of natural systems, especially at larger scales, and the nested scales involved, it is still challenging to adequately infuse results from these models into descriptions of matter fluxes and turnover at larger scales such as entire catchments (for example, in models like SWAT). Other, more conceptual modeling frameworks that link local processes to integral catchment response may emerge from the combination of enhanced local process understanding and improved integral descriptions of denitrification (for example, via proxies and refined landscape analysis) and may provide viable alternatives to currently existing models. To further develop, test, and refine such a new modeling framework we need adequate tracers and methods to describe the relevant key processes. The following sections provide examples of these tracers and methods.

Atmospheric Gases for Assessing Aquifer Residence Times

Atmospheric gases of anthropogenic origin including chlorofluorocarbons (CFC-11, CFC-12, CFC-113), sulfur hexafluoride (SF6), and Krypton-85 provide excellent tools for RTD in aquifers. Their atmospheric concentrations are well known and relatively uniform. Once in the saturated zone and without any exchange with the atmosphere, they provide excellent tracers of the aquifer circulations (Cook and Solomon 1995; Cook and others 1995; Newman and others 2010; Reilly and others 1994). They give essential information on the aquifer residence times from a few years up to around 50 years. Concentration measurements and interpretation require however some expertise. Because of low concentration levels, samples should avoid any contact with the atmosphere at any time of the sampling and analysis (Ayraud and others 2008). For SF6, geogenic production may occur in crystalline rocks and may reduce its applicability (Busenberg and Plummer 2000). Derivation of the RTD from the measured concentration also requires some a priori assumptions on the site-relevant dispersion of the residence times (Eberts and others 2012; Leray and others 2012; Massoudieh and Ginn 2011; Massoudieh and others 2012; Waugh and others 2003).

Heat as a Tracer to Elucidate Groundwater/Surface Exchange Processes

Using heat as a tracer has become a popular tool to characterize and quantify groundwater-surface water exchange processes (Constantz 2008). Natural temperature differences between ground and surface water can be used to qualitatively map exchange patterns (for example, Conant 2004), or to quantitatively invert flux rates based on numerical (for example, Constantz and others 2013) or simpler analytical solutions to the heat transport equation (for example, Hatch and others 2006; Schmidt and others 2006). Several analytical methods based on either an analysis of the amplitude damping and phase shift of transient temperature signals (Hatch and others 2006; Keery and others 2007) or a quasi-steady state evaluation of vertical temperature profiles (Schmidt and others 2006, 2007) have been proposed to quantify exchange fluxes. Despite their inherent simplifying assumptions these methods can provide reliable flux estimates over quite a broad range of field conditions (Anibas and others 2009; Lautz 2010; Munz and others 2011; Schornberg and others 2010) and allow a relatively simple means to characterize and quantify spatial exchange patterns (Anibas and others 2011). In addition easy-to-use software tools to estimate exchange fluxes based on these methods have become available (Swanson and Cardenas 2011; Gordon and others 2012). Spatially explicit flux estimates can in turn be linked to patterns of nitrogen turnover to identify and explain biogeochemical hotspots at the surface water-groundwater interface (Krause and others 2009, 2013; Briggs and others 2013). For even smaller-scale (<1 m) assessments of fluxes in streambeds active heat-pulse methods have been developed (Lewandowski and others 2011; Angermann and others 2012a), which allow quantification of small-scale spatial patterns of hyporheic exchange, which are superimposed onto the larger-scale river-aquifer exchange patterns (Bhaskar and others 2012; Angermann and others 2012b). At larger stream reach scale, substantial progress in identifying spatial patterns of groundwater–surface water exchange fluxes has been made by applying distributed fiber-optic temperature sensors (FO-DTS, Selker and others 2006). FO-DTS analyzes the backscatter of a laser signal propagating through a fiber-optic cable of up to several kilometers length, in this case ideally installed in the streambed or at the streambed—surface water interface. Using the seasonally variable difference in groundwater and surface water temperatures, groundwater upwelling is identified by hot or cold anomalies as long as end-member temperatures are significantly different. Although FO-DTS has been applied successfully for delineating groundwater–surface water exchange fluxes at multiple scales, its capabilities to provide quantitative predictions of fluxes remain limited and accuracies strongly depend on the existence of a suitable signal strength and signal size as well as correct calibration procedures (Rose and others 2013).

Isotopic Methods to Trace Sources of N and Denitrification

One of the major challenges of current research on the functioning of the continental environment is to develop integrative approaches allowing scale changes. Isotopic biogeochemistry is an important integrating tool. The basic idea is that the isotopic composition of a chemical species (δ15N for nitrogen) at a definite location reflects (i) its various sources and (ii) processes that affect its concentration (Figure 6). Despite some interpretation difficulties, especially in case of sewage treatment plants or manure application (Bedard-Haughn and others 2003), this method allows detection of the importance of denitrification activity in landscapes. δ15N signatures of plants and organisms are increasingly being used to identify the sources of N in aquatic environments, and to identify sites where extensive N cycling or transformations are occurring (Udy and Bunn 2001). It assumes that the δ15N of the riparian or aquatic vegetation will reflect its source, that is, the in-stream dissolved organic N or the groundwater. This may have been fractionated due to biogeochemical processes in transit. The N may also undergo small fractionation during assimilation within the plant (Udy and Bunn 2001) although Mariotti and others (1982) and Fry (1991) suggest that nitrate uptake by terrestrial vegetation appears to fractionate minimally or not at all (Clément and others 2003). Moreover, it was shown that diffusion and advection (physical processes) do not affect the isotopic composition of nitrate (Mariotti and others 1982; Semaoune and others 2012). Thus, any change in δ15N along flow paths can result either from a mixing of two different water bodies or from nitrification and/or denitrification. Moreover, a concomitant increase of δ15N and δ18O in the nitrate confirms the role of denitrification compared to mixing origin (Kendall1998; Sébilo and others 2006). When process rates are limited by nitrate diffusion through the water-sediment interface, the isotopic fractionation associated with denitrification is very poor. For this reason, only riparian denitrification (vs. benthic or in stream denitrification) generates a significant isotopic anomaly which leads an increase of the heavier fraction in the isotopic composition of nitrate (15N and 18O) with the decrease of nitrate concentration (Sébilo and others 2003).

Figure 6
figure 6

Geochemical proxies of the denitrification proximal drivers.

Full size image

Identify Anoxic Sites by Carbon Isotope Measurements of Dissolved Organic Carbon

The 13C/12C ratio of DOC could be used as a guide to quantify the proportion of DOM (and of waters) coming from wetlands, potential denitrification hotspots (Figure 6). Indeed, distinctively different carbon isotope signature can be expected for the DOM coming from uplands and adjacent wetlands because in the oxidative environment of the uplands, oxidative processes dominate during the decomposition of plant materials. Due to isotopic fractionation during those processes, residues are increasingly enriched in the heavier carbon isotopes (13C) as the lighter 12C will be preferentially involved in chemical reactions (for example, Wynn and others 2006). In contrast, wetland soils are characterized by anoxic conditions. The lack of oxygen results in an incomplete decomposition of organic material by anaerobic bacteria. Carbon compounds are preserved to a higher degree and keep their original (plant) isotopic signature. Therefore, the δ13C of soil organic matter in wetland soils is anticipated to be lighter than those of upland soils (Schaub and Alewell 2009). Thus, the δ13C of DOM in stream and river water provides a potentially extremely powerful tool to quantify the overall degree of interaction of drainage waters with wetland domains.

Another important distinction with the landscape can be between extended wetlands with a larger proportion of direct precipitation inputs compared to groundwater interactions, and fens where there is a greater proportion of groundwater. The differences between precipitation and groundwater chemistry have created different biogeochemical environments, and different chemical inputs to the suboxic to anoxic organogenic matrices bordering the surface water network. These differences lead to differences in the absorbance properties of the runoff water which can be used to determine the ratio of extended wetlands to forested fens (Berggren and others 2007; Laudon and others 2011). Given the complex nature of organic carbon, there are other potential tracers of C source environment that can be marshaled to the purpose of delineating dynamic contributions from different catchment environments. The basic starting point for this work with carbon character tracers is to measure the soil organic matter and DOM along upslope-riparian zone transects. These measurements, together with analysis of the DOM isotopic composition along water flow paths, can be used to quantify interaction of surface and subsurface water flows with denitrification hotspots. Although in-stream processing of carbon can eventually alter these carbon tracers, the focus of the new framework on headwater networks means that the time available for such processes is measured in hours to a day or two, limiting the scope for these processes to obscure the signal from these tracers of catchment DOM origin.

Rare Earth Elements

Marked negative Cerium (Ce) anomalies in groundwater are due to oxidation of Ce3+ to Ce4+ and subsequent secondary precipitation of cerianite (Figure 6). However, the Ce behavior in organic-rich waters is not completely controlled by redox processes. Indeed, the surface properties of dissolved organic matter are able to complex the Lanthanide series (rare earth elements: REE), including Ce, which has been shown to inhibit the development of negative Ce anomaly in waters (Dia and others 2000; Gruau and others 2004; Davranche and others 2005; Pourret and others 2010). Therefore, the extent of negative Ce anomaly development in ground and stream water can be used to quantify the degree of interaction of these waters with organic-rich zones in the basins, that is, with those areas where the electron acceptors necessary for the denitrification process to develop predominantly occur. More specifically, stream water showing no or insignificant negative Ce anomalies will be indicative of drainage basins where the ratio between potential denitrification landscape units (that is, wetland zones, thick organic-rich soil horizons,…) and overall water flux is high (high denitrification potential), while stream showing the reverse situation (that is, deep negative Ce anomalies) will be indicative of the opposite, that is, of basins where the ratio between potential denitrification units and effective water flux is low.

Conclusion

Small catchments constitute more than 80% of the drainage area of large river basins; thus they provide the right spatial scale for effective intervention to achieve water quality goals. Yet, an appropriate framework and methods to quantify the relationship between landscape structure/use and nitrogen fluxes/retention/removal is still lacking. We use an analysis of the existing frameworks related to denitrification as the basis for proposing a step forward in coupling landscape structures to denitrification at the small drainage basin scale. We propose to combine the landscape structure and the dynamic patterns of flow through that landscape pattern arrangement which produce the exposure times in anoxic zones and chemical inputs to denitrifying environments. In this context, exposure time distribution and Damköhler ratios provide an efficient means to evaluate and compare the denitrification capacity of different landscape units. Systematically combining local, process-based, and catchment-scale, integral assessments of denitrification capacity may ultimately lead to new modeling concepts to quantify catchment-scale denitrification. Integration of existing frameworks with new tools and methods offers the potential for significant breakthroughs in the quantification and modeling of denitrification in small drainage basins. This can provide a basis for improved protection and restoration of surface water and groundwater quality. Focusing on the hydrogeochemical architecture of small drainage basins can help place them in the center of monitoring and management issues related to protection and restoration of water quality.