Dissolved oxygen and inorganic carbon dynamics in a high-energy coastal environment near Victoria BC’s untreated municipal sewage outfalls

Dissolved oxygen and inorganic carbon dynamics in a high-energy coastal

environment near Victoria BC’s untreated municipal sewage outfalls

By Jeremy Krogh B.Sc., University of Victoria, 2014 A Thesis Submitted in Partial Fulfillment of the Requirements for the Degree of MASTER OF SCIENCE in the School of Earth and Ocean Sciences © Jeremy Krogh, 2017 University of Victoria All rights reserved. This thesis may not be reproduced in whole or in part, By photocopy or other means, without the permission of the author.

Supervisory Committee

Dissolved oxygen and inorganic carbon dynamics in a high-energy coastal

environment near Victoria BC’s untreated municipal sewage outfalls

By Jeremy Krogh B.Sc., University of Victoria, 2014 Supervisory Committee Dr. Roberta Hamme, Co-Supervisor School of Earth and Ocean Sciences Dr. Debby Ianson, Co-Supervisor School of Earth and Ocean Sciences Dr. Kenneth Denman Member School of Earth and Ocean Sciences Christopher Lowe

Abstract

Wastewater disposal often has deleterious impacts on the receiving environment. Low levels of dissolved oxygen are of particular concern. Here I investigate the impacts on dissolved oxygen and carbon chemistry of screened municipal wastewater in the marine waters off Victoria, B.C., Canada. I analyzed data from a series of undersea moorings, ship-based monitoring, and underwater remotely-operated vehicle video. I used these observations to construct a two-layer box model of the nearfield receiving environment. Despite the lack of more advanced treatment, dissolved oxygen levels near the outfalls are well above the commonly used 63 µmol kg-1 _{hypoxic threshold and} that the outfalls’ impact on water column oxygen is likely less than a few µmol kg-1_. Likewise, dissolved inorganic carbon is not elevated and pH not depressed compared to the surrounding region. Strong tidal currents and the cold, well-oxygenated waters of Victoria’s marine environment give these waters a high assimilative capacity for organic waste.

Table of Contents

Supervisory Committee ... ii Abstract ... iii Table of Contents ... iv List of Tables ... vi List of Figures ... viii Acknowledgments ... xii Key Words ... 1 Highlights ... 2 1 Introduction ... 3 1.1 Study Site ... 4 1.2 Regional Physical Oceanography ... 5 2 Methods ... 9 2.1 Observational ... 9 2.1.1 Vertical Profiles ... 9 2.1.2 Inorganic Carbon and Total Alkalinity ... 9 2.1.3 Mooring Network ... 10 2.1.4 Ferry Box Underway System ... 11 2.1.5 Remotely Operated Underwater Vehicle (ROV) ... 11 2.2 Data Analysis methods ... 12 2.2.1 Time Series Analysis ... 12 2.2.2 Water Mass Endmembers ... 13 2.3 Model Methods ... 13 2.3.1 Physical Model ... 14 2.3.2 Oxygen Depletion ... 20 3 Results and Discussion ... 25 3.1 Observational ... 25 3.1.1 Vertical Profiles ... 25 3.1.2 Water Mass Analysis ... 32 3.1.3 Seabed conditions at the outfalls ... 39 3.1.4 Benthic community ... 39 3.1.5 Inorganic Carbon and Total Alkalinity ... 40 3.2 Model ... 41 4 Summary and Conclusions ... 46 Bibliography ... 48 Appendices ... 54 Appendix A: Macaulay Point Seabed Conditions ... 54 Appendix B: ROV Dive ... 60

Appendix E: Vertical Profiles ... 73 Oxygen ... 73 Salinity ... 74 Temperature ... 75 Sewage Plume Detection ... 76 Location of Outfall Monitoring Stations ... 79 Appendix F: Moorings ... 80 Locations ... 80 Deployments ... 80 Mooring Calibrations ... 82 Appendix G: Nutrient Data ... 83 Appendix H: Summer profiles from Haro Strait and JdF ... 84

List of Tables

Table 1: Modeled quantities. ... 18 Table 3: Biological and anthropogenic forcings. ... 19 Table 4: Prescribed values of state variables at the model’s boundaries. ... 20 Table 5: Summer (day of year 138 - 225) mean and standard deviations for the two endmember water masses. * Ferry data from the summer of 2015 was only available for July. ... 34 Table 6: The mixing ratios as predicted based on the source water properties of table 4 and the conditions observed at the outfall mooring. Due to non-complete data sets the Haro Strait observations from the summer of 2015 are applied to 2013 and 2014. Likewise, the JdF-W 2016 data are taken to be the same as the 2015 observations. Only 2015 has complete data from both the Haro Strait Surface and JdF-W. ... 35 Table B1: Species count at Macaulay Point ROV dive. ... 61 Table B2: Species Count at Clover Point ROV dive. ... 62 Table E1: Macaulay Point monitoring station locations. ... 79 Table E2: Clover Point monitoring station locations. ... 79 Table F1: Location and depth of sea bottom moorings. ... 80 Table F2: The deployment/recovery dates and serial numbers of the instruments deployed at the Macaulay outfall mooring. ... 80 Table F3: The deployment/recovery dates and serial numbers of the instruments deployed at the JdF-East mooring. ... 80 Table F4: The deployment/recovery dates and serial numbers of the instruments deployed at the Boundary Passage moorings. ... 81 Table F5: The deployment/recovery dates and serial numbers of the instruments deployed at the JdF mooring. ... 81 Table F6: The deployment/recovery dates and serial numbers of the instruments

Table F7: The calibration dates and drift amounts for the oxygen and CTD sensors used as part of the study. Drift values less than 1 indicate a decrease in sensitivity, 0.9799 indicates a decline of 2.01%. ... 82 Table H1: Location of select IOS/DFO stations, all of which are located along or near the thalweg. ... 84 Table H2: Summer (June and September cruises) conditions at station 56 in Boundary Passage between 2002 and 2015. ... 86 Table H3: Summer (June and September cruises) conditions at station 102 at the western entrance to JdF between 2002 and 2015. ... 86

List of Figures

Figure 1: (a) Map of the region, the colored stars show the locations of sea bottom moorings, the hashed area is a zone of high surface nutrients where primary production is light limited year round (Mackas & Harrison 1997). (b) Close-up map of the sewage outfalls with the black dots showing the location of monitoring stations where DIC/TA samples were collected and vertical profiles taken. (c) A cross-section of the average flow and approximate depth along the thalweg. The outfall mooring is not located along the thalweg, but is still attached to the seabed, north of the thalweg nearer to shore. ... 7 Figure 2: (a) Cartoon of the sewage plume and processes affecting T, S, O2 and DIC the outfall within the physical structure of the mass balance model. Red arrows show processes that consume oxygen and green arrows show process that produce DIC. Modelled quantities are listed along the left side of each physical layer and in Table 1. In the model, these quantities are perfectly mixed within their respective layers. Daily averaged net ocean currents are indicated for the long direction, cross direction currents are not shown but range from 0.9 – 3.5 km day-1 (25th and 75th percentiles). (b) Schematic of physical and biological modelled fluxes in the lower layer of the model showing the flow of O2and consumption of organic matter. Solid arrows represent volume fluxes that affect the O2concentration while dotted arrows represent biological fluxes (consumption). Temperature and Salinity (not shown) only experience physical fluxes, i.e.; advection, vertical mixing and dilution by sewage.. ... 17 Figure 3: Time series of dissolved oxygen, temperature, and salinity from the four-year mooring dataset filtered with a 25-hour running mean to remove semi-diurnal and diurnal tidal effects. The horizontal dashed magenta line indicates a hypoxic threshold of 63 µmol kg-1_{. All moorings are anchored to the bottom with} instruments approximately 6 m off the seabed, with exact depths of each mooring shown in figure 1b and listed in Table F1. ... 26 Figure 4: The left column shows histograms (probability density functions) of hourly oxygen observations at the outfall mooring, while the right column shows the same information but expressed as percent saturation on the bottom x-axis and partial pressure of O2 on the upper x-axis. The means and standard deviations are in units of µmol kg-1 on the left column and percent saturation on the right column. The fours seasons are spring (a,b, day of year 77-137), summer (c,d, 138-225), fall (e,f, 256-314), and winter (g,h 315-76) (Bylhouwer et al., 2013). ... 30

Figure 5 (a) Hourly salinity and oxygen at the Macaulay Point mooring during July 1-15 2015 (b) black circles are the same data as (a) with a linear fit while gray circles are all data from the summer (day of year 138-225) of 2015. (c) black circles show the oxygen residual after the subtraction of the linear best fit from (b) plotted against the sewage flow rate as measured at the Macaulay pump station, again the gray circles show the residual after the seasonal best fit line (not shown) is subtracted. ... 32 Figure 6: (a) The two source water masses (green and yellow) from July 2015 with the outfall mooring water in the middle (red). Mooring data points (red outfall and yellow JdF-W) are hourly data while the ferry data (green) is one point per crossing of the ferry data area (Figure 1a). (b) The fraction of deep JdF water for the three parameters (T, S, O2) at the outfall mooring. ... 36 Figure 7: Hypoxia (<63 µmol kg-1_{) risk at the outfalls and surrounding region from the} potential declines of oxygen in the source water masses assuming the ratio of source water masses does not change in the future. The yellow area represents the combination of conditions required to give short term (hourly) hypoxia during summer flood tides (85% deep JdF water). The orange area represents the combination of conditions required to give summer seasonal hypoxia (60% deep JdF water). Current summer means and standard deviations from the surface of Haro Strait and deep JdF (Table 4) are marked by the cross in the upper left. The horizontal and vertical scales are different. ... 37 Figure 8: (a) Comparison of TA and DIC collected near the outfall sites (colored diamonds) with past observations from JdF (black dots) and Haro Strait (gray stars). (b) Calculated pH from CO2SYS. (c) Aragonite saturation state calculated using CO2SYS (Heuven et al. 2009). ... 41 Figure 9: Sensitivity of final model result to various parameters. The vertical red arrow represents my base case (Table 2 and 3), a conservative yet plausible scenario of current sewage flow rates, calm ocean currents, a deep trapping depth, and high estimates of sedimentation. The green arrow shows approximately the oxygen deficit if all POC and sediment organic carbon was removed. All parameters are varied independently in relation to the base case. ... 42 Figure B1: The rocky bottom near the Clover Point outfall as observed by the ROV. .... 63 Figure B2: The soft sediment bottom near the Macaulay Point outfall. Note the low visibility. ... 64 Figure B3: Diffuser port at Clover Point outfall showing considerable amounts of growth. ... 64

Figure B4: A Rock fish, circled, finds shelter immediately beneath a Clover Point diffuser port, undeterred by the sewage. ... 65 Figure B5: A Clover Point diffuser showing the large amount of particulate matter being discharged. ... 65 Figure C1: The power spectra as measured by the Macaulay Point Mooring. Tidal cycles dominate for all parameters including oxygen. The sewage flow rate has strong semi-diurnal and diurnal cycles as well corresponding to the daily routine of human activities. ... 67 Figure C2: Example time series of July 2015 oxygen (blue) from the outfall mooring and the sewage flow rate (green) from the Macaulay pump station. ... 68 Figure C3: Example time series of July 2015 oxygen residual (blue, observed oxygen minus salinity-predicted-oxygen) from the outfall mooring and the sewage flow rate (green) from the Macaulay pump station. ... 68 Figure D1: Daily net summer cross currents at the outfall mooring. The solid red line is the median value and the dashed red lines show the 25th (left) and 75th (right) percentiles……….. 69

Figure D2: Daily net summer long currents at the outfall mooring. The solid red line is the median value and the dashed red lines show the 25th (left) and 75th (right) percentiles……….. 70 Figure D3: Spring (doy 77 - 137) currents at outfall mooring………. 70 Figure D4: Summer (doy 138 – 225) currents at outfall mooring……….71 Figure D5: Autumn (226 - 314) currents at outfall mooring……….72 Figure D6: Winter (315 - 76) currents at outfall mooring……….. 72 Figure E1: Spring (left) and summer (right) oxygen profiles from outfall monitoring stations. ... 73 Figure E2: Fall (left) and winter (right) oxygen profiles from outfall monitoring stations. ... 74 Figure E3: Spring (left) and summer (right) salinity profiles from outfall monitoring stations. ... 74

Figure E4: Fall (left) and winter (right) salinity profiles from outfall monitoring stations. ... 75 Figure E6: Fall (left) and winter (right) temperature profiles from outfall monitoring stations. ... 76 Figure E7: Macaulay Point, Station 41, April 8 2013, 11:37am The sewage layer is near 40 m depth with another possible sewage layer near 48 m depth. ... 77 Figure E8: Macaulay Point, Station 42, July 22 2013 2:16pm. Below 52 m is likely a sewage influenced layer. ... 77 Figure E9: Macaulay Point, Station 0 (directly above sewage diffuser), November 3rd 2015, 3:02pm. Two sewage influenced layers are probably occurring, the more pronounced layer at 37 m with a second less pronounced layer at 48 m. ... 78 Figure E10: Clover Point, station 42, November 19th _{2015, 12:10pm. Two week} influenced sewage layers are likely, one near 52m and a second at 57 m. ... 78 Figure G1: Data collected by the Institute of Ocean Sciences at station ADCP (48.233, -123.300), roughly 10 km east of JDF-East. Nutrient samples were analyzed according to Barwell-Clarke and Whitney (1996). ... 83 Figure H1: Scatter of temperature and salinity (left) and oxygen and salinity (right) from around the Haro Strait and JdF region in September 2013. ... 85 Figure H2: Scatter of temperature and salinity (left) and oxygen and salinity (right) from around the Haro Strait and JdF region in September 2014. ... 85 Figure H3: Scatter of temperature and salinity from around the Haro Strait and JdF region in June 2014 (R) and June 2013 (L). ... 86

Acknowledgments

This project would not have been possible without the support of many people. I would first like to thank my co-supervisors Roberta Hamme, Debby Ianson, and Christopher Lowe. Without their help and encouragement this research would not have been possible. I would also like to acknowledge the contributions from Ocean Networks Canada, in particular Richard Dewey and Ken Denman, for their support and encouragement before this project even started. I would also like to thank CRD-Marine Programs for being such a wonderful sponsor to this research and for the helpfulness and willingness of marine programs staff to make data available, assist in sampling, take me along on their sampling, and answer my many questions. The experiences I have gained working with the CRD’s marine program have been immensely helpful, and I have learned a great deal. Debby’s funding from the Department of Fisheries and Ocean Canada, Aquatic Climate Change Adaptation Service Program (ACCASP) paid for the skillful analysis of my carbon samples by Marty Davelaar at the Institute for Ocean Sciences. Funding for the project was provided by a NSERC-IPS award and the by CRD-Marine Programs. Field observational support was provided by Ocean Networks Canada.

Key Words

Hypoxia, screened sewage, time series, ocean acidification, dissolved oxygen, anthropogenic, Victoria BC

Highlights

• Victoria B.C., a city of nearly 350,000 inhabitants, discharges screened wastewater. • Effluent is rapidly diluted by currents, and cold temperatures slow respiration. • Impacts on dissolved oxygen are limited to a few µmol kg-1_{; conditions are not} hypoxic. • Declining oxygen in upwelled water from offshore could put Victoria at risk of hypoxia. • Higher levels of treatment are unlikely to significantly increase dissolved oxygen.

1 Introduction

The discharge of minimally-treated municipal wastewater into the coastal marine environment can have disastrous impacts, and much regulatory attention is paid to the quality and quantity of the effluent (Canadian Fisheries Act). However, the properties of the receiving environment are critical in determining the ultimate impacts of the waste (Dinn et al., 2012; Chapman et al., 1996). If the waste is discharged close to shore or in a constrained environment, the public health risks are considerable. If discharged into a nutrient-limited environment, as is often the case for freshwater (Lake Winnipeg, Schindler et al., 2012) or sheltered marine ecosystems (Baltic Sea, Ronnberg and Bonsdorff, 2004), eutrophication can result, and the frequency of harmful algal blooms can increase. The rapid growth of aquatic plants can poison wildlife, contaminate drinking water, and deplete oxygen in lower levels of the water column (Smith and Schindler, 2009; Tchobanoglous, 1991). However, if discharges are into a non-nutrient-limited environment, impacts on oxygen are generally smaller. The organic load within the sewage itself can still be considerable, and, if the respiration of sewage organic material occurs in a small area, significant oxygen depletion can occur within that zone. As with all respiration, inorganic carbon is released, lowering the pH of water and worsening ocean acidification in the local marine environment (Cai et al., 2011). Low oxygen, combined with lower pH, and warmer water temperatures has the potential to act synergistically to the detriment of many marine organisms (Haigh et al., 2015).

Chemical contamination of the receiving environment with metals and persistent organic pollutants is of universal concern (Tchobanoglous, 1991), and, with the recent advent of more precise analytical methods, environmental contamination by pharmaceuticals and personal care products (PPCPs) has now garnered attention (Boxall et al., 2012). However, these contaminant concerns are often secondary to the eutrophic and oxygen-depleting impacts of domestic sewage (Tchobanoglous, 1991).

1.1 Study Site

Victoria, the capital city of British Columbia, Canada, is situated on the southern tip of Vancouver Island and is surrounded by the waters of Juan de Fuca (JdF) and Haro Straits (Fig 1a). The area served by greater Victoria’s sewage system is lightly industrialized and home to nearly 350,000 inhabitants (Statistics Canada, 2012). The region has come under increasing public and international pressure to upgrade its sewage treatment system, which currently consists of only preliminary sewage treatment (6mm screen) before allowing the wastewater to be discharged into the marine environment of JdF Strait (CRD, 2014). The region is served by two major sewage outfalls located at Clover and Macaulay Points (Figure 1a inset). The outfalls extend 1.1 and 1.8 km offshore, respectively, to depths of 62 m and 67 m, where wastewater is discharged through a series of small diffuser ports along the final 150 m of the outfall pipe (Figure 1a inset) designed to

with Clover experiencing slightly stronger currents (Chandler, 1997). The dominant bottom current direction, as determined by mooring studies, is to the southeast, matching a general pattern of seafloor sediment contamination (Dinn et al., 2012; Chandler, 1997; CRD, 2014; Chapman et al., 1996). Water temperatures are cool, varying seasonally between 7 and 10oC (Chandler, 1997). The seabed around the Clover Point outfall is rocky and sedimentation is minimal (Diaz, 1992; CRD, 2014; Markovic, 2003). Around the Macaulay outfall, the seabed is made up of a mixture of sediment size classes (Diaz, 1992; Dinn et al., 2012a; Krepakevich and Pospelova, 2010; Markovic, 2003) ranging from fines to gravel. The area around the Macaulay outfall is also littered with garbage from when the area served as greater Victoria’s municipal garbage dump between 1908 and 1958 (Ringuette, 2009). Some of this garbage (mainly glass, ceramics, and leather) is now buried a few centimeters below the surface (Shirley Lyons, scientific officer CRD, personal communication), indicating at least patchy sedimentation in the area (Appendix A).

1.2 Regional Physical Oceanography

The body of water that receives Victoria’s wastewater, JdF, is a long (150 km), deep (250 m), and wide (25 km) channel connecting the Pacific Ocean to the catchment basins of the Strait of Georgia (SoG) to the North-East and Puget Sound to the South (Figure 1a). Collectively these three main basins form an area known as the Salish Sea. The innermost basins, SoG and Puget Sound, receive large volumes of freshwater runoff

(LeBlond, 1983), and together the three basins form an estuarine complex, with JdF acting as an outer estuary (Thomson et al., 2007). The largest source of fresh water enters the complex via the SoG from the snow melt fed Fraser River (1000 – 9000 m3 s-1), whose flow peaks in early summer (May-June) and reaches its annual minimum in late winter, with a secondary peak occasionally occurring in late fall due to heavy rains (Morrison et al., 2002). Puget Sound, meanwhile, receives approximately 2200 m3 s-1 of runoff, with no strong seasonal cycle (Thomson et al., 2007). The total volume transported into and out of JdF by estuarine circulation changes little throughout the year. Seawards flows vary between 105,000 m3s-1 and 118,000 m3s-1 in the winter and summer respectively, while landward flows are stable near 87,000 m3s-1. The mean transition depth between seaward and landward flows is 60 m (Thomson et al., 2007), nearly the exact same depth as Victoria’s sewage outfalls.

Figure 1: (a) Map of the region, the colored stars show the locations of sea bottom moorings, the hashed area is a zone of high surface nutrients where primary production is light limited year round (Mackas & Harrison 1997). (b) Close-up map of the sewage outfalls with the black dots showing the location of monitoring stations where DIC/TA samples were collected and vertical profiles taken. (c) A cross-section of the average flow and approximate depth along the thalweg. The outfall mooring is not located along the thalweg, but is still attached to the seabed, north of the thalweg nearer to shore. Along the western shelf off Vancouver Island, summer winds bring upwelling (Crawford and Thomson, 1991) and allow hypoxic, nutrient-, and carbon-rich shelf-slope waters to enter the landward-flowing bottom layer of JdF (Davenne and Masson, 2001). The shelf-slope water moves south-east along the bottom of JdF as far as the Victoria sill before Ferry%Data%Zone Juan de Fuca Strait Victoria Sill Haro Strait Strait of Georgia Shelf Slope JdF-West (225 m) JdF (182 m) JdF-East (112 m) Outfall (55 m) Boundary Passage (224 m) Summer Upwelling 60 m Fraser River Approximate Mooring Locations (c) (b) (a) High Nitrate (10 – 30µM)

being significantly altered by mixing with surface layers originating in the SoG (Pawlowicz, 2001). In winter, downwelling winds predominate along the outer coast, and oxygen concentrations rebound in the bottom layers of JdF (Davenne and Masson, 2001). Important to the region’s oxygen dynamics is Haro Strait (Figure 1) that connects JdF with SoG. Haro Strait is a complex of shallow sills and narrow channels that induce mixing and allow the Strait of JdF to operate as a separate estuarine circulation cell (Figure 1b). Surface waters enter Haro Strait’s north-eastern end from the SoG, while deep waters from JdF enter Haro Strait’s southern end. Mixing of these water masses injects nutrients from deep Jdf water into the surface layers of Haro Strait, much of which, subsequently, flows into the surface layers of JdF, creating a large zone where primary production is not nutrient limited (Mackas and Harrison, 1997). At the same time, the vigorous mixing allows large amounts of oxygen from the atmosphere to diffuse into the water and be mixed throughout the water column (Ianson et al., 2016). Mixing is so complete that stratification breaks down, forcing the majority of water entering Haro Strait from JdF to be mixed into the surface and exit Haro Strait in the outward-flowing surface layer. Thus, much of the JdF landward flowing water never reaches the SoG (Pawlowicz et al., 2007).

2 Methods

2.1 Observational

2.1.1 Vertical Profiles Ship-based (MV John Strickland) monitoring near the outfalls was carried out on five days each season between fall 2011 and fall 2016. Stations were located approximately 100 meters from the sewage diffusers (black dots, Figure 1c). On each sampling day, four stations were sampled based on the conditions at the time in order to maximize the chance of plume detection (i.e. in the downstream current direction as predicted by an oceanographic model, Hodgins, 2006). Vertical profiles (Appendix E) were collected by Seabird SBE-19V conductivity, temperature, and pressure (CTD) sensor (± 0.0007 on the PSS-78 scale, ± 0.005 °C, ± 0.1 %) equipped with an SBE-43 oxygen sensor (±2 % of saturation). Niskin bottles attached to a rosette were used to collect water samples from three depths (roughly surface ~5 m, mid-depth ~40-50 m, and ~5 m above the bottom). Water samples for fecal coliform analysis were stored on ice and processed within 24 hours. Fecal coliforms were enumerated using 0.45 µm membrane filters placed in mFC growth medium and incubated at 44.5oC for 24 hours (CRD, 2014). 2.1.2 Inorganic Carbon and Total Alkalinity

On four occasions (July 22nd 2015, November 9th 2015, July 26th 2016 and January 24th 2017) discrete dissolved inorganic carbon (DIC, ±2 µmol kg-1) and total alkalinity (TA, ±4

select profiles (Figure 1c) and were analyzed according to Dickson et al. (2007). Values of pH and aragonite saturation were calculated using the MATLAB® version of CO2SYS (Van Heuven et al., 2009). Nutrient concentrations (phosphate and silicic acid) required for these calculations were estimated from their relationship with salinity in the region (Appendix G), with the exception of January 24th 2017 for which discrete nutrient data (Barwell-Clark and Whitney, 1996) were collected along with the DIC/TA. 2.1.3 Mooring Network In February of 2013, a bottom-mounted mooring (instruments approximately 6 meters off the seabed) was placed approximately 200 meters downstream (southeast) of the less energetic (Chandler, 1997) Macaulay Point sewage outfall (Figure 1c). The mooring was equipped with a Sea-Bird SBE-37SMP CTD (± 0.0006 PSS-78, ± 0.002 °C, ± 0.1 %), SBE-63 oxygen optode (greater of ±3 µmol kg-1 or ±2 %), and a Nortek Aquadopp acoustic current meter (±1 cm s-1_{, ±2}o_{), making this mooring identical to four others} placed around the region by Ocean Networks Canada (www.oceannetworks.ca) in the fall of 2012 (Figure 1a). Individual mooring deployments lasted anywhere from four to twelve months (Appendix F). Combined, these deployments make up a four-and-a-half-year time series spanning November 2012 to October 2016, with moorings still deployed.

2.1.4 Ferry Box Underway System In addition to the mooring data set, I make use of surface measurements collected by an underway system aboard the BC Ferries vessel Spirit of Vancouver Island as it crosses the ‘Ferry Data Zone’ box at the north end of Haro Strait (Figure 1a). The system collects data on temperature and salinity with a SeaBird SBE 21 SeaCAT Thermosalinograph (± 0.01 °C, ±0.005 PSS-78) and oxygen with an Anderaa optode 3835 (greater of ±5 uM or 5%) from an intake port approximately 3 meters (depending on ferry load) below the surface. Data from the ferry system were obtained from www.oceannetworks.ca. 2.1.5 Remotely Operated Underwater Vehicle (ROV) In August 2012, an ROV conducted an engineering inspection of both outfalls. The ROV traced the same route several times but did not maintain a constant speed or distance above the bottom. All individual organisms visible in the video were counted (Appendix B). Due to the nature of the dive and the low visibility, some double counting likely occurred. Furthermore, because the ROV did not maintain a constant distance above the bottom nor a constant speed, the smallest countable organism was not constant with time. The video also provided qualitative information on the physical attributes of the seabeds near the outfalls as well as on the behaviour of wastewater particles in the cold marine environment.

2.2 Data Analysis methods

The ferry data (Section 2.1.4) and the moored current meter data (Section 2.1.3) were pre-processed (converted from raw voltages) by ONC and obtained from www.oceannetworks.ca. All other electronic data were processed using the manufacturer’s software and calibrations. Mooring data from the various deployments were merged together and binned to 1-hour intervals. Linear corrections were applied to correct oxygen sensor drift (~1 % yr-1) as quantified by routine manufacturer calibrations (Appendix F). Sensor drift of the CTDs was not measurable, and no correction was applied. 2.2.1 Time Series Analysis Because sub-surface oxygen concentrations in the region are lowest during the summer upwelling season (Section 3.1.1), I focused my analysis on summer. I define summer based on the timing of upwelling near the western entrance of JdF, which typically spans year-day 138 – 225 (Bylhouwer et al., 2013). Hourly mooring data (T, S, O2) were filtered with a 25-hour running mean to remove tidal effects (Section 3.1.1). In addition, power spectral analysis (Welch, 1967) was carried out on the data from the outfall mooring (Appendix C) and on the time series of sewage discharge over the same time period.

2.2.2 Water Mass Endmembers Two primary water mass endmembers were identified from the outfall mooring dataset (Section 3.1.1). The first is cold, salty, and oxygen-poor water originating along the shelf slope and advected with the estuarine return flow in deep JdF (Masson, 2006). I assign the deep JdF water mass the average summer properties (T, S, O2) observed at 225 m by the mooring located at the western entrance of JdF (JdF-West, Figure 1a). I define the second water mass as a hybrid of deep JdF and surface SoG water (Section 3.1.2), and I assign it the summer averages (T, S, O2) observed by the ferry system (Section 2.1.4) as it crosses the northern edge of Haro Strait (Figure 1a).

2.3 Model Methods

To further assess and quantify the impact on oxygen and DIC near the outfall, a two-box box model was created. The large amount of observational data available (Section 2.1) allowed us to keep the physical model simple by fixing the boundary conditions to observations. The model simulates salinity (S), temperature (T), and sewage-derived dissolved organic carbon (DOC) in the upper and lower boxes (Figure 2). Dissolved oxygen (O2) and sewage particulate organic carbon (POC) are simulated in the lower box only. Model results were always compared to a base case (Tables 3 and 4) simulation with no added sewage. The difference in dissolved oxygen between the two model runs is the oxygen deficit attributable to the sewage.

Model parameters listed in Tables 1 - 3 are fixed for each model run but varied one at a time between runs (within the ranges stated) to assess the model’s sensitivity. The model is run to steady state in MATLAB® using a Runge-Kutta solver with a dynamic time step (Shampine and Reichelt, 1997). Directly solving the steady state equations analytically yields the same results as the MATLAB® solver. Like all box models, both boxes are perfectly mixed. 2.3.1 Physical Model The water column (60 m in total) is divided into an upper and a lower box. The transition depth between boxes is the trapping depth (DTrap), defined as the depth

above which the majority of sewage does not rise. DTrap is different and much shallower

than the transition depth between landward and seaward flowing layers in JdF that occurs around 60 m (Thomson et al. 2007). Past bacteriological observations (Section 2.1.1) show that the sewage plume is normally trapped below 30 m depth (DTrap = 30

m)(CRD, 2015), but at times can be trapped below 40 m (DTrap = 40 m) or, rarely (<2 %

of observations), it can rise to the surface (DTrap = 0 m)(CRD, 2014; Lorax, 2009).

Near-surface (5 m) fecal bacteria observations are used in conjunction with observed concentrations of fecal bacteria within the sewage to calculate the amount of sewage dilution and, thus, the vertical transport of sewage (i.e. vertical mixing) into the surface layers. Despite fecal bacteria being a non-conservative tracer of sewage (die-off rates can exceed 90% in 12 hours (Lorax, 2012)), they are the only measured parameter that

is consistently able to detect the sewage plume. The amount of vertical mixing is fixed for each model run; however, it is varied over a considerable range for sensitivity analysis (Section 3.2). Within the lower box, I simulate a reflective bottom sediment layer (i.e. instantaneous oxidation of all organic material) and a pool of slow-moving suspended particulate organic carbon (POC, Figure 2) equivalent to particles trapped in a slow-moving bottom boundary layer (Thomsen, 2002). The reflective sediment layer simulates a steady state of organic deposition, where the input rate equals the respiration rate (Bianucci et al. 2011). This method of modeling sediment does not take into account any burial which may be occurring. The pool of suspended sediment is respired with the same rate constant as that used for the dissolved pool, 0.7 day-1. The lower box is flushed with water that has properties (T, S, O2, Table 4) roughly matching the summer means as observed by instruments on the outfall mooring (Section 3.1.2). The upper box experiences the same advective flushing as the lower box, but the water is warmer and fresher to match observations of the upper 30 m of the water column (Section 2.1, Appendix E). Because tidal currents loop back on themselves the rate of flushing is determined by the summer net daily transport (daily sum of hourly current velocity observations) as observed by the Macaulay Point mooring (Appendix D). By doing this the model is effectively subjected to only the non-tidal residual currents, slower than if the tidal currents are included, this is a

conservative assumption and increases any oxygen deficits. The boxes are aligned such that the primary axis of the summer currents (south-east, 110o_{) is parallel with the long} axis of the boxes. Salinity is converted from PSS-78 to kg m-3 _{because the model is} volume based; this conversion was done using the TEOS-2010 Gibbs-SeaWater Oceanographic toolbox (McDougall and Barker, 2011) and depends slightly on the density of the water, but for the conditions near the outfall the conversion is nearly 1:1.

The upper box’s oxygen concentration, O2 upper, is fixed at 175 µmol kg-1, a low but

frequently observed value at 30m depth, in order to simulate a worse case for low oxygen. The fixed oxygen concentration in the upper box acts as a source of oxygen to the lower box through vertical mixing (vz). Oxygen sources of primary production and atmospheric gas exchange are implicitly modeled by the fixed oxygen concentration in the upper box. The spatial domain of the model is constrained to 800 m long (x) by 200 m across (y), such that (1) the most seriously polluted sediment, which has been shown to be within several hundred meters of the outfalls (Diaz, 1992; Dinn et al., 2012; CRD, 2014), is within the spatial domain of the model, and (2) so that the water column dilution of sewage is minimal due to the small spatial scale. By constraining the domain for these reasons, I simulate the largest reasonable oxygen deficit. Larger domains allow further dilution of the effluent and therefore smaller oxygen deficits, while much smaller

Figure 2: (a) Cartoon of the sewage plume and processes affecting T, S, O2 and DIC the outfall within the physical structure of the mass balance model. Red arrows show processes that consume oxygen and green arrows show process that produce DIC. Modelled quantities are listed along the left side of each physical layer and in Table 1. In the model, these quantities are perfectly mixed within their respective layers. Daily averaged net ocean currents are indicated for the long direction, cross direction currents are not shown but range from 0.9 – 3.5 km day-1 ₍₂₅th _and 75th percentiles). (b) Schematic of physical and biological modelled fluxes in the lower layer of the model showing the flow of O2and consumption of organic matter. Solid arrows represent volume fluxes that affect the O2concentration while dotted arrows represent biological fluxes (consumption). Temperature and Salinity (not shown) only experience physical fluxes, i.e.; advection, vertical mixing and dilution by sewage. 400 Surface Mooring Water Column O2Demand DO fixed at 175μmol/kg Dtrap (0 – 40m) 60 -DTr ap x (800m) + DIC + DIC vz(1-10 m/day) vx (4.5 - 9.2 km/day) Sediment O2Demand + DIC Suspended Sediment (POC) O2Demand vx (4.5 - 9.2 km/day) O2 env Senv lower Tenv lower Senv upper Tenv upper O2 lower DOClower POCsuspended Slower Tlower Supper Tupper DOCupper Sewage (F) Oxygen Advection (Vx) Vertical Mixing (VZ) POCSuspended Sewage DOC Sewage (F) Respiration of DOC Sediment O2Demand Advection (Vx) + Sewage (F) Vertical Mixing (VZ) Advection (Vx) + Sewage (F) Advection (Vx) + Sewage (F) Respiration of POC (a) (b)

Table 1: Modeled quantities.

Name Simulated Prescribed units

Dissolved Oxygen O2 lower O2 env lower, O2 sewage, O2 upper

mmol m-3

Salinity Supper, Slower Senv upper, Senv lower, Ssewage kg m-3

Temperature Tupper, Tlower Tenv upper, Tenv lower, Tsewage oC

Organic Carbon

(dissolved sewage) DOCupper, DOClower -- mmol m

-3

Organic Carbon Particles (sewage particles)

POCsuspended -- mmol m-3

Table 2: Physical dimensions and parameters of model.

Name Symbol Units Value Range Base Case

Horizontal long dimension x m 800 800 Horizontal cross dimension y m 200 200 Trapping Depth (measured down from the surface) DTrap m 0 – 40 40 Daily net long direction advective transport vx m day -1 4500 – 9200 4500 Daily net cross direction advective transport vy m day -1 900 – 3400 900 Daily net vertical transport vz m day -1 1 -10 5

Table 3: Biological and anthropogenic forcings.

Name Symbol Units Value

(Range) Base Case Reference Sewage Flow Rate F m 3 _day -1 40,000 – _60,000 50,000 CRD, 2014 Sewage carbonaceous biochemical oxygen demand CBOD mmol m-3 7600 -- CRD, 2014 Sewage biochemical oxygen demand BOD mmol m-3 8500 8500 CRD, 2014 Sewage Temperature Tsewage o_{C 20 20 CRD, 2014} Sewage Salinity Ssewage kg m-3 0 0 CRD, 2014

Sewage Oxygen Osewage µmol

kg-1 0 0 CRD, 2014 Molar ratio DIC : O2 r -- 0.7 - 0.9 -- CRD, 2014 Temperature correction constant 𝜃 -- 1.047 - 1.135 1.135 Tchobanoglous, 1991 Sewage decay rate k day -1 _{0.23 – 0.7 0.7 Tchobanoglous,} 1991 DOC Fraction of sewage BOD a -- 0.80 - 0.94 0.80 3.1.3, Tchobanoglous, 1991 POC Fraction of sewage BOD b -- 0.05 - 0.15 0.15 3.1.3, Tchobanoglous, 1991 OCsediment Fraction of BOD c -- _{0.01 -} 0.05 0.5 Dinn et al., 2012b Fractional speed of POC d -- 0.1 - 0. 5 0.1 3.1.3 Temperature and salinity in the lower box are simulated according to equations 1 and 2 "#$%%&' "( =

*+$%%&' #&,-.#$%%&' /01(#345&'.#$%%&')

789$%%&' (1)

"#345&'₌*+345&'#&,-.#345&' /01 #$%%&'.#345&' /:(#;&5<=&.#345&')

where J represents either temperature or salinity concentration and F is the sewage flow rate. Subscript env indicates the environmental concentrations (Table 4) of salt, temperature, and oxygen outside of the model domain while subscript sewage indicates the properties of sewage.

𝐻𝐴@AABC = 𝑥𝑣F𝐷HCIA + 𝑦𝑣L𝐷HCIA (2.1)

𝐻𝐴_98MBC = 𝑥𝑣_F(60 − 𝐷_HCIA) + 𝑦𝑣_L(60 − 𝐷_HCIA) (2.2) 𝑉𝑀 = 𝑥𝑦𝑣S (2.3) 𝑣𝑜𝑙_@AABC = 𝑥𝑦𝐷_HCIA (2.4) 𝑣𝑜𝑙_98MBC = 𝑥𝑦(60 − 𝐷_HCIA) (2.5) Although salinity does not impact dissolved oxygen, it was included for completeness and as a check on dilution of fresh wastewater in the marine environment. Table 4: Prescribed values of state variables at the model’s boundaries.

Parameter Units Upper Layer Lower Layer

Tenv oC 12.3 9.5 Senv kg m-3 29.5 31.5 O2 env mmol m-3 175 150 2.3.2 Oxygen Depletion Oxygen depletion in the water column due to a sewage outfall occurs in two steps. First, the dilution of the ambient environmental oxygen with anoxic sewage (the sewage system has elevated levels of H2S) creates an immediate drop in dissolved oxygen. Following this initial decline in nearfield oxygen concentration, organic material (DOC

the following days to weeks according to first order kinetics with rate constant k (common to both DOC and POC for my model) and a temperature-correcting coefficient θ (Tchobanoglous, 1991). Compared to natural sources of marine organic carbon raw sewage is labile with respiration rates up to 0.7 day-1_{, compared with about 0.005 day}-1 for the semi-labile fraction of marine DOC and 0.2 day-1 for marine POC (Ianson and Allen, 2002). The total amount of oxygen consumed from the ambient environment is the product of the sewage flow rate (F) and the sewage’s biochemical oxygen demand (BOD). Because the currents at the study site are so strong, dilution is able to continue indefinitely, so, unlike a river system which can see maximum oxygen deficits well downstream of the wastewater discharge (Tchobanoglous, 1991), I expect the largest oxygen depletion to be close to the outfall where water column concentrations of sewage are highest and where sewage particles are able to settle to the seafloor in significant quantities. These particles, which I simulating using a reflective boundary, enrich the sediment within about 600 m of the outfalls (Diaz, 1992; CRD, 2014) and therefore enlarge the sediment’s oxygen demand. Within the sewage itself I simulate three types of biodegradable organic carbon to best match qualitative observations from the ROV dive (Section 3.1.3) and what is known about the physical nature of the waste. The first type (DOCupper and DOClower) is

(POCsuspended) consists of small sewage particles which sink and become suspended in

the bottom boundary layer, moving along the seabed at a much-reduced speed, which I model via a fractional reduction in the ambient current speed (d). Lastly, the largest and densest organic carbon particles (OCsediment) quickly sink to the bottom of the seabed

and are not re-suspended. Because I assume steady state, OCsediment is instantaneously

oxidized in a reflective boundary layer (Bianucci et al., 2011). I do not simulate any burial or accumulation of organic carbon on the seabed to further maximize simulated oxygen deficits. Natural marine pools of DOC and POC are not simulated because they are assumed to be unaffected by the sewage and thus have no impact on the size of the oxygen deficit caused by the sewage (i.e. in both model runs with and without sewage the natural sources would be the same and thus cancel out). Dissolved oxygen, which is fixed in the upper box (175 µmol kg-1), is simulated in the lower box according to equation 3 "VW 345&' "( = eqn 2 − 𝑘𝜃 H.^_ _𝐷𝑂𝐶 98MBC + 𝑃𝑂𝐶c@cABd"B" −Ve₇₈₉;&fgh&hi 345&' (3) where 𝑂𝐶_cB"jkBd( = 𝑐 ∗ 𝐹 ∗ 𝐵𝑂𝐷 (3.1)

Dissolved organic carbon, (DOClower and DOCupper) are simulated in both the upper and

lower boxes according to equations 4 and 5 respectively.

"pVe$%%&' "( = pVe345&'∗01 789$%%&' −

pVe$%%&'∗ *+$%%&'/01

789$%%&' − 𝑘𝜃 (H$%%&'.^_)𝐷𝑂𝐶 @AABC (5) Suspended particles are only present in the lower box and are simulated according to equation 6. "rVe;$;%&,f&f "( = s∗:∗qVp 789345&' −

"∗*+345&'∗rVe;$%&,f&f

789345&' − 𝑘𝜃 (H345&'.^_)𝑃𝑂𝐶 c@cABd"B" (6) Primary treatment (settling tanks and grease skimmers) is able to reduce the BOD of wastewater by 25-40% (Tchobanoglous, 1991), so it is reasonable to assume that up to 40% of the Victoria’s wastewater BOD load is in the form of suspended sewage particles (POCsuspended) or sediment (OCsediment). However, because the ambient environment

around Victoria’s sewage outfalls is much higher energy than a settling tank, I assume that a relatively larger fraction of the wastewater’s BOD (at least 80%) will be in the functionally dissolved phase (DOC). Again, for the purposes of this study, I define dissolved based on the behaviour of any sewage particles (i.e. if it moves freely with the currents and does not sink it is dissolved), not filter size. Although the ROV video does not allow for the quantification of settling rates, it does confirm that some settling of sewage particles is occurring. Using the reasoning above and the ROV observations, I limit my simulations to POCsediment making up 5-15% of the total BOD load and true

organic material only making up at most 20% of the total sewage organic load, these two fractions have much longer residence times near the outfalls and therefore exert a disproportionately large impact on both the predicted oxygen deficit and DIC (Section 3.2).

3 Results and Discussion

3.1 Observational

3.1.1 Vertical Profiles Very few of the vertical profiles collected near the sewage diffusers show physical changes in the water column that indicate a sewage plume. Upon visual examination of nearly 850 profiles, only six profiles showed any evidence of a sewage plume layer, that is a layer that is fresher, warmer, more turbid, and lower in oxygen than the layers directly above or below it. From within the subset of profiles that did show plume effects (Appendix E), the largest declines in oxygen were close to sensor accuracy, at about 2 µmol kg-1_{, a drop of this magnitude would be expected to depress pH by} roughly 0.005 units (assuming a DIC enrichment of 1.7 µmol kg-1 _{and fixed TA of 2125}

µmol kg-1 _{at a salinity of 30 and a temperature of 10}o_{C). Despite the lack of oxygen}

deficits in the vertical profiles, bacteriological data collected at the same time as the electronic profiles frequently shows high fecal coliform (up to 6,000 CFU 100mL-1_{) in the} mid and deeper layers of the water column. These high bacterial counts confirm that sampling occurred in the correct downstream direction (CRD, 2014). Post screening but before the sewage is discharged, fecal coliform levels can approach 6,000,000 CFU per 100mL (CRD, 2014); the marine environmental concentrations observed at the monitoring stations represent dilution ratios of several thousand to one, if they can be quantified at all (i.e. many samples from the upper water column have bacteria counts comparable to natural levels making dilution calculations impossible). Dilution rates are highly variable, with the lowest observed dilution rates of around 500:1 occurring about once per year (~0.5% of observations).

3.1.2 Moorings The mooring dataset shows, for the first time, the full seasonal cycle of dissolved oxygen, not only at the outfall site but also in JdF. Previous mooring deployments in the region were not equipped with oxygen sensors. Although the general seasonal cycle has long been known from ship surveys (Herlinveaux and Tully, 1961; Johannessen et al., 2014; Masson, 2006), the mooring dataset (Figure 3) reveals significant temporal detail that was previously unknown. Figure 3: Time series of dissolved oxygen, temperature, and salinity from the four-year mooring dataset filtered with a 25-hour running mean to remove semi-diurnal and diurnal tidal effects. The horizontal dashed magenta line indicates a hypoxic threshold of 63 µmol kg-1. All moorings are anchored to the bottom with instruments approximately 6 m off the seabed, with exact depths of each mooring shown in figure 1b and listed in Table F1.

The outer moorings (JdF-W and JdF, orange and green) are heavily influenced by large scale summer coastal upwelling that occurs along the western coast of North America.

lower temperatures, lower oxygen, and higher salinities during the summer (Figure 3). The inner moorings, located at Boundary Passage (BNDYP, gray) and at the outfall site (red), show less influence of upwelled waters and more of waters originating in the SoG, with oxygen minima still occurring in summer. The mooring located at the eastern entrance to JdF (JdF-E, blue) shows a mixture of both the inner and outer mooring characteristics. In winter, downwelling predominates along the outer coast, causing temperature and oxygen to rebound at JdF-W. Intense winter storms combined with lower Fraser River discharge produces periodic estuary flow reversals (Thomson et al., 2007), causing surface waters to flow landwards and deep waters to flow towards the Pacific. These reversals likely cause the periodic drops in salinity seen in winter at JdF-W and JdF (Figure 3). At the outfall site (red line, Figure 3), the mooring provides temporal resolution that ensures that low oxygen events are not simply missed by comparatively sparse ship-based monitoring. When compared with the other moorings (Figure 3), the outfall site has the most oxygen throughout the time series (mean of 183 µmol kg-1 _{compared to} BNDYP’s mean of 172 µmol kg-1_{) and never approaches the commonly used hypoxic} threshold of 63 µmol kg-1 _{shown by the dashed magenta line in Figure 3. However,} hypoxic thresholds should be treated with caution as every species has differing limits of low oxygen tolerance, with some species requiring much more oxygen than 63 µmol kg-1 _{while others can tolerate less (Hofmann et al., 2011; Vaquer-Sunyer and Duarte.,} 2008). Acoustic ocean current observations by the outfall mooring match with those of

past mechanical current meters (Chandler, 1997), showing a dominant current direction to the south-east (110o_{), particularly in summer (Appendix D). High bacterial} concentrations in the lower layers of the water column confirm that not all of the sewage plume rises significantly off the bottom. The dominant direction is part of the larger scale estuarine exchange and suggests that at least some of the sewage may be carried within the bottom layer around Victoria’s waterfront and towards Haro Strait. This observation is contrary to the assumption that the sewage outfalls are shallow enough to ensure that all of the sewage is immediately carried toward the Pacific Ocean in the outward flowing surface layer (Chapman, 2006). However, even if some sewage is able to enter Haro Strait from JdF, the majority of it would likely be mixed into the surface layers and subsequently advected back to the Pacific Ocean ( Pawlowicz et al., 2007) (Section 1.2). The high temporal resolution of the outfall mooring reveals the highly variable nature of bottom layer oxygen near the outfall in all seasons, seen in the width of the frequency distributions (Figure 4). During the summer, when dissolved oxygen is lowest (Figure 3 and 4) and the water column is the most stratified (as determined from vertical profiles, Appendix E), it is common to observe dramatic shifts in oxygen. A single tidal cycle can cause oxygen to vary by nearly 100 µmol kg-1 (Figure 5a). Despite these dramatic changes, oxygen levels were always observed to be well above the frequently used 63 µmol kg-1 _{(2mg L}-1_{) hypoxic (Figure 4) threshold. However, because cold water can hold}

(0.07% of the whole time series), fall below 30 percent saturation (pO2 6.3kPa, Figure 4). Saturation states below 30 percent, especially if persistent, have been shown to have negative impacts on many marine organisms (Hofmann et al., 2011). Autumn (Figure 4c) has a slightly lower mean than summer (149 compared to 157 µmol kg-1) but has less negative skewness resulting in higher minimum concentrations and no observations below 30 percent saturation. Spring and winter (Figure 4a and 4d) have the highest seasonal mean concentrations and saturation states.

Figure 4: The left column shows histograms (probability density functions) of hourly oxygen observations at the outfall mooring, while the right column shows the same information but expressed as percent saturation on the bottom x-axis and partial pressure of O2 on the upper x-axis. The means and standard deviations are in units of µmol kg-1 on the left column and percent saturation on the right column. The fours seasons are spring (a,b, day of year 77-137), summer (c,d, 138-225), fall (e,f, 256-314), and winter (g,h 315-76) (Bylhouwer et al., 2013). 50 100 150 200 250 0 0.05 0.1 0.15 0.2 Mean: 206 STD: 24 30 40 50 60 70 80 90 0 0.05 0.1 0.15 0.2 Mean: 71 STD: 8 50 100 150 200 250 0 0.05 0.1 0.15 0.2 Mean: 157 STD: 26 30 40 50 60 70 80 90 0 0.05 0.1 0.15 0.2 Mean: 55 STD: 10 50 100 150 200 250 0 0.05 0.1 0.15 0.2 Mean: 149 STD: 18 30 40 50 60 70 80 90 0 0.05 0.1 0.15 0.2 Mean: 53 STD: 7 50 100 150 200 250 0 0.05 0.1 0.15 0.2 Mean: 208 STD: 22 30 40 50 60 70 80 90 0 0.05 0.1 0.15 0.2 Mean: 71 STD: 7 14.7 12.6 10.5 8.4 6.3 16.8 18.9 14.7 12.6 10.5 8.4 6.3 16.8 18.9 14.7 12.6 10.5 8.4 6.3 16.8 18.9 14.7 12.6 10.5 8.4 6.3 16.8 18.9 Oxygen Saturation (%) Dissolved Oxygen (μmol/kg)

pO2(kPa) (a) (b) (c) _(d) N o rm a lize d O ccu rr e n ce N o rm a lize d O ccu rr e n ce N o rm a lize d O ccu rr e n ce N o rm a lize d O ccu rr e n ce Spring Summer Autumn Winter (e) (f) (g) (h)

Oxygen and salinity both follow a tidal pattern (Figure 5a, Appendix C) at the outfall; this pattern can be best seen looking at a short period of the hourly time series (Figure 5a). Flood tides bring higher salinity and lower oxygen, while ebb tides bring lower salinities and higher oxygen. A scatter of the same data (black points are data between July 1-15 2015, Figure 5b) shows the strong linearity of the salinity to oxygen relationship with an r2 of 0.92. The entire summer of 2015 (gray points Figure 5b) shows more scatter, but the strong linearity remains (r2 of 0.88). The discharge rate of sewage is highly periodic with two peaks (morning and evening) and two lows (mid-day and overnight) occurring throughout the day (Appendix C). If the sewage was causing periodic draw downs in oxygen, a non-linear relation would likely result. From the salinity oxygen relation, the linear best fit (Figure 5b, seasonal best fit line not shown) was removed and the remaining oxygen residuals were compared to the sewage flow rate (Figure 5c) as measured at the Macaulay Point sewage pump station. Again, if the sewage was having a significant impact on oxygen, negative residuals (i.e. oxygen observations lower than the oxygen concentration predicted based on salinity observations) would correspond with high sewage flow rates, but this coupling is not observed (Figure 5c).

Figure 5 (a) Hourly salinity and oxygen at the Macaulay Point mooring during July 1-15 2015 (b) black circles are the same data as (a) with a linear fit while gray circles are all data from the summer (day of year 138-225) of 2015. (c) black circles show the oxygen residual after the subtraction of the linear best fit from (b) plotted against the sewage flow rate as measured at the Macaulay pump station, again the gray circles show the residual after the seasonal best fit line (not shown) is subtracted. 3.1.2 Water Mass Analysis To compare the conditions observed at the outfall mooring to what would be expected if no sewage outfalls were present, I determined the properties of source water masses that influence conditions at the outfall mooring. The strong linear relations between salinity and oxygen (Figure 5b) and between salinity and temperature (Figure 6a) indicate that only two source water masses are mixing in the area near the outfall to produce the observed variability. Masson (2006) conducted an optimum multiparameter water mass analysis of the entire region (SoG, Haro Strait, and JdF). Her analysis showed that the area just offshore from the outfalls was well defined by two source water masses, deep Juan de Fuca water and surface Fraser River plume water. Following Masson’s work, I defined a salty, cold, low oxygen water mass as the bottom water observed at the JdF-W mooring. Averaging (Section 2.2.2) over each summer

01-Jul 08-Jul 15-Jul

100 120 140 160 180 200 220 240 Oxygen ( 7 mol kg -1) (a) 30 30.5 31 31.5 32 32.5 33 33.5 Salinity 29 30 31 32 33 34 Salinity 80 100 120 140 160 180 200 220 240 Dissolved oxygen ( 7 mol kg -1) July 1-15 2015 R2 _{= 0.92} (b) r2 0 1000 2000 3000 4000 Flow at Macaulay Pump Station (m3 _hr-1₎ -40 -30 -20 -10 0 10 20 30 40 50

Oxygen Residual at Macaulay Mooring (

7

mol kg

-1) (c)

Summer 2015 July 1 - 15 2015

(2013, 2014 & 2015) my results for the deep JdF water mass closely match those found by Masson (Table 5). I found the other end member, that Masson (2006) defined as the Fraser River plume source water, to be too cold and fresh (Figure 6a) in my analysis to be the other source of water observed at the outfall. Instead, I define a unique end member using surface waters along the western edge of Haro Strait (summers of 2015 and 2016, ferry data zone, Figure 1a), which itself is made up of a mixture of SoG, Fraser River plume, and deep JdF waters. Despite the data not being from the main channel of Haro Strait, they do fall along a straight line in temperature-salinity space connecting the deep JdF water with the Macaulay bottom water (Figure 6a). Masson’s Fraser River plume end member is likely too cold and fresh because the outfalls are relatively nearshore features (1.1 and 1.8 km from the coast) and so may not experience the same conditions as those along the thalwag. The water that reaches the outfalls does not come from the middle of Haro Strait, but rather is dragged along the western coast of Haro Strait where it is slowed, warmed, and assimilates more oxygen as it moves towards Victoria (Appendix H).

Table 5: Summer (day of year 138 - 225) mean and standard deviations for the two endmember water masses. * Ferry data from the summer of 2015 was only available for July. Mean Summer JdF-W 2013 2014 2015 2016 Masson 2006 Salinity (PSS-78) 33.91 ± 0.03 33.90 ± 0.03 33.77 ± 0.04 -- 33.9 Temperature (o_{C) 6.3 ± 0.1 6.7 ± 0.2 6.7 ± 0.3 -- 6.4} Oxygen (µmol/kg) 69 ± 7 73 ± 8 84 ± 4 -- 78 Mean Summer Haro Strait Surface (Ferry Data) 2013 2014 2015* 2016 Salinity (PSS-78) -- -- 28.2 ± 0.7 27.8 ± 0.7 -- Temperature (oC) -- -- 13.5 ± 1.4 13.1 ± 1.1 -- Oxygen ç) -- -- 291 ± 30 291 ± 26 -- Mean Outfall Mooring Summer Conditions 2013 2014 2015 2016 -- Salinity (PSS-78) 31.6 ± 0.6 31.6 ± 0.6 31.6 ± 0.6 31.6 ± 0.6 -- Temperature (o_{C) 8.9 ± 0.7 9.3 ± 0.7 9.5 ± 0.9 9.4 ± 0.7 --} Oxygen (µmol/kg) 153 ± 25 167 ± 30 164 ± 23 166 ± 23 -- Interannual variability was generally low, with most seasonal means falling within the range defined by the respective standard deviations. However, conditions at JdF-W did change significantly. Salinity declined by 0.14 (mean seasonal variability 0.02), oxygen rose by 15 µmol kg-1 _{(mean seasonal variability 6 µmol kg}-1_{) and temperature rose by}

0.4oC (mean seasonal variability of 0.2oC). A regional scale warming event known as the ‘blob’ (Bond et al., 2015) likely influenced the mooring network post-November 2014 and contributed to the general warming trend seen through the time series (Figure 3). It is therefore likely that the summer mean temperatures in the surface of Haro Strait were slightly cooler in 2013 and 2014 compared to those observed in 2015 and 2016.

Using the water masses defined in table 5 and assuming that these are the only two source water masses, the average summer mixing fraction of deep JdF water at the outfall mooring was calculated for each of the three measured variables according to equation 7. 𝑂𝑢𝑡𝑓𝑎𝑙𝑙 = 𝑥𝐽𝑑𝐹 + 1 − 𝑥 𝐻𝑎𝑟𝑜 (7) 𝑥 =V@(|I99.*IC8 #":.*IC8 (7.1) where JdF, Haro, and Outfall represent T, S, or O2 at their respective locations, and x is the deep JdF mixing fraction. Differences in mixing fractions between the three variables would indicate non-conservative behaviour (such as oxygen being consumed by sewage discharge) or incorrect end member properties; however, they all agreed within their error (Table 6). Table 6: The mixing ratios as predicted based on the source water properties of table 4 and the conditions observed at the outfall mooring. Due to non-complete data sets the Haro Strait observations from the summer of 2015 are applied to 2013 and 2014. Likewise, the JdF-W 2016 data are taken to be the same as the 2015 observations. Only 2015 has complete data from both the Haro Strait Surface and JdF-W. Mean Summer Deep Juan de Fuca Fraction at the Outfall Mooring 2013* 2014* 2015 2016* Salinity (PSS) 0.59 ± 0.05 0.60 ± 0.05 0.62 ± 0.05 0.63 ± 0.05 Temperature (oC) 0.64 ± 0.08 0.62 ± 0.09 0.59 ± 0.10 0.58 ± 0.10 Oxygen (µmol/kg) 0.62 ± 0.07 0.57 ± 0.08 0.61 ± 0.06 0.60 ± 0.06 The overall agreement between the mixing fractions indicates that no large oxygen deficit is occurring at the outfall mooring; however, the resolution of this method is not good enough to rule out the possibility of a modest oxygen deficit (< 15 µmol kg-1),

Figure 6: (a) The two source water masses (green and yellow) from July 2015 with the outfall mooring water in the middle (red). Mooring data points (red outfall and yellow JdF-W) are hourly data while the ferry data (green) is one point per crossing of the ferry data area (Figure 1a). (b) The fraction of deep JdF water for the three parameters (T, S, O2) at the outfall mooring. In addition to seasonal mean mixing fractions, I used the hourly resolution of the outfall mooring to calculate the water mass fractions on an hourly basis (Figure 6b). These calculations show that there is generally good temporal agreement between temperature, salinity and oxygen (Figure 6b). For brief periods during flood tides, the fraction of deep JdF water found at the outfall can exceed 85%. At these times, the outfall area is most at risk of hypoxia. For example, if oxygen in deep JdF water were to fall below 23 µmol kg-1 _{(where the extension of the current condition axis intersects the} yellow zone, Figure 7), while surface Haro Strait maintains its current level of oxygenation, which is likely given the energetic mixing in the area (Ianson et al., 2016), short term (several hours or less) flood tide hypoxia (<63 µmol kg-1_{) at the outfall site} and surrounding area would be expected (yellow area, Figure 7). Declines in bottom 25 26 27 28 29 30 31 32 33 34 Salinity 6 8 10 12 14 16 18 Temperature (C)

July 2015 Water Masses

(a)

X FR Masson 2006

Outfall Mooring Surface Haro JdF-W

Jul−01 Jul−04 Jul−08 Jul−11 Jul−15 0.4 0.5 0.6 0.7 0.8 0.9 1 Deep JdF Fraction

Outfall Water Mass

(b) Salinity

Temperature Oxygen

and Peña, 2013). Although the drivers of this oxygen decline are not fully understood, and may be oscillatory in nature (Crawford and Peña, 2016), if oxygen concentrations were to continue to decline at the current rate of 0.83 µmol kg-1 y-1, suggested by Crawford and Peña (2013), oxygen in the deep water of JdF could fall below 23 µmol kg -1 _{by the 2050’s. Oxygen concentrations this low in deep JdF would put the near bottom} environment around the outfalls at risk of short term hypoxia, but, not as a result of the outfalls. Figure 7: Hypoxia (<63 µmol kg-1_{) risk at the outfalls and surrounding region from the potential declines of oxygen in} the source water masses assuming the ratio of source water masses does not change in the future. The yellow area represents the combination of conditions required to give short term (hourly) hypoxia during summer flood tides (85% deep JdF water). The orange area represents the combination of conditions required to give summer seasonal hypoxia (60% deep JdF water). Current summer means and standard deviations from the surface of Haro Strait and deep JdF (Table 4) are marked by the cross in the upper left. The horizontal and vertical scales are different. Climate change will likely raise the surface temperatures in Haro Strait (Mote and Salathé, 2010) and in so doing reduce the amount of oxygen in these waters 0 10 20 30 40 50 60 70 80

Dissolved Oxygen in Deep Juan de Fuca (7 mol kg-1) 50 100 150 200 250 300 350

Dissolved Oxygen in Haro Surface (

7

mol kg

-1)

Future Hypoxia Risk

Current Conditions

Flood Tide Hypoxia

(Johannessen, and Macdonald, 2009). Assuming no change to the mean summer oxygen saturation, currently 108% (Table 5) at the ferry crossing (calculated with Garcia and Gordon, 1992, 1993), a +2o_{C increase in surface temperature in Haro Strait would} result in a decrease in Haro Strait surface oxygen concentration from 291 to 279 µmol kg-1_{. Under such conditions, deep JdF water needs to maintain oxygen levels above 25}

µmol kg-1 _{to avoid the risk of flood tide hypoxia (63 µmol kg}-1_{) at the outfalls compared}

to 23 µmol kg-1 with no warming, calculated using equation 7. With +4oC of surface Haro Strait warming (a further decline of surface dissolved oxygen to 269 µmol kg-1) deep JdF water needs to maintain oxygen levels above 27 µmol kg-1 to avoid risk of flood tide hypoxia. The direct impacts of Haro Strait surface warming are comparatively small, such that the oxygen concentration at the outfalls is insensitive to Haro Strait temperature. Seasonally, the average summer deep JdF fraction hovers around 60% (Table 7). Even in the most extreme case where the deep JdF water goes anoxic, 40% of the water at the outfall sites comes from the surface of Haro Strait, which has a summer seasonal average oxygen concentration of 291 µmol kg-1_{. Thus, to have a mean oxygen} concentration at the outfalls below 63 µmol kg-1 _{(lower left corner, Figure 7), oxygen in} the surface waters of Haro Strait would need to fall over 130 µmol kg-1_{, an extremely} unlikely scenario. Even with a changing climate, the powerful tides in Haro Strait will maintain its high levels of oxygenation (Ianson et al., 2016; Johanneseen et al., 2014).

3.1.3 Seabed conditions at the outfalls The ROV outfall inspection provided useful information on both what species were present at the outfall site as well as the physical conditions of the seabed (Figures B1 and B2). The video clearly shows significant amounts of particulate matter being discharged from the sewage diffusers (Figure B5) and much of this particulate matter coming to rest on the seabed surrounding the outfall pipe. As the ROV approaches the seabed and the sewage outfall pipe, it can be clearly observed that many of these particles are not well consolidated and move slowly with the current. 3.1.4 Benthic community The video of the benthic community provides further biological evidence for the lack of hypoxia near the Macaulay and Clover point outfalls. Around the diffusers of the Macaulay outfall, spotted ratfish (Hydrolagus colliei) are present in high concentrations. Sunflower stars (Pycnopodia helianthoides), staghorn sculpin (Leptocottus armatus), English sole (Parophrys vetulus), and rock fish (Sebastes spp) are also present in significant numbers (Appendix B, Table B1). At the same time, no squat lobster (Munida quadrispina), slender sole (Lyopsetta exilis), or bacterial mats were observed, all of which are local species known to inhabit hypoxic environments (Chu and Tunnicliffe, 2015). Species present near Clover Point were different from those seen in such high abundance at Macaulay, but again no hypoxic-tolerant species were seen (Table B2). The difference in speciation between the two outfall sites is likely due to the rocky

bottom habitat at Clover (Figure B1) compared to the mud and gravel habitat at Macaulay (Figure B2, Castro and Huber, 2008). 3.1.5 Inorganic Carbon and Total Alkalinity DIC and TA samples (colored diamonds, Figure 8a) collected from the outfall monitoring stations (Figure 1c) overlap with previous observations from JdF and Haro Strait (Ianson et al., 2016). As the organic carbon within the sewage is respired it would be expected to raise DIC and slightly lower TA (Emerson and Hedges, 2008); however, my observations show no evidence of such a shift. This result is also consistent with the finding by Johannessen et al. (2015) that municipal sewage outfalls do not significantly contribute to the regional (SoG and JdF) carbon budget. Nearly all of the samples taken in November 2015 (blue), July 2016 (red), and January 2017 (cyan) were below a pH of 7.8 and corrosive to aragonite (ΩA < 1). The presence of corrosive near-surface waters is not unprecedented in the area, and similar observations have been made in the SoG (Ianson et al., 2016). The July 2015 (green) samples stand out in this regard as they show aragonite supersaturation throughout the water column, likely due to summer biological draw-down of DIC.