Interactive comment on “ The sensitivity of estuarine aragonite saturation state and pH to the carbonate chemistry of a freshet-dominated river ”

I must admit that I am conflicted in making a recommendation of this paper. First, this is my favorite subject and I like the approach of a combination of data and modeling (also a combination of numerical model and simple mixing model). However I do not think the combination is successful. My dissatisfaction is further aggravated by the frustration in reading it at every step. First, the paper doesn’t present much new data. I believe most new data and the numerical model were published in their two earlier publications (authors really need to say what is new here). Second, the coupled physical and biological ROMS model (first part of the paper) and the simple mixing (second part) are not very compatible. I like that they want to extrapolate the study to a global discussion. However their “transport” of the Fraser conditions to other rivers are clearly


Introduction
Estuaries support productive ecosystems (Cloern et al., 2014) and significant human populations (Cloern et al., 2015).Critical trophic links within many of these ecosystems may be negatively impacted by increases in dissolved inorganic carbon (DIC) and reduced pH associated with ocean acidification (Haigh et al., 2015).Those organisms using the calciumcarbonate mineral form aragonite in their external hard parts (e.g., some molluscs) are especially vulnerable since oceanic CO 2 uptake lowers the aragonite saturation state ( A ) of seawater (Waldbusser et al., 2015).While carbonate system dynamics and/or acid-base chemistry have been studied extensively in marine (e.g., Jiang et al., 2015) and freshwater (e.g., Schindler, 1988) environments, less is known in the diverse estuaries where these two zones meet (Salisbury et al., 2008;Hagens and Middelburg, 2016;Cai et al., 2017).
Estuarine systems are complex, cover large salinity ranges that challenge current measurement techniques and are generally under-sampled (Ianson et al., 2016;Cai et al., 2017).The carbonate chemistry in the rivers that feed these estuaries may also be exceptionally variable (Cai et al., 1998), ranging from rivers with low DIC to total alkalinity (DIC : TA) ratios and high pH (> 8) like the Mississippi River (Hu and Cai, 2013) to blackwater rivers that have low pH (< 5) (Cai et al., 1998;de Fátima et al., 2013).Along a salinity gradient there may exist maximum sensitivity zones that are especially vulnerable to acidification.The strength and location of these zones depend on the river endmember carbonate chemistry (Hofmann et al., 2009;Hu and Cai, 2013;Xue et al., 2017).This sensitivity is generally expected to increase as the ocean Published by Copernicus Publications on behalf of the European Geosciences Union.Figure 1.Maps of (a) the Fraser River watershed showing the major geologic belts (Wheeler et al., 1991), and (b) the lower Fraser River and Strait of Georgia showing the model location (red star), Environment and Climate Change Canada (ECCC) meteorological (Sandheads and YVR, magenta symbols) and river gauging (Fraser River at Hope and Englishman River at Parksville, green symbols) stations used to force the model, and Fisheries and Oceans Canada (DFO) sampling stations (yellow circles) and the ECCC Fraser River Water Quality Buoy (cyan circle) used to constrain Fraser River TA and pH, respectively.absorbs CO 2 and becomes warmer (Hagens and Middelburg, 2016).
In addition, seasonality is often strong and single rivers and estuaries may experience highly variable conditions in space and time (Hellings et al., 2001;de Fátima et al., 2013;Waldbusser and Salisbury, 2014;Hunt et al., 2014;Voss et al., 2014;Ianson et al., 2016;Xue et al., 2017).Carbonate and silicate weathering are major sources of river DIC and TA in most rivers (Meybeck, 1987;Amiotte Suchet et al., 2003).Carbonates are concentrated globally in the northern mid-latitudes (Amiotte Suchet et al., 2003), and several carbonate-rich, mid-latitude watersheds demonstrate high TA and large TA flow dependence, such as the Changjiang and Mississippi rivers (Cai et al., 2008).TA flow dependence is also common among low-TA, low-latitude rivers like the Amazon (Richey et al., 1990) and Congo (Wang et al., 2013) rivers, although organic carbon contributions can complicate this behavior.In the Congo River for example, TA is flow-dependent but DIC is persistently high year-round (Wang et al., 2013).On top of natural variability, many estuaries also experience heavy anthropogenic pressure (Frankignoulle et al., 1996;Zhai et al., 2007;Cai et al., 2017), making them particularly vulnerable to ocean acidification (Cai et al., 2011(Cai et al., , 2017) ) and in some cases the subject of intensive management and policy initiatives (Fennel et al., 2013).
In the present study, we determine the sensitivity of A and pH in a large, mid-latitude, fjord estuary (Strait of Georgia, Canada) to changes in freshwater TA and DIC using a quasi one-dimensional biogeochemical model.This model mechanistically parameterizes horizontal two-dimensional estuarine flow as a function of river flow.Few carbonate data exist in the river plume region and fewer still at the river mouth, motivating us to perform these broad analyses.We establish a freshwater TA range and flow-dependent variability for the system by extrapolating TA observations from the Fraser River plume region to zero salinity.We use autonomous pH measurements from an Environment and Climate Change Canada (ECCC) mooring near the Fraser River mouth with our estimated TA to constrain freshwater DIC, since we lack reliable direct DIC observations in the Fraser River.From the model estuarine A and pH results across 18 river chemistry scenarios and 12 recent hydrological annual cycles, we identify regions of the freshwater carbonate chemistry range, expressed as both DIC : TA and DIC-TA, that produce enhanced estuarine pH and A sensitivity to this freshwater chemistry.We characterize these regions in terms of (past) low DIC and (future) high DIC freshwater carbonate chemistry and show how conditions in this temperate estuary deviate from theoretical mixing curves due to the strong local seasonal biological cycles.We further discuss the implications of these results for future climate with emphasis on the anticipated changes in hydrological cycles.

Study area
The Strait of Georgia (Fig. 1) is a large (∼ 6800 km 2 , > 400 m deep), temperate, semi-enclosed, fjord-like estuarine sea with strong seasonal stratification, productivity, and carbonate chemistry cycles (Moore-Maley et al., 2016;Ianson et al., 2016).This high productivity supports abundant populations of shellfish, finfish, and other higher organisms that may be sensitive to pH and A anomalies (Haigh et al., 2015).The Fraser River, the primary freshwater source, drains approximately 238 000 km 2 with seasonally variable discharge (∼ 800 to 12 000 m 3 s −1 at Hope, ECCC data, http://wateroffice.ec.gc.ca) due to summer snow/ice melt and lack of dams throughout most of the watershed.This large freshwater flux is partially contained by narrow passages and tidal mixing over sills (Pawlowicz et al., 2007), and thus imparts a significant freshwater influence on the strait, especially compared with regions where large rivers meet the ocean directly, such as the nearby Columbia River plume region (Roegner et al., 2011).These same coastal and topographic features create long residence times, causing carbon to accumulate and making the strait DIC-rich relative to the open ocean (Ianson et al., 2016) despite a strong, seasonal DIC upwelling signal over the outer shelf (Bianucci et al., 2011).The large freshwater footprint, together with the abundance of previous circulation (e.g., LeBlond, 1983;Pawlowicz et al., 2007), ecology (e.g., Masson and Peña, 2009;Allen and Wolfe, 2013), and acidification studies (e.g., Moore-Maley et al., 2016;Ianson et al., 2016) make the strait an ideal system for investigating the response of pH and A to freshwater carbonate chemistry in a complicated estuarine setting.
The seasonal progression of productivity in the Strait of Georgia begins with a characteristic spring phytoplankton bloom (Allen and Wolfe, 2013) followed by a shallow (∼ 20 m) surface layer of productivity throughout the summer that transitions into weaker fall blooms before returning to the background winter state (Harrison et al., 1983).Below this productive surface layer, the DIC-rich intermediate basin is persistently aragonite undersaturated.This signal is mixed to the surface throughout winter, but summer productivity maintains high pH and aragonite supersaturation in the top 20 m of the water column (Ianson et al., 2016).The strength and timing of this seasonal DIC cycle are strongly linked to local wind, irradiance, and freshwater forcing, the latter of which maintains the strongest influence of the three during summer (Moore-Maley et al., 2016).
The Fraser River watershed spans four distinct geologic belts (Fig. 1a) that transition from the carbonate-rich Foreland Belt to the silicate-rich Coast Belt (Cameron and Hattori, 1997).Carbonate and silicate weathering thus dominate the watershed (Voss et al., 2014); carbonate weathering generally produces TA faster than silicate weathering (Meybeck, 1987).Observed Fraser River TA varies strongly throughout the watershed, but generally accumulates along the Foreland Belt, decreases along the Coast Belt, and is highest at the low-flow stage (Voss et al., 2014).Fraser River TA thus appears to be produced primarily by carbonate weathering in the upper watershed, diluted by low-TA seaward tributaries, and flow-dependent.There are no data in the Fraser River region to date that determine organic acids and bases that may contribute significantly to TA, as they do in some coastal areas (Cai et al., 1998;Koeve and Oschlies, 2012;Kim and Lee, 2009;Hernández-Ayón et al., 2007) and rivers (Hunt et al., 2011;Kennedy, 1965)  mated to be as low as 10 % of TA in the Congo (Wang et al., 2013) and Kennebec rivers (Hunt et al., 2014).

Overview
In order to resolve the primary productivity in the strongly stratified Strait of Georgia, it is necessary to have fine vertical resolution, but it is not necessary to model the whole water column.Given the summer depth of the chlorophyll maximum of 5 m (Peña et al., 2016) and typical mixing layer depths of 1-7 m (Collins et al., 2009), we use 0.5 m vertical resolution and model the upper 40 m.The vertical model is located directly to the northwest of the Fraser River plume region (Fig. 1b).The region of the Fraser plume is dominated by estuarine circulation and wind-mixing.The dynamics have been well studied (e.g., Pawlowicz et al., 2007), giving us the information to effectively parameterize higherdimensional processes with a one-dimensional (1-D) model (Collins et al., 2009;Allen and Wolfe, 2013).The benefits of a 1-D model are the quick run times that allow us to simulate many parameter variations over multiple annual cycles (Moore-Maley et al., 2016).
Accurately simulating mixing of the stratified plume with the waters below is required to reproduce the biology and chemistry of the plume.To this end, the physical model is based on the KPP mixed layer model which includes the impacts of winds and heat/cooling on currents and mixing (Large et al., 1994).We add baroclinic pressure gradients and estuarine circulation to the model.More specifically, we add (1) freshwater and freshwater tracers to the mixing layer of the model, (2) a vertical upwelling due to entrainment and (3) a seaward advective loss to conserve volume flux.All three are defined in terms of the total freshwater discharge from the Fraser River and other small rivers (Collins et al., 2009).
The model represents a column of water with radius about one tidal excursion.It uses a 15 min time step.We explicitly model in situ temperature, ITS-90 (Preston-Thomas, 1990), and practical salinity, PSS-78 (UNESCO, 1981), as physical state variables.
DIC and TA are both explicitly modeled and are coupled to the biological growth and remineralization cycles (Moore-Maley et al., 2016).Transfer of CO 2 across the airsea interface in the surface grid cell is parameterized according to Fick's first law of diffusion (Sarmiento and Gruber, 2006), using gas transfer coefficient (Nightingale et al., 2000), Schmidt number (Wanninkhof, 1992), and K 0 solubility coefficient (Weiss, 1974) parameterizations.Model pH (total scale) and A are calculated from model DIC, TA, dissolved silica, temperature, salinity, pressure, and estimated phosphate using the CO2SYS program (Lewis and Wallace, 1998) and full salinity range K 1 and K 2 constants (Millero, 2010).Phosphate is roughly approximated (±1 µmol kg −1 , Riche, 2011) from model nitrate using the Redfield N : P ratio.Calcium ion concentrations, required for A calculation, are approximated by a linear regression to salinity (Riley and Tongudai, 1967), and by the mean observed calcium ion concentration near the Fraser River mouth, 350 µmol kg −1 (Voss et al., 2014), where salinity < 1.

Initialization and forcing
We use profiles of temperature, salinity, fluorescence, chlorophyll a, nitrate, and dissolved silica (1999 through 2012) measured near the model site (Pawlowicz et al., 2007;Masson, 2006;Masson and Peña, 2009;Peña et al., 2016) (Diane Masson, personal communication, 2014) to initialize the model (see Moore-Maley et al., 2016).Model runs are initialized in fall and run through a full year and then beyond until the end of the following December.Our analysis starts at the beginning of the year following the initialization date, which is always a longer period than the 30-day spin-up period.Initial phytoplankton, zooplankton, and detritus concentrations are determined according to Moore-Maley et al. (2016).Since few DIC and TA data are available, we use a representative fall profile (11 September 2011, Ianson et al., 2016) to initialize model DIC and TA.Model pH and A are not sensitive to initial carbonate chemistry conditions after the spin-up period.Time averages (fluorescence, chlorophyll a, nitrate, dissolved silica) and annual fits (temperature, salinity, DIC, TA) of the initialization data near 40 m are used to set the 40 m boundary conditions (Moore-Maley et al., 2016).Average model 40 m boundary conditions are summarized as the model seawater endmember in Table 1.
The model is forced at the surface (Allen and Wolfe, 2013) by wind stress calculated from hourly wind speed and direction observed at Sandheads weather station, and by heat fluxes derived from cloud fraction, air temperature and relative humidity observed at Vancouver International Airport (Fig. 1b; ECCC observations, http://climate.weather.gc.ca/ historical_data/search_historic_data_e.html; last access: 13 June 2018).Total freshwater flux (volume/time) into the Strait of Georgia is prescribed (Allen and Wolfe, 2013) using daily river discharge measurements obtained by ECCC (https://wateroffice.ec.gc.ca/search/historical_e.html; last access: 13 June 2018) in the Fraser River at Hope and in the Englishman River at Parksville (Fig. 1b).Englishman River discharge is used in this study as a proxy for the contribution of small, rainfall-dominated rivers to the freshwater budget of the strait (Collins et al., 2009).Heat and nutrient fluxes due to freshwater are prescribed (Moore-Maley et al., 2016) as concentration × flux.Model freshwater endmember concentrations are summarized in Table 1.

Evaluation
Previous studies using the model have evaluated it against physical, biological and chemical data.The vertical profiles of density and in particular the depth of the halocline are well represented (Collins et al., 2009).The model captures interannual variability in the biology and in the physics driving the biological variability -the model accurately predicts the timing of the spring phytoplankton bloom (Allen and Wolfe, 2013;Allen et al., 2016).The large seasonal variation of the carbon cycle is captured with some underestimation of DIC in the summer due to over-productivity (Moore-Maley et al., 2016), a common problem with coupled models in the Strait of Georgia (e.g., Peña et al., 2016).Evaluation of pH and A in the model shows both systematic and non-systematic errors (Moore-Maley et al., 2016).The variations in pH observations are well captured by the model, with a correlation coefficient of 0.80.There is a positive bias which is higher at high pH.The root-mean-squared error (RMSE) for pH is 0.16, but removing the systematic error (by fitting a line between the model and observations) gives a non-systematic or scatter RMSE of 0.06.One-third of this discrepancy can be explained by mismatches between model time and space.The evaluation is similar for A , a correlation of 0.79, a RMSE of 0.51 with a scatter RMSE of 0.18, and average mismatches due to time and space of 0.05.Thus, the nonsystematic error is sufficiently small to support the model's use for the process studies described here.Note that the systematic bias means that the model is conservative, overpredicting both pH and A .

Sensitivity analysis
In order to test the sensitivity of the biogeochemical model to changes in freshwater chemistry, we select ranges of freshwater TA and DIC endmembers based on observed TA in the Fraser plume and pH near the Fraser River mouth.We ex-trapolate estuarine TA data (Ianson et al., 2016, Table S1) sampled throughout the water column, following modern sampling and analysis procedures (Dickson et al., 2007), from seven recent (2010-2016) sampling campaigns near the Fraser River plume (Fig. 1b) to zero salinity using linear regression in order to establish freshwater TA endmembers across multiple seasons (Fig. 2a).We only consider profiles that include at least one TA sample below a salinity of 20 to allow sufficient river influence.DIC can vary at constant pH; however, DIC : TA cannot (given constant temperature and salinity).Thus rather than prescribing freshwater DIC directly, we define freshwater DIC : TA endmembers based on observations from the ECCC Fraser River Water Quality Buoy (Fig. S1 in the Supple-ment) moored approximately 10 km upstream of the river mouth (Fig. 1b).Buoy pH was measured potentiometrically using a regularly inspected (bimonthly to monthly), hullmounted, YSI ADV6600 multisensor and recorded hourly from 2008 until present (Ethier and Bedard, 2007).We calculate DIC : TA from buoy pH using the CO2SYS program (Lewis and Wallace, 1998) and full salinity range K 1 and K 2 constants (Millero, 2010).
There are uncertainties associated with the extrapolated freshwater TA endmembers due to our assumptions that TA is conservative across such a wide salinity range and that organic alkalinity contributions in the Fraser plume region are small.There are also significant uncertainties in the freshwater pH observations given that potentiometric pH measurements in freshwater are generally no more precise than 0.2 units (Covington et al., 1983).A thorough error analysis of these data is intractable given the paucity of data and our lack of additional parameters such as pCO 2 .Instead, given these uncertainties, we consider the extrapolated freshwater TA endmembers and the freshwater DIC : TA endmembers calculated from the buoy pH record to represent approximate seasonal ranges of freshwater carbonate chemistry rather than absolute values.
The freshwater TA endmembers span an approximate range between 500 and 1000 µmol kg −1 (Fig. 2).We thus define five constant freshwater TA (hereinafter TA f ) scenarios: three at the minimum, center and maximum of the extrapolated endmember range (500, 750 and 1000 µmol kg −1 , respectively) and two more to represent extremes beyond our estimated range at 250 and 1250 µmol kg −1 (Table 2).The extrapolated endmembers also vary seasonally, demonstrating significant hysteresis with respect to Fraser River discharge at Hope (low-pass filtered with a 40-day cutoff, Q filt , Fig. 2b) and a positive correlation with dQ filt /dt (Fig. 2c).We thus add an additional flow-dependent TA f scenario based on a linear regression of the extrapolated TA endmembers to dQ filt /dt given by where TA 0 = 750 µmol kg −1 , Q 0 = 840 m 3 s −1 and t 0 = 86400 s (Table 2).This flow-dependent TA f converges to TA 0 when dQ filt /dt = 0 and increases or decreases when dQ filt /dt becomes either positive or negative (e.g., Fig. 2c).Buoy pH is seasonally variable (Fig. S1) and changes are likely driven primarily by biological productivity since the correlation with river discharge is weak and summer river warming would tend to drive pH seasonal cycles in the opposite direction.The complete buoy pH record follows a Gaussian distribution with a median of approximately 7.5 and 1st and 99th percentiles at approximately 7.1 and 7.9, respectively (Fig. S1).Over the annual model freshwater temperature range (2.5-19.3• C) these pH values, in order from highest to lowest, correspond to average DIC : TA values of 1.032, 1.089, and 1.226.We use these freshwater DIC : TA values (hereinafter DIC f : TA f ) as our respective Low Carbon, Med Carbon, and High Carbon freshwater endmember cases (Table 2).The Low Carbon case is typical of present-day high TA rivers with low DIC : TA such as the Mississippi (Cai, 2003) and Changjiang (Zhai et al., 2007) rivers.We suggest that this scenario represents a lower limit in the Fraser River and may represent past chemistry (lower pCO 2 ).Likewise, we consider the High Carbon case to represent possible future pCO 2 increases but still within the range of present-day Fraser River pH observations (buoy pH).This High Carbon scenario still has a lower DIC f : TA f than low pH, weakly buffered rivers such as the Kennebec (Hunt et al., 2014) and Congo (Wang et al., 2013) rivers.These three DIC f : TA f cases combined with our six TA f cases produce 18 individual river endmember chemistry scenarios that we use in the biogeochemical model (summarized in Table 2).

Sensitivity to physical forcing
Freshwater forcing exerts strong control over model biology (nutrients, light, phytoplankton, zooplankton) and carbonate chemistry (DIC, TA, pH, A ), particularly during summer (Moore-Maley et al., 2016).Although the model is forced by a combination of the Fraser River and local rainfalldominated rivers, the Fraser accounts for most of the summer freshwater signal.The Fraser River flow record at Hope during the 12-year period from 2001 through 2012 is seasonally and interannually variable in terms of freshet size, timing, and duration, with the smallest freshet (in 2010) just under 6000 m 3 s −1 and the largest freshet (in 2012) just under 12 000 m 3 s −1 (Fig. 3a).Model salinity (3 m averaged) is strongly inversely related to this freshwater signal and reaches a minimum of approximately 15 in 2010 and 5 in 2012 (Fig. 3b).Winter rainfall pulses and their effect on winter salinity are evident particularly in February 2005.However, several winter salinity dips appear without noticeable corresponding pulses in the Fraser record.These events are likely forced by the scaled Englishman River record (not shown).Overall, aside from a handful of stronger winter salinity dips, a salinity threshold of ∼ 20 appears to separate the summer high flow period from the lower background flow regime.Dec Model DIC : TA, pH and A (3 m averaged) during 2010 and 2012 all demonstrate strong seasonal variability between winter and summer (Fig. 4) which can be attributed to the seasonal cycle of productivity that is characteristic of the region and persistent in the model (Moore-Maley et al., 2016).Prior to the summer, differences between the two years arise primarily due to variable wind and irradiance affecting the onset and termination of spring blooms.During summer, interannual variability of Fraser River discharge drives larger changes between the two years.In 2010, throughout the summer, biological drawdown of DIC remains strong, keeping DIC : TA low and pH and A high (red lines, Fig. 4).However, in summer 2012 DIC drawdown weakens due to river shading caused by river turbidity, thus increasing DIC and DIC : TA and significantly lowering pH and A (black lines, Fig. 4).Meanwhile, superimposed on top of this physicalbiological response, the effects of river chemistry begin to emerge.During the summer, model DIC : TA and pH vary across the six TA f cases (Fig. 4a and b, gray/pink lines) by approximately half of the total difference between years (black−red) at peak freshet.This variability across TA f cases is not evident for model A .
To investigate the impact of flow-dependent TA f , we consider the Medium Carbon, flow-dependent TA f case.The strongest TA f variations (∼ 100-400 µmol kg −1 ) occur dur-ing the summer (Fig. 5a) due to the rapid rising and falling discharge rates associated with the Fraser River freshet (Fig. 3a).Positive TA f deviations during rising water at the beginning of the freshet cause an overall decrease in pH and A (3 m averaged), although the timing of these events can vary depending on the freshet timing.For example, model pH and A decrease by 0.1 and 0.07, respectively, in June 2007 and May 2008 (Fig. 5b and c) due to the large TA f increase (Fig. 5a), and associated DIC f increase, during the rapid early progression to freshet in those years (Fig. 3a).In contrast, the freshet in 2012, although larger, progresses more slowly (Fig. 3a) and the model pH and A decreases are thus smaller and occur at different times in the season (Fig. 5b and c).Following the freshet, negative TA f deviations during falling water cause an overall increase in pH and A , although the increase is not as strong as the prefreshet decrease for pH (Fig. 5b).While these changes are significant (e.g., pH = 0.1 and A = 0.07 at times), they are small relative to the differences across the range of TA f scenarios (Fig. 4) and practically insignificant relative to the differences between low and high flow years (e.g., 2010 and 2012, Fig. 4).Furthermore, the model pH and A deviations caused by flow-dependent TA f are moderately symmetric about the freshet, so that time-averaged results will likely be similar to the results of a constant TA f scenario.Thus flow  dependence is likely less important than other factors when considering freshwater chemistry in this system; however, it could play an important role on daily to weekly timescales.

Sensitivity to river chemistry
In order to examine the response of the model across the range of TA f for all years and across all three freshwater carbon (DIC f : TA f ) scenarios, similar to our analysis of 2010 versus 2012 at Med Carbon (DIC f : TA f = 1.089,Fig. 4), it is useful to look at the model results as time averages.We identified a model salinity of 20 to be an approximate threshold separating the summer flow regime from the background flow during the rest of the year.We thus time average our results over the period where salinity < 20 using the same three metrics (3 m averaged DIC : TA, pH, A ) as in Fig. 4 and the same flow year color scheme as in Fig. 3 to produce a comprehensive summary of our 216 model runs (18 freshwater chemistry scenarios × 12 years of freshwater, wind and meteorological forcing) at low salinity (Fig. 6).The sensitivity of the model to differences in river discharge between years (e.g., red and black lines, Fig. 4) is clear in these low salinity time averages as the spread of points along the vertical axis at a given TA f , with 2010 and 2012 again highlighted (red and black stars, Fig. 6).The sensitivity of the model to TA f (e.g., pink and gray lines, Fig. 4) is also clear as the trend along the horizontal axis for a given year in each panel (Fig. 6).Time-averaged (salinity < 20) model DIC : TA increases with increasing TA f in all three river carbon cases (Fig. 6ac), causing a decreasing trend in model pH (Fig. 6d-f).These trends are weak in the Low Carbon scenario such that the response to TA f is weaker than the response to freshwater flow (Fig. 6a and d).Conversely in the High Carbon scenario, the increase in model DIC : TA and decrease in model pH with TA f are significantly stronger than the corresponding responses to freshwater flow (Fig. 6c and f).These results show that a high carbon freshwater endmember produces a stronger estuarine carbonate chemistry response than a low carbon freshwater endmember.The response of model A to freshwater chemistry is more complicated.First of all, the dominance of freshwater flow over freshwater chemistry in determining model A between 2010 and 2012 (Fig. 4c) is robust throughout the 18 freshwater chemistry scenarios as evident by the large vertical spread at a given TA f relative to the trends along the horizontal axis (Figs.6g-i).Secondly, model A increases with TA f in the Low Carbon scenario despite increasing model DIC : TA and decreasing model pH (Figs. 6a, d and g), because increasing DIC also increases the carbonate ion concentration in this case.Then in the High Carbon scenario, model A reverses its sensitivity to TA f to follow the model pH trend (Fig. 6f  and i).In this case the shift within the DIC pool toward dissolved CO 2 causes a sufficient reduction in the carbonate ion, so that increasing DIC no longer increases [CO 2− 3 ].

Mechanisms influencing sensitivity
To conceptualize the results, especially the non-intuitive ones (such as low pH at high TA f ), we calculate conservative mixing curves between model freshwater and seawater endmembers (Table 1) for six different river chemistry scenarios, specifically, at low TA f (solid lines) and high TA f (dashed lines) across the Low (magenta), Med (black) and High (yellow) Carbon scenarios (Fig. 7a-c).We compared model results (daily averages in salinity space) for the largest freshet year (2012) with these curves.There are strong similarities below S ∼ 15 in the way that estuarine DIC : TA, pH, and A are affected by river chemistry.For example, increasing carbon shifts both theoretical curves and model points in the same direction and by about the same amount (e.g., Fig. 7a-c).The pH shift between high and low TA f curves at  (d-f) pH, and (g-i) A under the (left column) Low Carbon, (center column) Med Carbon, and (right column) High Carbon scenarios between all constant TA f cases (horizontal axis) for each year in the full 12-year period (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012).The salinity < 20 averages (circles) are plotted in order of increasing freshet size, with the smallest (2010) and largest (2012) years highlighted as red and black stars, respectively, as in Fig. 3. Model DIC : TA (pH) generally increases (decreases) with increasing TA f .Model A shares this trend under the High Carbon scenario, but follows an opposite trend under the Low Carbon scenario.
High Carbon is about 0.2 for both theory and model (Fig. 7b).These variations with freshwater chemistry changes are analogous to the trends discussed in Sect.3.2 (Fig. 6).
The mixing curves illustrate that increasing TA f within a given DIC f : TA f scenario results in a greater contribution of DIC to the estuary relative to TA.This carbon increase is seen most clearly in the difference of the DIC and TA mixing curves, DIC−TA (Fig. 7d), but is also evident by comparing the DIC : TA curves between the low and high TA f cases (Fig. 7a).For DIC : TA, the endmember ratios are the same in both TA f cases, but the ratio within the estuary is different.The extra DIC in the estuary is what causes pH and A to decrease with increasing TA f .These decreases are strongest in the High Carbon scenario (e.g., Fig. 6f and i) because DIC f : TA f is larger, and thus the excess DIC contribution to the estuary is higher between TA f scenarios (yellow curves, Fig. 7d).
Model A responds differently to TA f under each of the three freshwater carbon scenarios, reversing its trend between Low Carbon and High Carbon.The mixing curves demonstrate a similar pattern (Fig. 7c) where the Med Car-bon scenario appears to define a threshold at which A diverges from its similarity to pH in response to TA f .The responses of model pH and A to freshwater carbonate chemistry changes in this system are thus asynchronous.We attribute this asynchrony to higher estuarine carbonate in the high TA f case which is able to buffer effectively against the equilibria-driven carbonate loss that accompanies rising estuarine DIC : TA, but only in the Med and High Carbon scenarios.In the Low Carbon scenario, the equilibria-driven carbonate loss is too strong.
While the trends of the model with changing freshwater scenarios are similar to those of the endmember mixing curves, the striking differences highlight the importance of biology and gas exchange in reducing unfavorable carbonate chemistry conditions in the Fraser River-Strait of Georgia system.The model converges tightly to the seawater endmember for most of the winter, but estuarine carbon decreases significantly during the spring phytoplankton bloom (green arrows, Fig. 7a).The model remains carbon deficient relative to the mixing curves throughout the year, only partially converging toward the freshwater endmember during  Theoretical conservative mixing curves between the model freshwater and seawater endmembers (Table 1) calculated using CO2SYS are also shown (lines) for the corresponding DIC f : TA f and TA f cases as indicated by color and line style, respectively (see legend).All runs demonstrate significant divergence from the theoretical curves during spring and summer while converging toward the freshwater endmember at low salinity during the freshet and toward the seawater endmember at high salinity during winter, as indicated by the annotations in (a).Differences between TA f cases in the model are otherwise similar to the conservative mixing cases; however, the effect of DIC f : TA f in the model is reduced relative to the mixing cases.S2) at the (a, c) low and (b, d) high ends of the reported ranges.The dashed lines are the DIC : TA 1 : 1 lines.The legend is in order of increasing latitude.Cyan symbols represent tropical watersheds, black symbols represent urbanized/polluted watersheds, yellow symbols represent temperate watersheds, magenta symbols represent watershed type proxies, and white symbols represent Arctic watersheds.The Fraser River DIC and TA endmembers used in this study (Table 2) are shown as a red star with error bars to represent the range of DIC f : TA f scenarios.DIC data are unavailable for the Columbia River and thus calculated using reported TA and pCO 2 ranges.DIC and pCO 2 data are unavailable for the Arctic rivers and DIC is thus calculated from reported TA and pH ranges.TA data are unavailable for the Maipo and Biobio rivers and thus calculated using DIC and pCO 2 .For most rivers, DIC : TA can be significantly greater than 1 (c, d) despite falling visually close to the 1 : 1 line (a, b).Uncertainty is not shown and not always reported.
the freshet.The model then retraces its carbon-deficient path toward higher salinity and only converges back to the seawater endmember after fall phytoplankton blooms terminate (blue arrows, Fig. 7a).Most importantly, while biology and gas exchange do not completely buffer the effects of freshwater chemistry in the model, they do shift the system toward a significantly lower carbon state than would be found under mixing alone.Fig. 8).Pollution drives part of this variability (black symbols, Fig. 8), as does latitude, in large part from the presence of carbonate rocks in many temperate drainage basins (Amiotte Suchet et al., 2003) -the tropical Amazon and Congo rivers have exceedingly low TA and DIC by contrast (cyan symbols, Fig. 8).In some rivers seasonal variability of DIC and TA far exceeds that of the Fraser.This variability is particularly high in the polluted rivers and also the Yukon.However, overall, the range of DIC : TA in most rivers is near or within the range determined for this study.The tropical rivers are an exception.The Congo River (cyan square, Fig. 8c and d), for example, has nearly constant DIC throughout the year, while TA drops to less than half of the DIC concentration at peak river discharge (Wang et al., 2013).Another exception is cold glacial meltwater (magenta square), which, even at a DIC : TA more than twice as high as our High Carbon scenario, is pCO 2 -undersaturated relative to the atmosphere and continues to take up DIC as it mixes with seawater (Meire et al., 2015).Seasonal biological productivity may also regulate DIC in some temperate rivers, at least near the river mouth, as it does in the Fraser River.

Comparison
In the Columbia River estuary, DIC drawdown is strong in the spring, but is replaced by net heterotrophy during the summer and fall.This early transition to heterotrophy makes the Columbia estuary a net annual source of CO 2 to the atmosphere (Evans et al., 2013) despite maintaining a relatively low DIC : TA ratio (yellow square, Fig. 8).The Fraser estuary by contrast has a significantly higher DIC : TA ratio, but maintains a seasonally persistent pCO 2 undersaturation relative to the atmosphere because of strong productivity throughout most summers.
In the coming decades, weathering rates and mean river TA are unlikely to change significantly due to climate (Riebe et al., 2001).In contrast, increasing atmospheric pCO 2 will increase DIC in the ocean and most likely in rivers as well.River DIC may be influenced by many additional variables, including changes in freshwater flow, human pressures (local anthropogenic inputs) and anticipated increases in river temperature.For example, the relatively high DIC : TA ratio of the present-day glacial DIC and TA endmember (magenta square, Fig. 8c, d) is in equilibrium with the atmosphere at a river temperature of 0 • C (Meire et al., 2015), but if this pure glacier water were to experience a 10 • C increase during its passage to the ocean, then outgassing would decrease river DIC by about 10 % (∼ 9 µmol kg −1 ) if the meltwater remained in atmospheric equilibrium.However, its pH would stay about the same (slight increase) in this scenario.More generally, ocean pH will become more sensitive to changes in dissolved CO 2 as the oceans become warmer and higher in overall CO 2 because of weakening buffer capacity (Hagens and Middelburg, 2016).Were this trend to be significant in the Strait of Georgia, it would further exacerbate the increased estuarine pH and A sensitivity that we observe at high river DIC in the biogeochemical model.
In addition to freshwater chemistry, the large seasonal freshwater flux from the Fraser River exerts a significant influence on estuarine pH, and particularly A , in the Strait of Georgia, with the highest flow years producing the most acidic and corrosive conditions.The strait is relatively acidic throughout the water column compared to its primary conduit to the Pacific Ocean, the Juan de Fuca Strait, as well as the British Columbia continental shelf during upwelling (Ianson et al., 2016).The surface Strait of Georgia is also relatively acidic compared with oceanic extrapolations of freshwater endmembers from the otherwise similar, glacially fed temperate Corcovado estuary in Chile (Torres et al., 2011), despite the higher extrapolated TA in the Fraser relative to the Corcovado estuary's primary freshwater source, the Puelo River.This excess background DIC may contribute significantly to the persistent surface aragonite undersaturation observed (Ianson et al., 2016) and modeled (Moore-Maley et al., 2016) in the strait during winter and during large summer freshets.
A warming climate has (Zhang et al., 2001) and will continue (Morrison et al., 2014) to reduce the peak freshet flows of glacial, temperate rivers and move the freshet timing earlier, which may reduce the severity of summer aragonite undersaturation in the Strait of Georgia associated with large freshets.Winter flows may also increase, however, which could increase the sensitivity of the estuary to river chemistry during these low flow times beyond our model (currently insensitive) predictions.A changing freshet climatology may affect the significance of freshwater TA flow dependence as well.While many different studies consider flow-dependent river TA to anticorrelate with discharge (e.g., Wang et al., 2013;Evans et al., 2013), the strong Fraser River TA hysteresis demonstrated by our data links the strongest Fraser River TA fluctuations to rapid pulses in river discharge.A future shift to smaller freshets will thus likely reduce the already weak influence of flow-dependent freshwater TA in the model, supporting the use of fixed freshwater TA endmembers in large-scale acidification projection models (e.g., Volta et al., 2016).However, the present study highlights the importance of choosing a freshwater TA endmember carefully.

Conclusions
The Fraser River is a large, free-flowing, glacially fed river that exerts a strong physical and chemical influence on a neighboring, semi-enclosed estuarine sea, the Strait of Georgia, that is DIC-rich relative to the open ocean.Based on recent data, we find that the Fraser River is moderately alkaline (TA = 500-1000 µmol kg −1 ) but appears to carry a significant DIC load (high DIC : TA) even relative to many world rivers.TA appears to vary systematically with river flow, but does not display a simple dilution relationship.Rather it exhibits a strong hysteresis such that TA is a function of the change with respect to time in river flow.
We examined the sensitivity of A and pH in the Strait of Georgia to Fraser River carbonate chemistry by summarizing the results of a coupled biogeochemical model across 12 hydrological cycles and 18 freshwater TA and DIC : TA combinations.Model results show that estuarine pH and A are strongly sensitive to river flow and river carbonate chemistry at moderate to low model salinities (< 20), and generally decrease as river flow or river DIC increases.The primary reason for this sensitivity is the strong estuarine DIC regulation achieved both by physical river shading effects on primary productivity, and by the DIC contribution relative to TA from the freshwater endmember to the estuary.
Model DIC and pH are tightly coupled by the dominant influence of the carbonate equilibria; however, estuarine A responds asynchronously from DIC and pH to increasing freshwater TA at low and moderate freshwater DIC : TA (which we consider to represent past and present) scenarios as river DIC : TA approaches unity.Under these conditions the carbonate ion follows DIC and simple dilution/enhancement controls A .This estuarine A response reverses when freshwater DIC : TA becomes larger (> 1.1) because the carbonate ion makes up an increasingly smaller portion of the total DIC and the A dynamics are thus no longer controlled by physical dilution.In our highest freshwater DIC : TA scenario, A responds to estuarine DIC : TA, similarly to pH.
The Fraser River estuary is biologically productive, which modulates its sensitivity to river chemistry.In winter, productivity is low, so estuary chemistry nearly follows simple endmember mixing theory.Once the seasonal phytoplankton bloom occurs and there is significant biological drawdown of DIC, both estuarine pH and A increase away from the physical mixing line.Thus, even though the river flow is at a maximum during summer, estuarine sensitivity to river chemistry (in particular to DIC) is significantly reduced.Finally, the strong impact of Fraser River flow input on pH and A in the Strait of Georgia will be reduced as the Fraser watershed makes its expected transition to smaller, earlier freshets with climate change.However, if DIC : TA ratios increase in one or both of the river-ocean endmembers as anticipated with rising atmospheric CO 2 , then the increased sensitivity of estuarine pH and A associated with this endmember change may counteract the flow regime changes to some degree.
Code and data availability.Model source code, documentation, run environment scripts, and analysis scripts are available from the UBC Salish Sea Bitbucket repository (https://bitbucket.org/account/user/salishsea/projects/SOG; UBC Mesoscale Ocean and Atmospheric Dynamics Lab, 2018).Results files for the sensitivity experiments are hosted by the UBC Research Data Collection as part of the Abacus Dataverse Network (http://hdl.handle.net/11272/10606;Moore-Maley et al., 2018).Model initialization data are available from the Fisheries and Oceans Canada Institute of Ocean Sciences Data Archive (http://www.pac.dfo-mpo.gc.ca/science/ oceans/data-donnees/search-recherche/profiles-eng.asp;Fisheries and Oceans Canada, 2018).Carbon cruise data in the Fraser River plume cited from Ianson et al. (2016) are available from the NOAA Ocean Carbon Data System (OCADS) at https://www.nodc.noaa.gov/ocads/data/0173514.xml (Ianson et al., 2018).Additional new observations from two field campaigns are also available from OCADS and can be found using the Ocean Carbon and Acidification Data Portal (https://www.nodc.noaa.gov/oceanacidification/stewardship/data_portal.html; NOAA, 2018).All other data presented in this article are available from their cited URL locations.
B. L. Moore-Maley et al.: Estuarine A and pH sensitivity to river chemistry Figure 2. (a)Observed TA versus salinity (S) (circles) and linear regressions (dashed lines) from seven Department of Fisheries and Oceans Canada (DFO) Strait of Georgia cruises in the Fraser River plume(Ianson et al., 2016, Table  S1in the Supplement), (b) TA extrapolated to S = 0 versus low-pass filtered (40-day cutoff) Fraser River discharge (Q filt ) at Hope (circles) plotted in sequence of year day (yd) as indicated by the dashed line and labels, and (c) TA extrapolated to S = 0 versus Q filt /dt (circles) with the flow-dependent freshwater TA parameterization used in this study (Eq. 1) overplotted (dashed line).Cruise ID numbers (legend) begin with the sampling year.Each cruise contains at least one data point at S < 20.Error bars (b, c) are the 95 % confidence intervals associated with the extrapolation to S = 0 using linear regression (a).
used in previous studies involving this model (e.g.,Moore-Maley et al., 2016) to maintain consistency with previous work and to ensure a wide range of climatological forcing regimes.With 18 possible combinations of TA f and DIC f : TA f , we ran the model a total of 216 times.Since the seasonality of DIC and TA in the Strait of Georgia is surface-intensified(Moore- Maley et al., 2016), we use surface (3 m average) DIC : TA ratio, pH and A as our primary model sensitivity metrics.

Figure 3 .
Figure 3. (a) Fraser River discharge at Hope versus year day and (b) model top 3 m averaged salinity shown for each year in the 12-year period considered in this study.The Fraser River flow regime is characterized by a prominent summer freshet associated with summer glacial melt and smaller rainfall-dominated winter spikes.Model salinity reflects this overall pattern.Individual years are plotted in order of increasing freshet (peak) size.The smallest (2010) and largest (2012) flow years are highlighted in red and black, respectively.

Figure 4 .
Figure 4. Time series of model top 3 m averaged (a) DIC : TA, (b) pH, and (c) A for all TA scenarios (Table2) during 2010 (red/pink curves) and 2012 (black/gray) curves under the Med Carbon scenario (DIC f : TA f = 1.089), colors consistent with Fig.3.The bold curves show model results at TA f = 750 µmol kg −1 and the lighter curves show model results at all TA f .All three quantities demonstrate a seasonal cycle consistent with biologically driven DIC decreases in summer, increases in winter, and variations in spring associated with the spring bloom season.A strong mid-summer DIC increase is present in the 2012 results when river discharge is strong, but absent in the 2010 results when river discharge is weaker.

Figure 5 .
Figure5.Results of the flow-dependent TA f case for each year in the 12-year record: (a) the flow-dependent TA f parameterization (Eq.1), and (b, c) differences between flow-dependent and constant (750 µmol kg −1 ) TA f runs at (DIC : TA) f = 1.089 for (b) 3 m averaged pH and (c) 3 m averaged A .Years of record are plotted in order of increasing freshet size, with the lowest (2010) and highest (2012) years highlighted in red and black, as in Fig.3.Since flow dependence is proportional to dQ filt /dt, the strongest effects surround the freshet generally with pH and A decreases prior to and increases after the freshet.

Figure 7 .
Figure 7. Daily model averages in salinity space during the largest freshet year (2012) of 3 m averaged (a) DIC : TA, (b) pH, (c) A , and (d) DIC-TA under the Low Carbon (DIC f : TA f = 1.032, magenta circles), Med Carbon (DIC f : TA f = 1.089, black circles), and High Carbon (DIC f : TA f = 1.226, orange circles) scenarios during two constant TA f cases: TA f = 500 µmol kg −1 (filled) and TA f = 500 µmol kg −1 (open).Theoretical conservative mixing curves between the model freshwater and seawater endmembers (Table1) calculated using CO2SYS are also shown (lines) for the corresponding DIC f : TA f and TA f cases as indicated by color and line style, respectively (see legend).All runs demonstrate significant divergence from the theoretical curves during spring and summer while converging toward the freshwater endmember at low salinity during the freshet and toward the seawater endmember at high salinity during winter, as indicated by the annotations in (a).Differences between TA f cases in the model are otherwise similar to the conservative mixing cases; however, the effect of DIC f : TA f in the model is reduced relative to the mixing cases.

Figure 8 .
Figure 8. DIC and TA literature values for selected world rivers (TableS2) at the (a, c) low and (b, d) high ends of the reported ranges.The dashed lines are the DIC : TA 1 : 1 lines.The legend is in order of increasing latitude.Cyan symbols represent tropical watersheds, black symbols represent urbanized/polluted watersheds, yellow symbols represent temperate watersheds, magenta symbols represent watershed type proxies, and white symbols represent Arctic watersheds.The Fraser River DIC and TA endmembers used in this study (Table2) are shown as a red star with error bars to represent the range of DIC f : TA f scenarios.DIC data are unavailable for the Columbia River and thus calculated using reported TA and pCO 2 ranges.DIC and pCO 2 data are unavailable for the Arctic rivers and DIC is thus calculated from reported TA and pH ranges.TA data are unavailable for the Maipo and Biobio rivers and thus calculated using DIC and pCO 2 .For most rivers, DIC : TA can be significantly greater than 1 (c, d) despite falling visually close to the 1 : 1 line (a, b).Uncertainty is not shown and not always reported.

Table 1 .
. Carbonate alkalinity is esti-Model freshwater and seawater endmembers used to prescribe river fluxes and values at the 40 m boundary.

Table 2 .
Biogeochemical model freshwater endmember cases run for each year in the 2001-2012 year range.