Predicting Unsteady Pollutant Removal in Green Stormwater Infrastructure with Transit Time Distribution Theory

In this paper, we explore the use of unsteady transit time distribution (TTD) theory to model pollutant removal in biofilters, a popular form of nature-based or “green” stormwater infrastructure (GSI). TTD theory elegantly addresses many unresolved challenges associated with predicting pollutant fate and transport in 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 with or without first-order decay. The solution is calibrated and validated with breakthrough data from 17 simulated storm events (+/bromide as a conservative tracer) at a field-scale biofilter test facility in Southern California. 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 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 treatment credit attributable to GSI, are inexorably linked to the age distribution of water stored and released by these systems.


1
• A solution is derived from transit time distribution theory to model the removal of 29 stormwater pollutants in green stormwater infrastructure 30 • The solution is calibrated and validated with data from 17 simulated storm events 31 at a field-scale test facility in Southern California 32 • The solution reproduces measured breakthrough concentrations, provided that 33 lateral exchange with the surrounding soil is taken into account 34 35 breakthrough concentrations, provided that lateral exchange with the surrounding soil is 48 accounted for. At any given time, according to theory, more than half of water in storage is 49 from the most recent storm, while the rest is a mixture of penultimate and earlier storms. 50 Thus, key management endpoints, such as the treatment credit attributable to GSI, are 51 inexorably linked to the age distribution of water stored and released by these systems. 52

Introduction 69
Green stormwater infrastructure (GSI) provides many benefits beyond the retention and 70 detention of urban stormwater flows (Walsh et  and submerged zone is predicted with the one-dimensional ADE coupled with one or more 142 Figure 1. (a) Schematic diagram of 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), (ii) and (iii)); outflow by evapotranspiration, gravitational drainage through the filter media, and shortcircuiting ((iv), (v), (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).
experiments described later, which involved simulating a sequence of back-to-back 252 storms. Accordingly, for those experiments we approximated ET with an hourly time 253 series of reference crop potential evapotranspiration (cPET) following FAO guidelines 254  and based on measurements at, or nearby, the field site together with 255 plant-specific traits (details in Text S1, Supporting Information (SI)). 256

Numerical Implementation 257
The water balance bucket model was solved by substituting into equation (1a) the above 258 expressions for infiltration (equation (2)) and gravitational discharge (equation (3b)), 259 along with hourly estimates of cPET (Section 2.1.3). The model was then forced with 260 timeseries (sampling frequency ~1 min -1 ) of measured stormwater inflow (Section 3) and 261 numerically integrated following the procedure described in Text S2 (SI). These 262 simulations yielded ~1 min -1 timeseries of infiltration, storage and gravitational discharge 263 over the 17 simulated storm events described in Section 3. 264

Numerical Solution of the Richards Equation 265
To calibrate the gravitational discharge term (Section 2.1.2) and as a check on the bucket 266 model predictions described above, ~1 h -1 time series of infiltration, storage and 267 gravitational discharge were also simulated with the one-dimensional Richards equation 268 The conservation equation's dependent variable, age-ranked storage [L], represents 290 the area-normalized volume of water stored in the biofilter media control volume at any 291 time with ages or younger. Age-ranked storage is defined mathematically as the 292 product of the area-normalized volume of stored water, , and the cumulative 293 distribution function (CDF) for the fraction of stored water with ages less than or equal to 294 ; i.e., the stored water's residence time distribution (RTD), (equation (4b)). As 295 the age variable, , becomes large, the RTD's CDF tends to unity and the age-ranked 296 storage function collapses to the area-normalized volume of water in storage: 297 . The boundary condition (equation (4c)) ensures that no water stored in 298 the control volume has an age less than . The initial condition (equation (4d)) implies 299 that, at time , the volume of "original" water in storage, , has a single age, , 300 where the Heaviside function is denoted by . As applied to biofilters, equation (4a) 301 equates the change of age-ranked storage of water in the biofilter media (left hand side) to 302 the infiltration of stormwater of age (first term on right hand side); outflow of water 303 by gravitational discharge (second term) and ET (third term) with age distributions 304 and , respectively; and aging of water in storage (fourth term). 305 The two CDFs appearing in the outflow terms, and , represent the 306 fraction of water leaving the biofilter as gravitational discharge and ET with ages or less 307 at time . The backward arrows on these CDFs indicate they are "backward TTDs"; i.e., 308 they represent the age distribution of water leaving the biofilter at time . A corresponding 309 set of forward TTDs can be written for the "life expectancy" of water parcels entering the 310 biofilter at time, . The relationship between forward and backward TTDs is given by 311 Niemi's Theorem (Niemi, 1977;Benettin et al., 2015a;. Under unsteady 312 hydrology, the backward TTDs for gravitational discharge and ET are not necessarily 313 equal, nor are they necessarily equal to the RTD of water in storage (Botter et al., 2011). 314

Ranked StorAgeSelection (rSAS) Function 315
As written, equation (4a) is mathematically ill posed because it consists of a single 316 equation with three unknown functions: , , and . This closure 317 problem can be resolved by introducing a new CDF, the ranked StorAgeSelection (rSAS) 318 function, [-], which maps the fraction of outflow with ages less than or equal to 319 (i.e., the CDF form of the backward TTD for discharge or ET) to the fraction of age-320 ranked water in storage with that age or younger "selected" for outflow by either drainage According to equation (6a) age-ranked storage (left hand side) is influenced by the 343 evolving age distribution of both "original" water in storage at time (first term on 344 right hand side) and "young water" that infiltrates during storm events (second term). This 345 solution was numerically integrated (details in Text S4 (SI)) to yield ~1 min -1 timeseries of 346 age-ranked storage in the biofilter, after substituting bucket model simulations for 347 infiltration, , storage, , and gravitational discharge, (Section 2.1). 348

Age Structure of Stored Water in the Biofilter 349
Under uniform selection the backward TTDs for gravitational discharge and ET are equal, 350 and equal to the RTD of water in storage (compare with equation (4b)) : 351 (7a) 352 The 5 th , 50th, and 95 th percentile ages of water in storage and outflow at any time, , can 353 be obtained from equation (7a) by numerically solving the following implicit equations for 354 water age: , , and . The age-355 ranked storage's probability density function (PDF) can be calculated from equation (7a) 356 by differentiation where the symbol denotes the Dirac delta function: 357 The mean age in storage and outflow immediately follows by taking the first moment of 360 the PDF for age-ranked storage: 361 362 (7c) 363 Further details on the derivation and numerical implementation of equations (7b) and (7c) 364 are described in Text S5 (SI). 365

A TTD Theory for Solute Fate and Transport through a Biofilter 366
The concentration of a reactive or non-reactive (i.e., conservative) solute in water leaving 367 the biofilter by gravitational discharge, , can be calculated by convolving the PDF of 368 the backward TTD (equation (7b)) with the concentration of solute, , that entered 369 the biofilter at time, , and exited the biofilter as gravitational discharge at time 370 and age : 371 (8a) 372 Despite its simplicity, this convolution integral incorporates a rich set of processes, 373 including unsteadiness in the biofilter's water balance (e.g., time-varying inflows, 374 outflows, and storage, through the time-evolution of the backward TTD), unsteadiness in 375 the solute concentration entering the biofilter from the ponding zone (through the 376 dependence of on the inflow time, ) and any time-dependent reactions that 377 no bromide was present in the biofilter's original water), setting (because bromide is 393 assumed to be conservative) and using the distributive property of integration, we arrive at 394 the following expression for bromide concentration in water leaving the biofilter by 395 gravitational discharge (details of derivation in Text S6 (SI)): 396 (9b) 397 In deriving equation (9b) we have assumed that plants in the biofilter take up bromide and 399 water in the same proportion, which may not be the case in practice (e.g., if a solute is 400 excluded from plant uptake its pore fluid concentration will increase over time by in situ 401 evaporative concentration (Bertuzzo et al., 2013; Harman, 2015)). However, ET represents 402 a very small portion of the the overall water balance for the field experiments described 403 later (Section 4) and hence in situ evaporative concentration can be neglected in our case. 404 Time series (~1 min -1 ) of bromide breakthrough concentration were simulated with 405 equation (9b) following the numerical procedures described in Text S6 (SI). 406 2018 we discovered that, during construction, a ca., 5 cm diameter hole had been drilled 430 through the base of the cinderblock wall separating our test cell from the adjacent test cell 431 (and through the wall separating the adjacent test cell from the next test cell and so on) to 432 accommodate a buried irrigation pipe. A substantial fraction of stormwater added to our 433 test cell (approximately 50%) laterally exfiltrated to the adjacent cell through this hole. 434

Orange County Public Works (OCPW) Biofilter Test Facility
While not part of our original design, this feature made for a more realistic field 435 experiment, as most operational biofilters undergo at least some degree of subsurface 436 exfiltration (e.g., Brown & Hunt, 2011). Indeed, our exfiltration rate of ~50% is close to 437 the stormwater volume reduction design goal for GSI of 67% (Davis, 2008). To model 438 lateral exfiltration, gravitational discharge from the biofilter was routed as follows. A 439 fixed fraction, , was assigned to the underdrain and the rest, , to lateral exfiltration 440   (Table S1, SI). Across all 17 storms, 508 ET was a minor component of the water budget (< 0.3% of the ~1400 L added per storm, 509 Table S1, SI). Thus, the difference between these inflow and outflow volumes either went 510 to increasing storage or lateral exfiltration to the adjacent test cell (see Text S7, SI). 511 The fraction of inflow volume recovered at the outflow tank is inversely correlated 512 with each storm's antecedent dry period (R 2 =0.82, Figure S1, SI), consistent with the 513 hypothesis that at least some of the unrecovered water goes to storage.  ). Thus, these two parameters are robust to changes in the sequence 542 and length of antecedent dry periods as well as changes in saturated hydraulic conductivity 543 (within each simulated storm sequence, the saturated hydraulic conductivity declined over 544 time, see Text S3, SI). Our power-law exponent is also concordant with values inferred by 545 keeping with the small storage volume of our biofilter (e.g., compared to the volume of 550 water stored in a catchment). In summary, these results support the hypothesis that 551 Kirchner's power-law relationship (equation (3b)) applies at the scale of a single biofilter.

TTD Theory Predictions for Bromide Transport 560
To characterize the transport of solute through the experimental biofilter, we spiked a 561 subset of experimental storms with bromide as a conservative tracer. In 2018, we adopted 562 a semi-periodic study design involving, on each day, a bromide-free "flushing" storm in 563 the morning (orange arrows in Figure 4a) and a bromide-spiked "tracer" storm in the 564 afternoon (black arrows in Figure 4a). By the second day of the storm sequence, the 565 normalized bromide breakthrough curves (BTCs) settled into a periodic pattern, oscillating 566 between 0.3 and 0.6 during the morning and afternoon storms, respectively 567 (black dots in lower graph, Figure 4a). Here, the variable represents the measured 568  tendency to overshoot bromide measurements implies it is oversampling young water; i.e., 579 the predicted bromide BTC contains too much bromide-free water during the bromide-free 580 morning storm, and too much bromide-spiked water during the bromide-spiked afternoon 581 storm. 582 One possible explanation is that the uniform storage selection function, which 583 underpins our model (see equation (9b) and discussion thereof), oversamples young water 584 for gravitational discharge. Alternatively, the storage selection function is fine but there is 585 not enough old water in storage to select from. While the former explanation cannot be 586 ruled out (indeed, the science of selecting rSAS functions is an active area of current 587 research (Harman, 2019)), the latter explanation is compelling for several reasons. First, 588 the void volume of our biofilter (~900 L) is less than the volume of water flowing into the 589 biofilter with each experimental storm (~1400 L). Therefore, as far as the model is 590 concerned, the biofilter has very limited capacity to store older water from penultimate 591 and older storms. Second, a substantial fraction (>50%) of the inflow volume leaves our 592 biofilter by lateral exfiltration. As noted in Section 4.1, some of this exfiltrated water may 593 eventually return to the outflow tank and thereby increase the effective volume that solutes 594 experience as they transit through the system. 595 To test the last hypothesis-that the effective volume for solute transport is larger 596 than the biofilter's physical volume-we split the measured bromide data from 2018 into a 597 calibration period (the first bromide-spiked storm, storm #2) and a validation period (all 598 other storms). We then inferred a value of the biofilter's void volume by minimizing the 599 root-mean square error (RMSE) over the calibration period ( Figure S3

Age Distribution of Water Leaving the Biofilter by Gravitational Discharge 624
What can TTD theory tell us about the age structure of water leaving the biofilter by 625 gravitational discharge? By selecting a uniform rSAS function for our model (Section 626 2.3.2), the backward TTDs for gravitational discharge and ET are equal to the RTD of 627 water stored in the biofilter. Thus, under uniform storage sampling, the age distribution of 628 water in storage is equal to the age distribution of water leaving the biofilter as 629 gravitational discharge. 630 During the 2018 experiments, predictions for the median age of water stored in the 631 biofilter (equation (7a) and discussion thereof) follows a semi-periodic pattern, increasing 632 linearly with time between storms (as water stored in the biofilter ages) and rapidly 633 declining to near zero during storm events (as incoming stormwater, of age h, fills the 634 biofilter, Figure 5a). The 5 th and 50 th (median) age percentiles overlap but the 95 th age 635 percentile is much older, indicating that the age distribution is positively skewed (Ang & 636 Tang, 2007). The 5 th and 50 th percentiles overlap because, at any time , more than 50% 637 of water stored in the biofilter is from the most recent storm with an age roughly equal to 638 the antecedent dry period. The 95 th percentile age is much older because the rest of water 639 in storage (i.e., water not from the last storm) is from penultimate and earlier storms. 640 For the simulations presented in Figure 5a we arbitrarily set the initial age of 641 "original" water (i.e., water that was initially present in the biofilter at time, ) at 642 50 h. Until the fourth storm, this original water constituted more than 5% of water stored 643 in the biofilter, as evidenced by the upward slope of the grey band in Figure 5a (the 644 upward slope reflects the fact that that original water in storage is aging linearly with 645 . After the fourth storm, the original water's contribution to storage drops below 5%, 646 as evidenced by a steep drop in the 95 th percentile around 28 h (Figure 5a). Thus, four 647 storms were required to flush out 95% of the original water, even though more than 50% 648 of water in the biofilter, at any given time, is from the last storm. A similar pattern is 649 evident for the set of experiments conducted in 2019 (Figure 5b). Across both years the 650 mean age is 5 to 20 hours older than the median age, consistent with a positively skewed 651 age distribution (Ang & Tang, 2007). 652 represents the fraction of water in storage that is younger than the oldest water from the 661 storm indicated. The lower bound of the same color band represents the fraction of water 662 in storage that is younger than the oldest water from the next storm, and so on. 663

Mapping out the Contribution of Past Storms to Present
The influence of biofilter hydrology on the age structure of stored water (and by 664 implication the age structure of water leaving the biofilter by gravitational drainage under 665 uniform sampling) is striking. During the Filling Phase of each storm (e.g., Storm #3 in 666 Draining Phase, the age structure of water in storage does not change (i.e., during this 677 phase all boundaries in Figure 6 are horizontal lines). Figure 6 also vividly illustrates the 678 two key attributes of biofilter storage discussed previously: at any time about >50% of 679 water in storage is from the most recent storm while the rest is a mixture of penultimate 680 and earlier storms. 681

Age Structure and Water Quality 682
The age structure of water in storage has significant implications for the treatment credit 683 attributable to, and the pollution exported by, GSI. For example, during the 2019 storm 684 sequence we included a "worst case" scenario (from a water quality perspective) by using 685 a 50:50 mixture of stormwater and raw sewage for one of the storm events, Storm #3. With the unsteady water balance in hand, we next solved the age conservation 745 equation under the assumption that stored water is randomly selected for outflow 746 regardless of its age (i.e., we adopted the uniform rSAS function). From this solution 747 explicit expressions were derived for the mean age of water in storage (and leaving the 748 biofilter as ET and gravitational discharge), various age percentiles, as well as the 749 breakthrough concentration of a solute with or without first-order reaction (equations (8c) 750 and (9b)). When compared to bromide breakthrough measured during our field 751 experiments, we find the model over samples young water, either because the uniform 752 rSAS function oversamples young water in storage, or because there is simply not enough 753 old water in storage to sample from (Benettin et al., 2013;. 754 Given the magnitude of lateral exfiltration in our system, it is unlikely that water 755 entering the outflow tank was selected exclusively from water stored within the physical 756 boundaries of the biofilter test cell. Indeed, when we allow the volume of the biofilter to 757 be a free variable, the inferred volumes are 70% to 196% larger than the the biofilter's 758 void volume, consistent with the hypothesis that exfiltration increases the effective storage 759 experienced by solutes as they transit through the system. The concordance between 760 predicted and measured bromide breakthrough concentrations improves dramatically after 761 taking this extra storage into account (Figures 4b and 4d) Table S1 page 2 Introduction. This Supporting Information includes a table (Table S1) summarizing water volume recovery and bromide mass recovery for storm simulation experiments in 2018 and 2019. Table S1. Total inflow volume delivered to the biofilter's ponding zone, total volume recovered in the outflow tank, total ET (estimated by integrating the hourly cPET from the beginning of the storm to the beginning of the next storm (or in the case of 2018 Storm 10 and 2019 Storm 7, integrating hourly cPET over ten hours starting from the beginning of the storm)), and mass of bromide added and recovered for all storms simulated in 2018 and 2019. The inflow bromide mass column is marked not applicable (N/A) for storms that were not spiked with bromide, and the outflow bromide mass column is marked N/A for any storms that were conducted prior to the first bromide spiked storm of each year's sequence. Introduction. This Supporting Information includes a detailed description of how we estimated crop potential evapotranspiration (Text S1). Next, we describe our numerical evaluation of the biofilter bucket hydrology model (Text S2) and describe how we estimated saturated hydraulic conductivity (Text S3). We also provide detailed derivations (and illustrate their numerical application) for age-ranked storage under uniform selection (Text S4), the age structure of water in the biofilter (Text S5), and TTD theory predictions for reactive solute transport (Text S6). We then describe how the fraction of inflow that undergoes lateral exfiltration was estimated using

Text S1. Estimating Crop Potential Evapotranspiration (cPET)
Crop Potential Evapotranspiration (cPET) was estimated outside of Hydrus using the Penman-Monteith equation ; FAO-56): (1-1) Here, / = net radiation (calculated, not measured, as a function of albedo, latitude and longitude, and (2-1a) (2-1b) Terms on the righthand side were determined as follows: • Inflow to the ponding zone, , was measured based on the change in weight of the inflow tank (Section 3.1, main text); • The maximum storage, m, was calculated from the depth and porosity of the experimental biofilter, or inferred by minimizing the root mean square error between model predicted and measured bromide concentrations (Section 4.4, main text); • The saturated hydraulic conductivity, , was determined from in situ measurements and measurements of peak outflow (Text S3); • The gravitational discharge parameters, m and , were estimated by calibrating Kirchner's recessional streamflow model to Hydrus 1D simulations (Section 4.2, main text).
• cPET was estimated from local environmental measurements using FAO guidelines (Text S1).
• Finally, the initial storage, m, is equal to Hydrus 1D default values for the field capacity of loamy sand (assuming a depth and porosity of m and ).
Numerical integration of equation (2-1a) was carried out using the NDSolve command in Mathematica as follows (in this particular implementation there are two saturated hydraulic conductivities: ksat1 (which applies from t=0 to t=tswitch) and ksat2 (which applies from t=tswitch to t=tend), and sIC is the initial storage):

Text S3. Estimating Saturated Hydraulic Conductivity from Measured Biofilter Outflow
During the set of experiments conducted during the summer of 2018, peak flow rates measured at the outflow tank systematically declined from approximately 0.11 m h -1 (first two storms) to 0.09 m h -1 (next four storms) to 0.08 m h -1 (final four storms). For the set of experiments conducted in 2019, the peak flow rates declined from approximately 0.12 m h -1 (first two storms) to 0.08 m h -1 (final five storms).

Text S4. Derivation and Implementation of TTD Theory Solutions Under Uniform Selection.
Derviation of the Solution for Age-Ranked Storage. Substituting the uniform rSAS function into the conservation equation for age-ranked storage (equations (5b) and (4a), respectively, in the main text) we arrive at the governing partial differential equation to be solved where the superscript "U" denotes a uniform rSAS function: (4-1a) (4-1b) (4-1c) (4-1d) (4-1e) (4-1f) The boundary condition (equation (4-1f)) ensures that no water volume stored in the biofilter has an age less than or equal to zero, . For the initial condition, equation (4-1d), all of the volume of water (4-5b) Evaluation of equation (4-5a) requires estimating the functions and (see equations (4-3b) and (4-3c)) and the integral term on the right hand side of equation (4-5a). We consider these in turn. As illustrated in Figure 4-1, the function, , is highly oscillatory, taking on large values during storm events and approaching zero between storms. This oscillatory behavior implies that standard numerical algorithms for integrating are slow to converge and prone to numerical error (results not shown). We resolved this problem by expressing the integral as a rate equation, which was then numerically integrating with the NDSolve command in Mathematica; e.g., the code below took 32 microseconds to evaluate. Likewise, our 2019 time series for the function is highly oscillatory due to diurnal variation in solar radiation and hour-to-hour variations in vapor pressure deficit, cloud cover, and wind ( Figure 4-2). We therefore evaluated the function, (equation (4-3c)), using the same approach described above; namely, we expressed the integral as a rate equation which was numerically integrated with the NDSolve command in Mathematica; e.g., the Mathematica code illustrated below took 34 microseconds to evaluate.
Because infiltration is also highly oscillatory (taking on large positive values during storms, and equal to zero in the periods between storms) the integral term appearing on the right hand side of equation (4-5a) was evaluated using the same approach described above (i.e., expressing the integral as a rate equation and then integrating with the NDSolve command in Mathematica). The Mathematica code below took 88 ms to execute.
The Mathematica code below was used to evaluate the overall solution for age-ranked storage: Testing the Exact Solution for Age-Ranked Storage. It is easily verified that Equation (4-5a) satisfies the initial condition: . Next we demonstrate that this solution collapses to total storage when the water age, , is set to a value larger than the oldest water in storage: . At any time, , the oldest water in the biofilter is the initial age of the original water in the biofilter (at time ) plus the time that has elapsed since the start of the experiment: . The solution for ageranked storage was used to calculate the volume of water in the biofilter with age ; i.e., one hour than the oldest water in the biofilter, assuming (Figure 4-3).  (4-5a), for the choice of water age older than any water in storage, , where the initial age is . The corresponding probability density function (PDF) can be obtained by differentiating the CDF form of the RTD with respect to the age variable: (5-1) Combining equations (4-5a) and (5-1) we arrive at the following expression for the PDF form of the biofilter storage RTD under uniform sampling: Taking into consideration the two lower limits for (see equation (4-5b)) and applying Leibnitz integration rule, the second term in the bracket can be simplified: Thus, the PDF form of the biofilter's RTD can be written as follows (compare with equation (7b) in the main text): (5-2) The first moment of equation (5-2) yields the expected value, or mean, of the age distribution for water in storage under uniform sampling: (5-3a) (5-3b) The first integral on the righthand side of equation (5-3b) can be simplified after applying the combing property of the Dirac Delta function: Furthermore, the Heaviside function appearing in the second integral can be removed as follows: Back substituting these results into equation (5-3b) we arrive at the final expression for the mean RTD under uniform storage selection (equation (7c) in the main text): (5-3c) Equation (5-3c) was numerically implemented by expressing the integral as a rate equation (see earlier) in Mathematica as follows:

Text S6. Uniform Selection Predictions for Solute Concentration in Gravitational Discharge
Our solution for age-ranked storage under uniform selection can also be used to calculate the concentration of a reactive solute in the biofilter outflow, , by convolving the PDF form of the backward TTD with the concentration of solute, , that entered the biofilter at time, , and left the biofilter as discharge at age : If the solute in question undergoes first-order reaction (with rate constant ) in the biofilter, the function takes on the following form: (6-1b) Substituting equations (5-2) and (6-1b) into equation (6-1a) we arrive at the following solution for solute concentration in water discharged from the biofilter: From the combing property of the Dirac Delta, and letting represent the solute concentration in "original water", the first integral simplifies a follows: The Heaviside function appearing in the integrand of the second integral can be removed: Thus, solute concentration discharged from the biofilter, assuming uniform sampling, first-order decay in the biofilter, and no solute mass in the original water, becomes equal to equation (6-2) (compare with equation (8c) in the main text): (6-5a) (6-5b) The Mathematica code for implementing equation (6-5a) involves two steps; first numerically evaluating the integral and then evaluating the concentration in gravitational discharge from the biofilter at any time:

Text S7. Experimental Approaches for Estimating Lateral Exfiltration from the Biofilter Test Cell
A summary of the process we used to estimate key parameters for the Hydrus 1D and bucket model water balance calculations is presented in Figure S4.  Table S1) that not all inflow to the biofilter was routed to the outflow tank. We therefore let the fixed constant, , represent the fraction of total gravitational discharge routed to the outflow tank: . We determined using two different approaches: (1) a steady-state experiment conducted during the summer of 2018; and (2) an analysis of the fraction of inflow volume recovered at the outflow tank for transient experimental storms simulated during the summers of 2018 and 2019. These two approaches are described next.
In the first approach, we conducted an experiment in July 2018 (using the same tank and scale experimental set-up described in Section 3.1 main text) in which we supplied a steady inflow of water to the ponding zone of the biofilter (10.1 L/min; the inflow tank was modified to maintain constant head) and simultaneously measured the outflow rate at the outflow tank. The outflow rate eventually stabilized at 5.6 L/min, implying a recovery of .
The second approach involved, for each experimental storm in 2018 and 2019, calculating the total volume collected at the outflow tank and dividing it by the total volume of inflow delivered to the biofilter to obtain the storm's "volume recovery fraction". The volume recovery fraction values thus obtained increase inversely with antecedent dry period (ADP) (R 2 =0.82, Figure S1), presumably because more of the inflow water went to increasing storage during storm events with longer ADPs. By this logic, extrapolating the fractional volume recovery back to an ADP of zero hours should yield the fraction of inflow routed to the outflow tank (under the assumption that none of the inflow water is going toward increasing storage when the ADP is zero). When this exercise is carried out for the set of storm experiments conducted during the summer of 2018 and 2019 we obtain a y-intercept of ( Figure S1), which compares closely with the average volume recovery estimated during the steady-state experiments described above ( , see red horizontal line in Figure S1).

Text S8. Correcting the Time Stamps of Measured Outflow Rate and Bromide Breakthrough
Water leaving our biofilter through the underdrain had to flow by gravity down a manifold to an underground sump where it was pumped into the outflow tank (see Figure 1c in the main text). Our modeling framework was applied to a control volume that included the biofilter media. Consequently, in comparing model predictions measurements of volume and bromide concentration at the outflow tank, we needed to account for the time it takes water to transit from the biofilter underdrain to the outlet tank. We estimated this travel time experimentally by setting up a steady inflow of bromide-free water to the biofilter (using the constant head approach described in the last section). Once steady-state flow had been achieved we poured ~1L of bromide solution (NaBr dissolved in water) into the lateral connecting the biofilter's underdrain to the manifold. Water samples were collected at the outflow collection tank every minute for 14 minutes and measured for bromide concentration using an Orion bromide electrode (Thermo Fisher Scientific, Waltham, MA). From the bromide breakthrough curve, we estimated that the travel time from the biofilter underdrain to the outlet tank is about 3 -12 minutes. As shown in Figure S5, our simplified infiltration model tended to over-estimate infiltration rates into the biofilter during the Second, averaged across all 2018 experimental storm events, ~ 50% of the water added to the biofilter was recovered at the outlet tank (see Text S6). As evapotranspiration (ET) was a relatively minor component of the overall water balance (<1%, see Table S1) we suspected that water was escaping the biofilter test cell into adjacent cells or the surrounding soil; i.e., exfiltration was occurring. Following the 2018 experiments and about five months prior to the 2019 experiments, we removed all plants from the biofilter, excavated the media and discovered that, when the test cell was originally constructed, the contractor had drilled a ca., 5 cm diameter hole into the base of the cinderblock wall separating our test cell from the adjacent cell to accommodate a buried irrigation water supply pipe. The hole was partially sealed (to the extent possible while leaving the irrigation supply pipe in place), the media was replaced, and the biofilter was replanted with a native sedge, Carex spissa, recommended for biofilters in southern California (County of San Diego, 2014).
Third, the operation of the control valve differed over the two years. In 2018, storms consisted of tap water and the discharge rate from the inflow tank through the control valve was continuously adjusted with a real-time control algorithm based on continuous measurements of water level in the inflow tank.
This set-up could not be replicated in 2019, because the 50:50 mixture of stormwater and sewage applied during the third storm required vigorous mixing of the inflow tank, which interfered with the water level sensor. Consequently, during the 2019 experiments we operated the control valve manually, according to a schedule that roughly conformed to the design storm hydrograph. The storm hydrographs (estimated by differentiating the measured decline in weight of the inflow tank over time) generated using this approach compare reasonably well to the design storm, although the peak flow is reduced in both years, and the storm hydrographs in 2019 are smoother and more reproducible from storm-to-storm ( Figure S2).