Rapid warming and salinity changes in the Gulf of Maine alter surface ocean carbonate parameters and hide ocean acidification

A profound warming event in the Gulf of Maine during the last decade has caused sea surface temperatures to rise to levels exceeding any earlier observations recorded in the region over the last 150 years. This event dramatically affected CO2 solubility and, in turn, the status of the sea surface carbonate system. When combined with the concomitant increase in sea surface salinity and assumed rapid equilibration of carbon dioxide across the air sea interface, thermodynamic forcing partially mitigated the effects of ocean acidification for pH, while raising the saturation index of aragonite (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varOmega_{AR}$$\end{document}ΩAR) by an average of 0.14 U. Although the recent event is categorically extreme, we find that carbonate system parameters also respond to interannual and decadal variability in temperature and salinity, and that such phenomena can mask the expression of ocean acidification caused by increasing atmospheric carbon dioxide. An analysis of a 34-year salinity and SST time series (1981–2014) shows instances of 5–10 years anomalies in temperature and salinity that perturb the carbonate system to an extent greater than that expected from ocean acidification. Because such conditions are not uncommon in our time series, it is critical to understand processes controlling the carbonate system and how ecosystems with calcifying organisms respond to its rapidly changing conditions. It is also imperative that regional and global models used to estimate carbonate system trends carefully resolve variations in the physical processes that control CO2 concentrations in the surface ocean on timescales from episodic events to decades and longer.


Introduction
Global ocean acidification (OA) proceeds as rising CO 2 levels in the atmosphere (CO 2ðatmÞ ) lead to higher oceanic carbon dioxide concentrations via uptake across the air-sea interface. In surface ocean chemistry, the term ''carbonate system'' refers to a combination of species produced by the equilibria The uptake of atmospheric carbon over time perturbs the carbonate system such that there is an increase in surface ocean CO 2 (CO 2ðaqÞ ) with a concomitant reduction in surface pH. Since the beginning of the industrial revolution, the world's surface ocean has decreased by about 0.1 pH units ), and further reductions on the order of 0.2-0.3 pH units are expected by 2100 ). Reduction in surface ocean pH due to increasing CO 2ðatmÞ is not the only effect on the carbonate system; additionally, OA causes reductions in carbonate ion (CO 2À 3 ) concentration and in the saturation states of calcium carbonate minerals (X) (Bates et al. 2014). While there is still debate on the direct role that CO 2À 3 availability and its proxy X play in shell development (Bach 2015;Jokiel 2016;Bednaršek et al. 2016), it is widely stated that reductions in CO 2À 3 and X represent a stressor to a variety of marine invertebrates that fix their shells or skeletons from calcium carbonate (e.g. Waldbusser et al. 2015). How these chemical changes will propagate into marine ecosystems is a subject of growing concern.
The potential threat by OA on marine organisms is of critical importance in the Gulf of Maine (GOM), where much of the value from fishery landings originates from potentially susceptible organisms such as lobsters (Homarus americanus) and sea scallops (Argopecten irradians), (Cooley and Kite-Powell 2009;Ekstrom et al. 2015). Declining pH, CO 2À 3 concentrations, and the saturation state of the shell building mineral aragonite (X AR ), can affect organisms in a variety of ways. The most cited impacts accrue through decreased rates of calcification, particularly during the larval phases of growth (e.g. Barton et al. 2012;Waldbusser et al. 2013), or via increased respiration that can consume energy required for mobility or reproduction . OA can also affect an organism's immune response, organ development, and olfactory discrimination (Ekstrom et al. 2015;Munday et al. 2009). Despite this knowledge, the effects of OA on individual species and community ecology are not well understood. While a number of studies have been initiated to fill this knowledge void, much of the work continues to be done within the context of controlled experimental studies rather than within functioning ecosystems. As such, it is difficult to assess a species' response to multiple stressors or mitigative factors that may occur in the natural environment (Breitberg et al. 2015).
In addition to OA, the carbonate system in coastal waters is affected by varying fluxes of Dissolved Inorganic Carbon ðDIC : P ðCO 2ðaqÞ þ HCO À 3 þ CO 2À 3 ÞÞ, total alkalinity (TA), and nutrients derived from local or remote sources. These processes are collectively known as coastal acidification and include acidic river discharge (Salisbury et al. 2008), atmospheric fluxes of acidic and alkaline compounds occurring predominantly in coastal regions (Doney et al. 2007), and coastal eutrophication. The latter is attributable to land-and atmospherically-derived nutrient fluxes that promote intense autotrophic production with subsequent CO 2ðaqÞ evolution and pH reduction via heterotrophic respiration (Cai et al. 2011).
Physical processes in the GOM (e.g., strong tides, wind-driven mixing, coastal currents) and large annual ranges in sea surface temperature (SST) and salinity generate significant thermodynamic variability in the carbonate system from diurnal to annual timescales. Of particular importance is the pronounced annual cycle of CO 2ðaqÞ , whereby disequilibrium with the atmosphere is partially balanced by an air sea flux of DIC (Shadwick et al. 2010;Vandemark et al. 2011). The GOM is therefore an ideal region to investigate how the effects of varying SST and salinity on the carbonate system relate to long-term trends driven by OA. In a climatological study using data from 1950 to 2013, the GOM records an annual SST range of 15:5 C and salinity range of 2.2 (Richaud et al. 2016). The annual range alone elicits a significant change in the carbonate system. For example, using approximate mean GOM salinity (32.2), mean TA (2184 lmol kg À1 ) and an atmospherically equilibrated seawater surface of pCO 2 (presently $ 400 latm), temperature alone produces an annual change of 0.013 in pH and 1.06 in X AR .
In this work we use simple data-driven decomposition models to explore how the carbonate system is affected by variability in SST, salinity, and its covarying carbonate parameter, TA. We show that over timescales of 5-10 years, such changes can partially mitigate or overwhelm the effects of OA. The variability imposed by thermodynamic changes is described in several texts (e.g. Butler 1991;Stumm and Morgan 1996). All carbonate dissociation and calcite mineral solubility constants are dependent on, and are modeled using temperature, and thus a change in temperature will alter the relative proportions of carbonate species.
Experimental determinations at constant TA, salinity, and DIC show that the sensitivity of the partial pressure of CO 2ðaqÞ (pCO 2ðaqÞ to temperature variability is 4:23% C À1 , a relationship that is consistent over a majority of the world's oceans (Takahashi et al. 1993). Similarly, at typical GOM SSTs, assuming no exchange of inorganic carbon with the atmosphere, and a constant mean GOM TA, salinity, and DIC, pH decreases in a nonlinear fashion at a rate of 0.015-0.017 U C À1 . However given the rapid equilibration timescales indicative of the adjacent Northwestern Atlantic (1-4 months; Jones et al. 2014), in reality the system would lose DIC during warming, which would mitigate the decrease, or even raise pH. For contrast, at constant mean GOM TA, and salinity, with pCO 2ðaqÞ held at 400 latm, pH increases slightly ( $ 0:001 U C À1 ) due to the release of DIC.
As water temperature increases, the equilibrium in Eq. 1 shifts to favor increases in CO 2À 3 and X AR (Dickson and Millero 1987). A dependency of X AR on temperature is also due to the influence of temperature on the apparent solubility product (Ksp). The Ksp for aragonite decreases by $ 0:4% C À1 (Mucci 1983). Over the range of observed GOM SSTs, assuming rapid equilibration, this translates into a X AR change of 0.05-0.09 U C À1 . Salinity variability has a more modest effect on the carbonate system of seawater by changing the ionic strength of the solution (Harris 2010), which decreases activity coefficients and, as a result, the values of pH, and X AR . Tracking salinity also enables regional modeling of the quasi-conservative TA (e.g. Lee et al. 2006;Cai et al. 2010) that can be used as one of two parameters needed to resolve the full carbonate system (Wolf-Gladrow et al. 2007).
Within the GOM, the decade spanning 2005-2014 was characterized by an extreme warming event with average SST increasing by over 0:2 C y À1 . Satellite observations of SST within the GOM show that the region was warming at a faster rate than 99% of the global ocean (Pershing et al. 2015), with the highest average annual values exceeding over 150 years of observations held in NOAA's Merged Land-Ocean Surface Temperature Analysis database (Vose et al. 2012). The warming initiated ecosystem changes that included a northward (or deeper) shift in the distributions of many planktonic and nektonic organisms as they sought out suitable temperatures (Nye et al. 2009;Fogarty et al. 2012). Enhanced warming was accompanied by an increase in salinity that is consistent with a change in water mass distribution related to a retreat of the Labrador Current and a northerly shift of the Gulf Stream as described in Saba et al. (2016) and Grodsky et al. (2017).
We describe the relationship between physical variability and the carbonate system in the GOM during a time span over which large changes in SST and salinity were observed. Our main emphasis seeks to quantify the effects that recent extreme environmental conditions played in driving changes in pH and X AR and to highlight the difficulties that variability in physical drivers may cause in resolving longer-term trends in OA. The work is structured as follows. ''Description of study region, methods, and data'' section describes the observations of SST, salinity, and carbonate parameters in the GOM, as well as the model framework and algorithms used to determine carbonate system variability. ''Results'' section examines the causes of variability in the carbonate parameters by decomposing the model into components that highlight atmospheric versus physical influences. We discuss the implications of the recent heating event and sub-decadal variability in terms of potential ecosystem stress and the longer-term trends in OA in ''Discussion'' section.

Site description
The GOM is a productive temperate continental shelf sea bounded by Cape Cod to the south and Nova Scotia to the northeast (Fig. 1). It is well known for its large semidiurnal tides and their resulting impact on mixing, and also for the high commercial value of its fish and shellfish landings. It is separated from the open northwest Atlantic by both Georges and Browns Banks. Considerable control on seasonal to interannual circulation patterns is exerted via shelf-sea exchange through the northheast Channel (NEC) which separates Georges and Browns Banks from the fresher coastal source waters on the Scotian Shelf (Townsend 1991;Pringle 2006;Hetland and Signell 2005;Geyer et al. 2004;Feng et al. 2016). The GOM has an average tidal range [ 3:0 m, and experiences large seasonal amplitudes in surface salinity (Geyer et al. 2004), net primary productivity (O'Reilly et al. 1987), SST and pCO 2ðaqÞ (Vandemark et al. 2011).
The key circulation feature impinging on the region is the Maine Coastal Current (MCC), which flows counterclockwise and delivers freshwater and constituents from the northeast along the coast and into the Gulf (Pettigrew et al. 2005).

Data preparation
We base our study on a combination of data sources that include output from a physical General Circulation Model (GCM) for salinity; satellite-derived SST and chlorophyll, and pCO 2ðaqÞ observed in the GOM. The carbonate system is modeled using thermodynamic equilibrium equations as described in the Handbook of Methods for Analysis of the Various Parameters of the Carbon Dioxide System in Seawater (Dickson and Goyet 1994): with the apparent dissociation constants where brackets represent the stoichiometric concentrations of the chemical species and [CO 2ðaqÞ *] represents the sum of the combined concentrations of CO 2ðaqÞ and H 2 CO 3 . A full description of the carbonate system requires salinity, temperature, pressure, and two carbonate parameters-in this study TA and pCO 2ðaqÞ . Total scale pH and X AR are subsequently estimated using the CO2SYS package (Lewis et al. 1998). The K1 and K2 constants are based on Mehrbach et al. (1973) and refitted by Dickson and Millero (1987), with the borate-to-salinity ratio of Uppström (1974). Alkalinity modifications by phosphate and silicate are assumed to be negligible. The saturation state with respect to aragonite is defined as where Ca 2þ and CO 2À 3 are the molar concentrations of calcium and carbonate ions in solution, and K sp is the solubility product of the mineral aragonite that is modeled as a function of in situ temperature, salinity, and pressure.

Data and algorithms
Our analyses are performed on time series data that begin in September 1981 and run through December 2014. All spatially resolved datasets are averaged over the GOM domain to generate integrated system-wide time series with a monthly temporal resolution. SST is based on the 1=4 daily Optimum Interpolation Sea Surface Temperature Version 2 (dOISSTv2), which combines observations from satellites, ships, and buoys into a blended product on a unified grid. The main data sources are NOAA's Advanced Very High Resolution Radiometer (AVHRR) 7-19 satellites. The resulting product typically has an average root-meansquare error (RMSE) of 0:3 C compared to buoy data (Banzon et al. 2016;Reynolds and Chelton 2010).
Salinity is derived from the Northeast Coastal Ocean Forecast System (NECOFS). NECOFS is an integrated atmosphere-ocean model forecast system The horizontal grid has a resolution (measured by the length of the longest edge of a triangular cell) that varies from 0.3 to 15 km over the entire domain. The resolution in the Bay of Fundy ranges from 0.5 km inside inlets to 2.0 km along the coast and 4.0 km in the interior of the Bay. NECOFS is a product of the Northeast Regional Coastal Ocean Observation System (NERACOOS), with support from the Massachusetts Fishery Institution and the MIT Sea Grant Program. The salinity time series is validated using data from five NERACOOS moorings in the GOM (Fig. 1). We use monthly averages for each buoy and perform a point-to-point comparison with each nearest grid cell in the model. We note that uncertainty of salinity (Table 1) can be considered conservative due to spatial and temporal mismatches between buoy and modeled data.
Time series of chlorophyll (Chl) are made using data from the NASA satellites' Coastal Zone Color Scanner (CZCS 1981(CZCS -1986 Zonal averages (40 N-50 N) of marine atmospheric boundary layer mole fraction of CO 2 (xCO 2 ) from 1981 to 2014 were acquired from the NOAA Earth System Research Laboratory GLOBALVIEW-CO2 database (Masarie and Tans 1995). Atmospheric CO 2 partial pressure (pCO 2ðatmÞ ) was estimated as pCO 2ðatmÞ ¼ XCO 2 (SLP-pH 2 O). Sea level pressure (SLP) was set to 1 standard atmosphere (1013.25mb), and water vapor pressure (pH 2 O) was calculated as a function of mean monthly GOM SST according to the empirical formula described by Cooper et al. (1998).
TA is modeled as a function of salinity and dOISSTv2 temperature data according to the North Atlantic model of Lee et al. (2006, Fig. 2). Model output was compared to all surface TA data found in the National Centers for Environmental Information (n ¼ 2140) producing an RMSE of 11:4 lmol kg À1 (Table 1). pCO 2ðaqÞ algorithm pCO 2ðaqÞ is modeled by deriving a relationship between SST, salinity, Chl, day of year, and observed pCO 2ðaqÞ based on a modification of methods presented by Signorini et al. (2013). The Signorini et al. algorithm performed reasonably well with modeled CO 2ðaqÞ showing an RSME 34:6 latm for the GOM region. Our algorithm is a multiple linear regression of a similar form, i.e.: where K 0 is the solubility for CO 2ðaqÞ and Chl is satellite chlorophyll. K 0 and Chl are system-wide mean values for the data set. yd 0 is the fit of a third order Fourier function that models a daily climatology of GOM pCO 2ðaqÞ based on day of year. The slope term represents the annual derivative of the atmospheric pCO 2ðaqÞ estimates for the North Atlantic region. These data have been smoothed by 3 months to account for the estimated time required for CO 2ðaqÞ to equilibrate into the surface ocean in our region (Jones et al. 2014).
Training data for the multi regression model of pCO 2ðaqÞ are taken from the Surface Ocean CO 2ðaqÞ Atlas (Bakker et al. 2016) within our domain. In an effort to avoid extreme pCO 2ðaqÞ values found in (Color figure online) regional river plumes (Cai et al. 2010;Salisbury et al. 2009), we remove all pCO 2ðaqÞ data measured at salinities \30:5. Further, to avoid disproportionate effects from extreme outliers, we remove data below and above the 1 and 99 percentiles respectively. The resulting data set contained 200 784 values. These data were averaged by month prior to regression analysis. Model coefficients, results, and uncertainty are found in Fig. 2 and Table 1.
Our approach differs from that of Signorini et al. (2013) in three ways. First, we substitute SST and salinity with K 0 , a function of temperature, salinity, and pressure (Weiss 1974). This approach provides a better linear relationship compared to the non-linearity that changing temperature imposes on CO 2ðaqÞ solubility. Second, instead of a simple sine function for day of year, we use yd 0 , which represents a better approximation of the annual pCO 2ðaqÞ cycle. Finally, we base the trend in CO 2ðatmÞ on observations instead of the linear increase of 1:68 latm y À1 used in Signorini et al. (2013). While the average linear slope determined here is similar (1:78 latm y À1 ), our approach accounts for the non-linear increase in CO 2ðatmÞ over time (Fay and Mckinley 2013).

Monte Carlo simulations
The uncertainties of X AR and pH are estimated via 10,000 Monte Carlo simulations using SST, salinity, TA, and pCO 2ðaqÞ data with accompanying uncertainty (standard deviation) values (Table 1). Prior to analysis, single-sample Kolmogorov-Smirnov tests were performed on each variable to verify the data were normally distributed.

Anomalies and sensitivity analysis
To better elucidate long-term trends and short-term variability obscured by the seasonal cycle, we use 2004 as a reference to which we compare all other data. This year is noteworthy, as it represents the beginning of the recent extreme warming event (Pershing et al. 2015). The modeled time series of TA, pCO 2ðaqÞ , salinity, and SST are used to construct time series of X AR and pH. Anomalies are calculated by removing mean monthly values in 2004 from all other years, month by month. We use the resulting time series to estimate the sensitivity of X AR and pH to variability imposed by OA, the effects of variable TA, and the combined effects of SST and salinity variability, the latter which are intended to track the variability of carbon dioxide solubility. For this work, we narrowly define OA as the changes to the carbonate system brought about by the increase in sea surface pCO 2ðaqÞ that are directly attributable to increasing CO 2ðatmÞ .
To understand the relative contribution of changes in pCO 2ðaqÞ , TA and combined SST and salinity on pH and X AR variability, we perform the decompositions shown in Eqs. 6 and 7. The partial derivatives quantify changes in pH and X AR attributable to incremental changes in pCO 2ðaqÞ , TA, and P, the last, which represents the combined temperature and salinity effects intended to track changes in CO 2ðaqÞ solubility (Eq. 8), with f representing the parameter of interest (i.e. pH or X AR ).
The sensitivity to each variable is calculated by holding each of the other components constant (

Results
The annual cycle of carbonate parameters Our results are affected by the combination of seasonal cycles, interannual variability, and decadal-scale changes. These different scales are evident in Fig. 2 (e.g. salinity), where full time series of the data used for model parameterization are presented. The seasonal cycles of these variabiles are detailed in Fig. 3, where the upper right panel highlights the annual cycle of SST, ranging from a low of \4 C in February and March to a mean high of 18 C in August. Salinity (upper left) has a range of $ 1:5 from a low of 31.0 in July to a high of 32.5 in November. The cycle of satellite Chl (bottom right) shows a peak in April, which corresponds to the maximum rates of net primary production (e.g. Behrenfeld and Falkowski 1997;Friedrichs et al. 2009). While not apparent in the climatology, a smaller, shorter duration, peak corresponding to a fall bloom, is often observed during September or October (e.g. Riley 1947;Townsend and Spinrad 1986;Thomas et al. 2003). This combination of seasonal variability in salinity, SST, and net production controls much of the annual cycle of pCO 2ðaqÞ in the GOM and is largely consistent with the seasonal dynamics described by Vandemark et al. (2011), who have suggested that changes in pCO 2ðaqÞ are controlled by the countervailing annual cycles of solubility and net biological process. The biological ''new year'' begins with an intense spring phytoplankton bloom that typically extends from March to May. By removing DIC, the spring bloom generates the most significant perturbation to pCO 2ðaqÞ throughout the year. Following the bloom, surface water warms, promoting an increase in pCO 2ðaqÞ . Changes in K 0 due to the summer warming are expected to increase X AR and decrease pH. The annual salinity cycle, which controls TA, imposes a modest modulation of pCO 2ðaqÞ via buffering of the carbonate system, leading to slightly higher pCO 2ðaqÞ values when the salinity is low and vice versa. A fall bloom may also cause a slight reduction in pCO 2ðaqÞ via uptake of DIC. As winter approaches, the water becomes well mixed, entraining deeper waters with higher DIC concentrations to the surface, raising the pCO 2ðaqÞ , while lowering X AR and pH (Wang et al. 2017).

Anomaly analyses
Our anomaly analyses (Fig. 4) show that pH (lower panel) exhibits a markedly different behavior than X AR (top panel). pH declines at an average rate of 0:0018 y À1 from 1981 to 2015, primarily in response to the OA signal imposed on the carbonate system by increasing atmospheric pCO 2ðaqÞ . The effect of OA on X AR is partially obscured because of their greater sensitivity to variations in SST and salinity. Both parameters show significant interannual variability caused by differential timing and magnitude of the seasonal cycle in SST, salinity, and net community production relative to 2004. It is notable that X AR shows a higher interannual variability (often [ 5%, lX AR ¼ 1:96) than pH (( 1%, lpH ¼ 8:06). The long-term trend in pH is not entirely consistent over the full time series, but varies on decadal scales with patterns in the SST and salinity data (see Fig. 4). For example, the time periods 1981-1987, 1995-2001, and 2004-2014 show considerable deviation from the long-term trend 1981-2015 driven by OA (i.e. pH ¼ À 0:018 and X AR ¼ À 0:065 per decade, Table 3, See OA). The decadal variability is more pronounced for X AR , which is dominated by the variation in TA and combined SST and salinity, with less dependence on the OA signal.

Sensitivity analyses
By modeling the carbonate system while holding certain terms constant, we are able to investigate in greater detail how specific processes affect carbonate parameters. Figure 5 shows results from these analyses. We note the annual amplitude of each curve is partially affected by our choice to hold variables relative to each month of 2004. While OA imparts a signal into each parameter, the effect is most significant for pH (Table 3). The trends in TA perturbations are significant when the entire time series is viewed (Table 3), but are even more pronounced when the data are divided into sub-decadal periods from 1991 to  (Table 4).
The combined effects of SST and salinity on pH produce small but significant changes over all timescales considered (Tables 3, 4) with the sign of pH Uncertainty is estimated as confidence intervals at p\0:05 and presented in parentheses. Note that the decadal change attributable to temperature and salinity variability is positive Fig. 5 The sensitivity of X AR (top panels) and pH (bottom panels), to OA (left panels), variable TA (center panels) and combined effects of variable SST and salinity (right panels). Slope information for various time intervals can be found in Table 4 Table 4 The response of X AR and pH to perturbations by acidification, variable TA, and combined variability of SSS and SST over specific time ranges with high variance in SST and salinity Confidence intervals are shown in parentheses as estimated at p\0:05. n/a indicates that the slope is not significant change being negative or positive depending primarily on the direction of the SST change. The combined effects of SST and salinity cause large variability in X AR at all time intervals except for 1991-1998 (Table 4). We note that the curves describing the sensitivity of X AR to TA and combined salinity and SST share similarities for two reasons. First, the TA values and their corresponding buffering capacity are dependent on salinity and, to a lesser degree, temperature (Lee et al. 2006). Second, as discussed in the introduction, decadal variability in SST and salinity often covary due to the interplay between fresher water entering from the north, and warmer waters from the south, with warming accompanied by an increase in salinity (and alkalinity) and vice versa. Both warming and increased alkalinity will increase X AR . Using mean annual values, we assess the respective contribution by acidification and other effects to recent trends in pH and X AR . Figure 6 shows a subset of the sensitivity analyses that focus on the extreme warming and salinity event spanning from 2004 to 2014. The figure demonstrates the magnitude of the OA signal relative to all other effects considered here (TA, SST and salinity). During this time OA is clearly the dominating factor affecting the trend in pH, whereas X AR is mainly influenced by the other factors (cf. Table 4). These distinctions can also be see in annual mean of the time series in Fig. 7.

Discussion
Over sufficiently long timescales it is understood that the oceanic carbonate system will respond proportionally to increasing DIC brought about by OA. Several recent studies have presented results showing OA to be the predominant driver of trends in pH (Lauvset et al. 2015), pCO 2ðaqÞ (Fay and Mckinley 2013;Turi et al. 2016), and X AR (Jiang et al. 2015). However, trends in these parameters are often found to be larger or smaller than that expected by OA alone. Sensitivity analyses by these investigators demonstrate that processes such as net warming, variable salinity, and upwelling affect the rate of change over time and space. This is particularly true in coastal regions (Turi et al. 2016), where upwelling and freshwater inputs can affect the variability of the carbonate system and increase the length of time series needed to detect statistically distinct differences in the rate of change (Fay and Mckinley 2013;Tjiputra et al. 2014).
Because of the climatological, geographical, and hydrological characteristics found in the GOM and its surrounding watersheds, we expect the OA signal to be modulated by physical forcing. This forcing can occur over sub-daily timescales such as the case in which heat flux or mixing of water masses rapidly change the distribution of carbonate parameters. Figure 6 shows an example occurring at the decadal scale where the OA signal dominates the variability of pH, but for X AR , the OA signal is overwhelmed by other factors. Over this time period, thermodynamic forcing partially mitigated the effects of OA for pH, while X AR increased by an average of 0.14 U. We note that physical forcing leading to a change in surface pCO 2ðaqÞ would also alter the disequilibrium of pCO 2ðaqÞ with the atmosphere, initiating a change in the rate of the net CO 2ðaqÞ flux. Although we could find no published information describing the timescales of mixed layer equilibration of CO 2ðaqÞ in the GOM, we assume that it is similar to or less than the adjacent Western Atlantic that exhibits timescales of 1-4 months due to the relatively high gas transfer velocities and moderate mixed layer depths found in the GOM (Jones et al. 2014;Galbraith et al. 2015). Such short equilibrium timescales would mean that increases in surface pCO 2ðaqÞ due to warming are accompanied with rapid removal of DIC across the air-sea interface while pCO 2ðaqÞ remains higher than pCO 2ðatmÞ . The opposite would be true during cooling. In the case of heating (increased pCO 2ðaqÞ ), the removal of DIC would raise the ratio of TA:DIC and in turn, modulate the change in pH and X AR to higher values. Thus we speculate that thermodynamic forcing with subsequent gas exchange plays a major role in regulating the status of the carbonate system in the GOM.
Our results demonstrate that pH and X AR in the GOM respond to OA, but are also sensitive to thermodynamic (net heating) and chemical (TA) variability. The magnitude of response to each process is dependent on the parameter and the conditions to which the parameters are subjected. For example, pH demonstrates significant correlation with OA over the full time series (r 2 ¼ À 0:95, p\0:001) (Fig. 4, Table 3). However, pH is less correlated with OA during periods with high overall variability (Table 4), owing to greater influences in variability of SST, salinity and TA. X AR is also affected by OA, but is far more sensitive to relative changes in temperature and salinity than pH (Tables 3, 4). Over the entire time range, and for each parameter, the effect on the carbonate systems by TA variability alone generates significant decdal-scale trends (Table 3). Since the distributions of TA closely follow salinity patterns, TA is primarily controlled by decadal scale salinity anomalies (Petrie and Drinkwater 1993)  linked to the upstream discharge of the Saint Lawrence River and the behavior of Arctic water masses (Khatiwala et al. 1999), as well as inputs of warm salty slope waters originating from the South (Saba et al. 2016;Townsend et al. 2015). This control has important implications in that discharge from local watersheds to the GOM has been increasing over the last four decades (Huntington et al. 2016;Huntington and Billmire 2014), and with higher precipitation than under current conditions expected in the future (Rawlins et al. 2012), local discharge could continue to increase. The possibility of fresher surface waters in the GOM from both regional and remote sources implies that TA could decrease together with the TA:DIC ratio, and as a consequence, pH and X AR would also decrease. The sensitivity to temperature shown by X AR highlights the degree to which interannual warming and cooling events influence longer-term trends. The last decade's warming and increasing salinity event has had a dramatic effect on X AR that is [ 2:5 times greater than that of OA alone (Fig. 6, Table 4). Put another way, the changes in salinity and SST linked to the recent event actually served to mitigate the equivalent of 25 years of decline from OA. We note that if this event had instead contributed fresher and cooler waters, the effect would have reversed and greatly exacerbated the impact of OA. While the perturbations in salinity and SST over the last decade are extreme in our data set, smaller SST and salinity fluctuations capable of affecting the consecutive interannual means of X AR ( [ AE 5%) can be found throughout the time series (Fig. 7). By contrast, such variability is an order of magnitude greater than the interannual variability in the offshore surface waters of the Atlantic and Pacific Oceans (Jiang et al. 2015).

The hidden signal of ocean acidification
While OA is clearly decreasing X AR and pH, it is difficult to observe this signal without sufficiently long time series (Fay and Mckinley 2013;McKinley et al. 2017;Henson et al. 2016). Because the strong seasonality and vigorous physical processes serve to attenuate the OA signal in the GOM, we seek to understand how long a time series of observations would need to be in order to observe the expression of OA. One can evaluate the degree to which total variability obscures the OA signal by estimating the time it would take for the current trend in OA to emerge from the background variability caused by biological and physical processes affecting the carbonate system. One simple approach, termed Time of Emergence (ToE), is designed to detect biogeochemical signals in the context of noise imposed by other natural processes (Keller et al. 2014). ToE is defined as where S is the absolute value of the slope of the trend of interest and r is the variance (standard deviation) of the observed time series. Doubling r implies that the signal emerges through twice the observed variability and that ToE will be estimated at the 95% confidence level. Assuming an average positive trend of 1:78 latm pCO 2ðaqÞ y À1 based on increasing CO 2ðatmÞ , our results for the OA simulations provide annual linear trends (S) for pH (À 0:0018 y À1 ), and X AR (À 0:0065 y À1 ), over the entire time range (Table 3). We calculate ToEs for pH and X AR based on imposition of the calculated slopes and r for the entire time series, as well as partial time series in 11-year increments starting in January 1982. After eliminating the partial year of 1981, 11-year intervals represent 1/3 of the dataset. When calculated over the entire time series using annual averaged data (cf. (Keller et al. 2014)), we find our estimate of ToE for pH in the GOM is longer (approximately 20 years vs 14 years). We speculate that the main reason for this difference is the result of using higher fidelity data to model pH, which imparts greater variance in the ToE estimate.
When based on monthly data and its full variance, longer ToE values are estimated (Table 5). We also find that the ToE for X AR is considerably longer than The standard deviation (r) of each time series is indicated in parentheses for pH, particularly for time intervals with high variance. Indeed, given the variability seen over the last 34 years, application of Eq. 9 to the data suggests it could take up to a century of observations for an OA signal to emerge in X AR . The differences in ToEs for different intervals are consistent with our findings that variability in pH is dominated by the OA signal, while X AR responds more strongly to perturbations arising from both OA and thermodynamic variability that operates over interannual to decadal timescales. Such results point to difficulties in trend analyses when using inadequately long time series or data of poor temporal resolution that are averaged annually. Further, the ToE estimate used here cannot provide a date at which the data set will reveal that OA is actually dominating the biogeochemical signal, but instead can only indicate the time it would take for the imposed slope to emerge through the variance. For example, trend analyses performed with the last interval of our data set (2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014) may conclude that X AR is steadily increasing, while over the longer time span it is actually decreasing in response to OA. The fact that the ToEs estimated for the GOM are longer than for than many ocean regions described in Keller et al. (2014) highlights the importance of understanding how local physical and biological processes affect the carbonate system. In our data X AR shows an annual range between 0.7 and 1.05 in our time series, which is more than twice as large as what is typically observed at open ocean sites such as the Hawaii Ocean Time Series ). The annual cycle of SST and salinity (with its attendant TA change), can account for most of the observed variability (see Fig. 2, top panels). When extreme events like the one experienced over the last decade are added to the seasonality in SST and salinity, the range of X AR values experienced by calcifying ecosystems becomes even larger. An important finding from these analyses centers on the need for long-term sustained observations in order to establish OA treands and to distinguish the different effects that chemical, physical and biological processes have on the observed signals and trends. Presently such time series are rare, with few (e.g. Hawaiian Ocean Time Series; Bermuda Atlantic Time Series) possessing the required data holdings necessary to resolve the drivers of carbonate system variability.

Potential significance to ecosystems
Questions arise concerning how rapid changes in the carbonate system may affect calcifying ecosystems. Many organisms have life histories that start in the water column prior to settlement in the benthos, and these early stages may have distinct sensitivities to surface carbonate parameter thresholds (Waldbusser and Salisbury 2014;Waldbusser et al. 2015). For example, Salisbury et al. (2008) have shown slower shell formation and growth in larvae of the commercially harvested GOM clam species Mercenaria mercenaria when X AR \1:6.
Sutton et al. (2016) investigated present-day coastal X AR distributions using pH and pCO 2ðaqÞ data from several buoyed assets, including the Coastal Western GOM Mooring, located at 43:02 N, 70:54 W. They also modeled preindustrial values, making assumptions about pCO 2ðatmÞ values of the past. One finding was that during preindustrial times, X AR never dropped below the 1.6 threshold at the GOM site. However, in present day conditions the threshold is exceeded in the coastal GOM 11-31% of the time from December through April, with peak exposure to low X AR in February and March. While the average SSTs during this time are typically below the 11 C necessary to initiate clam spawning (Ropes and Stickney 1965), continued warming in combination with OA could create conditions where larval shellfish are exposed more frequently to suboptimal X AR , leading to less growth. While it is beyond our scope to speculate whether the recent events could have affected ecosystem function, we note that during the last decade, surface ocean ecosystems in the GOM have been exposed to an X AR range that contains nearly the entire envelope of values observed over the 34-year time series (exposure range since 2004 = 1.2; time series max-min = 1.3). To assess whether an ecosystem or species is at risk or aided by such events, it is important to characterize the drivers and ranges of chemical conditions over periods of low and high variance, and to better account for organismal responses to such conditions. Our work highlights the importance of long-term monitoring of coastal ecosystems, especially those in ecosystems like the GOM that experience high variability at multiple timescales. Only time series data taken at regular intervals and over decades enable us to overcome signal-to-noise issues, thereby allowing the identification of extreme events and the separation of signals into those affected by physical, biological, and OA processes. Because long-term observation assets are costly to deploy and maintain, it is incumbent upon the ocean carbon modeling communities to continue efforts to resolve variations in the physical processes that control CO 2ðaqÞ on timescales from episodic events to decades and longer.