Predicting Solute Transport Through Green Stormwater Infrastructure With Unsteady Transit Time Distribution Theory

In this study, we explore the use of unsteady transit time distribution (TTD) theory to model solute transport in biofilters, a popular form of nature‐based or “green” storm water infrastructure (GSI). TTD theory has the potential to address many unresolved challenges associated with predicting pollutant fate and transport through these systems, including unsteadiness in the water balance (time‐varying inflows, outflows, and storage), unsteadiness in pollutant loading, time‐dependent reactions, and scale‐up to GSI networks and urban catchments. From a solution to the unsteady age conservation equation under uniform sampling, we derive an explicit expression for solute breakthrough during and after one or more storm events. The solution is calibrated and validated with breakthrough data from 17 simulated storms at a field‐scale biofilter test facility in Southern California, using bromide as a conservative tracer. TTD theory closely reproduces bromide breakthrough concentrations, provided that lateral exchange with the surrounding soil is accounted for. At any given time, according to theory, more than half of the water in storage is from the most recent storm, while the rest is a mixture of penultimate and earlier storms. Thus, key management endpoints, such as the pollutant treatment credit attributable to GSI, are likely to depend on the evolving age distribution of water stored and released by these systems.

GSI called biofilters, also known as bioretention systems or rain gardens. As illustrated in Figure 1a, these vertically oriented systems filter water through planted soil or sand-based media and are easily integrated into the urban landscape over a range of scales (Grant et al., 2013;Roy-Poirier et al., 2010;Wong, 2006). Their possible features include: (1) a ponding zone that retains water prior to infiltration; (2) biological components including upright vegetation and naturally colonizing soil invertebrates and microorganisms; (3) engineered filter media (sand, sandy loam, or loamy sand with or without media amendments [e.g., biochar; Boehm et al., 2020;Mohanty & Boehm, 2014]); (4) a coarse sand transition layer; (5) a drainage layer consisting of coarse sand or fine gravel which can be lined or unlined and with or without an underdrain; (6) an overflow structure that releases excess storm water; and (7) a raised outlet to facilitate the formation of a permanently wet "submerged zone" (Clar et al., 2004;A. P. Davis et al, 2009;H. Kim et al., 2003;Payne et al., 2015;Rippy, 2015).
The water quality benefits attributable to GSI are often quantified based on the fraction of storm water pollutants (measured on a concentration or mass basis) removed during laboratory or field challenge experiments (Bedan & Clausen, 2009; A. P. Davis et al., 2009;Feng et al., 2012;Hatt et al., 2009;Kranner et al. 2019;L. Li & Davis, 2014;Y. Li et al., 2012;Parker et al., 2017;Ulrich et al., 2017). Much of this research has focused on the link between system design and pollutant removal, for example how the choice of plant species and the presence or absence of a submerged zone influences the removal of nutrients (e.g., H. Kim  transport processes that may influence solute transport through the lined biofilter used in our field experiments, including lateral infiltration and exfiltration with the surrounding soil (i-iii); outflow by evapotranspiration, gravitational drainage through the filter media, and shortcircuiting (iv-vi); diffusive exchange of water and solutes between pore spaces and micro-porosity within individual grains or organic material (vii), and uptake of water and solutes by plant roots (viii). (b) Photographs of the experimental set-up at the OCPW Low Impact Development Demonstration Facility (left), and the programmable controller and valve (top and bottom right). (c) Detail of the storm simulation experiments, including the vertical and horizontal positions of inflow and outflow tanks, control valve, biofilter test cell, and sampling set-up in the outflow tank (not to scale). OCPW, Orange County Public Works.

(a) (b) (c)
and how media amendments influence the removal of microbial contaminants and heavy metals (e.g., Y. Li et al., 2016;Mohanty & Boehm, 2014;L. Zhang et al., 2010). The effects of transient unsaturated flow, a defining feature of biofilters and GSI generally, are less often considered. Occasionally, transient unsaturated flow is indirectly acknowledged through experimental designs that incorporate an antecedent dry period between storm water dosing events (Chandrasena, Pham, et al., 2014;Payne et al., 2014). Similarly, biofilter design guidelines often recommend the inclusion of a submerged zone so that a portion of storm water passing through the biofilter spends a longer time undergoing treatment (e.g., nitrogen removal by denitrification) between storms (LeFevre et al., 2015;Payne et al., 2015). Yet, a detailed understanding of how transient unsaturated flow influences contaminant removal remains elusive.
Part of the problem is that transient unsaturated flow imposes severe challenges for predictive modeling. The Richards equation, which describes transient unsaturated flow through porous media, can be solved to estimate time varying flow and saturation through biofilters in one-, two-, or three-dimensions (e.g., using the numerical package Hydrus [Simunek et al., 2008]). These solutions can be coupled to the advection-dispersion equation (ADE) and one or more hypothesized pollutant removal mechanisms, to estimate pollutant removal in biofilters and analogous vadose zone systems (Massoudieh et al., 2017;Radcliffe & Simunek, 2010;Simunek et al., 2008). While this approach may work well as a theoretical exercise (Alikhani, et al., 2020;Dos Santos et al., 2013) and under highly controlled conditions in the laboratory (Al-Mashaqbeh & McLaughlan, 2012;Behroozi et al., 2020;Henrichs et al., 2009;Horel et al., 2015;Saiers & Lenhart, 2003;Trenouth & Gharabaghi, 2015) or field (Boivin et al., 2006;Massoudieh et al., 2017), its general application is limited by the information required (boundary conditions, soil hydraulic properties, root profiles, and etc.) and the fact that biofilters, like the catchments they are nested within, are "complex, heterogeneous, and poorly characterized by direct measurement" (Kirchner, 2009).
Alternatively, mass balance over a control volume drawn around the biofilter media-so-called "bucket models"-can track the temporal evolution of soil moisture and solutes in these systems. Daly et al. (2012) used a bucket model to derive the probability density of water volume stored in a biofilter based on a stochastic description of rainfall together with biophysical models for gravitational drainage and evapotranspiration (ET). The power of bucket models lies in their simplicity; an enduring challenge has been how to leverage their output into estimates of pollutant removal. Daly et al. (2012) bridged this gap by relating the pore fluid total nitrogen concentration in a biofilter to degree of soil saturation (estimated from the bucket model) on the premise that low saturation levels are associated with "a reduction in nitrogen plant uptake and denitrification with a consequent accumulation of nitrogen in the filter media that is then washed out during the next inflow event." While clever, this approach does not address the more general problem of predicting pollutant removal as storm water passes through a biofilter or other GSI. Further, Daly et al.'s bucket model was implemented at a daily time step, and consequently within-storm processes (such as pollutant breakthrough curves) cannot be resolved. Randelovic et al. (2016) proposed a hybrid approach in which the unsteady water balance is solved with a bucket model while pollutant removal in the filter media and submerged zone is predicted with the one-dimensional ADE coupled with one or more pollutant removal mechanisms. Their framework accurately predicts micropollutant removal in field-scale storm water biofilters (Randelovic et al., 2016) but at the cost of significant complexity-the model consists of 25 coupled equations and requires the specification of 32 variables. Refinements of this model continue to be published (Shen et al., 2018;K. Zhang et al., 2019) and incorporated into practice oriented GSI design software (e.g., MUSIC [eWater Ltd., 2020]) (reviewed in Jefferson et al., 2017 andC. Li et al., 2017).
In this paper we propose and test an entirely new approach for predicting unsteady solute transport through GSI: unsteady transit time distribution (TTD) theory. TTD theory combines the simplicity of bucket models with the temporal resolution and physical insights provided by process-based models of solute fate and transport in transiently unsaturated flow systems. The theory was developed by hydrologists to characterize the myriad transport pathways and timescales associated with water and solute transport through catchments to streams  but has since been applied to a diverse array of environmental problems (e.g., Metzler et al., 2018;Smith et al., 2018). In place of a mass conservation equation (such as the ADE), TTD theory is premised on an unsteady conservation equation for the age distribution of water entering, stored in, and leaving a system. These age distributions, in turn, encode all information needed to estimate solute breakthrough concentrations, including the time history of inflows (e.g., in the case of a biofilter, the magnitude and timing of storm events) and the transformation reactions that occur as water flows through the system along diverse flow paths. In short, TTD theory elegantly links hydrologic processes to water quality outcomes (Hrachowitz et al., 2016).
There are at least two reasons why TTD theory is a potentially important advance over current approaches for modeling pollutant removal in GSI: (1) parsimony and (2) extensibility. The unsteady TTD model presented later requires the specification of just one unknown parameter, the effective size of the biofilter (taking lateral exchange of water and solutes with the surrounding soil into account). TTD models can be linked in series and parallel (Bertuzzo et al., 2013;Hrachowitz et al., 2016) to represent GSI networks and hydrological response units (e.g., hillslopes, groundwater, wetlands, and rivers) in urban catchments. Thus, TTD theory directly addresses a significant limitation with existing GSI modeling frameworks; namely, their "uncertainty in simulating the propagation of flows through pathways such as storm water networks, pervious runoff, and subsurface flows" (C. Li et al., 2017).
The paper is organized as follows. We begin by developing the TTD modeling framework needed for GSI applications (Section 2). A field-scale test of the TTD theory is then described (Section 3) followed by experimental and modeling results (Section 4), a discussion of water quality implications (Section 5), and conclusions (Section 6). All variables are defined in Table 1.

Modeling Framework
The application of TTD theory to GSI entails three steps: (1) a control volume is drawn around the feature of interest, in our case the soil media component of a biofilter; (2) an unsteady water balance is performed over the control volume, taking into account time-varying inflows, outflows, and change in water storage; and (3) the age distribution of water in the control volume, and in water leaving the control volume by gravitational drainage and ET, is estimated from TTD theory, along with any water quality metrics (e.g., solute breakthrough curves) of interest. In this section we begin by describing two different approaches for preparing the unsteady water balance, including a simple bucket model (Section 2.1) and a numerical solution of the Richards equation coupled with soil hydraulic functions (Section 2.2). We then solve the age conservation equation (Section 2.3.3) and derive a set of expressions for the age distribution's central tendency and spread (Section 2.3.4) and the breakthrough concentration of a conservative solute (Section 2.3.5).

Bucket Model Water Balance
Equation 1a is an unsteady water balance over a control volume drawn around the biofilter media, where the variable t [T] is time and the functions S(t) [L], J(t) [L T −1 ], Q(t) [L T −1 ], and ET(t) [L T −1 ] represent, respectively, the time-varying volume of water in storage, infiltration rate of water into the biofilter from the ponding zone, gravitational discharge of water out of the biofilter, and ET across the top boundary of the biofilter (Daly et al., 2012). All volumes and flows are normalized by the biofilter's surface area.
The initial condition (Equation 1b) stipulates that the area-normalized water in storage at time t = 0 is S=S 0 [L]. To solve Equation 1a we must first specify the storage dependence of all terms on the right-hand side. These are discussed in turn.

S T (T, t)
Age-ranked storage, equal to area-normalized volume of water stored in the biofilter media at any time t with ages T or younger L P RTD (T, t) CDF for the fraction of stored water with ages less than or equal to T at time t (the unsteady RTD of water in storage) - This simple expression approximates the three phases of infiltration (Williams et al., 1998) as follows. Infiltration equals inflow during the Filling Phase, which begins when storm water first enters the ponding zone and infiltration is dominated by capillary forces: S(t) < S max , J(t) = I(t). Infiltration equals the saturated hydraulic conductivity during the Transition Phase as the biofilter approaches full saturation: S(t) = S max , J(t) = K sat . During this phase, water level in the ponding zone rises whenever inflow exceeds the media's saturated hydraulic conductivity. Infiltration is zero during the Draining Phase, which commences once inflow has ceased and the ponding zone has drained: S(t) < S max , J(t) = I(t) = 0. While process-based models of infiltration are available (e.g., Green & Ampt, 1911), Equation 2 is consistent with the field observations described later (see Section 4) and its sole variables (K sat and S max ) are easily measured biofilter design parameters (Le Coustumer et al., 2009Coustumer et al., , 2012Payne et al., 2015;Peng et al., 2016).

Dependence of Gravitational Discharge on Storage
Kirchner (2009) posited that streamflow out of a catchment can be represented by a single nonlinear function of the catchment's water storage, Q(t) = f(S(t)). One such functional relationship derives from the power-law recessional model for streamflow where the prefactor, a, and exponent, b, are empirical constants: When coupled with an unsteady water balance over the catchment, Kirchner demonstrated that Equation 3a can be manipulated to yield an algebraic expression for streamflow as a function of storage (Equation 14 in Kirchner, 2009   The constants appearing in Equations 3a and 3b are related as follows: The new variables, S min and g, are emergent properties of the transient unsaturated flow field; that is, they must be determined empirically based on experimental observations or numerical solutions of the Richards equation.

Dependence of ET on Storage
ET also depends nonlinearly on water storage, but only when storage falls below a critical value known as the incipient water stress (Allen et al., 1998;Daly et al., 2012). Above the incipient water stress, ET approaches a maximum rate (set by local environmental conditions, including wind speed, vapor pressure deficit, temperature, and plant-specific characteristics) called potential ET. While biofilters often operate at or below the incipient water stress (Hess et al., 2019) this was not the case for the experiments described later, which involved a sequence of back-to-back simulated storms. Accordingly, in this study we approximated ET with an hourly time series of reference crop potential ET (cPET) following FAO guidelines (Allen et al., 1998) and based on measurements at, or nearby, the field site together with plant-specific traits (details in Text S1 of the Supporting Information).

Numerical Implementation
The water balance bucket model was solved by substituting into Equation 1a the above expressions for infiltration (Equation 2) and gravitational discharge (Equation 3b), along with hourly estimates of cPET (Section 2.1.3). The model was then forced with timeseries (sampling frequency ∼1 min −1 ) of measured storm water inflow and numerically integrated (details in Text S2 and Figure S1). These simulations yielded ∼1 min −1 timeseries of infiltration, storage and gravitational discharge over the 17 simulated storm events described in Section 3.

Numerical Solution of the Richards Equation
To calibrate the gravitational discharge expression (Equation 3b) and as a check on the bucket model predictions described above, ∼1 min −1 time series of infiltration, storage and gravitational discharge were also simulated with the one-dimensional Richards equation ( ). The exception was saturated hydraulic conductivity, K sat , which was estimated from measurements of peak discharge and in situ measurements with a modified Philip-Dunne Infiltrometer (Text S3).

Solving the Age Conservation Equation
The age distribution of water in the control volume surrounding the biofilter's soil media is governed by the following age conservation equation (Botter et al., 2011;Harman, 2015): The conservation equation's dependent variable, age-ranked storage S T (T, t) [L], represents the area-normalized volume of water stored in the biofilter media control volume at any time t with ages T or younger. Age-ranked storage is defined mathematically as the product of the area-normalized volume of stored water, S(t), and the cumulative distribution function (CDF) for the fraction of that stored water with ages less than or equal to T; that is, the stored water's residence time distribution (RTD), P RTD (T, t) (Equation 4b).
As the age variable, T, becomes large, the RTD's CDF tends to unity and the age-ranked storage function collapses to the area-normalized volume of water in storage: ) ensures that no water stored in the control volume has an age less than T = 0. The initial condition (Equation 4d) implies that, at time t = 0, the volume of "original" water in storage, S 0 , has a single age, ). As applied to biofilters, Equation 4a equates the change of age-ranked storage of water in the biofilter media (left-hand side) to the infiltration of storm water of age T = 0 (first term on right-hand side); outflow of water by gravitational discharge (second term) and ET (third term) with age distributions   the biofilter as gravitational discharge and ET with ages T or less at time t. The backward arrows on these CDFs indicate they are "backward TTDs"; that is, they represent the age distribution of water leaving the biofilter at time t. A corresponding set of "forward TTDs" can be written for the life expectancy of water parcels entering the biofilter at time, t i . The relationship between forward and backward TTDs is given by Niemi's Theorem (Benettin, Rinaldo, & Botter, 2015;Harman, 2015;Niemi, 1977). Under unsteady conditions and depending on the nature of the storage selection function (see next section), the backward TTDs for gravitational discharge and ET are not necessarily equal, nor are they necessarily equal to the RTD of water in storage (Botter et al., 2011). , which maps the fraction of outflow with ages less than or equal to T (i.e., the CDF form of the backward TTD for discharge or ET) to the fraction of age-ranked water in storage with that age or younger "selected" for outflow by either gravitational drainage or ET (Botter et al., 2011;Harman, 2015):

Ranked StorAge Selection Function
In principle, the functional form of the SAS function can be calculated by averaging the ADE for solute transport over the control volume (Benettin et al., 2013;Rinaldo et al., 2015). For the purposes of this study, however, we adopted a "uniform SAS" function, under the assumption that water in storage has an equal probability of being selected for outflow regardless of its age (Harman, 2015): Uniform SAS functions often apply to systems, such as ours, that are far from well-mixed (Benettin, Bailey, et al., 2015;Benettin et al., 2013;Bertuzzo et al., 2013;Danesh-Yazdi et al., 2018;M. Kim et al., 2016;Rodriguez et al., 2018).

Exact Solution for Age-Ranked Storage under Uniform Selection
Under uniform sampling, the age conservation equation (Equation 4a) can be solved exactly for certain choices of initial and boundary conditions (Bertuzzo et al., 2013;Botter et al., 2011). Equation 6a is one such solution that satisfies the initial and boundary conditions presented earlier (Equations 4c and 4d, see Text S4 for derivation); the superscript "U" denotes that the solution is premised on the choice of a uniform storage selection function.
According to Equation 6a age-ranked storage (left-hand side) is influenced by the evolving age distribution of both "original" water in storage at time t = 0 (first term on right-hand side) and "young water" that infiltrates during storm events (second term). This solution was numerically integrated (details in Text S4) to yield ∼1 min −1 timeseries of age-ranked storage in the biofilter, after substituting bucket model simulations for infiltration, J(t), storage, S(t), and gravitational discharge, Q(t) (Section 2.1).

Age Structure of Stored Water in the Biofilter
Under uniform selection the backward TTDs for gravitational discharge and ET are equal, and equal to the RTD of water in storage (compare with Equation 4b) (Harman, 2015): The fifth, 50th, and 95th percentile ages of water in storage and outflow at any time, t, can be obtained from Equation 7a by numerically solving the following implicit equations for water age: S T (T 0.05 ,t)/S(t) = 0.05, S T (T 0.5 ,t)/S(t) = 0.5, and S T (T 0.95 ,t)/S(t) = 0.95. The age-ranked storage's probability density function (PDF) can be calculated from Equation 7a by differentiation where the symbol  denotes the Dirac delta function:

Water Resources Research
The mean age in storage and outflow immediately follows by taking the first moment of the PDF for ageranked storage: Further details on the derivation and numerical implementation of Equations 7b and 7c are described in Text S5.

A TTD Theory for Solute Fate and Transport through a Biofilter
The concentration of a nonreactive and nonadsorbing (i.e., conservative) solute in water leaving the biofilter by gravitational discharge, C Q (t), can be calculated by convolving the PDF of the backward TTD (Equation 7b) with the concentration of solute, C J (t i ,T), that entered the biofilter at time, t i = t − T, and exited the biofilter as gravitational discharge at time t and age T (Harman, 2015): Despite its simplicity, this convolution integral incorporates a rich set of processes, including unsteadiness in the biofilter's water balance (e.g., time-varying inflows, outflows, and storage, encoded through the time-evolution of the backward TTD) and unsteadiness in the solute concentration entering the biofilter from the ponding zone (through the dependence of C J (t i ,T) on the inflow time, t = t i ) (Harman, 2015). Combining Equations 7b and 8a, we arrive at the following solution for the concentration of a conservative solute in water discharged from the biofilter, where C 0 is the concentration of solute present in the original water stored in the biofilter at time, t = 0: For the experiments described later, a subset of 17 simulated storms were tagged with bromide, which we assumed behaved conservatively (Levy & Chambers, 1987). The inflow concentration for these storms can be expressed as follows, where C J,m , t m,s , and t m,e are the m-th storm's bromide concentration, start time and end time, respectively, and the sum is taken over all N storms: Substituting Equation 9a into Equation 8b, setting C 0 = 0 (because, in our experiments, no bromide was present in the biofilter's original water), and using the distributive property of integration, we arrive at the following expression for bromide concentration in water leaving the biofilter by gravitational discharge (details in Text S6): In deriving Equation 9b we have assumed that plants in the biofilter take up bromide and water in the same proportion, which may not be the case in practice (e.g., if a solute is excluded from plant uptake its pore fluid concentration will increase over time by in situ evaporative concentration [Bertuzzo et al., 2013;Harman, 2015]). However, in the field experiments described later (Section 4) ET represents a very small portion of the overall water balance and hence in situ evaporative concentration can be neglected in our case. Time series (∼1 min −1 ) of bromide breakthrough concentration were simulated with Equation 9b following the numerical procedures described in Text S6.

Orange County Public Works Biofilter Test Facility
Field-scale biofilter challenge experiments were carried out at the Orange County Public Works (OCPW) low impact development demonstration facility located in the City of Orange, Orange County, California. Experiments were conducted in a biofilter test cell (approximately 2.4 × 1.5 × 0.6 m deep) built by a local contractor with previous GSI construction experience (Tobo Construction, Figure 1b). The test cell, which was lined and outfitted with an underdrain, consisted of a concrete slab floor and four cinderblock walls extending approximately 0.5 m above the filter media surface to create a ponding zone (Figures 1b and 1c).
The filter media consisted of sand (65%), sandy loam (20%), and compost (15%) (v/v basis). In January 2017, the media was planted with a European gray sedge, Carex divulsa tumulicola. When we conducted our first set of experiments in June 2018 the plant community also included opportunist ruderal weed species (e.g., the common dandelion). While our biofilter is open to the atmosphere, no measurable rainfall occurred during the set of experiments described below.
The research team retrofitted the biofilter with an upstream 1,890 L "inflow" tank (Custom Roto-Molding, Inc., Caldwell, ID) which drained by gravity through a programmable control valve (Sigma Controls, Inc., Perkasie, PA) to the biofilter's ponding zone (Figure 1b). The weight of the inflow tank was monitored continuously at ∼10 Hz (WinWedge, TAL Technologies Inc., Philadelphia, PA) with a calibrated industrial scale (PCE-SW 3000N Pallet Scale, PCE Americas Inc., Jupiter, FL). The weight measurements were lowpass filtered (J. C. Davis, 2002), differentiated, and divided by the density of water to yield ∼1 min −1 estimates for the volumetric discharge of water entering the ponding zone (these measurements correspond to the inflow term, I(t), in Equation 2). Following the experiments conducted in the summer of 2018 we discovered that, during construction, a ca., 5 cm diameter hole had been drilled through the base of the cinderblock wall separating our test cell from the adjacent test cell (and through the wall separating the adjacent test cell from the next test cell and so on) to accommodate a buried irrigation pipe. A substantial fraction of storm water added to our test cell (approximately 50%) laterally exfiltrated to the adjacent cell through this hole. While not part of our original design, this feature made for a more realistic field experiment, as most operational biofilters undergo at least some degree of subsurface exfiltration (e.g., Brown & Hunt, 2011). Indeed, our exfiltration rate of ∼50% is close to the storm water volume reduction design goal for GSI of 67% (A. P. Davis, 2008). To account for lateral exfiltration, gravitational discharge estimated from the bucket model (Section 2.1) was routed as follows. A fixed fraction, α, was assigned to the underdrain and the rest, 1−α, to lateral exfiltration to the adjacent test cell (Figure 1c). The fraction α was estimated using several independent experimental methods (see Text S7 and Figure S6).
Water exiting the biofilter through the underdrain flowed by gravity through a buried manifold to an underground sump and from there was periodically pumped (Model 98 Sump Pump, Zoeller Pump Company, Louisville, KY) up to an "outflow" tank sitting on a calibrated industrial scale at ground level (identical to the inflow tank set up, Figure 1c). A timeseries (∼1 min −1 ) of volumetric discharge entering the outflow tank was estimated from high frequency (10 Hz) measurements of the tank's weight following the same procedure described above for the inflow tank (corresponding to the Q(t) term in Equations 1a and 3b). Time series of volumetric discharge and bromide concentration measured at the outflow tank were time-shifted backwards by 30 min to account for the transit of water from the biofilter's underdrain to the outflow tank and the overly fast response of the bucket model to storm events (Text S8).

Experimental Storm Hydrograph
Municipal separate storm sewer system (MS4) permit requirements for the Santa Ana Region (where our experiments were carried out) stipulate that new development and significant re-development projects include storm water control measures sufficient to capture runoff volume generated by the 85th percentile storm, which at this field site corresponds to 2.1 cm of storm water depth (volume per unit catchment area) over a 24-h period (OCPW, 2017). With this regulatory requirement in mind, we designed our experimental storm hydrograph as follows: (1) seven 24-h rainfall events were selected from measurements at an onsite rain gauge over the years 2011-2016; (2) an average hyetograph was constructed from these seven events after aligning peaks and standardizing the total 24-h rainfall depth to 2.1 cm; and (3) a design storm hydrograph was calculated from the average hyetograph by the Rational Method (Brooks et al., 2013) assuming a unit runoff coefficient (corresponding to 100% imperviousness) and a catchment area of 82.3 m 2 . The corresponding biofilter-to-catchment area ratio (4.5%) is typical for urban landscapes in Southern California (Ambrose & Winfrey, 2015). The design storm hydrograph was subsequently programmed into the control valve to regulate the flow of water into the biofilter's ponding zone during each simulated storm (Figure 1b) (see Text S9 and Figure S8 for a comparison of measured and design storm hydrographs).

Bromide Tracer Experiments
A sequence of simulated experimental storms (each of which conformed to the storm hydrograph described in Section 3.2) were discharged to our experimental biofilter over a 5-day period in the summer of 2018 (June 25-29) and again over a 5-day period in the summer of 2019 (June 1-5). Ten storms were simulated in 2018, one in the morning and another in the afternoon on each day. The afternoon storms were spiked with bromide (final concentration ∼50 Brmg/L) while the morning storms were bromide free. The storm sequence in 2019 consisted of: (1) two bromide-free storms on the first day, one in the morning and one in the afternoon; (2) 2 days later a single bromide-spiked storm in the morning (final concentration of 124 Brmg L −1 ); and (3) over the following 2 days two bromide-free storms per day, one in the morning and one in the afternoon. Three replicate 40 mL samples were collected from the inflow tank before each simulated storm. Outflow samples were collected as follows. Water entering the outflow tank from the sump was directed into a continuously overflowing 5 L bucket affixed to the top of the tank (the bucket overflowed into the tank, see Figure 1c). Water in the bucket was continuously subsampled (

Lateral Exfiltration and the Effective Volume of the Biofilter
Each experimental storm discharged roughly the same volume of water (1,400-1,500 L) to the biofilter's ponding zone over 1 to 2 h. The volume of water captured in the outflow tank varied by storm, from 378 to 751 L (25%-49% of the inflow volume) for the 10 experiments conducted in 2018, and from 266 to 654 L (21%-46% of the inflow volume) for the seven experiments conducted in 2019 (Table S1). Across all 17 storms, ET was a minor component of the water budget (<0.3% of the ∼1400 L added per storm, Table S1). Thus, the difference between these inflow and outflow volumes either went to increasing storage or lateral exfiltration to the adjacent test cell (see Text S7).
The fraction of inflow volume recovered at the outflow tank is inversely correlated with antecedent dry period (R 2 = 0.82, Figure S6), consistent with the hypothesis that at least some of the unrecovered water goes to storage. Extrapolating the fractional water recovery back to an antecedent dry period of zero hours (under the premise that the change in storage should be zero in this case) we estimate that, in both 2018 and 2019, approximately α = 46% of the water added to the biofilter is routed to the outflow tank while 1 − α = 54% is lost to lateral exfiltration (Text S7). These results are consistent with loss rates we independently measured under steady-state flow conditions (Text S7) and the observed wetting of biofilter media in the adjacent test cell (data not shown), along with previously published modeling studies (Browne et al., 2008;J. G. Lee et al., 2015) and field measurements (Winston et al., 2016) that indicate exfiltration is a dominant mechanism for water volume reduction in GSI. At our site, some of the exfiltrated water and solute may eventually find its way back to the outflow tank, for example by recirculating back through our biofilter test cell (mechanism [iii], Figure 1a) or transiting along another subsurface route to the buried collection manifold (indeed, the test cell adjacent to our biofilter also had an underdrain that could have contributed flow and solute to the manifold and, ultimately, to the outflow tank). Thus, exfiltration can potentially increase the effective volume of our biofilter during solute transport. We return to this idea in Section 4.4.

Power-Law Model for Gravitational Discharge
Kirchner's power-law model for gravitational discharge (Equation 3b) could not be evaluated from our inflow and outflow measurements, because lateral exfiltration precluded accurate estimates of the gravitational discharge and water storage terms of the equation. Instead, hourly time series for these two quantities were numerically simulated, with Hydrus 1D, over the 17 experimental storms in 2018 and 2019 (Section 2.2). Consistent with Kirchner's power-law relationship, when normalized and plotted on a log-log basis, the Hydrus-generated time series of discharge and storage collapse to a single line for Q(t) > 0.05K sat (Figure 2). Inferred values of the power-law exponent and minimum storage are the same, within error, across both years (2018: g = 4.99 ± 0.01 and S min = −3.6 × 10 −20 ± 1.9 × 10 −4 m; 2019: g = 5.00 ± 0.01 and S min = −3.4 × 10 −21 ± 2.2 × 10 −4 m). This result indicates these parameters are robust to changes in the sequence and length of antecedent dry periods as well as to changes in saturated hydraulic conductivity (within each simulated storm sequence, the saturated hydraulic conductivity declined over time, see Text S3). Substituting g = 5 into Equation 3d yields a recessional exponent of b = 1.8, which is toward the flashy end of the allowable range, b ª ¬ º ¼ 1 2 , (Kirchner, 2009), consistent with the small storage volume of the biofilter. Our inferred value for the power-law exponent is also concordant with exponent values inferred by Bertuzzo et al. (2013) for gravitational drainage from the vadose zone of a 46 km 2 catchment in Switzerland (compare g = 5 with the posterior distribution for the exponent c in their Figure 4).
On the other hand, the inferred value for minimum storage, S min = 0, is at odds with the expectation that some residual water (e.g., equal to the area normalized volume of water in the biofilter media at field capacity [Hillel, 2003]) should be present in the biofilter following the cessation of gravitational discharge. Two points are relevant here. First, for the set of Hydrus 1D simulations used in the inference step, water storage in the biofilter never fell below the soil media's field capacity due to the large volume of water associated with each storm (compared to the biofilter's void volume) and the back-to-back nature of the storm sequences in 2018 and 2019 (Section 4.1). Thus, one plausible explanation for the inferred value, S min = 0, is that the biofilter was never dry enough to observe an influence of minimum storage on gravitational discharge. Second, the inferred value of the power-law exponent, g = 5, implies that discharge declines very rapidly with storage. For example, substituting g = 5 and S min = 0 into Equation 3b we predict that gravitational discharge will be ≤0.05% of the media's saturated hydraulic conductivity (Q ≤ 0.0005 * K sat ) when storage falls below the media's field capacity (S < S fc = 0.054m for our biofilter). Put another way, the large power-law exponent ensures that there will be very little gravitational discharge once storage drops below field capacity, even if the minimum storage is set to zero. For this reason and consistent with other field applications of Equation 3b (e.g., Bertuzzo et al., 2013), we adopted the inferred value, S min = 0, in the bucket model and TTD simulations presented below.   To compare simulated and measured discharge, the former was multiplied by α = 0.46 (to account for the fraction of water that is laterally exfiltrated) and the latter was backward time shifted by 30 min. cPET, crop potential evapotranspiration.

(a) (b)
biofilter media) is consistent with our field observations and the predicted gravitational discharge rates closely match measurements at the outflow tank (light blue curves, bottom panels of Figures 3a and 3b).

TTD Theory Predictions for Bromide Transport
In 2018, our semi-periodic study design consisted of a bromide-free "flushing" storm in the morning (orange arrows in Figure 4a) and a bromide-spiked "tracer" storm in the afternoon (black arrows in Figure 4a) of each day. By the second day of the storm sequence, the normalized bromide breakthrough curves (BTCs) settled into a periodic pattern, oscillating between  ,1 / Q J C C 0.3 and 0.6 during the morning and afternoon storms, respectively (black dots in lower graph, Figure 4a). Here, the variable C Q (t) represents the measured bromide concentration at the outlet tank and the variable C J,1 = 48.7 mg L −1 represents the initial bromide concentration measured in the first afternoon tracer storm (across all five tracer storms the initial bromide concentrations were C J,m = 48.7, 48.8, 49.3, 48.3, and 44.4 mg L −1 ). The bromide BTC predicted by TTD theory also follows a periodic pattern (solid curve, bottom graph, Figure 4a), but the model consistently under and overpredicts measured bromide concentrations in the morning and afternoon storms, respectively. These model predictions were calculated from Equation 9b after specifying the timing and initial bromide concentrations associated with each bromide-spiked storm, and running bucket model simulations for J(t), S(t), and Q(t) after setting the maximum storage term (see Equation 3b) equal to the void volume of the biofilter, S max = 0.246 m.
The TTD model's tendency to overshoot bromide measurements implies it is oversampling young water; that is, the predicted bromide BTC contains too much bromide-free water during the bromide-free morning storm, and too much bromide-spiked water during the bromide-spiked afternoon storm. One possible explanation for this result is that the uniform storage selection function, which underpins our model (see Equation 9b and discussion thereof), oversamples young water for gravitational discharge. Alternatively, the storage selection function is fine but there is not enough old water in storage to select from. While the former explanation cannot be ruled out (indeed, the science of selecting SAS functions is an active area of current research [Harman, 2019]), the latter explanation is compelling for several reasons. First, the void volume of our biofilter (∼900 L) is less than the volume of water flowing into the biofilter with each experi-PARKER ET AL.   mental storm (∼1,400 L). Therefore, as far as the model is concerned, our biofilter has very limited capacity to store older water from penultimate and older storms. Second, a substantial fraction (>50%) of the inflow volume leaves our biofilter by lateral exfiltration. As noted in Section 4.1, some of this exfiltrated water may eventually return to the outflow tank and thereby increase the effective volume that solutes experience as they transit through the system.
To test the last hypothesis-that the effective volume for solute transport is larger than the biofilter's physical volume-we split the measured bromide data from 2018 into a calibration period (the first bromide-spiked storm, storm #2) and a validation period (all other storms). We then inferred a value of the biofilter's void volume by minimizing the root-mean square error (RMSE) over the calibration period ( Figure S9). The optimal volume thus obtained (S max = 0.73 m, RMSE = 0.041) is about two times larger than the physical volume of the biofilter (S max = 0.246 m) consistent with the hypothesis that a substantial fraction of the water that exfiltrated into the adjacent biofilter test cell eventually returned to the outflow tank. When the inferred value of S max = 0.73 m is substituted back into the bucket model and the hydrologic water balance is recomputed, Equation 9b closely tracks the bromide BTC over the validation period (storms #3 through #10, bottom graph in Figure 4b).
Application of TTD theory to the 2019 storm sequence yields similar results. If the maximum volume of the biofilter is set equal to its physical volume (S max = 0.246 m), Equation 9b consistently underpredicts bromide breakthrough during the four bromide-free flushing storms (storms #4, 5, 6, and 7, bottom graph in Figure 4c). However, when the effective volume is raised to S max = 0.42 m (obtained by minimizing the RMSE for storms #4-7, see Figure S9; RMSE = 0.31) the model's performance improves markedly (Figure 4d). The void volume inferred from the 2018 experiments (S max = 0.73 m) is about 74% larger than the void volume inferred from the 2019 experiments (S max = 0.42 m), perhaps reflecting changes in the biofilter test cell over the 2 years; for example, after the 2018 experiments the hole at the base of the test cell was partially sealed and the media was replanted (see Text S9 and discussion in Section 5.4 below).

Age Distribution of Water Leaving the Biofilter by Gravitational Discharge
By selecting a uniform SAS function for our model (Section 2.3.2), the backward TTDs for gravitational discharge and ET are equal to the RTD of water stored in the biofilter. Thus, under uniform storage sampling, the age distribution of water in storage is equal to the age distribution of water leaving the biofilter as gravitational discharge.
During the 2018 experiments, predictions for the median age of water stored in the biofilter (Equation 7a and discussion thereof) follows a semi-periodic pattern, increasing linearly with time between storms (as water stored in the biofilter ages) and rapidly declining to near zero during storm events (as incoming storm water, of age T = 0 h, fills the biofilter, Figure 5a). The 5th and 50th (median) age percentiles overlap but the 95th age percentile is much older, indicating that the age distribution is positively skewed (Ang & Tang, 2007). The fifth and 50th percentiles overlap because, at any time t, more than 50% of water stored in the biofilter is from the most recent storm with an age roughly equal to the antecedent dry period. The 95th percentile age is much older because the rest of water in storage (i.e., water not from the last storm) is from penultimate and earlier storms.
For the simulations presented in Figure 5a we arbitrarily set the initial age of "original" water (i.e., water that was initially present in the biofilter at time, t = 0) at T 0 = 50 h. Until the fourth storm, this original water constituted more than 5% of water stored in the biofilter, as evidenced by the upward slope of the gray band in Figure 5a (the upward slope reflects the fact that that original water in storage is aging linearly with time). After the fourth storm, the original water's contribution to storage drops below 5%, as evidenced by a steep drop in the 95th percentile around t = 28 h (Figure 5a). Thus, four storms were required to flush out 95% of the original water, even though more than 50% of water in the biofilter, at any given time, is from the last storm. A similar pattern is evident for the set of experiments conducted in 2019 (Figure 5b). Across both years the mean age is 5-20 h older than the median age, consistent with a positively skewed age distribution (Ang & Tang, 2007).

Mapping out the Contribution of Past Storms to Present Storage
TTD theory also allows us to determine the relative contribution of all past storms to water stored in a biofilter at any time, t. If the n-th storm begins at time, t = t b,n , then the fraction, f n (t) [−], of water in storage with that age or younger can be estimated from the RTD's CDF (see Section 2.3.4) (Benettin et al., 2017;Kirchner, 2016;Lutz et al., 2018): We applied Equation 10 to all seven storms simulated in 2019, along with the original water present in the biofilter at time, t = 0 ( Figure 6). The upper bound of each color band in the figure represents the fraction of water in storage that is younger than the oldest water from the storm indicated. The lower bound of the same color band represents the fraction of water in storage that is younger than the oldest water from the next storm, and so on.
The influence of biofilter hydrology on the age structure of stored water (and by implication the age structure of water leaving the biofilter by gravitational drainage under uniform sampling) is striking. During the Filling Phase of each storm (e.g., Storm #3 in Figure 6) new water entering the biofilter from the ponding zone rapidly dominates the age distribution of water in storage for two reasons: (1) the new water fills up portions of storage that were previously dry; and (2) the new water displaces older water, driving it out of the biofilter as gravitational drainage. The first mechanism explains why different storms initially dominate storage to different degrees. For example, the assumed volume of original water in storage at time t = 0 was relatively small in these simulations (S 0 = 0.092 m, corresponding to the biofilter's field capacity) which explains why Storm #1 very quickly constituted more than 80% of the biofilter's storage (orange band in Figure 6). Because the initial moisture content for the first storm event was set equal to the soil media's field capacity, this last example may be more representative of the situation, during long antecedent dry periods, where biofilters are irrigated only to the extent necessary to maintain soil moisture at field capacity (Ambrose & Winfrey, 2015).
Under a uniform SAS, all water parcels (regardless of their age) have an equal probability of being selected for outflow by gravitational discharge or ET. This explains why, during the Draining Phase, the age structure PARKER ET AL.
10.1029/2020WR028579 Top panels indicate the infiltration (red curve) and gravitational discharge (black curve) rate for each storm sequence. For these simulations, we adopted inferred values of maximum storage: S max = 0.73 and 0.42 m for 2018 and 2019, respectively. The age of original water was arbitrarily set to T 0 = 50 h and its volume was taken as the product of initial moisture content (s 0 = 0.22, corresponding to field capacity) and maximum storage:

(a) (b)
of water in storage does not change (i.e., during this phase all boundaries in Figure 6 are horizontal lines). Figure 6 also vividly illustrates the two key attributes of biofilter storage discussed previously: at any time about >50% of water in storage is from the most recent storm while the rest is a mixture of penultimate and earlier storms.

Age Structure and Water Quality
The age structure of water in storage has significant implications for the treatment credit attributable to, and the pollution exported by, GSI. For example, during the 2019 storm sequence we included a "worst case" scenario (from a water quality perspective) by making one of the experimental storm events, Storm #3, a 50:50 mixture of storm water and raw sewage. From the thickness of the horizontal portions of the green band in Figure 6 (divided by two to account for the storm's 50:50 sewage:storm water mixture) we can infer that the fraction of water in storage that is raw sewage declined over time as follows: 35% (Storm #3), 11% (Storm #4), 4% (Storm #5), 1.2% (Storm #6), and 0.5% (Storm #7). Raw sewage harbors very high concentrations of human fecal bacteria (e.g., in the range of 10 6 E. coli mL −1 [Garcia-Aljaro et al., 2018]). Therefore, even after the biofilter has been flushed with four sewage-free storms, the E. coli concentration in gravitational drainage could still be as high as 5,000 mL −1 -more than enough bacteria to close beaches if the biofilter drained to a recreational lake or river (US EPA, 2018). This example assumes that bacteria behave conservatively which is rarely the case ; C. M. Lee et al., 2006). Indeed, the tendency of biofilters to retain older water might enhance their treatment performance, for example by increasing the opportunity for bacterial removal in the rhizosphere by die-off, grazing by predacious protozoa, and other processes (Chandrasena et al., 2017;Surbeck et al., 2010).
PARKER ET AL.

10.1029/2020WR028579
18 of 24 More generally, the retention and release of older water could impact the quality of water released from these systems either positively or negatively, depending on interstorm pollutant transformation mechanisms. For example, between storms the biofilter media's organic material can be respired by resident bacteria, potentially leading to the liberation of ammonium (by ammonification) and nitrate (by nitrification) (Canfield et al., 2010). Thus, water retained in the biofilter from penultimate and older storms could serve as a perpetual source of nitrate that is exported during storm events-a pattern often observed in practice (e.g., Hatt et al., 2009;McPhillips et al., 2018). On the other hand, if anaerobic conditions develop between storms (as is likely to occur if the biofilter contains a submerged zone [H. Kim et al., 2003]) nitrate may be further transformed to harmless N 2 gas and, potentially, the potent greenhouse gas N 2 O (McPhillips et al., 2018) by denitrification.
Evaporative concentration is another aging process that could have potentially significant water quality implications. In northern climates, GSI can serve as a conduit through which deicing salts from road runoff are transferred into shallow groundwater and surface waters, thereby contributing to inland freshwater salinization (Snodgrass et al., 2017). Depending on the degree to which vegetation in these systems fractionate salt during transpiration, evaporative concentration of older water in storage could increase the salinity of water released from GSI, making an already serious problem worse.

Limitations of the Uniform Selection Function
A key assumption of the TTD model derived and tested in this study is that all water in storage is randomly selected for outflow regardless of its age; that is, a uniform SAS function was adopted. A primary benefit of this approach is that the resulting age conservation equation is linear and can be solved exactly in the Laplace domain (see Text S4), leading to explicit formulae for the age distribution of water leaving the system by gravitational drainage and the breakthrough concentration of a conservative solute (Equations 7a and 8b, respectively). A potential disadvantage of this approach is that the physical processes by which water and solute are selected for outflow may not be consistent with the choice of a uniform selection function; that is, by adopting the uniform SAS function, we may be misrepresenting the underlying physics of water and solute transport through the biofilter.
On the one hand, our model closely reproduces measured bromide breakthrough curves with only a single fitting parameter-the effective maximum storage. On the other hand, inferred values for this parameter declined over the two sets of experiments conducted here, from 0.73 m in 2018 to 0.42 m in 2019. These maximum storage volumes are two to three times the actual storage volume of the biofilter (0.25 m), consistent with the idea that some of the water exfiltrating from our test cell may have found its way back to the outflow tank, either by recirculating through our test cell or by flowing directly to the collection manifold and, from there, to the outflow tank. We can point to several changes in the experimental system that occurred over the two years (e.g., the "leak" was partially sealed and the biofilter was replanted) but these changes might have altered the effective volume of the biofilter, the way in which stored water in the biofilter was sampled for outflow, or both. In particular, we cannot rule out the possibility that the apparent decrease in the maximum effective volume from 2018 to 2019 is an artifact of assuming a uniform SAS function for both years. This equifinality problem underscores a key research priority going forward: to delineate under what conditions the uniform SAS function is appropriate for modeling the transport and transformation of solutes within GSI. While the analytical results presented in this study are valid only for a uniform SAS function, software tools are available for numerically integrating the age conservation equation under nonuniform storage selection as well (e.g., Benettin & Bertuzzo, 2018).

Conclusions
Unsteady TTD theory directly links the hydrology and treatment performance of GSI. Its practical application therefore requires, as a first step, delineation of the unsteady water balance over the GSI element of interest. In this study, we demonstrate that this first step can be accomplished with a simple bucket model that tracks time varying infiltration, storage, ET and gravitational discharge over a control volume drawn around the biofilter soil media, which in our case was lined with an underdrain open to the atmosphere.
A set of expressions was proposed and tested for the storage-dependence of water fluxes across the control volume surface, including: (1) an empirical relationship for infiltration that toggles between the inflow of storm water (when the biofilter media is partially unsaturated) and the saturated hydraulic conductivity (when the biofilter media approaches full saturation); (2) cPET for ET; (3) Kirchner's power-law model for gravitational discharge (Equation 3b); and (4) a simple flow routing scheme that directs a constant fraction, α, of the gravitational discharge to the underdrain and, from there, to the outflow tank.
Generalizing the hydrologic bucket model beyond the experimental system described here may require modifying one or more of these four relationships, for example by adopting a process-based model (such as the Green-Ampt equation [Green & Ampt, 1911]) for infiltration, accounting for the reduction of ET that occurs as moisture content falls below the incipient water stress (Hess et al., 2019;Zhao et al., 2013), tailoring Kirchner's gravitational outflow expression (Equation 3b) by adjusting the minimum and maximum storage values (S min and S max ), and directly measuring, modeling or controlling water loss to the surrounding media (thus eliminating the need for a flow routing scheme).
In our implementation of Kirchner's gravitational outflow model we set the minimum storage value to zero, S min = 0, for two reasons: (1) our biofilter was never dry enough to warrant the specification of a value for this parameter; and (2) the very large power-law exponent inferred for Kirchner's gravitational discharge expression (g = 5, Equation 3b) ensures that gravitational drainage will be near zero once storage falls below the field capacity of the biofilter media. On the other hand, depending on context, the minimum storage can be equated to: (1) the biofilter media's area-normalized field capacity (e.g., in settings with long antecedent dry periods where supplemental irrigation is applied to maintain plant communities [Ambrose & Winfrey, 2015]); and (2) the area-normalized volume of the submerged zone for biofilters that include this engineered feature (Brown & Hunt, 2011;H. Kim et al., 2003).
In the context of our TTD theory, the maximum storage parameter, S max , represents the area-normalized void volume of soil media through which solute and water circulate before being transferred to the next hydrologic response unit; for example, shallow groundwater in the case of an unlined biofilter or an MS4 drainage network in the case of a lined biofilter (Grant et al., 2013). Our study suggests that this maximum storage may exceed the biofilter's physical void volume, due to lateral exfiltration of water and solutes into the surrounding soil. Roughly speaking, from the results presented here, it appears that a volume loss of around 55% (close to the design goal for GSI of 67% [A. P. Davis, 2008]) increases the effective maximum storage volume by a factor of two to three.
With the unsteady water balance in hand, we next solved the age conservation equation under the assumption that stored water is randomly selected for outflow regardless of its age (i.e., we adopted the uniform SAS function). From this solution explicit expressions were derived for the mean age of water in storage (and leaving the biofilter as ET and gravitational discharge), various age percentiles, as well as the breakthrough concentration of a conservative solute (Equations 8b and 9b). When compared to bromide breakthrough data measured during our field experiments, we find the model over samples young water, either because the uniform SAS function oversamples young water in storage, or because there is simply not enough old water in storage to sample from (Benettin et al., 2013;Harman, 2015).
Given the magnitude of lateral exfiltration in our system, it is unlikely that water entering the outflow tank was selected exclusively from water stored within the physical boundaries of the biofilter test cell. Indeed, when we allow the maximum storage volume of the biofilter to be a free variable, the inferred volumes are 70%-196% larger than the biofilter's void volume, consistent with the hypothesis that exfiltration increases the effective storage experienced by solutes as they transit through the system. The concordance between predicted and measured bromide breakthrough concentrations improves dramatically after taking this extra storage into account (Figures 4b and 4d). Remarkably, the final model-which includes both the unsteady water balance over the biofilter media (Equation 1a) and the convolution integral for solute breakthrough (Equation 9b)-has only one fitting parameter: the effective volume of the biofilter.
Looking forward, TTD theory could inform the design of GSI elements, for example by allowing engineers to tailor passive features of the porous media component of these systems-such as their depth, porosity and saturated hydraulic conductivity-to achieve target mixtures of young and old water in storage and outflow. Opportunities also exist to control water quality "on the fly," for example by embedding TTD theory into control algorithms that dynamically actuate inflows and outflows using so-called smart storm water technology (Kerkez et al., 2016). The parsimony, extensibility and predictive power of TTD theory opens up new possibilities for modeling and optimizing pollutant removal at the scale of individual biofilters, as well as GSI networks and the urban catchments in which they are embedded (Hrachowitz et al., 2016).

Data Availability Statement
All data used in this study are publicly available (http://www.hydroshare.org/resource/b7187928719446bbaf 30d181787efdad).