On the role of climate modes in modulating the air-sea CO<sub>2</sub> fluxes in Eastern Boundary Upwelling Systems

Abstract. The air-sea CO2 fluxes in Eastern Boundary Upwelling Systems (EBUS) vary strongly in time and space with some of the highest flux densities globally. The processes controlling this variability have not yet been investigated consistently across all four major EBUS, i.e., the California (CalCS), Humboldt (HumCS), Canary (CanCS), and Benguela (BenCS) Current Systems. In this study, we diagnose the physical and biological mechanisms that contribute to historical (1920–2015) CO2 flux variability in these regions using simulation results from the Community Earth System Model Large Ensemble (CESM-LENS), a global coupled climate model ensemble that is forced by historical and RCP8.5 radiative forcing. Differences between simulations can be attributed entirely to internal climate variability. We find that the deviations from the ensemble mean, i.e., the anomalous CO2 fluxes, in the CalCS and CanCS are strongly affected by modes of variability associated with atmospheric subtropical gyres: the North Pacific Gyre Oscillation (NPGO) and the North Atlantic Oscillation (NAO), respectively. The CalCS (CanCS) has anomalous uptake (outgassing) of CO2 during the positive phase of the NPGO (NAO). The HumCS is mainly affected by El Niño Southern Oscillation (ENSO), with anomalous uptake of CO2 during an El Niño event. Variations in dissolved inorganic carbon (DIC) and sea surface temperature (SST) are the major contributors to these anomalous CO2 fluxes, and are generally driven by changes to gyre circulation, upwelling, the mixed layer depth, and biological processes. A better understanding of the sensitivity of EBUS CO2 fluxes to modes of climate variability may improve our ability to predict the ocean–atmosphere carbon cycle in EBUS, which are particularly susceptible to ocean acidification.



Introduction
The four major Eastern Boundary Upwelling Systems (EBUS) occur at the eastern edges of subtropical gyres in the Atlantic and Pacific oceans -the California (CalCS), Humboldt (HumCS), Canary (CanCS), and Benguela (BenCS) Current Systems. These regions are characterized by seasonal or permanent equatorward winds that cause upwelling due to both offshore Ekman transport as well as wind stress curl-driven Ekman suction within the first 200 km of the coastline (Chavez and Messié, 2009). 20 Upwelling delivers deep waters with respired nutrients to the surface, fueling primary production and ultimately supporting fisheries that are highly productive with respect to the small surface area they cover (Ryther, 1969). Upwelled waters also However, single realizations provide a limited sample size of internal variability and confound the impacts of external forcing (e.g., land use change, fossil fuel emissions, volcanic eruptions) with internal climate variability. The former is problematic, because ENSO events evolve differently in the tropics (Capotondi et al., 2015) and mid-latitude atmospheric noise can obscure the tropical-extratropical connections associated with climate modes such as ENSO, causing a diversity of responses in EBUS (Deser et al., 2017(Deser et al., , 2018. The latter makes it difficult to isolate the internal component of variability in CO 2 fluxes from the 25 seasonal cycle and anthropogenic and other external forcing. One solution to this problem is to use a single-model ensemble that is derived by introducing perturbations to the initial state of the climate system. This gives rise to a set of realizations with unique representations of internal climate variability and gives one access to many hundred ENSO events. By performing historical experiments with increasing atmospheric CO 2 rather than a long control simulation, we can account for variability in air-sea flux of both natural CO 2 (the component of ocean CO 2 in equilibrium with pre-industrial atmospheric CO 2 ) and 30 anthropogenic CO 2 .
In this study, we utilize output from the single-model Community Earth System Model "Large Ensemble" (CESM-LENS; Kay et al., 2015;Lovenduski et al., 2016) to identify major modes of climate variability that are associated with anomalous CO 2 fluxes in the major EBUS. We expand on this by investigating the physical and biological drivers that underpin these anomalies. The single-model ensemble allows for an assessment of both the forced signal and the CO 2 flux response to internal 35 3 Biogeosciences Discuss., https://doi.org /10.5194/bg-2018-415 Manuscript under review for journal Biogeosciences Discussion started: 21 September 2018 c Author(s) 2018. CC BY 4.0 License. climate variability. Since the simulations are forced with historical CO 2 emissions, each member accounts for variability in both natural and anthropogenic CO 2 . Furthermore, the availability of 34 simulations allows us to find statistically robust relationships between anomalous CO 2 fluxes and internal climate variability. 5 We utilize monthly mean output from 34 members of the CESM-LENS (Kay et al., 2015;Lovenduski et al., 2016), which is derived from the Community Earth System Model, version 1, with the Community Atmosphere Model, version 5 (CESM1(CAM5); Hurrell et al., 2013). Along with standard atmosphere, ocean, land, and sea ice components, the CESM-LENS simulations include land and ocean biogeochemistry. The ocean biogeochemical component of CESM1(CAM5) is the Biogeochemical Elemental Cycling (BEC) model, which has three phytoplankton functional groups and tracks the cycling of C, N, P, Fe, Si, 10 and O in the ocean. Further information on the implementation of BEC in CESM1 can be found in Moore et al. (2013) and Lindsay et al. (2014). The ocean component is the Parallel Ocean Program, version 2 (Smith et al., 2010) and has a nominal 1 o horizontal resolution with vertical resolution of 10 m through the upper 250 m, thereby resolving the Ekman layer. Due to the coarse horizontal resolution, neither curl-driven nor coastal upwelling is directly resolved, but both are represented in the model. A more detailed discussion of coastal upwelling in the CESM-LENS for the CalCS in particular can be found in (Brady 15 et al., 2017).

Model Configuration
The 34 ensemble members of CESM-LENS were generated using round-off level (order 10 −14 K) level perturbations in the initial atmospheric temperature. (Kay et al., 2015). This generates an ensemble of simulations that diverge solely due to the influence of internal variability. Each ensemble member was forced with a common external forcing: historical radiative forcing from 1920 to 2005 and RCP8.5 radiative forcing from 2006 to 2100. In contrast to a pre-industrial control simulation, 20 CESM-LENS includes the addition of anthropogenic carbon to air-sea CO 2 fluxes. Anthropogenic CO 2 is explicitly separated from natural CO 2 by computing air-sea CO 2 fluxes relative to a constant 284.7 ppm pre-industrial atmosphere and subtracting it from CO 2 fluxes responding to the evolving atmosphere under historical and RCP8.5 radiative forcing.

Upwelling Regions
Upwelling regions in each EBUS span approximately the 10 o latitude of most active upwelling as defined by Chavez and Messié 25 (2009); we shifted the CanCS upwelling domain north by 9 o to capture the more intense upwelling off the Western Sahara in CESM1 (Table 1). Following Turi et al. (2014), the EBUS upwelling regions span from the coastline to 800 km offshore. The black outlines in Figure 1e-h display these regions. The CESM-LENS ensemble mean incorporates the seasonal cycle and any long-term anthropogenic trends for a given variable. We therefore removed the ensemble mean from each simulation to create a time series of anomalies produced solely by internal climate variability. Note that for CO 2 flux, these anomalies represent 30 internal variability in the contemporary air-sea flux of CO 2 , or the combined variability of natural and anthropogenic CO 2 .

Climate Indices
Throughout this paper, we regress anomaly time series of variables from each EBUS onto climate indices. Climate index time series were available for each ensemble member for Nino3, the PDO, NAO, Atlantic Multidecadal Oscillation (AMO), and Atlantic Meridional Overturning Circulation (AMOC) through the Climate Variability Diagnostics Package (CVDP;Phillips et al., 2014). We computed the NPGO index for each simulation following Di Lorenzo and Mantua (2016), and these indices 5 are available through the CESM-LENS project page on the NPGO.

Statistical Analysis and Model Equations
Air-sea CO 2 fluxes in CESM are computed following the parameterization of Wanninkhof (2014): where k represents the gas transfer velocity (dependent on the wind speed squared), K 0 the solubility of CO 2 in seawater, and 10 pCO o 2 and pCO a 2 the partial pressures of CO 2 in the surface ocean and atmosphere, respectively. We use a linear Taylor expansion to quantify the relative contribution of each variable to the overall CO 2 flux anomaly in response to internally generated variability following Lovenduski et al. (2007) and Turi et al. (2014), where ∂F ∂U and ∂F ∂pCO oc 2 are determined from the model equations and mean values in each EBUS. ∆'s represent the linear 15 regression of the given variable's anomalies onto a climate index. The contributions from ∆pCO oc 2 is further decomposed into DIC, Alk, SST, and salinity terms: Because DIC and Alk can be diluted by freshwater fluxes, we introduce salinity-normalized DIC (sDIC) and Alk (sAlk), 20 Due to the significance of sDIC anomaly contributions to the total CO 2 flux anomaly in EBUS, we approximate the mechanisms controlling sDIC anomalies following Lovenduski et al. (2007), where J circ , J bio , and J ex represent the sources and sinks of sDIC from circulation, biology, and CO 2 flux anomalies integrated over the upper 100m, respectively. See Lovenduski et al. (2007, their Equation 4) for additional details on these terms and their 25 Appendix B for the computation of J bio in particular.
To compensate for autocorrelation that is characteristic of climate indices and is also introduced from smoothing, we replace the t-statistic sample size N with an effective sample size, N ef f : where r 1 and r 2 are the lag-1 autocorrelation coefficients of the two time series being correlated (Bretherton et al., 1999;Lovenduski and Gruber, 2005). N ef f represents the number of statistically independent measurements.

Model Evaluation
CESM-LENS air-sea CO 2 flux climatologies and pCO 2 seasonal cycles were compared to the SOM-FFN (

CO 2 Flux Climatology
The mean state of Pacific EBUS CO 2 fluxes is particularly well-modeled in CESM-LENS ( Figure 1). The CESM-LENS 10 captures the meridional gradient of poleward uptake and equatorward outgassing of CO 2 in the CalCS (Figure 1e). In this and the other EBUS, the coarse resolution model is not capable of resolving the narrow band of nearshore CO 2 outgassing found in high resolution model solutions (e.g., Turi et al., 2014;Fiechter et al., 2014). In the HumCS, the model depicts the strong outgassing that is characteristic of a tropical upwelling system ( Figure 1f)  comm. with RJ Small, 2018). Further, the BenCS only improves to a 5 o C bias at 0.5 o atmospheric resolution (Gent et al., 2010).
This bias is likely driven by the fact that the Angola-Benguela Front is simulated too far south, in addition to deficiencies in upwelling and meridional transport that are caused by unrealistic alongshore wind stress structure (Small et al., 2015). Because these deficiencies are specific to the BenCS, we will only discuss its representation of the pCO 2 seasonal cycle in Section 3.2, and its internal variability in CO 2 fluxes in Section 4.1, but will not perform a full analysis on its connections to larger-scale 30 climate variability.

pCO 2 Seasonal Cycle
CESM-LENS simulates the pCO 2 seasonal cycle for the EBUS with similar accuracy to its depiction of the mean state -the Pacific systems are generally well-modeled, while larger biases exist in the Atlantic regions. Beginning with the CalCS, CESM-LENS nearly perfectly matches the SOM-FFN in both amplitude and phase ( Figure 2a). The system exhibits its maximum pCO 2 (and thus CO 2 outgassing) in August, and its minimum pCO 2 (and thus CO 2 uptake) in April. We further decomposed 5 the model seasonal cycle into its thermal component (driven by the seasonality of SST) and its non-thermal component (driven by the seasonality of factors such as DIC, ALK, and salinity) following Takahashi et al. (2002). This decomposition suggests that the phase of the CalCS seasonal cycle is determined by thermal (solubility) effects, with its amplitude modulated by nonthermal factors ( Figure 2a). These non-thermal factors are almost entirely driven by the seasonal cycle of DIC (not shown), which is characterized by photosynthetic uptake of CO 2 in the summer and fall, coinciding with upwelling season. These

Internal Variability in Upwelling Systems
We emphasize the magnitude of internal variability in EBUS CO 2 fluxes in Figure Table 1). The BenCS is uniquely exposed to variability from the Southern Ocean and Agulhas Current (Reason et al., 2006). The HumCS likely has intense variability due to its proximity to the tropical Pacific Ocean and thus rapid communication with ENSO (e.g., Colas et al., 2008;Montes et al., 2011).
All four systems have statistically significant trends toward a weaker CO 2 source or greater CO 2 sink over 1920-2015 due 15 mainly to the invasion of anthropogenic carbon ( Figure 4; Table 1). Note that in the CanCS and BenCS, there is a trend toward more outgassing of natural CO 2 , which is compensated for by the relatively large uptake tendency of anthropogenic CO 2 (Table 1). The long-term trend forces the HumCS, CanCS, and BenCS to act as intermittent seasonal sinks by 2015 in some realizations due to the combination of the long-term trend and internal variability ( Figure 4). The BenCS and CanCS have the largest uptake of anthropogenic CO 2 over the historical period, although the HumCS has a relatively large uptake of both 20 natural and anthropogenic CO 2 (Table 1). For the HumCS and BenCS, the contemporary CO 2 flux trend is on the order of the magnitude of their seasonal cycles over the course of 96 years. The CanCS is a unique case, where the contemporary trend is more than double the magnitude of its seasonal cycle (Table 1).
The magnitude of internal variability is greater than that of the seasonal cycle for the majority of systems. The non-seasonal component of the total variability (the sum of the seasonal and internal components) is 59% for the HumCS, 73% for the 25 CanCS, and 56% for the BenCS (Table 1). Only the CalCS has a stronger seasonal cycle of CO 2 flux than internal variability, but the non-seasonal component still accounts for 33% of the total variability in this system (Table 1). Perhaps for the CalCS, as well as the other EBUS, more significant internal variability would be captured at a higher resolution that resolves coastal upwelling, such as in Turi et al. (2014, their Figure 8c). Lastly, internal variability in CO 2 fluxes tends to be phase-locked with the seasonal cycle, as the peak magnitudes of internal variability track the ridges and troughs of the seasonal component

Climate Variability and CO 2 Fluxes
Our primary goal for each EBUS was to identify the mode of climate variability most strongly associated with its CO 2 flux anomalies. We correlated area-weighted anomalies from the black boxes in Figure 1e-h for each simulation with every grid cell globally for a set of predictor variable anomalies: SST, sea level pressure (SLP), 10m wind speed, and wind stress curl.
We then assessed the ensemble mean of the correlations to determine the mode of climate variability associated with the given 5 global spatial pattern.

California Current
The correlation between CalCS CO 2 flux anomalies and SSTa yields a map suggestive of Pacific Decadal Variability, due to the  Figure 5d). Thus, we highlight the NPGO as the major mode of climate variability associated with anomalous CO 2 flux in the CalCS. 15 We find that the CalCS has a single-signed response to the NPGO with anomalous uptake of CO 2 during a positive event, intensifying the mean state of the system as an uptake site ( Figure 6). The direct regression of ∆F onto the NPGO results in an anomalous uptake of 0.10 mol m −2 yr −1 (Table 2), which is roughly 24% of the long-term historical flux mean of -0.42 mol m −2 yr −1 . The primary contributions to this uptake anomaly come from variations in SST and sDIC, which are mainly driven by changes to offshore gyre dynamics. As the oceanic expression of the North Pacific Oscillation, the NPGO is 20 associated with intensified geostrophic circulation, resulting in increased transport in both the California and Alaskan Coastal Currents (Di Lorenzo et al., 2008). Further, a positive NPGO leads to enhanced upwelling-favorable winds south of Cape Mendocino ( Figure S1). During a positive NPGO event, the entire CalCS cools, increasing the CO 2 solubility of the system ( Figure 6, Figure S1). Although downwelling is enhanced offshore, increased DIC from the Alaskan Gyre leads to a tendency for outgassing of CO 2 ( Figure S1). Nearshore south of Cape Mendocino, increased upwelling of DIC-enriched waters is 25 compensated for by photosynthesis, leading to near-zero CO 2 flux anomalies ( Figure S2). North of Cape Mendocino, weakened upwelling and low subsurface DIC anomalies lead to enhanced uptake anomalies ( Figures S1, S2). Because the system-wide contributions of SST and sDIC to the anomalous flux nearly balance each other, minor contributions from wind, salinity, sAlk, and freshwater flux push the system in favor of anomalous uptake ( Figure 6). The CalCS has the largest relative ensemble spread in sDIC and sAlk ( Figure 6; Table 2). This is potentially because of inter-simulation variability in the response of CalCS 30 dynamics to the NPGO due to atmospheric noise (as in the case of ENSO in Deser et al., 2017Deser et al., , 2018 which can directly alter the biogeochemical properties of source waters that feed the region (Pozo Buil and Di Lorenzo, 2017). Although the linear Taylor expansion approximates a CO 2 flux anomaly nearly half that of the direct regression of ∆F onto the NPGO, it is still of the same sign. This discrepancy is due to the influence of higher-order and cross-derivative terms that we did not account for in our linear approximation.
We also performed this analysis for the CalCS response to a 1σ positive (warm) phase of the PDO (Figure 7a and b). Every ensemble member displayed a dipole response to the PDO (not shown), with anomalous uptake in the nearshore region south of Cape Mendocino, and anomalous outgassing elsewhere in the domain (Figure 7c). This was the only case in which we found  Figure S3). In coordination with reduced subsurface DIC throughout the region, these changes in circulation contribute toward 15 anomalous uptake of CO 2 throughout the CalCS. Note that the nearshore decomposition in Figure 7a has a y-axis range four times smaller than that of the offshore decomposition. This slight uptake anomaly is the result of a delicate balance of minor terms, where the sDIC reduction slightly outweighs the warming effect. On the other hand, the offshore region has contributions from SST and sDIC that are as much as triple the magnitude as that for the NPGO (Table 2). Despite the sDIC reduction being larger than the SST term, the reduced sAlk is substantial enough to cause a slight outgassing anomaly offshore (Figure 7b).

20
The direct response of winds to the NPGO and PDO plays a negligible role in influencing anomalous CO 2 flux in the CalCS (Table 2). Although ∆U in response to the NPGO and PDO is on the order of the HumCS and CanCS, ∂F ∂U is 3-10 times smaller than the other systems. ∂F ∂U is based on the climatological mean U, ∆pCO 2 , and Schmidt number. The CalCS has the smallest mean ∆pCO 2 of the EBUS -just 0.2µatm. This causes CO 2 flux in the system to be relatively insensitive to fluctuations in the wind.

Humboldt Current
Global correlations between the HumCS CO 2 flux anomalies and SSTa display ENSO as a major influencer, with regions of high correlation focused around the equatorial Pacific. Correlations between HumCS CO 2 flux anomalies and the Nino3 index resulted in an r-value of -0.40 ± 0.04 (Figure 5e). Similar results were found for the Nino3.4 index (-0.38 ± 0.04) and the Nino4 index (-0.36 ± 0.05). We chose the Nino3 index as our primary predictor of HumCS CO 2 flux anomalies, since it is 30 more eastern-focused and thus captures the stronger spatial correlations closest to the HumCS (Figure 5b).
We present the results of a linear Taylor expansion for HumCS CO 2 flux anomalies regressed onto a 1 o El Niño in Figure 8 (Eq. 4). We find that the HumCS responds with a nearly single-signed CO 2 uptake anomaly, resulting in a weakening of the climatological outgassing (Figure 8b). Although there is a small region in the northern HumCS that responds with an Biogeosciences Discuss., https://doi.org /10.5194/bg-2018-415 Manuscript under review for journal Biogeosciences Discussion started: 21 September 2018 c Author(s) 2018. CC BY 4.0 License. outgassing anomaly, it is nowhere near as coherent across the ensemble as was the spatial dipole response of the CalCS to the PDO (Figure 7c). The direct regression of ∆F onto the Nino3 index results in an anomalous uptake of 0.49 mol m −2 yr −1 , which is approximately 18% of the long-term historical mean of 2.8 mol m −2 yr −1 . As in the case of the CalCS, the two major terms contributing to the uptake anomaly are sDIC and SST, which are in opposition to one another (Figure 8a). We would anticipate this to be the case, as an El Niño event induces warming along the HumCS as well as reduces the efficacy of 5 upwelling for bringing carbon-and nutrient-rich waters to the surface due to the presence of an anomalously deep thermocline ( Figure S4; Strub et al., 1998). While upwelling-favorable winds tend to decrease along Chile (outside of our study region) during an El Niño event, the coastal wind response along Peru is variable, ranging anywhere from slight downwelling-favorable to slight upwelling-favorable anomalies (Wyrtki, 1975;Enfield, 1981;Huyer et al., 1987).
In CESM-LENS, the HumCS experiences a warming of 0.7 o C for a 1 o El Niño, which results in an outgassing tendency of  (Table 2). However, sDIC in the system is reduced by 13.2 mmol m −3 for the same event, which translates to a large uptake contribution of 0.8 mol m −2 yr −1 (Table 2). This is an enormous change in sDIC, which is partially driven by the high subsurface DIC bias in the east equatorial Pacific in CESM1 (see Lovenduski et al., 2015, their Figure 2). The large sDIC reduction is due to weakened upwelling and a deepening of the thermocline by advected warm waters from the equatorial Pacific ( Figure S4). Although the passage of coastally trapped waves are important in the HumCS, they are not resolved in 15 the coarse resolution CESM-LENS. Lastly, there is a minor outgassing anomaly of 0.06 mol m −2 yr −1 in response to a slight intensification of wind magnitude during El Niño (Table 2). Despite the significant contributions of wind speed, SST, and sAlk toward outgassing, the large reduction in sDIC drives an uptake anomaly that weakens the HumCS outgassing during an El Niño event. Note that these are the same mechanisms that are responsible for reduced outgassing in the tropical belt in response to an El Niño .

Canary Current
Correlations between the CanCS CO 2 flux anomalies and SLPa globally reveal a region of high positive r-values just northwest of Africa. This encircles the climatological position of the Azores High, the atmospheric subtropical gyre which forces the CanCS. The climate index that most directly captures variability in the Azores High is the NAO, and will thus be considered the main mode of climate variability that modulates anomalous CO 2 flux in the CanCS. We find modest correlations of 0.28 25 ± 0.03 between annually smoothed CanCS CO 2 flux anomalies and the NAO (Figure 5f). Relatively lower correlations are expected between anomalous EBUS CO 2 fluxes and atmospheric indices, as the atmosphere is noisier than the more slowly evolving ocean.
Grid cell correlations between CanCS CO 2 flux anomalies and the NAO are displayed in Figure 9b. The CanCS has a nearly single-signed response of increased outgassing during the positive phase of the NAO. The direct regression of ∆F phase of the NAO, a stronger Azores High leads to intensified alongshore winds and thus more vigorous upwelling ( Figure   S5). This brings up additional deep cold water which in turn increases the CO 2 solubility of the system, tending toward an uptake anomaly of 0.15 mol m −2 yr −1 (Table 2). On the other hand, the increased sDIC from intensified upwelling is double the magnitude of the SST contribution, despite the increased biological activity which reduces some of the physical sDIC input, leading to an outgassing anomaly of 0.33 mol m −2 yr −1 . This large sDIC response is driven both by a high ∆sDIC 5 of 3.9 mmol m −3 per 1σ NAO as well as the fact that the CanCS has the highest ∂F ∂sDIC of the major EBUS. Increased winds of 0.3 m s −1 per 1σ NAO lead to a significant outgassing pressure of 0.05 mol m −2 yr −1 . This is due both to a large system sensitivity, ∂F ∂U , to changes in wind and a high wind anomaly in response to the NAO. Ultimately, intensified winds and an anomalous increase in sDIC due to enhanced upwelling counteracts the solubility effects of colder SSTs and increased photosynthetic uptake of sDIC. This leads to the highest relative CO 2 flux anomaly of any system, with a 21% increase in 10 outgassing per 1σ NAO event relative to the long-term mean.

Conclusions and Discussion
We find that the seasonal cycle of pCO 2 is single-peaked and its phase is driven by thermal (solubility) effects in the CalCS, CanCS, and BenCS, while non-thermal effects modulate its amplitude. The seasonal cycle of pCO 2 in the HumCS is dualpeaked and driven by an interchanging importance between thermal and non-thermal factors. Although CESM-LENS models 15 the CalCS and HumCS pCO 2 seasonal cycle well, the CanCS seasonal cycle is damped and the BenCS seasonal cycle entirely out of phase. Variations in sDIC and SST in response to large modes of climate variability exert the most influence on anomalous CO 2 fluxes in the CalCS, CanCS, and HumCS. Further, these terms always act in opposition to one another in the portions of variability associated with the climate indices explored in this study. Secondary to these terms are wind speed and sAlk.
Although their contributions do not rival those of SST and sDIC in magnitude, they act to further reinforce anomalies or to tip 20 the balance toward outgassing or uptake when sDIC and SST are of equal magnitude. In all systems, salinity and freshwater fluxes have negligible contributions toward the total CO 2 flux anomaly. CalCS and CanCS CO 2 flux anomalies are associated mainly with climate modes related to the subtropical anticyclonic gyres, with their mean states (uptake for the CalCS and outgassing for the CanCS) intensified during phases of enhanced gyre circulation in the NPGO and NAO, respectively. On the other hand, CO 2 flux anomalies in the HumCS are mostly driven by ENSO, due to its close proximity to the equatorial Pacific. 25 We find that the HumCS has weakened outgassing during El Niño due to reduced upwelling and a deepening of the thermocline, which are similar mechanisms to the equatorial Pacific's response to El Niño . The major EBUS are often lumped together in studies due to their similarities -they are all characterized by their presence on the eastern flank of subtropical gyres, their Ekman dynamics associated with this positioning which leads to coastal and curl-driven upwelling, their productive fisheries, and in the case of our study, their high variability in CO 2 fluxes. However, we show in this study that 30 their position in terms of latitude and ocean basin as well as their coastal geometry leads to unique physical and biogeochemical responses to climate variability. In particular, despite variations in sDIC being a leading contributor to CO 2 flux anomalies, the drivers of these sDIC anomalies differ between EBUS. Biogeosciences Discuss., https://doi.org /10.5194/bg-2018-415 Manuscript under review for journal Biogeosciences Discussion started: 21 September 2018 c Author(s) 2018. CC BY 4.0 License.
Because sDIC anomalies are so critical to the EBUS response to internal climate variability, we further decomposed these anomalies into their absolute (Table 3) and relative ( Figure 10) contributions from circulation (e.g., upwelling, advection, diffusion), biology (photosynthesis, respiration, and calcium carbonate formation and dissolution), and CO 2 fluxes. Changes in biological activity in the HumCS in response to ENSO contribute to nearly 50% (8.611 TgC yr −1 ) of the total sDIC tendency anomaly ( Figure 10, Table 3). During an El Niño (La Niña) enhanced respiration (photosynthesis) pumps sDIC into (out of) the 5 upper water column. While biology is the major contributor to sDIC tendency anomalies in the HumCS, circulation changes play a leading role in the CalCS. In response to the NPGO, circulation anomalies in the CalCS contribute to roughly 47% (1.467 TgC yr −1 ) of the total sDIC tendency anomaly ( Figure 10, Table 3). It is difficult to assess the relative contributions of J circ and J bio in the CanCS, as the large ensemble spread in J bio drives a highly uncertain J circ , which is computed as the anomaly of the other three terms ( Figure 10). In all three systems, CO 2 flux anomalies play an important role in modifying 10 the sDIC tendency, with an ensemble mean relative contribution greater than 25% ( Figure 10). Note, however, that significant (moderate) spatial variability exists for circulation and biology contributions in the HumCS (CanCS), while contributions are relatively uniform throughout the CalCS ( Figure S7).
It is important to note that anthropogenic climate change will likely modify our findings over the coming decades in a number of ways. The long-term addition of anthropogenic carbon to the surface ocean causes EBUS to become greater sinks for CO 2 . 15 This trend shifts the mean state of the EBUS, causing historical outgassing sites (the HumCS, CanCS, and BenCS) to become intermittent net sinks of CO 2 under the influence of internal variability by the end of the historical period (Figure 4). Projecting to 2100 under RCP8.5 forcing, we find that the CanCS and BenCS become net sinks for CO 2 due to the anthropogenic trend, with mean values of -0.43 mol m −2 yr −1 and -0.04 mol m −2 yr −1 over 2016-2100, respectively. The addition of anthropogenic CO 2 into the surface ocean can also intensify the magnitude of the seasonal cycle of CO 2 flux. This is due to the increased 20 concentration of CO 2 in the surface ocean as well as the reduction in the ocean's buffer capacity, which makes pCO 2 more sensitive to seasonal fluctuations in DIC and alkalinity (Landschützer et al., 2018, and references therein). This effect has been shown in observations (Landschützer et al., 2018) and is also projected in climate models (Kwiatkowski and Orr, 2018), but is only seen significantly in the CalCS and CanCS, with an approximate 37% and 30% increase in the CO 2 flux seasonal cycle over 2016-2100, respectively. Negligible changes to the seasonal cycle occur in the HumCS and BenCS.

25
External forcing will also alter the dynamics of the EBUS (Bakun et al., 2015;, potentially inducing changes to alongshore winds (Narayan et al., 2010;Sydeman et al., 2014;Oerder et al., 2015;Wang et al., 2015) and upper ocean stratification (Di Lorenzo et al., 2005;Rykaczewski and Dunne, 2010;Oerder et al., 2015), which will in turn influence the rate and efficacy of upwelling. It is unlikely that changes to vertical transport will be detectable until mid-century, as the anthropogenic signal is obscured by internal variability (Brady et al., 2017). However, 30 externally forced trends in stratification are likely to emerge much quicker (Henson et al., 2016). The biogeochemical signature of waters feeding upwelling (e.g., O 2 , CO 2 , and nutrient concentrations) will likely also be modified due to these dynamical changes as well as to changes to ocean ventilation and source water pathways (Rykaczewski and Dunne, 2010). Externally forced changes in physical and biogeochemical properties of EBUS will likely alter the relative contributions of variables to anomalous CO 2 fluxes, as approximated by Equation 4 and shown in Figures 6-9 and Table 2. It is also possible that modes of climate variability will change in response to anthropogenic forcing. Model projections and long-term observations have suggested that the intensity, frequency, or variance associated with ENSO (e.g., Timmermann et al., 1999;Cai et al., 2014Cai et al., , 2015, the NAO (Kuzmina et al., 2005), and the NPGO (Sydeman et al., 2013) has changed significantly in recent decades or will change over the next century. While not observed in our historical modeling study, modifications to modes of climate variability associated with the major EBUS could directly influence the magnitude of internally generated anomalies in CO 2 5 fluxes in the future. Further investigation of the impacts of anthropogenic climate change on CO 2 fluxes in EBUS is necessary for the community.
This study serves as a starting point toward better understanding how internal climate variability modulates CO 2 fluxes in the major EBUS. Here, we only present the leading mode of climate variability associated with anomalous CO 2 fluxes in the HumCS and CanCS, and the leading two modes for the CalCS. Because of this, we explain approximately 16% of the 10 total CO 2 flux variance in the HumCS and 8% in the CanCS. Since we analyzed statistically orthogonal modes (the PDO and NPGO) in the CalCS, we were able to explain as much as 31% of the total CO 2 flux variance. It is difficult to explain a large chunk of the remaining variance for the EBUS, as other modes of internal variability are physically dependent on one another.
Further, locally generated atmospheric noise in the EBUS can contribute substantially to the total modeled CO 2 flux variance.
Previous studies that relate anomalous CO 2 fluxes to internal climate variability have explained similar amounts of variance. Our study is limited by our use of a coarse single model ensemble, which carries uncertainty due to structural biases in the representation of the climate system and biogeochemistry, the ensemble's ability to accurately simulate the magnitude 20 and frequency of internal climate variability, and processes that occur at a finer scale than the grid resolution. The community would benefit from future studies involving multiple single-model ensembles, which would reduce uncertainty due to structural biases, such as in the dynamics of the BenCS and the elevated sub-surface DIC concentration in the east equatorial Pacific.
Due to model resolution, we do not directly resolve the coastal upwelling process which induces vigorous outgassing within the first O(10km) of the coastline. This problem could be mitigated by nesting high-resolution EBUS simulations within a 25 coarser global ensemble or by using regional mesh refinement techniques. This would allow the remote propagation of climate variability into the EBUS, while avoiding the high computational cost of running multiple high-resolution global simulations.
In particular, the BenCS requires significant attention. We find pronounced internal variability in CO 2 fluxes in the BenCS in CESM-LENS that warrants investigation in a high-resolution model specific to the BenCS. We anticipate that these results and further investigation of the relationship between internal climate variability and anomalous CO 2 fluxes in EBUS will be 30 useful for the rapidly developing subseasonal to decadal prediction community. Skillful prediction of climate variability, such as ENSO, the NPGO, and NAO, could be linked directly to anomalous fluxes of CO 2 in the major EBUS. As these systems are naturally sensitive to the undersaturation of calcium carbonate, these predictions could also aid in detecting and managing the onset of seasonal ocean acidification.        in the CESM-LENS for each of the four studied upwelling systems (a-d). The ensemble mean yields both the seasonal cycle (black) and the forced trend (red). Gray shading shows the bounds of the maximum and minimum realizations due to internal variability, but also contains signal from the seasonal cycle and forced trend. Table 1 displays Nea rsho re Offs hore Figure 7. As in Figure 6, but in response to the PDO for a nearshore region (a) and offshore region (b). Note that the offshore decomposition (b) has a y-axis range four times that of the nearshore decomposition (a).

29
Biogeosciences Discuss., https://doi.org /10.5194/bg-2018-415 Manuscript under review for journal Biogeosciences    Figure 10. Relative contributions of anomalous CO2 flux, physical circulation, and biology to the approximated sDIC tendency anomaly for the CalCS (yellow), HumCS (orange), and CanCS (red) in response to the NPGO, Nino3, and NAO, respectively. Error bars show the 1 standard deviation spread of the ensemble. For a given simulation, the absolute values of the three components sum to 100%. These values were approximated using Equation 5. Table 1. Statistics for CO2 fluxes in the CalCS, HumCS, CanCS, and BenCS from 1920-2015. The seasonal component is computed as the standard deviation of the ensemble monthly climatology after removing a fourth-order polynomial fit to remove the long-term trend.
The internal component is computed as the ensemble mean standard deviation of the anomalies. The trend is computed as the first-order ordinary least squares regression coefficient for the contemporary (total), anthropogenic, and natural CO2 fluxes. The intercept is derived from this linear regression. The non-seasonal component represents the fraction of total variability (seasonal plus internal) to which the internal component contributes.