Including high-frequency variability in coastal ocean acidification projections

Assessing the impacts of anthropogenic ocean acidification requires knowledge of present-day and future environmental conditions. Here, we present a simple model for upwelling margins that projects anthropogenic acidification trajectories by combining high-temporal-resolution sensor data, hydrographic surveys for source water characterization, empirical relationships of the CO2 system, and the atmospheric CO2 record. This model characterizes CO2 variability on timescales ranging from hours (e.g., tidal) to months (e.g., seasonal), bridging a critical knowledge gap in ocean acidification research. The amount of anthropogenic carbon in a given water mass is dependent on the age; therefore a density–age relationship was derived for the study region and then combined with the 2013 Intergovernmental Panel on Climate Change CO2 emission scenarios to add density-dependent anthropogenic carbon to the sensor time series. The model was applied to time series from autonomous pH sensors deployed in the surf zone, kelp forest, submarine canyon edge, and shelf break in the upper 100 m of the Southern California Bight. All habitats were within 5 km of one another, and exhibited unique, habitatspecific CO2 variability signatures and acidification trajectories, demonstrating the importance of making projections in the context of habitat-specific CO2 signatures. In general, both the mean and range of pCO2 increase in the future, with the greatest increase in both magnitude and range occurring in the deeper habitats due to reduced buffering capacity. On the other hand, the saturation state of aragonite (Ar) decreased in both magnitude and range. This approach can be applied to the entire California Current System, and upwelling margins in general, where sensor and complementary hydrographic data are available.


Introduction
It has become increasingly apparent that upwelling systems, including the California Current System (CCS), are particularly vulnerable to anthropogenic ocean acidification due to their unique physical and chemical traits (Feely et al., 2008(Feely et al., , 2010;;Gruber et al., 2012;Hauri et al., 2013a, b).Upwelled waters have been isolated from the atmosphere and are naturally elevated in CO 2 from remineralization of organic matter; depending on the age of the upwelled water mass it may also contain anthropogenic carbon (Harris et al., 2013;Sabine et al., 2002).Recent observations estimate the saturation horizon with respect to aragonite (depth at which

Y. Takeshita et al.: Unique habitat acidification projections
Oregon border during a strong upwelling event (Alin et al., 2012;Feely et al., 2008;Harris et al., 2013).Furthermore, the rate of acidification (i.e., pH yr −1 ) is expected to be significantly higher along upwelling margins than observed in the surface open ocean (Bates et al., 2014;Gruber et al., 2012;Hauri et al., 2013b;Leinweber and Gruber, 2013;Rykaczewski and Dunne, 2010) due to the reduced buffering capacity of seawater at higher levels of CO 2 (Frankignoulle, 1994).This effect has caused parts of the CCS to venture beyond the envelope (defined as mean ±1 SD) from modeled preindustrial Ar conditions (Hauri et al., 2013b).This is a concern because organisms may need to survive outside of the environmental conditions to which they are acclimatized, and the evolutionary potential for key ecological species to adapt to such rapid and unprecedented changes is poorly understood.For example, a significant decrease in calcareous benthic organisms was observed along a natural pH gradient near a cold volcanic CO 2 vent (Hall-Spencer et al., 2008;Kroeker et al., 2011).However, some calcareous organisms such as limpets seemed to have adapted to higher CO 2 levels compared to corals and mussels in the same system (Rodolfo-Metalpa et al., 2011).Furthermore for upwelling margins, CO 2 covaries with other environmental stressors such as temperature and O 2 (Reum et al., 2015), making predictions more difficult due to potential nonlinear synergistic effects (Frieder et al., 2014).
A critical component in making accurate impact assessments of ocean acidification is the development of robust, ecosystem-specific projections of future CO 2 conditions (Andersson et al., 2013;Cai et al., 2011;Feely et al., 2009Feely et al., , 2010;;McNeil and Matear, 2008;Sunda and Cai, 2012).Development of surface, open-ocean acidification projections has been relatively straightforward, as they rely on welldefined chemical principles of CO 2 equilibrium at the airsea interface (Byrne et al., 2010;Lauvset and Gruber, 2014).Models become more complicated when attempting to resolve biological and physical processes that contribute significantly to the natural variability in the system.For example, biologically mediated "enhanced acidification" was identified in the northern Gulf of Mexico, causing significantly faster rates of acidification than the open ocean (Cai et al., 2011).On many tropical coral reefs, seasonal patterns in CO 2 are minimal, whereas the dominant frequency of variability occurs on diel and tidal frequencies (Hofmann et al., 2011).On upwelling margins, both biological and physical processes contribute to the observed natural variability in carbonate conditions (Fassbender et al., 2011).
One approach to develop region-specific ocean acidification projections is to apply an eddy-resolving regional ocean model system (ROMS) coupled with a biogeochemical component, as has been developed for the CCS (Gruber et al., 2011(Gruber et al., , 2012;;Hauri et al., 2013b).Such models have highlighted the importance of capturing physical and biological processes in highly dynamic upwelling systems.The model simulations show complex spatiotemporal variability (Hauri et al., 2013b), and predict that the frequency of "corrosive" upwelling events will intensify in both magnitude and duration by the year 2050 (Hauri et al., 2013a).However, the eddy-resolving ROMS project pCO 2 and saturation state conditions on a 5 km grid, whereas many marine animals experience the environment on the scale of centimeters to meters.In addition, regional models can largely resolve eventscale (weeks) and seasonal features, but cannot capture fluctuations on diel to tidal timescales, which can be the dominant frequency of variability in many near-shore environments (Duarte et al., 2013;Frieder et al., 2012;Hofmann et al., 2011).Due to this discrepancy in both space and time, numerical models tend to underestimate or entirely miss the high-frequency variability that exists for the microclimate of organisms.
To transition from region-specific to habitat-specific ocean acidification projections, high-temporal-resolution data from autonomous chemical sensors deployed across many habitat types can be used to directly quantify the full range of present-day carbonate conditions (Harris et al., 2013;Hofmann et al., 2011Hofmann et al., , 2014;;Martz et al., 2014;Sutton et al., 2014).The CCS supports many ecosystems that are of great ecological and economic value.In particular, there is great habitat and species diversity near-shore on the shelf, which includes a large number of commercially important invertebrates and fishes.Habitats in the region include bays and estuaries, rocky and sandy intertidal, eelgrass beds and kelp forests, sub-tidal reefs, canyons, and extensive sandy sea floor.There are many endangered and harvested benthic organisms that inhabit only one or a subset of these habitats.
Sensor data provide key observations for the mechanistic understanding of the controls on environmental conditions, and are particularly needed for coastal marine environments where complex physical and biological processes influence the observed variability.For example, sensor data from a near-shore kelp forest in the southern CCS revealed that local biological feedbacks and episodic upwelling events were the dominant drivers of CO 2 variability, with pCO 2 fluctuating by 600 µatm at 17 m water depth (Frieder et al., 2012).This scale of variability associated with near-shore environments is not captured by regional model simulations, but is most relevant for organisms living inside the kelp forest.
Here, we present an anthropogenic ocean acidification model to project CO 2 chemistry into the future by combining autonomous chemical sensor data, regional empirical relationships for the CO 2 system (Alin et al., 2012), hydrographic data, and the atmospheric CO 2 record (Keeling et al., 2005).This model was applied to four habitats ranging from the surface to 100 m water depth and all within 5 km of each other in the Southern California Bight (SCB).Each site showed distinct CO 2 variability signatures and acidification trajectories, highlighting the importance of interpreting ocean acidification projections in the context of present and future habitat-specific CO 2 signatures.Implications for future ocean acidification research are discussed.

Study sites
Moored autonomous sensor packages SeapHOx or SeaFET (Bresnahan et al., 2014) were deployed at four depths (4, 17, 30, and 88 m) within several distinct habitats on the San Diego continental shelf for 1 year starting June 2012 (Fig. 1).All sensors were deployed near the seafloor; the three shallowest sensors were deployed within 3 m of the bottom, and the deepest sensor was moored 12 m above the seafloor (Table 1).
A SeaFET was deployed at the Ellen Browning Scripps Pier 2 m above the benthos as a part of the Scripps Ocean Acidification Real-time Monitoring Program.The sensor was located approximately 400 m from the shore in the surf zone.Weekly discrete samples for total alkalinity (TA) and dissolved inorganic carbon (DIC) were taken alongside the sensors for calibration and quality control following standard protocols (Dickson et al., 2007).The sensor was serviced every 1-2 months to remove biofouling organisms.
The La Jolla kelp forest is part of the South La Jolla State Marine Reserve and is characterized by a dense population of Macrocystis pyrifera.The chemical variability in this ecosystem is strongly influenced by regional physical processes (e.g., upwelling and stratification) and local biological feedbacks (e.g., production and respiration).A SeapHOx was deployed at 17 m in the southern portion of the kelp forest, 3 m above the bottom.The reader is referred to Frieder et al. (2012) for further details on site and deployment description.
The La Jolla canyon is a submarine canyon plunging from approximately 20 m to a depth of 1000 m within several kilometers from shore.A SeapHOx was deployed over a sandy bottom at the southern canyon edge at 30 m depth within the Matlahuayl state marine reserve; the O 2 sensor malfunctioned and thus is not included.The water in the La Jolla canyon is characterized by higher salinity and lower temperature, O 2 , and pH (data not shown).Tidal energy in submarine canyons is significantly amplified (Swart et al., 2011), bring-ing deep water from the canyon to the canyon edge.Therefore, physical forcings are the dominant drivers for chemical variability at this site (Navarro et al., 2013).
The Del Mar buoy was first deployed in 2006 at 100 m off of Del Mar in northern San Diego at the shelf break, and has provided continuous time series data (e.g., temperature, salinity, oxygen, and current) at discrete depths (Frieder et al., 2012;Send and Nam, 2012).A SeaFET sensor was deployed on the mooring at 88 m in 2011, and has provided a near-continuous time series of pH since.Colocated sensors include temperature, salinity (SBE 37), and dissolved oxygen (O 2 ; Aanderaa optode).Water at this depth is isolated from the atmosphere and below the euphotic zone, and thus influenced primarily by upwelling and tidal dynamics.

Cruise data
Hydrographic data were collected aboard R/V Melville during the student-led San Diego Coastal Expedition cruises in June/July and December 2012 (Fig. 1).The SCB is characterized by relatively weak (compared to the northern CCS) but nearly year-round upwelling.However, there is a clear seasonal cycle based on climatological data, where upwelling intensifies generally between April and August, with the maximum occurring in May (Bograd et al., 2009).The cruises therefore corresponded with upwelling (June/July) and non-upwelling (December) seasons.Water samples were collected at stations ranging from > 100 km from shore at 1200 m water depth, to within 5 km from shore at 30 m depth.Discrete samples were analyzed for O 2 , pH, and DIC; duplicate samples were collected during every cast.
Discrete samples for O 2 were collected and analyzed by titration using a custom-built system (Martz et al., 2012).The titrant was standardized prior to and after each cruise using KIO 3 standard solutions prepared in house (Fisher, lot 105 595); no detectable drift was observed for either cruise.Precision was ±0.6 µmol kg −1 (duplicate n = 62, 1 SD), and the accuracy was estimated to be ±0.5 % because KIO 3 standards were not recrystallized (Emerson, 1999).
Samples for DIC and pH were collected in 150 or 250 mL Pyrex serum bottles (13 mm neck) following standard procedures (Dickson et al., 2007).However, rather than leaving headspace, the bottle was filled completely, and a gray butyl stopper was inserted to prevent gas exchange; samples were analyzed within 4 h of collection.
DIC samples were analyzed using a custom-built system based on an infrared analyzer (LI-COR 7000) similar to systems built by others (Friederich et al., 2002;O'Sullivan and Millero, 1998).The DIC measurements were calibrated using certified reference materials provided by the Dickson Lab at SIO by applying a gain correction (slope) and assuming an offset of zero (intercept).The reference materials were stored in CO 2 -impermeable bags (3 L Scholle DuraShield ® ) and were measured frequently throughout the cruise.The stability of the reference material in the bag was verified  (Clayton and Byrne, 1993) using an automated system (Carter et al., 2013).The pH is reported on the total hydrogen ion concentration scale.The indicator dye (m-cresol purple, ACROS lot A0264321) was used as received from the manufacturer without further purification.An offset was applied based on measurements in certified Tris buffer provided by the Dickson Lab.The precision and accuracy of the measurements were estimated to be ±0.0015(duplicate n = 86, 1 SD) and ±0.02 (Liu et al., 2011), respectively.TA and pCO 2 were calculated using CO2SYS (van Heuven et al., 2011), with pH and DIC as inputs and carbonic acid dissociation constants from Mehrbach et al. (1973) refit by Lueker et al. (2000).

Sensor data
The SeapHOx and SeaFET sensor packages utilize a modified Honeywell Durafet III pH combination electrode for high-frequency pH measurements (Martz et al., 2010).These sensor packages have been successfully deployed in ecosystems worldwide (Frieder et al., 2012;Hofmann et al., 2011;Kroeker et al., 2011;Martz et al., 2014;Price et al., 2012), and have been shown to have excellent stability in seawater for months to years (Bresnahan et al., 2014).The SeapHOx is an integrated sensor package that also consists of an Aanderaa 3835 oxygen optode and a Sea-Bird SBE 37 conductivity-temperature sensor all plumbed into a pumped flow stream; the SeaFET measures pH using a passively flushed cell.Sampling frequencies were 1 h −1 or greater at all depths.
All pH measurements were calibrated based on discrete TA and DIC samples taken alongside the sensor, at minimum at the beginning and end of each sensor deployment (n > 4 for every site), as recommended by the best practices (Bresnahan et al., 2014).The resolution of the pH measurements is better than 0.0005 pH, stability is estimated to be better than 0.005, and accuracy is estimated to be ±0.015.Sensors were removed periodically for maintenance, but all were deployed for > 50 days during both the upwelling and relaxation season.
At the surf zone (surface waters), a constant TA value of 2240 µmol kg −1 was assumed, since discrete TA samples showed low variability (2240 ± 7 (1 SD) µmol kg −1 , n = 57).For the three subsurface sensors, TA was estimated (TA est ) using a regional empirical relationship developed for the CCS, with temperature and salinity as inputs (Alin et al., 2012); an offset of +8 µmol kg −1 was applied to TA est based on comparisons to discrete samples collected (root mean squared error (RMSE) = 6 µmol kg −1 , n = 25).This offset was a persistent feature over multiple years (2010)(2011)(2012), thus most likely reflecting a regional surface TA influence that is not incorporated in the empirical relationship developed for the whole CCS.DIC, pCO 2 , and Ar were calculated using CO2SYS (van Heuven et al., 2011) with pH sensor data and TA est as inputs.Uncertainty for the calculated DIC, pCO 2 , and Ar is pH-dependent but is on average estimated to be ±13 µmol kg −1 , ±25 µatm, and ±0.04, respectively.The daily range of sensor data was calculated by first high-pass-filtering the data with a 36 h window and then taking the difference between the daily maximum and minimum.The mean daily range was then calculated by averaging the resultant time series.

Approach
The carbonate conditions were modeled by decreasing or increasing DIC while using TA conditions from 2012.Modeled projections were made for preindustrial times and for year t, where t ranges between 2012 and 2100.The model presented here is based on the C * t approach (Gruber et al., 1996), but instead of using tracers to estimate the age of the water mass (e.g., CFCs), we used the atmospheric CO 2 record as a quasiage tracer.The age of the water mass ranged between 0 and 50 years in this study region.Although this approach must be used with caution, we demonstrate that our estimates are in good agreement with previously published anthropogenic carbon inventory estimates using age tracer measurements in this region (Feely et al., 2008;Sabine et al., 2002).In this model, it was assumed that ocean acidification is due to anthropogenic CO 2 invasion through the air-sea interface alone.We also assumed that both the path of a particular water mass between the subduction and upwelling site and the rate of remineralization processes remain unchanged.Sensitivity to these assumptions is explored in the Discussion.
The DIC of the modeled year t (DIC t ) is calculated by where DIC 2012 is the DIC observed in 2012, and DIC anth is the additional anthropogenic CO 2 that the water mass would have absorbed since 2012.Different formulations for DIC anth were used for surface waters (i.e., above the seasonal mixed layer depth, defined here as σ θ ≤ 25.2 kg m −3 ) and subsurface waters (σ θ > 25.2 kg m −3 ), and are outlined below.
For surface waters, DIC anth was calculated as the difference in surface DIC between year t and 2012.Surface DIC was calculated by assuming atmospheric equilibrium with TA = 2240 µmol kg −1 (based on water samples from the Scripps pier) and using pCO 2,atm projection under the 2013 IPCC RCP6.0 scenario (Hijioka et al., 2008).Although large deviations from equilibrium conditions are often observed in the coastal ocean due to upwelling and biological production (Hales et al., 2005), the mean pCO 2 calculated from sensor data at the surf zone was 394 ± 43 (1 SD) µatm (Table 2), suggesting that the surface water at the study site was near atmospheric equilibrium.
For subsurface waters, DIC anth was quantified as the increase in DIC due to anthropogenic CO 2 when the water parcel was last in contact with the atmosphere.The mass balance of DIC for subsurface waters is where DIC • is the preformed DIC, and DIC bio is the DIC added by remineralization processes in the ocean interior.DIC • can be expressed as the sum of DIC if it were in equilibrium with the atmosphere (DIC eq ) and the degree of air-sea disequilibrium due to slow gas exchange kinetics and biological processes ( DIC diseq ): Since anthropogenic CO 2 only enters the ocean at the surface, the increase in DIC eq represents the anthropogenic ocean acidification signal, DIC anth , assuming DIC diseq is invariant with time.However, in order to use this approach, the age of the water parcel must first be quantified, as this determines the pCO 2,atm with which it was last in contact.
The age of the water parcel was established by combining Eqs. ( 2) and (3): where the superscript denotes the year at which the water parcel was last at the surface (i.e., equal to 2012 minus age of the water mass).The age of the water mass was calculated by comparing the atmospheric CO 2 record to the pCO 2,atm that is necessary to generate DIC was calculated from Eq. ( 4).The DIC anth for subsurface waters was modeled for each projection year as a linear function of σ θ , and the surface and subsurface DIC anth were connected assuming a two-endmember linear mixing between σ θ 25.2 and 25.5 kg m −3 to prevent step changes (Fig. 2).

Calculation of DIC bio
DIC bio was quantified following formulations in Sabine et al. (2002): where AOU (apparent oxygen utilization) = (O 2,sat − O 2,obs ), TA • is the preformed alkalinity, and the r's are the elemental remineralization ratios (Anderson and Sarmiento, 1994).
The oxygen saturation concentration (O 2,sat ) was calculated using the equations in Garcia and Gordon (1992), and TA • was estimated based on historical near-surface TA data in the Pacific (Eq. 3 in Sabine et al., 2002).Phosphate concentrations necessary to estimate TA • were not directly measured but instead estimated from a regional empirical relationship using historical data (Supplement); the uncertainty in estimating phosphate using this approach propagates to an error in TA • of 4 µmol kg −1 .

Estimation of DIC diseq
Making accurate estimates of DIC diseq is important because it is a source of large uncertainty for anthropogenic carbon inventory calculations (Matsumoto and Gruber, 2005).
Traditionally, age of the water mass is quantified using tracers such as CFCs and then the DIC diseq is subsequently calculated (Gruber et al., 1996;Sabine and Tanhua, 2010).However, such tracer measurements were not made for this study.Alternatively, we estimated DIC diseq based on θ and S data to overcome this limitation (Sabine et al., 2002).The mean θ and S between σ θ of 25.5 and 26.5 kg m −3 were 10.0 • C and 33.9, respectively, resembling water type 1e in Sabine et al. (2010) with a corresponding DIC diseq = −6.24µmol kg −1 , the value used in this study.

Calculation of the age of water parcel
In order to estimate the age of the water mass, we use Eq. ( 4) to calculate DI C

2012−age eq
, which is the DIC of the water parcel that was in equilibrium with the atmosphere when it was last at the surface (i.e., equal to 2012 -age of the water mass).Therefore the age of the water mass can be calculated by comparing the atmospheric CO 2 record to the pCO 2,atm that is necessary to generate DIC 2012−age eq .The latter was calculated from the fugacity of CO 2 of the water mass when it was last in contact with the atmosphere at the time of subduction (f CO 2012−age 2,eq ), assuming 100 % relative humidity and a barometric pressure of 1 atm (Dickson, 2007).The year that the water parcel subducted was determined by matching the calculated CO 2,atm to the mean annual CO 2,atm record (Keeling et al., 2005); the age is the difference between 2012 and the calculated year (Fig. 3).A relationship between σ θ and the age was established by fitting a secondorder polynomial to the subsurface data (n = 186, R 2 = 0.92) and assuming the age of the surface water (σ θ < 25.2) to be 0 (Fig. 3).The nonzero age of the water that appears around σ θ = 24.4kg m −3 corresponds to the shallow oxygen maximum layer that formed during the summer.However, since this density range is still shallower than the seasonal mixed layer, its age was considered to be 0. The age of the water ranged between 0 and 50 years between σ θ of 25.2 and 26.5 kg m −3 .

Estimation of preindustrial DIC
In order to calculate the preindustrial DIC, Eq. ( 3) is written as where DIC anth represents the anthropogenic carbon present in the water parcel in 2012, and DIC prein eq is the DIC of the water parcel if it were in equilibrium with pCO 2,atm = 280 µatm.Combining Eqs. ( 2) and ( 7) and rearranging gives Calculated DIC anth as a function of σ θ is shown in Fig. 4.
Note that the values calculated here are in good agreement with published values using age tracers (Feely et al., 2008;Sabine et al., 2002) for higher σ θ but are significantly higher at lower densities.This is because the literature values were quantified using offshore subsurface waters, whereas our study region is near the coast along an upwelling margin, where subsurface waters are brought near the surface and are thus affected by surface processes.The agreement at higher density where surface influence is minimal demonstrates that the model presented here is capable of making accurate estimates of anthropogenic CO 2 .Furthermore, the comparison illustrates the importance of incorporating surface influence when making acidification projections in shallow, coastal ecosystems.Preindustrial DIC (DIC prein ) was calculated by subtracting DIC anth from DIC observed in 2012.Preindustrial pCO 2 , Ar , and pH were calculated using DIC prein and TA conditions from 2012.

Carbonate chemistry variability observed in 2012
The results are presented using pH, pCO 2 , or Ar , since pH was directly measured and pCO 2 and Ar are commonly used as stress indicators for respiration (Brewer and Peltzer, 2009) and calcification (Langdon et al., 2010), respectively.Across all four sites in 2012, pCO 2 increased with depth.In 2012, the mean pCO 2 in the surf zone was near atmospheric equilibrium (394 µatm), while the mean pCO 2 at 88 m was 878 µatm (Table 2) and reached a maximum of 1270 µatm.
The variability in pCO 2 also increased with depth (indicated by the SD of the time series), which was only 43 µatm in the surf zone but was 149 µatm at 88 m depth.The mean (±1 SD) Ar decreased with depth; the mean Ar in the surf zone was 2.4 ± 0.25, whereas it was 1.05 ± 0.18 at 88 m.Undersaturated conditions ( Ar < 1) were observed 48 % of the time at 88 m in 2012 but were not observed at other sites.However, unlike pCO 2 , the variability in Ar , indicated by the SD, decreased with depth (0.25 at the surface to 0.18 at 88 m) (Table 2; see Discussion).The mean pH decreased with depth; the mean pH in the surf zone was 8.05, whereas it was 7.73 at 88 m.The variability in pH increased with depth until 30 m, but decreased at 88 m (Table 2).
Distinct, habitat-specific CO 2 signatures were observed at the four deployment sites (Figs. 5, 6 and 7).Here, we define habitat-specific CO 2 signatures as how CO 2 conditions varied in that habitat, regardless of biological or physical origin.In the surf zone, the conditions were near atmospheric equilibrium, with intrusions of higher pCO 2 waters through internal tidal bores, a common feature observed in shallow, upwelling environments (Booth et al., 2012;Pineda, 1991); temperature and pH were correlated during these events (Supplement Fig. S2).This leads to a high, mean daily range of CO 2 conditions (e.g., 96 µatm, 0.46, and 0.085 for pCO 2 , Ar , and pH, respectively) (Table 3).However, the signature from the internal bores usually only lasted several hours, and remained at near atmospheric equilibrium for the majority of the time.Higher occurrence of tidal bores was observed during the spring and summer months relative to winter (Fig. 5), consistent with previous observations (Pineda, 1991).The mean diel range of pCO 2 at the surf zone was significantly higher than measurements made by a surface mooring located offshore in the SCB (Leinweber et al., 2009).
The mean (±SD) CO 2 conditions in the kelp forest and the canyon edge were similar, and the SDs for Ar and pH were the highest among the sites (Table 2).However, the timescales of the variability were different, indicating that distinct processes control the CO 2 conditions in these two habitats.For example, the mean daily range for all the variables was significantly higher at the canyon edge compared to the kelp forest (Table 3).Submarine canyons are known to amplify tidal energy (Navarro et al., 2013;Swart et al., 2011)  occurred on semi-diurnal and diurnal cycles, indicative of tidal forcing.Temperature and pH were correlated on these shorter timescales (Fig. S2), further supporting the fact that the variability was dominantly driven by intrusion of cold, deep waters from the canyon.While tidal forcings and daily biological production are drivers for carbonate chemistry in the La Jolla kelp forest, the largest variability occurred on event timescales (Frieder et al., 2012); event timescales are defined as longer than a day but shorter than several weeks.For example, pH, pCO 2 , and Ar regularly changed by up to 0.3, 250 µatm, and 1.3 on event timescales, more than 3 times the mean daily range.Variability on event timescales is due to a combination of changing water mass, stratification, and biological respiration (Frieder et al., 2012).In addition, a clear seasonal pattern was observed at the canyon edge, where higher pCO 2 and lower pH and Ar were observed during the spring and summer months (upwelling season) and lower pCO 2 and higher pH and Ar were observed during the fall and winter (relaxation season).Due to incomplete data coverage, a seasonal trend at the kelp forest and pier sites could not be discerned.The largest seasonal change among the four sites for Ar (∼ 1) was observed at the canyon edge; pCO 2 differed by roughly 200-300 µatm between the two seasons.
The shelf break experienced the highest mean CO 2 conditions and the highest and lowest SD for pCO 2 (149 µatm) and Ar (0.18), respectively (Table 2); the SD of pH (0.070) was lower than the kelp forest (0.083) or the canyon edge (0.075).Variability on tidal, event, and seasonal timescales were observed at this site (Figs. 5, 6, and 7), as has been previously reported for oxygen (Send and Nam, 2012).In general, upwelling on event timescales led to greater changes in pCO 2 and pH than on tidal frequencies (Figs. 5 and 6).The largest variability for all parameters occurred between the seasons, where changes in pCO 2 , pH, and Ar were approximately 350 µatm, 0.2, and 0.5, respectively.The close proximity of these four sites demonstrates the wide variety of habitat-specific CO 2 signatures that exist over a small spatial scale, especially in near-shore environments.

Modeled carbonate chemistry
Each habitat showed distinct trends in both modeled mean and variability in pCO 2 , pH, and Ar owing to increased levels of anthropogenic DIC ( DIC anth ) (Table 2).For example, the mean pCO 2 at the surf zone (4 m), canyon edge (30 m), and shelf break (88 m) increased by 225, 435, and 738 µatm, respectively, from 2012 to 2100; this drastic difference in increased mean pCO 2 is driven by different buffer factors due to depth differences among the sites.The increase in variability (i.e., SD) was also larger at 88 m (97 µatm) compared to the surf zone (37 µatm), although the largest increase occurred at the canyon edge at 30 m (126 µatm) (Table 2).Similar trends were observed for mean pH, where the largest mean decrease in pH occurred at 88 m water depth (0.26).However, the SD increased into the future for the three shallowest sites, whereas the SD decreased at the shelf break (88 m).In contrast, the largest decrease in the mean Ar occurred at the surface relative to the deeper sites, whereas the decrease in range was equivalent across all depths.
The measured and modeled time series for pCO 2 and Ar at the shelf break for the year 2012 and 2100 are shown in Figs. 8 and 9.The variability in pCO 2 increases on both seasonal and tidal timescales; the seasonal amplitude increases from approximately 350 to 650 µatm and the mean daily range increases from 110 to 325 µatm by 2100.This greater variability is in addition to an increase in mean pCO 2 of > 700 µatm.On the other hand, the variability in Ar on both seasonal and shorter timescales decreases.Furthermore, the shelf break is projected to experience undersaturated waters over 90 % of the time by 2060, compared to 48 % in 2012.Similar patterns were observed in the kelp forest as well, where both the mean conditions and variability in pCO 2 increased and Ar decreased (Fig. 10).The largest variability at the kelp forest occurred on timescales of days to weeks, and high-frequency (< 1 day) variability was significantly smaller than at the shelf break.Therefore benthic organisms at the kelp forest would experience elevated CO 2 conditions for prolonged periods of time, with only intermittent exposure to near-atmospheric conditions.
Preindustrial pCO 2 and Ar were compared to conditions observed in 2012 (Table 2).At most sites, the observed pCO 2 , pH, and Ar in 2012 were already outside of their preindustrial variability envelopes (defined as mean ±1 SD), which is consistent with results from a previous ROMS simulation in the CCS (Hauri et al., 2013b).These results suggest that all habitats studied here have left, or are about to leave, the pCO 2 , pH, and Ar conditions that were experienced during preindustrial times.This is significant as organisms at these sites are now surviving in conditions that are significantly different than the conditions under which their ancestors evolved.
The modeled habitat-specific pCO 2 and Ar conditions for preindustrial, 2012, 2060, and 2100 are shown in Fig. 11.The histograms represent the full range of carbonate conditions at each habitat that was captured by the sensors, which includes both the seasonal and high-frequency variability.The shape of each distribution skews towards more "corrosive" conditions at all sites as the model steps forward into the future.This translates to not only increases in mean pCO 2 but also greater extremes and amount of time spent in extremes.
The projected pCO 2 and Ar envelopes (mean ±SD) at each habitat throughout the modeled period are shown in Fig. 12.An increasing rate of change in pCO 2 ( pCO 2 yr −1 ) is observed, whereas Ar tends to decrease at a relatively constant rate.The rate of increase in pCO 2 is higher than the projected atmospheric CO 2 increase at all subsurface sites.This indicates that as ocean acidification progresses, the effects due to elevated pCO 2 are more likely to become exacerbated with increasing depth.Mean Ar is projected to be < 1 at the shelf break by 2020 and leave the 2012 variability envelope around 2070.

Changes in the buffer factors
The general patterns of the acidification trajectories presented here can be explained by changing buffer factors of seawater, as deeper sites are more strongly influenced by CO 2 -rich upwelled waters.The buffer factors pCO 2 , pH , and CO 3 , are defined as representing the change in each carbonate parameter with respect to a change in DIC (Frankignoulle, 1994).The ef- fect of temperature on is small (< 10 %) between 0 and 15 • C for the DIC and TA values observed here; thus subsequent values were calculated assuming a temperature of 10 • C, TA = 2240 µmol kg −1 , and salinity = 33.5 (Fig. 13).The ability for seawater to buffer changes in pCO 2 diminishes under higher concentrations of DIC.For example, pCO 2 increases from 1.6 to 3.3 at the surface between 2012 and 2100 under the RCP6.0 scenario.However, since deeper waters are naturally elevated in DIC, this effect is more pronounced at the shelf break: pCO 2 increases from 6.2 to 12.3 during the same time interval.This explains why the surf zone had the lowest mean increase in pCO 2 (225 µatm) despite having the highest increase in DIC (82 µmol kg −1 ) out of all of the sites.The shelf break, on the other hand, had the highest increase in pCO 2 (737 µatm) while having the smallest increase in mean DIC (77 µmol kg −1 ) during the same time period.Furthermore, the increase in variability with depth can be explained as well, as the same biological and physical forcings on tidal to seasonal cycles cause a larger change in pCO 2 .
Changes in CO 3 can explain the patterns for Ar , since [Ca 2+ ] and K SP remain unchanged.Unlike pCO 2 , | CO 3 | decreases at higher concentrations of DIC (Fig. 13b); | CO 3 | decreases from 0.62 to 0.57 at the surface, and 0.49 to 0.3 at 88 m between 2012 and 2100.This change in CO 3 explains the decrease in both the rate and range of Ar as anthropogenic CO 2 continues to infiltrate the ocean.
The pH follows a parabolic shape, where there is a maximum decrease in pH per DIC added (Fig. 13c).In a pure carbonate solution, this maximum occurs when DIC = TA, but in seawater it occurs at slightly lower DIC (Frankignoulle, 1994); this maxima occurs at DIC = 2225 µmol kg −1 using the parameters listed above.Therefore we would expect to see a similar trend to pCO 2 for pH in which greatest changes occur at depth relative to the surface as long as the mean DIC is lower than this threshold.This condition is only met at the shelf break (88 m) near the end of the century; thus, as expected, a greater decrease in mean pH and increase in variability (i.e., SD) was observed with depth (Table 2).One exception was observed where the SD decreased as ocean acidification progressed at the shelf break.This is because an increased proportion of time is spent at greater DIC where pH is past its maxima, leading to a smaller variability in pH under the same changes in DIC.It is important to note that the buffer factor of H + ( H+ ) follows a similar pattern to pCO 2 , where it continues to increase as DIC increases (Fig. 13d).Therefore the rate of increase of [H + ] will continue to increase as ocean acidification progresses, and thus biological responses to [H + ] may become exacerbated in the future.

Observed and modeled carbonate chemistry variability
The carbonate conditions presented here are consistent with previous studies.For example, at the shelf break, Ar had a strong seasonal cycle, and undersaturated waters were observed almost continuously throughout the upwelling season (Nam et al., 2015), and remained supersaturated for the rest of the year.This is in good agreement with previous hydrographic surveys in this region, where aragoniteundersaturated waters were observed as shallow as 60 m during the beginning of the upwelling season (Feely et al., 2008) but were not observed in the upper 100 m at the end of the upwelling season in this region (Bednaršek et al., 2014).Furthermore, estimates based on empirical equations showed a similar seasonal pattern in Ar at 88 m, where undersaturated waters were observed every upwelling season (Alin et al., 2012).However, undersaturated waters were not observed in the upper 30 m, unlike northern parts of the CCS where undersaturated conditions are repeatedly observed at the surface during the upwelling season (Bednaršek et al., 2014;Feely et al., 2008;Harris et al., 2013).Due to these traits, the southern portion of the CCS is commonly considered less vulnerable to ocean acidification compared to its northern counterpart.However, our results demonstrate that Ar as low as 1.3 is routinely observed in the kelp forest (17 m), demonstrating the imminent threat of anthropogenic ocean acidification to the southern CCS.
The subsurface habitats characterized in this study routinely experience Ar conditions that have been shown to have non-lethal chronic effects on various bivalve larvae between Ar of 1.2 and 2.0 (Barton et al., 2012;Gaylord et al., 2011;Gazeau et al., 2011;Hettinger et al., 2012;Waldbusser et al., 2015).However, the length of exposure to these unfavorable conditions varies between habitats.For example, the organisms in the kelp forest would be exposed to low-Ar conditions for days to weeks, whereas large tidal variability at the canyon edge could result in periodic exposure to low-Ar conditions on the order of hours.Therefore the effects of low Ar will largely depend on the reproductive timing and environmental variability that occur on event to seasonal timescales; the effects of exposure on various timescales are poorly understood.Such events are expected to become more severe in the future (Hauri et al., 2013a) and thus could lead to an increased rate of failed recruitment of bivalves and other keystone organisms (Byrne et al., 2013).
It may be surprising that the mean diel range of pH was the smallest at the kelp forest (Table 3), as one might expect a large diel cycle driven by photosynthesis and respiration in a highly productive kelp forest.This is most likely because the sensor was deployed near the benthos, below the most productive region of the forest.Frieder et al. (2012) observed significantly larger diel pH variability closer to the surface (7 m depth) compared to near the bottom (17 m depth; same as this study), demonstrating that the biologically driven diel cycle diminishes with increasing depth within the canopy.Therefore it is important to keep in mind that the results presented here are not reflecting kelp forest production dynamics but rather conditions that are experienced by benthic dwelling organisms inside the kelp forest.
The trajectories are sensitive to the choice of the emission scenario (Fig. 14).Trends are similar at all depths; thus only the mean pCO 2 and Ar projections at the shelf break are shown in Fig. 14.The highest emission scenario (RCP8.5)diverges from the two intermediate scenarios around 2030, while the lowest emission scenario (RCP2.6)diverges around 2050.The two intermediate scenarios (RCP4.5 and RCP6.0) do not diverge significantly until 2070.The delayed response to different atmospheric CO 2 trajectories occurs because upwelled waters have spent several decades since they were last in contact with the atmosphere (Feely et al., 2008).Therefore the anthropogenic ocean acidification trajectory for the SCB is already determined for the next several decades, and any mitigation due to changing CO 2 emissions will be delayed.
The results presented here are site-specific, and they do not necessarily reflect conditions at all kelp forests, canyon edges, and shelf breaks.However, if sensor pH data and corresponding regional hydrographic surveys are available, then a DIC anth -σ θ relationship can be established for that region and applied to the sensor data.For example, this approach can potentially expanded to many regions for the CCS, using the North American Carbon Program West Coast Cruise (Feely et al., 2008) and the DIC diseq for the Pacific Ocean  (Sabine et al., 2002).If similar data exist, then this approach can be expanded to other upwelling margins as well.
The SCB experiences a steady but weaker degree of upwelling compared to the northern regions of the CCS, where upwelling events are more pronounced (Bograd et al., 2009).These regions could experience more extreme conditions regularly, as well as significantly higher variability in carbonate conditions (Harris et al., 2013).However, such dynamics are poorly understood, and more high-frequency observations of carbonate parameters along this system are needed.Source water properties must be characterized through hydrographic surveys.Alternatively, for regions where such data for source waters are not available, sensor data can be combined with either global circulation model or ROMS outputs.This approach will alleviate the cost associated with characterizing source waters and, to a large degree, incorporate processes such as interannual variability, decadal changes in source water properties, and reduced ventilation.It is critical that inorganic carbon sensors (e.g., pH or pCO 2 ) are colocated with basic physical oceanographic measurements (e.g., T and S) to determine source water properties especially for subsurface deployments.

Model assessment
The sensitivity of the projected carbonate conditions to the assumptions made in the model is explored here.For example, temperatures observed in 2012 were used to parameterize the model.Sea surface temperature has increased over the past century due to climate change (Smith et al., 2008), and is expected to continue.This will affect the CO 2 equilibrium concentration (DIC eq ), but the effects are small and will reduce DIC eq by only several µmol kg −1 .Both pCO 2 and Ar are dependent on in situ temperature; the effects on Ar are negligible ( / T < 0.01 • C −1 ), whereas pCO 2 / T increases at higher pCO 2 levels, and can be as large as 60 µatm • C −1 at the end of the century, compared to 30 µatm • C −1 at present day at the shelf break.These temperature dependencies will affect the mean conditions, but the magnitude of the variability will be relatively unaffected.However, it should be noted that this simple error analysis does not include any biological feedbacks that increased temperature or CO 2 may induce.For example, phase shifts from kelp-dominated to algal turfs might be an outcome of sea surface warming and acidification (Connell and Russell, 2010), with implications for habitat-scale biogeochemical cycling.Likewise, higher temperatures may increase remineralization rates along the path of the subducted water (Rivkin and Legendre, 2001), further enhancing acidification.TA conditions from 2012 were used to calculate pCO 2 and Ar for all years.Changes in TA affect the buffer factors of seawater; thus, alterations in TA distribution will either speed up or slow down the progression of ocean acidification.However, trends in TA along the CCS on decadal timescales are unknown due to insufficient data.Reduced ventilation in high-latitude seas, altered precipitation patterns, and changes in surface calcification and water-column dissolution rates would all lead to changes in upwelled TA conditions (Fassbender et al., 2011;Lee et al., 2006).Quantifying these processes is difficult and out of the scope of this study.Nevertheless, to demonstrate the magnitude of the uncertainty due to TA, pCO 2 was projected for the year 2100 with a +20 µmol kg −1 bias to TA.The effects were strongly dependent on depth: mean pCO 2 was reduced by approximately 240, 130, and 70 µatm at 88 m, 30 m, and the surface, respectively.
Finally, the model presented here projects future carbonate conditions by assuming that the dynamics that control the variability at each habitat (e.g., seasonal and episodic upwelling events, internal waves and tides, and biological production and respiration) remain the same as 2012 conditions, and it does not account for any variability that occurs on interannual to decadal timescales.For example, changes in O 2 and pH on the continental shelf associated with interannual climate events, such as El Niño, have been observed (Nam et al., 2011).However, since 2012 did not correspond with a strong El Niño or La Niña phase, we believe that it was not strongly biased by such events.Furthermore, recent evidence suggests that the proportion of Pacific Equatorial Waters in the California Undercurrent has been increasing over the past several decades, thus modifying the source water properties for upwelled waters onto the continental shelf (Bograd et al., 2015).Since waters of equatorial origin observed between 100 and 500 m are elevated in DIC and lower in O 2 (Bograd et al., 2015), it is expected that the SCB will experience higher levels of acidification than predicted from this study if this redistribution of water masses of equatorial origin continues.However, at this time, we lack observations with sufficient longevity to predict how climate variability on interannual to decadal timescales might modify the acidification trajectory over the course of the next century.Sustained, high-frequency time series of inorganic carbon parameters are required to elucidate such effects.

Implications for ocean acidification research
In order to properly assess the impacts of anthropogenic ocean acidification through laboratory manipulation experiments, the control and experimental conditions should accurately reflect the study organism's present-day and future habitat conditions (McElhany and Busch, 2013;Reum et al., 2015).The most common control treatment used in ocean acidification experiments for organisms found in the CCS was a pCO 2 value of ∼ 400 µatm, reflecting atmospheric conditions (compiled by Reum et al., 2015).However, our sensor data showed that all subsurface habitats had significantly greater pCO 2 relative to the atmosphere (Table 2).For example, the mean pCO 2 at the kelp forest is about 100 µatm greater than the atmosphere and routinely experiences conditions of more than 300 µatm above atmospheric.Therefore utilizing atmospheric pCO 2 conditions for control treatments will necessarily underestimate the baseline pCO 2 for organisms collected from subsurface habitats.
Recent studies that incorporate natural variability into ocean acidification experiments observed modified responses relative to constant conditions (Dufault et al., 2012;Frieder et al., 2014).However, the effect of natural variability on organismal response to ocean acidification, especially through various life stages is still poorly understood.Our model results demonstrate that variability trajectories are also habitatspecific.For example, in the kelp forest, the variability, approximated by the SD, was 93 µatm in 2012, whereas this increased to 202 µatm in 2100 (Table 2).Furthermore, despite having similar mean CO 2 conditions, the largest variability was observed on event timescales in the kelp forest, whereas the dominant variability occurred on tidal and seasonal cycles at the canyon edge.Therefore future ocean acidification studies investigating the effect of natural variability should not only incorporate increasing magnitude into their experimental design but also consider variability patterns on appropriate timescales.
Temperature and O 2 were tightly correlated with carbonate parameters across habitats and various timescales (daily to seasonal) in this study (Fig. 15); similar correlation has been documented across the CCS in general (Reum et al., 2014(Reum et al., , 2015)).These parameters can potentially act as additional stressors (Padilla-Gamiño et al., 2013) or stress reliefs (Gooding et al., 2009) for ocean acidification.However, laboratory experiments incorporating the effects of temperature (Gooding et al., 2009;Padilla-Gamiño et al., 2013) and O 2 (Frieder et al., 2014) have just started to be explored for the CCS, and no studies have been conducted that incorporate all three variables in their experimental design.Future studies investigating the synergistic effects of O 2 , temperature, and CO 2 should establish experimental conditions based on environmental data (Fig. 15).Although development of systems that can manipulate individual parameters is challenging, important strides have been made to make such experimental setups accessible to the community (Bockmon et al.,  2013).The development of habitat-specific ocean acidification models provides a link between environment and laboratory to facilitate interpretations of physiological responses to elevated CO 2 in the context of current and future environmental conditions.
Discerning habitat-specific CO 2 signatures could lead to the discovery of local populations that are more tolerant of future CO 2 conditions.For example, large high-frequency variability in CO 2 could lead to a greater capacity for physiological and phenotypic plasticity, as organisms are routinely exposed to a wide range of CO 2 .The embryos of Doryteuthis opalescens, an important fishery species in California, can tolerate low pH and O 2 , perhaps due to the fact that they routinely experience a wide range of pH and O 2 (Navarro, 2014).Furthermore, such environmental conditions may be conducive for the existence of highly CO 2 -tolerant subpopulations, allowing for adaptation to buffer some of the negative effects of ocean acidification (Hofmann and Todgham, 2010).Alternatively, these populations could be living near critical biological thresholds, as has been suggested for the thermal of some organisms living in the intertidal (Somero, 2002).A massive failure in an oyster hatchery in Oregon was linked to upwelling of high-CO 2 waters during a critical life stage of oyster larvae (Barton et al., 2012), indicating the existence of CO 2 thresholds for some marine organisms (Bednaršek et al., 2014).However, such thresholds may be dependent on species, life stage, and/or environmental history.As we begin to realize which populations of species and life stages are living near acidification thresholds versus those that exhibit acidification tolerance, implementation of habitat-specific acidification models can be used as a tool to aid protection, management, and remediation efforts of critical marine habitats now and in the future.

Conclusions
Here we have presented habitat-specific carbonate chemistry projections for four coastal habitats along an upwelling margin.The projections were generated by combining highfrequency sensor measurements, a regional empirical rela-tionship for TA, hydrographic survey data to quantify sourcewater properties of upwelled waters, and the atmospheric CO 2 record.Even though the four habitats were within 5 km of one another, distinct habitat-specific variability signatures and acidification trajectories were observed.These results reveal the existence of highly variable CO 2 signatures within a small geographic area, and the potential for discoveries of habitats that could act as refugia from ocean acidification.Changes in the buffer factors largely explained the observed patterns; however, local biological feedbacks could also produce a large acidification signal.In all habitats studied, carbonate conditions have left, or are leaving, preindustrial variability envelopes.Model projections suggest that anthropogenic ocean acidification will continue to progress in the CCS and other upwelling margins over the next several decades regardless of any changes in CO 2 emissions; any impacts from reduced emissions will only be observed midcentury and beyond.This demonstrates the urgency of the situation, and this delayed response must be taken into account when assessing the impacts of ocean acidification and developing mitigation and monitoring strategies.
The Supplement related to this article is available online at doi:10.5194/bg-12-5853-2015-supplement.

Figure 1 .
Figure 1.Map of study region.Hydrographic stations (black dots) and sensor deployment sites (black squares) are shown.Initials are as follows: CE, canyon edge; SB, shelf break; SZ, surf zone; and KF, kelp forest.

Figure 2 .
Figure 2. DIC anth as a function of σ θ for certain modeled years (indicated above line) using the IPCC RCP6.0 projection.The three regimes used in this model -surface, mixing, and subsurface -are labeled.

Figure 4 .
Figure 4. DIC anth as a function of σ θ .The calculated values and the fit are represented by black circles and a red line, respectively.The blue line shows DIC anth using the formulations from Feely et al. (2008).

Figure 5 .
Figure 5.Time series of sensor pH between June 2012 and June.The two lower panels are month-long snapshots.pH is reported at in situ conditions.

Figure 6 .
Figure 6.Time series of pCO 2 calculated from sensor pH and TA est between June 2012 and June 2013.The two lower panels are monthlong snapshots.pCO 2 is reported at in situ conditions.

Figure 7 .
Figure 7. Time series of Ar calculated from sensor pH and TA est between June 2012 and June 2013.The two lower panels are monthlong snapshots.Ar is reported at in situ conditions.

Figure 8 .
Figure 8. Observed pCO 2 in 2012 (black) and modeled pCO 2 using the IPCC RCP 6.0 scenario for the year 2100 (red) at the Del Mar buoy (88 m) over an annual cycle (top).A close-up for the month of December is shown in the bottom panel.Note that the range, but not the absolute values, of the vertical axes for each figure is the same.

Figure 9 .
Figure 9. Observed Ar in 2012 (black) and modeled Ar using the IPCC RCP6.0 scenario for the year 2100 (red) at the shelf break (88 m) over an annual cycle (top).A close-up for the month of December is shown in the bottom panel.Dashed lines represent Ar = 1.Note that the range, but not the absolute values, of the vertical axes for each figure is the same.

Figure 10 .
Figure 10.Observed (black) and modeled (red) pCO 2 (top) and Ar (bottom) at the La Jolla kelp forest (17 m).Modeled values correspond to projected values in 2100 using the IPCC RCP6.0 scenario.Note that the range, but not the absolute values, of the vertical axes for each figure is the same.

Figure 14 .
Figure 14.Projections of mean pCO 2 (left) and Ar (right) at the shelf break (88 m) based on four projections from the Fifth Assessment of the IPCC.

Figure 15 .
Figure 15.pCO 2 as a function of O 2 (left) and temperature (right) from the kelp forest and shelf break.Data observed in 2012 and projected for 2100 are plotted.

Table 2 .
Mean ±SD of modeled carbonate parameters at in situ conditions for preindustrial, 2012, 2060, and 2100 using the RCP6.0 projection at each habitat.

Table 3 .
Mean daily range of carbonate parameters at in situ conditions for 2012 at each habitat.