Pairing High-Frequency Data with a Link-Node Model to Manage

SUPPLEMENTAL MATERIALS FOR:

Pairing high-frequency data with a link-node model to manage

dissolved oxygen impairment in a dredged estuary

by

Mary Kay Camarilloa*, Gregory A. Weissmanna, Shelly Gulatia,

Joel Herrb, Scott Sheederb, and William T. Stringfellowa,c

aEcological Engineering Research Program, School of Engineering & Computer Science, University of the Pacific, 3601 Pacific Avenue, Stockton, CA 95211, USA.

bSystech Water Resources, Inc., 1200 Mount Diablo Boulevard, Suite 102, Walnut Creek, CA 94596, USA.

cEarth Sciences Division, Lawrence Berkeley National Laboratory, 1 Cyclotron Road, Berkeley, CA 94720, USA.

*Corresponding author: Email: , Phone: +1-209-946-3056, FAX: +1-209-946-3086

Prepared for submission to:

Environmental Monitoring and Assessment

S1.  Model Formulation

The SJR-Link-Node model has been previously described (Chen and Tsai 2002a, 2002b; Sheeder and Herr 2013), but is summarized here for clarity. The model was not developed as part of this project. Rather, an established model was used that had previously undergone peer-review and had been tested in the SJR Estuary.

The SJR-Link-Node model simulates flow between the nodes as described by Eqn. 1 where U is velocity (ft s-1), t is time (s), X is horizontal distance (ft), g is gravity acceleration (ft s-2), H is water surface elevation (ft), and n is Manning’s friction factor (unitless). In Eqn. 1, flow acceleration is a function of momentum, head differential, and friction loss.

dUdt=-UdUdX-gdHdX-nUU / (1)

Once integrated, Eqn. 1 is transformed into Eqn. 2 that is used for numerical simulation.

Ut=Ut-1-UdUdX+gdHdX+nUt-1Ut-1 / (2)

Once the velocity, U(t), is calculated using Eqn. 2, it is then used to determine the flow rate, Q (ft3 s-1), as shown in Eqn. 3, where Xa is the cross sectional area (ft2) of the channels (links).

Q=UXa / (3)

Once flows for the individual links are determined, they can be summed at the nodes and used to determine the resulting water surface elevation, H(t), as shown in Eqn. 4, where Sa is the surface area (ft2) of each node.

Ht=Ht-1-QSa∆t / (4)

Calculation of Eqn. 1 – 4 yields model output consisting of velocity and flow rate for the links and water surface elevation for the nodes that is then used to calculate node volumes. Once the flow rate is established within the model, it is used to calculate the concentrations of water quality parameters using a mass balance approach (Eqn. 5) where V is the water volume (ft3) at the nodes, C is the concentration (mg L-1), and α is a diffusion coefficient (ft2 s-1).

dVCdt=QC+αdCdX-AND-Sinks+Sources / (5)

The term AND is a dimensionless anti-numerical dispersion term (discussed below). Sinks are losses of individual water quality parameters resulting from the decay, uptake, and diversion of the individual water quality parameters, while sources are gains from waste discharges, chemical transformations, and biological growth. Most of the biological and chemical processes within the SJR-Link-Node model—contributing to the sources and sinks of DO—are represented by first-order reactions (Table S1 in Section S2).

As shown in Eqn. 5, the volume of the individual nodes, V, and concentration of water quality parameters, C, can change over time. Eqn. 5 can be transformed into Eqn. 6 for calculation by the model. Values for V and Q, necessary to calculate Eqn. 6, are obtained using Eqn. 3 and 4.

VdCdt=-CdVdt+QC+αdCdX-AND-Sinks+Sources / (6)

The output of Eqn. 6 is the changing rate of concentration (dC dt-1) that is then used to calculate C(t) as shown in Eqn. 7. In Eqn. 6, the upstream concentrations are used to calculate mass flow in the links.

Ct=Ct-1+dCdt∆t / (7)

The anti-dispersion term, AND, present in Eqn. 5 and 6, was introduced to prevent masses from advancing too quickly from one node to another, which can be problematic in link-node models (Chen and Tsai 2002a). The model approach is based on the UPWIND scheme for numerical dispersion as described by Roache (1972). The term AND is calculated using Eqn. 8 and 9.

AND=12u∆x1-c / (8)
c=u∆t∆x / (9)

Eqn. 6 is used for each water quality parameter included within the model where each parameter has different sinks and sources. Water quality parameters simulated in SJR-Link-Node are those affecting DO: carbonaceous biochemical oxygen demand (CBOD, mg L-1), DO (mg L-1), ammonia-N (mg L-1), nitrate-N (mg L-1), coliform (MPN (100 mL)-1), total dissolved solids (mg L-1), phosphate-P (mg L-1), chlorophyll-A (μg L-1), pheophytin-a (μg L-1), volatile suspended solids (VSS, mg L-1), total suspended solids (TSS, mg L-1), and sediment (mg L-1). For DO, the sinks are CBOD decay, nitrification, sediment oxygen demand, algal respiration, and decay of detritus (represented by VSS). In addition to mass input from adjacent waters (point sources and tributary inflows), DO sources include algal photosynthesis and reaeration from the atmosphere. Loss of DO to the atmosphere—which occurs when the water is supersaturated with DO—can also occur. The aeration rate was calculated using the empirical equation of O’Connor and Dobbins (O’Connor and Dobbins 1958) with an added term for the effects of wind, as shown in Eqn. 10, where Ka is the aeration coefficient (d-1), θ is a coefficient used to correct for temperature effects (unitless), W is the wind speed (m s-1), and D is the water depth (ft).

Ka=12.9U12D32θT-20+0.15W2D / (10)

Using the aeration coefficient, Ka, calculated in Eqn. 10, the flux of oxygen across the water surface, (F, in lb d-1) is calculated as shown in Eqn. 11 where DOs is the dissolved oxygen solubility at temperature T (°F). An additional term, a linear aeration adjustment factor, is used to control the atmospheric oxygen flux, F.

F=KaADOsT-DO / (11)

Algae growth is a function of light intensity in the water column, which is calculated using Beer’s Law, per Eqn. 12, where IZ (W m-2) is the light intensity at depth Z (ft), I is the surface light intensity (W m-2), and λ is the light extinction coefficient (ft-1). Because SJR-Link-Node is a 1-D model, light intensity is averaged over the entire depth of the water column. Since the photic zone is only present in the top 0.5 – 1.0 m of the water column and the channel depth is 10.7 m, the result is that average light intensity in the DWSC is minimal. To account for photosynthesis in the top of the water column where light is available, phytoplankton growth was simulated separately in the top 0.6 m of the water column.

IZ=Iexp-λZ / (12)

Since light intensity is a function of depth, the growth rate of phytoplankton is also depth-dependent, as shown in Eqn. 13, where GZ is the growth rate at depth Z and K is the half-saturation coefficient of light (W m-2).

GZ=IZK+IZ / (13)

Integrating Eqn. 13, the depth-averaged, light-limited algal growth rate (Gd) can be calculated using Eqn. 14.

Gd=1λDlnK+IK+Ie-λd / (14)

The model includes four types of settleable solids: pheophytin (dead phytoplankton), detritus (measured as VSS), inorganic solids, and sand. These settled solids can be scoured from the bottom surface and re-suspended. Scouring occurs when the flow velocity exceeds a critical value, VCR (m s-1) as shown in Eqn. 15. Once the critical scouring velocity has been exceeded, the scouring rate (S, kg s-1) can be calculated using Eqn. 16 where K is a calibration parameter, AW is the area of the wetted channel bed, VAVE is the average flow velocity (m s-1), and b is a model fitting parameter (unitless).

VCR=25∙0.65gd0.8D0.2 / (15)
S=KAWVAVE-VCRb / (16)

Heavier sand undergoes bed load transport, which is a function of the shear velocity V* (m s-1) as shown in Eqn. 17. The shear velocity can be used to calculate shear stress Y (unitless), as shown in Eqn. 18. The critical shear stress can determined from the Shield Diagram (Graf 1971), as shown in Eqn. 19, which is a function of Reynold’s Number NR (unitless), as shown in Eqn. 20. The sand density, defined by the specific gravity γ (unitless), also affects sediment transport.

V*=gDS / (17)
Y=V*2γ-1gd / (18)
YCR=fNR / (19)
NR=V∙dν / (20)

The Yalin equation, as shown in Eqn. 21 and 22, is used to calculate bed load transport of sand where ρw is the density of water (kg m-3), and W is the wetted perimeter of the channel (m). The eroded sand in excess of Tf is re-deposited in the riverbed. Only the excess remains in suspension.

Tf=PSγρWdV∙W / (21)
Ps=0.635YYCR-11-ln1+2.45ρ-0.40∆YCR2.45ρ-0.40∆YCR / (22)

The SJR-Link-Node model simulates natural spring and neap tide cycles. The model is configured to allow water and mass to flow from one node to another via the links, simulating flux and exchange resulting from tidal action. The effect of the tidal exchange is described in Eqn. 23 where E (unitless) is an exchange coefficient. A tidal exchange of 3% was used for every hour outside of the tidal boundary and approximately 10% over a tidal cycle.

Ci=Ci+1E+1-ECo / (23)

S2.  Model Coefficients

The parameters used in formulating the SJR-Link-Node model are shown in Table S1:

Table S1. San Joaquin River Link-Node model parameters.

Parameter / Value / Units /
Diffusion coefficient / 600 / ft2 s-1
Advection weighting factor / 1.0 / unitless
Exchange coefficient / 0.97 / unitless
Atmospheric turbidity factor / 1.5 / unitless
Evaporation coefficient A / 0.0 / unitless
Evaporation coefficient B / 3.60E-9 / unitless
Bank shading factor / 0.0 / unitless
Soluble CBOD decay coefficient / 0.30 / d-1
Ultimate CBOD CBOD5-1 ratio / 2.54 / mg mg-1
Aeration adjustment factor / 1.80 / unitless
Nitrification rate coefficient / 0.10 / d-1
Oxygen demand of ammonium / 4.57 / mg mg-1
Ammonia preference by algae factor / 0.60 / unitless
Sediment ammonia release rate / 0.000864 / g ft-2 d-1
Light extinction rate / 0.45 / ft-1
Algae self-shading factor / 0.236 / L ft-1 mg-1
Maximum growth rate of algae / 3.2 / d-1
Mortality and predation rate of algae / 0.05 / d-1
Respiration rate of algae / 0.15 / d-1
Maximum settling rate of algae / 2.46 / ft d-1
Nutrient effect on algae settling / 0.8 / unitless
Algae half-saturation constant for light / 0.0004 / kcal ft-1 s-1
Algae half-saturation constant for N / 0.01 / mg L-1
Algae half-saturation constant for P / 0.001 / mg L-1
Algae photo oxygenation demand / 1.6 / mg DO (mg algae)-1
Algae respiration oxygen demand / 1.6 / mg DO (mg algae)-1
N content of algae / 0.15 / mg N (mg algae)-1
Chlorophyll-a content of algae / 0.04 / mg (mg algae)-1
Pheophytin-a decay rate / 0.1 / d-1
Detritus decay rate / 0.04 / d-1
N content of detritus / 0.15 / mg N (mg VSS) -1
Oxygen demand of detritus / 1.6 / mg DO (mg VSS) -1
P content of detritus / 0.009 / mg P (mg VSS) -1
Coliform decay rate / 1.0 / d-1
SOD / g ft2 d-1
Theta values for temperature correction
BOD decay rate / 1.04 / unitless
SOD decay rate / 1.04 / unitless
Aeration rate / 1.02 / unitless
Nitrification rate / 1.08 / unitless
Sediment ammonium release rate / 1.02 / unitless
Algae growth rate / 1.03 / unitless
Algae mortality rate / 1.07 / unitless
Algae respiration rate / 1.03 / unitless
Pheophytin-a decay rate / 1.02 / unitless
Detritus decay rate / 1.02 / unitless
Coliform decay rate / 1.05 / Unitless

aBOD: biochemical oxygen demand; CBOD: carbonaceous BOD; CBOD5: 5-day CBOD; DO: dissolved oxygen; N: nitrogen; P:phosphorus; SOD: sediment oxygen demand

S3.  ANOVA Analyses

Tables S2 and S3 contain the results of ANOVA analysis performed (using JMP software) on model DO output (relative to the model baseline) and flow data from the Garwood Bridge Station.

Table S2. Results of a one-way ANOVA with Tukey HSD for San Joaquin River Link-Node model scenariosa,b.

"DWSC" scenario / "SJR" scenario / "RWCF" scenario / "Tribs" scenario
Month / Connecting letter / Mean DO improve-ment (mg/L) / Month / Connecting letter / Mean DO improve-ment (mg/L) / Month / Connecting letter / Mean DO improve-ment (mg/L) / Month / Connecting letter / Mean DO improve-ment (mg/L)
Nov / A / 0.584 / May / A / 1.175 / Sep / A / 0.17 / Mar / A / 0.165
Dec / A / 0.583 / Mar / B / 0.818 / Aug / A / 0.164 / Feb / B / 0.14
Oct / B / 0.556 / Feb / B / 0.8 / Jan / B / 0.139 / Apr / C / 0.104
Sep / B / 0.545 / Apr / C / 0.763 / Dec / C / 0.123 / Sep / D / 0.094
Jan / C / 0.515 / Oct / C / 0.741 / Feb / C / 0.123 / Jan / D / 0.094
Feb / C / 0.51 / Jun / D / 0.69 / Jul / D / 0.106 / Aug / E / 0.089
May / C / 0.509 / Nov / E / 0.631 / Oct / E / 0.098 / Jul / F / 0.073
Aug / C / 0.508 / Jan / F / 0.568 / Nov / F / 0.088 / May / G / 0.06
Mar / D / 0.455 / Jul / G / 0.497 / Jun / G / 0.064 / Nov / G, H / 0.058
Jul / E / 0.435 / Dec / H / 0.446 / Mar / G / 0.064 / Jun / H / 0.055
Jun / E / 0.434 / Aug / I / 0.415 / May / H / 0.047 / Oct / I / 0.043
Apr / F / 0.389 / Sep / I / 0.409 / Apr / H / 0.044 / Dec / I / 0.042

aModel calibrated using data from 2005-2010.

bSimulations based on: 1) Natural channel depth instead of dredged depth (“DWSC”), 2) Elimination of upstream mass loads (“SJR”), 3) Elimination of the Regional Wastewater Control Facility mass loads (“RWCF”), and 4) Elimination of mass loads from urban tributaries (“Tribs”).

Table S3. Results of a one-way ANOVA with Tukey HSD for observed net flow in the San Joaquin River, 2005 – 2010.

Month / Connecting letter / Mean net flow (m3 s-1)
May / A / 92.3
Apr / A / 90.1
Jun / B / 63.8
Mar / C / 48.6
Feb / D / 31.6
Jan / D,E / 29.2
Oct / E / 26.3
Dec / F / 21.2
Jul / G / 16.3
Sep / G, H / 14.3
Nov / H, I / 10.3
Aug / I / 7.8

S4.  References

Chen, C. W., & Tsai, W. (2002a). Final Report: Improvements and Calibrations of Lower San Joaquin River DO Model, Revised. CALFED 2000 Grant, CALFED 99-B16, DWR 4600000989. San Ramon, CA: Systech Engineering, Inc.