Evaluation of Nitrogen and Phosphorus Removal from a Denitrifying Woodchip Bioreactor Treatment System Receiving Silage Bunker Runo ﬀ

Featured Application: Denitrifying woodchip bioreactors are an effective treatment technology for agriculturaltiledrainage.Thisworkfurtherexpandstheapplicationtotreatleachateandstorm-driven runoff from silage storage bunkers. Additional applications of this treatment technology may effectively treat other agricultural wastewaters containing high nutrient concentrations. Abstract: Leachate and storm-driven runo ﬀ from silage storage bunkers can degrade receiving water bodies if left untreated. This study evaluated a novel treatment system consisting of three treatment tanks with a moving-bed bioﬁlm reactor and paired side-by-side denitrifying woodchip bioreactors for the ability to reduce inﬂuent nutrient mass loads. Flow-based samples were taken at four locations throughout the system, at the inﬂow to the ﬁrst tank, outﬂow from the tanks prior to entering the woodchip bioreactors, and from the outﬂows of both bioreactors. Samples were analyzed for concentrations of nitrogen (N) and phosphorus (P) species. Inﬂow concentrations were reduced from the bioreactor outﬂows by an average of 35% for total N (TN) and 16% for total P (TP) concentrations on a storm event basis. The treatment system cumulatively removed 76% of the TN mass load, 71% of the nitrite + nitrate-N (NO 2 − + NO 3 − -N) load, 26% of the TP mass load, and 19% of the soluble reactive P load, but was a source of ammonium-N, based on the monitoring of 16 storm events throughout 2019. While the system was e ﬀ ective, very low NO 2 − + NO 3 − -N concentrations in the silage bunker runo ﬀ entered the bioreactors, which may have inhibited denitriﬁcation performance.


Introduction
The discharge of nitrogen (N) and phosphorus (P) from agricultural production areas due to storm-driven runoff can degrade water quality and lead to eutrophication, cause fish kills, and decrease the esthetic value of water bodies [1]. One specific harmful production area wastewater originates from silage storage bunkers on confinement dairy farms. To store feed for livestock throughout the year, plant material, typically chopped corn or hay, is packed in horizontal concrete floored and walled bunkers and covered with polyethylene tarps. This process, known as ensiling, promotes anaerobic fermentation of soluble carbohydrates by producing organic acids, which preserve proteins in the feed and inhibit bacterial growth that can cause spoilage [2]. Silage leachate is a byproduct from the ensiling process, which is extremely high in nutrient content, consisting of average concentrations of over 3500 mg L −1 for total N (TN), 800 mg L −1 for total P (TP), 5-day biochemical oxygen demand (BOD) of over 60,000 mg L −1 , and a pH between 3.5 and 5 [3]. Controlling the moisture content of ensiled forage by optimizing timing for harvesting plant material is recommended to limit concentrated leachate production [4], but appropriate harvest timing is often difficult to achieve in the Northeastern United States, due to unfavorable drying conditions and a limited harvest period.
During storm events, exposed silage particles and leachate on impervious storage bunker surfaces mix with precipitation and are transported as runoff [4]. While dilution with rainwater yields lower nutrient concentrations in silage bunker runoff more than the concentrated leachate, nutrients are still a major concern if released into surface waters [3][4][5][6]. Runoff concentrations from silage bunkers generally depend on storm size and intensity, concentration of undiluted silage leachate, and seasonality [4]. Best management practices to reduce leachate losses and subsequent environmental degradation are recommended or required on many dairy farms [3].
Management options include intercepting silage leachate so that it can be disposed of properly (e.g., stored, aerated, and diluted to be applied as fertilizer to hay lands or other feed crops during the growing season) [7]. There are minimal data available describing silage leachate composition changes after storage, and more research is needed in this area [3]. For storm-driven silage runoff, vegetative treatment areas (VTAs) are often recommended. VTAs are vegetated systems located downslope of a pollution source designed to treat runoff and provide a buffer zone to the surrounding environment [8]. In a New York study of three VTAs receiving silage bunker runoff, Faulkner [6] found significant ammonium-N (NH 4 + -N) mass removal at all sites and one VTA had 40% soluble reactive P (SRP) mass removal. Incoming nitrate-N (NO 3 − -N) concentrations were low, and some sites experienced a net export, presumably due to nitrification processes within the VTA. Whereas VTAs are generally an economically viable and low maintenance option for farms, several concerns still exist. These issues include surface runoff remaining untreated [9], vegetative burning due to high acidity in the runoff [10], neglected operation and maintenance demands [6], and subsurface leaching that can lead to groundwater contamination [9]. Accordingly, there is a need to explore alternative remediation practices for silage bunker leachate and storm-driven runoff.
In this research study two technologies for managing silage bunker runoff were evaluated. First, a moving-bed biofilm bioreactor (MBBR) is considered, which is a tank-based wastewater treatment technology that incorporates engineered carriers to provide surface area for biofilm growth and which are suspended through aeration [11]. The goal of this technology is to increase the removal of organic matter and nutrients. MBBR systems have simple and compact designs, can operate under varying temperatures, and sustain the growth of microorganisms [12]. Incorporating inoculant at the startup of bioreactors is typically done to establish populations of nitrifying bacteria [13,14]. When designed to accommodate influent wastewater characteristics, sufficient hydraulic retention time [15], and appropriate ratio between carrier volume and tank volume [12], MBBR systems have been shown to effectively treat municipal wastewater [11,16], industrial wastewaters [17], and dairy processing wash waters [18].
The main component of this system is a denitrifying woodchip bioreactor; a passive treatment technology consisting of a bed of saturated woodchips that provide anerobic conditions and a carbon (C) source for microbial denitrification, in which nitrate (NO 3 − ) is converted to inert dinitrogen gas (N 2 ) or, in some cases, nitrous oxide (N 2 O) [19]. Woodchip bioreactors have been shown to effectively improve water quality from wastewater sources containing elevated N concentrations [20]. While most bioreactor studies focus on NO 3 − removal, reductions of other water quality parameters such as P [21], BOD [22], and fecal indicator bacteria (Escherichia coli) have also been reported [23]. Designs that consider variables such as optimal NO 3 − loading rates [24] and hydraulic retention time [25] have been successful in a number of settings treating agricultural tile drainage [19,26], dairy farm wastewaters, [27], aquaculture effluent [28], and septic tank effluent [23]. This system is the first application of a denitrifying bioreactor used to treat silage leachate or silage bunker runoff.
To address the lack of reliable treatment technologies for silage storage bunker runoff, this study evaluated the nutrient removal performance of an innovative treatment system including treatment tanks with MBBR, followed by denitrifying woodchip bioreactors, on a research farm located in Vermont, USA. The objectives of this study were to:

1.
Characterize the N and P composition of storm-driven silage storage bunker runoff.

2.
Evaluate concentration reductions and nutrient mass removal of the two components of the treatment system (i.e., the treatment tanks with MBBR, and the denitrifying woodchip bioreactors) during storm events.

3.
Consider the design of the system and its effectiveness in treating silage storage bunker runoff in order to offer recommendations for future designs.

Study Site Description
A treatment system for silage bunker runoff ( Figure 1) was installed in 2017 at the University of Vermont Paul R. Miller Research Complex in South Burlington, VT, a small dairy farm with 45 milking cows and an equine center used for research and teaching. The site is located within the Potash Brook watershed, which drains to Lake Champlain. The drainage area of the system is 2767 m 2 consisting of silage storage bunkers, where varying amounts of corn and hay is ensiled and stored under plastic tarps throughout the year, as well as a paved asphalt entry area for heavy equipment. Silage leachate and silage bunker runoff flow downslope away from the stored silage and are directed to the treatment system. A diversion berm in the asphalt prevents other production area runoff from entering the treatment system. The system consists of a particle separating screen assembly, a sequence of three treatment tanks for containment of low flows, settling and aeration, woodchip bioreactors, and an earthen infiltration basin ( Figure 2).
Appl. Sci. 2020, 10, x FOR PEER REVIEW  3 of 17 tanks with MBBR, followed by denitrifying woodchip bioreactors, on a research farm located in Vermont, USA. The objectives of this study were to: 1. Characterize the N and P composition of storm-driven silage storage bunker runoff. 2. Evaluate concentration reductions and nutrient mass removal of the two components of the treatment system (i.e., the treatment tanks with MBBR, and the denitrifying woodchip bioreactors) during storm events. 3. Consider the design of the system and its effectiveness in treating silage storage bunker runoff in order to offer recommendations for future designs.

Study Site Description
A treatment system for silage bunker runoff ( Figure 1) was installed in 2017 at the University of Vermont Paul R. Miller Research Complex in South Burlington, VT, a small dairy farm with 45 milking cows and an equine center used for research and teaching. The site is located within the Potash Brook watershed, which drains to Lake Champlain. The drainage area of the system is 2767 m 2 consisting of silage storage bunkers, where varying amounts of corn and hay is ensiled and stored under plastic tarps throughout the year, as well as a paved asphalt entry area for heavy equipment. Silage leachate and silage bunker runoff flow downslope away from the stored silage and are directed to the treatment system. A diversion berm in the asphalt prevents other production area runoff from entering the treatment system. The system consists of a particle separating screen assembly, a sequence of three treatment tanks for containment of low flows, settling and aeration, woodchip bioreactors, and an earthen infiltration basin ( Figure 2).   Storm-driven runoff from the bunker area enters the treatment system ( Figure 2) through a stainless-steel screen assembly that has a series of three screens with mesh holes of decreasing size (7.2 cm 2 , 3 cm 2 , 0.79 cm 2 ) to block silage particles and other large solids transported in the runoff. After the screen assembly, runoff and leachate flow through a series of three treatment tanks that overflow into each subsequent tank during storms. Tank 1 has a storage capacity of 7.57 m 3 and is intended for BOD removal, settling of organic solids, and to prevent clogging of the downstream system components. Tank 2 is a 3.76-m 3 MBBR tank equipped with a regenerative blower, which operates every other hour. This tank is oxygenated to allow for reduction of BOD and N transformations; specifically, mineralization of organic N to NH4 + and nitrification of NH4 + to NO3 -. In May 2019, cylindrical media carriers were added to fill 40% of the Tank 2 size by volume (1514.2 L) which each have a diameter of 2.5 cm, height of 0.4 cm, and a projected surface area of 800 m 2 /m 3 (AnoxKaldness TM -Veolia Water Technologies AB). Tank 2 was inoculated with 18.9 L of a liquid suspension containing ammonia oxidizing microbial strains (MICROCAT-XNL, Bioscience, Inc., Allentown, PA, USA) for establishment of nitrifying bacteria populations. A third 3.76-m 3 tank, Tank 3, is for settling of any remaining solids.
Paired side-by-side denitrifying woodchip bioreactors referred to as East Bioreactor (EB) and West Bioreactor (WB) are each 9.1 m by 12.2 m with a depth of 1.4 m. The woodchip bioreactors are continuously saturated, the bottom, side walls, and center barrier splitting the two bioreactors are lined with 1.1-mm ethylene propylene diene monomer (EPDM) and contain mixed hardwood woodchips (approximate size of individual chips: 5 cm × 5 cm × 0.6 cm). The species composition of the woodchips was 60% Ash (Fraxinus sp.), 20% Yellow Birch (Betula alleghaniensis), and 20% Silver Maple (Acer saccharinum), determined by a local woodchip supplier, and was based on available feedstock species at the time of construction. The bioreactors receive effluent that is split after leaving Tank 3 through two sets of pipes distributed across the top perimeter of each bioreactor. The effluent percolates vertically into the woodchips in each bioreactor. The bioreactors remain saturated and internal water depth is governed by an adjustable water level control structure. When the woodchip bioreactors discharge, their outlets flow through parallel 100-mm diameter underdrain pipes into an earthen infiltration basin. Sampling occurred at four locations: before Tank 1 (System Inflow), after Tank 3 (Tank Outflow, prior to entering the woodchip bioreactors), and at the outflows of both EB and WB ( Figure 2).
The system has three distinct runoff flow paths (Figure 2), depending on the flow rate, which allowed for bypass of some treatment components during larger storm events to protect the system. Storm-driven runoff from the bunker area enters the treatment system ( Figure 2) through a stainless-steel screen assembly that has a series of three screens with mesh holes of decreasing size (7.2 cm 2 , 3 cm 2 , 0.79 cm 2 ) to block silage particles and other large solids transported in the runoff. After the screen assembly, runoff and leachate flow through a series of three treatment tanks that overflow into each subsequent tank during storms. Tank 1 has a storage capacity of 7.57 m 3 and is intended for BOD removal, settling of organic solids, and to prevent clogging of the downstream system components. Tank 2 is a 3.76-m 3 MBBR tank equipped with a regenerative blower, which operates every other hour. This tank is oxygenated to allow for reduction of BOD and N transformations; specifically, mineralization of organic N to NH 4 + and nitrification of NH 4 + to NO 3 − . In May 2019, cylindrical media carriers were added to fill 40% of the Tank 2 size by volume (1514.2 L) which each have a diameter of 2.5 cm, height of 0.4 cm, and a projected surface area of 800 m 2 /m 3 (AnoxKaldness TM -Veolia Water Technologies AB). Tank 2 was inoculated with 18.9 L of a liquid suspension containing ammonia oxidizing microbial strains (MICROCAT-XNL, Bioscience, Inc., Allentown, PA, USA) for establishment of nitrifying bacteria populations. A third 3.76-m 3 tank, Tank 3, is for settling of any remaining solids. Paired side-by-side denitrifying woodchip bioreactors referred to as East Bioreactor (EB) and West Bioreactor (WB) are each 9.1 m by 12.2 m with a depth of 1.4 m. The woodchip bioreactors are continuously saturated, the bottom, side walls, and center barrier splitting the two bioreactors are lined with 1.1-mm ethylene propylene diene monomer (EPDM) and contain mixed hardwood woodchips (approximate size of individual chips: 5 cm × 5 cm × 0.6 cm). The species composition of the woodchips was 60% Ash (Fraxinus sp.), 20% Yellow Birch (Betula alleghaniensis), and 20% Silver Maple (Acer saccharinum), determined by a local woodchip supplier, and was based on available feedstock species at the time of construction. The bioreactors receive effluent that is split after leaving Tank 3 through two sets of pipes distributed across the top perimeter of each bioreactor. The effluent percolates vertically into the woodchips in each bioreactor. The bioreactors remain saturated and internal water depth is governed by an adjustable water level control structure. When the woodchip bioreactors discharge, their outlets flow through parallel 100-mm diameter underdrain pipes into an earthen infiltration basin. Sampling occurred at four locations: before Tank 1 (System Inflow), after Tank 3 (Tank Outflow, prior to entering the woodchip bioreactors), and at the outflows of both EB and WB ( Figure 2).
The system has three distinct runoff flow paths ( Figure 2), depending on the flow rate, which allowed for bypass of some treatment components during larger storm events to protect the system. The three paths are referred to as 'low flow', 'high flow', and 'extreme events' and are designed to accommodate flow up to a certain intensity ( Figure 2).

•
Low flow (leachate flows and storm intensities of up to 45.7 mm hr −1 ): behind the screen assembly, the lowest inlet pipe in the flow diversion structure directs leachate and runoff to enter the treatment tanks.

•
High flow (storm intensities up to 71.1 mm hr −1 ): runoff enters the high flow inlet in the flow diversion structure. In this case, runoff bypasses the treatment tanks and is sent directly to the surface of the woodchip bioreactors.

•
Extreme events (storm intensities greater than 71.1 mm hr −1 ): runoff enters a third inlet, bypasses the treatment tanks and the bioreactors, and sent directly into the infiltration basin.
For this study, the treatment system was monitored from June to October 2019. A tipping bucket rain gauge (Onset RG3 Hobo Rain Gauge Data Logger, Bourne, MA, USA), installed at the site, measured rainfall depths throughout the sampling period.
Routine maintenance was performed to prevent malfunctioning of the system. Weekly inspections were conducted during the growing season and subsequent maintenance involved removal of solids from the screen assembly and area around the flow diversion structure, removal of plant growth from the bioreactors, and addressing any potential flow blockages in the pipe network.

Storm Event Sampling
Collection of flow-based water samples occurred at four locations throughout the system ( Figure 2) using automated water samplers (Teledyne ISCO 6712, Lincoln, NE, USA). Autosamplers were equipped with water-level measurement modules, which were set to record water level on a one-minute interval. Water-level measurements were used to calculate flow, as described below. To prevent drift, modules were calibrated at the zero-level point at the start of and throughout the sampling season.
Inflow samples were taken from the inlet pipe to Tank 1 (System Inflow, Figure 2), where a compound weir (Thel-Mar, Boone, NC, USA) and a water level sensor (Teledyne ISCO 730 Bubbler Flow Module, Lincoln, NE, USA) measured stage. Flow rates were determined using stage-discharge values provided by the weir manufacturer. A second autosampler took samples from the effluent of Tank 3 (Tank Outflow, Figure 2). This autosampler was triggered via a communications cable to take samples concurrently with the first autosampler.
Bioreactor outflow samples were collected individually from both water level control structures (WB and EB, Figure 2). As the outflow from each bioreactor filled the sumps, effluent exited via a 15.2-cm diameter spillway equipped with an orifice plate containing a 7.6-cm diameter orifice, the top of which was the reference zero point for the water level sensor. Flow rate (q) was calculated from level readings from the attached water level sensor using a table of values developed from plotting H (autosampler level readings) vs. q; under weir (Equation (1)) and orifice (Equation (2)) conditions.
Each autosampler collected water samples on a predetermined flow-based interval. When compositing multiple samples across a storm hydrograph, flow-based sampling is preferred to limit sample processing [29]. Samples of 100 mL were taken after a specific volume of runoff passed through the sampling location. Each autosampler had four 3.7 L sampling bottles and a sampling interval was chosen to accommodate up to the anticipated runoff quantity for this drainage area for a 10-year 24-h storm in South Burlington, VT calculated using the Soil Conservation Service Curve Number equation [30]. Within 24 h of a storm event, samples were collected, transported to the laboratory, and composited in the laboratory to prepare for analysis.

Water-Quality Analysis
Samples were analyzed for concentrations of TN, combined nitrite + nitrate-N (NO 2 − +NO 3 − -N), NH 4 + -N, TP, and SRP. Samples for N species analysis were preserved with sulfuric acid, SRP samples were filtered through 0.45-µm nylon mesh filters and samples for TN and TP were prepared through persulfate digestion. Nutrient concentrations were analyzed using flow injection analysis instruments (Lachat QuickChem8000 AE, Hach Inc., Loveland, CO) using standard methods (listed in Supplementary Materials Table S1) [31]. System Inflow samples during the 2018 season were analyzed by the Vermont State Agriculture and Environmental testing laboratory. During the 2019 season, all samples were analyzed by the University of Vermont Agriculture and Environmental testing laboratory.
Any concentration values below the instrument's detection limit were substituted with one-half the detection limit [32,33].

Nutrient Mass Load Calculations and Removal Efficiency
For each sampling location, total nutrient mass loads were determined based on runoff volumes and measured event mean concentration (EMC) (Equation (3)). Volumes were calculated using flow rates recorded each minute by autosampler flow modules and summing the total volume over the sampling interval.
Total Mass = Q * EMC, where, Q = flow volume (L) and EMC = event mean concentration of analyte (mg L −1 ). Nutrient loads were used to determine the mass load removals and mass removal efficiency using the equation [32]: where, M i = initial mass entering the system (g), M f = final mass leaving the system (g). For positive values nutrients are retained, for negative values there is an export of nutrients. One storm event was removed from the analysis due to flow measurement instrumentation error.

Bioreactor Flow Distribution Challenges
For most storm events, a portion of the storm flow volume bypassed the tanks (i.e., entered the high flow path, Figure 2) and went directly to the woodchip bioreactors. This was likely due to issues including clogging of the inflow pipe, backup of flow through the tanks, or rapid influent flow rate during high-intensity storm events. It was determined that summed volumes from bioreactor outflow autosamplers, with a correction for added effective precipitation on the open bioreactors, more accurately accounted for storm event runoff volumes that passed through the entire system. Additionally, the uneven split of flow between the bioreactors, which may have occurred at the flow diversion structure or at the split in the distribution system just past the Tank Outflow, led to the WB receiving greater runoff volume than EB.

Statistical Analysis
Event mean concentrations and nutrient mass loads were evaluated between the four sampling locations, System Inflow, Tank Outflow, WB Outflow, and EB Outflow. Each monitored storm event was considered a replicate for statistical purposes [34]. These data violated both normality and equal variance assumptions, therefore a Wilcoxon signed rank test, a non-parametric analog to the paired t-test, was used to make comparisons between median influent and effluent nutrient concentrations of the treatment system components. For all storm events where nutrient mass loads were calculated, System Inflow and Combined Bioreactor Outflow nutrient mass loads were compared via Wilcoxon signed rank tests. Results were evaluated at 95% confidence (α < 0.05), with p values between 0.05 and 0.10 considered marginally significant. All models were run using R statistical software version 4.0 [35].

Silage Storage Bunker Runoff Nutrient Composition
To characterize the composition of silage bunker runoff, storm event samples were collected at the System Inflow from June through November 2018 and June through October 2019. Mean concentrations from flow-weighted samples for each nutrient are presented in Table 1.

Storm Events and Flow Rates
During the months of June through October 2019, a total of 16 storm events were monitored at all autosampler points in the system (Figure 3). The mean storm depth was 18.78 mm and the storm depths ranged from 5.33 mm to 52.58 mm. The average flow rate during storm events leaving the WB was 0.28 L s −1 and 0.17 L s −1 from the EB.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 7 of 17 and 0.10 considered marginally significant. All models were run using R statistical software version 4.0 [35].

Silage Storage Bunker Runoff Nutrient Composition
To characterize the composition of silage bunker runoff, storm event samples were collected at the System Inflow from June through November 2018 and June through October 2019. Mean concentrations from flow-weighted samples for each nutrient are presented in Table 1.

Storm Events and Flow Rates
During the months of June through October 2019, a total of 16 storm events were monitored at all autosampler points in the system (Figure 3). The mean storm depth was 18.78 mm and the storm depths ranged from 5.33 mm to 52.58 mm. The average flow rate during storm events leaving the WB was 0.28 L s −1 and 0.17 L s −1 from the EB.

Nitrogen
Total N concentrations were highest at the System Inflow, and generally decreased as the runoff was treated in the system (Figure 4). Ammonium-N concentrations increased within the treatment tanks (i.e., between System Inflow and Tank Outflow), decreased from the WB (i.e., between Tank Outflow and WB Outflow) and were generally similar or increased within the EB (i.e., between Tank Outflow and EB Outflow) (Figure 4). All means, sample sizes, and standard deviations are shown in Supplementary Materials Table S2

Nitrogen
Total N concentrations were highest at the System Inflow, and generally decreased as the runoff was treated in the system (Figure 4). Ammonium-N concentrations increased within the treatment tanks (i.e., between System Inflow and Tank Outflow), decreased from the WB (i.e., between Tank Outflow and WB Outflow) and were generally similar or increased within the EB (i.e., between Tank Outflow and EB Outflow) (Figure 4) 4.74-40.7) mg L −1 for NH4 + -N, and 0.025 (0.025-7.7) mg L −1 for NO2 -+NO3 --N. Compared to the System Inflow, NH4 + -N median concentration in Tank Outflow increased, and that difference was found to be statistically significant (p < 0.001).
The WB was most effective at reducing N concentrations. Median concentrations measured at the WB Outflow were 16.61 (range: 12.15-52.34) mg L −1 for TN, 5.14 (2.75-39.55) mg L −1 for NH4 + -N, and 0.025 (0.025-0.13) mg L -1 for NO2 -+NO3 --N. The median effluent concentrations of TN and NH4 + -N at the EB Outflow were both reduced relative to Tank Outflow (i.e., inflow to the WB) and those differences were found to be statistically significant (p < 0.001, p = 0.002, respectively).
The EB Outflow median TN concentrations were higher than the WB Outflow. At the EB Outflow, the median TN concentration was 39.74 (range: 17.7-92.34) mg L −1 , the NH4 + -N median concentration was 24.5 (range: 1.88-75.5) mg L −1 , and the NO2 -+NO3 --N median concentration was 0.025 (range: 0.025-0.089) mg L −1 . Differences in concentrations compared to Tank Outflow were not found to be significant.

Phosphorus
Overall, the treatment system was successful at reducing TP and SRP concentrations, but P reduction was most effective within the bioreactors rather than the tanks ( Figure 5 and Supplementary Materials Table S2 The WB Outflow median TP concentration was 12.1 (range: 9.72-19.5) mg L −1 and 10.27 (range: 7.8-18.16) mg L −1 for SRP. There was a reduction for both median TP and SRP concentrations compared to the Tank Outflow measurements, and those difference were found to be statistically significant (p = 0.002, p < 0.001, respectively). Median concentrations from the EB Outflow were 13.1 (range: 11.0-27.3) mg L −1 for TP and 11.58 (range: 6.75-21.6) mg L −1 for SRP, and there was not a statistically significant difference compared to Tank Outflow concentrations.

Phosphorus
Overall, the treatment system was successful at reducing TP and SRP concentrations, but P reduction was most effective within the bioreactors rather than the tanks ( Figure 5 and Supplementary Materials Table S2 The WB Outflow median TP concentration was 12.1 (range: 9.72-19.5) mg L −1 and 10.27 (range: 7.8-18.16) mg L −1 for SRP. There was a reduction for both median TP and SRP concentrations compared to the Tank Outflow measurements, and those difference were found to be statistically significant (p = 0.002, p < 0.001, respectively). Median concentrations from the EB Outflow were 13.1 (range: 11.0-27.3) mg L −1 for TP and 11.58 (range: 6.75-21.6) mg L −1 for SRP, and there was not a statistically significant difference compared to Tank Outflow concentrations.

Analysis of All Storms
Cumulative mass load removals and removal efficiencies (RE %) were calculated from summed System Inflow and Combined WB and EB Outflow mass loads for 15 storm events during the 2019 sampling period to assess overall system performance ( Table 2). The system removed 76% of the influent TN mass load. On a storm-by-storm basis comparing TN inflow and outflow loads, the difference was statistically significant (p = 0.002). The treatment system had a net export of NH4 + -N mass load over the season, and the difference between inflow and outflow loads was significant (p = 0.004). There was very little NO2 -+NO3 -N measured throughout the system, but overall there was a reduction in mass load and a statistically significant difference in inflow and outflow loads (p = 0.005). Moderate mass load reductions for TP and SRP occurred within the system, however those differences were not found to be statistically significant. Table 2. Cumulative mass removal of nitrogen (N) and phosphorus (P) mass loads and overall removal efficiency (RE %) in 2019 for 15 storm events where loads were calculated, between System Inflow and Combined Bioreactor Outflow (i.e., from East and West Bioreactor Outflows). Negative values indicate an export of nutrient. TN = total N, NH4 + -N = ammonium-N, NO2 -+NO3 --N = nitrite + nitrate-N, TP = total P, SRP = soluble reactive P. There were five storm events in 2019 for which tank bypass did not occur and the full runoff volume from each storm event was treated in the treatment tanks with MBBR prior to entering the woodchip bioreactors ( Figure 6). For all storms, TN mass load decreased between System Inflow and

Analysis of All Storms
Cumulative mass load removals and removal efficiencies (RE %) were calculated from summed System Inflow and Combined WB and EB Outflow mass loads for 15 storm events during the 2019 sampling period to assess overall system performance ( Table 2). The system removed 76% of the influent TN mass load. On a storm-by-storm basis comparing TN inflow and outflow loads, the difference was statistically significant (p = 0.002). The treatment system had a net export of NH 4 + -N mass load over the season, and the difference between inflow and outflow loads was significant (p = 0.004). There was very little NO 2 − +NO 3 − -N measured throughout the system, but overall there was a reduction in mass load and a statistically significant difference in inflow and outflow loads (p = 0.005). Moderate mass load reductions for TP and SRP occurred within the system, however those differences were not found to be statistically significant.

Nutrient Mass Loads for Storms with No Tank Bypass
There were five storm events in 2019 for which tank bypass did not occur and the full runoff volume from each storm event was treated in the treatment tanks with MBBR prior to entering the woodchip bioreactors ( Figure 6). For all storms, TN mass load decreased between System Inflow and Combined Bioreactor Outflow (i.e., sum of WB Outflow and EB Outflow), however there was generally an increase in TN mass load within the treatment Tanks ( Figure 6).
Combined Bioreactor Outflow (i.e., sum of WB Outflow and EB Outflow), however there was generally an increase in TN mass load within the treatment Tanks ( Figure 6). These storm events were analyzed separately (Table 3). Cumulative nutrient mass loads and average removal efficiencies were calculated between the System Inflow and Tank Outflow, Tank Outflow and WB, and Tank Outflow and EB. Nutrient mass loads from individual storm events are shown for these locations. The treatment tanks with MBBR on average for each storm event increased nutrient mass loads, and overall most mass removal occurred within the WB (Table 3). Table 3. Cumulative nutrient mass load removal and mean removal efficiency for each part of the system (i.e., treatment tanks with MBBR, West Bioreactor (WB), and East Bioreactor (EB)) for five storm events where no bypass of treatment tanks occurred.  These storm events were analyzed separately (Table 3). Cumulative nutrient mass loads and average removal efficiencies were calculated between the System Inflow and Tank Outflow, Tank Outflow and WB, and Tank Outflow and EB. Nutrient mass loads from individual storm events are shown for these locations. The treatment tanks with MBBR on average for each storm event increased nutrient mass loads, and overall most mass removal occurred within the WB (Table 3). Table 3. Cumulative nutrient mass load removal and mean removal efficiency for each part of the system (i.e., treatment tanks with MBBR, West Bioreactor (WB), and East Bioreactor (EB)) for five storm events where no bypass of treatment tanks occurred.

Silage Storage Bunker Runoff Composition
This study contributes to the understanding of nutrient composition of silage bunker runoff, which is rarely reported in the literature. Leachate composition is highly variable depending on time that elapsed since ensiling and the initial crop moisture content, which in turn leads to variation in runoff composition once it is diluted with precipitation [4]. Total N is comprised of dissolved species, NH 4 + -N and NO 2 − +NO 3 − -N and dissolved and particulate organic N. Studies reviewed by Gebrehanna [3] suggest that NO 2 − +NO 3 − -N may only make up 0.02% of the N content of silage leachate with most in the organic form, which is consistent with silage bunker runoff concentrations in this study (Table 1). Faulkner [6] reported silage bunker runoff EMC concentrations for NH 4 + -N of 59 mg L −1 and NO 3 -N of 4 mg L −1 , both higher than averages reported in this study (Table 1). For P, SRP contributes over 75% of the TP concentration on average (Table 1). Mean P concentrations were less than concentrations measured by Holly [5], for which TP ranged from 26 to 71 mg L −1 and SRP from 20 to 67 mg L −1 and the EMC of SRP from Faulkner [6] was 37 mg L −1 .

Treatment Tanks with Moving-Bed Biofilm Bioreactor (MBBR)
This treatment system was designed to facilitate conditions for N removal via microbial transformations. The MBBR system in Tank 2 is the first step in this process, as it is designed for conditions for nitrification, where biofilm reactors and aeration provide oxygen and nitrifying microbial populations, for the conversion of NH 4 + to NO 3 − . The treatment tanks were evaluated as an entire system since sampling occurred at the inflow to Tank 1 and outflow from Tank 3. A slight decrease in the average TN concentration and a significant increase in average NH 4 + -N concentrations from the System Inflow to Tank Outflow, suggest that ammonification, conversion of organic N to NH 4 + , likely occurred in the tanks ( Figure 4). The average NO 2 − +NO 3 − -N concentration entering the bioreactors from the tanks was 1.66 mg L −1 , however it is notable that for 11 of 15 monitored storm events, the NO 2 − +NO 3 − -N concentration at the Tank Outflow was below 0.05 mg L −1 (the analyzing instrument's detection limit). In 2019, there were four consecutively monitored storms in mid-summer where Tank Outflow sample concentrations were between 3 and 7.7 mg L −1 ; an increase from System Inflow, suggesting there might have been a time where nitrification successfully occurred in the treatment tanks. It is likely that nitrification performance was inhibited by an inadequate hydraulic retention time during storm events, which is further discussed below.
It is important to note that monitoring of this system occurred in conjunction with storm events, when each component of the system overflows to the next component and autosampling is triggered by flow. In other words, the system is not monitored at times between storms. Generally, nutrient concentrations and mass loads were higher in Tank Outflow samples leaving the treatment tanks, compared to the System Inflow (Figures 4 and 6, Table 3), which may be attributed to conditions that occur in between storm events when the system was not monitored. This could potentially come from unmonitored silage leachate (i.e., not mixed with storm-driven runoff), which flows into Tank 1. Leachate concentration would be diluted within the treatment tanks with MBBR and also contribute to an increase in concentration measured in Tank Outflow during the following storm event.

Bioreactors
Cumulatively over the sampling season, TN mass load was reduced by 76% (Table 2) Table S2). This overall decrease indicates that N removal via denitrification occurred within the beds, preceded by nitrification in the aerobic zones of the bioreactor, which has been described as a possible removal mechanism of NH 4 + -N in in other woodchip bioreactor studies [21,25]. However, there was likely mineralization of organic N to NH 4 + -N, due to the increase in median NH 4 + -N concentration between the Tank Outflow and WB and EB Outflows, which was not consistently followed by nitrification and denitrification processes.
Most woodchip bioreactor studies sample influent and effluents on regular intervals when flow rates are relatively constant, however treatment of storm flows within bioreactors is thought to be less effective due to disturbance of steady-state conditions [36]. A meta-analysis of bioreactor studies found that performance was significantly lower with a hydraulic retention time (HRT) of less than 6 h compared to longer times between 6 and 20 h and >20 h [20]. For this system, the steady-state HRT within each bioreactor during storm events ranged from 29 to 309 hr and HRT between storm events varied from 1 day to 12 days. This variation in HRT and flow rates through the system may have inhibited performance by disrupting stable microbial populations, especially at rapid flow rates. Future work could modify outflow structures to control flow rates and HRT.
Bioreactor performance was limited by low NO 2 − +NO 3 − -N concentrations entering the woodchip bioreactors from the tanks (Figure 4), which has also been reported in other studies [27,37]. A consequence of low NO 2 − +NO 3 − -N concentrations in the bioreactors can be undesirable redox reactions such as sulfate (SO 4 2− ) reduction, which produces hydrogen sulfide (H 2 S). There is the potential for these reactions to occur in woodchip bioreactors if complete NO 2 −− +NO 3 − -N removal is achieved and the NO 2 − +NO 3 − -N concentration falls below 0.5-1 mg L −1 [22]. However, SO 4 2− concentrations were not measured in this study, so future work could investigate if sulfate reduction occurs in this system due to low NO 2 −− +NO 3 − -N concentrations that were measured (below 0.05 mg L −1 to 0.13 mg L − ). The release of nitrous oxide (N 2 O) due to incomplete denitrification is a further 'pollution swapping' concern, which should be investigated for this system. Generally, bioreactor studies have found low levels of N 2 O production [21,38], but more work on greenhouse gas emissions from bioreactors is needed [20].

Treatment Tanks with MBBR
Throughout the sampling season, moderate removal occurred in the entire system for both TP and SRP, as measured in both mass loads and concentration reductions (Table 2 and Figure 5). Analysis of the subset of storms without tank bypass storms suggest that most of the load removal occurred in the bioreactors, as opposed to the treatment tanks with MBBR (Table 3 and Figure 6). Phosphorus concentrations and loads generally increased from System Inflow to Tank Outflow locations. This supports the hypothesis that undiluted silage leachate enters the tanks between monitored storm events, increasing nutrient concentrations in Tank 1, and subsequently the concentrations entering the bioreactors from Tank Outflow.

Woodchip Bioreactors
Woodchip bioreactors are typically designed to target N, so P removal reported in woodchip bioreactor studies is highly variable. Some studies found moderate removal [21], while others observed a net P export [39,40]. The exact removal mechanism for TP and SRP in the woodchip bioreactors ( Figure 5, Tables 2 and 3) in this study is uncertain, but possibilities include sorption to woodchips and microbial immobilization of P. High removal of TP concentrations was observed in a woodchip bioreactor treatment system, which reduced influent concentrations by 71%, but most reduction was observed from a settling tank [41]. Robertson [42] found a high TP removal rate in pilot-scale woodchip biofilters with high influent TP concentrations from stream flow. Removal rates increased at higher loading rates, a promising finding for storm event flow. Soupir [40] reported that non-sterile woodchips, i.e., with microbial growth, removed more dissolved P than sterile woodchips. Incorporating a P-binding additive into the bioreactor, such as drinking water treatment plant residuals (DWTR) has been successful in removing both TP and SRP via sorption onto these materials, however removal efficiencies were greater for SRP [43]. The aforementioned study also compared pilot scale woodchip bioreactors amended with DWTR and without, and for the non-amended woodchip bioreactors authors saw an average of 19% (TP) and 12% (SRP) removal efficiency, which is comparable to the average RE observed in the WB from this study ( Table 3). Removal of P in these non-DWTR bioreactors was thought to have been due to physical filtration [43]. A possibility to further enhance removal of P in this system would be to incorporate a P-binding additive into one woodchip bioreactor, which could be compared with this paired system.
Notably, System Inflow P concentrations (Table 1) were much higher than other woodchip bioreactor studies mentioned above. The treatment system from Choudhury [41] received vegetable wash water with observed TP concentrations from 5 to 10 mg L −1 . Dissolved P concentrations were 1.12 mg L −1 for the woodchip sorption experiment from Soupir [40]. In the lab study from Gottschall [43] influent concentrations were 0.6 mg L −1 and 0.2 mg L −1 for TP and dissolved P, respectively. Higher influent concentrations in this system may have an impact on P sorption dynamics, however further consideration of these processes is needed.

Bioreactor Performance: East Bioreactor (EB) Versus West Bioreactor (WB)
Superior nutrient removal performance was observed from the WB compared to EB for both concentration reductions and mass load removals (Figures 4-6 and Table 3). As designed, the bioreactors were intended to be replicates of one another; however, observations during the study period made it clear that as-built features resulted in unequal hydraulic loading. Notably, greater runoff volumes flowed through WB compared to EB for 14 out of the 15 storms analyzed for nutrient mass loads, and on average WB received 28% more runoff during each event than EB. Since the WB outflow sampling structure had larger volumes of water passing through it compared to the EB outflow, a potential explanation is that increased water mounding occurred in the WB than EB, before drainage could occur. This additional mounding would have resulted in a greater volume of the top woodchip layers fluctuating between saturated and unsaturated conditions.
The fluctuating water level in the WB is comparable to woodchip bioreactor studies on drying-rewetting cycles (DRW), where a bioreactor is drained and air is allowed to fill the once saturated spaces between the woodchips, overall increasing nitrate removal [44,45]. These conditions in soil systems are known to enhance decomposition of organic matter, making C and N available, which can accelerate rates of soil metabolic processes [46][47][48]. Results from a DRW event was first reported in a woodchip bioreactor column study from Christianson [49]. Researchers saw increased nitrate removal following a 96-hr drying time and the capacity for dissolved P removal in the woodchips also increased. It is thought that DRW cycles cause decomposition of the woodchips to occur in the presence of oxygen, which increases dissolved organic C (DOC) and can then be used by microbes in denitrification. In this study where more woodchips in WB are saturated during storm flow and then subsequently dry out, the DRW may explain a difference in performance, leading to more DOC available for rapid denitrification. For future studies, modification of the outlet structures from each bed could be used to control outflow and retention time to investigate whether more water mounding in the WB does indeed occur. While woodchips have been shown to be an effective C source for NO 3 − removal for at least 15 years of continuous bioreactor operation [50], decomposition of woodchips is a concern especially if DRW cycles may enhance this process.

Conclusions
The design of this treatment system is novel, as it is the first reported combined MBBR system and woodchip bioreactor study treating silage bunker runoff. As a whole, the treatment system successfully reduced TN mass load by 76% throughout the study duration. Most reduction of TN took place in the bioreactors, compared to the tanks, despite low nitrate concentrations entering the bioreactors. DRW may have been the primary driver for N transformations and removal. While the system's design specifically targeted treatment of N species, moderate P removal was also observed; 26% of the TP mass load and 19% of SRP mass load over the sampling season. Overall, the system was an effective treatment option for high concentrations of N and P in the influent silage bunker runoff evaluated during the third year of operation. Other treatment options or treatment tank designs should be considered prior to runoff entering the bioreactor to maximize the transformation of N into NO 3 − . To increase hydraulic retention time and nitrification performance of the treatment tanks, future designs could consider a decreased flow rate through the system that incorporates substantial storage up-stream of the treatment tanks. Simple hydraulic control (e.g., an orifice) could then be used to slowly release this stored runoff through the treatment system in between storm events at a rate that would allow for full nitrification to occur. Such modifications would likely improve performance but could also increase installation costs. Longevity of the woodchips should be monitored on an annual basis to ensure a readily available source of C for denitrification. Issues of unequal flow splitting and extra unintended bypass to the high flow path around treatment tanks highlight the importance of proper operation and maintenance of silage bunker runoff treatment systems. Areas for solids settling that can easily and frequently be cleaned out are critical to effective performance [6]. In other water-quality contexts, where N is already in the form of NO 3 − , this combined MBBR and woodchip bioreactor system may also offer a promising solution.