You may notice some images loading slow across the Open Collections website. Thank you for your patience as we rebuild the cache to make images load faster.

Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

The inorganic carbonate chemistry of the southern Strait of Georgia Moore-Maley, Benjamin Lee 2014

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.

Item Metadata


24-ubc_2015_february_mooremaley_benjamin.pdf [ 3.59MB ]
JSON: 24-1.0167635.json
JSON-LD: 24-1.0167635-ld.json
RDF/XML (Pretty): 24-1.0167635-rdf.xml
RDF/JSON: 24-1.0167635-rdf.json
Turtle: 24-1.0167635-turtle.txt
N-Triples: 24-1.0167635-rdf-ntriples.txt
Original Record: 24-1.0167635-source.json
Full Text

Full Text

The Inorganic Carbonate Chemistry of the SouthernStrait of GeorgiabyBenjamin Lee Moore-MaleyB.Sc. Biology, The Evergreen State College, 2007A THESIS SUBMITTED IN PARTIAL FULFILLMENTOF THE REQUIREMENTS FOR THE DEGREE OFMaster of ScienceinTHE FACULTY OF GRADUATE AND POSTDOCTORALSTUDIES(Oceanography)The University of British Columbia(Vancouver)December 2014© Benjamin Lee Moore-Maley, 2014AbstractA one-dimensional, biophysical, mixing layer model was modified to hindcast pHand aragonite saturation state (ΩA) in the southern Strait of Georgia. The modelskill in predicting spring phytoplankton bloom timing in previous studies was akey factor in its selection. Dissolved inorganic carbon (DIC) and total alkalinity(TA) were added as state variables, coupled to the existing nitrogen-based biologi-cal equations. Additional processes determined to be important to the system suchas air-sea gas exchange and nutrient-limited excess carbon uptake were also im-plemented. pH and ΩA could then be calculated from DIC and TA. Modeled DIC,TA, pH, and ΩA were evaluated against data collected between 2003 and 2012.Modeled and observed quantities agreed except in some summers, with surfacedisagreement driven primarily by plume variability and subsurface disagreementdriven primarily by model overproductivity. Model outputs demonstrated a near-surface seasonal cycle characterized by low pH and ΩA in winter and high pH andΩA in summer. In order to evaluate the sensitivity of model pH and ΩA to lo-cal forcing quantities, the model was run in one year increments over the periodfrom 2001 through 2012. For each year, each of three forcing records (wind speed,freshwater flux, cloud fraction) was shifted across all possible years during thesame period for a total of 432 experimental runs. When regressed against springwind speed, model surface pH demonstrated a clear, negative correlation. Modelspring ΩA demonstrated a negative correlation to cloud fraction. Summer pH andΩA were most sensitive to freshwater flux, both showing negative correlations.Model pH and ΩA sensitivity to freshwater TA and pH were also evaluated overthe same period using a set of realisitic freshwater chemistry scenarios determinedfrom observations in the Fraser River. Model pH and ΩA demonstrated oppositeiicorrelations to freshwater TA with sensitivities at opposite extremes of freshwaterpH. The sensitivity results identify important links between local processes and thecarbonate chemistry in the southern Strait of Georgia, and perhaps provide somesimple forecasting tools to be tested in the future.iiiPrefaceThis project was initiated by my supervisors Dr. Susan Allen and Dr. DebbyIanson. Dr. Allen previously adapted the model developed in (Collins et al., 2009)to include multiple classes of nutrients, phytoplankton, and grazers, and JeremySklad performed the preliminary tuning to the STRATOGEM dataset. Dr. Allenperformed further tuning using the genetic algorithm and developed the northernreturn parameterization, and as such co-wrote appendix 2 and wrote appendix 3.I implemented dissolved inorganic carbon and total alkalinity as state variablesand the air-sea CO2 gas exchange module under the guidance of Dr. Allen andDr. Ianson. I also implemented the carbonate system calculations module usingFortran code from Dr. Ianson and from Lewis and Wallace (1998). Additionally,I refit the model mesozooplankton and tuned the grazing preference parameters totry and reduce the persistent summer overproductivity in the model. Dr. Iansonorganized the sampling behind the majority of unpublished carbon data presentedin this thesis, although I did personally perform the sampling during the summer2011 carbon cruises. I have based some analyses in this thesis on additional datagenerously provided by Dr. Sophia Johannessen and Dr. Robie Macdonald at theInstitute of Ocean Sciences.Following the guidance of Dr. Allen and Dr. Ianson, I performed the literaturereview contained in chapter one. I evaluated the model carbonate system parame-ters against the IOS carbon dataset, and analyzed the Fraser River carbonate chem-istry data to obtain the proposed freshwater chemistry scenarios. I also developedand performed the sensitivity studies presented in chapters two and three of thisthesis and evaluated the results. I performed the majority of the runs using Pythonautomation scripts developed by Doug Latornell at UBC, and processed and plot-ivted the results using my own Matlab scripts. Chapters 2 and 3 will submitted asmanuscripts for publication with myself, Dr. Allen, and Dr. Ianson as authors.A third manuscript will also be submitted containing an analysis of the Strait ofGeorgia carbon data amassed during recent IOS cruises. Data from this pendingpublication are referenced in this thesis as: D. Ianson, manuscript in preparation,2014.vTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixGlossary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiiAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvi1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 Modeling study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112.1 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112.1.1 The SOG physical model . . . . . . . . . . . . . . . . . . 132.1.2 The SOG biology model . . . . . . . . . . . . . . . . . . 142.1.3 The carbonate chemistry model . . . . . . . . . . . . . . 182.1.4 Initialization and boundary conditions . . . . . . . . . . . 212.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 242.2.1 Distribution and seasonal variability . . . . . . . . . . . . 242.2.2 Evaluation . . . . . . . . . . . . . . . . . . . . . . . . . 262.2.3 Sensitivity studies . . . . . . . . . . . . . . . . . . . . . 35vi2.3 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 412.3.1 Mechanisms . . . . . . . . . . . . . . . . . . . . . . . . 412.3.2 Spatial distribution . . . . . . . . . . . . . . . . . . . . . 522.3.3 Comparison to other regions . . . . . . . . . . . . . . . . 552.3.4 Limitations . . . . . . . . . . . . . . . . . . . . . . . . . 573 Freshwater analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . 593.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 593.2 Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 643.2.1 Sources . . . . . . . . . . . . . . . . . . . . . . . . . . . 643.2.2 Uncertainty . . . . . . . . . . . . . . . . . . . . . . . . . 643.2.3 Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . 663.3 Sensitivity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 723.3.1 Scenarios . . . . . . . . . . . . . . . . . . . . . . . . . . 723.3.2 Model and evaluation . . . . . . . . . . . . . . . . . . . . 743.3.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . 753.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 763.4.1 Mixing . . . . . . . . . . . . . . . . . . . . . . . . . . . 763.4.2 Parameterization for model . . . . . . . . . . . . . . . . . 803.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 814 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89A Model Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103B Model Tuning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111C Northern Return Flow . . . . . . . . . . . . . . . . . . . . . . . . . . 115D Residence Time . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119E Correlation Lag . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123viiList of TablesTable 2.1 Primary datasets used for initialization and boundary condi-tions, tuning, and evaluation of the model. . . . . . . . . . . . 21Table 2.2 Initialization dates and profile sources for each model base run. 23Table 2.3 Coordinates for the model and sampling locations. . . . . . . . 27Table 2.4 CTD casts used to estimate the spatial variability of T and S(figure 2.6a and b) for the given sampling dates from the IOScarbon dataset. . . . . . . . . . . . . . . . . . . . . . . . . . . 27Table 2.5 Summary of model sensitivity experiments. . . . . . . . . . . . 36Table 3.1 Datasets used to assess the range and variability of freshwaterTA and pH near the Fraser River estuary. . . . . . . . . . . . . 65Table 3.2 SoG freshwater TA scenarios based on fits between Fraser Riverand SoG TA data and 15-day backaveraged total parameterizedfreshwater discharge, Qt. . . . . . . . . . . . . . . . . . . . . 73Table A.1 Biological model parameters. . . . . . . . . . . . . . . . . . . 104Table A.2 Tuned biological model parameters. . . . . . . . . . . . . . . . 106Table A.3 Grazer parameters, broken down by grazer class. . . . . . . . . 107Table A.4 Waste fractions from natural mortality, excretion, and sloppygrazing, broken down by waste destination. . . . . . . . . . . . 108Table D.1 Difference thresholds for T , S, NO−3 , chl a, DIC, and TA usedfor determining convergence between runs. . . . . . . . . . . . 121viiiList of FiguresFigure 2.1 Map of the Salish Sea including important waterways, localmetropolitan areas, hydrometric river gauging sites, samplingstations, and the location of the model used in this study. . . . 12Figure 2.2 Biological scheme depicting modeled state variables and pre-scribed mesozooplankton. . . . . . . . . . . . . . . . . . . . 16Figure 2.3 Conceptual diagram of electroneutral nutrient uptake in phyto-plankton. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19Figure 2.4 Boundary conditions for DIC, TA, and S. . . . . . . . . . . . 22Figure 2.5 Four years of daily observed wind speed cubed and total fresh-water discharge, and modeled mixing layer depth, S, chl a,DIC, pH, and ΩA. . . . . . . . . . . . . . . . . . . . . . . . . 25Figure 2.6 Profiles of modeled and observed T , S, and NO−3 at DFO Northand DFO South for ten carbon cruise dates. . . . . . . . . . . 28Figure 2.7 Linear regression of all DIC and TA versus S observed at DFONorth and DFO South for the ten carbon sampling dates inspring through fall. . . . . . . . . . . . . . . . . . . . . . . . 29Figure 2.8 Profiles of modeled and observed DIC, TA, pH, and ΩA at DFONorth and DFO South for ten carbon cruises in spring throughfall. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31Figure 2.9 Modeled versus observed pH for ten carbon cruises in springthrough fall. . . . . . . . . . . . . . . . . . . . . . . . . . . . 33Figure 2.10 Modeled versus observed ΩA for ten carbon cruises in springthrough fall. . . . . . . . . . . . . . . . . . . . . . . . . . . . 34ixFigure 2.11 Scatterplots of model pH3m versus observed U backaveraged30 days from 15 Apr, model pH3m versus observed Qt backav-eraged 30 days from 10 Jul, model spring ΩAsur > 1 onset ver-sus observed CF backaveraged 30 days from 6 Mar, and modelsummer ΩAsur < 1 duration versus observed Qt backaveraged30 days from 25 Jul. . . . . . . . . . . . . . . . . . . . . . . 37Figure 2.12 Timeseries of mean, standard error, and minimum and maxi-mum 30-day backaveraged pH3m, and of 30-day backaveragedpH3m correlation against 30-day backaveraged U , Qt, CF, andPPint. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39Figure 2.13 Correlation timeseries of spring ΩAsur > 1 onset, summerΩAsur <1 duration, and winter ΩAsur < 1 onset versus U , Qt, CF, andPPint. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40Figure 2.14 Timeseries of mean, standard error, and minimum and maxi-mum 30-day backaveraged ΦmixDIC, ΦfreshDIC , and ΦPPDIC, and 30-daybackaveraged ΦmixDIC, ΦfreshDIC and ΦPPDIC correlation againt 30-daybackaveraged U and Qt. . . . . . . . . . . . . . . . . . . . . 44Figure 2.15 Timeseries of modeled ΦEN and ΦGPPN from May through mid-August 2010 and 2011. . . . . . . . . . . . . . . . . . . . . . 46Figure 2.16 Spring ΩAsur > 1 onset versus spring bloom onset for all 432experimental runs, colormapped by chl a biomass. . . . . . . 47Figure 2.17 2010/2011 ratios of modeled Ca2+ and CO2−3 . . . . . . . . . . 48Figure 2.18 Winter ΩAsur < 1 onset results from sensitivity experiment 1versus observed U backaveraged 30 days from 13 Nov for allruns, backaveraged Qt < 2,500 m3 s−1, Qt < 2,500 m3 s−1 andCF < 0.75 and Qt < 2,500 m3 s−1 and CF > 0.75. . . . . . . 49Figure 3.1 Maps of the Fraser River watershed and estuary with geologicbelts, sampling sites, and the model site. . . . . . . . . . . . . 60Figure 3.2 Hydrographs showing mean, standard deviation and minimum/-maximum observed discharge for the Fraser River at Hope,BC and the Chilliwack River at Chilliwack, BC over the 1913-2013 Environment Canada daily record. . . . . . . . . . . . . 62xFigure 3.3 TA at various locations along the Fraser River and in selectFraser River tribuaries sampled at low-stage flow between 14to 27 October 2010, mid-stage flow between 28 July to 13 Au-gust 2009, and high-stage flow between 26 May to 7 June 2011. 67Figure 3.4 Timeseries of Fraser River TA at Hope, BC and the Fraser Es-tuary and Fraser River pH at the Environment Canada WaterQuality Buoy in the Fraser estuary. . . . . . . . . . . . . . . . 68Figure 3.5 TA versus S regressions from observations on several cruisesin the southern SoG and Fraser River estuary. . . . . . . . . . 70Figure 3.6 TA versus 15-day backaveraged parameterized total dischargefor Hope, BC, the Fraser estuary, and the SoG (extrapolated toS = 0). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71Figure 3.7 Model sensitivity results for nDIC3m, pH3m, and ΩAsur < 1 du-ration averaged over 1 May to 31 October and 31 October to 1May. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76Figure 3.8 Theoretical two-endmember mixing curves of DIC and TA,pH, and DIC:TA for TA f resh of 500.0 µeq kg−1 and 1300.0 µeqkg−1 at pH f resh of 7.4, 7.7, and 8.0. . . . . . . . . . . . . . . 78Figure B.1 Monthly averages of 1-3mm copepods, 3-5mm copepods, andeuphausiids from 20 years of full water-column zooplanktonsampling in the Strait of Georgia. . . . . . . . . . . . . . . . 113Figure C.1 Observed spice contours from 0 to 20 m along the Strait ofGeorgia southwestern boundary from the July 2004 STRATO-GEM cruise. . . . . . . . . . . . . . . . . . . . . . . . . . . 116Figure C.2 Contoured difference between modeled and observed T as afunction of STRATOGEM cruise and depth before and afterincluding northern flow. . . . . . . . . . . . . . . . . . . . . 117Figure D.1 Model spinup times for T , S, NO−3 , chl a, DIC, and TA usingeach available set of initialization profiles from the STRATO-GEM dataset. . . . . . . . . . . . . . . . . . . . . . . . . . . 122xiFigure E.1 Timeseries of the lag between 30-day backaveraged pH3m and30-day backaveraged U , Qt, CF, and PPint that maximizes R2,and the corresponding optimized correlation of regression be-tween the same metrics. . . . . . . . . . . . . . . . . . . . . 124Figure E.2 Timeseries of the lag between 30-day backaveraged ΦmixDIC andΦPPDIC and 30-day backaveraged U , and the lag between the 30-day backaveraged ΦfreshDIC and ΦPPDIC and 30-day backaveraged Qtthat maximizes R2. Timeseries of the corresponding optimizedcorrelation of regression between the same metrics are included. 125xiiGlossaryUBC University of British ColumbiaIOS Institute of Ocean SciencesEC Environment CanadaSoG Strait of GeorgiaVICC Vancouver Island Coastal CurrentCCS California Current SystemC-09 Collins et al. (2009)AW-13 Allen and Wolfe (2013)CO2SYS carbonate system calculation algorithm (Lewis and Wallace, 1998)STRATOGEM UBC-based 2002-2005 SoG sampling program (table 2.1)IOS water properties IOS quarterly sampling program (table 2.1)IOS carbon IOS carbon sampling program (tables 2.1 and 3.1)IOS estuary IOS low-salinity carbon sampling study (table 3.1)EC Hope EC 1979-1998 TA sampling program (table 3.1)EC estuary EC 2004-2009 TA sampling program (table 3.1)EC buoy EC 2008-present pH sampling program (table 3.1)xiiiUBC estuary UBC low-salinity 1978-79 sampling study (table 3.1)T in situ temperatureS practical salinityΩA aragonite saturation state (equation 1.1)pCO2 seawater carbon dioxide partial pressureDIC dissolved inorganic carbonTA total alkalinityCA carbonate alkalinityDIN dissolved inorganic nitrogenDON dissolved organic nitrogenPON particulate organic nitrogenPP primary productivitychl a chlorophyll aU windspeedQt freshwater flux into the southern SoGCF cloud fractionpH3m 3 m averaged pHPPint depth-integrated primary productivityΩAsur surface aragonite saturation statenDIC3m dissolved inorganic carbon normalized to TA = 2,000 µeq kg−1DIC f resh freshwater endmember dissolved inorganic carbonDICsea seawater endmember dissolved inorganic carbonxivDICmix dissolved inorganic carbon along the freshwater-seawater mixing curveTA f resh freshwater total alkalinityTAsea seawater total alkalinityTAmix total alkalinity along the freshwater-seawater mixing curvepH f resh freshwater pHpHsea seawater pHpHmix pH along the freshwater-seawater mixing curveΦmixDIC dissolved inorganic carbon vertical mixing flux to surface (depth chosen togive maximum flux)ΦfreshDIC dissolved inorganic carbon freshwater dilution fluxΦPPDIC dissolved inorganic carbon biological flux (integrated over the model do-main)ΦEN net nitrogen entrainment flux over the model domainΦGPPN gross nitrogen uptake integrated over the model domainxvAcknowledgmentsFirst and foremost I would like to thank my advisors Dr. Susan Allen and Dr.Debby Ianson for your patience and support throughout this work. You have bothbeen so genuinely invested in my progress, and your input along the way has beeninvaluable to the way I now approach research. It has truly been a pleasure towork with both of you. I would also like to thank my third committee member, Dr.Kristin Orians, for your consultations and insight not only on this written work butduring seminars as well.I would like to thank Doug Latornell for his technical support with running andmaintaining the model, Dr. Rich Pawlowicz, Dr. Mark Halverson, and Dr. OlivierRiche for their help navigating through all the data, and Dr. Kristina Brown for herinsights on carbonate chemistry early on in my project. Additionally, I am gratefulto both the Institute of Ocean Sciences and the MEOPAR network for their fundingand support, and to Dr. Sophia Johannessen and Dr. Robie Macdonald for the useof their data. I would also like to acknowledge my fellow students, particularlyJessica Spurgin, Kate LeSouef, Ben Scheifele, Tara Howatt, Chuning Wang, Jennyand Drew Snauffer, Michal Hodal, and Megan Wolfe.I would also like to thank my friends and family, my parents Belle and Larryand my sister Jane especially, thank you for your patience and support. And finallyGenna, you have shared the weight of this project with me for over two years.Thank you for being there every step of the way.xviChapter 1IntroductionHow do we translate our understanding of ocean acidification to coastal and estu-arine environments? Acidification, defined here as persistent decreasing pH of amarine system due to anthropogenic carbon accumulation, can be clearly demon-strated in the pelagic ocean where spatial and temporal variability of pH are on theorder of the mean trend. Evidence of acidification in the coastal and estuarine zonesof the world oceans is less clear due to the larger spatial and temporal variabilityencountered there. As these zones are disproportionately productive relative to thepelagic ocean and thus relied upon by human civilizations for natural resources, itis important to understand acidification in this new context.It is now well-accepted that oceanic uptake of anthropogenic atmospheric CO2decreases ocean pH (Brewer, 1978; Caldeira and Wickett, 2003; Doney et al.,2009). Measureable declines in ocean pH have been observed at Station ALOHA(Dore et al., 2009), at the Bermuda Atlantic Time series Study stations (Bates,2007), at the European station for time series at the Canary Islands (Gonza´lez-Da´vila et al., 2007), and in the north Pacific Ocean (WOCE P16N line; Byrneet al., 2010). Further mean declines of 0.2 to 0.3 for the surface ocean are pro-jected by 2100 based on IPCC projections for future CO2 emissions (Caldeira andWickett, 2005; Orr et al., 2005).Decreasing seawater pH reduces the concentration of carbonate and thus thecalcium carbonate (CaCO3) saturation state. The CaCO3 saturation state, Ω, is1defined here as the rearranged solubility equilibrium equationΩ=[Ca2+][CO2−3 ]Ksp(1.1)where Ksp is the solubility equilibrium constant. The two most common mineralcompositions of CaCO3 are calcite and aragonite; the saturation state of aragonite,ΩA, is lower than that of calcite, ΩC, as aragonite is about twice as soluble (Mucci,1983; Morse and Mackenzie, 1990). Seawater that is undersaturated with respectto aragonite (ΩA < 1) is potentially corrosive to organisms that use aragonite intheir shells and other hard parts (Feely et al., 2012a).Following the earliest research into the effects of high pCO2 on calcareous or-ganisms, a wealth of studies have emerged demonstrating a variety of biologicalimpacts. Notably, CaCO3 dissolution has been demonstrated in coccolithophores(calcite; Riebesell et al., 2000), pteropods (aragonite; Bednarsek et al., 2012; Lis-chka and Riebesell, 2012; Roger et al., 2012), and oysters (aragonite larval stage;Barton et al., 2012) during high (> 800 µatm) pCO2 seawater treatments. Thephysiology of shelled organisms is complex, and not all responses to high pCO2are negative. Increased growth and calcification rates have been demonstrated incoccolithophores (Riebesell and Tortell, 2011) and other favorable impacts havebeen observed in certain species of bivalves (Ries et al., 2009). However, the pre-dominant outcomes from studies over a large suite of organisms have been over-whelmingly negative (Kroeker et al., 2010).The primary habitats for the majority of the organisms for which negative acid-ification impacts are anticipated lie in the coastal and estuarine zones of the world’soceans. Coastal upwelling zones are also more vulnerable to low pH and ΩA ex-tremes than the open ocean. Source waters in coastal upwelling originate fromintermediate water masses where DIC accumulates naturally and is not ventilatedat the surface. The added DIC from surface anthropogenic CO2 that has been ex-ported to these water masses combined with the naturally high DIC concentrationdrives an enhanced low pH signal in upwelled waters (Feely et al., 2008). As of2004, the depth of the aragonite saturation horizon (ΩA = 1) in the temperate NorthAtlantic Ocean had not changed much since preindustrial times (Feely et al., 2004).In the North Pacific Ocean however, anthropogenic CO2 invasion combined with2decadal oscillations in strength of the North Pacific Subtropical Gyre have causedthe aragonite saturation horizon to shoal by 1 to 2 m yr−1 (Feely et al., 2012b).The reduction in the ΩA of Pacific Intermediate Water (PIW) has strong impli-cations for the California Current upwelling system. Waters undersaturated withrespect to aragonite have been observed within 50 m of the surface at multiple sitesalong the California Current (Feely et al., 2008). Central California Current watersbelow 60 m are projected to become permanently undersaturated with respect toaragonite by 2030 (Gruber et al., 2012).Estuarine habitats for marine calcifiers that are connected to the California Cur-rent upwelling are particularly vulnerable to reductions in pH and ΩA. In NetartsBay, Oregon, near-surface aragonite-undersaturated waters have been observed toflow over the estuarine sill and into the bay during upwelling events. Large-scalefluctuations in carbonate chemistry due to the inflow of upwelled waters were ob-served on the order of days (Barton et al., 2012). At the Whiskey Creek shellfishhatchery in Netarts Bay, Pacific oyster larvae demonstrated reduced growth andsurvival when exposed to DIC-rich upwelled waters (Barton et al., 2012). To thenorth in Puget Sound and Hood Canal, Washington, deep inflows of aragonite-undersaturated, low-pH waters were evident from casts taken during the 2008 up-welling season (Feely et al., 2010).The largest variations in carbonate chemistry throughout the California Currentare natural. Seasonal variations are attributed to biological processes (e. g. nutrientuptake and remineralization) and contributions from terrestrial runoff in additionto changes in upwelling strength (Hauri et al., 2013; Ianson and Allen, 2002). Inthe Southern California Bight (Point Conception to La Jolla, CA), pH and ΩAempirically-derived from T and dissolved O2 data collected during the 2008-2010CalCOFI program ranged from 7.6 to 8.15 and 0.75 to 3.0, respectively, above125 m (Alin et al., 2012). Greater shelf width and nutrient fluxes in the northernreach of the California Current enhance the role of biology in carbonate chemistryvariability. On the west coast of Vancouver Island, upwelling does not stronglypenetrate the inner shelf inshore of roughly the 100 m isobath (Thomson, 1981)and DIC fluxes there are dominated by primary productivity near the surface andremineralization in the subsurface (Bianucci et al., 2011). The presence of thenutrient-rich Vancouver Island Coastal Current and the relatively small volume of3this isolated shelf water column make DIC fluxes especially large (Ianson et al.,2003). Bianucci et al. (2011) did not find ΩA < 1 anywhere over the VancouverIsland shelf during a modeling study of the 1993 upwelling season, however, ΩAdata calculated from T and dissolved O2 demonstrate a 40% annual probability ofencountering ΩA < 1 in the top 60 m (Ianson, 2013).The inland basins that share the Vancouver Island shelf as their primary pointof connection with the Pacific Ocean bear special mention. They are collectivelynamed the Salish Sea (figure 2.1). The Salish Sea behaves as a large estuary andreceives PIW primarily through estuarine return flow sourced from the polewardCalifornia Undercurrent (Freeland and Denman, 1982), although shelf water asshallow as 75 m can also feed this inflow (Hickey et al., 1991). Like the VancouverIsland shelf, carbonate chemistry in the Salish Sea is highly variable spatially andtemporally. In the sea’s southern basins, Puget Sound and Hood Canal, pH calcu-lated from winter and summer 2008 DIC and TA observations ranged from 7.32 inthe deep waters of the nearly hypoxic Hood Canal basin to 8.25 near the surface ofthe main basin (Feely et al., 2010). Spatial variability near the surface was high,with pH as low as 7.77 over the silled region at the mouth of Puget Sound duringsummer 2008. In the Strait of Georgia, surface pH and ΩA are generally about 8.0and 1.5, respectively, although pH > 8.4 and ΩA > 3 have been observed duringblooms (D. Ianson, manuscript in preparation, 2014). In the deep Strait of Geor-gia, pH values as low as 7.48 have been reported accompanied by ΩA < 0.5. In thebrackish Fraser River plume region, ΩA can be less than 0.4 (D. Ianson, manuscriptin preparation, 2014). Feely et al. (2010) calculated the pH change due to anthro-pogenic DIC in Puget Sound to be as high as 0.11 in surface and as high as 0.06 inthe deep basin. Though these estimates are similar to mean estimates for the openocean, they are small relative to the local natural variability.The Salish Sea can be separated into three main basins: Puget Sound to thesouth, the Strait of Juan de Fuca to the west, and the Strait of Georgia to the north(figure 2.1). For this study, I will focus on the largest of the three basins, the Straitof Georgia (SoG hereinafter). The SoG is semi-enclosed with a primary route ofexchange with the Pacific Ocean to the south through a network of passages, thelargest of which is Haro Strait. Exchange to the north through Johnstone Strait isgenerally regarded as small (LeBlond, 1983; Pawlowicz et al., 2007) but a mod-4eling study by Crean et al. (1988) suggests that surface layer outflow to the northaccounts for up to 17% of the SoG freshwater budget. The SoG is ∼25 km acrossand 220 km along the thalweg. The central basin is ∼400 m deep and significantlyshallower toward Haro Strait.It is important to consider the SoG in the context of the other basins. The SoGis the most stratified of the three basins due to its proximity to the Fraser River out-flow, which drains ∼238,000 km2 of temperate British Columbia, Canada. FraserRiver discharge is snowmelt-dominated, characterized by low winter outflow andhigh summer outflow with a peak summer freshet that ranges from ∼5,000 to12,000 m3 s−1 (Environment Canada data, The Fraser River is the largest freshwater contribution to the SoG, particu-larly during the summer (Waldichuk, 1957). Rainfall-dominated rivers do makelarge freshwater contributions during winter months.The rivers that feed the SoG drive an estuarine circulation causing a relativelyfresh upper layer outflow and deep return flow to pass through the Strait of Juande Fuca (Masson and Cummins, 2004). Brackish estuarine outflows exiting theSoG pass over restrictive sills to the north and south of Haro Strait (figure 2.1).These sills are sites of enhanced tidal mixing such that deep estuarine inflow ismodified with upper layer estuarine outflow before entering the SoG (Pawlowiczet al., 2007). Similar estuarine circulation and tidal mixing occur over sills inAdmiralty Inlet between the Strait of Juan de Fuca and Puget Sound (Geyer andCannon, 1982). As a result of tidal mixing, the Strait of Juan de Fuca is lessstratified, less productive, and higher in surface nutrients than the SoG or PugetSound (Masson and Pen˜a, 2009).Estuarine circulation combined with the silled bathymetry and tidal mixing inHaro Strait create a three-layer structure in the SoG (Pawlowicz et al., 2007). Thesurface layer is a brackish mixture of river discharge and entrained seawater fromthe intermediate layer. Pawlowicz et al. (2007) define the depth range of this layeras the surface to 50 m. Below 50 m, the intermediate layer is created by estuarinereturn flow through Juan de Fuca Strait that mixes with brackish surface water inHaro Strait. The density of the estuarine inflow depends on the properties of bothinflowing and outflowing water, and on the intensity of mixing over the sills. Gen-erally, the density after mixing causes this intrusion to penetrate slightly below the5sill depth (∼150 m). The deep layer (>200 m, Pawlowicz et al., 2007) is isolatedfrom exchange beyond the SoG for most of the year. Only when intrusion densitiesin Haro Strait are high enough can they penetrate the deep layer. Renewal eventsgenerally take place twice a year in pulses during the neap tides when mixing overthe sills is weakest. Late spring renewal is characterized by inflow of relativelycold, fresh, nutrient-poor shelf water and late summer renewal is characterized byrelatively warm, salty, nutrient rich PIW (Masson, 2002).Wind mixing is a strong driver of exchange between the surface and intermedi-ate water masses in the SoG. Wind patterns are driven primarily by the prevailingoceanic winds associated with the North Pacific High and Aleutian Low pressurecenters over the northeast Pacific Ocean (Thomson, 1981). Winter winds blowtypically from the southeast (Aleutian Low circulation) and summer winds are typ-ically from the northwest (North Pacific High circulation), although complicatedtopography and differential heating between land and sea frequently produce local-ized winds that deviate from this idealized behavior (Thomson, 1981). Regardless,windspeed is generally weaker in the summer than in the winter, following thespring transition between prevailing wind patterns (Collins et al., 2009).The SoG is a productive ecosystem (∼280 g C m−2 yr−1, Harrison et al., 1983;∼ 223± 53 g C m−2 yr−1, Riche, 2011), and estuarine circulation and wind mix-ing are coupled to the local ecosystem dynamics. Phytoplankton in the SoG aregenerally either light-limited or N-limited with respect to growth (Harrison et al.,1983). While nitrate is supplied year round via estuarine return flow to the inter-mediate layer (Mackas and Harrison, 1997), light-limitation persists during winterdue to poor solar irradiance and wind-induced mixed layer deepening. In springwhen wind relaxes and solar irradiance increases, freshwater stratification allowsphytoplankton to remain in the euphotic zone initiating a characteristic spring di-atom bloom. Following the bloom, nitrate is nearly depleted from the surface andweaker, N-limited growth of flagellates and other small phytoplankton persists un-til fall when increased winds once again mix nitrate-rich water to the surface andfall diatoms can bloom (Harrison et al., 1983). Despite surface N-limitation, itis still estimated that approximately half of annual production occurs in summer(Riche, 2011).Consumption of DIC by phytoplankton reduces pCO2 near the surface, as6demonstrated for the west coast of Vancouver Island (Ianson et al., 2003). In theSoG, Evans et al. (2012) observed surface pCO2 throughout the basin to be wellbelow the atmospheric concentration during July and August 2010. Observationsfor pCO2 in October/November 2008 and in February/March 2009 were all abovethe atmospheric concentration with the exception of February/March 2009 near theFraser River plume. A model hindcast of the 2009 spring diatom bloom date sug-gests the bloom did not occur until after April 1 (Allen and Wolfe, 2013), howeversatellite observations found evidence of a bloom as early as late February (Goweret al., 2013). The low spring 2009 pCO2 observations near the Fraser River plumemay be indicative of the early productivity suggested in these satellite chlorophyllobservations.The roles of wind and freshwater in the seasonal productivity of the SoG arecomplicated by competing mechanisms. Although wind mixing brings nitrate-richwater to the surface, it also deepens the mixing layer (St. John et al., 1993). Thussome period of stratification is required to keep phytoplankton in the euphotic zone.A modeling study by St. John et al. (1993) found that wind had little effect on ni-trate entrainment during peak river discharge due to the presence of the freshwaterlens. However, Yin et al. (1997a) found that wind mixing and estuarine entrain-ment both increased primary productivity in August 1991. The timing of the springbloom is sensitive to wind and cloud fraction because of light limitation (Collinset al., 2009; Yin et al., 1997b), however, spring productivity is further complicatedby the seasonal zooplankton cycle. There is an ontogenic zooplankton migration ofprimarily Neocalanus plumcrus copepods in spring that is independent of bloomtiming (Harrison et al., 1983). Late blooms can be rapidly attenuated by zoo-plankton grazing (Yin et al., 1997b), and copepod nauplii may miss early bloomsaltogether (El-Sabaawi et al., 2009).The Fraser River plume can have an impact on the distribution of phytoplank-ton. St. John et al. (1993) demonstrated spring blooms to be stronger near theplume. However, using fluorescence measurements from an instrumented ferry,Halverson and Pawlowicz (2013) found integrated chlorophyll to be much lowerin the plume than outside, and the depth of the chlorophyll maximum was muchshallower in the plume. They suggested light limitation due to river turbidity asa possible mechanism. Stockner et al. (1979) found extinction coefficients to be73-fold higher in the summer plume than in winter during high-salinity conditions.In addition, a modeling study in Howe Sound found that turbidity in the SquamishRiver plume reduced PP by 35% (Stockner et al., 1977).Tools used in previous ecosystem modeling studies in the Salish Sea rangefrom small scale 1-dimensional (1-D) near-surface biophysical coupled models(Collins et al., 2009; Allen and Wolfe, 2013) to mesoscale 3-dimensional (3-D)biogeochemical-physical coupled models (Khangaonkar et al., 2012; Giddings et al.,2014). Since the early modeling work of Crean et al. (1988) and Foreman et al.(1995), a wealth of 3-D modeling tools have been implemented for the regionincluding POM (Masson and Cummins, 2004), ROMS (Sutherland et al., 2011;Snauffer et al., 2014), FVCOM (Khangaonkar et al., 2012) and NEMO (N. Soon-tiens, manuscript in preparation, 2014). While 3-D models generally capture thespatial extent of the system, their computational expense limits grid spacing to rel-atively coarse resolution. Alternatively, 1-D models allow for fine grid spacing dueto their small domain. Collins et al. (2009) and Allen and Wolfe (2013) demon-strated the importance of fine vertical grid spacing in accurately predicting springdiatom bloom timing in the southern SoG.1-D ocean boundary layer models differ primarily in their parametrizations ofturbulent mixing, the most common being the K profile parametrization (KPP) firstorder turbulence closure scheme of Large et al. (1994), the Mellor Yamada 2.5(MY2.5) second order turbulence closure scheme of Mellor and Yamada (1982),and the generic length scale (GLS) scheme of Warner et al. (2005); bulk mixingschemes (e. g. Price et al., 1986) also provide a simpler alternative. Several 1-D vertical ecosystem models using compartmentalized nutrients, phytoplankton,zooplankton and detritus (NPZ and NPZD) have been developed for the northeastPacific Ocean at Ocean Weather Station Papa using the MY2.5 (Kawamiya et al.,1995; Denman and Pen˜a, 2002) and KPP (Jeffery, 2002) mixing parametrizations.The KPP scheme is generally considered to simulate stronger and thus more real-istic mixing than the MY2.5 scheme (Large, 1998), and the KPP model of Jeffery(2002) is the template for the SOG model used in Collins et al. (2009).While carbonate system models were used primarily to estimate carbon fluxesin earlier studies (e. g. Ianson and Allen, 2002) a wide range of carbonate systemmodels have been developed recently for the purpose of predicting pH and ΩA,8particularly in the Baltic and North Seas (Kuznetsov and Neumann, 2013; Artioliet al., 2012; Gypens et al., 2011) and in the California Current Ecosystem (Hauriet al., 2013; S. Siedlecki, manuscript in preparation, 2014). Carbonate systemmodels are generally based on DIC and TA equations coupled to a biophysicalmodel, however statistically-derived empirical models for pH and ΩA from tracerssuch as T and O2 have proved useful in the California Current as well (Lee et al.,2006; Juranek et al., 2009; Alin et al., 2012; Lara Espinosa, 2012), although suchmodels are generally not well-suited for estuaries (Reum et al., 2014).This thesis aims to address the seasonal variability of pH and ΩA in the southernSoG in the context of local scale physical forcing. Previous studies (e. g. Collinset al., 2009; Allen and Wolfe, 2013; Yin et al., 1997a,b) have established that localforcing such as windspeed, cloud fraction, and freshwater input can impact theSoG ecosystem biology directly (e. g. light availability) and indirectly throughphysical processes like mixing, stratification, and estuarine circulation. As manyof the physical and biological fluxes of DIC and TA are inherently linked to theseecosystem processes, I will show that pH and ΩA in the near-surface southern SoGare strongly influenced by local forcing as well.I begin by building on the modeling work of Collins et al. (2009). The model is1-D vertical and is implemented for a station just outside of the Fraser River plumein the central southern SoG. The 1-D domain allows for fine vertical grid spacingat a low computational expense. The fine grid spacing makes the model a skilledpredictor of mixing layer depth, and the low computational expense is ideal forsensitivity studies where high numbers of model runs are required. Mixing layerdepth is an important physical feature of this study because it is directly influencedby wind and freshwater discharge, and it is so closely connected to phytoplanktongrowth and vertical mixing.Near-surface chlorophyll maxima and surface-intensified freshwater fluxes makethe near-surface layer the most highly variable over time with respect to carbonatechemistry. DIC and TA collected since 2010 in the southern SoG do not demon-strate high variability over time below 40 m (D. Ianson, manuscript in preparation,2014). Similarly, pH below about 40 m in the main basin of Puget Sound wassimilar between February and August 2008 (Feely et al., 2010). A study of locally-forced variability in the carbonate system would ideally focus on the near-surface9layer of the SoG down below the euphotic zone and the maximum mixing layerdepth. Thus, rather than resolve the entire water column to the basin’s sediments, Ichoose only the top 40 m at the model site as the domain and consider deeper watercolumn properties in terms of their 40 m boundary conditions and the parameter-ized estuarine circulation.In the next chapter I will describe in detail the implementation of the modelincluding the unpublished modifications that were made by Dr. Susan Allen inorder to resolve ecosystem dynamics following the spring bloom, and the carbonatechemistry additions that I made in order to model pH and ΩA. I will also describethe evaluation of the model using the existing data for the southern SoG (D. Ianson,manuscript in preparation, 2014), and I will present high-resolution model outputsof pH and ΩA for the top 40 m. In order to address the primary research question,I will describe and present results from sensitivity experiments of near-surface pHand ΩA to local windspeed, cloud fraction, and freshwater input. I will discussthese results in the context of seasons, spatial variability, and the oceanography ofneighboring basins.In chapter three I will describe the freshwater chemistry of the southern SoG,beginning with a discussion of the current literature and several unpublished datasetsin the Fraser River. I will then present a set of plausible freshwater TA and pH sce-narios that span the range of observed river values for both quantities. Next I willdescribe and present results from sensitivity experiments similar to those in chaptertwo, but for the sensitivity of near-surface pH and ΩA to freshwater TA and pH. Iwill discuss the outcomes of these experiments in the context of theoretical estuar-ine mixing of DIC, TA, and pH, and following this discussion I will propose onefreshwater pH and TA scenario as a suitable selection for the model described inchapter two.10Chapter 2Modeling study2.1 MethodsI employ a 1-dimensional (1-D) biophysical coupled near-surface (0 to 40 m) mix-ing model in order to investigate the behavior of pH and aragonite saturation state(ΩA) in the southern Strait of Georgia (SoG; figure 2.1a). Earlier versions of themodel have been used to successfully hindcast the spring diatom bloom in the SoG(Collins et al., 2009, hereinafter C-09; Allen and Wolfe, 2013, hereinafter AW-13).The predictive skill of the model stems from three critical attributes: fine verti-cal grid-spacing (0.5 m), high-resolution local forcing data (hourly to daily), andrigorous tuning to the Strait of Georgia Ecosystem Monitoring Program (STRATO-GEM) dataset (Pawlowicz et al., 2007; Riche, 2011; table 2.1). The STRATOGEMdataset is comprised of 48 multiparameter cruises in the southern SoG that span theperiod between April 2002 and June 2005.For the purposes of this study, I have adapted the C-09 model as modified byAW-13. Part of this adaptation includes the addition of new nutrient, photosynthe-sizer, and grazer classes in order to resolve the community dynamics following thespring bloom season. The central piece of this adaptation involves the coupling ofdissolved inorganic carbon (DIC) and total alkalinity (TA) as state variables to theexisting physical and biological variables, and implementing new parametrizationsnecessary to model these quantities (e. g. air-sea gas exchange).11−150−150−150−150 30’  125oW  30’  124oW  30’  123oW  30’  122oW   48oN  30’   49oN  30’   50oN 020406080 kmFraser RiverVictoriaParksville Hope  Strait of GeorgiaStrait of Juan de FucaDepth (m)−450−400−350−300−250−200−150−100−500 165oW  150oW  135oW  120oW   40oN   48oN   56oN CanadaUSAPacificOcean−150−150−150 124oW  45’  30’  15’  123oW  50’   49oN  10’  20’  30’ 010203040 kmVancouverNHaro StraitPuget SoundJohnstoneStraitHoweSoundModel siteSandheadsDFO NorthDFO 3 DFO 4Model siteDFO 5DFO SouthDFO 7DFO 8DFO 9DFO 1YVRWaterQualityBuoyTexadaIslandFigure 2.1: (upper) Map of the Salish Sea including important waterways,local metropolitan areas, hydrometric river gauging sites, and the loca-tion of the model used in this study (red square). Depths are contouredin blue (m). (lower) Locations of sampling stations. Sandheads windstation, YVR airport, and the Fraser River water quality buoy (greentriangles) are maintained by Environment Canada. DFO 1-9 (black cir-cles) are locations of CTD casts performed regularly by Fisheries andOceans Canada (DFO North and DFO South included; Masson, 2006).Additionally, DFO North and DFO South are sampling sites for DICand TA (D. Ianson, manuscript in preparation, 2014).122.1.1 The SOG physical modelThe physical model uses a K-profile parametrized (KPP) vertical diffusion schemeadapted from Large et al. (1994). Baroclinic pressure gradients in the SoG are sim-ulated by prescribing shoalings (depressions) in isopycnals due to the divergence(convergence) that results from horizontal transport in an enclosed basin (C-09).Three-dimensional (3-D) estuarine circulation is parameterized as an upward en-trainment velocity and outward advective flux, both defined in terms of the totalfreshwater discharge (Qt) from the Fraser River and other small rivers (C-09). Theadvective loss arises due to water column convergence, since the upward entrain-ment velocity increases with depth. The parameterization of estuarine circulationfrom 3-D to 1-D was done using Knudsen relations (C-09). In the region of themodel site (figure 2.1a), the strongest horizontal gradients are cross-strait due tothe movement of the Fraser River plume away from the river mouth (Halversonand Pawlowicz, 2008). Since cross-strait advection is primarily driven by estu-arine circulation while wind and tidal currents are mainly along-strait (LeBlond,1983), wind and tidal advection are not explicitly modeled. Net transport by windand tides is also low since they both frequently reverse (LeBlond, 1983).A realistic model of the ocean surface boundary layer is important in determin-ing primary productivity (PP) and phytoplankton abundance. The Critical DepthHypothesis (Sverdrup, 1953) posits that blooms occur when the turbulent bound-ary layer depth shoals to a critical depth where phytoplankton light-limited growthexceeds losses. More recently, Behrenfeld (2010) has suggested that blooms areinstead initiated upon cessation of boundary layer deepening which allows phyto-plankton growth to exceed losses, briefly decoupling phytoplankton growth fromzooplankton grazing. The surface boundary layer is generally regarded as themixed layer and quantified using a vertical density gradient. However, the verticaldensity gradient criteria frequently fails to predict the depth of turbulent mixing. Inthe open ocean, the diurnal heating cycle produces a dynamic layer of active mix-ing that is shallower than the density-defined mixed layer (Brainerd and Gregg,1995). In the SoG, a mixed layer is often undetectable due to continuous restrat-ification from freshwater inputs, however a dynamic mixing layer is well-defined(C-09). The present study uses a bulk Richardson number criterion to calculate the13mixing layer depth (Large et al., 1994). I use the term “mixing layer” rather than“mixed layer” throughout this work to emphasize the distinction.I explicitly model in situ temperature (t; ITS-90, Preston-Thomas, 1990) andpractical salinity (Sp; PSS-78, UNESCO, 1981) as state variables. A new ther-modynamic equation of state has been developed that defines the thermodynamicproperties of seawater in terms of conservative temperature and absolute salinity(TEOS-10; IOC et al., 2010). However, adaptation to the new standard will taketime in the literature, and the preferred carbonate system constants are currentlystill parameterized in terms of t and Sp. Since I make extensive use of carbonatesystem parameterizations throughout this study, I present the model results and anycomplimentary datasets in terms of t (hereinafter T ) and Sp (hereinafter S) for thebenefit of the reader.2.1.2 The SOG biology modelThe biological model (figure 2.2) is an adaptation of the nutrient phytoplanktonzooplankton detritus (NPZD) scheme developed to predict spring phytoplanktonbloom timing in the southern SoG (C-09). Ecosystem dynamics following thebloom were not well constrained in the original implementation. Since I am inter-ested in annual cycles, I use the model of C-09 with updates to the cloud model,freshwater flux parametrization, 40 m T boundary conditions, and river properties(AW-13) and the following additions to the biological model.The C-09 model used a single class of phytoplankton based on the centricdiatom Thalassiosira nordenskioldii, and was tuned to best reproduce the springbloom timing. However, centric diatoms tend to dominate in spring and fall onlywhile flagellates, primarily nanoflagellates (< 20 µm), are particularly abundantduring the summer (Harrison et al., 1983). The heterotrophic ciliate Myrionectarubra is also a strong contributor to photosynthetic biomass. Although not a phy-toplanker, M. rubra retains functional chloroplasts during grazing and uses themto perform photosynthesis. The ciliate is prevalent in the SoG throughout the year,particularly in summer (Harrison et al., 1983).In the present study, the diatom class has been modified from Thalassiosiranordenskioldii to represent a wider range of diatoms, and the other two dominant14primary producers, M. rubra, and nanoflagellates, have been added as state vari-ables to the model (figure 2.2). Despite their abundance on the west coast of Van-couver Island (Pen˜a et al., 1999), calcifying phytoplankton (e. g. coccolithophores)are absent from satellite observations in the Strait of Georgia (J. Gower, personalcommunication, 2014) and are thus not represented in the model. Because ammo-nium (NH+4 ) uptake is important on an annual scale in the southern SoG (Mackasand Harrison, 1997), I model NH+4 as an additional nitrogen (N) source to the sin-gle nutrient of the C-09 model, nitrate (NO−3 ). Despite using a single diatom fortheir primary producer class, C-09 did not explicitly model dissolved silica (Si) asa potentially limiting nutrient. Si loading is high in the Fraser River relative to theSoG (Swain, 2007), and terrestrial runoff provides an important source of Si fordiatoms in the SoG. However, near-limiting Si values were observed in the surfaceSoG in spring during STRATOGEM (Riche, 2011). I therefore model Si as statevariable in order to allow for Si limitation (figure 2.2).I have refit the prescribed mesozooplankton of C-09 to a new dataset (Mackaset al., 2013), and have added microzooplankton as an explicitly-modeled state vari-able. Waste quantities have been diverted from grazing, natural mortality, and ex-cretion to new state variables for dissolved organic N (DON), particulate organic N(PON), and biogenic silica detritus. Remineralization of organic matter becomes asignificant source of nutrients over seasonal timescales.Nutrient uptakeThe biological model is N-based and all biological quantities are expressed in µmolN L−1 except Si-containing compounds, and DIC and TA. Uptake of nutrients byphytoplankton is dependent on light and nutrient availability, the least abundantbeing the limiting quantity and determining the uptake rate kinetics (Liebig’s Law).Light-limited uptake is parameterized in terms of the incoming photosyntheticallyactive radiation (PAR) to give an adapted Steele’s scheme (Steele, 1962) with anextended high light range suitable for multiple organisms (equation A.1).N-limited uptake in the presence of a single N-based nutrient can be modeledby the Michaelis-Menten (MM) or Monod equation. To determine uptake ratesfor multiple N sources, I use the substitutable model of O’Neill et al. (1989, in15diatoms M. rubra flagellatesmicrozooplanktonmesozooplanktonDONPONH2SiO4bSiNO3– NH4+Figure 2.2: Biological scheme. All variables are modeled using state equa-tions (appendix A) except for mesozooplankton, which is prescribedusing a Gaussian distribution. Dotted arrows represent sinking losses.Jeffery, 2002; equations A.2a and A.2b). This model reduces to MM kinetics inthe absence of a second N source. If N is not limiting, NO−3 and NH+4 uptake mustbe assigned based on the total light (or Si) limited N uptake. Following Jeffery(2002), NH+4 is used first until the uptake capacity (equation A.4a) is reached,and then NO−3 uptake fulfills the remaining N demand (equation A.4b). Si-limiteduptake in modeled diatoms follows MM kinetics (equation A.3).Photosynthesis, grazing and wasteNutrient uptake rates are calculated separately for each of the model’s three classesof photosynthesizers: the centric diatoms, the nanoflagellates, and M. rubra. Thegoverning equations that describe the concentrations of each class follow the con-ceptual model:∂P∂ t = growth−mortality−grazing+ physical transport (2.1)16where growth, mortality and grazing are proportional to the concentration of pho-tosynthesizer P, and the constant of proportionality for growth is defined as thenutrient uptake rate discussed above (equation A.6). A sinking loss term is in-cluded for diatoms and a positive grazing term is added to represent heterotrophyin M. rubra (figure 2.2).Grazing is modeled using a Holling type II response (MM kinetics; Fasham,1995), and is calculated separately for each photosynthesizer and for each of themodel’s three grazers: M. rubra, microzooplankton, and mesozooplankton (equa-tions A.7a and A.7b). The governing equation for microzooplankton is based onequation 2.1 where P is the concentration of microzooplankton and the proportion-ality constant for growth is defined as the microzooplankton grazing rate (equa-tion A.8).The biology loop is closed by mesozooplankton. Average mesozooplanktonbiomass is fit to a summer bloom Gaussian model (equation A.9) using a 20-yeartimeseries of zooplankton sampling in the SoG (Mackas et al., 2013). Mesozoo-plankton are assumed to choose their depth based on the profile of their prey, sotheir vertical distribution in the model is parametrized in terms of the total abun-dance of prey, including PON (figure 2.2), and the average mesozooplankton con-centration (equation A.10).Waste N from natural mortality, grazing, and excretion from all phytoplanktonand zooplankton flows to pools of PON and DON (figure 2.2). Natural mortal-ity and excretion rates are fixed and fluxes from both processes are proportionalto biomass (equations A.11a, A.11b, and A.11c). Grazing waste is the result ofsloppy feeding and is proportional to grazing (equation A.12). PON and DON areremineralized to NH+4 at fixed rates (equation A.13a), and NH+4 can be used as anN source for phytoplankton or remineralized to NO−3 (equation A.13b), completingthe biological loop.Nutrient uptake, grazing, and natural mortality rates are temperature depen-dent. Implementation of the temperature dependence term based on Eppley (1972)is described in C- The carbonate chemistry modelCarbonate system parametersThe carbonate system quantities DIC and TA are modeled as state variables. DICis treated as a nutrient and a waste product, and I assume that nutrient uptake inprimary producers and remineralization from organic waste follows the RedfieldC:N:P ratio (Redfield et al., 1963) except under nutrient limitation where I explic-itly model excess DIC uptake. I allow uptake and remineralization of charged par-ticles to modify TA by assuming that the two processes occur such that electroneu-trality is maintained within the organism (figure 2.3; Zeebe and Wolf-Gladrow,2001; Brewer and Goldman, 1976).∂DIC∂ t =−(UN +UPC−RN)RC:N + physical transport (2.2a)∂TA∂ t =UNO−3−2RNO−3 −UNH+4 +RNH+4 +(UPO4−RPO4)RN:P+ physical transport (2.2b)where U represents uptake, R represents remineralization, PO4 represents dissolvedorthophosphate (HPO2−4 and PO3−4 ), and RX:Y is the Redfield ratio of element Xto element Y (equations A.14a and A.14b). UPC represents DIC uptake that is un-coupled to uptake of DIN under conditions of nutrient limitation (equation A.15).In the absence of available DIN, carbon uptake does not contribute to chlorophyllbiomass and instead passes directly to the DOC pool. The result is that the C:chlratio increases under nutrient limitation while the phytoplankton biomass C:N ratioremains constant (Ianson and Allen, 2002).DIC and TA can be mixed and advected like any other model scalar. An addi-tional pathway for DIC exists through air-sea gas exchange. Transfer of CO2 acrossthe air-sea interface in the surface grid cell is parametrized according to Fick’s sec-ond law of diffusion (Sarmiento and Gruber, 2006), using the transfer coefficientparametrization of Nightingale et al. (2000), the Schmidt number of Wanninkhof(1992), and the solubility coefficient of Weiss (1974).Model pH and ΩA are calculated from model DIC, TA, Si, T , S, and pressureusing the CO2SYS program (Lewis and Wallace, 1998). PO4 is roughly approxi-18mated from model NO−3 using the Redfield ratio of N:P. While not a robust assump-tion, differences in carbonate chemistry between the NO−3 -coupled PO4 scenarioand a PO4 = 0 scenario were negligible. The K1 and K2 equilibrium constants ofLueker et al. (2000) are recommended for seawater carbonate system calculations(i. e. Dickson et al., 2007). Since low S values (S < 19) are frequently encounteredin the SoG and in the model (C-09), I use the K1 and K2 equilibrium constants ofMillero (2010) which are suitable for a wide range of S (1 to 50). Differences be-tween calculations of model pH and ΩA using the Lueker et al. (2000) and Millero(2010) constants are negligible where S > 19.  NO3–  (0 eq)NH4+  (0 eq)Cell Seawatermesozsplanszkp 1 eq mol–1t????e??lanszkp–1 eq mol–1? ?? ?n?zsplanszkp 1 eq mol–1 uptake ?e?eᔅsplanszkp????lanszkpNo Effect ??m llll m? ?? 1 eq mol–1??? llll ?? ? ? ?–1 eq mol–1mesoe?eᔅse??–2 eq mol–1DONanszkp ?pంฆoz?e∅seഎDOPOH – 2  OH – OH –H+ 3  OH – OH –HPO42–  (1 eq)PO43–  (2 eq)SiO(OH)3–  (1 eq)HCO3–  (1 eq)NH4+ (0 eq) OH –H+NH4+ (0 eq) NO3– (0 eq)H+PO43– (2 eq)3   H+Figure 2.3: Conceptual diagram of electroneutral nutrient uptake in phyto-plankton. H+ and OH− contribute -1 and 1 eq mol−1 to TA, respectively.DON is dissolved organic nitrogen and DOP is dissolved organic phos-phate; neither contributes to TA. Based on Zeebe and Wolf-Gladrow(2001).19Model river chemistryI assigned the model freshwater carbonate chemistry parameters following a de-tailed analysis of datasets from Environment Canada sampling programs in theFraser River estuary and ∼ 150 km upstream at Hope, BC (, and a set of ten TA casts,five from de Mora (1981) and five from D. Ianson (manuscript in preparation,2014), in the southern SoG near the model site (figure 2.1; see chapter 3). Theestuary TA data are consistently lower than data from the Hope dataset, but higherthan extrapolated endmembers (S = 0) from the ten SoG casts. It is not surpris-ing that TA varies between the three sampling locations. Fraser Valley agriculturaland urban wastewater runoff enter the Fraser River downstream of Hope (Mackasand Harrison, 1997) and may affect TA. Suspended particles affecting TA may alsoprecipitate and sink on entering the SoG (de Mora, 1983).The model is sensitive to freshwater TA, and so the selection of a suitablemodel freshwater TA value is important. However, in addition to the wide rangeof river TA between sampling, this selection is complicated by the presence oforganic matter. Organics are present throughout estuarine systems and can intro-duce error into the approximation of carbonate alkalinity from TA (Hunt et al.,2011). In chapter 3, I investigate the model sensivity to six freshwater TA sce-narios based on the datasets at Hope, in the Fraser estuary, and in the southernSoG in order to support my selection. I define the model freshwater TA as themean value (873.0 µmol kg−1) from the 2004 through 2008 Fraser River estuarydataset. I chose the mean Fraser estuary value because the data lie in betweenthe other datasets, and because the estuary sampling sites are closest to the FraserRiver mouth. I did not consider TA in other rivers due to the lack of available data.Waldichuk (1957) estimated the Fraser River to contribute ∼77% of the annualtotal freshwater runoff to the SoG based on observations for 1950, however thisestimate was much lower (∼55 to 60%) from December through April. Still, theFraser River is the largest freshwater source to the SoG year-round.Freshwater DIC and TA are assumed to be coupled in the carbonate system tothe mean observed Fraser River pH (7.77) determined from hourly pH observationsbetween 2008 and 2013 at Environment Canada’s Fraser River Water Quality Buoy20(; fig-ure 2.1b). The mean buoy pH is used along with the model freshwater TA, Tand S to calculate the model freshwater DIC (CO2SYS) using the K1 and K2 con-stants of Millero (1979). The Fraser River datasets and model freshwater carbonatechemistry values are discussed in detail in chapter Initialization and boundary conditionsI make use of three datasets for the initialization, 40 m boundary conditions, tun-ing, and evaluation of the model. The STRATOGEM dataset (Pawlowicz et al.,2007; Riche, 2011; table 2.1) contains approximately monthly profiles of T , S,fluorescence, chl a, NO−3 , and Si taken at the model site between April 2002 andJune 2005. For years beyond the STRATOGEM dataset, I use similar data col-lected at nearby locations (DFO North, figure 2.1) during a quarterly samplingprogram maintained by the Institute of Ocean Sciences (IOS) in Sidney, BC (IOSwater properties dataset; Masson, 2006; Masson and Pen˜a, 2009; D. Masson, un-published data, 2014; table 2.1). Where I require carbonate chemistry data, I useDIC, TA, NO−3 , PO4, Si, T , and S from ten IOS carbon cruises in the southern SoG(IOS carbon dataset; D. Ianson, manuscript in preparation, 2014; table 2.1).Dataset Source(s)STRATOGEM Pawlowicz et al., 2007; Riche, 2011IOS water properties Masson, 2006; Masson and Pen˜a, 2009; D. Masson, pers. com., 2014IOS carbon D. Ianson, manuscript in preparation, 2014Table 2.1: Primary datasets used for initialization and boundary conditions(section 2.1.4), tuning (appendix B), and evaluation (section 2.2.2) of themodel.40 m boundariesModel T , S, and chl a fluorescence at 40 m were fit to seasonal curves (C-09; AW-13) and model NO−3 and Si at 40 m were assigned constant values (C-09) usingthe STRATOGEM dataset. Model NH+4 at 40 m was set to 1.0 µM to reflect a lowbackground concentration. Dilution is a strong driver of TA in the SoG, and the21variation of TA with S at 40 m would ideally be present in the boundary conditions.However, the sampling frequency within the IOS carbon dataset is too low to reflectthis seasonal variation. Instead, I calculated 40 m TA from 40 m S (figure 2.4) usinga linear regression of TA and S data from the IOS carbon dataset. I observed the40 m pH from the IOS carbon dataset to be approximately constant, and thus usedthe mean value to calculate the 40 m DIC boundary condition (figure 2.4) alongwith the seasonal fits of TA, T , and S and the mean 40 m values of PO4 and Si(CO2SYS) using the Millero (2010) K1 and K2 constants.J F M A M J J A S O N D J204020502060207020802090210021102120DIC [µmol kg−1] or TA [µeq kg−1]  29.329.429.529.629.729.829.930S p [PSS−78]DICTASpFigure 2.4: Boundary conditions for DIC (solid), TA (dash), and S (thick).InitializationModel runs are initialized in fall on the date of available observations (table 2.2)and run through a full year and then beyond until the following December 31st.Fall was selected to allow adequate spinup time before late winter dynamics be-come important for the spring bloom. State variables are initialized using profilesof T , S, chl a, NO−3 , and Si from either the STRATOGEM dataset or the IOS waterproperties dataset. Phytoplankton biomass is converted from chl a using a NO−3 -to-chl a conversion factor of 1.6 determined from the STRATOGEM dataset. NH+4is initialized to the 40 m boundary value and set to zero in the mixing layer. Micro-22zooplankton biomass is assumed to be equal to the nanophytoplankton biomass.PON and biogenic Si detritus are estimated to be 0.2 of the microphytoplanktonbiomass and DON is 0.1 of PON. Semi-labile DON is generally higher than PONin the SoG (Johannessen et al., 2008), however model DON is not sensitive tothe initialization value and exceeds model PON after spin-up. IOS carbon profilesfrom the 11 Sept, 2011 cast date are used to initialize DIC and TA. Runs with DICand TA initialized using each of the ten IOS carbon cruises converged within 30days.Base year Initialization date Source2001 31 Aug 2000 Masson, 20062002 19 Sep 2001 Masson and Pen˜a, 20092003 8 Oct 2002 Pawlowicz et al., 2007; Riche, 20112004 9 Oct 2003 Pawlowicz et al., 2007; Riche, 20112005 19 Oct 2004 Pawlowicz et al., 2007; Riche, 20112006 15 Sep 2005 Masson and Pen˜a, 20092007 14 Sep 2006 Masson and Pen˜a, 20092008 4 Oct 2007 Masson and Pen˜a, 20092009 11 Sep 2008 D. Masson, pers. com., 20142010 18 Sep 2009 D. Masson, pers. com., 20142011 1 Oct 2010 D. Masson, pers. com., 20142012 12 Sep 2011 D. Masson, pers. com., 2014Table 2.2: Initialization dates (section 2.1.4) and profile sources for eachmodel base run (section 2.2.3).ForcingThe model is forced at the surface with hourly wind speed (U) and direction, cloudfraction (CF), air temperature and relative humidity observations from Environ-ment Canada ( Hourly 10 m U andwind direction observations were obtained at Environment Canada’s Sandheads CSweather station located near the model site at the edge of the Fraser River sedimentbank (figure 2.1). Hourly CF, relative humidity, and air temperature observationswere obtained at the Environment Canada weather station at Vancouver Interna-tional Airport (YVR, figure 2.1). U and wind direction have been transformed intovelocity components and rotated such that the u velocity is oriented cross-strait23toward 35◦NE and the v velocity is oriented along-strait toward 305◦NW. Windstress is calculated according to Large and Pond (1982). Solar irradiance is deter-mined by positioning the sun according to time of day and year and filtering thedownward radiation using a cloud model (AW-13).Total freshwater flux (Qt, volume/time) into the SoG is prescribed using dailyriver discharge measurements obtained by Environment Canada ( in the Fraser River at Hope and in the EnglishmanRiver at Parksville (figure 2.1). Englishman River discharge is used in this studyas a proxy for the contribution of small, rainfall-dominated rivers to the freshwaterbudget of the SoG (C-09). A total freshwater flux profile, Ft, is parameterized interms of Qt, model S, and the model mixing layer depth to represent the freshwaterflux into the model domain (AW-13).2.2 Results2.2.1 Distribution and seasonal variabilityModel S varies seasonally (∼5 to 30) in the top 10 m due primarily to changesin U and Qt (figure 2.5a and b). Surface S can be quite low (< 5) in May to Au-gust, particularly during large Qt. The surface freshwater lense during this periodreaches a maximum depth of about 5 m, below which lies a strong halocline. Thedepth of the mixing layer also varies seasonally (∼1 to 20 m) and is shallowestduring peak Qt (figure 2.5a). The mixing layer depth appears to be best describedby a combination of U3 and Qt. Strong winds and a weak halocline in Septemberthrough April coincide with a relatively deep mixing layer (> 5 m, figure 2.5a andb). Strong wind events in May through August are associated with mixing-layer-deepening and increased surface S (e. g. summer 2009), however, wind cannotmaintain a mixing layer as deep in summer as it can in winter due to increasedstratification.Background chl a biomass is low (< 1 mg m−3, figure 2.5c). An elevated near-surface signal (> 3 mg m−3, top 10 m) appears between February and April andpersists until sometime in October, except during large freshets when near-surfacechl a biomass is reduced to background concentrations. The high chl a signal24U3  [103 m3 s−3 ]Q t [104 m3 s−1 ]  01Depth (m)(a)20100U3 (observed) Qt (observed) MLD (model)(b) 2929 2521252125212521010203040S P [PSS−78]51525Depth (m)(c) 39399 3393939010203040Chl a [mg m3]051015(d)2221. DIC [103 µmol kg−1]11.52(e)7.77.788.388.388.388.32009 2010 2011 2012 2013010203040pH [total scale]7.788.3Figure 2.5: Four years (1 Jan 2009 through 31 Dec 2012) of (a) daily ob-served wind speed cubed (U3, light grey), daily total freshwater dis-charge (Qt, solid line), and modeled mixing layer depth (dark grey),(b) model S, (c) model chl a, (d) model DIC, and (e) model pH. U3 isobtained from hourly wind speed observations by Environment Canadaat Sandheads (figure 2.1; section 2.1.4) that is low-pass filtered ( f = 1day−1) and cubed. Qt is obtained from Environment Canada observa-tions of Fraser River discharge at Hope and Englishman River dischargeat Parksville (C-09; figure 2.1). Plots (b)-(e) are temporally averagedusing a 7-day running mean. The bold black contour (e) is the ΩA = 1horizon.25moves away from the surface in May and until September the chl a maximum isfound at ∼10 m depth. Masson and Pen˜a (2009) observe a similar summer 10 mchl a maximum in the southern SoG. The shift of the phytoplankton distributionaway from the surface in May directly follows the depletion of surface NO−3 , andthe depth of the chlorophyll maximum is consistent with the maximum depth ofNO−3 depletion (NO−3 < 1 µM) during this time.Model DIC is high (> 2,000 µmol kg−1) at the 40 m boundary relative tothe surface (figure 2.5d). Near-surface dilution or drawdown of DIC is evidentin the top 10 m during the times of reduced S and elevated chl a biomass. LikeS and chl a, DIC also follows a seasonal cycle in the top 10 m with the high-est (∼1,900 µmol kg−1) values occuring in December to April and the lowest(∼1,000 µmol kg−1) values occuring in June and July. Because of the high 40 mDIC boundary condition throughout the year (∼2,050 to 2,090 µmol kg−1, fig-ure 2.4), the vertical gradient of DIC is strongest during periods of low surfaceDIC (April through October).Model pH ranges from 7.6 near the 40 m boundary to 8.4 near the surface whensignificant chl a biomass is present (figure 2.5c and e). Near-surface pH variesseasonally with lower (∼7.8) values occurring in October to sometime in Februaryor March and a high (8.0 to 8.4) pH signal beginning in March and persistingthrough September (figure 2.5e). Water below about 10 m is generally low in pH(< 8.0), and water below about 20 m is generally undersaturated with respect toaragonite (ΩA < 1). The aragonite saturation horizon (ΩA = 1) is present at∼20 mduring most of the year, except from November through February and during large(> 6,000 m3 s−1) freshets when it shoals to the surface (figure 2.5a and e). Strongwind events appear to mix the water column across the chemocline, particularly inMarch and April before peak Qt onset. This mixing brings DIC-rich water to thesurface and lowers the surface pH (e. g. spring 2010, figure 2.5a, d and e).2.2.2 EvaluationI evaluate the model results against ten profiles of DIC, TA, pH and ΩA from theIOS carbon dataset. The profiles span all seasons and cover a three year period from2010 through 2012 with a December 2003 cast included. I additionally use profiles26of T , S and NO−3 from these casts to evaluate the model physics and biology in thecontext of the carbonate system.Station Latitude Longitude Dataset(s)Model site 49.125 N 123.558 W STRATOGEMDFO North 49.164 N 123.548 W IOS water properties, IOS carbonDFO South 49.030 N 123.437 W IOS water properties, IOS carbonDFO 1 49.144 N 123.613 W IOS water propertiesDFO 3 49.201 N 123.438 W IOS water propertiesDFO 4 49.227 N 123.344 W IOS water propertiesDFO 5 49.003 N 123.501 W IOS water propertiesDFO 7 49.055 N 123.372 W IOS water propertiesDFO 8 48.974 N 123.301 W IOS water propertiesDFO 9 49.016 N 123.221 W IOS water propertiesTable 2.3: Coordinates for the model and sampling locations (figure 2.1).Datasets (table 2.1) that contain casts at each respective location are listedin the far right column.CTD date Carbon sampling dates North South 1 3 4 5 7 8 97 Apr 2012 2 Apr 2012, 7 Apr 2012 F X X X X23 Jun 2011 21 Jul 2010, 23 Jun 2011, 16 Jul 2012 X F X X X X X11-12 Sep 2011 7 Aug 2010, 11 Sep 2011, 12 Sep 2011 F F X X X X X4 Dec 2003 4 Dec 2003, 31 Oct 2010 F X X X X X XTable 2.4: CTD casts used to estimate the spatial variability of T and S (fig-ure 2.6a and b) for the given sampling dates from the IOS carbon dataset(table 2.1). Locations of CTD casts (see table 2.3) are identified byX andF. Bold text indicates carbon profiles taken during the same cruise as theaccompanying CTD casts, and F indicates carbon sampling locations.The location and horizontal structure of the Fraser River plume is strongly in-fluenced by wind, rivers, and tides (Royer and Emery, 1982). Since the 1-D modelparametrizes freshwater input in terms of Qt only and the model site is fixed (Eule-rian rather than Lagrangian), it is not possible to capture movement of the plume,leading to disagreement at times between the model and observations at DFO Northand DFO South. To estimate horizontal spatial variability in T and S in the plumevicinity I use additional CTD casts from the IOS water properties dataset taken atlocations near the Fraser River mouth (figure 2.1a; table 2.3) during four of the IOS27carbon cruises (table 2.4). The difference between the minimum (or maximum) andmean observed values across all stations at each depth is interpreted as the lower(or upper) range of horizontal variability in the plume region during each of thefour cruises. The four T and S horizontal variability profiles are assigned to the tenT and S carbon cruise profiles by season, one per cruise (table 2.4). I estimate tem-poral variability as the minimum and maximum modeled values at a given depth ina 7-day window centered on the date of comparison with an observed profile.Depth (m)(a) 2 Apr 2012010203040 10 15 207 Apr 201210 15 2023 Jun 201110 15 2016 Jul 201210 15 20 T [oC]21 Jul 201010 15 207 Aug 201010 15 2011 Sep 201110 15 2012 Sep 201110 15 2031 Oct 201010 15 204 Dec 200310 15 20Depth (m)(b) 0102030405 15 25 5 15 25 5 15 25 5 15 25 SP [PSS−78]5 15 25 5 15 25 5 15 25 5 15 25 5 15 25 5 15 25Depth (m)(c)  010203040 5 15 25   5 15 25 5 15 25 5 15 25 NO3−1 [umol kg−1]5 15 25 5 15 25 5 15 25 5 15 25 5 15 25 5 15 25Model DFO North DFO SouthFigure 2.6: Profiles of modeled and observed (a) T , (b) S, and (c) NO−3 atDFO North and DFO South for ten carbon cruise dates (spring throughfall, IOS carbon dataset). The grey envelope around the solid modelcurve represents the model temporal variability at the given depth withina 7 day window. The errorbars represent the spatial variability in T andS estimated from CTD data (table 2.4). The thick box indicates wheresurface T and S are overestimated and underestimated, respectively, bythe model. The thin box indicates where 5 m NO−3 is underestimated bythe model.28Observed horizontal T and S variability range from near 0◦C and ∆S ∼ 0 on 7April 2012 and 4 Dec 2003 and below 20 m to ∼ 7◦C at 10 m on 11-12 Sep 2011(see error bars, figure 2.6a) and ∆S ∼ 15 at 5 m on 23 Jun 2011 (see error bars,figure 2.6b). Temporal T and S variability are highest near the surface on 11-12September 2011 and 23 Jun 2011, respectively (∼ 3◦C and ∆S ∼ 4). The greatestspatial and temporal S variability generally occurs during periods of strong FraserRiver discharge (freshet season, June-July). High spatial variability in S at 5 m inthe June 2011, July 2012, and July 2010 profiles (∆S ∼ 7 to 15, figure 2.6b) canbe interpreted as abrupt changes in the halocline depth with increasing/decreasingproximity to the river plume.aDICaTA = 65.8 = 45.4DIC or TA [103  µmol kg−1 ]2 Apr 201211.52 aDICaTA = 56.6 = 40.57 Apr 2012aDICaTA = 41.7 = 40.723 Jun 2011  aDICaTA = 50.6 = 48.216 Jul 2012aDICaTA = 77.3 = 56.421 Jul 2010aDICaTA = 71.8 = 47.17 Aug 201010 20 3011.52 aDICaTA = 71.8 = 49.111 Sep 201110 20 30aDICaTA = 78.8 = 42.6SP [PSS−78]12 Sep 201110 20 30aDICaTA = 46.1 = 46.131 Oct 201010 20 30aDICaTA = 58.3 = 51.04 Dec 200310 20 30DIC TA DICfit TAfitFigure 2.7: Linear regression of all DIC (solid) and TA (open) versus S ob-served at DFO North (square) and DFO South (circle) for the ten carbonsampling dates (IOS carbon dataset) in spring through fall. The regres-sion slopes (aDIC and aTA, solid and dashed lines, respectively) are usedto evaluate the model, specifically by estimating spatial variability ofDIC and TA (figures 2.8a and 2.8b) from the S variability in figure 2.6b.Overall, modeled T and S agree with observed T and S within the ranges ofvariability, particularly below 5 m and at the surface in September through April.Disagreement in T and S can be seen in 40% of the comparisons at the surfaceand is focused in July and August. The model underestimates surface S on 16 Jul2012, 21 Jul 2010, and 7 Aug 2010 by ∼3 to 7 (thick box, figure 2.6b) suggesting29a fresher surface or stronger plume influence in the model than at the samplinglocation. The model also overestimates T by ∼1 to 3◦C on those dates (thick box,figure 2.6a), which supports a stronger model river influence. Additionally, themodel overestimates surface T on 23 Jun 2011 by ∼ 6◦C and surface S on 31 Oct2010 by ∼11.Observed DIC and TA are to first order linear with respect to S for all ten pro-files (figure 2.7), indicating that dilution is the strongest driver of changes in bothquantities. Biological uptake and remineralization generally cause deviations froma linear relationship to S, particularly in DIC (Ianson et al., 2003). The linear re-gressions of DIC and TA with S (figure 2.7) are therefore approximations. Theslopes of these regressions, one per cast for each DIC and TA cast to S (figure 2.7),are used to scale the spatial S variability (figure 2.6b) to the two quantities (fig-ure 2.8a and b). Since these slopes are seasonally variable due to biological activityand changes to endmember concentrations, I use a unique regression slope for eachscaling (e. g. 2 Apr 2012 in figure 2.7 is used to estimate DIC and TA variabilityfrom S variability for 2 Apr 2012 in figure 2.6b). Model surface DIC and TA areunderestimated (∼ 65 to 185 µmol kg−1 and ∼65 to 100 µeq kg−1, respectively)on two of the three dates (16 Jul 2012 and 7 Aug 2010) where model surface S isunderestimated (thick box, figure 2.8a and b). These low estimates are consistentwith the closer plume proximity discussed earlier, however, the low model surfaceS estimate on 21 Jul 2010 (thick box center, figure 2.6b) is not accompanied by lowmodel surface DIC or TA estimates.On 21 Jul and 7 Aug 2010, DIC is underestimated by the model (∼ 90 to120 µmol kg−1) where TA is not (thin box, figure 2.8a and b). This DIC disagree-ment is located below the surface at ∼ 5 m, and is not consistent with model/datadisagreement in either T or S. Instead, these instances where the model underes-timates DIC but not TA occur during the times of greatest model NO−3 underpre-diction (thin box, figure 2.6c). These low (∼ 8.5 to 10 µM) model NO−3 estimatesare indicative of residual overproductivity in the tuned biology model and are cen-tered at 5 m due to the presence of a subsurface chl a maximum. The magnitudeof the DIC disagreement is ∼ 1.6 to 1.8 times larger than the product of the NO−3disagreement and the Redfield C:N ratio, suggesting that overproductivity couldaccount for the majority of the low model DIC estimates at 5 m on these dates.30Depth (m)(a) 2 Apr 20120102030401.0 2.07 Apr 20121.0 2.023 Jun 20111.0 2.016 Jul 20121.0 2.0 DIC [103 µmol kg−1]21 Jul 20101.0 2.07 Aug 20101.0 2.011 Sep 20111.0 2.012 Sep 20111.0 2.031 Oct 20101.0 2.04 Dec 20031.0 2.0Depth (m)(b) 0102030401.0 2.01.0 2.01.0 2.01.0 2.0 TA [103 µeq kg−1]1.0 2.01.0 2.01.0 2.01.0 2.01.0 2.01.0 2.0Depth (m)(c) 010203040 7.7 8 8.3 7.7 8 8.3 7.7 8 8.3 7.7 8 8.3 pH [total scale]7.7 8 8.3 7.7 8 8.3 7.7 8 8.3 7.7 8 8.3 7.7 8 8.3 7.7 8 8.3Depth (m)(d)  010203040 1.0 2.0   1.0 2.0 1.0 2.0 1.0 2.0 ΩA1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0Model DFO North DFO SouthFigure 2.8: Profiles of modeled and observed (a) DIC, (b) TA, (c) pH, and(d) ΩA at DFO North and DFO South for ten carbon cruises in springthrough fall (IOS carbon dataset). pH and ΩA are calculated from DICand TA (section 2.1.3). The grey envelope around the solid model curverepresents the model temporal variability at the given depth within a7 day window (minimum and maximum values). Errorbars representthe spatial variability in DIC and TA estimated from S (figure 2.6b, fig-ure 2.7), and in pH and ΩA estimated from DIC (a) and TA (b) andT and S (figure 2.6a and b). The thick box indicates where surface Tand S are overestimated and underestimated, respectively, by the model(figure 2.6a and b). The thin box indicates where 5 m NO−3 is underes-timated by the model (figure 2.6c).31High (∼ 9 µM) NO−3 disagreement also occurs on 16 Jul 2012 at 5 and 10 m (fig-ure 2.6c), however, the model overestimates NO−3 in this case and the disagreementdoes not appear to affect DIC relative to TA (figure 2.8a and b).Spatial variability of pH and ΩA calculated from DIC and TA observationsis determined (CO2SYS) using the spatial variability limits of DIC and TA (fig-ure 2.8a and b) and T and S (figure 2.6a and b). Low (or high) model estimatesof surface S (or T ) (thick box, figure 2.6a and b) that indicate a stronger plumeinfluence are associated with high model estimates of surface pH (∼ 0.07 to 0.27)and ΩA (∼ 0.0 to 0.9) on 16 Jul 2012, 21 Jul 2010, and 7 Aug 2010 (thick box, fig-ure 2.8c and d). The high model estimates of pH and ΩA during an elevated riverinfluence are unexpected since the Fraser River is generally low in pH relative tothe surface SoG with ΩA < 1. However, estuarine mixing where freshwater pH isless than seawater pH does not produce a linear mixing curve of pH with S, andpH maxima can occur at as low as S = 8 or as high as S = 17 (figure 3.8; Hofmannet al., 2009; Whitfield and Turner, 1986).Low model estimates of NO−3 at 5 m (thin box, figure 2.6c) that indicate over-productivity are associated with high model estimates of pH (∼ 0.0 to 0.2) and ΩA(∼ 0.6 to 1.1) on 21 Jul 2010 and 7 Aug 2010 (thin box, figure 2.8c and d). Con-trary to the behavior of pH and ΩA with increasing freshwater influence, the highmodel estimates of pH and ΩA attributed to overproductivity are intuitive resultsof reducing model DIC relative to TA. The model overestimates pH and ΩA mostseverely on 21 Jul 2010 (0.27 and 0.9 in the surface and 0.20 and 1.1 at 5 m, re-spectively). The collective disagreement of surface T and S and 5 m NO−3 is alsomost severe on that date.Direct comparisons of modeled and observed pH (figure 2.9) and ΩA (fig-ure 2.10) independent of depth further illustrate the relative importance of disagree-ment with respect to physics (T and S) and biology (NO−3 ). The model tendency tooverestimate pH and ΩA is clear at model pH > 8.0 (figure 2.9a and c) and modelΩA > 1.5 (figure 2.10a and c) and is accompanied by low model estimates of S andNO−3 . I evaluate model error due to physics (figures 2.9b and 2.10b) by recalculat-ing model pH and ΩA with subtracted model disagreement in T and S (figure 2.6aand b) and corrections to DIC and TA estimated from S regressions (figure 2.7).The physics corrections appear to improve the pH and ΩA agreement for some327.67.888.28.4Model pH(a) (b)∆SP [PSS−78]−8−6−4−2024687.4 7.6 7.8 8 8.2 8.4 pHCalculated pH(c)7.4 7.6 7.8 8 8.2 8.4 8.6Calculated pH(d)∆NO 3 [µmol kg−1 ]−8−6−4−202468Figure 2.9: Modeled versus observed pH for ten carbon cruises in springthrough fall (IOS carbon dataset). pH is calculated from DIC and TA(section 2.1.3). Color indicates S disagreement (model - observations)(a, b) and NO−3 disagreement (c, d). Model error due to physics (b) isevaluated by recalculating model pH with subtracted model disagree-ment in T and S (figure 2.6a and b) and corrections to DIC and TAestimated from S regressions (figure 2.7). Model error due to biology(d) is evaluated by recalculating model pH with subtracted correctionsto DIC and TA estimated from the Redfield ratio and electroneutral up-take assumptions (figure 2.3). Vertical errorbars represent the modeltemporal variability within a 7 day window (minimum and maximumvalues). Horizontal errorbars represent the spatial variability estimatedfrom DIC and TA (figure 2.8a and b) and T and S (figure 2.6a and b).Observations are from DFO North and DFO South. The 1:1 line is alsoshown (dash).330123Model ΩA(a) (b)∆SP [PSS−78]−8−6−4−2024680 1 2 30123Model ΩACalculated ΩA(c)0 1 2 3Calculated ΩA(d)∆NO 3 [µmol kg−1 ]−8−6−4−202468Figure 2.10: Modeled versus observed ΩA for ten carbon cruises in springthrough fall (IOS carbon dataset). ΩA is calculated from DIC and TA(section 2.1.3). Color indicates S disagreement (model - observations)(a, b) and NO−3 disagreement (c, d). Model error due to physics (b) isevaluated by recalculating model ΩA with subtracted model disagree-ment in T and S (figure 2.6a and b) and corrections to DIC and TAestimated from S regressions (figure 2.7). Model error due to biology(d) is evaluated by recalculating model ΩA with subtracted correctionsto DIC and TA estimated from the Redfield ratio and electroneutral up-take assumptions (figure 2.3). Vertical errorbars represent the modeltemporal variability within a 7 day window (minimum and maximumvalues). Horizontal errorbars represent the spatial variability estimatedfrom DIC and TA (figure 2.8a and b) and T and S (figure 2.6a and b).Observations are from DFO North and DFO South. The 1:1 line is alsoshown (dash).34points with low S estimates, however it is not clear whether the corrections im-prove the overall pH and ΩA agreement. I evaluate model error due to biology(figures 2.9d and 2.10d) by recalculating model pH and ΩA with subtracted correc-tions to DIC and TA estimated from the Redfield ratio and electroneutral uptakeassumptions (figure 2.3). The biology corrections clearly improve the pH and ΩAagreement at several points with low NO−3 estimates, and they appear to improvethe overall agreement as well.Despite the model-data disagreement in pH and ΩA discussed, the model gener-ally reproduces observations of the two quantities consistently across the ten sam-pling dates from the IOS carbon dataset. Where model-data disagreement arises,the tendancy of the model to overestimate the freshwater influence and primaryproductivity during the summer leads to low predictions of near-surface DIC andhigh predictions of near-surface pH and ΩA. Near-surface pH and ΩA are overes-timated by as much as 0.27 and 1.1, respectively. The two quantities are rarely,and never severely, underestimated by the model. Because low pH and aragonitesaturation state are unfavorable for the SoG ecosystem, I consider the high modelestimates during disagreement to represent a best case scenario.2.2.3 Sensitivity studiesNear-surface (top 10 m) pH and ΩA appear to follow the seasonal cycle of chl a(figure 2.5c and e). The appearance of elevated chl a biomass (> 3 mg m−3) be-tween February and April represents the onset of the spring diatom bloom, thetiming of which is determined primarily by U and CF averaged from December toMarch and is weakly sensitive to Qt averaged over the same period (AW-13). Uand Qt can also impact pH and ΩA directly (e. g. March 2010 wind events appear tomix pH across depth, June/July 2011 and 2012 freshets occur with surface ΩA < 1,figures 2.5a and e).I designed a series of sensitivity runs in order to quantify the importance ofU , Qt, and CF on the seasonal distribution of pH and ΩA. The model base runswere initialized between September and October from 2001 to 2012 on the date ofthe available sampling cast, unique for each year (table 2.2), and run through theend of the following year. For each base run, 3 sets of 12 experimental runs were35performed (table 2.5). The first set, experiment 1, cycles through 12 separate years(2001-2012) of U data while keeping Qt and CF the same as the base run year.Experiments 2 and 3 use the same years of Qt and CF data, respectively. The 3experiments (12 base runs × 3 sets of 12 experimental runs each) total 432 runs. Ichose these 432 runs to create a wide range of realistic scenarios for the southernSoG.Experiment U Qt CF Runs1 2001-2012 Base year Base year 1442 Base year 2001-2012 Base year 1443 Base year Base year 2001-2012 144Table 2.5: Summary of model sensitivity experiments. Columns indicate theyear(s) of data used for each forcing quantity. Each experiment is per-formed for all base years (table 2.2) totaling 432 runs.The standard metrics I define for evaluation against the varied forcing quanti-ties are pH averaged from the surface to 3 m (pH3m), the onset of spring surfacearagonite supersaturation (ΩAsur > 1), the onset of winter surface aragonite under-saturation (ΩAsur < 1), and the duration of summer surface ΩAsur < 1. These metricsare also evaluated against the depth-integrated (0 to 40 m) mid-day PP (PPint) al-though PP is generally concentrated in the top 10 m (figure 2.5). I selected a 30-daybackaveraging window based on the biological residence time (appendix D) in or-der to standardize pH3m, PPint, and forcing quantities.The correlation between a model metric and a forcing quantity using a specificbackaveraging window across all relevant sensitivity runs provides the basis for theevaluation of these experiments. I quantify the correlation using linear regression.Figure 2.11 shows four of these regressions during times when the correlation isrelatively strong (r2 > 0.5 except figure 2.11c). pH3m is anticorrelated to U duringmid-March to mid-April (figure 2.11a) and anticorrelated to Qt during mid-Juneto mid-July (figure 2.11b). Spring ΩAsur > 1 onset is positively correlated to CFaveraged over February (figure 2.11c) and summer ΩAsur < 1 duration is positivelycorrelated to Qt averaged over July (figure 2.11d). The four plots suggest thatincreased U and CF during February, March, and early April and increased Qtduring June and July all act to reduce near-surface pH and ΩA.366 7 8 988. 15U (m s−1)pH3mR2 = 0.57p < 0.01(a)6 8 1088. 10Qt (103 m3 s−1)pH3mR2 = 0.93p < 0.01(b)0.5 0.6 0.7 0.8 0.9455055606570Mar 6CFSpring ΩA sur > 1 onset (yearday) R2 = 0.49p < 0.01(c)4 6 8 10010203040Jul 25Qt (103 m3 s−1)Summer Ω A sur < 1 duration (days) R2 = 0.81p < 0.01(d)Figure 2.11: Scatterplots of (a) experiment 1: model pH3m versus observedU backaveraged 30 days from 15 Apr, (b) experiment 2: model pH3mversus observed Qt backaveraged 30 days from 10 Jul, (c) experiment3: model spring ΩAsur > 1 onset versus observed CF backaveraged 30days from 6 Mar (CF only), and (d) experiment 2: model summerΩAsur < 1 duration versus observed Qt backaveraged 30 days from 25Jul (Qt only). Experimental runs sharing the same U (a), Qt (b and d),and CF (c) timeseries are averaged and errorbars represent the stan-dard deviation (n = 12) of each average (e. g. all base runs using the2005 U record are summarized as a point with errorbars). Linear re-gression is performed using the unaveraged experimental runs, and pis determined from a Student’s t-test.3 m averaged pHIt is difficult to characterize the relationship between any two quantities through-out the year using only single specific windows (e. g. figure 2.11). Additionally,37these snapshots do not address possible lag time between a forcing event and theresponse of a model metric. Thus, I also consider a timeseries of R2 values, eachpoint representing a regression between pH3m and either U , Qt, or CF for a par-ticular backaveraging window determined by the yearday (figure 2.12b). To helpme determine times when biological interactions are important, I further consider atimeseries of the correlation between pH3m and PPint (figure 2.12b). I also considerthe variability of the annual cycle of 30-day backaveraged pH3m using the standarddeviation and minimum and maximum values across the 432 sensitivity runs (fig-ure 2.12a). A variable lag between the averaging windows for the independent(i. e. forcing, PPint) and dependent (i. e. pH3m) quantities is chosen to maximizethe value of R2 (appendix E). The lag designates the number of days that the inde-pendent quantity backaverage precedes the dependent quantity backaverage. Lagperiods are limited to be less than or equal to 30-days. Those that exceeded 30 daysdid not significantly increase the value of R2.The strong anticorrelation between pH3m and U (R2 > 0.5) during mid-Marchto mid-April and between pH3m and Qt (R2 > 0.9) during mid-June to mid-July(figure 2.11a and b) correspond to peaks in the timeseries of pH3m correlationagainst U and Qt (figure 2.12b), although the values of R2 in the timeseries arehigher due to the lag. The width of each peak represents the period during whichthe relationship is strong (∼ 2 to 3 months). Lag periods between the forcingquantities and pH3m are generally short (∼ 5 to 10 days for U , Qt, and CF) dur-ing periods of strong correlation (figure E.1a). Anticorrelation between pH3m andU is relatively strong (R2 > 0.5) throughout the year except between June andAugust, and reaches a maximum of 0.8 in March and April. Anticorrelation be-tween pH3m and Qt is isolated to June through early August, and a period of weakpositive correlation between the two quantities occurs in mid-April to mid-May(figure 2.12b). Anticorrelation between pH3m and CF is weak (R2 < 0.2) except inSeptember through October (R2 > 0.6, figure 2.12b). Additional sensitivity exper-iments in which air temperature was varied resulted in even lower and generallynon-significant strengths of correlation (not shown).The seasonal cycle of near-surface pH discussed in section 2.2.1 appears acrossall sensitivity runs, however there are two periods, February through April (spring)and mid-May through mid-August (summer), where pH3m variability is high (total38(a) pH3m variabilitypH3m  7.888.28.4Mean Std Error Min/MaxR2(b) pH3m correlation  J F M A M J J A S O N D J00. Qt CF PPintFigure 2.12: Timeseries of (a) mean (dash), standard error (light gray), andminimum and maximum (dark gray) 30-day backaveraged pH3m, and(b) 30-day backaveraged pH3m correlation against 30-day backaver-aged U (bold), Qt (dash dot), CF (solid) and PPint (patch). Regres-sions against U , Qt, and CF (e. g. figure 2.11) are performed usingthe set of runs from experiments 1, 2 and 3, respectively (table 2.5),and regressions against PPint are performed using the total set of 432experimental runs. pH3m is backaveraged from the indicated date onthe x-axis and the U , Qt, CF, and PPint backaverages precede the in-dicated date by a lag which maximizes R2 (figure E.1b). Red linesand patches indicate anticorrelation while black lines and gray patchesindicate positive correlation.range > 0.3) relative to the rest of the year (figure 2.12a). The spring period of highpH3m variability occurs with the maximum anticorrelation between pH3m and U(R2 > 0.8). Likewise, the summer period of high pH3m variability occurs with themaximum anticorrelation (R2 > 0.9) between pH3m and Qt (figure 2.12b). Thereis also a third period of increased pH3m variability (∼ 0.2) in October throughNovember (autumn) associated with a rise in anticorrelation between pH3m and U(figure 2.12a and b).There are two time periods during the year when pH3m is strongly correlated(R2 > 0.8) to PPint (figure 2.12b). The earliest, February through April, nearlyoverlaps the spring period of peak anticorrelation between pH3m and U . Consider-39ing the role wind has been shown to play in the timing of the spring phytoplanktonbloom (C-09), this overlap suggests that the effect of U on pH3m here is at least par-tially biologically mediated (e. g. wind attenuation of spring productivity). Sim-ilarly, the strong correlation between pH3m and PPint from June through Augustduring summer anticorrelation between pH3m and Qt (figure 2.12b) suggests thatpeak outflow negatively affects PPint which in turn lowers pH3m. A period of ele-vated correlation (R2 ∼ 0.6) between pH3m and PPint from mid-September throughNovember (figure 2.12b) is associated with the autumn period of increased pH3mvariability (figure 2.12a).There are considerable lag times that follow the onset of the peaks in correla-tion of pH3m to PPint (up to 30 days, figure E.1a). Most of the backaveraged pH3mvalues during these periods are actually best correlated with backaveraged PPintnear the beginning of the period of strong correlation. This persisting correlationindicates that PP is not a continuous driver of pH, but rather that the effects of bi-ological events, such as the spring bloom, on pH tend to last for periods as long as2 months.R2(a) ΩA sur > 1 onsetJ F M A00. (b) ΩA sur < 1 durationM J J A S O   (c) ΩA sur < 1 onsetN DU Qt CF PPintFigure 2.13: Correlation timeseries of (a) spring ΩAsur > 1 onset, (b) summerΩAsur < 1 duration, and (c) winter ΩAsur < 1 onset versus U (bold),Qt (dash-dot), CF (solid), and PPint (patch). Regressions against U ,Qt, and CF (e. g. figure 2.11) are performed using the set of runs fromexperiments 1, 2 and 3, respectively (table 2.5), and regressions againstPPint are performed using the total set of 432 experimental runs. Redlines and patches indicate anticorrelation while black lines and graypatches indicate positive correlation.40Surface ΩAThe correlation between spring ΩAsur > 1 onset and CF (figure 2.11c) can be seen atthe R2 peak in February and March (figure 2.13a). Similarly, the strong (R2 > 0.8)correlation between summer ΩAsur < 1 duration and Qt (figure 2.11d) can be seenas a peak in the correlation timeseries in June through July (figure 2.13b). Thecorrelation between spring ΩAsur > 1 onset and CF is overlapped by a nearly iden-tical correlation between spring ΩAsur > 1 onset and PPint (figure 2.13a). Likewise,the correlation between summer ΩAsur < 1 duration and Qt overlaps with a similarcorrelation between summer ΩAsur < 1 duration and PPint (figure 2.13b). These cor-relations to PPint suggest that the indicated relationships between spring ΩAsur > 1onset and PPint (figure 2.13a) and summer ΩAsur < 1 duration and Qt (figure 2.13b)are biologically-mediated. These biological connections are supported by the sim-ilarities between model ΩA and model chl a described earlier (figure 2.5). Thewinter ΩAsur < 1 onset is only weakly correlated to U , Qt, CF, and PPint duringOctober and November (figure 2.13c).2.3 Discussion2.3.1 MechanismsSpring pH versus windIncreased U is strongly associated with reduced pH3m during the spring (Februarythrough April) period of high pH3m variability (figure 2.12a and b). I propose twomechanisms for this relationship. The first mechanism is strictly physical, althoughit relies on the prior biological drawdown of near-surface DIC (and increase inpH). Increased U during spring mixes DIC-rich water from below the chemoclineinto the near-surface, thus lowering pH3m. This mechanism is supported by twokey features of the pH3m variability timeseries (figure 2.12a). The low variabilityof pH3m prior to the spring season represents winter conditions where the watercolumn is well-mixed and the vertical DIC gradient is weak. Wind mixing of thewater column will have little effect on pH3m during this time. The low variability ofpH3m following the spring season arises when the effect of freshwater stratification41on pH3m becomes stronger than that of wind mixing (figure 2.12a). Since U andQt have competing effects on pH3m in May, pH3m is stabilized.The second mechanism is biological, in that it directly involves changes in PP(which are driven by physics). Increased U during spring will either deepen themixing layer depriving phytoplankton of light, or transport phytoplankton belowthe mixing layer all together. This mechanism is supported by the strong corre-lation between pH3m and PPint during spring (figure 2.12b). This mechanism isalso consistent with the low pH3m variability in December through February (win-ter) and in May (figure 2.12a) since winter is a time of low solar irradiance andfreshwater stratification in May is generally strong enough to resist wind mixing.In order to explore which mechanism is the primary driver of the relationshipbetween U and pH3m in spring, I examine DIC fluxes from vertical mixing and netPP. The DIC mixing flux, ΦmixDIC, is defined according to Fick’s first law evaluatedat depth dΦmixDIC =− K(z)∂DIC∂ z∣∣∣∣−d(2.3)where K(z) is the vertical diffusivity determined by the KPP model and d is thedepth that maximizes ΦmixDIC. I further require d to be greater than 5 m to ensuremixing of DIC from below the near-surface. The DIC biological flux, ΦPPDIC, isdefined as the depth-integrated difference of DIC uptake and remineralizationΦPPDIC =∫ 0−D[uptake− remin]dz (2.4)where D is the depth of the model domain. Biological processes, including rem-ineralization, in the model are near-surface intensified so integration over the modeldepth resolves the net fluxes without distorting their magnitude.The magnitude and variability of ΦmixDIC are about a factor of 2 larger than thoseof ΦPPDIC (figure 2.14a). A preliminary survey would suggest that vertical mix-ing has a stronger impact on near-surface DIC than biological activity. However,the correlations between ΦmixDIC and ΦPPDIC provide a more complete scenario. Bothfluxes, when backaveraged over a 30 day window and regressed against U , can bepresented in terms of their correlation coefficients as timeseries similar to those infigure 2.12 (figures 2.14c). In March, ΦmixDIC is anticorrelated to U suggesting, sur-42prisingly, that increased winds prevent net DIC transport to the near-surface (andreductions in near-surface pH) rather than enhancing it. However, the magnitudeof ΦPPDIC is also anticorrelated to U in March. This relationship suggests that thereason U seems to be preventing DIC transport rather than enhancing it is becausethe biological DIC uptake fluxes necessary to establish vertical DIC gradients areweak. The net DIC transport over these weak gradients is small as a result, thusexplaining the anticorrelation between ΦmixDIC and U . In April, the correlations arereversed.The turning point where R changes sign signifies a fundamental shift in thefunction of U on PP (ΦPPDIC). In March, increased U tends to suppress PP by deep-ening the mixing layer and delaying the spring bloom. In April, however, bloomshave already occurred in several of the model scenarios, and the role of U haschanged such that wind mixing provides nutrients and enhances PP. The function ofU on ΦmixDIC is determined by ΦPPDIC. In March, net vertical DIC mixing fluxes (ΦmixDIC)are small unless a strong vertical DIC gradient has been established by near-surfacebiological uptake. In April, a strong vertical DIC gradient has been established andΦmixDIC is more directly coupled to wind mixing. This interpretation indicates that bi-ological fluxes are primarily responsible for the relationship between pH3m and Uprior to April during the bloom, and mixing fluxes are primarily responsible for therelationship between pH3m and U from April onward following the bloom. In thecontext of equation 2.3, variations in the vertical gradient of DIC (∂DIC/∂ z) formthe primary driver of pH3m prior to April, and variations in the vertical diffusivity,K(z), form the primary driver of pH3m from April onward.Summer pH versus freshwaterIncreased Qt is strongly associated with reduced pH3m during the summer (mid-May through mid-August) period of high pH3m variability (figure 2.12a and b). Ipropose four potential mechanisms for this relationship: (1) freshwater dilution ofTA relative to DIC, (2) N limitation due to reduced entrainment velocity, (3) lightlimitation due to increased turbidity, and (4) loss of phytoplankton due to increasedestuarine advection. All four mechanisms are consistent with the negative correla-tion between Qt and pH3m during the freshet season (figure 2.12b). However, the43µmol m−2  s−1Spring  (a) ΦDIC variabilityΦDICmixΦDICPP−505Mean Std Error Min/MaxR2(c) ΦDIC U correlation  M A M00.ΦDICmixSummer(b) ΦDIC variabilityΦDICfreshΦDICPP(d) ΦDIC Qt correlation  J J A SΦDICfresh ΦDICPPFigure 2.14: (a,b) Timeseries of mean (dash), standard error (light patch), andminimum and maximum (dark patch) 30-day backaveraged (a) ΦmixDIC(gray) and ΦPPDIC (green), and (b) ΦfreshDIC (blue) and ΦPPDIC (green). (c,d)Timeseries of 30-day backaveraged ΦmixDIC (bold), ΦfreshDIC (dash-dot) andΦPPDIC (solid) correlation againt (c) 30-day backaveraged U , and (d)30-day backaveraged Qt. Regressions against U and Qt are performedusing the set of runs from experiments 1 and 2, respectively (table 2.5).ΦmixDIC, ΦfreshDIC , and ΦPPDIC are backaveraged from the given date and theU and Qt backaverages precede the indicated date by a lag which max-imizes R2 (figure E.2). Red lines and patches indicate anticorrelationwhile black lines and gray patches indicate positive correlation.strong summer correlation between pH3m and PPint would suggest mechanisms (2)through (4) since freshwater dilution of TA relative to DIC could occur indepen-dently of changes in PP.In order to determine which mechanism is driving the relationship betweenpH3m and Qt in summer (figure 2.12b), I examine an additional DIC flux attributedto freshwater dilution. Scalar fluxes from freshwater into the model are more accu-rately defined in terms of the parameterized total freshwater flux, Ft, rather than Qtsince Ft describes the freshwater received at the model location. Ft (volume/time)44is parametrized in terms of Qt, mixing layer S, and mixing layer depth (AW-13),and is approximately linear with Qt (ratio ∼ 2.0× 10−8) in the model surface atQt < 6,000 m3 s−1 and approximately constant at Qt > 6,000 m3 s−1. The DICfreshwater dilution flux, ΦfreshDIC , is defined as the product of Ft and the DIC circula-tion strength, DICb− DICrΦfreshDIC = Ft(DICb−DICr)(2.5)where DICb is the DIC 40 m boundary condition (constant, section 2.1.4) and DICris the freshwater DIC value (variable, section 2.1.3).Despite the large DIC loss due to freshwater dilution (figure 2.14b), the effectof ΦfreshDIC on pH3m is determined by the relative loss of DIC with respect to TA.The magnitude of ΦfreshDIC is therefore a weaker indicator of the relationship betweenpH3m and Qt than the correlation between ΦfreshDIC and Qt (figure 2.14d). I considerΦfreshDIC to be a potentially strong driver of summer pH3m variability if (a) ΦfreshDICdemonstrates high variability during summer (indicating changes driven primarilyby Qt), and (b) the correlation between ΦfreshDIC and Qt is strong. The variability ofΦfreshDIC between mid-May and mid-August is overall on the same order as that ofΦPPDIC (figure 2.14b), but the correlation between ΦfreshDIC and Qt reaches a lowpointduring the freshet season (figure 2.14d). The weak correlation between ΦfreshDIC andQt suggests that freshwater dilution of DIC and TA is not a strong driver of pH3mvariability. Conversely, the strong anticorrelation between ΦPPDIC and Qt betweenmid-May and mid-August (figure 2.14d) suggests that PP is the primary driver ofsummer pH3m variability.Finally, in order to identify the dominant biological mechanism driving therelationship between pH3m and Qt in summer, I compare the net entrainment fluxof DIN, ΦEN, (vertical entrainment - horizontal advection) over the model domainto depth-integrated gross N-uptake, ΦGPPN , during a small freshet year (2010, ∼6,000 m3 s−1 peak outflow) and a large freshet year (2011, ∼ 11,000 m3 s−1 peakoutflow). There is a dramatic reduction in ΦEN from 2010 to 2011 accompanied by asimilar reduction in ΦGPPN (figure 2.15). The time-integrated fluxes (June-July) forΦGPPN in 2010 are near those for ΦEN (∼ 86%). In other words, primary producersuse most of the entrained N. However, in 2011 summer time-integrated ΦGPPN is45−1−0.500.5µmol N ms−1 s−1(a) 2010May Jun Jul Aug−1−0.500.5µmol N ms−1 s−1(b) 2011  Net Entrainment ΦNUptake ΦNFigure 2.15: Timeseries of modeled ΦEN (dash) and ΦGPPN (solid) from Maythrough mid-August for (a) 2010 and (b) 2011.only∼ 37% of that for ΦEN. This comparison suggests that although phytoplanktonare generally N-limited during the freshet season, and that ΦEN is reduced at largeQt, phytoplankton are limited by another mechanism during large freshets. Thismechanism is most likely light limitation or estuarine advection of phytoplanktonor both; processes that are not easily separated in the model.Surface ΩAOf the three forcing quantities used in the sensitivity experiments, spring surfaceΩA > 1 onset is most strongly correlated to CF (figure 2.13a). The nearly identicalstrength of correlation (R2) demonstrated between CF and PPint (figure 2.13a) sug-gests that the effect of CF on spring surface ΩA > 1 onset is biologically-mediated(e. g. light-availability to phytoplankton). Rapid phytoplankton growth due tochanges in light-limitation in early spring is analagous to the spring bloom, andprevious studies have demonstrated spring bloom timing and CF to be correlated(AW-13; C-09). Next I consider how well-correlated spring bloom timing and46spring surface ΩA > 1 onset are to each other.The spring bloom timing is subject to definition. In the SoG, previous estimatesfrom modeling studies and satellite observations have defined the bloom in termsof both increases in phytoplankton concentration and decreases in surface nutrientabundance (AW-13). The spring bloom signifies a transition between nutrient-replete and nutrient-limited surface conditions. I choose a definition based on thistransition: the time at which the surface NO−3 is reduced below 15 µM for at least1 day. The value of 15 µM is low enough to ensure that the reduction is biological,and high enough to represent the onset of spring PP rather than the peak. I furtherdefine the bloom magnitude as the peak chl a biomass within a 15 day window ofthe bloom onset date.30 40 50 60 70 80 9030405060708090Spring bloom onset [yearday]Spring Ω A sur > 1 onset [yearday]  1:1Chl a [mg m−3 ]468101214161820Figure 2.16: Spring ΩAsur > 1 onset versus spring bloom onset for all 432experimental runs, colormapped by chl a biomass. The 1:1 line is alsopictured (solid black line).There is a weak correlation between spring ΩAsur > 1 onset and spring bloomtiming (figure 2.16). The two quantities approach a 1:1 relationship, however whenlow chl a values are excluded (< 10 mg chl a m−3), the spring ΩAsur > 1 onset pre-47cedes the spring bloom across all 432 sensitivity runs. This relationship indicatesthat spring PP leads to surface aragonite supersaturation more rapidly than NO−3depletion. Furthermore, the constant threshold of NO−3 depletion does not corre-late strongly with the spring ΩAsur > 1 onset, and the relationship between ΩA andNO−3 consumption is likely more complicated.The appearance of surface aragonite undersaturation during large freshets hasbeen addressed multiple times in this chapter. The strong correlation between sum-mer ΩAsur < 1 duration and Qt was demonstrated in figure 2.13b. As increased Qtreduces summer pH3m, it follows that increased Qt will reduce ΩAsur since pH influ-ences the carbonate equilibrium. The strong correlation between summer ΩAsur < 1duration and PPint that overlaps the correlation to Qt (figure 2.13b) supports theconnection to pH. However, it is equally plausible that freshwater dilution of DIC(and thus CO2−3 ) is reducing ΩAsur independently of pH. It is further possible thatfreshwater dilution of Ca2+ is also reducing ΩAsur .J F M A M J J A S O N D0.511.522.533.52010 / 2011  Ca2+CO32−Figure 2.17: 2010/2011 ratios of modeled Ca2+ (bold) and CO2−3 (dash).According to the definition of ΩA (equation 1.1), CO2−3 and Ca2+ diluted bythe same factor should contribute equally to reductions in ΩA. Model Ca2+ scaleslinearly with S on the modeled range when changes in T are small (Riley andTongudai, 1967), however, pH (and therefore CO2−3 ) does not (see chapter 3, fig-ure 3.8). To examine the relative dilution of Ca2+ and CO2−3 , I compute the ratiosof both quantities between a small freshet year (2010, ∼ 6,000 m3 s−1) and a largefreshet year (2011, ∼ 11,000 m3 s−1). Despite the non-linearities associated withCO2−3 dilution, both Ca2+ and CO2−3 appear to be diluted approximately equally48(∼ 3 fold) between the summers of 2010 and 2011 (figure 2.17). The change inCO2−3 may not be attributed solely to dilution, however, because of the reductionin PP.310320330340R2 = 0.28p < 0.01(a) All UΩ A sur < 1 Onset (yearday)R2 = 0.52p < 0.01(b) Qt < 2,500 m3 s−16 7 8310320330340R2 = 0.55p < 0.01(c) Qt < 2,500 m3 s−1 & CF < 0.75Ω A sur < 1 Onset (yearday)U (m s−1)6 7 8R2 = 0.55p < 0.01(d) Qt < 2,500 m3 s−1 & CF > 0.75U (m s−1)Figure 2.18: Winter ΩAsur < 1 onset results from sensitivity experiment 1 ver-sus observed U backaveraged 30 days from 13 Nov for (a) all runs,(b) backaveraged Qt < 2,500 m3 s−1, (c) Qt < 2,500 m3 s−1 and CF< 0.75 and (d) Qt < 2,500 m3 s−1 and CF > 0.75. Experimental runssharing the same U (e. g. the 2005 U record) are averaged and error-bars represent the standard deviation of each average. Linear regres-sion is performed using the unaveraged experimental runs, and p isdetermined from a Student’s t-test.Winter surface aragonite undersaturation is a persistent feature of the model(figure 2.5) and is supported by late October 2010 and December 2003 cruise data(figure 2.8). The presence of a well-mixed water column and low PP make thewinter near-surface properties similar to those of the intermediate layer (i. e. 40 mboundary conditions). While spring ΩAsur > 1 onset is well-correlated to CF, corre-49lations between the onset of winter ΩAsur < 1 and the three forcing quantities wereweak (< 0.2) across the sensitivity experiments (figure 2.13c).The correlation between winter ΩA < 1 onset and U improves dramatically(R2 2-fold) when runs with Qt > 2,500 m3 s−1 are excluded (figure 2.18a and b).The improvement in correlation indicates that increased fall U accelerates winterΩA < 1 onset, and that this relationship is complicated by strong fall freshwaterfluxes. Further separation by CF demonstrates the effect of CF on ΩA < 1 onset. Atlow U and Qt , experiment 1 runs with CF< 0.75 produce average onsets 5 - 10 dayslater than runs with CF > 0.75 (figure 2.18c and d). As the correlation betweenspring CF and ΩA > 1 onset is linked to spring productivity, the fall correlationsbetween U , CF, and ΩA < 1 onset are most likely linked to fall productivity. Herelight limitation and removal of phytoplankton by mixing attenuate fall blooms thusreducing surface ΩA, and just as in the spring relationship between pH3m and U ,enhanced DIC transport to the surface may be a possible driver of early winterΩA < 1 onset as well.Seasonal summaryThe onset of spring ΩAsur > 1 occurs sometime between mid-February and mid-March (day 45 to day 80, figure 2.16). Thinning of cloud cover between mid-January and April determines the onset timing (figure 2.13a), which is driven byPP but precedes the spring bloom timing (figure 2.16). This onset is the beginningof a ΩAsur > 1 signal that persists until anywhere from the end of October to thebeginning of December, unless interrupted by large Qt (e. g. summers 2011 and2012, figure 2.5).Anticorrelation between pH3m and U begins to exceed the background winteranticorrelation (R2 > 0.6) in mid-February, around the time of spring ΩAsur > 1 on-set (thick red curve, figure 2.12b), although the variability of pH3m starts to increaseas early as the beginning of February (figure 2.12a). Increases in U drive decreasesin pH3m from mid-February through April, however, these decreases occur by pre-dominantly different mechanisms depending on the time of year. From Februarythrough March, wind relaxation allows PP to decrease the near-surface DIC con-centration, increasing pH3m. Wind inhibition of PP is more strongly coupled than50wind mixing to increases in near-surface DIC during this time. The role of windchanges around the beginning of April. Since the strong vertical gradient of DICis already established in most model scenarios by the end of March, wind eventsin April to mid-May increase the near-surface DIC primarily by mixing rather thanby supressing PP.The variability of pH3m throughout the sensitivity experiments decreases inlate May (figure 2.12a), suggesting that none of U , Qt, or CF are strong drivers ofpH3m during that time. The variability of pH3m increases again in June and remainshigh through mid-August (figure 2.12a). The anticorrelation between pH3m andQt increases during this time and reaches a maximum from mid-June to mid-July(figure 2.12b). This summer relationship between pH3m and Qt is biologically-mediated, i. e. increases in Qt inhibit PP which in turn reduce pH3m. Since June-July 2011 time-integrated ΦGPPN is only ∼ 37% of time-integrated ΦEN compared to86% in summer 2010 (figure 2.15), N-limitation due to reduced upward entrain-ment is not likely. It is more likely that phytoplankton are either light-limited dueto river turbidity, advected out of the water column near the surface, or both.During some summers when Qt is large, a period of ΩAsur < 1 occurs. Theduration of this period is strongly correlated with Qt (figure 2.13b). The concurrentreduction of pH3m due to suppression of DIC uptake during large Qt suggests thatΩAsur < 1 may also be a result of PP suppression by Qt. However, freshwaterdilution of both Ca2+ and CO2−3 are strong enough during large Qt to reduce ΩAsurby as much as 3-fold (equation 1.1, figure 2.17), which is enough to create surfaceundersaturation by dilution alone.The variability of pH3m is weak from September through the remainder of theyear (figure 2.12a), indicating that U , Qt, and CF are not strong drivers of pH3m dur-ing this period. Variability of pH3m does increase slightly in late September, whichis consistent with a rise in anticorrelation between pH3m and U (figure 2.12b).Winds and CF in mid-October determine the winter ΩAsur < 1 onset, bringing thesystem once again back to a low pH, undersaturated, and well-mixed state.512.3.2 Spatial distributionFraser RiverThe summer relationships between Qt and pH3m (figure 2.12) and between Qtand ΩAsur (figure 2.13) demonstrate the importance of Fraser River discharge onthe carbonate chemistry of the southern SoG. However, the spatial distribution offreshwater during the freshet season is complicated by the variable position of theFraser River plume. Harrison et al. (1991) characterized freshwater in the southernSoG as two plumes: a riverine plume and an estuarine plume, the latter repre-senting the residual structure of the former after entraining deeper, higher-S wa-ter. Naturally, the estuarine plume covers a larger area than the riverine plume,generally representing the surface SoG south of Texada Island (figure 2.1). Theriverine plume covers a more localized area bounded approximately by the sam-pling stations in figure 2.1 (lower panel). Riverine plume boundaries are typicallypronounced, characterized by fronts and positioned between regions of contrastingvertical structure. Riverine plume location is driven by wind, tides, and the Cori-olis force, and will vary at a given discharge with these processes. The model issituated to represent either estuarine plume or riverine plume conditions dependingon Qt.Summer anticorrelation between Qt and ΦPPDIC is consistent with observed prop-erties in the Fraser River plume region. Halverson and Pawlowicz (2013) estimatedthe mean annual chl a maximum depth from ferry observations to be 3 m shallowerin the riverine plume than in the estuarine plume. They further estimated the meanannual riverine plume depth-integrated phytoplankton biomass per area to be 62-83% of that in the estuarine plume. Harrison et al. (1991) also found lower PP andless DIN drawdown in the riverine plume relative to the estuarine plume during aJuly 1987 cruise.The degree to which the model represents the riverine plume versus the estu-arine plume is determined by a light attenuation parameterization with a thresh-old Qt of ∼7,240 m−3 s−1 (C-09). Model outputs for 2011 and 2012, both large(> 11,000 m−3 s−1) freshet years, demonstrate shallow chl a maximum depths withchl a nearly absent below 5 m. These results are consistent with riverine plume con-52ditions. By contrast, model runs for 2009 and 2010, both smaller (< 7,000 m−3 s−1)freshet years, had pronounced chl a maximum depths just shallow of 10 m duringthe summer freshet. Despite a critical Qt threshold that causes the model to switchbetween estuarine plume conditions and riverine plume conditions, the relation-ships between Qt and pH3m and ΩAsur are gradual and support freshwater impactson carbonate chemistry at both ends of the discharge spectrum. Whether the modelis in the riverine plume or the estuarine plume, the summer sensitivity results holdfor both areas. Nevertheless, the freshwater impacts on summer carbonate chem-istry demonstrated by the sensitivity studies are likely to become less importantwith decreasing proximity to the riverine plume.Other basins and passagesThe Fraser River discharges into the southern reach of the southern SoG basin. Thesouthern basin extends to ∼100 km northwest of the rivermouth to Texada Island,where it joins the northern basin (figure 2.1a). Summer stratification due to theFraser River is important throughout the southern and northern basins for main-taining high primary productivity (Masson and Pen˜a, 2009). However, during thesummer the Fraser River plume appears to inhibit PP in the model. This inhibitionis likely due at least in part to light limitation from high turbidity near the riverplume, and may weaken with decreasing proximity to the river mouth. Weakeningof the inhibition of PP to the north of the plume is evident in recent observations.Mean annual chl a biomass is highest in the southern SoG basin west of 123◦ 30’W, and in the northern basin to the west of Texada Island despite decreasing strati-fication (figure 2.1a; Masson and Pen˜a, 2009). Water column carbonate chemistryis similar for both basins in winter (pH ∼ 7.8 and ΩA < 1 in the surface, respec-tively), however, surface pH and ΩA vary throughout the SoG during summer whenQt is large (D. Ianson, manuscript in preparation, 2014). Surface pH and ΩA nearthe plume (S < 20) are low (∼ 7.9 to 8.2 and ∼ 0.4 to 0.9, respectively) relativeto the rest of the SoG. In the southern basin north of the plume, surface pH andΩA are higher (∼ 8.2 and ∼ 1.8, respectively) and the ΩA = 1 horizon is present atabout 20 m. The highest observations of surface pH and ΩA are found in the north-ern basin, but the variability there is also highest (∼ 7.9 to 8.7 and ∼ 1.0 to 4.3,53respectively). Summer surface pH and ΩA do appear to increase from the FraserRiver plume region to the northern basin, supporting a reduction in the negativeinfluence of Qt on PP. However, highly variable observations in the northern basinpoint to the existence of additional mechanisms that control pH and ΩA that are notas strong in the southern basin.In the spring, wind appears to be the dominant driver of near-surface pH in thesouthern basin, and surface ΩA appears to be connected to wind-influenced PP. Uis strong throughout the SoG in winter, although instantaneous winds are variablebetween the northern and southern reaches of the southern basin. One feature of thesouthern reach that is absent from the north is the Squamish wind: a gap wind thatis funneled from the northern Canadian interior through narrow topography andarrives, often as a northerly gale, through Howe Sound (Thomson, 1981). Strongwind events could delay spring PP in the southern reach while allowing earlier PPto the north. An analysis of satellite imagery by Gower et al. (2013) has providedevidence that spring blooms are often ‘seeded’ by smaller blooms in the northerninlets, particularly Malaspina Inlet. However, Gower et al. (2013) also find evi-dence of elevated chl a near the Fraser River plume during some spring blooms,which is consistent with the location of maximum spring stratification in the SoG(Masson and Pen˜a, 2009). A proper climatology of winter winds in the SoG hasnot been published. However, if winter winds are stronger in the southern reachrelative to the rest of the SoG, they may occur with a stronger freshwater stratifyingeffect due to the proximity of the Fraser River. Stronger wind in the southern SoGwould then be accompanied by a water column that is more difficult to mix. Windsmight be weaker to the north but the water column would be mixed more easily.In the shallow passages that connect the SoG to its neighboring basins, tidesreplace wind as the dominant driver of vertical mixing. Neither wind nor fresh-water have a strong effect on stratification in these regions, which remains lowyear-round. Masson and Pen˜a (2009) found mean annual chl a levels to the southin Haro and Rosario Straits, and to the north near Johnstone Strait to be less thanhalf of those observed in the southern SoG basin. It follows that surface pH andΩA would be persistently lower in those passages relative to the SoG throughoutmost of the year, except in winter. Surface pH and ΩA observations are lower inHaro Strait than in the SoG (7.8 to 7.9 and 1.0 to 1.8, respectively), except near the54Fraser River plume where low surface ΩA (0.4 to 0.9) is driven by freshwater Ca2+and CO2−3 dilution (D. Ianson, manuscript in preparation, 2014). Because wind andfreshwater have a lesser impact on stratification and vertical mixing in these pas-sages, their effects on surface pH and ΩA are likely reduced as well. One exceptionmay be during peak discharge on neap tides when freshwater pulses pass throughHaro and Rosario Straits relatively unmixed (Masson and Cummins, 2004). Thesepulses could carry surface properties from the southern SoG, and thus communi-cate the effect of freshwater on pH and ΩA in the SoG to the south through HaroStrait and into the Strait of Juan de Fuca and Puget Sound.2.3.3 Comparison to other regionsThe model results are generally consistent with 2008 observations in Puget Sound.Feely et al. (2010) observed a well-mixed water column in February 2008 in thecentral basin of Puget Sound that was entirely undersaturated with respect to arag-onite. Summer surface pH observations in Puget Sound ranged from 7.77 to 8.25.The lower values tended to occur over shallow tidal mixing zones in AdmiraltyInlet, analogous to the tidal mixing and low surface pH observed in Haro Strait(Masson and Pen˜a, 2009; D. Ianson, manuscript in preparation, 2014). Water be-low 50 m in Puget Sound was not entirely undersaturated with respect to aragonitein August 2008, and these higher deep ΩA values may represent a fundamentaldifference between Puget Sound and the SoG.The strong seasonal cycles of pH and ΩA found in this study and in Feelyet al. (2010) appear to be unique to the semi-enclosed waters of the SoG and PugetSound relative to the Juan de Fuca Strait and the North American continental mar-gins. Unlike the SoG, winter surface pH in Juan de Fuca Strait is higher than insummer, and the winter ΩA = 1 horizon is significantly deeper than in summer (D.Ianson, manuscript in preparation, 2014). Surface pCO2 values measured fromunderway sampling systems in the Strait of Juan de Fuca between 2005 and 2011were supersaturated with respect to the atmospheric concentration year-round ex-cept in April (Evans et al., 2012). Low surface pH and ΩA and supersaturatedpCO2 in the Strait of Juan de Fuca reflect a year-round, weak stratification due totidal mixing near the basin head, and a resulting lack of PP (Masson and Pen˜a,552009).On the west coast of Vancouver Island, surface nutrients are supplied year-round by the Vancouver Island Coastal Current (VICC) which is fed by the buoyancy-driven estuarine outflow from Juan de Fuca Strait. Bianucci et al. (2011) deter-mined that nutrients in the VICC are responsible for enhanced PP over the shelfthat increases the upward vertical flux of DIC and alleviates ΩA < 1 in the lowerwater column. Undersaturated surface pCO2 with respect to the atmospheric con-centration in spring and summer (Evans et al., 2012; Ianson and Allen, 2002) isconsistent with the seasonal biological signal in the SoG, however periods of near-shore pCO2 supersaturation particularly near the Juan de Fuca canyon are evidenceof the summer upwelling signal (Ianson and Allen, 2002).In the California Current System (CCS), a modeling study by Hauri et al.(2013) found monthly means of surface pH for 2011 to range from 7.85 to 8.15 andmonthly means of surface ΩA for 2011 to range from 1.5 to 2.3. These ranges aresimilar to modeled surface pH and ΩA in the southern SoG during March throughNovember, however the seasonal cycle in surface pH and ΩA in the CCS was muchdifferent than the results for the SoG. The lowest surface pH and ΩA values oc-curred in May through September, which is consistent with seasonal upwelling,despite PP over the shelf. The only exception was a temperature-driven surface ΩAmaximum in the Southern California Bight during August through October.Water-column undersaturation with respect to aragonite in winter is a partic-ularly striking feature of the SoG that is not observed in the Juan de Fuca Strait,the west coast of Vancouver Island, or the CCS, although periodic surface ΩA < 1has been reported north of Cape Mendocino (Feely et al., 2008; Hauri et al., 2013).In the Strait of Juan de Fuca, the ΩA = 1 horizon can be as shallow as 8 m in thesummer but has not been observed to reach the surface. The mean winter ΩA = 1horizon depth of 50 to 60 m in the Strait of Juan de Fuca is a strong contrast to thesurface ΩA < 1 in the SoG (D. Ianson, manuscript in preparation, 2014). On thewest coast of Vancouver Island, Bianucci et al. (2011) did not find ΩA < 1 over theshelf when the VICC was present, and only removing the current from the modelled to ΩA < 1 near the sediments. Aragonite undersaturation is not absent fromthe Vancouver Island shelf, however, as ΩA data calculated from T and dissolvedoxygen demonstrate a 40% annual probability of encountering ΩA < 1 in the top5660 m (Ianson, 2013). Periods of surface ΩA < 1 in the CSS are geographicallyisolated and associated with the strongest upwelling events (i. e. summer; Feelyet al., 2008; Hauri et al., 2013).The location and relative isolation of the SoG with respect to coastal pro-cesses (e. g. upwelling and advection) combined with the strong local forcingfrom wind and river discharge are likely responsible for the carbon cycle extremes,particularly surface ΩA < 1 that occur there. The isolation from coastal advectiveprocesses and longer residence times (Pawlowicz et al., 2007), particularly in thedeep basins, allow local biological carbon fluxes to accumulate (Johannessen et al.,2008), amplifying sources of deep DIC. Local physical processes then drive an am-plified seasonal cycle of pH and ΩA relative to the Strait of Juan de Fuca and thecontinental shelf.2.3.4 LimitationsThere are three primary limitations of the model used in this study. The first, resid-ual over-productivity of phytoplankton, is clear from the model evaluation againstobservations (section 2.2.2). The second is the river flow parameterization, whichis determined by the flow record of only two rivers and TA and pH observationsfrom only the Fraser River. The third is the limited-depth, 1-D model domain,which parameterizes deeper processes as boundary conditions instead of resolvingthem.Despite rigorous tuning of a variety of parameters to the 4-year STRATOGEMdataset (appendix B), the model remains slightly overproductive during the sum-mer as indicated by surface NO−3 depletion (figure 2.6c). This excess NO−3 uptakeis coupled to DIC uptake that may cause the model to overestimate pH and ΩAby as much as 0.27 and 1.1, respectively (figure 2.8c and d). Despite this limi-tation within the biological model, pH and ΩA were not severely underestimatedduring any of the experimental runs, and only low model estimates of pH and ΩAwould generally be interpreted as unfavorable to marine calcifiers like shellfish(e. g. Kroeker et al., 2010). Thus this limitation does not produce unrealisticallyalarming results, but instead serves as a best case scenario.With the exception of the nutrient-limited, non-Redfield DIC uptake term, DIC57uptake and remineralization proceed in the model according to the Redfield ratio.However, the simple non-Redfield DIC parameterization employed here may notbe sufficient to reproduce the actual uptake ratios observed in the Strait of Georgia.Ianson et al. (2003) estimate uptake ratio uncertainties of±100% for the west coastof Vancouver Island. I did not investigate the effects of tuning this parameteriza-tion, however, since summer, near-surface model DIC is already deficient relativeto summer observations it is unlikely that increasing the strength of the excess DICuptake would improve the accuracy of model pH and ΩA.Although the Fraser River drainage basin is the largest source of freshwater tothe SoG, there are at least 7 other regional drainage basins flowing into the SoGas well (Morrison et al., 2012). While the drainage of these basins is small rela-tive to the Fraser, they may exert strong influence during the winter with increasedprecipitation or even during the summer as plumes migrate. I have attempted toaccount for the influence of these drainage basins by selecting one river (English-man) draining South Vancouver Island as a proxy. However, there are at least threecharacteristic streamflow regimes observed in rivers in the SoG and rivers withinone regime can demonstrate characteristics of another (Halverson and Fleming,2014). Furthermore, the lack of availability of carbonate chemistry data in otherSoG rivers forces me to base the DIC and TA contributions from those rivers ondata primarily from the Fraser River. Unlike the Fraser watershed which containsa significant carbonate weathering source, the smaller watersheds draining to theSoG are composed of primarily younger, Si-rich volcanic rock and are likely lowerin TA (chapter 3).The model parameterizes all processes below 40 m as boundary conditions.These boundary conditions are fit to observations at 40 m and thus allow the modelto reproduce observations in the near-surface. However, the relative importance ofDIC sources below 40 m remains unclear. DIC of upwelling origin is the dominantdriver of low pH and ΩA observations in the California Current region (Feely et al.,2008), however DIC-rich water entering the SoG through Juan de Fuca Strait isactually DIC-poor relative to SoG intermediate water on the same dnsity surfaces(D. Ianson, manuscript in preparation, 2014). Without resolving these deeper pro-cesses the SoG more explicitly, it is difficult to draw quantitative conclusions aboutlong term trends in SoG carbonate chemistry.58Chapter 3Freshwater analysis3.1 IntroductionThe pH of an estuarine system can be strongly affected by the carbonate chemistryof its freshwater sources. While this concept may seem intuitive for variations infreshwater pH, variations in freshwater dissolved inorganic carbon (DIC) and totalalkalinity (TA) at constant freshwater pH can also produce large changes in estuar-ine pH (Mook and Koene, 1975; de Mora, 1983; Whitfield and Turner, 1986; Hof-mann et al., 2009). In temperate rivers, DIC, TA and pH are seldom constant overseasonal timescales. Riverine TA is primarily a product of chemical weathering.Since weathering is generally kinetic-limited (e. g. Fraser River, Voss et al., 2014),variations in TA generally arise due to dilution (via precipitation and snowmelt) orthe confluence of tributaries. Seasonal changes in DIC and pH can be driven byprimary productivity and remineralization of terrestrial dissolved organic matter(DOM; Wallin et al., 2010).Of the freshwater sources to the Strait of Georgia (SoG), the Fraser River isthe largest, contributing ∼77% of the annual total freshwater runoff and ∼ 55 to60% between December and May (Waldichuk, 1957). The Fraser has been studiedpreviously with respect to dissolved and suspended constituents, many of whichcontribute to TA (de Mora, 1983; Voss et al., 2014; Cameron et al., 1995; Cameron,1996; Whitfield, 1983; Whitfield and Schreier, 1981). de Mora (1983) observedTA and pH in the Fraser River and in the SoG near the Fraser River mouth to vary59 128oW  126oW  124oW  122oW  120oW  118oW   48oN   50oN   52oN   54oN   56oN 200 km 30’  123oW  30’  122oW  30’  50’   49oN  10’  20’  30’ 50 kmGeologic BeltsForelandOminecaIntermontaineCoastFraser Sites1. Fitzwilliam2. McBride3. Hansard5. Lillooet6. Lytton7. Hope8. Ft. Langley9. VancouverTributary Sites10. Robson11. Bowron12. McGregor13. Willow14. Nechako15. Blackwater16. Quesnel17. Chilcotin18. Bridge19. Thompson20. Harrison21. Chilliwack22. PittEstuary Sites23. Gravesend24. Oak St.104. Stoner211134121761975720218222324 9131415 1618(b)(a)PacificOcean25. Buoy26. Model site2526Strait ofGeorgiaFigure 3.1: Maps of (a) the Fraser River watershed with geologic belts(Wheeler et al., 1991) and Voss et al. (2014) sampling sites, and (b)the Fraser Valley and estuary with Voss et al. (2014) and EnvironmentCanada (EC, 23-25) sampling sites and the model site. Hope (7) and theChilliwack River (21) are locations of EC river gauging.seasonally by as much as 500 µeq kg−1 and 0.5 pH units, respectively below asalinity (S) of 10. These fluctuations are on the order of the seasonal variability inthe surface SoG presented in chapter 2.The Fraser River follows its main course ∼1,375 km from the headwaters to60the estuary directly south of the city of Vancouver and drains ∼238,000 km2 or25% of mainland British Columbia by area (figure 3.1). The Fraser flows from itsheadwaters in the Rocky Mountains near the Alberta border to the northwest, join-ing with the Bowron and McGregor Rivers (figure 3.1-11,12), and then abruptlyturns south where it joins the Nechako River (figure 3.1-14). From the NechakoRiver confluence, the Fraser River flows south for ∼700 km joining several tribu-taries and cutting its way into the narrow Fraser Canyon near the confluence withthe Chilcotin River (figure 3.1-17). The Fraser joins its principle tributary (fig-ure 3.1-19), the Thompson River, along the canyon, and turns west at the town ofHope (figure 3.1-7) and into the Fraser Valley alluvial plain. The Fraser estuary isconsidered to begin near the city of New Westminster (figure 3.1-9), although theseaward tip of Annacis Island (figure 3.1-23) is the farthest upstream that the saltwedge has been observed.Along its course, the Fraser River passes through four geologic belts: the Fore-land, Omineca, Intermontaine, and Coast belts in order from east to west and byage. The four belts represent 6%, 22%, 58%, and 14% of the Fraser Basin by area,respectively (Cameron and Hattori, 1997). The Foreland belt, which encompassesthe Fraser headwaters and borders the first ∼ 500 km of the river, is composed pri-marily of carbonate-rich sedimentary rock while the other three belts have strongcomponents of silicate-rich mantle and crustal igneous rocks (Cameron and Hat-tori, 1997). Ca/Na and Mg/Na ratios throughout the river demonstrate a character-istic mixing line between carbonate weathering and silicate weather endmembers,with a shift to silicate weathering beginning downstream of the McGregor River(Voss et al., 2014). Fraser River TA is dominated by carbonate alkalinity (CA),and thus carbonate weathering in the Foreland Belt is the primary source of TA tothe river (Voss et al., 2014).Despite hydroelectric projects along select tributaries (e. g. Nechako and BridgeRivers), the main river course remains undammed. With the exception of MooseLake near the headwaters (figure 3.1-1), the river flows unrestricted to the SoG. Thelack of dams and lakes along the river and the cold winter climate of the moun-tainous basin produce a strong seasonality to the drainage volume of the river.Background flow is low (500 to 2,500 m3 s−1) in winter and increases to peakflow (4,000 to 15,000 m3 s−1) marked by a characteristic freshet occurring typi-61Q F [103  m3  s−1 ]  0246810121416MeanSt DevMin/MaxQ C [103  m3  s−1 ]J F M A M J J A S O N D J00. (b) Chilliwack(a) FraserFigure 3.2: Hydrographs showing mean (solid black), standard deviation(dark gray) and minimum/maximum (light gray) observed discharge for(a) the Fraser River at Hope, BC and (b) the Chilliwack River at Chill-iwack, BC over the 1913-2013 Environment Canada (EC) daily record( in June (figure 3.2a). Flow has been historically gauged at Hope, howeverdownstream tributaries (e. g. the Chilliwack, Harrison, and Pitt Rivers, figure 3.1-21,20,22) can make significant contributions to flow along this ∼ 150 km path tothe estuary (e. g. Chilliwack, figure 3.2b). Pawlowicz et al. (2007) estimated flowat Port Mann (figure 3.1-9) near the mouth to be∼ 20% higher than at Hope duringthe summer freshet and as much as double during winter when rainfall contribu-tions are large.Currently available observations and literature estimates of Fraser River TA62and pH are either insufficient or too inconsistent to provide a seasonal summary ofcarbonate chemistry near the mouth. While several studies in the Fraser River re-port TA and pH measurements (Voss et al., 2014; Cameron et al., 1995; Cameron,1996), de Mora (1983) has been the only study to discuss the two quantities indetail. However, this discussion focuses primarily on freshwater/seawater mix-ing and does not elaborate on the variability of freshwater TA and pH, of whichonly four sampling dates are presented. TA measurements reported in other stud-ies are generally consistent to within the variability observed by de Mora (1983),however, pH throughout the river is severely underestimated (< 7.0) by Cameronet al. (1995) relative to the other studies. Unpublished TA and pH monitoring pro-grams in the Fraser River also exist; both quantities were sampled semi-monthlyby Environment Canada at Hope between 1979 and 1999 (table 3.1). However,DIC sampling overlapped TA and pH in 1998, and pH calculated from DIC andTA using CO2SYS and the Millero (1979) freshwater K1 and K2 constants doesnot reproduce the observed pH (error > 1.0 units), indicating significant error inone or more of the three observed quantities. DIC observations in the Fraser Riverare generally scarce - the most comprehensive DIC sampling program to date wascarried out by Environment Canada at Hope between 1998 and 2006. The lack ofrobust DIC datasets and the high sensitivity of the carbonate system equilibrium toDIC make TA and pH more reliable indicators of river carbonate chemsitry.In order to fully describe the seasonal variability of TA and pH of river waterentering the estuary, an analysis of accurate, higher frequency datasets is needed.However, inherent uncertainties associated with measuring TA and pH in estuarineenvironments will complicate such an analysis. An ideal study would use multi-ple datasets to provide a range of variability with which to bound the uncertaintyin each measurement. The aims of this chapter are to (1) identify the range andvariability of TA and pH in the Fraser River freshwater contribution to the SoGusing data from recent sampling programs, (2) assess the sensitivity of southernsurface SoG pH and ΩA to (1) using the model described in chapter 2, and (3) de-termine realistic parametrizations of freshwater pH and TA for use in the chapter 2modeling study based on the outcome of (2). Assessing SoG surface pH and ΩAsensitivity to freshwater TA and pH is essential to providing error estimates for themodel, and those estimates are included in the outcomes of this chapter. I base my63analysis on observations from the Fraser River since it is the largest contributor offreshwater by volume to the southern SoG year-round, and is the most extensivelysampled with respect to TA and pH. I have chosen to focus this analysis on TA andpH rather than DIC as the two former quantities are more consistently sampled inthe lower Fraser River and more robust carbonate system metrics.3.2 Data3.2.1 SourcesThere are three published studies that report TA and pH in the Fraser River andFraser estuary (de Mora, 1983; Cameron et al., 1995; Cameron, 1996) and one ad-ditional study that reports TA only (Voss et al., 2014). While the spatial coverageof these studies is rich (headwaters to Vancouver in Voss et al. (2014), Cameronet al. (1995) and Cameron (1996); Vancouver to SoG in de Mora (1983); fig-ure 3.1), the temporal coverage is poor (maximum four sampling dates), althoughCameron (1996) does provide five-year averages from a higher frequency Environ-ment Canada dataset.The present study makes use of data from several sampling programs in theFraser River and southern SoG in addition to the literature studies. TA datasetsfrom monitoring programs at Hope, BC (EC Hope) and in the Fraser estuary (ECestuary) along the north and main arms were acquired from Environment Canadaalong with hourly pH data at a mooring in the Fraser estuary (EC buoy; table 3.1).TA measurements in the Fraser estuary and the southern SoG were obtained froma published study by researchers at the University of British Columbia (UBC es-tuary), and from several cruises conducted by Fisheries and Oceans Canada at theInstitute of Ocean Sciences (IOS estuary, IOS carbon; table 3.1).3.2.2 UncertaintyThe purpose of measuring TA as a carbonate system variable is ultimately to esti-mate CA. In seawater, TA is defined in terms of the proton donors and acceptors(including CA) present in high enough concentrations to influence the result be-yond the precision limit (Dickson, 1981). Most species are easily measured or64Dataset Location(s) Years Frequency Qty MethodEC Hope∗ 49.387 N 121.451 W 1979 - 1999 semi-monthly TA Environment Canada, 2006EC estuary∗ 49.202 N 123.122 W49.147 N 123.033 W2004 - 2009 semi-monthly TA Environment Canada, 2006EC buoy∗ 49.148 N 123.039 W 2008 - 2013 hourly pH Ethier and Bedard, 2007UBC estuary∗∗ Fraser and SoGsee de Mora, 19811978 - 1979 seasonal TA de Mora, 1983IOS estuary† 49.165 N 123.547 W -49.113 N 123.277 W2010 twice TA Dickson et al., 2007IOS carbon‡ see chapter 2, table 2.1 2010 - 2012 seasonal TA Dickson et al., 2007∗Environment Canada (∗∗de Mora, 1981†R. W. MacDonald and S. Johannessen, pers. com., 2014‡D. Ianson, manuscript in preparation, 2014Table 3.1: Datasets used to assess the range and variability of freshwater TAand pH near the Fraser River estuary.approximated (Lewis and Wallace, 1998), leaving CA to be calculated as an un-known. In many freshwater systems, dissolved constituents disregarded by theDickson (1981) TA definition, especially organics which are difficult to measure,occur in high enough concentrations to significantly increase or reduce the TA(Hunt et al., 2011). Unknown organics have been shown to influence TA in coastalzones as well, particularly during strong phytoplankton bloom conditions (Koeveand Oschlies, 2012; Kim and Lee, 2009; Herna´ndez-Ayo´n et al., 2007; Cai et al.,1998). The presence of these unknown, possibly organic, contributions to TA in-troduces uncertainty to the CA estimate of coastal/estuarine systems like the Fraserestuary and SoG.It was mentioned earlier that pH calculated from 1998 observations of DIC andTA (CO2SYS; Millero, 1979) at Hope, BC (EC Hope, table 3.1) failed to repro-duce observed pH. The inconsistencies between DIC, TA, and pH observations inthe EC Hope dataset could be partially explained by inaccurate CA estimates dueto the presence of unknown, possibly organic, river constituents. However, errorcould also be introduced during sample treatment as chemical preservation is notspecified in the standard operating procedures for DIC and TA sample collection(Environment Canada, 2005, 2006). Chemical preservation using HgCl2 has be-65come an oceanographic standard with respect to DIC and TA sampling (e. g. Dick-son et al., 2007; Dickson and Goyet, 1994) in order to restrict biological processesfrom influencing the two quantities after the sample is taken.Measurements of Fraser River pH are also problematic in the literature. Cameronet al. (1995) calculated pCO2 values throughout the river basin to be > 5,000 ppmusing potentiometric pH observations from a 1993 study. However, these pH val-ues are consistently low, as much as 0.5 pH units less than the 5 year averagesreported by Cameron (1996). Nearly all pH measurements from the Fraser Riverare potentiometric. It has been argued that precision of potentiometric pH measure-ments is generally strong enough for estuarine work (0.01 units; Dickson, 1993),however the instability of the liquid-junction potential makes these measurementsinherently uncertain (Covington et al., 1983) and accuracy can be easily compro-mised by instrument drift (UNESCO, 1988). Additionally, accurate measurementsrequire calibration in buffers of a similar ionic strength, and if the ionic strengthof the in situ medium changes, measurements can become less accurate (Dickson,1993).Despite inconsistent methods and potentially large uncertainty associated withthe described datasets (table 3.1), I assume that a comprehensive analysis of theavailable data will produce a range of plausible scenarios for each of the two quan-tities in the Fraser River outflow. I further assume that realistic scenarios for eachquantity will lie within that range. The sensitivity of the SoG carbonate chem-istry to the prescribed range will determine how important these uncertainties arein defining the SoG freshwater chemistry for the model in chapter AnalysisData from Voss et al. (2014) demonstrate the variability of TA throughout theFraser River basin across three different flow regimes (figure 3.3). TA accumu-lates along the main course from the headwaters to Hansard near the confluenceswith the McGregor and Bowron rivers (∼ 450 km downstream, figure 3.1). October2010 sampling indicates a gradual decrease in TA downstream of Hansard towardthe estuary, which is consistent with a shift between carbonate and silicate weath-ering (Voss et al., 2014). However, this decrease is not observed during summer66sampling when flow stage is higher. TA in the Fraser River decreases with increas-ing discharge between the low and mid-stage flow regimes (figure 3.3a). Largedecreases in TA with discharge are also evident in the Bowron, McGregor, Willow,Blackwater, and Bridge Rivers while most of the remaining tributaries maintainapproximately constant TA (figure 3.3b). TA appears to increase significantly withdischarge in only two of the tributaries (Robson and Chilcotin Rivers) and onlybetween mid and high-stage flow.1 2 3 4 6 8 95 7TA [103  µmol L−1 ]00.511.52101112 131415 16 17 18 19 20 22TA [103  µmol L−1 ]Distance from Fraser River source [km]  0 200 400 600 800 1000 1200 140000.511.52Oct 2010Jul/Aug 2009May/Jun 2011Fraser Sites1. Fitzwilliam2. McBride3. Hansard4. Stoner5. Lillooet6. Lytton7. Hope8. Ft. Langley9. VancouverTributary Sites10. Robson11. Bowron12. McGregor13. Willow14. Nechako15. Blackwater16. Quesnel17. Chilcotin18. Bridge19. Thompson20. Harrison22. Pitt(a)(b)Figure 3.3: TA at various locations (a) along the Fraser River and (b) inselect Fraser River tribuaries (figure 3.1) sampled at low-stage flow(1719 m3 s−1) between 14 to 27 October 2010 (black), mid-stage flow(3616 m3 s−1) between 28 July to 13 August 2009 (gray), and high-stage flow (9040 m3 s−1) between 26 May to 7 June 2011 (white).TA data are from Voss et al. (2014). Fraser River flow is calculatedfrom Environment Canada observations at the Hope, BC gauging sta-tion ( seasonal pattern in Fraser River TA can be observed in the EC Hope dataset(figure 3.4a). The seasonal pattern in the estuary is less clear (figure 3.4b), but thevariability is still large (∼700 to 1150 µmol kg−1). TA differences between the67TA [103  µeq kg−1 ]  (a)1980 1982 1984 1986 1988 1990 1992 1994 1996 19980.40.60.811.21.4HopeTA [103  µeq kg−1 ]  (b)2004 2005 2006 2007 2008 20090.40.60.811.21.4EstuarypH [NBS scale]  (c)2008 2009 2010 2011 2012 20137. 3.4: Timeseries of Fraser River TA at (a) Hope, BC (EC Hope) and(b) the Fraser Estuary (EC estuary) and (c) Fraser River pH at theEnvironment Canada Water Quality Buoy in the Fraser estuary (ECbuoy). White/gray patches span the duration between two consecu-tive freshets as recorded at Hope, BC (Environment Canada, Date ranges differ between thethree subplots. Errorbars in (b) represent standard error across replicates(n = 5 before 1 August 2007, n = 2 after). Daily pH values in (c) arebinned from hourly data with a mean standard deviation of ∼ 0.04. De-tails regarding each dataset are provided in table 3.1.estuary and Hope (∼200 µmol kg−1) are on the order of the change between Hopeand the upstream stations observed in the Voss et al. (2014) data (figure 3.3a). How-ever, the reach between Hope and the estuary is far downstream of the general shiftto silicate weathering identified by Voss et al. (2014). Instead, local processes arelikely responsible for either removing or diluting TA. River confluences seaward68of Hope may reduce TA through dilution, especially considering the large (17%to 50%) seasonal flow contributions of downstream tributaries (Pawlowicz et al.,2007). Observations of TA in the lower tributaries are sparse, however the PittRiver demonstrated markedly low TA (< 100 µeq L−1) between 26 May and 6 July2011 (Voss et al., 2014; figure 3.3). The Fraser River is rich in suspended sediments(Whitfield and Schreier, 1981), many of which could contribute to TA by cation ex-change of TA constituents at charged sites on particulate surfaces (Kennedy, 1965).de Mora (1983) suggested that settling of these sediments may remove TA near theFraser River mouth. Sedimentation is also likely seaward of Hope since the riverhas lost most of its elevation by that point, and thus sedimentation may removeTA between Hope and the estuary. Precipitation of CaCO3 may be responsible forremoving carbonate in systems with a high Ω, however, this process is unlikely inthe Fraser River since Ω is low throughout (de Mora, 1983).A third freshwater TA timeseries can be constructed from the IOS carbon, IOSestuary and UBC estuary datasets that cover the Fraser estuary and nearby SoGregions (S > 0) by extrapolating the S = 0 TA endmember from each cruise (fig-ure 3.5). For this analysis I only consider cruises that have at least one TA samplebelow S = 20. Each extrapolated endmember is an estimate of the combined TAfrom all of the freshwater sources to the sampling area, but because of its closeproximity and magnitude the Fraser River likely still dominates the signal. Theset of ten extrapolated endmembers range from ∼550 to 1150 µeq kg−1 which islower than the EC Hope data range but near the EC estuary data range. Strong,non-conservative TA behavior at low S was observed in four of the five cruisesfrom the UBC estuary dataset (de Mora, 1983), suggesting that extrapolated fresh-water TA may underestimate true freshwater TA. Since it is not clear whether thisnon-conservative behavior is associated with the CA or organic component of TA,I use the extrapolated endmembers for this analysis and consider TA to be approx-imately conservative in the UBC estuary, IOS carbon, and IOS estuary datasets.pH at the Environment Canada mooring in the Fraser estuary (EC Buoy) variesseasonally (figure 3.4c). There is no relationship between buoy pH and conduc-tivity (data not shown) which suggests that the seasonal pattern is an attribute ofthe river and not an influence of the southern SoG during weak discharge or strongtides. Buoy pH demonstrates a weakly hyperbolic relationship to river discharge690 5 10 15 20 25 30 350.40.60.811. [PSS−78]TA [103 µeq kg−1 ]  1978−041978−071978−161979−011979−072010−572010−732011−092011−602012−57Figure 3.5: TA versus S regressions from observations on several cruises inthe southern SoG and Fraser River estuary. Cruise ID numbers beginwith the year of sampling. Cruises prior to 2000 are from the UBC es-tuary dataset, 2010 cruises are from the IOS estuary dataset, and cruisespost-2010 are from the IOS carbon dataset (table 3.1). Each cruise plot-ted contains at least one datapoint at S < 20. Observations at S = 0generally deviate from conservative behavior (e. g. de Mora, 1983) andare excluded from these regressions.calculated at Port Mann. However, the weak correlation suggests that pH variabil-ity in the lower Fraser River is primarily due to local processes that influence DIC.Primary productivity, specifically, would be a classic driver of such variability, andthe cycle of buoy pH is consistent with the biologically-driven cycle modeled inchapter 2 (e. g. low in winter, high in summer). The buoy pH data demonstrate amean of∼7.7 with a 0.17 unit standard deviation and approximate 95% confidenceintervals at 7.4 and 8.0.The Voss et al. (2014) data suggest that flow rate is the dominant driver ofchanges in Fraser River TA. Voss et al. (2014) also observed significant anticor-relations to discharge for several elements in the Fraser River (Ca, Mg, Na, Cl,K, Sr, Ba, also SO2−4 ) that they sampled more frequently than TA. They founda rapid decline in concentration of each constituent with discharge near base flowrelative to peak flow in these relationships consistent with primarily kinetic-limited700.40.60.811.21.4 y = 1.1670x−0.1571pa = 0.000 pb = 0.000r2 = 0.589(a) HopeTA [103  µeq kg−1 ] y = −0.0241x + 0.9325p1 = 0.003 p0 = 0.000r2 = 0.090mean = 0.87(b) EstuaryTA [103  µeq kg−1 ]0 2 4 6 8 10 y = 0.0125x2 − 0.1653x + 1.1160p2 = 0.041 p1 = 0.011 p0 = 0.000r2 = 0.757(c) SoGTA [103  µeq kg−1 ]Qt [backaveraged over 15 days; 103 m3 s−1]  1978−041978−071978−161979−011979−072010−572010−732011−092011−602012−57Figure 3.6: TA versus 15-day backaveraged parameterized total discharge(Qt, section 2.1.4) for (a) Hope, BC (EC Hope, table 3.1), (b) the Fraserestuary (EC estuary, table 3.1), and (c) the SoG (extrapolated to S = 0,figure 3.5). Errorbars in (b) represent standard error across replicates(n = 5 before 1 August 2007, n = 2 after). Errorbars in (c) representstandard error of the S = 0 intercept calculated using linear regression.Solid symbols in (c) are the observed TA at S = 0, where available, fromthe corresponding open symbol cruises. Solid symbol errorbars are thestandard error across replicates (n = 2).71weathering followed by modest dilution. These observed relationships and the TAdata in figure 3.3 suggest that chemical weathering rates of TA in the Fraser Basinare likely flow-dependent with higher rates associated with low base flow levels.However, relationships to discharge are not always clear, even in the presence ofan overall significant correlation. Variability due to hysteresis, or a lag in tracerflux that produces differing discharge relationships between rising and falling flowstages, can complicate the relationship to discharge for some river constituents,as demonstrated by Voss et al. (2014) for PO3−4 and SiO2 and by Whitfield andSchreier (1981) for sediment concentration.The EC Hope and EC estuary TA datasets and the ten extrapolated freshwa-ter TA datapoints demonstrate significant negative correlations to discharge (fig-ure 3.6) similar to those observed for Ca, Mg, Na, Cl, K, Sr, Ba, and SO2−4 byVoss et al. (2014). The EC Hope data are best described by a power fit with thelow discharge tail removed (figure 3.6a) and the EC estuary data are best describedby a linear function (figure 3.6b). Interestingly, the extrapolated freshwater TAvalues increase with discharge when Qt > 7,000 m3 s−1, and are best describedby a second order polynomial with the tails removed (figure 3.6c). These regres-sions provide a high-resolution summary of the seasonal variability and extremesof Fraser River TA downstream of Hope, and their shapes are consistent withthe kinetic-limited weathering and dilution mechanisms identified by Voss et al.(2014). The three regressions demonstrate a total freshwater TA range of ∼500 to1,350 µeq kg−1.3.3 Sensitivity3.3.1 ScenariosFreshwater TAThe ranges of freshwater TA and pH and the potential relationships between TAvariability and discharge established from the described datasets (table 3.1) providethe basis for a set of reasonable scenarios with which to test the sensitivity ofsurface SoG pH and ΩA. As seasonal variations in freshwater TA are apparent72from the available data (figure 3.4), I use the regressions between the EC Hope, ECestuary, and extrapolated freshwater TA datasets and discharge (figure 3.6) as thethree flow-dependent freshwater TA scenarios (table 3.2). However, SoG surfacepH and ΩA may not be as strongly sensitive to the seasonality of freshwater TAas they are to the overall mean value. To test this case, I use constant freshwaterTA scenarios in addition to the flow-dependent scenarios. I propose the minimumand maximum observed freshwater TA values, 500 µeq L−1 and 1,350 µeq L−1respectively (figure 3.6a and b), and the mean value from the EC estuary dataset(870 µeq L−1, figure 3.6c) as an additional three plausible constant freshwater TAscenarios for a total of six scenarios (table 3.2).Scenario TA f resh (µeq kg−1) Qt limit Dataset(s)1 Minimum 500 UBC estuary, IOS estuary, IOS carbon2 Quadratic 125000Q2t −0.1653Qt +1116.0 > 11,000 UBC estuary, IOS estuary, IOS carbon3 Linear −0.0241Qt +932.5 EC estuary4 Mean 870 EC estuary5 Power 3454.4Q−0.1571t < 400 EC Hope6 Maximum 1350 EC HopeTable 3.2: SoG freshwater TA scenarios based on fits between Fraser Riverand SoG TA data and 15-day backaveraged total parameterized freshwa-ter discharge, Qt (figure 3.6). The Qt limit designates the range of Qt (m3s−1) where the scenario is invalid. TA f resh from the closest end of thevalid Qt range is used when the invalid range is reached.Freshwater pHThe EC buoy pH data follow a seasonal pattern (figure 3.4c). However, since theydo not demonstrate a significant correlation to discharge, I only consider constantfreshwater pH scenarios for the southern SoG, and use the approximate 95% con-fidence intervals of the buoy pH to represent the range of seasonal variability. Ipropose low (7.4), medium (7.7), and high (8.0) freshwater pH scenarios, wherethe medium scenario is the approximate mean buoy pH and the low and high sce-narios reflect the confidence intervals.73Freshwater DICFor simplicity, I assume the Fraser River to behave similarly to an S = 0 seawaterendmember. Part of this assumption means using the Dickson (1981) definition ofseawater TA, and thus regarding freshwater TA as essentially carbonate alkalinity(CA). Considering the data caveats discussed earlier (section 3.2.2), this assump-tion is not accurate. However, since quantifying the unknown contribution to TAin the Fraser River is beyond the scope of this study, my treatment of the FraserRiver as S = 0 seawater is necessary in order to assess the range of freshwater-driven uncertainty in marine carbonate chemistry forecasts for the SoG. Since Ihave made this assumption, I calculated freshwater DIC from pH and TA usingCO2SYS (Lewis and Wallace, 1998) and the K1 and K2 constants of (Millero,2010).3.3.2 Model and evaluationIn order to assess the sensitivity of near-surface pH and ΩA in the southern SoG tothe range of proposed freshwater TA and pH scenarios, I use the model described inchapter 2 and the base model run parameters summarized in section 2.2.3. I carriedout each set of base runs separately for the six freshwater TA scenarios (table 3.2)and three pH scenarios (7.4, 7.7, 8.0) for a total of 216 runs.I define three model metrics to assess the sensitivity of the model to freshwaterchemistry: 3 m averaged DIC (DIC3m), 3 m averaged pH (pH3m), and surface arag-onite undersaturation (ΩAsur < 1) duration. I normalized DIC3m to TA′= 2,000 µeqkg−1nDIC3m = DIC3mTA′TA3m(3.1)to show DIC trends independent of changes to river TA. As in chapter 2, I evaluateeach metric at the near-surface (3 m average) because the seasonality of pH andDIC is surface-intensified. I chose two 180-day averaging periods for nDIC3mand pH3m: one beginning 1 May and the other beginning 31 October. I considerΩAsur < 1 duration for each of the two periods as well. The two periods generallyrepresent summer and winter conditions, respectively, although the 180-day periodbeginning 31 October includes the spring bloom season.743.3.3 ResultsNot surprisingly, the strongest sensitivity occurs between 1 May and 31 Octo-ber, when freshwater discharge is highest (figure 3.7). Despite an increase in thefreshwater pH buffering capacity, raising the TA f resh appears to reduce pH3m (fig-ure 3.7c). This result is consistent with an increase in nDIC3m (figure 3.7a). Thesetrends are visible for the period from 31 October to 1 May as well (figure 3.7b andd) but with a smaller range of variability.Despite rising nDIC3m and declining pH3m, ΩAsur < 1 duration appears to de-cline with overall increasing TA f resh across the six variable and constant scenarios(figure 3.7e and f). The decline is particularly strong for TA f resh scenario 6 (con-stant maximum) at pH f resh 7.7 and 8.0 for both 180 day periods. In the case ofpH f resh 8.0 and TA f resh scenario 6 for the period beginning 1 May, the numberof runs exhibiting any duration of ΩAsur < 1 decreased from four to one of twelverelative to the other TA f resh scenarios at high pH f resh.The importance of flow dependence is most evident between scenarios 5 (power)and 6 (constant maximum) determined from the EC Hope dataset. Scenario 6produces notably lower pH3m (∼ 0.5 unit median reduction at pH f resh = 7.4, fig-ure 3.7c) and summer ΩAsur < 1 (nearly complete surface supersaturation, fig-ure 3.7e) compared to scenario 5. Scenario 6 also produces lower winter ΩAsur < 1(∼ 10-day reduction, figure 3.7f) relative to scenario 5 at pH f resh > 7.4. The differ-ence between model outcomes for scenarios 3 and 4 (linear fit and constant mean)is weak by comparison. However, an exception lies at pH f resh = 8.0 where sum-mer ΩAsur < 1 is reduced by ∼ 10 to 20 days in scenario 4 relative to scenario 3(figure 3.7e). There does not appear to be a strong difference in model outcomesbetween TA f resh scenarios 1 and 2 (constant minimum and quadratic fit).Overall, the sensitivity of pH3m and nDIC3m to the range of pH f resh is approx-imately equal to the TA f resh range sensitivity. However, decreasing the pH f reshdoes increase the sensitivity of pH3m and nDIC3m to TA f resh, especially in sum-mer. At a pH f resh 7.4, the range of pH3m and nDIC3m sensitivity to TA f resh hasnearly doubled relative to pH f resh 8.0. Changes in pH f resh have the opposite effecton ΩAsur < 1 duration. Lowering pH f resh appears to stabilize ΩAsur < 1 durationwhereas raising pH f resh strengthens the sensitivity to TA f resh.75nDIC [103  µmol kg−1 ]1 May to 31 Oct(a)pHfresh = 7.4 pHfresh = 7.7 pHfresh = 8.01.751.81.851.931 Oct to 1 May(b)1.751.81.851.9pH [total scale](c)7.988.18.28.3(d)7.988.18.28.31 2 3 4 5 6 1 2 3 4 5 6 1 2 3 4 5 6Ω A < 1 duration [days]5 5 5 5 5 5 4 4 4 4 4 4 4 4 4 4 4 1(e)02040601 2 3 4 5 6 1 2 3 4 5 6 1 2 3 4 5 6(f)120140160180Figure 3.7: Model sensitivity results for (a, b) nDIC3m, (c, d) pH3m, and (e,f) ΩAsur < 1 duration averaged over (a, c, e) 1 May to 31 October and(b, d, f) 31 October to 1 May. Each box indicates the median, standarddeviation about the mean, and confidence intervals across 12 runs (years2001 through 2012) at the indicated TA f resh (horizontal axis, table 3.2)and pH f resh (grey/white sections) scenarios, except in (e) where eachbox only contains runs that exhibited ΩAsur < 1 in summer (total labeledabove each box). TA f resh and pH f resh scenarios are in ascending orderof mean annual value from left to right.3.4 Discussion3.4.1 MixingHow do uncertainties in freshwater pH and TA affect pH in an estuary? Considerthe simple case of conservative mixing between model freshwater ( f resh) and sea-water from the 40 m boundary (sea). In the absence of additional processes, phys-ical dilution alone will drive linear mixing between high and low concentrations,76and thus the mixing curves (mix) for DIC, TA, T , Si, and PO3−4 will be linear withrespect to S (e. g. DIC and TA, figure 3.8a). pH can be calculated from thesequantities using the method for DIC calculation described in section 3.3.1. Cal-culated pH is not linear with S and demonstrates a characteristic minimum at lowS when pH f resh <= pHsea (figure 3.8b). Similar theoretical pH versus S curveshave been calculated by Mook and Koene (1975), de Mora (1983), Whitfield andTurner (1986), and Hofmann et al. (2009). Similar curves have also been observedto occur in the Scheldt River estuary on the Netherlands coast (Mook and Koene,1975; Hofmann et al., 2009) and locally in the Fraser estuary (de Mora, 1983).The pH versus S relationship that results from conservative mixing is clearlysensitive to pH f resh (figure 3.8b). However, despite sharing a common pH f resh, theTA f resh scenarios also demonstrate variability along the pH mixing curve. Specif-ically, increasing TA f resh from 500.0 µeq kg−1 to 1300.0 µeq kg−1 lowers pHmixand shifts the pHmix minimum, where present, to higher S. This behavior can be ex-plained by the difference between the DIC:TA ratio of the high TA f resh case and thelow TA f resh case along the mixing curve. High TA f resh produces a greater DIC:TAratio along the mixing curve than that produced by low TA f resh (figure 3.8c), caus-ing pHmix to be lower for high TA f resh than for low TA f resh. Similar behaviorwas demonstrated by Mook and Koene (1975) by increasing the ratio of TA f resh toTAsea.The strongest sensitivity of pHmix to TA f resh occurs at low pH f resh (pH f resh <pHsea, figure 3.8b). The sensitivity of pHmix to TA f resh at pH f resh > pHsea is rela-tively weak by comparison. The difference in sensitivity is again explained by thedifference in the DIC:TA ratios of the TA f resh cases along the mixing curve. AspH f resh increases, the freshwater DIC:TA ratio approaches the seawater DIC:TAratio and the difference between the low and high TA f resh DIC:TA ratios along themixing curve is small relative to the difference in ratios for the low pH f resh case(figure 3.8c). The difference in pHmix between the two TA f resh cases is thus smallerfor the high pH f resh case than for the low pH f resh case.The model sensitivity studies demonstrate the impact of freshwater carbonatechemistry on the marine environment beyond the simple bulk mixing case. Air-sea gas exchange, biological activity, KPP mixing, and estuarine circulation makethese studies more realistic, yet also more complicated. Hofmann et al. (2009)770.511.52103  µmol kg−1 DICTADICTA pHfresh = 7.4(a) [total scale](b)pHfresh = 7.4pHfresh = 7.7pHfresh = 8.00 5 10 15 20 25 3011.1Sp [PSS−78]DIC:TA(c)  TAfresh = 500.0 µeq kg−1TAfresh = 1300.0 µeq kg−1Figure 3.8: Theoretical two-endmember mixing curves of (a) DIC and TA,(b) pH, and (c) DIC:TA for TA f resh of 500.0 µeq kg−1 (solid) and1300.0 µeq kg−1 (broken) at pH f resh of 7.4 (black), 7.7 (dark grey),and 8.0 (light grey). Freshwater and seawater endmember DIC, TA,T , S, Si, and PO3−4 , unless otherwise specified, are defined accordingto the average freshwater and 40 m boundary values, respectively, de-scribed for the model in chapter 2. Mixing is conservative (linear) forthose quantities. pH and DIC f resh are calculated using CO2SYS (Lewisand Wallace, 1998) and the K1 and K2 dissociation constants of Millero(2010).found that adding a biogeochemical model to a prescribed physical mixing scenarioin the Scheldt estuary decreased the low S pH minimum and increased pH at higher78S. Biological activity created a net decrease in pH (net remineralization) while air-sea gas exchange resulted in a net pH increase (net outgassing). The 1-D modelI have used is positioned too far away from the Fraser River mouth to resolve theestuarine dynamics along the full S range, therefore I cannot make a comparisonsimilar to Hofmann et al. (2009) between my model results and the simple mixingcase.Despite the lack of a full S range, I do observe similar behavior in the mixingonly case and in the model sensitivity results. In both cases, increasing TA f reshreduces the estuarine pH despite a constant pH f resh (pHmix in figure 3.8b, pH3min figures 3.7c and d). The estuarine pH reduction is non-intuitive since TA isgenerally associated with higher pH. However, this result occurs because DICmix(DIC3m in model) increases relative to TAmix (TA3m in model) when TA f resh is in-creased (figure 3.8a). More generally, DIC f resh is a stronger driver of model pHthan TA f resh at constant pH f resh. This result is consistent with the high freshwa-ter DIC:TA ratio concept for the simple mixing case, where increasing TA f resh atconstant pH f resh increases the DIC:TA ratio along the mixing curve (figure 3.8c).The results in figures 3.7c and d also support the increased sensitivity of pHsea toTA f resh at low pH f resh (pH f resh < pHsea) demonstrated in the simple mixing case(figure 3.8b). This result again highlights the importance of the freshwater DIC:TAratio. Reducing pH f resh increases the freshwater DIC:TA ratio making the ratio be-tween TAmix and DICmix more sensitive to changes in TA f resh (figure 3.8c).The decline of the ΩAsur < 1 duration with decreasing pH3m (figure 3.7c throughf) is counterintuitive since model pH and ΩA are generally coupled through pH-driven changes in [CO2−3 ]. However, the response of [CO2−3 ] in this study is moresensitive to changes in total DIC than shifts in the equilibrium point of the carbon-ate system. Although low pH shifts the carbonate equilibrium toward [HCO−3 ], theadded [CO2−3 ] from increasing total DIC is enough to create a net increase in ΩAand decrease in ΩAsur < 1 duration. More generally, although increases in DIC f reshreduce model pH at constant pH f resh, the increase in model [CO2−3 ] drives increasesin ΩA.The increase in [CO2−3 ] caused by increasing DIC has additional implicationsfor the ΩA sensitivity to pH f resh. For the low pH f resh scenario, the combined ef-fect of increasing [CO2−3 ] and the seawater DIC:TA ratio (figure 3.8c) stabilizes79the ΩAsur < 1 duration to changes in TA f resh. As pH f resh increases, the seawaterDIC:TA ratio approaches 1 and changes in [CO2−3 ] dominate. The result is to cre-ate a stronger sensitivity of the ΩAsur < 1 duration to changes in TA f resh at higherpH f resh (figure 3.7e and f). On the other hand, pH3m and nDIC3m are most sensitiveto TA f resh at low pH f resh and least sensitive at high pH f resh (figure 3.7a through d).3.4.2 Parameterization for modelAcross all scenarios proposed here (table 3.2; 216 total runs), nDIC3m averagedover the 180-day period beginning 1 May varies by > 50 µmol kg−1 (figure 3.7aand b). This variation amounts to a∼ 50 µmol kg−1 change in DIC in the top 3 m atTA′3m = 2,000 µeq kg−1, and a 0.1 to 0.2 unit change in pH. These effects will oftenbe stronger for shorter averaging periods. The summer ΩAsur < 1 duration rangesfrom a maximum of ∼ 50 days to a minimum of ∼ 1 day across all scenarios, andthe total number of runs exhibiting any period of summer ΩAsur < 1 ranges from 5to 1 (figure 3.7e). At high pH f resh and TA f resh, the summer ΩAsur < 1 is virtuallynon-existent (1 day duration in only 1 run).The range of winter ΩAsur < 1 duration is ∼ 20 days across all scenarios (fig-ure 3.7f). Considering winter is a time of relatively weak freshwater influence, thestrength of the winter ΩAsur < 1 duration trend is striking. However, the winter av-eraging period includes the month of April when flow stage begins to increase dra-matically. Additionally, winter rainfall can produce large freshwater pulses whichmay strengthen the river influence on ΩAsur . Furthermore, the winter averagingperiod includes the spring bloom season which likely accounts for much of thewinter ΩAsur < 1 variability, although the spring bloom timing is unlikely to affectthe trend in winter ΩAsur < 1 duration with TA f resh since model phytoplankton arenot chemically-stressed and TA and DIC do not influence model growth rates.The multiple effects of pH f resh complicate the selection of a suitable freshwaterpH value for the model. At the low end of pH f resh, the ΩAsur < 1 duration is rela-tively stable over the range of TA f resh scenarios, but nDIC3m and pH3m are at theirmost sensitive to TA f resh (figure 3.7). At the high end of pH f resh it is the ΩAsur < 1duration that is most sensitive to TA f resh.The TA f resh scenarios based on the EC estuary dataset, scenarios 3 (linear fit)80and 4 (constant mean), are the strongest candidates for defining the model freshwa-ter TA. This dataset is representative of the river mouth which is the approximatephysical location of the freshwater endmember along the estuarine mixing curve(figure 3.8), and the data were sampled recently (i. e. beginning in 2004). ModelnDIC3m, pH3m, and ΩAsur < 1 are not especially sensitive to the differences betweenscenarios 3 and 4 (figure 3.7) except in the case of summer ΩAsur < 1 duration athigh pH f resh (figure 3.7e). The model metrics are more sensitive to the differencebetween scenarios 3 and 4 and scenario 2 (quadratic fit), however scenario 2 is apoor candidate since the freshwater TA values are extrapolated rather than mea-sured and are generally much lower (∼ 200 µeq kg−1, figure 3.6) than those mea-sured in the river (i. e. EC Hope and EC estuary datasets). Since the model metricsare not sensitive to the difference between the constant and flow-dependent TA f reshscenarios based on the EC estuary dataset, I use the constant scenario (scenario 4)as the model TA f resh in chapter 2 for simplicity.The EC buoy dataset provides a satifactory summary of pH near the rivermouth. The buoy is ideally located and the high frequency sampling allows fora temporal resolution in pH otherwise unavailable in the region. Despite the inher-ent inaccuracy of the pH instrument, the lack of more robust pH measurements inthe estuary deem this dataset to be the best alternative. I select a mean buoy pH of7.7 as the model pH f resh for simplicity. The variability of pH f resh at the selectedmodel TA f resh may produce a ∼ 0.1 unit change in pH3m during the summer (fig-ure 3.7c) and a ∼ 5-day change in ΩAsur < 1 duration during summer and winter(figure 3.7e and f).3.5 ConclusionsThe aims of this chapter were to estimate the range and variability of TA and pHin the Fraser River freshwater contribution to the SoG, to assess the sensitivityof model SoG pH and ΩA to this range and variability, and to determine realisticfreshwater pH and TA parameterizations for use in the chapter 2 modeling study.The datasets used in this study demonstrate a total plausible freshwater TArange of∼500 to 1,350 µeq kg−1 near the Fraser River mouth, and an observed pHrange of ∼7.4 to 8.0 at the Environment Canada Fraser River Water Quality Buoy.81Both TA f resh and pH f resh are seasonally-variable with the highest values observedin winter and summer, respectively. Significant correlations between TA f resh ob-servations and parameterized total freshwater discharge in the southern SoG nearthe Fraser River mouth support a flow-dependent mechanism to describe the sea-sonality of TA f resh. While pH f resh may also demonstrate flow dependence nearpeak river flow, observations generally span the range of variability at most flowstages, and pH f resh seasonality is likely determined by primarily local processes(e. g. biological productivity).The range of freshwater TA and pH scenarios produce an average (1 May to Oct31) near-surface model pH range of∼8.1 to 8.3, median winter ΩAsur < 1 durationsof between 140 and 160 days (total range 120 to 185 days), and median summerΩAsur < 1 durations, when present, of between 0 and 25 days (∼50 days maxi-mum for pH f resh = 7.4, TA f resh = 530 µeq kg−1). Increasing TA f resh decreasesmodel pH by increasing model DIC relative to TA. The two-endmember mixingdiscussion supports this mechanism by demonstrating an increase in the estuarineDIC:TA ratio at high TA f resh. However, increasing TA f resh also decreases modelsurface ΩA < 1 duration by raising the ΩA. Despite the high model DIC relative toTA, ΩA is more strongly sensitive to the net rise in estuarine CO2−3 resulting fromthe bulk freshwater DIC increase in this case. At low pH f resh, model pH is moresensitive to changes in TA f resh due to the high freshwater DIC:TA ratio. However,since model surface ΩA < 1 duration sensitivity is determined by a different mech-anism, this quantity is actually stabilized at low pH f resh by simultaneous additionand removal of CO2−3 .I select the mean EC estuary TA f resh scenario (scenario 4, table 3.2), andpH f resh = 7.7 as the freshwater chemistry parameters for the model used in chap-ter 2. The combined sensitivity of model pH and surface ΩA < 1 duration isstrongest at the extreme scenarios and weakest near the mean values I have se-lected. Near-surface model pH and ΩA do not demonstrate sensitivity to differ-ences between the constant and flow-dependent EC estuary TA f resh scenarios (3and 4) at pH f resh = 7.7 so I choose the constant scenario for simplicity.82Chapter 4ConclusionsI set out in this thesis to investigate the role of local physical forcing in determiningthe seasonal variability of near-surface pH and ΩA in the southern Strait of Georgia(SoG). In order to evaluate this connection, I developed a carbonate system modulerepresentative of the system that could be coupled to a collaboratively-developed1-D biophysical model of the southern SoG. The model required tuning and evalu-ation, and a freshwater sensitivity analysis in order to provide realistic estimates ofpH and ΩA at the study site. Using the tuned model with the appropriate freshwaterchemistry parameterizations, I tested the sensitivity of near-surface pH and ΩA to acomprehensive set of local forcing scenarios representative of the period between2001 and 2012.The model generally reproduces observations of pH and ΩA from the top 40 mof the southern SoG, however disagreement arises in the near-surface during sum-mer, primarily in 2010. Where surface values disagree, the model predicts a fresher,warmer surface associated with a strong river influence. Since these conditionslikely represent a nearby station closer to the river plume, the model may be thoughtof as Langrangian rather than representative of only the defined coordinates. Wheresubsurface values disagree, the model underestimates DIC due to overproductivity.Both cases result in high estimates for model pH and ΩA relative to observations.Since only low pH and ΩA are considered to be potentially unfavorable to theecosystem, the model represents a best case scenario.The selection of freshwater endmember concentrations of DIC and TA can af-83fect near-surface pH and ΩA in the model. Between 2001 and 2012 over the rangeof observed freshwater TA and pH in the southern SoG and Fraser River, six-monthaverages (1 May to 31 October) of near-surface model pH vary between 8.1 to 8.3,median winter surface ΩA < 1 durations vary by 20 days, and median summersurface ΩA < 1 durations vary by 25 days. Average model ΩA duration and pHstabilize at opposite ends of the observed freshwater pH range because changes inthe freshwater and estuarine DIC:TA ratios affect each quantity differently. Modelfreshwater pH is thus defined as the mean observed value of 7.7 to optimize sta-bility for both quantities. Since the model does not demonstrate sensitivity to thedifference between constant and flow-depended freshwater TA scenarios for theFraser River mouth at pH 7.7, the mean Fraser River estuary value defines themodel freshwater TA.Near-surface model pH and ΩA follow robust seasonal cycles associated withthe local ecosystem dynamics. Within these cycles, near-surface pH exhibits twoperiods of increased variability (> 0.2 pH units) when evaluated over a wide rangeof realistic forcing scenarios. During each period, pH demonstrates strong nega-tive correlations to two local forcing quantities across these scenarios: wind speedin spring and freshwater flux in summer. Analysis of water-column DIC fluxesreveals two mechanisms of spring pH regulation by wind. In early spring windmixing primarily inhibits phytoplankton growth, weakening or delaying bloomsand keeping pH low. In late spring most blooms have already happened increasingambient near-surface pH significantly, and net wind mixing of DIC to the surfaceacross the newly-established chemocline dominates. Nitrate flux analysis revealsthat phytoplankton are primarily light-limited, as opposed to nutrient-limited, insummer during peak freshwater flux. This finding suggests that increased turbidityis the primary pH regulation mechanism in summer rather than the decreased es-tuarine entrainment velocity (and upward nutrient transport) associated with peakriver discharge. These findings do not rule out increased advection of phytoplank-ton as a mechanism for reduced summer pH, which is likely to contribute to thecorrelation between summer pH and freshwater flux as well.Winter surface aragonite undersaturation (ΩA < 1) is a persistent feature of themodel. Duration ranges from ∼70 to 130 days over the range of forcing scenar-ios, and varies by ∼40 days at onset and 30 days at termination. The relationships84between onset and local forcing quantities are weak, but negative correlations towind speed and cloud fraction point to light limitation, removal of phytoplanktonby mixing, and enhanced DIC transport to the surface as possible drivers of earlywinter ΩA < 1 onset. These correlations improve at low river discharge, suggestingthat other mechanisms associated with higher river discharge weaken the influenceof wind and cloud cover. Strong correlation between termination (spring surfaceΩA > 1 onset) and cloud fraction indicates surface light availability to phytoplank-ton as the dominant factor in terminating winter ΩA < 1. Allen and Wolfe (2013)also found cloud fraction to be important in determining spring diatom bloom tim-ing, however I do not observe a strong relationship between spring ΩA > 1 onsetand spring bloom timing. Perhaps rethinking the spring bloom definition couldstrengthen this relationship in future research.Summer surface ΩA < 1 occurs during years with large (> 5,500 m3 s−1) Julyaverage Fraser River flows. When present, the duration reaches as long as 50 daysover the range of forcing scenarios. With respect to local forcing, duration demon-strates the strongest correlation to the July-average freshwater flux. Comparisonsbetween summer 2010 (small freshet) and summer 2011 (large freshet) model sur-face Ca2+ and CO2−3 concentrations indicate that both ions are deficient by as muchas 3-fold during large freshets. Since ΩA is directly proportional to both ions, fresh-water dilution is the dominant driver of reductions in surface ΩA, although CO2−3reductions associated with biologically-sourced pH decreases in summer may beimportant as well.The model results demonstrate that pH and ΩA extremes in the context of sea-sonal timing are common. For instance, near-surface model pH and ΩA increaseduring spring, but these increases can occur gradually or rapidly at any point be-tween early February and early April. Likewise, low near-surface model pH andΩA are uncharacteristic of summer, but can arise nontheless during strong freshwa-ter fluxes. Low extremes in this context might be observed during years with strongspring winds and cloud cover accompanied by weak, late spring phytoplanktonblooms, or during years with strong Fraser River freshets.It is possible that in other parts of the SoG, the observed range of pH and ΩAmay exceed the range determined by this study for the model site. Such obser-vations would rely on a different balance of fluxes than represented in the model.85High chlorophyll-a observations in the northern SoG relative to the southern basin(Masson and Pen˜a, 2009) are consistent with abnormally high (variable) surfacepH and ΩA observations to the north (D. Ianson, manuscript in preparation, 2014).In protected areas near shorelines, abnormally low pH observations may arise dueto benthic processes (e. g. enhanced remineralization) associated with DIC pro-duction (Johannessen et al., 2008)The majority of shellfish impact studies report the greatest vulnerability to beduring the larval stage (Hettinger et al., 2012; Kroeker et al., 2013; Parker et al.,2013). By contrast, wild adult shellfish in the SoG cannot be especially vulnerableto the observed low extremes of pH and ΩA because they experience the full rangeof observed conditions during their lifespan without persistent population decline.The spawning season for most commercial clam varieties in the SoG occurs be-tween May and August, during which time each animal will spend 2 to 3 weeks inthe larval stage before settling (Quayle, 1970). In most years, near-surface modelpH exceeds 8.0 units and ΩA exceeds saturation during this season. It is unlikelythat wind-driven low pH and cloud-driven low ΩA would persist into the spawningseason since both scenarios appear to weaken by the end of April. However, lownear-surface summer pH and ΩA driven by large freshwater fluxes could create un-favorable growth conditions for new shellfish larvae during some years. This effectwill be much less severe at large distances from the Fraser River plume, so shell-fish in the northern SoG are unlikely to be affected by corrosive, brackish water.Additionally, aragonite undersaturation due to large freshwater fluxes is associatedwith severe freshening events, and most shellfish living in these zones may alreadybe adapted to fresher, undersaturated environments (e. g. (Parker et al., 2012)).Aquaculture shellfish (oysters, clams, scallops) may experience more severeimpacts from local carbonate chemistry than resident shellfish due to variable seedtiming and depth practices used by the industry. Seeding events taking place beforethe local spawning season may risk exposure to unfavorable low pH and ΩA condi-tions. Additionally, seeding stocks too deep may unintentionally expose vulnerablelarvae to shallow, undersaturated waters. While some seed is stored directly in theSoG prior to seeding (D. Ianson, personal communication, 2014), seed temporar-ily stored in land facilities are also vulnerable to local carbonate chemistry as thesefacilities use SoG water for their storage containers. Many facilities use liming86practices for their incoming water to mitigate low pH and ΩA, however there arestill many other aspects of these land facilities that do not reflect the natural envi-ronment and may reduce fitness.Quantifying the long term trends of near-surface SoG pH and ΩA due to an-thropogenic CO2 forcing and climate change is beyond the scope of this study.However such trends are plausible given the model setup and sensitivity results.Direct anthropogenic CO2 forcing sourced in inflowing PIW (Feely et al., 2008,2012a) would be parameterized in the model as a change in the 40 m boundarycondition. Upwelling source waters were DIC-rich relative to Puget Sound basinwaters during a summer 2008 cruise in Juan de Fuca Strait and Puget Sound (Feelyet al., 2012a). However, upwelling source waters are consistently DIC-poor rela-tive to deep SoG waters on similar density surfaces based on recent observationsin Juan de Fuca Strait and the SoG (D. Ianson, manuscript in preparation, 2014),and so the relative strength of oceanic DIC is a driver of SoG acidification remainsunclear.Climate forcing could potentially drive changes in the local forcing climatol-ogy, altering seasonal wind, cloud, and river flow behavior. Hindcasts from 1968 to2010 using a diatom-only version of the 1-D model in this study demonstrate thatspring bloom timing became more variable after 1990, and progressively later after2005 (Allen and Wolfe, 2013). The authors attributed this variability to changesin the timing of the spring transition of prevailing wind patterns over the NorthPacific. Similar large scale to local scale effects could produce more highly vari-able and potentially longer periods of low near-surface pH and aragonite under-saturation in the model used here. Fraser River flow patterns are also forecastedto change with shifts in snow patterns in the Fraser River watershed, with futureflow profiles rising earlier in the season and peaking at smaller freshets (Morrisonet al., 2002). Such flow would produce overall higher near-surface pH and ΩA inthe model by increasing stratification in spring without adding light limitation insummer. However Collins et al. (2009) suggested that increased freshwater flux inspring increases phytoplankton advection, reducing bloom strength.Long term changes in near-surface pH and ΩA seasonality could disrupt thecurrent SoG ecosystem. Given the predominantly negative effects of low pH andΩA on secondary producers and higher trophic levels determined from recent meta-87analyses (Kroeker et al., 2013; Parker et al., 2013) and the stronger sensitivity tolow pH and ΩA observed in a variety of organisms at early growth stages (e. g.oysters; Hettinger et al., 2012; Barton et al., 2012), changes in seasonal carbonatechemistry could expose the vulnerable spawning stage of local organisms to unfa-vorable conditions. Resident shellfish populations could be particularly affected,with the decline in benthic structure negatively impacting the benthic community(Gutie´rrez et al., 2003).The next step in this research is to investigate the potential long-term changesto SoG carbonate chemistry. This study provides a starting point and the modelprovides a starting tool. I have described the strengths and limitations of thismodel and their implications in its use for long term projections. With furtherdevelopment, this model is suited to be used in a long-term study. However, asbiophysical-coupled modeling continues to advance in the region, the role of 3-Dmodels in addressing these long-term trends will be important as well.88BibliographyS. R. Alin, R. A. Feely, A. G. Dickson, J. M. Herna´ndez-Ayo´n, L. W. Juranek,M. D. Ohman, and R. Goericke. Robust empirical relationships for estimatingthe carbonate system in the southern California Current System and applicationto CalCOFI hydrographic cruise data (20052011). J. Geophys. Res., 117:C05033, 2012.S. E. Allen and M. A. Wolfe. Hindcast of the timing of the spring phytoplanktonbloom in the Strait of Georgia, 19682010. Prog. Oceanogr., 26:81–87, 2013.Y. Artioli, J. C. Blackford, M. Butenscho¨n, J. T. Holt, S. L. Wakelin, H. Thomas,A. V. Borges, and J. I. Allen. The carbonate system in the North Sea:Sensitivity and model validation. J. Mar. Sys., 102-104:1–13, 2012.U. Bamstedt and K. Karlson. Euphausiid predation on copepods in coastal waterof the Northeast Atlantic. Mar. Ecol. Prog. Ser., 172:149–168, 1998.A. Barton, B. Hales, G. G. Waldbusser, C. Langdon, and R. A. Feely. The Pacificoyster, Crassostrea gigas, shows negative correlation to naturally elevatedcarbon dioxide levels: implications for near-term ocean acidification effects.Limnol. Oceanogr., 57(3):698–710, 2012.N. R. Bates. Interannual variability of the oceanic CO2 sink in the subtropicalgyre of the North Atlantic Ocean over the last 2 decades. J. Geophys. Res., 112:C09013, 2007.N. Bednarsek, G. A. Tarling, D. C. E. Bakker, S. Fielding, E. M. Jones, H. J.Venables, P. Ward, A. Kuzirian, B. Le´ze´, R. A. Feely, and E. J. Murphy.Extensive dissolution of live pteropods in the Southern Ocean. Nature Geosci.,5(12):881–885, 2012.M. J. Behrenfeld. Abandoning Sverdrup’s Critical Depth Hypothesis onphytoplankton blooms. Ecol., 91(4):977–989, 2010.89L. Bianucci, K. L. Denman, and D. Ianson. Low oxygen and high inorganiccarbon on the Vancouver Island Shelf. J. Geophys. Res., 116:C07011, 2011.B. Bolin and H. Rodhe. A note on the concepts of age distribution and transit timein natural reservoirs. Tellus, 25(1):58–62, 1973.K. E. Brainerd and M. C. Gregg. Surface mixed and mixing layer depths. DeepSea Res. I, 42(9):1521–1543, 1995.P. G. Brewer. Direct observation of the oceanic CO2 increase. Geophys. Res.Lett., 5(12):997–1000, 1978.P. G. Brewer and J. C. Goldman. Alkalinity changes generated by phytoplanktongrowth. Limnol. Oceanogr., 21(1):108–117, 1976.R. H. Byrne, S. Mecking, R. A. Feely, and X. Liu. Direct observations ofbasin-wide acidification of the North Pacific Ocean. Geophy. Res. Lett., 37:L02601, 2010.W.-J. Cai, Y. Wang, and R. E. Hodson. Acid-base properties of dissolved organicmatter in the estuarine waters of Georgia, USA. Geochim. Cosmochim. Ac., 62(3):473–483, 1998.K. Caldeira and M. E. Wickett. Anthropogenic carbon and ocean pH. Nature,425:365–365, 2003.K. Caldeira and M. E. Wickett. Ocean model predictions of chemistry changesfrom carbon dioxide emissions to the atmosphere and ocean. J. Geophys. Res.,110:C09S04, 2005.E. M. Cameron. Hydrogeochemistry of the Fraser River, British Columbia:seasonal variation in major and minor components. J. Hydrol., 182:209–225,1996.E. M. Cameron and K. Hattori. Strontium and neodymium isotope ratios in theFraser River, British Columbia: a riverine transect across the Cordilleranorogen. Chem. Geol., 137:209–225, 1997.E. M. Cameron, G. E. M. Hall, J. Veizer, and H. R. Krouse. Isotopic andelemental hydrogeochemistry of a major river system: Fraser River, BritishColumbia, Canada. Chem. Geol., 122:149–169, 1995.R. W. Campbell and J. F. Dower. Depth distribution during the life history ofNeocalanus plumchrus in the Strait of Georgia. J. Plankton Res., 30:7–20,2008.90A. K. Collins, S. E. Allen, and R. Pawlowicz. The role of wind in determining thetiming of the spring bloom in the Strait of Georgia. Can. J. Fish. Aquat. Sci.,66:1597–1616, 2009.A. K. Covington, P. D. Whalley, and W. Davison. Procedures for themeasurement of pH in low ionic strength solutions including freshwater.Analyst, 108:1528–1532, 1983.P. B. Crean, T. S. Murty, and J. A. Stronarch. Mathematical modeling of tides andestuarine circulation: the coastal seas of southern British Columbia andWashington State. In Lecture notes on coastal and estuarine studies,volume 30, page 471. Springer-Verlag, 1988.P. Cugier, A. Me´nesguen, and J.-F. Guillaud. Three-dimensional (3d) ecologicalmodelling of the Bay of Seine (English Channel, France). J. Sea. Res., 54:104–124, 2005.S. J. de Mora. Manganese chemistry in the Fraser Estuary. PhD thesis,University of British Columbia, 1981.S. J. de Mora. The distribution of alkalinity and pH in the Fraser Estuary.Environ. Technol. Lett., 4:35–46, 1983.K. Denman. Modelling planktonic ecosystems: parameterizing complexity. Prog.Oceanogr., 57:429–452, 2003.K. Denman and A. Pen˜a. The response of two coupled one-dimensional mixedlayer/planktonic ecosystem models to climate change in the NE subarcticPacific Ocean. Deep Sea Res. II, 49:5739–5757, 2002.A. G. Dickson. An exact definition of total alkalinity and a procedure for theestimation of alkalinity and total inorganic carbon from titration data. Deep SeaRes., 28(6):609–623, 1981.A. G. Dickson. The measurement of sea water pH. Mar. Chem., 44:131–142,1993.A. G. Dickson and C. Goyet. Handbook of methods for the analysis of the variousparameters of the carbon dioxide system in sea water. DOE (U.S. Departmentof Energy), 2nd edition, 1994. ORNL/CDIAC-74.A. G. Dickson, C. L. Sabine, and J. R. Christian. Guide to Best Practices forOcean CO2 Measurements. PICES Special Publication 3, 2007.91S. C. Doney, V. J. Fabry, R. A. Feely, and J. A. Kleypas. Ocean acidication: theother CO2 problem. Annu. Rev. Mar. Sci., 1:169–192, 2009.J. E. Dore, R. Lukas, D. W. Sadler, M. J. Church, and D. M. Karl. Physical andbiogeochemical modulation of ocean acidification in the central North Pacific.PNAS, 106(30):12235–12240, 2009.R. El-Sabaawi, J. F. Dower, M. Kainz, and A. Mazumder. Interannual variabilityin fatty acid composition of the copepod Neocalanus plumchrus in the Strait ofGeorgia, British Columbia. Mar. Ecol. Prog. Ser., 382:151–161, 2009.Environment Canada. Summary of standard operating procedure for: carbon -TIC/TOC combustion-infrared, acidification, sparging. Environment CanadaChemistry Section, Pacific Environmental Science Centre, 2005. v2.4.Environment Canada. Summary of standard operating procedure for: alkalinityby electrometric titration. Environment Canada Chemistry Section, PacificEnvironmental Science Centre, 2006. v3.4.R. W. Eppley. Temperature and phytoplankton growth in the sea. Fish. Bull. Nat.Ocean. Atmos. Adm., 70(4):1063–1085, 1972.A. Ethier and J. Bedard. Development of a real-time water quality buoy for theFraser River Estuary. AXYS Technologies Inc., Sidney, British Columbia,2007.W. Evans, B. Hales, P. G. Strutton, and D. Ianson. Sea-air CO2 fluxes in thewestern Canadian coastal ocean. Prog. Oceanogr., 101(1):78–91, 2012.M. J. R. Fasham. Variations in the seasonal cycle of biological production insubarctic oceans: a model sensitivity analysis. Deep-Sea Res. Pt. I, 42(7):1111–1149, 1995.R. A. Feely, C. L. Sabine, K. Lee, W. Berelson, J. Kleypas, V. J. Fabry, and F. J.Millero. Impact of anthropogenic CO2 on the CaCO3 system in the oceans.Science, 305:362–366, 2004.R. A. Feely, C. L. Sabine, J. M. Herna´ndez-Ayo´n, D. Ianson, and B. Hales.Evidence for upwelling of corrosive ”acidified” water onto the continentalshelf. Science, 320:1490–1492, 2008.R. A. Feely, S. R. Alin, J. Newton, C. L. Sabine, M. Warner, A. Devol,C. Krembs, and C. Maloy. The combined effects of ocean acidification, mixing,and respiration on pH and carbonate saturation in an urbanized estuary. Estuar.Coast. Shelf. Sci, 88:442–449, 2010.92R. A. Feely, T. Klinger, J. A. Newton, and M. Chadsey. Scientific Summary ofOcean Acidification in Washington State Marine Waters. NOAA OAR SpecialReport, 2012a.R. A. Feely, C. L. Sabine, R. H. Byrne, F. J. Millero, A. G. Dickson,R. Wanninkhof, A. Murata, L. A. Miller, and D. Greely. Decadal changes in thearagonite and calcite saturation state of the Pacific Ocean. Glob. Biogeochem.Cycles, 26(3):2011GB004157, 2012b.M. G. G. Foreman, R. A. Walters, R. F. Henry, C. P. Keller, and A. G. Dolling. Atidal model for eastern Juan de Fuca Strait and the southern Strait of Georgia. J.Geophys. Res., 100(C1):721–740, 1995.H. J. Freeland and K. L. Denman. A topographically controlled upwelling centreoff southern Vancouver Island. J. Mar. Res., 40(4):1069–1093, 1982.W. R. Geyer and G. A. Cannon. Sill processes related to deep water renewal in afjord. J. Geophys. Res. Oceans, 87(C10):7985–7996, 1982.S. N. Giddings, P. MacCready, B. M. Hickey, N. S. Banas, K. A. Davis, S. A.Siedlecki, V. L. Trainer, R. M. Kudela, N. A. Pelland, and T. P. Connolly.Hindcasts of potential harmful algal bloom transport pathways on the PacificNorthwest coast. J. Geophys. Res. Oceans, 119:2439–2461, 2014.M. Gonza´lez-Da´vila, J. M. Santana-Casiano, and E. F. Gonza´lez-Da´vila.Interannual variability of the upper ocean carbon cycle in the northeast AtlanticOcean. Geophys. Res. Lett., 34:L07608, 2007.J. Gower, S. King, S. Statham, R. Fox, and E. Young. The Malaspina Dragon: anewly-discovered pattern of the early spring bloom in the Strait of Georgia,British Columbia, Canada. Prog. Oceanogr., 115:181–188, 2013.N. Gruber, C. Hauri, Z. Lachkar, D. Loher, T. L. Fro¨licher, and G.-K. Plattner.Rapid progression of ocean acidification in the California Current System.Science, 337:220–223, 2012.J. L. Gutie´rrez, C. G. Jones, D. L. Strayer, and O. O. Iribarne. Mollusks asecosystem engineers: the role of shell production in aquatic habitats. Oikos,101:79–90, 2003.N. Gypens, G. Lacroix, C. Lancelot, and A. V. Borges. Seasonal and inter-annualvariability of airsea CO2 fluxes and seawater carbonate chemistry in theSouthern North Sea. Prog. Oceanogr., 88:59–77, 2011.93M. J. Halverson and S. Fleming. Complex networks, streamflow, and hydrometricmonitoring system design. Hydrol. Earth Sys. Sci. Disc., 11:13711–13744,2014.M. J. Halverson and R. Pawlowicz. Estuarine forcing of a river plume by riverflow and tides. J. Geophys. Res., 113:C09033, 2008.M. J. Halverson and R. Pawlowicz. High-resolution observations of chlorophyll-abiomass from an instrumented ferry: influence of the Fraser River plume from2003 to 2006. Cont. Shelf Res., 59:52–64, 2013.P. J. Harrison, J. D. Fulton, F. J. R. Taylor, and T. R. Parsons. Review of thebiological oceanography of the Strait of Georgia: pelagic environment. Fish.Aquat. Sci., 40:1064–1094, 1983.P. J. Harrison, P. J. Clifford, W. P. Cochlan, K. Yin, M. A. St. John, P. A.Thompson, M. J. Sibbald, and L. J. Albright. Nutrient and plankton dynamicsin the Fraser River plume, Strait of Georgia, British Columbia. Mar. Ecol.Prog. Ser., 70:291–304, 1991.C. Hauri, N. Gruber, M. Vogt, S. C. Doney, R. A. Feely, Z. Lachkar,A. Leinweber, A. M. P. McDonnell, M. Mu¨nnich, and G.-K. Plattner.Spatiotemporal variability and long-term trends of ocean acidification in theCalifornia Current System. Biogeosci., 10:193–216, 2013.J. M. Herna´ndez-Ayo´n, A. Zirino, A. G. Dickson, T. Camiro-Vargas, andE. Valenzuela-Espinoza. Estimating the contribution of organic bases frommicroalgae to the titration alkalinity in coastal seawaters. Limnol. Oceanogr.Methods, 5:225–232, 2007.A. Hettinger, E. Sanford, T. M. Hill, A. D. Russell, and K. N. S. Sato. Persistentcarry-over effects of planktonic exposure to ocean acidification in the Olympiaoyster. Ecology, 93:2758–2768, 2012.B. M. Hickey, R. E. Thomson, H. Yih, and P. H. LeBlond. Velocity andtemperature fluctuations in a buoyancy-driven current off Vancouver Island. J.Geophys. Res., 96(C6):10507–10538, 1991.A. F. Hofmann, J. J. Middleburg, K. Soetaert, and F. J. R. Meysman. pHmodelling in aquatic systems with time-variable acid-base dissociationconstants applied to the turbid, tidal Scheldt estuary. Biogeosci., 6:1539–1561,2009.94C. W. Hunt, J. E. Salisbury, and D. Vandemark. Contribution of non-carbonateanions to total alkalinity and overestimation of pCO2 in New England and NewBrunswick rivers. Biogeosci., 8:3069–3076, 2011.D. Ianson. The increase in carbon along the Canadian Pacific coast. In J. R.Christian and M. G. G. Foreman, editors, Climate Trends and Projections forthe Pacific Large Aquatic Basin, pages 57–66. Can. Tech. Rep. Fish. Aquat.Sci. 3032, 2013.D. Ianson and S. E. Allen. A two-dimensional nitrogen and carbon flux model ina coastal upwelling region. Global Biogeochem. Cycles, 16(1):11–1–11–16,2002.D. Ianson, S. E. Allen, S. L. Harris, K. J. Orians, D. E. Varela, and C. S. Wong.The inorganic carbon system in the coastal upwelling region west of VancouverIsland, Canada. Deep Sea Res. I, 50:1023–1042, 2003.IOC, SCOR, and IAPSO. The international thermodynamic equation of seawater2010: Calculation and use of thermodynamic properties. IntergovernmentalOceanographic Commission, Manuals and Guides No. 56, UNESCO (English),2010.N. Jeffery. Modelling a phytoplankton dichotomy in the eastern subarctic Pacific:Impact of atmospheric variability, iron surface flux, and life cycle dynamics ofthe calanoid copepods, Neocalanus spp. PhD thesis, University of BritishColumbia, 2002.S. C. Johannessen, G. Potentier, C. A. Wright, D. Masson, and R. W. Macdonald.Water column organic carbon in a Pacific marginal sea (Strait of Georgia,Canada). Mar. Environ. Res., 66:S49–S61, 2008.L. W. Juranek, R. A. Feely, W. T. Peterson, S. R. Alin, B. Hales, K. Lee, C. L.Sabine, and J. Peterson. A novel method for determination of aragonitesaturation state on the continental shelf of central oregon using multi-parameterrelationships with hydrographic data. Geophys. Res. Lett., 36:L24601, 2009.M. Kawamiya, M. J. Kishi, Y. Yamanaka, and N. Suginohara. Anecological-physical coupled model applied to station papa. J. Oceanogr., 51:635–664, 1995.V. C. Kennedy. Minerology and cation-exchange capacity of sediments fromselected streams. Prof. Paper 433-D, U.S. Geol. Surv., 1965.95T. Khangaonkar, B. Sackmann, W. Long, T. Mohamedali, and M. Roberts.Simulation of annual biogeochemical cycles of nutrient balance, phytoplanktonbloom(s), and DO in Puget Sound using an unstructured grid model. OceanDynamics, 62:1353–1379, 2012.H.-C. Kim and K. Lee. Significant contribution of dissolved organic matter toseawater alkalinity. Geophys. Res. Lett., 36:L20603, 2009.W. Koeve and A. Oschlies. Potential impact of DOM accumulation on f CO2 andcarbonate ion computations in ocean acidification experiments. Biogeosci., 9:3787–3798, 2012.K. J. Kroeker, R. L. Kordas, R. N. Crim, and G. G. Singh. Meta-analysis revealsnegative yet variable effects of ocean acidification on marine organisms. Ecol.Lett., 13:1419–1434, 2010.K. J. Kroeker, R. L. Kordas, R. N. Crim, I. E. Hendriks, L. Ramajo, G. G. Singh,C. M. Duarte, and J.-P. Gattuso. Impacts of ocean acidification on marineorganisms: quantifying sensitivities and interaction with warming. Glob.Change Biol., 19:1884–1896, 2013.I. Kuznetsov and T. Neumann. Simulation of carbon dynamics in the Baltic Seawith a 3D model. J. Mar. Sys., 111-112:167–174, 2013.A. Lara Espinosa. Determination of the acidification state of Canadian Pacificcoastal waters using empirical relationships with hydrographic data. Master’sthesis, University of Victoria, 2012.W. G. Large. Modeling and parametrizing the ocean planetary boundary layer. InE. P. Chassignet and J. Verron, editors, Ocean Modeling and Parameterization,pages 81–120. NATO Science Series, 1998.W. G. Large and S. Pond. Sensible and latent heat flux measurements over theocean. J. Phys. Oceanogr., 12:464–482, 1982.W. G. Large, J. C. McWilliams, and S. C. Doney. Oceanic vertical mixing: areview and a model with a nonlocal boundary layer parametrization. Rev.Geophys., 32(4):363–403, 1994.P. H. LeBlond. The Strait of Georgia: functional anatomy of a coastal sea. Can. J.Fish. Aquat. Sci., 40:1033–1063, 1983.K. Lee, L. T. Tong, F. J. Millero, C. L. Sabine, A. G. Dickson, C. Goyet, G.-H.Park, R. Wanninkhof, R. A. Feely, and R. M. Key. Global relationships of total96alkalinity with salinity and temperature in surface waters of the worlds oceans.Geophys. Res. Lett., 33:L19605, 2006.E. R. Lewis and D. W. R. Wallace. Program Developed for CO2 SystemCalculations. DOE (U.S. Department of Energy), 1998. ORNL/CDIAC-105.S. Lischka and U. Riebesell. Synergistic effects of ocean acidification andwarming on overwintering pteropods in the Arctic. Glob. Change Biol., 18:3517–3528, 2012.T. J. Lueker, A. G. Dickson, and C. D. Keeling. Ocean pCO2 calculated fromdissolved inorganic carbon, alkalinity, and equations for K1 and K2: validationbased on laboratory measurements of CO2 in gas and seawater at equilibrium.Mar. Chem., 70:105–119, 2000.D. L. Mackas and P. J. Harrison. Nitrogenous nutrient sources and sinks in theJuan de Fuca Strait/Strait of Georgia/Puget Sound estuarine system: assessingthe potential for eutrophication. Estuarine Coastal Shelf Sci., 44:1–21, 1997.D. L. Mackas, M. D. Galbraith, D. Faust, D. Masson, K. Young, W. Shaw,S. Romaine, M. Trudel, J. F. Dower, R. W. Campbell, A. R. Sastri, E. A. B.Pechter, E. Pakhomov, and R. W. El-Sabaawi. Zooplankton time series fromthe Strait of Georgia: results from year-round sampling at deep water locations,1990-2010. Prog. Oceanogr., 115:129–159, 2013.D. Masson. Deep water renewal in the Strait of Georgia. Estuar. Coast. Shelf.Sci., 54:115–126, 2002.D. Masson. Seasonal water mass analysis for the Straits of Juan de Fuca andGeorgia. Atmos. Ocean, 44(1):1–15, 2006.D. Masson and P. F. Cummins. Observations and modeling of seasonal variabilityin the Straits of Georgia and Juan de Fuca. J. Mar. Res., 62:491–516, 2004.D. Masson and A. Pen˜a. Chlorophyll distribution in a temperate estuary: TheStrait of Georgia and Juan de Fuca Strait. Estuarine Coastal Shelf Sci., 82:19–28, 2009.G. L. Mellor and T. Yamada. Development of a turbulence closure model forgeophysical fluid problems. Rev. Geophys. Space Phys., 20(4):851–875, 1982.F. J. Millero. Thermodynamics of the carbonate sytem in seawater. Geochim.Cosmochim. Ac., 43(10):1651–1661, 1979.97F. J. Millero. Carbonate constants for estuarine waters. Mar. Fresh. Res., 61:139–142, 2010.W. G. Mook and B. K. S. Koene. Chemistry of dissolved inorganic carbon inestuarine and coastal brackish waters. Estuar. Coast. Mar. Sci., 3:325–336,1975.J. Morrison, M. C. Quick, and M. G. G. Foreman. Climate change in the FraserRiver watershed: flow and temperature projections. J. Hydrol., 263:230–244,2002.J. Morrison, M. G. G. Foreman, and D. Masson. A method for estimating themonthly freshwater discharge affecting British Columbia coastal waters.Atmos. Ocean, 50(1):1–8, 2012.J. W. Morse and F. T. Mackenzie. Geochemistry of sedimentary carbonates. InDevelopments in Sedimentology, volume 48, page 707. Elsevier, 1990.A. Mucci. The solubility of calcite and aragonite in seawater at various salinities,temperatures, and one atmosphere total pressure. Am. J. Sci., 283(7):780–799,1983.P. D. Nightingale, G. Malin, C. S. Law, A. J. Watson, P. S. Liss, M. I. Liddicoat,J. Boutin, and R. C. Upstill-Goddard. In situ evaluation of air-sea gas exchangeparametrizations using novel conservative and volatile tracers. GlobalBiogeochem. Cycles, 14(1):373–387, 2000.R. V. O’Neill, D. L. DeAngelis, J. J. Pastor, B. J. Jackson, and W. M. Post.Multiple nutrient limitations in ecological models. Ecol. Modell., 46:147–163,1989.J. C. Orr, V. J. Fabry, O. Aumont, L. Bopp, S. C. Doney, R. A. Feely,A. Gnanadesikan, N. Gruber, A. Ishida, F. Joos, R. M. Key, K. Lindsay,E. Maier-Reimer, R. Matear, P. Monfray, A. Mouchet, R. G. Najjar, G.-K.Plattner, K. B. Rodgers, C. L. Sabine, J. L. Sarmiento, R. Schlitzer, R. D.Slater, I. J. Totterdell, M.-F. Weirig, Y. Yamanaka, and A. Yool. Anthropogenicocean acidification over the twenty-first century and its impact on calcifyingorganisms. Nature, 437:681–686, 2005.L. M. Parker, P. M. Ross, W. A. O’Connor, L. Borysko, D. A. Raftos, and H. O.Po¨rtner. Adult exposure influences offspring response to ocean acidification inoysters. Glob. Change Biol., 18(1):82–92, 2012.98L. M. Parker, P. M. Ross, W. A. O’Connor, H. O. Po¨rtner, E. Scanes, and J. M.Wright. Predicting the response of molluscs to the impact of oceanacidification. Biology, 2:651–692, 2013.R. Pawlowicz, O. Riche, and M. Halverson. The circulation and residence time ofthe Strait of Georgia using a simple mixing-box approach. Atmos. Ocean, 45(4):173–193, 2007.A. Pen˜a, K. Denman, S. E. Calvert, R. E. Thomson, and J. R. Forbes. The seasonalcycle in sinking particle fluxes off Vancouver Island, British Columbia. DeepSea Res. II: Topical Studies in Oceanography, 46(11-12):2969–2992, 1999.H. Preston-Thomas. The international temperature scale of 1990 (ITS-90).Metrologia, 27(107):1990, 1990.J. F. Price, R. A. Weller, and R. Pinkel. Diurnal cycling: Observations and modelsof the upper ocean response to diurnal heating, cooling, and wind mixing. J.Geophys. Res. Oceans, 91(C7):8411–8427, 1986.D. B. Quayle. The Intertidal Bivalves of British Columbia. British ColumbiaProvincial Museum: Handbook, 1970.A. C. Redfield, B. H. Ketchum, and F. A. Richards. The influence of organisms onthe composition of seawater. In M. N. Hill, editor, The composition ofsea-water: comparative and descriptive oceanography, pages 26–77. JohnWiley and Sons, New York, 1963.J. C. P. Reum, S. R. Alin, R. A. Feely, J. Newton, M. Warner, and P. McElhany.Seasonal carbonate chemistry covariation with temperature, oxygen, andsalinity in a fjord estuary: implications for the design of ocean acidificationexperiments. PLoS ONE, 9(2):e89619, 2014.O. Riche. Time-dependent inverse box-model for the estuarine circulation andprimary productivity in the Strait of Georgia. PhD thesis, University of BritishColumbia, 2011.U. Riebesell and P. Tortell. Effects of ocean acidification on pelagic organismsand ecosystems. In Ocean Acidification, pages 99–121. Oxford UniversityPress, 2011.U. Riebesell, I. Zondervan, B. Rost, P. D. Tortell, R. E. Zeebe, and F. M. M.Morel. Reduced calcification of marine plankton in response to increasedatmospheric CO2. Nature, 407:364–367, 2000.99J. B. Ries, A. L. Cohen, and D. C. McCorkle. Marine calcifiers exhibit mixedresponses to CO2-induced ocean acidification. Geol., 37(12):1131–1134, 2009.J. P. Riley and M. Tongudai. The major cation/chlorinity ratios in sea water.Chem. Geol., 2:263–269, 1967.L. M. Roger, A. J. Richardson, A. D. McKinnon, B. Knott, R. Matear, andC. Scadding. Comparison of the shell structure of two tropical Thecosomata(Creseis acicula and Diacavolinia longirostris) from 1963 to 2009: potentialimplications of declining aragonite saturation. ICES J. Mar. Sci., 69(3):465–474, 2012.L. Royer and W. J. Emery. Variations of the Fraser River plume and theirrelationship to forcing by tide, wind and discharge. Atmos. Ocean, 20(4):357–372, 1982.J. L. Sarmiento and N. Gruber. Ocean Biogeochemical Dynamics. PrincetonUniversity Press, 2006.E. L. Snauffer, D. Masson, and S. E. Allen. Modelling the dispersal of herring andhake larvae in the Strait of Georgia for the period 20072009. Fish. Oceanogr.,23(4):375–388, 2014.M. A. St. John, S. C. Marinone, J. A. Stronach, P. J. Harrison, J. Fyfe, and R. J.Beamish. A horizontally resolving physical-biological model of nitrateconcentration and primary productivity in the Strait of Georgia. Can. J. Fish.Aquat. Sci., 50:1456–1466, 1993.J. H. Steele. Environmental control of photosynthesis in the sea. Limnol.Oceanogr., 7:137–150, 1962.J. G. Stockner, D. D. Cliff, and D. B. Buchanan. Phytoplankton production anddistribution in Howe Sound, British Columbia: a coastal marineembayment-fjord under stress. J. Fish. Res. Board Can., 34(7):907–917, 1977.J. G. Stockner, D. D. Cliff, and K. R. S. Shortreed. Phytoplankton ecology of theStrait of Georgia, British Columbia. J. Fish. Res. Board Can., 36(6):657–666,1979.S. L. Strom, E. L. Macri, and M. B. Olson. Microzooplankton grazing in thecoastal Gulf of Alaska: variations in top-down control of phytoplankton.Limnol. Oceanogr., 52(4):1480–1494, 2007.100V. Stuart and S. C. Pillar. Diel grazing patterns of all ontogenetic stages ofEuphausia lucens and in situ predation rates on copepods in the southernBenguela upwelling region. Mar. Ecol. Prog. Ser., 64:227–241, 1990.D. A. Sutherland, P. MacCready, N. S. Banas, and L. F. Smedstad. A model studyof the Salish Sea estuarine circulation. J. Phys. Oceanogr., 41:1125–1143,2011.H. U. Sverdrup. On conditions for the vernal blooming of phytoplankton. J. Cons.Int. Explor. Mer., 18:287–295, 1953.L. G. Swain. Water quality assessment of Fraser River at Hope (1979 - 2004).Technical report, Environment Canada, 2007.R. E. Thomson. Oceanography of the British Columbia coast. Can. Spec. Publ.Fish. Aquat. Sci. 56, 1981.UNESCO. The practical salinity scale 1978 and the international equation of stateof seawater 1980. Technical report, UNESCO technical papers in marinescience No. 36, 1981.UNESCO. Methodology for oceanic CO2 measurements. Technical report,UNESCO technical papers in marine science No. 65, 1988.B. M. Voss, B. Peucker-Ehrenbrink, T. I. Eglinton, G. Fiske, Z. A. Wang, K. A.Hoering, D. B. Montluc¸on, C. LeCroy, S. Pal, S. Marsh, S. L. Gillies,A. Janmaat, M. Bennett, B. Downey, J. Fanslau, H. Fraser,G. Macklam-Harron, M. Martinec, and B. Wiebe. Tracing river chemistry inspace and time: dissolved inorganic constituents of the Fraser River, Canada.Geochim. Cosmochim. Ac., 124:283–308, 2014.M. Waldichuk. Physical oceanography of the Strait of Georgia, British Columbia.J. Fish. Res. Board Can., 14(3):321–486, 1957.M. Wallin, I. Buffam, M. O¨quist, H. Laudon, and K. Bishop. Temporal and spatialvariability of dissolved inorganic carbon in a boreal stream network:Concentrations and downstream fluxes. J. Geophys. Res., 115(G2):G02014,2010.R. Wanninkhof. Relationship between wind speed and gas exchange over theocean. J. Geophys. Res., 97:7373–7382, 1992.J. C. Warner, C. R. Sherwood, H. G. Arango, and R. P. Signell. Performance offour turbulence closure models implemented using a generic length scalemethod. Ocean Modelling, 8:81–113, 2005.101R. F. Weiss. Carbon dioxide in water and seawater: the solubility of a non-idealgas. Mar. Chem., 2:203–215, 1974.J. O. Wheeler, A. J. Brookfield, H. Gabrielse, J. W. H. Monger, H. W. Tipper, andG. J. Woodsworth. Terrane map of the Canadian Cordillera. Map 1713A, Geol.Surv. Can., 1991. scale 1:200,000.M. Whitfield and D. R. Turner. The carbon dioxide system in estuaries - aninorganic perspective. Sci. Total Environ., 49:235–255, 1986.P. H. Whitfield. Regionalization of water quality in the upper Fraser River basin,British Columbia. Limnol. Oceanogr., 17(9):1053–1066, 1983.P. H. Whitfield and H. Schreier. Hysteresis in relationships between discharge andwater chemistry in the Fraser River basin, British Columbia. Limnol.Oceanogr., 26(6):1179–1182, 1981.K. Yin, R. H. Goldblatt, P. J. Harrison, M. A. St. John, P. J. Clifford, and R. J.Beamish. Importance of wind and river discharge in influencing nutrientdynamics and phytoplankton production in summer in the central Strait ofGeorgia. Mar. Ecol. Prog. Ser., 161:173–183, 1997a.K. Yin, P. J. Harrison, R. H. Goldblatt, M. A. St. John, and R. J. Beamish. Factorscontrolling the timing of the spring bloom in the Strait of Georgia estuary,British Columbia, Canada. Mar. Ecol. Prog. Ser., 161:173–183, 1997b.R. E. Zeebe and D. Wolf-Gladrow. CO2 in Seawater: Equilibrium, Kinetics,Isotopes. Elsevier Oceanography Series, 65, 2001.102Appendix AModel EquationsNutrient uptakeUnder light limitation, the nutrient uptake rate for phytoplankton class i is given byuic = Rimax(T )[1− exp(−IPAR0.33Iiopt)][1.06exp(−IPAR30Iiopt)](A.1)where Rimax(T ) is the maximum growth rate for phytoplankton class i, IPAR is theavailable phytosynthetically active radiation, and Iiopt is the optimum light intensityfor class i.The N-limited uptake rate for phytoplankton class i is determined using thesubstitutable model of O’Neill et al. (1989, in Jeffery, 2002)uino = Rimax(T )(κi[NO−3 ]Kinh +κi[NO−3 ]+ [NH+4 ])(A.2a)uinh = Rimax(T )([NH+4 ]Kinh +κi[NO−3 ]+ [NH+4 ])(A.2b)where Kinh is the half-saturation constant for NH+4 for phytoplankton class i and κiis the ratio Knh/Kno for class i (κ < 1 as NH+4 is always preferred).103Symbol Definition Unit Value SourceKs Si half-saturation constant* µM N 0.6 tunedSi:N Si:N ratio* - 2.0 tunedchl:N Chlorophyll:N ratio µg (µmol N)−1 1.6 tunedchlr Chlorophyll reduction factor µM C d−1 0.2 Ianson and Allen (2002)Mesozooplankton fitZw Winter concentration µM N 0.0369 fit from Mackas et al. (2013)a Summer peak magnitude µM N 0.339 fit from Mackas et al. (2013)b Summer peak timing yearday 217.0 fit from Mackas et al. (2013)c Summer peak width d 114.0 fit from Mackas et al. (2013)Remineralization ratesrNH+4 NH+4 remineralization rate µM N d−1 0.0346 Denman (2003)rDON DON remineralization rate µM N d−1 0.1987 Jeffery (2002); Cugier et al. (2005)rPON PON remineralization rate µM N d−1 0.1987 Ianson and Allen (2002)rbSi bSi remineralization rate µM N d−1 0.2799 Cugier et al. (2005)Sinking ratesW1 Nutrient replete sinking rate* m d−1 0.5 Collins et al. (2009)W2 Nutrient depleted sinking rate* m d−1 1.2 Collins et al. (2009)wPON PON sinking rate m d−1 6.912 Jeffery (2002)wbSi bSi sinking rate m d−1 6.912 set to PONTable A.1: Biological model parameters. *diatoms only.The Si-limited uptake rate (silicic acid, H2SiO4) is given by the MM equationuis = Rimax(T )[H2SiO4]Kis +[H2SiO4](A.3)where Kis is the half-saturation constant for silicic acid for phytoplankton class i.uis is nonzero only for diatoms.The N uptake rates for each photosynthesizer class i are assigned based onthe limiting nutrient (or light). The rates are parsed such that ammonium is usedfirst until the uptake capacity, uinh, is reached, and then nitrate uptake fulfills theremaining N demand. The final N uptake terms becomeU iNH+4={min(uic,uis,uinh)when N is not limitinguinh when N is limiting(A.4a)104U iNO−3={min(uic,uis)−U iNH+4when N is not limitinguino when N is limiting(A.4b)The time rates of change of nitrate and ammonium concentrations are given by∂ [NO−3 ]∂ t =−3∑i=1(U iNO−3Pi)+RNO−3 −∂∂ z(we[NO−3 ])+V ([NO−3 ]) (A.5a)∂ [NH+4 ]∂ t =−3∑i=1(U iNH+4Pi)+RNH+4 −RNO−3 −∂∂ z(we[NH+4 ])+V ([NH+4 ]) (A.5b)where RX is the biological remineralization of nutrient X and Pi is the phytoplank-ton class i biomass concentration. The second to last term is the advective loss dueto estuarine circulation where we is the entrainment velocity, and the last term, V ,is the turbulent diffusion determined by the KPP model.PhytoplanktonThere are three classes of photosynthesizers included in the SOG biology model:diatoms, flagellates, and Myrionecta rubra. The time rate of change of photosyn-thesizer class i biomass is given by∂Pi∂ t =(U iNO−3+U iNH+4+ εiGi(Pf )−Mi)Pi−n∑j=1min[GiT j(Pi),Gij(Pi)]+∂∂ z (wsPi)−∂∂ z (wePi)+V (Pi) (A.6)where M is the natural mortality for phytoplankton class i, n is the number ofgrazers for Pi, and GiT j and Gij are the grazing loss terms for class i by the jthgrazer. The first vertical gradient term is the loss due to sinking (6= 0 for diatomsonly). The sinking velocity ws is parametrized in terms of a nutrient replete termW1 and a nutrient depleted term W2 (Collins et al., 2009).An additional grazing term εiGi(Pf ) is added to represent heterotrophy in M.rubra where εi is the M. rubra grazing efficiency and Pf is the biomass concentra-105Symbol Definition Unit Diatoms M. rubra Flagellates Microzoo PONPhytoplanktonRmax Maximum growth rate µM N d−1 5.184 2.160 2.160 - -Iopt Optimum light intensity W m−2 42.0 37.0 10.0 - -Knh NH+4 half-saturation constant µM N 2.0 0.5 0.1 - -κ Knh/Kno - 1.0 0.5 0.3 - -M Natural mortality µM N d−1 0.229 0.233 0.216 - -Microzooplankton grazingKµ Half-saturation constant µM N 1.0 1.0 1.0 0.5 2.0αµ Predation threshold µM N 0.3 0.5 0.4 0.3 0.6δµ Preference - 0.28 0.17 0.28 0.18 0.09Mesozooplankton grazingKm Half-saturation constant µM N 0.2 1.0 0.4 1.2 0.4αm Predation threshold µM N 0.0 0.2 0.0 1.0 0.0δm Preference - 0.25 0.2 0.2 0.2 0.15Table A.2: Tuned biological model parameters. Phytoplankton parametersare broken down by photosynthesizer class. Microzooplankton graz-ing and Mesozooplankton grazing parameters are broken down by preyclass.tion of M. rubra prey (flagellates). εiGi(Pf ) is nonzero only for M. rubra.ZooplanktonZooplankton grazing is modeled using a Holling type II response (Fasham, 1995).GiT j = ϒ j(Z j(t,z)(FT−αT j)KT j +FT j−αT j)( δ ijFiδˆ · Fˆ)(A.7a)Gij = ϒ j(Z j(t,z)(Fi−α ij)Kij +Fij −α ij)(A.7b)where ϒ j is the maximum grazing rate for predator j, α ij is the predator j predationthreshold for food source Fi, αT j is the predator j predation threshold for the totalfood abundance FT, Kij is the predator j half-saturation constant for Fi, KT j is thepredator j half-saturation constant for FT, δ ij is the preference of predator j for Fi,and δˆ j · Fˆ is the sum of products of food abundances and predator j preferences.106Symbol Definition Unit M. rubra Microzoo Mesozoo Sourceϒ Maximum ingestion µM N d−1 0.6912 1.296 1.201 tuned*KT Half-saturation constant µM N 1.0 1.25 1.0 tunedαT Predation threshold µM N 0.4 0.5 0.5 tunedε Grazing efficiency - 0.6 0.6 0.6 Cugier et al. (2005)M Natural mortality µM N d−1 - 0.0605 0.1987 tunedE Excretion µM N d−1 - 0.0043 0.1201 setTable A.3: Grazer parameters, broken down by grazer class. *Mesozoo-plankton maximum ingestion is from Collins et al. (2009).There are 2 zooplankton classes in the SOG model. The time rate of change ofmicrozooplankton biomass is given by∂Zµ∂ t = εµm∑i=1min[GiTµ(Zµ ,Fi),Giµ(Zµ ,Fi)]−(Mµ +Eµ)Zµ−n∑j=1min[GµT j(Zµ),Gµj (Zµ)]−∂∂ z(weZµ)+V (Zµ) (A.8)where εµ is the grazing efficiency, Eµ is the excretion rate, m is the number offood sources and Fi is the concentration of the ith food source. Following recentobservations, microzooplankton graze diatoms as well as flagellates (Strom et al.,2007).Average mesozooplankton biomass is fit to an annual Gaussian using a 20-yeartimeseries of zooplankton sampling in the Strait (Mackas et al., 2013; appendix B).Z¯m = Zw +aexp[−(D−b)2c2](A.9)where Z¯m is the average mesozooplankton biomass over the model domain, Zw isthe average winter biomass, and D is the yearday. The coefficients a, b, and c de-termine the magnitude, timing, and width, respectively, of each mesozooplanktonpeak.Mesozooplankton are assumed to choose their depth based on the profile oftheir prey, so their vertical distribution in the model is parametrized in terms of thetotal abundance of phytoplankton and the average mesozooplankton concentration.107Symbol Definition NH+4 DON PON Ref SourceσM Mortality waste fraction 0.000 0.475 0.475 0.050 setσE Excretion waste fraction 1.0 0.0 0.0 0.0 Jeffery (2002)σ rG M. rubra grazing waste fraction 0.0 0.5 0.5 0.0 setσµG Microzoo grazing waste fraction 0.0 0.5 0.5 0.0 setσmG Mesozoo grazing waste fraction 0.03 0.20 0.20 0.57 setTable A.4: Waste fractions from natural mortality, excretion, and sloppy graz-ing, broken down by waste destination. The sum of fractions for eachparameter add to 1. Diatom waste from natural mortality and grazingcontributes an additional whole portion to biogenic silica (bSi) detritus.Zm = Z¯mPP¯(A.10)where P is the total biomass profile of mesozooplankton prey and P¯ is the watercolumn average of P.DetritusThere are 3 classes of detritus in the SOG model: particulate organic nitrogen(PON), dissolved organic nitrogen (DON), and biogenic silica (bSi). The timerates of change of detritus concentrations are given by∂PON∂ t =WPON−n∑j=1min [GT j(PON),G j(PON)]− rPONPON+∂∂ z (wPONPON)−∂∂ z (wePON)+V (PON) (A.11a)∂DON∂ t =WDON− rDONDON−∂∂ z (weDON)+V (DON) (A.11b)∂bSi∂ t =WbSi− rbSibSi+∂∂ z (wbSibSi)−∂∂ z (webSi)+V (bSi) (A.11c)where rX is the remineralization rate for X and WX is the X fraction waste fromnatural mortality, grazing, and excretion processes. Sinking of PON and bSi isrepresented by the first vertical gradient terms in equations A.11a and A.11c wherewPON and wbSi are the respective sinking speeds. The waste term WX is the sum of108the X fractions of all waste-generating processesWX = σˆMX ·Mˆ+n∑j=1[σˆ jGX · Gˆ j(1− ε j)]+ σˆEX · Eˆ (A.12)where σˆYX is the array of X waste fractions from process Y , Yˆ is the array of wastevalues, and n is the number of grazers.Remineralization to dissolved inorganic nitrogen (DIN) from DON and PONfollows Jeffery (2002). DON and PON are first converted to ammonium at therates rDON and rPON. Ammonium is then converted to nitrate, and this conversionis proportional to [NH+4 ]2. The remineralization terms areRNH+4 = rDONDON+ rPONPON (A.13a)RNO−3 = rNH+4 [NH+4 ]2 (A.13b)Carbonate system parametersThe model equations for DIC and TA are∂DIC∂ t =[RNH+4 −3∑i=1((U iNO−3+U iNH+4+U iPC)Pi)]RC:N−∂∂ z(weDIC)+V (DIC) (A.14a)∂TA∂ t = (RP:N +1)3∑i=1(U iNO−3Pi)+(RP:N−1)(3∑i=1(U iNH+4Pi)−RNH+4)−2RNO−3 −∂∂ z(weTA)+V (TA) (A.14b)where U iPC is the uncoupled DIC uptake by phytoplankton class i under N limitationand RX:Y is the Redfield ratio of element X to element Y.The uncoupled DIC uptake term U iPC is given byU iPC = chlr(uic−UiNO−3−U iNH+4)(A.15)where chlr is the reduction factor for phytoplankton chlorophyll relative to carbon109(Ianson and Allen, 2002).Transfer of CO2 across the air/sea interface in the surface grid cell is parametrizedaccording to Fick’s second law of diffusionΦ=−kw(U10)K0∆pCO2 (A.16)where kw is the transfer coefficient as a function of the 10 m wind speed U10(Nightingale et al., 2000) and K0 is the CO2 solubility of Weiss (1974).110Appendix BModel TuningMost of the physical model parameters used in this study are as described in C-09.Several of those that differ from C-09 have been re-tuned to the STRATOGEMtemperature (T ) and salinity (S) dataset (table 2.1) more recently in AW-13. Es-tuarine circulation has three fitting parameters: the vertical velocity, the amountof freshwater added and the depth distribution of added freshwater. The verticalvelocity is set based on an entrainment calculation using Knudsen’s relations, theobserved link between Fraser River flow and surface S, and dynamics (C-09). Theamount of freshwater added is fit based on the STRATOGEM S observations, andthe depth distribution is set based on the observed halocline depth (AW-13). Theamount of mixing due to internal waves below the mixing layer was tuned to give agood match between model and observed halocline strength (C-09). The additionof flow from the northern SoG and its tuning is explained in Appendix C.The biology model was re-tuned to the STRATOGEM dataset following the ad-dition of multiple nutrient, phytoplankton, and zooplankton classes using a geneticalgorithm. Parameters chosen to tune were:• Table A.1 Chlorophyll:N ratio, Si:N ratio (diatoms), Si half-saturation con-stant (diatoms)• Table A.2 All parameters• Table A.3 All parameters except grazing efficiency, excretion, and mesozoo-plankton maximum ingestion.111The goodness of fit or cost function was constructed from the mean absolute errorbetween the model and observations for the following quantities:• Nutrients NO−3 and NH+4 concentration timeseries (optimized to within a10-day window) and total range• Phytoplankton Diatom, M.rubra, and flagellate biomass timeseries (opti-mized to within a 10-day window) and total range• Spring bloom Timing as defined by AW-13A good convergence was not found. It was determined that the amount of meso-zooplankton was too low to allow biomasses and nutrients in the model to matchthe observations, particularly during the summer. Zooplankton sampling duringthe STRATOGEM project yielded uncharacteristically low numbers for the Straitof Georgia ecosystem perhaps due to the years of sampling. The resulting lowmodel background mesozooplankton appears to give rise to phytoplankton over-productivity following the spring bloom.To improve the prescribed model mesozooplankton, I constructed a new pa-rameterization based on monthly averages from a 20-year timeseries in the SoG(Mackas et al., 2013). I only considered copepods larger than 1 mm and euphausi-ids (figure B.1 circles) as other classes of zooplankton were either more repre-sentative of model microzooplankton or present in such low abundance as to bedisregarded. The 1 to 3 mm copepod class is generally representative of Calanussp., Pseudocalanus sp., and Metridia pacifica whereas copepods greater than 3mm are dominated by Neocalanus plumchrus (Mackas et al., 2013). N. plum-chrus demonstrates an overwintering phase in the SoG from approximately Junethrough January during which it does not migrate to the surface to feed (Campbelland Dower, 2008). I assumed copepods greater than 3 mm sampled during thesemonths to be primarily overwintering N. plumchrus and excluded them.Euphausiids are generally omnivorous and carnivorously feed on primarilyother mesozooplankton. There is no mechanism for mesozooplankton self-grazingsince mesozooplankton are not explicitly modeled, thus I estimated the percentageof carnivorous euphausiids and excluded them from the mesozooplankton parame-terization. Bamstedt and Karlson (1998) estimate carnivorous feeding in Meganyc-112J F M A M J J A S O N D00.511.5YeardayBiomass [µM N]  Zm = 0.37 + 0.34exp[−(D − 217)2/1162]     r2 = 0.68STRATOGEMfitMediumCopepods1to3mmLargeCopepods3to5mmEuphausiidsTotalModified TotalFigure B.1: Monthly averages of 1-3mm copepods (open circles), 3-5mmcopepods (black circles), and euphausiids (gray circles) from 20 yearsof full water-column zooplankton sampling in the Strait of Georgia(Mackas et al., 2013). An annual least-squares gaussian fit (solid thickline) to the combined total (open squares) with estimated overwinter-ing copepods and carnivorous euphausiids removed (black squares)replaces the existing depth-average model mesozooplankton (brokenthick line) previously developed by C-09 using the STRATOGEM zoo-plankton data.tiphanes norvegica on the Swedish Skagerrak coast to be 40% or less. Stuart andPillar (1990) make a similar estimate (40%) for Euphausia lucens in the SouthAfrican Benguela current when phytoplankton abundance is high, however at lowphytoplankton abundance they estimate carnivory to be as high as 85%. The rangesof carnivory between the two studies are significantly different, and not surpris-ingly so considering the difference in euphausiid species and study areas. I esti-mated carnivory of SoG euphausiids to be 45% which is consistent with Stuart andPillar (1990) at moderately high phytoplankton abundance.The new mesozooplankton parameterization is based on a least-squares gaus-sian fit (figure B.1 solid bold line) to the sum of copepods and euphausiids (fig-113ure B.1 open squares) with overwintering copepods and carnivorous euphausiidssubtracted (figure B.1 filled squares). Compared to the C-09 and AW-13 parame-terization fit from the STRATOGEM data (figure B.1 dash bold line), the new pa-rameterization prescribes a similar mesozooplankton biomass in spring, but higherbiomass during the rest of the year with a peak concentration in late summer. Thisnew parameterization improved the model-data summer productivity agreement,but did not eliminate the model overproductivity.In order to further increase the grazing pressure on summer phytoplankton, Ireduced the mesozooplankton grazing preference for microzooplankton from 0.3to 0.2 and increased the mesozooplankton preference for flagellates from 0.1 to 0.2.The impact of this shift in grazing preference was two-fold. The added pressureon flagellates reduced their summer biomass, and the reduced pressure on micro-zooplankton increased their capacity for phytoplankton grazing. The final tuningof the biological model significantly improves the agreement with summer produc-tivity data, however the model still underpredicts summer NO−3 in the near-surface(figure 2.6c).114Appendix CNorthern Return FlowAs water enters the Strait of Georgia (SoG) from the Fraser River, it flows out in abuoyant plume spreading both north and south. As it spreads north, denser water,originally at the surface in the northern SoG, is subducted below the plume. Thusthe water entrained into the plume is not only intermediate layer intrusions fromthe Haro Strait region as was implicitly assumed in C-09 and AW-13.In summer, this subducted surface water is warm and low in nutrients. Evi-dence of this watermass can be seen in June 2004 STRATOGEM data (table C.1)as a southward flow of spicey (high temperature, high salinity for a given density)water from the northern SoG along the eastern coast of Vancouver Island (fig-ure 2.1). This water is near the surface (approximately 0-12 m) south of TexadaIsland but has subducted (approximately 8-18 m) by the time it has reached Gabri-ola Island (directly northwest of model site; figure 2.1). Water above is brackish,likely of eroded Fraser River plume origin and probably flowing northward. Theabsence of this intrusion in the C-09 model can be clearly seen in comparisons withthe STRATOGEM data (Figure C.2a).To include this surface-influenced water from the north, one needs parametersfor the temperature and the nutrient levels. However, in a one-dimensional modelsuch values are not available. Thus I assume that the northern values of temperatureand nutrients are correlated linearly with the southern values. Then I advect-in overa prescribed depth range, a lagged temperature/nutrient signal from the north.115LongitudeDepth [m]  125 W 40’ 20’ 124 W 40’ 20’5101520Spice−10−9−8−7−6−5−4−3Figure C.1: Observed spice contours from 0 to 20 m along the Strait of Geor-gia southwestern boundary from the July 2004 STRATOGEM cruise.Calculate the lagged temperature signal (nutrients are done similarly)Tin = (1−α)Tin +αTs (C.1)where α = dt/τ where dt is the model time step, τ = 2 days is the delay timescale,Ts is the surface temperature from the model and Tin is the time-delayed tempera-ture.Advect the temperaturedT (z) = dtSnwu (Tin−T (z))exp(− [z− zc]2W 2)(C.2)where dT is the northern return flow increment in the temperature equation, Sn =0.0206 m−1 is the strength of the northern flow, wu is the upwelling velocity, T isthe insitu temperature from the model, z is depth, zc = 23 m is the center of thenorthwater intrusion, W = Wu = 12 m is z is above zc and W = W` = 10 m if z isbelow zc.The coefficients for the fit (τ , Sn, zc, Wu, W`) were determined using a geneticalgorithm. Reasonable ranges for each parameter were chosen, 20 random sets ofparameters were generated, the fit to the data was determined, the top 6 were cho-116Depth [m]∆T (ModelNR off − Obsv.)(a)2003 2004 2005010203040∆T (ModelNR on − Obsv.)  (b)2003 2004 2005∆T [o C]−6−4−20246Figure C.2: Contoured difference between modeled and observed T as afunction of STRATOGEM cruise and depth (a) before and (b) afterincluding northern flow.sen, combined and mutated to generate 20 new sets and this was repeated until anoptimum parameter set was determined. Adding the lagged temperature signal tothe model greatly reduced the excess summer cooling seen at 5-10 m depth in theoriginal model (Figure C.2b) and reduced the root-mean-square error in tempera-ture over the 40 m domain from 1.3 to 0.74 ◦C.The values obtained using the genetic algorithm are close to what was expectedtheoretically. The time scale is slightly shorter; the initial guess was 5 days. Theshorter timescale implies much of the returned water is coming from closer thanoriginally considered. The depths are as expected, however they are deeper thanthe anomalies seen in Figure C.2. Because the water column is being upwelled,anomalies appear higher in the water column than the depth where the northernreturn water originally enters the system. The strength parameter is the horizon-tal velocity of the return flow, normalized by the upwelling flow speed and thehorizontal lengthscale of the subduction (approximately 50 km). The horizontalflow speed is expected to be smaller than, but on the order of, the estuarine flowspeed, v and by conservation of mass: v/wu = O(R/D) where R is the distancefrom the river to S3 (approximately 25 km) and D is the depth of the surface flow(approximately 12 m). This gives an estimate for Sn of less than but of order of0.04 m−1.The genetic algorithm was run only on the temperature data. However, the117influence of the northern return flow is also implemented for the nutrients (nitrate,ammonium and silicon), DIC and dissolved oxygen. Adding the northern returnflow improves the nitrate agreement with the data to a lesser degree (not shown).118Appendix DResidence TimeIn order to compare large sets of model runs, it is necessary to define model metricsthat condense outputs with high vertical and temporal resolution into single datapoints. Temporal averages are particularly useful here, but the choice of temporalaveraging window is not arbitrary. I choose to define these temporal averages interms of the model residence time, or the average time a parcel spends in the modeldomain before being removed or transformed.Consider the concentration θ of a scalar in a box of constant volume V . If aninflow F of concentration θs is added (and thus a mandatory outflow F of concen-tration θ is removed), the time rate of change of θ is given by∂θ∂ t =FV(θs−θ) (D.1)where V/F is the flushing time, or physical residence time, of the system (Bolinand Rodhe, 1973). For the biophysical model discussed in chapter 2, it is moreuseful to describe residence time as the convergence between scalar quantities intwo runs, θ i and θ j, which differ in initial conditions. In the case of the model,residence time is not solely described by flushing, but is instead determined by aset of sources and sinks summarized as the multivariable function F .∂θ i∂ t =F (. . . ,θi) (D.2a)119∂θ j∂ t =F (. . . ,θj) (D.2b)If we consider θ j to be a perturbation of θ i (θ j = θ i +∆θ ) then we can expandequation D.2b as a Taylor series.∂θ j∂ t =∂ (θ i +∆θ)∂ t =F (. . . ,θi)+∆θ ∂F (. . . ,θi)∂θ+12∆θ 2 ∂2F (. . . ,θ i)∂θ 2 + · · · (D.3)If we assume ∆θ to be small, which is a valid assumption where θ nears conver-gence, we can disregard the higher-order terms to produce an expression for ∆θ .∂ (∆θ)∂ t = ∆θ∂F (. . . ,θ i)∂θ (D.4)We now have an expression in the form of equation D.1∂ (∆θ)∂ t =∆θτ (D.5)where τ is the residence time of θ . Equation D.5 has an analytical solution∆θ = ∆θ0 exp{−tτ}(D.6)where ∆θ0 is the difference in θ between model runs at t = 0 and ∆θ <∆θ0/e is thethreshold for convergence. Thus the time required to reach convergence betweentwo models runs with difference initial conditions is the model residence time, orthe model spin-up time.In order to quantify the model spinup time, I initialized model runs (θ j) onall STRATOGEM sampling dates according to section 2.1.4 and tested them forconvergence against a control run (θ i) initialized on the first STRATOGEM sam-pling date. I chose six quantities to examine: two for model physics (T and S),two for the biology (NO−3 and chl a), and two for the carbonate chemistry (DICand TA). I used the water column maximum at each timestep and calculated the15-day running mean of that timeseries for the convergence tests to account for120Quantity Threshold UnitsT 0.2 degCS 0.2NO−3 1.0 µMChl a 0.5 mg m−3DIC 15.0 µmol kg−1TA 15.0 µeq kg−1Table D.1: Difference thresholds for T , S, NO−3 , chl a, DIC, and TA used fordetermining convergence between runs.high-frequency deviations in ∆θ .The threshold ∆θ0/e is not a sufficient convergence criterion for many sets ofmodel runs as ∆θ frequently grows or oscillates after reaching that threshold duringspinup. Thus, I adopt a more conservative definition of convergence for ∆θ thatI have verified graphically. I define convergence to occur when ∆θ falls below apre-defined threshold value for 30 days or more. I have defined threshold valuesfor T , S, NO−3 , chl a, DIC, and TA (table D.1).Model residence times are approximately constant for S, NO−3 , DIC, and TAacross the STRATOGEM initialization profiles with the exception of isolated spikes(figure D.1). Residence times are markedly longer for 2002 initializations, butthese times appear to be anomalies relative to the years that follow. Disregardingthe anomalies and spikes, these times average between 20 and 40 days (less forS). By contrast, residence times for T and chlorophyll a are seasonally-variableranging from the order of a week to 2 months (figure D.1a and d). Seasonality ofT residence time is a product of the northern return physics as similar runs withthe northern return parameterization disabled did not demonstrate similar seasonalvariability. Seasonality of chlorophyll a residence time arises due to oscillations inphytoplankton populations driven by grazing.Since most residence times are between 20 and 40 days and approximatelyconstant, I chose to use 30-day averages to compare model outputs and observedforcing quantities across multiple runs.121050100150 (a) TSpinup Time (days) (b) S050100150 (c) NO3−Spinup Time (days) (d) Chl a2003 2004 2005050100150 (e) DICSpinup Time (days)2003 2004 2005(f) TAFigure D.1: Model spinup times for (a) T , (b) S, (c) NO−3 , (d) chl a, (e) DIC,and (f) TA using each available set of initialization profiles from theSTRATOGEM dataset. Points are plotted according to the initialization(profile) date.122Appendix ECorrelation LagThe sensitivity of model pH, ΩA, and DIC fluxes to local forcing parameters andprimary productivity (PP) is evaluated using linear regression between two timeaverages. Determining the relationship between the averaging windows of the in-dependent and dependent quantities of the regression is not arbitrary, as delayedimpacts of changes in windspeed and freshwater fluxes in the water column couldoccur on the order of days to weeks. Collins et al. (2009) demonstrate strong corre-lations between spring bloom timing and averages of windspeed and cloud fractionoccurring more than a month earlier.In the present study, a lag must be determined for a correlation between anytwo quantities that both have variable time windows. pH3m and windspeed, for ex-ample, are two such quantities, while ΩAsur < 1 duration is not as it is a fixed modelrun metric. Since I evaluate the sensitivity of pH3m (and fluxes of DIC) through-out the entire year, I determine a unique lag for each day, defined as the numberof days between the independent and dependent 30-day averaging windows thatproduces the largest correlation coefficient (R2). I restrict the lag to a maximum of30 days as correlations do not significantly improve beyond this limit. Lag valuesthat accompany the pH regressions in figure 2.12 and the DIC flux regressions infigure 2.14 are presented in figures E.1 and E.2.123Lag (days)(a) pH3m lag0102030R2(b) pH3m correlation  J F M A M J J A S O N D J00. Qt CF PPintFigure E.1: Timeseries of (a) the lag between 30-day backaveraged pH3mand 30-day backaveraged U (bold), Qt (dash-dot), CF (solid), and PPint(patch) that maximizes R2, and (b) the corresponding optimized corre-lation of regression between the same metrics (as seen in figure 2.12).Lags denote the shift between the end of the forcing or PPint backaver-aging window and the plotted date. Lags are restricted to 30 days. Redlines and patches indicate anticorrelation while black lines and graypatches indicate positive correlation.124Lag (days)(a) ΦDIC U lag0102030(b) ΦDIC Qt lagR2(c) ΦDIC U correlation  M A M00.ΦDICmix(d) ΦDIC Qt correlation  J J A SΦDICfresh ΦDICPPFigure E.2: Timeseries of the lag between (a) 30-day backaveraged ΦmixDIC(bold) and ΦPPDIC (patch) and 30-day backaveraged U , and (b) 30-daybackaveraged ΦfreshDIC (dash-dot) and ΦPPDIC (patch) and 30-day backav-eraged Qt that maximizes R2. Timeseries of the corresponding opti-mized correlation of regression between the same metrics (as seen infigure 2.14) are included in (c) and (d). Lags denote the shift betweenthe end of the forcing backaveraging window and the plotted date. Lagsare restricted to 30 days. Red lines and patches indicate anticorrelationwhile black lines and gray patches indicate positive correlation.125


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items