Sizing an open-channel woodchip bioreactor to treat nitrate from agricultural tile drainage and achieve water quality targets

Nitrate contamination of surface and groundwater is a pervasive issue in areas where intensive agriculture uses fertilizer to increase crop yields (Wang & Li 2019). Management practices that match fertilizer and water applications to meet crop needs will help reduce the infiltration of nitrate to groundwater and runoff to surface water; however, 100% efficiency of nutrient uptake is not possible and about 80% recovery of nitrate by the crop is the practical upper limit (Harter et al. 2012). As nitrate is highly water soluble, any excess nitrate not taken up by crops is conveyed to ground or surface water where excess concentrations can cause environmental and human harm. California's Environmental Protection Agency sets limits for contaminants based on assessments of beneficial uses of a waterbody. These beneficial uses include drinking water, recreation, wildlife habitat, fisheries and other uses, with allowable limits for nitrate concentration varying between 1 and 10 mg N L⁻¹ on California's central coast (CCRWQCB 2017). Regional Water Quality Control Boards are implementing a phased and progressive approach to agricultural compliance with achievement of water quality objectives (WQOs) at the level of individual farms (SWRCB 2018; CCRWQCB 2020). This will put California growers in the position of needing to implement management practices both to increase the efficiency of fertilizer use and to remove excess nutrients that the crop cannot absorb. A relatively low-cost and focused method of removing excess nitrate in runoff is denitrification in woodchip bioreactors.

Tile drain systems are often installed beneath agricultural fields to dewater perched groundwater and direct that water to agriculture ditches, as occurs on lower elevation fields on California's central coast. Because these tile-drained systems are typically pumped, a woodchip bioreactor at the edge of the ditch could intercept the pumped water prior to discharging it to the ditch. Although woodchip bioreactors (WBR) have been designed for treating tile drain runoff in farmlands of the Midwest, the design and sizing concept in that region is based on gravitational flow through a subsurface structure filled with woodchips and a weir that allows for controlling retention time (Schipper et al. 2010a; Christianson et al. 2012). For fields with perched water tables, Hartz et al. 2017 proposed that an above-ground WBR sized at 61 m by 16.8 m by 1.8 m would treat a 200-acre field producing 250 m³ of water with concentrations of Nitrate-N of 60 mg L⁻¹ being removed to achieve the drinking water standard of 10 mg L⁻¹. This work raises the need for further validation and development of a model for WBR sizing. On California's central coast, where land prices are high, farmers pursuing a WBR solution need to meet regulatory nitrate WQOs while compromising the least possible area of arable land. A simple effective method to determine the minimum necessary size of a WBR treating tile drain systems, such as the sump pumped systems, has not been previously developed and should be based on locally determined WBR nitrate decay rates to provide effective design guidance.

Correctly sizing a WBR to meet an outlet concentration target involves understanding the hydrology of the woodchip channel, the reaction kinetics of nitrate removal in the woodchip substrate and nitrate concentrations in the WBR's input to determine the retention time required to meet WQOs (Kadlec & Wallace 2009; Schipper et al. 2010b). Reaction kinetics remain under investigation for WBRs. Findings from the literature suggest zero order (with a constant rate of nitrate removal independent of incoming nitrate concentration), first order (with a constant proportion of incoming to outgoing nitrate removal) or Michaelis–Menten kinetics (first order at low input concentrations and zero order at high input concentrations). Ghane et al. (2015) found that Michaelis–Menten best described their results, whereas Leverenz et al. (2010) and Chun et al. (2010) found that first-order kinetics best described removal and Robertson (2010) concluded that zero-order reaction kinetics were the best fit. In a review of WBR performance, Schipper et al. (2010a) suggested that denitrification likely follows Michaelis–Menten kinetics; however, because NO₃ concentrations entering WBRs are typically high, zero-order kinetics are functionally appropriate. The different findings of these different studies indicate the need for further evaluation in order to gain confidence in a sizing model.

A further complication in sizing a WBR arises because denitrification is a microbial process, thus reaction rates vary with temperature (Schipper et al. 2010a). Temperature dependence is exemplified by the many seasonal studies of woodchip bioreactors demonstrating lower reaction rates in winters compared with summers (Robertson & Merkley 2009, Schipper et al. 2010b; Miller 2014; Pfannerstill et al. 2016). The Arrhenius equation has been widely used to describe the temperature dependence of microbial growth and reaction rates (Peleg et al. 2012; Noll et al. 2020). In its simplified form, the Arrhenius equation has been used for modeling both wetland and bioreactor nitrate removal rates: k_v = k * θ ^ (T – To) (Kadlec & Wallace 2009; Metcalf & Eddy Inc. et al. 2014; Halaburka et al. 2017). If WBR temperatures are variable, as is normally the case in field settings, kinetic modeling should include the temperature dependence of the reaction rate, for which the Arrhenius equation is ideally suited.

In addition to reaction chemistry, understanding hydrology is critically important to sizing an effective treatment system; however, the importance of hydrology is frequently neglected and this can lead to design inefficiencies (Kadlec & Wallace 2009). A model commonly used for hydrology and residence time distributions in environmental transport modeling is the tanks-in-series (TIS) model (Clark 2009; Kadlec & Wallace 2009; Krone-Davis et al. 2013). Because nitrate removal depends on treatment time, a TIS model accounts for the different residence times a given parcel of water experiences. Tanks-in-series models have been used for successfully modeling nitrate removal in both surface flow and subsurface flow wetlands (Kadlec & Wallace 2009); however, this use has not been extended to WBRs. WBRs have more typically been modeled using Darcy's Law or Forchhemier's equation with dispersion coefficients (Ghane et al. 2015; Halaburka et al. 2017). The advantage of the TIS model is that it has a single parameter and transport is modeled with a simple gamma function (Clark 2009).

Clearly, there is a great deal that can be refined in the use of WBRs, but we concentrate here on a better understanding of the reaction kinetics and hydrology of WBRs in order to develop a simple approach that can be used for estimating size to meet a targeted outlet concentration. Using data on inputs and outputs from a field-scale WBR over a period that included the entire crop cycle, we modeled the process of denitrification in the hope of improving our understanding of this valuable agricultural practice and providing a basis for responsible and efficient WBR sizing.

• 1.

  To determine whether zero- or first-order reaction kinetics best describe woodchip bioreactor performance.

• 2.

  To determine the Arrhenius constant and reduction rate coefficient values that best describe nitrate removal rates for woodchip bioreactors scaled to be applicable to the climatic and agricultural context of California's Central Coast region.

• 3.

  To develop a simple approach to estimate woodchip bioreactor size for pumped tile-drained systems to achieve outlet WQOs.

The Central Coast Wetlands Group (CCWG) at Moss Landing Marine Laboratories and California State University Monterey Bay (CSUMB) built a multi-channel experimental water treatment system in 2016 near Castroville, California on a private utility company property adjacent to agricultural land. The experimental system is comprised of a series of 12 parallel channels of approximately equal size lined with impermeable plastic designed to treat agricultural runoff collected from the Castroville Ditch. The Castroville Ditch receives water from 1080 acres of irrigated farmland, predominantly under vegetable production that includes artichokes, Brussel sprouts, lettuce, and celery. Water from the ditch was pumped into the experimental channels using one of the many regional pumps installed to lift tailwater from near sea level up to elevations from which gravity flow toward the ocean proceeds. Given the widespread existing deployment of these pumps systems, there is great potential to achieve in-line treatment by installing reactors immediately below the pumps. Without the pumps, there would not be enough hydraulic head to drain the tailwater or achieve meaningful flow through reactors in this flat terrain. The channels were sized to be able to remove a substantial amount of nitrate from a typical 6.5 acre farm block and to be able to be installed in a relatively compact footprint below pump stations and alongside fields. Reactor length was maximized within the available area to promote maximum residence time. Reactor depth was maximized to reduce site footprint, but constrained to fit vertically between typical pump head elevation and downstream tailwater elevation. The resulting dimensions for each channel were 23 m long by 1.7 m wide by 0.45 m deep. The 12-channel system allows for three replicates of as many as three treatments and one control. However, for the purposes of this study, we were interested in modeling the performance of a single woodchip-filled channel.

The experimental channels receive tile-drained agricultural water pumped from the Castroville Ditch into a sediment basin immediately upstream of the channels with a retention time of approximately 1–2 days. Water from the sediment basin gravity drains into a plastic (HDPE) pipe that acts as a manifold for distributing water to each of the 12 channels (Figure 1). Inlet structures to the individual channels are engineered V notches that can be raised and lowered to change the flow rate and retention time of water independently in each. Water flows by gravity through each of the channels and then exits an outlet drain into a downstream treatment wetland. The multi-chamber flow system is managed by the CCWG and visited approximately once per week to measure flow rates, remove biofouling and make V notch adjustments to achieve the desired hydraulic retention time. Additional water parameters were monitored quarterly using a calibrated YSI multi-parameter sonde to measure dissolved oxygen (inlet range: 0.52–6.34 mg L⁻¹, mean 1.47 mg L⁻¹; outlet range: 0.36–12.30 mg L⁻¹, mean 6.58 mg L⁻¹), pH (inlet range: 7.00–8.01 mg L⁻¹, mean 7.64 mg L⁻¹; outlet range: 6.85–7.82 mg L⁻¹, mean 7.51 mg L⁻¹) and conductivity (inlet range: 367–3,788 μS cm⁻¹, mean 2,736 μS cm⁻¹; outlet range: 542–3,840 μS cm⁻¹, mean 2,563 μS cm⁻¹).

Woodchips to fill the channel were purchased from the local waste management facility and consisted predominantly of pine and oak with size variation between 0.2 to 2 cm wide by 3 to 7 cm long. Eucalyptus and pressure-treated wood were excluded from the source materials. Woodchips were added to the channels following their construction in 2016 and then supplemented in March 2018 when noticeable shrinkage had occurred. The woodchip addition represented less than 10% of the original volume. Drained porosity of the woodchips was measured at 50%.

Nitrate data were collected every 2 hours from 12/13/18 through 1/14/20 except for several data gaps when the bioreactor or Hydra monitoring array were undergoing maintenance or repairs. A custom-built water monitoring system (Hydra) integrated with a Seabird SUNA V2 optical nitrate sensor measured water temperature and nitrate concentrations of the inlet source water and the 12 outlets of the multi-channel system every two hours. The Seabird nitrate probe is cleaned weekly and calibrated every 6 months. Nitrate and temperature data from Hydra sensors were downloaded weekly, reviewed, subjected to quality assurance protocols, and saved. Hydra nitrate concentration measurements were compared with weekly grab samples analyzed at the MLML Nutrient Analytical Laboratory to verify accuracy. The laboratory measures nitrate concentration using colorimetric analysis by a Lachat Quickchem 8000 flow injection autosampler in accordance with standards in the Irrigation and Nutrient Management Quality Assurance Project Plan (QAPP).

From the 13-month dataset, we selected data from a number of time periods when the bioreactor channels and the water monitoring system were functioning properly. We chose long periods (>140 hours) of continuous data so that we could compare inflow and outflow over many data points. Periods with variability in nitrate inlet concentration and water temperature were selected to assess reduction rates under a range of conditions and to support the calculation of decay rate parameters for the model (see below). We omitted periods with inaccurate or incomplete nitrate or temperature data or with high variability in inflow rate. We also omitted data when outlet concentrations were zero, as inclusion of these points would not provide information about decay rate.

where is a gamma distribution of hydraulic residence times, with shape and scale parameters k_R and θ_R. N is the hypothetical number of TIS that represents hydrology and is the average hydraulic residence time water receives treatment in the WBR. (More detailed information on this model can be found in Krone-Davis et al. 2013.)

The reference temperature T_o can be set at any value, however the convention for wetlands is 20 °C (Kadlec & Wallace 2009). If a different value for T_o is adopted, this will influence the value of θ, thus T_o is necessarily a selected model parameter used in the model while varying k and θ to evaluate goodness of fit.

We implemented the model using numerical simulation with a two-hour time step for each inlet parcel of water entering the bioreactor, such that fractions of this inlet parcel exited the wetland according to the hydrological distribution function (Equation (2)) and each parcel was thus present in the bioreactor undergoing degradation for the length of time represented by this distribution. We calculated the nitrate removal for each outlet parcel, using the first-order and zero-order decay rate equations over the amount of time that the parcel was present in the WBR. We included temperature dependency in the modeled nitrate removal. We eliminated early outlet model data (1.7 times τ) to allow the model time steps to fill so that representative outlet water was nearly complete (>95%), and we incorporated a prior concentration into the model of the average concentration over the time interval to account for the remaining 5%. To determine predicted outlet nitrate concentration, we summed the partial nitrate contribution of all the parcels of water exiting the WBR at the same time step in the model.

We performed a tracer study of the woodchip-filled Channel 1 on 8/17/19 by adding 160 g of bromide and collecting outlet water at two-hour intervals for future lab analysis. Bromide concentrations collected at the outlet were analyzed using inductively coupled plasma mass spectrometry (Agilent Technologies, 7900 ICP-MS). Outlet flow rate of Channel 1 on 8/19/20 was measured as 0.27 L sec⁻¹. The recovery rate of bromide tracer in Channel 1 was calculated by multiplying outlet bromide concentration less background bromide concentration by flow rate and using the trapezoidal method to calculate mass removed for each two-hour period. The bromide recovery rate in Channel 1 was 108%. We performed a similar tracer test on Channel 5; however, flow rate inconsistencies did not allow evaluation of the tracer test results.

For both reaction types (zero and first order) we identified a range of feasible parameter values for zero and first-order decay rate constants and for Arrhenius θ based on the literature (Table 2). To determine the parameter values that best fit the measured data, we used the field data collected for each of the four run time intervals and ran the model over the feasible range of values for decay rate and Arrhenius theta, separately for the first-order and zero-order models. For each model run we calculated RMSE. As we approached the numeric regions where RMSE was lowest, we used more finely graded parameter values to determine best fit to the outlet nitrate concentration values.

A comparison of model runs over the four time periods showed the zero-order reaction rate model matched the outlet measured nitrate concentration values better than did the first-order reaction rate. The lowest RMSE was 1.145 for the zero-order model and 2.057 for the first-order model. The value for zero-order decay rate that minimized RMSE was 17.5 g-N m⁻³ day⁻¹ an Arrhenius θ of 1.12. The value for first-order decay rate that minimized RMSE were 0.47 day⁻¹ and Arrhenius θ of 1.08. During the time period with the lowest nitrate input concentration, when Michaelis–Menten kinetics might predict better performance for first-order reduction, the zero-order model still outperformed the first-order model (Figure 3). Model performance of the zero-order model showed a good fit of model predictions to the measured outlet concentrations for all four time periods as depicted in Figure 4.

The estimated relationship between WBR performance and temperature is quite strong as exemplified by its exponential form (Figure 5). Using the Arrhenius equation over the range of temperatures observed during our model runs, the decay rate (k_v0) at 23 °C is increased by 7.7 times relative to k_v0 at 5 °C. The temperature dependence is also expressed as a Q₁₀ value, the change in reaction rate for a 10 degree change in temperature; in this case Q₁₀ = 1.76.

We reviewed RMSE over the model runs and graphed them in three dimensions against the parameters we were evaluating (decay rate and Arrhenius theta) to determine that there were no local nadir. The resulting graphics showed a smooth monotonic function, reinforcing our confidence in the parameter values at the lowest RMSE (Figure 6).

Zero-order decay best described nitrate removal compared with first order in the WBR over the range of concentrations and temperatures evaluated in our field study. Several other studies have compared zero-order, first-order and Michaelis–Menten kinetics. Ghane et al. (2015) found that neither zero nor first-order rates explained their data and found that Michaelis–Menten kinetics was a better solution for a small denitrification bed. They determined a K_M value of 7.2 mg N L⁻¹ at V_max/2, where V_max is a zero-order rate constant that represents the maximum removal rate, which was achieved when the inlet nitrate concentration was >30 mg N/L. Halaburka et al. (2017) used a multi-species reactive transport model that included DO and DOC along with nitrate and found that denitrification in a woodchip-filled PVC column at constant temperature with concentrations >2 mg N L⁻¹ could be effectively modeled using zero-order kinetics. Miller (2014), using an Akaike information criterion (AIC) approach and including combined and independent zero-order and first-order models, found that nitrate removal was best supported by the model that contained both zero-order and first-order coefficients. Inlet nitrate concentrations in their study included a lower range (0.39–84.9 mg L⁻¹) compared with ours (6.9–60.8 mg L⁻¹). Leverenz et al. (2010) concluded that first-order kinetics best described their removal rates at inlet nitrate concentrations varying between 9–21 mg N L⁻¹ and temperatures between 10 and 30 °C. However, they noted the need for additional studies to better characterize removal kinetics. In a mesocosm study, Schmidt & Clark (2013) found that denitrification rates were unresponsive to changes in nitrate concentration, therefore matching zero-order kinetics. Schipper et al. (2010a) reviewed studies of WBRs of different sizes, reporting zero-order rate constants between 2 and 22 g N m⁻³ day⁻¹. A plausible reconciliation of findings is that Michaelis–Menten dynamics governs nitrate removal in WBRs, resulting in first-order performance at lower concentrations and zero-order performance at higher concentrations. We suggest that deployments of agricultural WBRs are going to be in response to nitrate concentrations that are above the range where first-order dynamics would be expected and therefor zero-order models are the best choice for estimating size in field applications.

We modeled hydrology and compared model outputs with measured outlet concentrations for the purpose of accurately determining rate coefficients. The TIS modeling approach uses a gamma distribution for the residence time density (RTD) function that describes the hydraulics of the channelized WBRs. The model performed well in depicting the bromide tracer test results that were used for determining the value of N as 7.8 (Figure 1). As portrayed in the four model runs (Figure 4), TIS disperses inlet concentrations in accord with the residence time distribution as they pass through the WBR, such that peaks in input are rounded off at the outlet. TIS accurately predicted this spread and modeled predictions closely paralleled measured values throughout the entire temperature range. TIS has been used for modeling wetlands with subsurface flow, which are similar to WBRs but use gravel rather than woodchips as the media. Typically, subsurface flow wetlands are characterized by 11 TIS (Kadlec 2009), indicating a comparatively lower dispersion or mixing than in the WBR. Although the TIS/gamma distribution approach has not been previously utilized for WBRs, the simplicity and effectiveness of this model makes it attractive.

First-order TIS models have been solved analytically for steady state flow and stable inlet nitrate conditions as described by Clark (2009); however, there is no analytical solution for TIS steady state zero-order model. Because hydrology has an important effect on reactor performance it should not be ignored in sizing calculations when achieving an outlet concentration is important (Kadlec 2009). In the absence of an analytical solution portraying zero-order removal, we developed a TIS Excel version of the model that can be used for estimating WBR size available at github.com and found by searching ‘woodchip bioreactors TIS model’. An example from our model is that with specified input conditions (inlet NO₃-N = 40 mg L⁻¹, water temperature = 18 °C, discharge = 2 gpm and desired outlet NO₃-N<10 mg L⁻¹) the WBR size would need to be 46 m³. Locally relevant values for k and θ can be substituted for those derived in our study, which we recommend for California's central coast. Field measures needed to run the model are inlet nitrate concentration, water temperature, discharge and desired outlet nitrate concentration. Computing size for at least two time periods is informative, one for highest nitrate concentration in runoff and second for highest flow rate.

If the reactor footprint is larger than can be accommodated on the field edge, WBR denitrification can be augmented by the addition of a second carbon source such as glycerol, methanol, acetic acid or vegetable oil to achieve a higher rate of nitrate removal and thereby reduce the volumetric size of the WBR (cf. Shrimali & Singh 2001; cf. Hagman et al. 2008; Hartz et al. 2017; Roser et al. 2018). Hartz et al.’s (2017) study of the addition of ethanol or glycerol to two field WBRs designed for an HRT of 2 days found that nearly complete removal of inlet nitrate concentration of 160 mg L⁻¹ was achieved if a second carbon source was added at the ratio of 1.4:1 (C applied: N denitrified on a mass basis). Roser et al. (2018) found that by spiking woodchips with acetate in a column study, nitrate removal was increased to a rate of 30 g N m⁻³ day⁻¹ at 5 °C and to 121 g N m⁻³ day⁻¹ at 15 °C. Achieving these higher removal rates would greatly reduce the footprint of the WBR. Added costs of carbon augmentation include designing and installing the carbon delivery system to modulate addition based on nitrate concentration and the cost of the carbon source. Regulatory requirements or permits associated with adding a secondary carbon source that could potentially contaminate water exiting the bioreactor would need to be investigated. These increased system and operational costs might offset the cost of land associated with requirements for a larger WBR relying on woodchips alone as the carbon source.

Sizing should also consider annual variation in irrigation runoff, aging of woodchips, variations in crops and fertilizer inputs, and sedimentation of the bioreactor resulting in diminishing performance. Other design considerations include resiliency to stormflows, groundwater interactions, and placement to ensure the farmer can make repairs and replenish woodchips. Finally, management decisions will include consideration of the risk and penalty of not achieving the target nitrate concentration, WBR cost, and alternative management practices that might be utilized to reduce nitrate.

The environmental and human benefits of removing nitrate from ground and surface water include restoring drinking water to safe limits (nitrate < 10 mg N L⁻¹), reducing the occurrence of harmful algal blooms (HABs), and preventing excessive algal growth in surface waters that can produce hypoxic states and change the ecology of the rivers, estuaries and bays (Wurtsbaugh et al. 2019). There are simultaneously negative consequences of the denitrification process that can occur in WBRs that can be partially managed by their design and operation. WBRs produce high levels of biological oxygen demand (BOD) resulting in low levels of dissolved oxygen during their first 1–4 months of employment and immediately following the addition of more woodchips (Cameron & Schipper 2010). To avoid negative impacts on downstream biology, this should be considered and mitigated if necessary. Nitrous oxide and methane are powerful greenhouse gases that can be produced in WBRs. Researchers have reported nitrous oxide emissions from denitrification beds ranging between 0.12 and 78.58 μg N m⁻² min⁻¹ (Warneke et al. 2011; Ghane et al. 2015). Nitrous oxide is produced when the denitrification process is incomplete because of high levels of nitrate and low amounts of bioavailable carbon, and thus are partially manageable (Aalto et al. 2020). When retention times are long and low redox potentials occur, hydrogen sulfide also can be produced in WBRs (Robertson et al. 2009; Shih et al. 2011). WBRs can additionally produce methyl mercury when sulfide reducing conditions develop, however this process can be suppressed by maintaining nitrate concentrations >0.5 mg N L⁻¹ (Shih et al. 2011). WBRs should ideally be sized to achieve the needed nitrate removal and designed to avoid or minimize these negative consequences. Due to variations in inlet nitrate concentrations and flow rates, the addition of weirs, bypasses or other control features may sometimes be necessary to control for these conditions. Similarly, designing WBRs in parallel or series will allow for operators to take one or multiple sections of the treatment system off line if they are not needed or require maintenance.

This study corroborates the findings and assumptions of other WBR research that denitrification is explained by zero-order kinetics over the range of inlet concentrations that is of concern in agricultural runoff. Because we evaluated a high number of inlet-outlet concentration data points (n = 264) that included both winter and summer conditions we were able to estimate the values for nitrate decay (k₀ = 17.5 g N m⁻³ day⁻¹) and reaction rate temperature dependency (θ = 1.12) over a substantial variety of conditions, thus increasing our confidence in the derived values. We developed an Excel version of a TIS model that can be used to estimate WBR sizing where local decay rate parameters can be inputted along with field conditions for flow rate, water temperature and nitrate concentration.

Woodchip bioreactors have demonstrated the capacity for substantial nitrate removal from agricultural tile-drained water and runoff over a variety of designs, including subsurface flow, tile drain sidewall, in-streambed, and in-channel alternatives (Blowes et al. 1994; Jaynes et al. 2008; Robertson & Merkley 2009; Schipper et al. 2010a; Pfannerstill et al. 2016). WBRs have been more widely studied in the Midwest over the past 15 years as subsurface flow systems intercepting tile drain water. Considerable effort has been invested in methods for determining the size and design for these subsurface systems, and a practice standard has been adopted by the Natural Resources Conservation Service (NRCS), which is evidence of their utility (USDA NRCS 2009). Development of a design concept for sizing a pump assisted open-channel WBR to treat tile-drained water is rare (Hartz et al. 2017); yet pumped tile-drained systems are common on California's central coast. Being able to size a WBR can assist an individual grower or watershed collective attempting to meet regulatory goals for water quality. This same framework can provide a basis to estimate credits for nitrate removal in some regulatory frameworks or for nitrate trading strategies. Further design and operational considerations for WBRs will need to include an additional carbon source, varying hydrology, tile drain outfalls, variation in flow rate, proximity to neighborhoods, regional regulations, and other factors. By using existing drainage features and including design concepts that allow for operation and management of the system under the appropriate range of conditions, farmers and technical assistance providers on California's central coast can design relatively low-cost, simple systems that can be sized to achieve WQOs. The woodchip bioreactor TIS model available on line at github.com can provide an initial estimate of the size needed.

We thank David Winkler for help with 3-D graphics and general editing.

Financial support for research and preparation of the journal article for this grant work was funded by a USDA Conservation Innovation Grant (Agreement # 69-3A75-17-27), CDFA Specialty Crop Block Grant (Agreement # SCB14031), California Nonpoint Source Pollution Control Program Implementation Grant (Agreement # SWRCB-D 181330300), and the Anthropocene Institute. Funders were not involved with the study in any way other than financial support.

All relevant data are included in the paper or its Supplementary Information.