A monthly surface pCO₂ product for the California Current Large Marine Ecosystem

A common strategy for calculating the direction and rate of carbon dioxide gas (CO₂) exchange between the ocean and atmosphere relies on knowledge of the partial pressure of CO₂ in surface seawater (pCO_2(sw)), a quantity that is frequently observed by autonomous sensors on ships and moored buoys, albeit with significant spatial and temporal gaps. Here we present a monthly gridded data product of pCO_2(sw) at 0.25^∘ latitude by 0.25^∘ longitude resolution in the northeastern Pacific Ocean, centered on the California Current System (CCS) and spanning all months from January 1998 to December 2020. The data product (RFR-CCS; Sharp et al., 2022; https://doi.org/10.5281/zenodo.5523389) was created using observations from the most recent (2021) version of the Surface Ocean CO₂ Atlas (Bakker et al., 2016). These observations were fit against a variety of collocated and contemporaneous satellite- and model-derived surface variables using a random forest regression (RFR) model. We validate RFR-CCS in multiple ways, including direct comparisons with observations from sensors on moored buoys, and find that the data product effectively captures seasonal pCO_2(sw) cycles at nearshore sites. This result is notable because global gridded pCO_2(sw) products do not capture local variability effectively in this region, suggesting that RFR-CCS is a better option than regional extractions from global products to represent pCO_2(sw) in the CCS over the last 2 decades. Lessons learned from the construction of RFR-CCS provide insight into how global pCO_2(sw) products could effectively characterize seasonal variability in nearshore coastal environments. We briefly review the physical and biological processes – acting across a variety of spatial and temporal scales – that are responsible for the latitudinal and nearshore-to-offshore pCO_2(sw) gradients seen in the RFR-CCS reconstruction of pCO_2(sw). RFR-CCS will be valuable for the validation of high-resolution models, the attribution of spatiotemporal carbonate system variability to physical and biological drivers, and the quantification of multiyear trends and interannual variability of ocean acidification.

1 Introduction

The concentration of carbon dioxide gas (CO₂) in Earth's atmosphere has rapidly increased from about 280 parts per million in 1750 to over 400 parts per million today (Joos and Spahni, 2008; Dlugokencky and Tans, 2019). This rise in CO₂ concentration is a direct result of human activities such as fossil fuel combustion, deforestation, and agriculture (Ciais et al., 2014; Friedlingstein et al., 2020). The presence of human-produced or “anthropogenic” CO₂ in the atmosphere – along with other anthropogenic greenhouse gases – leads to planetary warming, with a disproportionate amount of heat (∼ 90 %) being absorbed by the ocean (von Schuckmann et al., 2020). About a quarter of annually produced anthropogenic CO₂ dissolves directly into the ocean (Friedlingstein et al., 2020), mitigating its warming potential. However, dissolved CO₂ reacts with seawater to form carbonic acid, which rapidly dissociates and acidifies (primarily) surface ocean environments (Caldeira and Wickett, 2003), with adverse effects for many marine organisms and ecosystems (Orr et al., 2005; Fabry et al., 2008; Pörtner, 2008; Doney et al., 2009, 2020). Closing the global carbon budget involves accurately estimating the amount of CO₂ taken up by the ocean (e.g., Hauck et al., 2020). A primary method for calculating the amount of CO₂ transferred to the ocean requires knowing the difference between the partial pressure of CO₂ in the atmosphere and surface seawater.

Compared to atmospheric CO₂ partial pressure (pCO_2(atm)), which can be determined with some certainty at a given location even without direct observations due to the well-mixed nature of the atmosphere, surface seawater CO₂ partial pressure (pCO_2(sw)) is more variable and therefore more difficult to constrain (Wanninkhof, 2014; Landschützer et al., 2014; Woolf et al., 2019). This variability is a result of ocean mixing, equilibration kinetics between the atmosphere and ocean, biological processes, and thermal effects on pCO_2(sw). Filling temporal and spatial data gaps in the observational coverage of pCO_2(sw) can therefore be challenging (Hauck et al., 2020; Fay et al., 2021) and a variety of strategies have been attempted over several decades (Takahashi et al., 1993; Rödenbeck et al., 2015), becoming even more prevalent and varied in the literature over time. Briefly, statistical interpolations (Takahashi et al., 1993, 2002, 2009; Rödenbeck et al., 2013, 2014; Jones et al., 2015; Shutler et al., 2016), multiple linear regressions (Schuster et al., 2013; Iida et al., 2015; Becker et al., 2021), machine-learning-based regression methods (Landschützer et al., 2013; 2014, 2016, 2018; Nakaoka et al., 2013; Zeng et al., 2014; Laruelle et al., 2017; Ritter et al., 2017; Gregor et al., 2017, 2018; Chen et al., 2019; Denvil-Sommer et al., 2019), and biogeochemical-model-based approaches (Valsala and Maksyutov, 2010; Majkut et al., 2014; Verdy and Mazloff, 2017) have been common tactics, each one with its own strengths and weaknesses. Recently, ensemble averages of multiple data- or model-based approaches have become popular options as well (Gregor et al., 2019; Lebehot et al., 2019; Fay et al., 2021).

One widely used machine-learning-based pCO_2(sw) gap-filling strategy relies on a two-step approach consisting of unsupervised clustering using a self-organizing-map (SOM) followed by construction of a feed-forward neural network (FFN) for each cluster (Landschützer et al., 2013). This SOM-FFN approach is well-established in the literature (Landschützer et al., 2013, 2014, 2015, 2016, 2018; Laruelle et al., 2017; Ritter et al., 2017; Denvil-Sommer et al., 2019) and is recognized as one of the most effective approaches for filling gaps in the observational pCO_2(sw) record (Rödenbeck et al., 2015). The SOM-FFN approach was recently applied to coastal ocean areas, resulting in the first globally continuous, multiyear data product of monthly coastal ocean pCO_2(sw) at 0.25^∘ resolution (Laruelle et al., 2017). Even more recently, that coastal product was combined with an updated open-ocean product (Landschützer et al., 2020a) to produce a uniform 12-month climatology of pCO_2(sw) across the coastal to open-ocean continuum (Landschützer et al., 2020b, c).

The data products provided by Laruelle et al. (2017) and Landschützer et al. (2020b) – hereafter L17 and L20, respectively – are important advancements toward characterizing pCO_2(sw) across the entire ocean domain for carbon budget analyses. Most data-based estimates of oceanic CO₂ uptake have considered only the open ocean (e.g., Landschützer et al., 2014; Iida et al., 2015; Denvil-Sommer et al., 2019; Gregor et al., 2019; Watson et al., 2020) or are based on coarse spatial representations of the coastal ocean (Rödenbeck et al., 2013). However, coastal ocean CO₂ uptake is estimated to be about 10 % of the open-ocean figure (Laruelle et al., 2010, 2014; Bourgeois et al., 2016; Roobaert et al., 2019; Chau et al., 2022), is far more spatially variable (Liu et al., 2010), and may be changing at a different rate relative to open-ocean CO₂ uptake (Laruelle et al., 2018). Therefore, augmenting global open-ocean pCO_2(sw) data products to include the coastal ocean is quite valuable (Fay et al., 2021). Despite the greater spatial coverage and temporal resolution offered by these new gap-filled pCO_2(sw) data products, significant challenges remain for accurately representing pCO_2(sw).

One of those challenges involves characterizing seasonal cycles in pCO_2(sw), particularly in the nearshore coastal ocean. Although the L17 product effectively captures pCO_2(sw) seasonality when averaged across relatively large coastal ocean regions, the authors assert that “the coastal SOM-FFN tends to systematically underestimate the amplitude of the seasonal pCO₂ cycle” in locations where they can make comparisons with direct observations. This result is logical given that (1) direct observations are made at discrete locations and times, whereas gridded products are averaged over some spatial area and time, which tempers extremes; and (2) fits obtained via least squares regressions or machine-learning methods generally tend to perform better when temporal and spatial variability is low and worse when variability is high (Landschützer et al., 2014), such as in the coastal ocean. However, this problem must be addressed if we hope to achieve realistic global representations of pCO_2(sw) seasonality, which are necessary for investigating the processes that drive this variability (Roobaert et al., 2019) and for ensuring the fidelity of future air–sea CO₂ flux projections (Hauck et al., 2020). Addressing carbon exchange in coastal margins has recently been highlighted as a fundamental and emerging research topic in ocean carbon research (Dai, 2021).

Here, we present a reconstruction of pCO_2(sw) (1998–2020) in a broad region of the northeastern Pacific that includes the California Current System (CCS), the surrounding open-ocean regions, and the highly variable continental shelf of the North American west coast spanning from southern Alaska to Baja California. We apply a random forest regression (RFR) approach (Breiman, 2001) to fill observational gaps, constraining pCO_2(sw) across the coastal to open-ocean continuum. We show that the RFR approach in the northeastern Pacific produces realistic monthly maps of surface pCO_2(sw) from 1998 to 2020 and that these maps reliably capture seasonal pCO_2(sw) variability in the coastal and open ocean.

We compare pCO_2(sw) values from our gap-filled product – RFR-CCS – to coastal ocean mooring measurements and other direct observations and to the available global-scale 0.25^∘ resolution SOM-FFN products in the region (i.e., L17 and L20). We speculate as to why nearshore seasonal cycles are better represented by RFR-CCS than by global-scale gap-filled products and discuss implications for how to best capture seasonal variability in global products going forward. We describe spatial and seasonal patterns in pCO_2(sw) revealed by RFR-CCS and discuss the physical and biological processes that likely produce those patterns. Finally, we compare air–sea CO₂ flux computed from RFR-CCS to that from a recently released CO₂ flux product (Gregor and Fay, 2021) and discuss the implications of sporadic sampling for calculations of CO₂ flux in the coastal ocean.

2 Methods

2.1 Sea surface $f_{{\mathrm{CO}}_{\mathrm{2}}}$ data acquisition and conversion to pCO₂

Sea surface CO₂ fugacity ($f_{{\mathrm{CO}}_{\mathrm{2}}\left(\mathrm{sw}\right)}$) data, along with ancillary variables, were obtained from the Surface Ocean CO₂ Atlas (SOCAT; Pfeil et al., 2013; Bakker et al., 2016) version 2021 (SOCATv2021) for latitudes between 15 and 60^∘ N and longitudes between 105 and 140^∘ W (hereafter referred to as “the study region”). SOCAT is an international effort to synthesize quality-controlled $f_{{\mathrm{CO}}_{\mathrm{2}}\left(\mathrm{sw}\right)}$ observations for the global surface ocean, and has released datasets of individual surface ocean $f_{{\mathrm{CO}}_{\mathrm{2}}\left(\mathrm{sw}\right)}$ observations and gridded values since 2011, with annual releases since 2015. SOCATv2021 contains nearly 30.6 million $f_{{\mathrm{CO}}_{\mathrm{2}}\left(\mathrm{sw}\right)}$ observations globally and over 1.4 million $f_{{\mathrm{CO}}_{\mathrm{2}}\left(\mathrm{sw}\right)}$ observations within the study region.

SOCAT data in the study region were filtered to retain $f_{{\mathrm{CO}}_{\mathrm{2}}\left(\mathrm{sw}\right)}$ observations with a measurement quality control (QC) flag of 2 (“good”) and dataset QC flags of A through D ($f_{{\mathrm{CO}}_{\mathrm{2}}\left(\mathrm{sw}\right)}$ accuracy of 5 µatm or better). This is identical to the QC procedure followed by the SOCAT team for producing gridded data products (Sabine et al., 2013; Bakker et al., 2016). SOCATv2021 provides ancillary variables along with $f_{{\mathrm{CO}}_{\mathrm{2}}\left(\mathrm{sw}\right)}$, including contemporaneous observations of sea surface temperature (SST) and sea surface salinity (SSS), as well as atmospheric pressure at the ocean surface (P_atm) from the National Centers for Environmental Prediction (NCEP) reanalysis; these values were used only for fugacity to partial pressure conversions (Eq. 1). Though SST and SSS are considered surface values, it is important to note that these are primarily underway measurements taken a few meters beneath the surface and that nontrivial differences in temperature and salinity may exist between the measurement depth and the surface (Robertson and Watson, 1992; Donlon et al., 2002; Goddijn-Murphy et al., 2015; Woolf et al., 2016; Ho and Schanze, 2020; Watson et al., 2020). Also, while SST and SSS are not assigned explicit QC flags in SOCAT, these parameters do undergo quality control checks during the calculation of $f_{{\mathrm{CO}}_{\mathrm{2}}\left(\mathrm{sw}\right)}$ (Lauvset et al., 2018).

Sea surface CO₂ fugacity represents CO₂ partial pressure corrected for the nonideality of CO₂ gas. It was converted to sea surface CO₂ partial pressure (pCO_2(sw)) following (Weiss, 1974)

$$\begin{array}{}\text{(1)}& p{\mathrm{CO}}_{\mathrm{2}\left(\mathrm{sw}\right)}=f_{{\mathrm{CO}}_{\mathrm{2}}\left(\mathrm{sw}\right)}\cdot \exp {\left(P_{\mathrm{atm}}{\displaystyle \frac{B+\mathrm{2}\cdot \mathit{\delta }}{R\cdot T}}\right)}^{-\mathrm{1}},\end{array}$$

where B and δ are virial coefficients, R is the ideal gas constant, and T is SST in Kelvin.

Figure 1 Annual mean pCO_2(sw) from the 0.25^∘ resolution gridded dataset computed as an average over the monthly climatology from 1998 to 2020 for each grid cell. The two extremes of the color bar can represent pCO_2(sw) values less than or greater than the color bar limits; the chosen range represents most of the values and emphasizes regional contrast.

2.2 Binning of pCO_2(sw) observations

Sea surface CO₂ partial pressure data were aggregated onto a 0.25^∘ latitude by 0.25^∘ longitude grid for each month from January 1998 to December 2020 using a bin-averaging procedure that consisted of computing the means (μ) and standard deviations (σ) of all observations of pCO_2(sw) included within each grid cell. Observations prior to 1998 were excluded as an increase in $f_{{\mathrm{CO}}_{\mathrm{2}}\left(\mathrm{sw}\right)}$ data coverage occurs around the start of 1998 and the first full year of SeaWiFS chlorophyll observations (which are used in our procedure to fill gaps in the pCO_2(sw) dataset) is 1998. For cases in which observations in a given grid cell originated from two or more platforms (e.g., cruises or autonomous assets), platform-weighted μ and σ were computed by first taking the means and standard deviations of all observations made by each platform, then taking the means of those values. This ensured that all platforms contributing observations to a given grid cell were weighted equally, mitigating unwanted biases toward high-resolution measurement systems (Sabine et al., 2013).

Table 1 Sources of data for interpolation of surface pCO₂. Chlorophyll a (Chl) and mixed layer depth (MLD) were log₁₀-transformed to produce a distribution of values that was closer to normal before constructing the regression model. Gaps in CHL data were filled by linear interpolation over time within each grid cell (see Appendix A). Month of the year was transformed by cosine and sine functions to retain its cyclical nature.

This bin-averaging procedure is identical to the one followed by the SOCAT team for producing monthly datasets for coastal regions with 0.25^∘ resolution as well as for open-ocean regions with 1^∘ resolution (Sabine et al., 2013; Bakker et al., 2016). However, here we produced a monthly gridded dataset with 0.25^∘ resolution for a region of the northeastern Pacific (15 to 60^∘ N, 105 to 140^∘ W) that spans both the coastal and open ocean. Means of pCO_2(sw) from this gridded dataset (averages over the monthly climatology from 1998 to 2020 for each spatial grid cell) are shown in Fig. 1. Some of the apparent fine-scale spatial variability in this bin-averaged map is not indicative of true environmental conditions but originates from the combination of large temporal variability within each grid cell and uneven sampling of each grid cell across and within years. This form of temporal variability is exactly the kind of spurious result that advanced pCO_2(sw) mapping techniques are intended to circumvent. Figure B1 shows the number of years containing an observation within each month of our gridded pCO_2(sw) dataset. Unsurprisingly, temporal coverage is highest close to the coast, especially in the summer months.

2.3 Predictor variable acquisition and processing

Of the 4 014 844 grid cells that represent the surface ocean gridded in three dimensions at 0.25^∘ resolution over 276 months (1998–2020) in the study region, just 1.25 % have an associated gridded pCO_2(sw) value. To fill gaps in this dataset, relationships between pCO_2(sw) and various predictor variables need to be determined. The predictor variables used in this study are primarily derived from satellite observations or reanalysis models due to the condition that they be resolved with temporal and spatial continuity across the study region and selected time span.

Predictor variables are intended to capture conditions that mechanistically influence pCO_2(sw) (e.g., SST and atmospheric pCO₂), serve as a proxy for mechanisms that influence pCO_2(sw) (e.g., sea surface chlorophyll), or, in the case of temporal and spatial information, constrain additional patterned variability not captured by the mechanistic variables alone. The chosen predictor variables for this study (Table 1) have all been used before for pCO_2(sw) gap-filling methods (e.g., Landschützer et al., 2014; Gregor et al., 2018; Denvil-Sommer et al., 2019; Watson et al., 2020); temporal and spatial predictors were included to ensure robust representation of pCO_2(sw) seasonal cycles (Gregor et al., 2017). Included in Table 1 are the sources of each dataset, the original resolutions of each dataset, and the steps that were taken to process each dataset. Appendix A provides more detail about the acquisition and processing of the driver variables and includes figures showing annual means of selected variables.

2.4 Construction of nonlinear relationships using random forest regression

We used the random forest regression approach (Breiman, 2001) to identify relationships between pCO_2(sw) and predictor variables in order to fill gaps in the gridded pCO_2(sw) dataset. This method averages the results from a number of decision and/or regression trees (i.e., a “forest”) built on bootstrapped replicates of the dataset – which individually have low bias and high variance – to produce a final regression model with reduced variance (Hastie et al., 2009). RFR is the machine-learning method of choice for this study as early testing showed better performance than the SOM-FNN method in the northeastern Pacific. Further, RFR is less computationally expensive than fitting a neural network and has been shown to produce results comparable to the SOM-FFN approach in terms of overall performance (Gregor et al., 2017). It should be noted, however, that the two approaches differ mechanistically and therefore adapt to variability within a training dataset in different ways. Finally, while RFR has been explored more frequently in recent years as a method of spatiotemporal pCO_2(sw) gap-filling both globally (Gregor et al., 2017, 2018) and regionally in the Gulf of Mexico (Chen et al., 2019), far fewer RFR-based pCO_2(sw) products exist than neural-network-based products. So, this study provides a good opportunity to further demonstrate the utility of RFR for producing monthly fields of pCO_2(sw), in this case on a regional scale in the northeastern Pacific.

Each decision tree within a random forest regression model is built on a different subset of the training dataset (that contains both the predictor variables and corresponding gridded pCO_2(sw) values). This subset is generated by bootstrapping, in which a random set of training data points is selected with replacement – meaning the same data point can be selected more than once (Breiman, 1996). The number of data points in the bootstrapped dataset is equal to a defined fraction (InBagFraction in Table 2) of the original dataset; however, a fraction equal to 1 does not mean the bootstrapped dataset is identical to the original dataset because selection is made with replacement. Since each regression tree is built on a different subset of the training data, it will contain somewhat different relationships between the predictor variables and the corresponding gridded pCO_2(sw) values.

Table 2 Model parameters for the random forest regression. Parameter names are the default property names for the MATLAB TreeBagger class.

The process of building a decision tree begins at the top “node” of the tree with the values of a single predictor variable being used to split that tree's bootstrapped subset of the training dataset into two smaller subsets (not necessarily of equal size) containing the most similar pCO_2(sw) observations (i.e., sets of pCO_2(sw) observations with the smallest variance among them). These subsets are then further divided into progressively smaller sets of similar observations until either the variance among the pCO_2(sw) observations in a node drops below a prescribed tolerance level or the number of observations in the node reaches the user-defined minimum (MinLeafSize in Table 2). To ensure that the algorithm does not always pick the same predictor variable (e.g., the one most highly correlated with pCO_2(sw) overall) for the split at every node, we limit it to choosing from a different random subset of the predictor variables (equal in number to NumPredictorsToSample in Table 2) at each node. This introduces another “random” element into the tree-building process. The random forest contains a large number of these regression trees (NumTrees in Table 2) each built on a different, random bootstrapped subsample of the training data. Once the random forest is built, a set of predictor variables can be provided to the model and the average of the pCO_2(sw) values provided by each regression tree is used as the pCO_2(sw) prediction for that particular set of inputs.

We constructed an RFR model using the MATLAB TreeBagger function with the predictor variables given in Table 1 and the parameters given in Table 2, along with gridded pCO_2(sw) values that were obtained as described in Sect. 2.1 and 2.2. To produce the northeastern Pacific random forest regressionpCO_2(sw) product (RFR-CCS) that is the main result of this work (Sharp et al., 2022; https://doi.org/10.5281/zenodo.5523389), the full dataset of gridded pCO_2(sw) values was used. For optimization and evaluation, subsets of the full dataset were used as described in the following sections.

2.5 Optimization of random forest regression model

The predictor variables used (Table 1) and the values for the model parameters (Table 2) were determined by iteratively optimizing the model performance. First, default model parameters were used to train an RFR model using a subset of the data for training (80 % of full dataset, distributed randomly across the space and time domains of interest) and a number of possible predictor variables: latitude, longitude, sea surface height, bottom depth, and those given in Table 1. During model selection, the generalization skill for the RFR model was assessed using a validation dataset comprised of 10 % of the full dataset, none of which was included in the training data. After the initial model fit, predictors with a “feature importance” (computed during the RFR fit) significantly lower than all other predictors were sequentially dropped (latitude, longitude, and sea surface height), and this did not substantially change the training or validation root mean squared error (RMSE). Remaining predictor variables were dropped one at a time for subsequent fits, and the goodness-of-fit and generalization skill of the model were assessed using the RMSE values calculated from applying the model to the training and validation datasets, respectively. The set of predictors with the lowest RMSE after dropping one predictor was carried into the next iteration. If removing a predictor did not increase the validation RMSE significantly, then that predictor was removed from the set of predictors (only bottom depth was dropped in this step). The final set of predictor variables is shown in Table 1.

Next, different values for model parameters (Table 2) were tried iteratively with the retained predictors to identify the optimal values, again by minimizing the RMSE of the validation dataset. Although lower values for the minimum terminal node size performed better in this analysis, additional testing indicated that retaining the default value of 5 was important to prevent overfitting. To determine the appropriate number of trees, we examined how the out-of-bag mean squared error changed as more and more trees were included in the random forest (up to 5000 trees) and selected a number of trees well past the point at which this error had stabilized (1200 trees). Finally, the remaining 10 % of the full dataset that was withheld from both the model training and model validation (i.e., the “test data”) was used to quantify the mapping uncertainties from the RFR approach (discussed further in Sects. 2.7 and 3.5). The predictor variable feature importances for the final RFR-CCS fit are given in Fig. B2.

2.6 Evaluation of random forest regression approach and resulting data product

Once predictor variables and model parameters were optimized, the skill of the RFR approach was further evaluated by splitting the full dataset into different subsets of training data and test data. Evaluation models (RFR-CCS-Evals) were constructed in three different ways: (1) by removing a random (20 %) subset of cruises and/or measurement platforms from the training data (repeated 10 times with different subsets removed each time; n=10), (2) by removing all observations from every fifth year from the training data (repeated five times such that data from every year was removed from one of the trials; n=5), and (3) by removing all moored autonomous pCO_2(sw) measurements (i.e., discrete time series sites primarily located in the coastal ocean) from the training data (n=1). The first two strategies were relevant for assessing bulk error statistics for the method applied across the region and the third strategy for evaluating the ability of the RFR to represent local seasonal variability without the use of high-temporal-resolution mooring data. These RFR-CCS-Evals are distinct model variants that are only used for assessment; the final RFR-CCS model uses all available training data.

Each data split for an RFR-CCS-Eval was applied directly to SOCATv2021 observations, before bin-averaging the data according to the procedure given in Sect. 2.2; as a result, a gridded training dataset and a gridded test dataset were produced from each split. Data splits were performed in this way to ensure that autocorrelation among measurements from a specific platform did not bias the error statistics. Each split was repeated n times, and error statistics (bias, RMSE, and R²) from comparing pCO_2(sw) predicted from the RFR-CCS-Eval models versus pCO_2(sw) from the gridded test dataset were averaged.

The final RFR-CCS data product was evaluated through comparisons with gridded pCO_2(sw) observations from SOCATv4 (Bakker et al., 2016) and pCO_2(sw) observations from surface ocean moorings (Sutton et al., 2019). Of the surface ocean moorings within the study site that are not located within an inland sea and have available data from all 12 months of the year, four (CCE2, NH10, Cape Elizabeth, Châ bá) are located within 40 km of shore and one (CCE1) is about 215 km from shore. RFR-CCS was also compared to global-scale gap-filled pCO_2(sw) products that are available in the region. Namely, we focused on the coastal multi-month pCO_2(sw) product from Laruelle et al. (2017; i.e., L17) and the combined coastal and open-ocean pCO_2(sw) climatology from Landschützer et al. (2020b; i.e., L20).

2.7 Uncertainty analysis

Uncertainty in pCO_2(sw) for each grid cell was calculated according to the approach used by Landschützer et al. (2014, 2018) and Roobaert et al. (2019), in which total uncertainty in pCO_2(sw) results from a combination of observational uncertainty, mapping uncertainty, and gridding uncertainty. Observational uncertainty (θ_obs) is uncertainty inherent to the original measurements of pCO_2(sw) evaluated as the average of reported uncertainties in the $f_{{\mathrm{CO}}_{\mathrm{2}}\left(\mathrm{sw}\right)}$ observations from our training dataset, which are flagged by SOCAT with a dataset QC flag of A or B ($f_{{\mathrm{CO}}_{\mathrm{2}}\left(\mathrm{sw}\right)}$ accuracy of 2 µatm or better) and of C or D ($f_{{\mathrm{CO}}_{\mathrm{2}}\left(\mathrm{sw}\right)}$ accuracy of 5 µatm or better); we weighted θ_obs by the number of observations assigned each flag. Mapping uncertainty (θ_map) is uncertainty contributed by the RFR mapping procedure and was evaluated as separate values for the coastal (< 400 km from shore) and open ocean (> 400 km from shore) using the mean of the root mean squared errors for a subset of test data (10 %) withheld from both the model training data (80 %) and model validation data (10 %) (see Sect. 2.5). Gridding uncertainty (θ_grid) is uncertainty attributable to aggregating observations into monthly 0.25^∘ resolution grid cells and was evaluated as separate values for the coastal and open ocean by taking the average unweighted standard deviation among pCO_2(sw) values within each grid cell in which two or more platforms were represented. Grid cells with mooring observations were excluded from the θ_grid calculation to avoid the high number of observations swamping the signal from other platforms. These three components were combined to obtain total pCO_2(sw) uncertainty (${\mathit{\theta }}_{p{\mathrm{CO}}_{\mathrm{2}}}$) applicable to each open-ocean grid cell and to each coastal grid cell:

$$\begin{array}{}\text{(2)}& {\mathit{\theta }}_{p{\mathrm{CO}}_{\mathrm{2}}}=\sqrt{{\mathit{\theta }}_{\mathrm{obs}}^{\mathrm{2}}+{\mathit{\theta }}_{\mathrm{map}}^{\mathrm{2}}+{\mathit{\theta }}_{\mathrm{grid}}^{\mathrm{2}}}.\end{array}$$

Whereas ${\mathit{\theta }}_{p{\mathrm{CO}}_{\mathrm{2}}}$ represents the uncertainty in pCO_2(sw) for a given grid cell in a given month, uncertainty averaged regionally or over time will not scale exactly with ${\mathit{\theta }}_{p{\mathrm{CO}}_{\mathrm{2}}}$ due to the spatial correlation of pCO_2(sw) values and the autocorrelation features of the model error (e.g., Landschützer et al., 2014).

2.8 Calculation of CO₂ flux

The flux of CO₂ across the ocean–atmosphere interface ($F_{{\mathrm{CO}}_{\mathrm{2}}}$) was calculated using a bulk formula:

$$\begin{array}{}\text{(3)}& F_{{\mathrm{CO}}_{\mathrm{2}}}=k_{\mathrm{w}}\times K_{\mathrm{0}}\times \mathrm{\Delta }p{\mathrm{CO}}_{\mathrm{2}},\end{array}$$

where k_w is the gas transfer velocity, K₀ is the CO₂ solubility constant, and ΔpCO₂ is the difference between CO₂ partial pressure in seawater and in the overlying atmosphere (pCO_2(sw)−pCO_2(atm)) The salinity- and temperature-dependent equations of Weiss (1974) were used to calculate K₀.

Gas transfer velocities were parameterized using a quadratic dependence on wind speed (Wanninkhof, 1992):

$$\begin{array}{}\text{(4)}& k_{\mathrm{w}}={\mathrm{\Gamma }}_{\mathrm{660}}\cdot U^{\mathrm{2}}\cdot \sqrt{\mathrm{660}}/Sc,\end{array}$$

where Γ₆₆₀ is a gas exchange coefficient normalized to Sc = 660, U² is the squared wind speed, and Sc is the Schmidt number for CO₂. Our calculations used Γ₆₆₀=0.276, which is a gas exchange coefficient that is specific to ERA5 reanalysis winds and scaled to a bomb-¹⁴C flux estimate of 16.5 cm h⁻¹ (Fay et al., 2021). Sc was calculated using the fourth-order polynomial fit of Wanninkhof (2014). U was obtained from ERA5 reanalysis (Hersbach et al., 2020). Flux calculations used monthly averages of squared 3-hourly wind speeds to retain the influence of the quadratic wind term (Fay et al., 2021).

3 Results and discussion

3.1 Evaluation by comparison to withheld data

As described in Sect. 2.6, training and test datasets were created by splitting the full dataset prior to bin-averaging. Evaluation models (RFR-CCS-Evals) were constructed by fitting RFR models using the various gridded training datasets. Values of pCO_2(sw) predicted by RFR-CCS-Evals were compared to corresponding values from gridded test datasets. Error statistics (bias, RMSE, and R²) averaged over the n sets of evaluation tests are given in Table 3. When RFR-CCS is compared against all the gridded observations used to construct it, error statistics are predictably strong (last row in Table 3), with a mean bias of 0.00 µatm and an RMSE of 13.33 µatm (R²=0.93). These error statistics demonstrate the ability of the RFR model to fit the training data; the evaluation tests provide insight into the model's ability to predict independent data.

Table 3 Error statistics for comparisons of predicted pCO_2(sw) from evaluation models versus gridded pCO_2(sw) from test datasets. The number of times each test was repeated is given by n; where n is greater than 1, different subsets of data were removed for each iteration of the test and error statistics are the mean of all iterations.

^∗ These statistics represent model training statistics (i.e., evaluated with the same data used to train the model) rather than model validation statistics.

Tests 1 and 2 are good indicators of the overall skill of RFR-CCS. The mean absolute bias for each of those tests is less than 2 µatm, and the RMSEs are near or below 30 µatm. These error statistics can be compared with those of L17, who obtained biases with a mean of 0.0 and RMSEs ranging from 20.5 to 53.1 µatm (mean of 39.2 µatm) for independent evaluations of coastal pCO_2(sw) values fit using the SOM-FFN method in 10 separate global subregions at 0.25^∘ resolution. For an open-ocean comparison, Denvil-Sommer et al. (2019) obtained an RMSE of 15.86 for an independent evaluation of pCO_2(sw) values fit using a similar neural network approach (LSCE-FFNN) for the subtropical North Pacific (18 to 49^∘ N) at 1^∘ resolution. The error statistics for our study region, which spans the coastal to open-ocean continuum on a finely resolved spatial grid, lie comfortably between those coastal and open-ocean comparison points.

Test 3 is a good indicator of how well the RFR approach is able to reproduce the values and seasonalities of coastal pCO_2(sw) at fixed locations when mooring data at a given location are not provided as training data, as each of the moorings makes continuous pCO_2(sw) measurements throughout the year and all but one of the mooring locations included in SOCATv2021 in this region are within 40 km of shore. The positive mean bias (8.39 µatm) suggests that RFR-CCS somewhat overestimates pCO_2(sw) at grid cells corresponding to mooring locations, but this is strongly influenced by high biases at the Cape Elizabeth and Châ bá mooring locations (Table B1). The relatively high RMSE (43.28 µatm) is a result of higher variability in coastal grid cells compared to the open ocean; this is confirmed by a comparison to the offshore CCE1 mooring (Table B1), where the RMSE from the mooring-excluded RFR-CCS-Eval is just 10.5 µatm.

Figure 2 Monthly values of pCO_2(sw) from mooring observations (black), RFR-CCS (blue), the mooring-excluded RFR-CCS-Eval model (orange), and L17 (green). The envelope around the black line equals the standard deviation of all mooring observations within each month, representing the natural variability of the 3-hourly mooring measurements; the envelopes around the blue and orange lines represent the RFR-CCS and RFR-CCS-Eval results plus 1 standard uncertainty (43.6 µatm; Sect. 3.5); the envelope around the green line represents the L17 data product plus the RMSE of an independent data evaluation in the province associated with CCE2 (52.5 µatm; Table 3 of Laruelle et al., 2017; Province P7).

Figure 2 provides an example of one coastal mooring record (CCE2, which is positioned on the shelf break off the coast of Point Conception, CA, at 34.324^∘ N, 120.814^∘ W) compared to pCO_2(sw) predicted in the corresponding grid cell (centered at 34.375^∘ N, 120.875^∘ W) by the mooring-excluded RFR-CCS-Eval model (Test 3) as well as the full RFR-CCS model. For comparison, pCO_2(sw) in the same grid cell provided by the L17 coastal product is also shown. At the CCE2 mooring location, RFR-CCS reproduces mooring-observed monthly pCO_2(sw) with a mean bias of −2.2 µatm and an RMSE of 16.1 µatm (R²=0.81). These error statistics are expected to be relatively favorable, as the RFR-CCS model is trained using mooring observations from CCE2. In contrast, the mooring-excluded RFR-CCS-Eval reproduces monthly mooring-observed pCO_2(sw) at CCE2 with a mean bias of −4.6 µatm and an RMSE of 28.9 µatm (R²=0.41). This can be compared to the L17 coastal pCO_2(sw) product, which reproduces monthly mooring-observed pCO_2(sw) at CCE2 with a mean bias of −44.2 µatm and an RMSE of 57.3 µatm (R²=0.06). Notably, the mooring-excluded RFR-CCS-Eval captures pCO_2(sw) variability at CCE2 more effectively than the L17 product, even though RFR-CCS-Eval was trained without mooring observations and the L17 training dataset (i.e., SOCATv4) includes CCE2 mooring observations through 2014. Similar results are obtained for comparisons to other mooring records (Table B1; Fig. B3), with RFR-CCS always producing the best error statistics (as expected) and RFR-CCS-Eval always producing a better R² than L17, indicating that coastal seasonality at mooring locations is better captured by our regional random forest regression model, even when mooring observations themselves are not included in the model training. This is an important conclusion, especially in light of the recommendation by Hauck et al. (2020) that the inclusion of coastal areas and marginal seas in pCO_2(sw) mapping methods will be critical for improving the ocean carbon sink estimate. If these areas are to be included, it is sensible to attempt to capture their unique modes of variability as accurately as possible.

Figure 3 Differences between annual means (a, c) and seasonal amplitudes (b, d) of pCO_2(sw) from RFR-CCS-clim versus the L20 climatology (a, b; RFR-CCS-clim – L20) and versus a climatological average of the L17 product (c, d; RFR-CCS-clim – L17).

3.2 Evaluation by comparison to global pCO_2(sw) products

Across the study area, values of pCO_2(sw) from RFR-CCS were compared against corresponding values from L17 and L20. For temporal compatibility with L17 and L20, a climatology of average monthly values from RFR-CCS spanning 1998 to 2015 (RFR-CCS-clim) was created for these comparisons. Figure 3 shows mapped differences in annual means and seasonal amplitudes (calculated as the maximum climatological pCO_2(sw) minus the minimum) of pCO_2(sw) between RFR-CCS-clim versus L20 (top panels) and RFR-CCS-clim versus a climatological average of L17 (bottom panels); monthly mean differences in are given in Fig. B4.

The most notable feature of the annual mean difference maps is that RFR-CCS-clim produces much higher annual mean pCO_2(sw) than both L17 and L20 in the nearshore coastal ocean and slightly higher pCO_2(sw) in the remainder of the study area. Similarly, RFR-CCS-clim produces much higher seasonal variability than both L17 and L20 in the nearshore coastal ocean, especially north of about 34^∘ N. On average, RFR-CCS-clim produces an area-weighted annual mean pCO_2(sw) that is greater than L17 by 19.0 µatm and L20 by 8.4 µatm, as well as an area-weighted seasonal amplitude that is greater than L17 by 13.0 µatm and L20 by 5.6 µatm.

3.3 Evaluation by comparison to gridded observations of pCO_2(sw)

Values of pCO_2(sw) from RFR-CCS, L17, and L20 were compared against the SOCATv4 gridded pCO_2(sw) data product. SOCATv4 was used in the development of the coastal L17 product, whereas SOCATv5 was used in the development of the open-ocean product for the merged L20 climatology, and SOCATv2021 was used in the development of RFR-CCS. Therefore, comparisons were made to both the gridded open-ocean observations (1^∘ resolution) and gridded coastal observations (0.25^∘ resolution) from SOCATv4 to include only data points that were available to the training of all three data products. To match the resolution of the gridded open-ocean observations from SOCATv4, aggregation from a 0.25^∘ resolution grid to a 1^∘ resolution grid was performed for RFR-CCS, RFR-CCS-clim, and L20. L17 was only compared to gridded coastal observations from SOCATv4 because the two are gridded to the same spatial resolution and cover the same coastally limited spatial domain.

Figure 4 shows two-dimensional histograms of bin-averaged differences between RFR-CCS-clim, L20, RFR-CCS, and L17, each compared against gridded observations from SOCATv4. For comparisons to climatological products (RFR-CCS-clim and L20), gridded SOCATv4 observations were averaged to a monthly climatology across 1998–2015 for consistency with the products. The regional RFR-CCS product and its climatology outperform both global SOM-FFN products: RFR-CCS-clim shows better agreement with gridded monthly means of observations from SOCATv4 than L20 (R²=0.85 versus R²=0.73), and RFR-CCS (within the coastally limited spatial domain of L17) shows better agreement with gridded observations from SOCATv4 than L17 (R²=0.96 versus $R^{\mathrm{2}}=\mathrm{0}.$61). In particular, the two global products (L20 and L17) struggle to match pCO_2(sw) values in the nearshore coastal ocean (within 100 km of the coast), indicated by dark blue cells in Fig. 4.

Mismatches between global pCO_2(sw) products and observations in the nearshore coastal ocean are not unexpected, as regional error statistics for reconstructed global pCO_2(sw) are typically larger than the global mean error statistics (Laruelle et al., 2017; Landschützer et al., 2020c), and it is generally more challenging to model pCO_2(sw) in environments with high temporal and spatial variability, such as in the nearshore coastal ocean (Landschützer et al., 2014). This result emphasizes the importance of carefully addressing nearshore pCO_2(sw) when constructing global products if one hopes to achieve an accurate representation of coastal ocean variability. This may be achieved (1) by using a greater number of model clusters for coastal ocean reconstructions (L17 uses just 10 biogeochemical clusters for the global coastal ocean), (2) by increasing the spatial and/or temporal resolution of pCO_2(sw) data products to better account for small-scale variability (Gregor et al., 2019), (3) by carefully accounting for mismatches between the temperature (and salinity) at which pCO_2(sw) is measured versus that at which it is reported in surface data products (Ho and Schanze, 2020; Watson et al., 2020), or (4) by taking an ensemble approach to pCO_2(sw) gap-filling to reduce errors overall, especially in undersampled regions (Gregor et al., 2019; Fay et al., 2021). Ultimately, it will be critical to continue to expand our observational capabilities by means of shipboard underway systems (Pierrot et al., 2009), uncrewed surface vehicles (Meinig et al., 2015; Sutton et al., 2021), biogeochemical Argo floats (Roemmich et al., 2019), moored buoys (Sutton et al., 2019), and other platforms, as well as to make strides toward incorporating these novel measurements into pCO_2(sw) gap-filling schemes (Gregor et al., 2019; Djeutchouang et al., 2022).

Figure 4 Two-dimensional histograms showing bin-averaged comparisons of pCO_2(sw) from (a) RFR-CCS-clim and (b) L20 to SOCATv4 gridded observations that have been averaged to a climatology, as well as comparisons of pCO_2(sw) from (c) RFR-CCS and (d) L17 to SOCATv4 gridded monthly observations in the coastally limited spatial domain of L17. Grid cells are color-coded by the average base-10 logarithm of distance from shore (km) of the observations included within each bin; the transparency of each grid cell is set by the relative number of observations within each bin.

Figure 5 Climatological mean pCO_2(sw) from five NOAA ocean moorings and the corresponding grid cells in RFR-CCS-clim, L20, and a climatological average of L17. Shading represents the standard deviation of all monthly values for each mooring or data product.

3.4 Evaluation by comparison to seasonal observations of pCO_2(sw) at ocean moorings

Values of pCO_2(sw) from RFR-CCS-clim, L17, and L20 were compared against monthly climatologies from mooring observations to evaluate how well each product captured seasonal variability at fixed time series sites. Figure 5 shows climatologies of mooring-observed pCO_2(sw) (each averaged over available years and normalized to their annual mean) compared to pCO_2(sw) from RFR-CCS-clim, L20, and climatological monthly averages of L17 (each normalized to their annual mean) in the grid cell corresponding to the mooring location. Overall, RFR-CCS-clim does a much better job of capturing the variability in mooring observations than either L17 or L20 (Table 4).

Table 4 Seasonal amplitudes of pCO_2(sw) (µatm) from mooring observations and corresponding grid cells of climatological averages (from 1998–2015) of RFR-CCS-clim, L17, and L20.

3.5 Uncertainty calculations

Three components comprised the estimate of uncertainty for pCO_2(sw) values from RFR-CCS: observational uncertainty (θ_obs), mapping uncertainty (θ_map), and gridding uncertainty (θ_grid). According to the procedure detailed in Sect. 2.7, θ_obs was calculated as 3.3 µatm, θ_map as 4.4 µatm for the open ocean and 35.3 µatm for the coastal ocean, and θ_grid as 3.7 µatm for the open ocean (n=268) and 25.1 µatm for the coastal ocean (n=889). These three components were combined to obtain total pCO_2(sw) uncertainty (${\mathit{\theta }}_{p{\mathrm{CO}}_{\mathrm{2}}}$) according to Eq. (2), resulting in ${\mathit{\theta }}_{p{\mathrm{CO}}_{\mathrm{2}}}$ equal to 6.6 µatm for the open ocean and 43.4 µatm for the coastal ocean. The open-ocean value determined through this analysis compares well with the grid-level uncertainty estimated in open-ocean grid cells by Landschützer et al. (2014), which ranged from 8.6 to 17.7 µatm for different regions. The large coastal uncertainty value emphasizes the high degree of variability in monthly pCO_2(sw) near ocean margins.

As noted in Sect. 2.7, uncertainties reported here are appropriate for a given grid cell (i.e., monthly 0.25^∘ latitude by 0.25^∘ longitude bin). Values averaged over time or over larger regions will have reduced pCO_2(sw) (and CO₂ flux) uncertainties due to the spatiotemporal correlation of pCO_2(sw) and the autocorrelation features of the model error (e.g., Landschützer et al., 2014).

3.6 Spatial and seasonal patterns of sea surface pCO₂

In the open ocean, relatively high pCO_2(sw) values can be observed off southern Baja California (Fig. 6a) and extending toward the northwest, especially during summer months and into autumn (Fig. 7) when higher sea surface temperatures drive higher pCO_2(sw) (Nakaoka et al., 2013). This area also corresponds to low chlorophyll (Fig. A2) and the lowest wind speeds across the study region (Fig. A4), suggesting that a lack of nutrient delivery from deep convection may be limiting biological production, also driving high pCO_2(sw). Relatively low open-ocean pCO_2(sw) values can be observed in the northern part of the study region from about 45 to 60^∘ N (Fig. 6a). Wintertime cooling drives low pCO_2(sw) in this area, though that effect is compensated for by dissolved inorganic carbon (DIC) brought to the surface by deep winter mixing (Ishii et al., 2014). Figure B5 illustrates competing effects between temperature and winds by displaying correlations between SST and pCO_2(sw), which are mainly positive below 50^∘ N, and between wind speed and pCO_2(sw), which are mainly positive above 50^∘ N.

Figure 6 Annual mean pCO_2(sw) (a) and the seasonal amplitude of pCO_2(sw) (b) from RFR-CCS. Also shown are annual mean pCO_2(sw) and the seasonal amplitude of pCO_2(sw) measured at ocean mooring locations.

In the summer, high biological production in the northern portion of the study region (Fig. A2) removes DIC, keeping pCO_2(sw) relatively low. This low-pCO_2(sw) region extends southward along the California coast to about 34^∘ N between both offshore and nearshore high-pCO_2(sw) waters. The southward extension of the low-pCO_2(sw) region is consistent with what we know about the dynamics of the CCS: a narrow band of nearshore waters is high in DIC in the spring and summer due to the direct effects of wind-driven upwelling (Fig. 7), but a wider band of waters farther offshore is lower in DIC due to drawdown by high biological production stimulated by nutrients delivered to the euphotic zone by upwelling (Hales et al., 2005; Fassbender et al., 2011; Fiechter et al., 2014; Turi et al., 2014).

In the coastal ocean, high pCO_2(sw) occurs in the central CCS (∼ 34 to ∼ 42^∘ N), with values of 400 µatm or greater beginning in April off Pt. Conception (34^∘ N) and propagating northward to around Cape Arago (43^∘ N) through October (Fig. 7). This corresponds to the latitudinal band of the CCS with the strongest and most consistent equatorward winds (Huyer, 1983), which induce upwelling of CO₂-rich subsurface waters by wind-driven Ekman transport very near the coast and wind-stress-curl-driven Ekman pumping farther offshore (Checkley and Barth, 2009). This nearshore band of high summertime pCO_2(sw) has been previously reported by observational (Hales et al., 2012) and modeling (Fiechter et al., 2014; Turi et al., 2014; Deutsch et al., 2021) studies. It corresponds to naturally low surface pH values and aragonite saturation states, which will be exacerbated by increasing atmospheric CO₂ concentrations (Gruber et al., 2012; Hauri et al., 2013), with likely deleterious effects for calcifying organisms (Feely et al., 2008).

Figure 7 Monthly mean pCO_2(sw) fields from RFR-CCS.

Relatively low coastal pCO_2(sw) values (340 µatm or lower) develop during April off the coasts of Oregon, Washington, and Vancouver Island and propagate northward toward southern Alaska through September (Fig. 7). Low summertime pCO_2(sw) in the northern CCS (∼ 42 to ∼ 50^∘ N) has been demonstrated before (Hales et al., 2005, 2012; Evans et al., 2011; Fassbender et al., 2018) and corresponds to the weaker and more variable equatorward winds in summer in the northern CCS (Checkley and Barth, 2009) as well as the effect of DIC drawdown by high primary productivity, which offsets upwelling-induced increases in pCO_2(sw). Primary productivity in the northern CCS can be enhanced relative to the rest of the CCS due to factors like riverine nutrient delivery and distribution, submarine canyon-enhanced upwelling, and physical retention of phytoplankton blooms (Hickey and Banas, 2008).

The coastal ocean from Vancouver Island northward is a high-pCO_2(sw) region from October to March (Fig. 7), which is broadly consistent with observations of high pCO_2(sw) in the western Canadian coastal ocean during autumn and winter (Evans et al., 2012, 2022). This high pCO_2(sw) is perhaps due to the influence of deep tidal mixing (Tortell et al., 2012) and wintertime light limitation of DIC drawdown by primary production. The northern coastal area shifts to a low-pCO_2(sw) region from April to September, again consistent with observations (Evans et al., 2012, 2022) and likely reflecting surface DIC drawdown by primary production in the region (Ianson et al., 2003).

The coastal ocean from the Southern California Bight (SCB) southward along Baja California (∼ 22 to ∼ 34^∘ N) shows relatively low pCO_2(sw) seasonality (Fig. 6b). In this region, pCO_2(sw) is generally lower than in offshore waters of the same latitude, which matches previous results well (Fig. 6a; Hales et al., 2012; Deutsch et al., 2021). One exception is directly off the southern tip of Baja California, where especially high summertime pCO_2(sw) is observed. This may in part reflect the tendency for wind-driven upwelling to bring significant amounts of CO₂-rich subsurface waters to the surface just south of major topographic features (Van Geen et al., 2000; Friederich et al., 2002; Fiechter et al., 2014). Coastal pCO_2(sw) within the Gulf of California (GoC) appears to be strongly influenced by thermally induced seasonal effects, though the lack of observational data coverage in the GoC within SOCATv2021 (Fig. 1), especially within the nonsummer months (Fig. B1), may mask more dynamic variability.

The seasonal amplitude of pCO_2(sw) (Fig. 6b) exhibits interesting variation in the central and northern CCS. Here, nearshore seasonality is extremely high due to dominant effects from upwelling and primary production; however, seasonality farther offshore is extremely low, likely due to compensating effects by thermally driven changes to pCO_2(sw) (high temperature in summer increases pCO_2(sw), low temperature in winter decreases pCO_2(sw)) and biologically or physically driven changes to pCO_2(sw) (high primary production in summer decreases pCO_2(sw), deep mixing in winter increases pCO_2(sw)). Elsewhere, a hotspot of high seasonality exists offshore around 40^∘ N, possibly due to thermal control of pCO_2(sw) without strong biophysical compensatory effects.

3.7 Carbon uptake in the RFR-CCS domain

A recently published data product (SeaFlux; Gregor and Fay, 2021) described by Fay et al. (2021) harmonizes calculations of global CO₂ flux by standardizing the areas covered by different global pCO_2(sw) products and by scaling the gas exchange coefficient to different wind products. As part of this procedure, the L20 climatology is used to fill spatial gaps in some of the pCO_2(sw) products. As we have demonstrated here, filling gaps with this climatology may result in an underestimate of the seasonal pCO_2(sw) cycle in certain locations, especially nearshore (Fig. 5). For comparison we calculate monthly CO₂ flux in our study region from SeaFlux and from RFR-CCS, resulting in the monthly climatologies shown in Fig. 8.

Figure 8 Monthly CO₂ flux per unit area for the nearshore coastal (a) and open-ocean (b) portions of the RFR-CCS domain calculated from the SeaFlux ensemble average (dotted black line; individual products in thin dotted lines) and from RFR-CCS (solid blue line). The grey shaded area represents variability in SeaFlux pCO_2(sw) products calculated as plus and minus 1 standard deviation. Also shown is the spatially distributed CO₂ flux per unit area calculated from RFR-CCS (c) and from SeaFlux (d), as well as the difference between them (e). Red in panels (c) and (d) indicates net release of CO₂ to the atmosphere, whereas blue indicates net uptake of CO₂; in panel (e), red indicates where the RFR-CCS CO₂ flux per unit area is greater and blue where the SeaFlux CO₂ flux per unit area is greater. The solid black line denotes the boundary between the nearshore coastal (a) and open ocean (b) calculated as 100 km from the coast. All calculations are performed using ERA5 winds and an identical gas exchange coefficient (Γ₆₆₀=0.276).

Overall, the SeaFlux ensemble (with ERA5 winds) suggests an oceanic uptake of 69.2 Tg C yr⁻¹ for the RFR-CCS domain between 1998 and 2019 (inclusive) compared to an uptake of 60.0 Tg C yr⁻¹ calculated from RFR-CCS. Of the excess 9.2 Tg C yr⁻¹ uptake from SeaFlux, 5.7 Tg C yr⁻¹ comes from the open ocean and 3.5 Tg C yr⁻¹ from the nearshore coastal ocean (within 100 km of the coast). Given that the nearshore coastal ocean only comprises about 9 % of the RFR-CCS region yet ∼ 38 % of the discrepancy, the discrepancy in coastal uptake is more significant on a per area basis than the open-ocean discrepancy, as can be observed visually in Fig. 8. This discrepancy may reflect more coastal outgassing captured by RFR-CCS than the SeaFlux ensemble, consistent with the annual mean differences shown in Fig. 3a and c. Still, the RFR-CCS results do lie within the variability of SeaFlux pCO_2(sw) products (Fig. 8a and b).

3.8 Effect of sporadic sampling on coastal CO₂ flux calculations

RFR-CCS includes pCO_2(sw) values for the coastal and offshore ocean in the northeastern Pacific that are representative of monthly conditions. However, air–sea carbon dioxide exchange, which is driven by the difference between oceanic and atmospheric pCO₂, operates on shorter timescales. It has been demonstrated in the past that inadequate sampling frequency can be a significant factor biasing CO₂ flux ($F_{{\mathrm{CO}}_{\mathrm{2}}}$) estimates (e.g., Monteiro et al., 2015).

To demonstrate this potential bias, Fig. 9 shows $F_{{\mathrm{CO}}_{\mathrm{2}}}$ at the CCE2 mooring over the course of 2015 (1) calculated from RFR-CCS monthly pCO_2(sw) matched to NOAA marine boundary layer monthly pCO_2(atm) and monthly averages of squared 3-hourly ERA5 winds (blue), (2) calculated from 3-hourly mooring measurements of pCO_2(sw) and pCO_2(atm) matched to squared 3-hourly ERA5 winds (grey), and (3) as the 1-standard-deviation envelope obtained by the following Monte Carlo process: assigning one randomly selected pair of 3-hourly mooring measurements of pCO_2(sw) and pCO_2(atm) from each month as the monthly values, matching them with squared 3-hourly ERA5 winds to calculate $F_{{\mathrm{CO}}_{\mathrm{2}}}$, and repeating this 100 000 times to obtain statistically meaningful values (green).

The $F_{{\mathrm{CO}}_{\mathrm{2}}}$ values provided by the 3-hourly mooring measurements are as close as possible to the true flux. Those provided by RFR-CCS are a best-case scenario for monthly flux approximations in the absence of continuous measurements (because the RFR model was trained on monthly mean pCO_2(sw) values from 3-hourly observations at CCE2). Those provided by the Monte Carlo analysis provided reasonable ranges of $F_{{\mathrm{CO}}_{\mathrm{2}}}$ that might be obtained from sporadic sampling of one measurement per month without the benefit of an advanced interpolation routine like RFR-CCS. The annual $F_{{\mathrm{CO}}_{\mathrm{2}}}$ calculated from RFR-CCS (−0.26 mol C m⁻² yr⁻¹) agrees fairly well with that from the 3-hourly mooring measurements (−0.18 mol C m⁻² yr⁻¹). The smaller uptake from the mooring measurements likely reflects the effect of transient outgassing events in the spring and summer, when positive ΔpCO₂ coincides with high wind speeds. The range from the Monte Carlo analysis (−0.00 to −0.36 mol C m⁻² yr⁻¹) highlights the variety of outcomes in calculated $F_{{\mathrm{CO}}_{\mathrm{2}}}$ that might result from sporadic sampling in the coastal ocean, representative of a region with no high-resolution mooring measurements that may be observed by a ship's underway system only a few times a year.

Figure 9 Hourly flux of CO₂ across the air–sea interface ($F_{{\mathrm{CO}}_{\mathrm{2}}}$) calculated from 3-hourly mooring observations (grey), monthly values from RFR-CCS (blue), and the 1-standard-deviation envelope of a Monte Carlo analysis (n=100 000) whereby one randomly selected 3-hourly mooring observation from each month is selected to represent that month (green). The bar chart on the right gives annual $F_{{\mathrm{CO}}_{\mathrm{2}}}$ based on 3-hourly mooring observations (−0.18 mol C m⁻² yr⁻¹) and RFR-CCS (−0.26 mol C m⁻² yr⁻¹), along with the uncertainty in annual $F_{{\mathrm{CO}}_{\mathrm{2}}}$ from mooring measurements based on the Monte Carlo analysis (−0.00 to −0.36 mol C m⁻² yr⁻¹).

In large portions of the open ocean, low temporal variability and high spatial correlation mean that the aliasing problem may be a relatively low-priority concern for calculations of $F_{{\mathrm{CO}}_{\mathrm{2}}}$ from sporadic pCO_2(sw) measurements (Bushinsky et al., 2019). However, the dynamic coastal ocean is dominated by processes that influence pCO_2(sw) and $F_{{\mathrm{CO}}_{\mathrm{2}}}$ on short spatial and temporal scales, making observational frequency a significant factor that can bias annual $F_{{\mathrm{CO}}_{\mathrm{2}}}$ calculations. This bolsters the case for the expansion and enhancement of coastal carbon observing systems even with pCO_2(sw) gap-filling methods, such as the one described here, at our disposal.

4 Data availability

The RFR-CCS data product (Sharp et al., 2022) is available as a NetCDF and MATLAB file at https://doi.org/10.5281/zenodo.5523389.

5 Code availability

MATLAB code used to process data and create figures included in this paper is provided at https://github.com/jonathansharp/RFR-CCS (last access: 1 April 2022) and https://doi.org/10.5281/zenodo.6484875 (Sharp, 2022). The majority of this code is also compatible with the open-source software GNU Octave.

6 Conclusions

This work presents a data product, called RFR-CCS, of surface ocean pCO₂ in the California Current System and surrounding ocean regions. RFR-CCS was constructed from pCO_2(sw) observations in the Surface Ocean CO₂ Atlas version 2021 (Bakker et al., 2016), which were related to predictor variables (Table 1) using a random forest regression approach. Validation exercises (Table 3) reveal that this approach is able to predict independent pCO_2(sw) values with a skill commensurate with expectations (mean bias near zero and RMSE ≈ 30 µatm), considering the highly variable coastal ocean comprises a large portion of the study region.

RFR-CCS captures variability in pCO_2(sw) in the northeastern Pacific, especially at coastal time series locations, more effectively than global-scale data products of pCO_2(sw). This is evident through comparisons to gridded monthly observations in the SOCAT database (Fig. 4), to monthly observations at fixed mooring sites (Figs. 2 and B3; Table B1), and to seasonal amplitudes of pCO_2(sw) measured at those moorings (Fig. 5; Table 4). The improvements made by RFR-CCS mainly represent the enhanced ability of regional data fits to capture local-scale variability compared to global data fits. Going forward, perhaps global-scale gap-filled pCO_2(sw) products that include a clustering step would benefit from the creation of a greater number of clusters in the coastal ocean, allowing for more robust reconstruction of local variability. Improvements detailed here may also be due to the flexibility of RFR in capturing multiple different length scales of variability (Gregor et al., 2017), which may make the method especially useful for regions that span both the coastal and open ocean. The CCS is also particularly data-rich, and this work demonstrates the excellent resolution of nearshore variability that can be achieved in gap-filled pCO_2(sw) products when coastal observing systems are sustained over time. Examination of the spatiotemporal distribution of pCO_2(sw) observations contained in the SOCAT database (Bakker et al., 2016; http://www.socat.info, last access: 7 March 2022) suggests that analyses similar to this one could be effective for other coastal regions around North America, the western North Pacific, and the eastern North Atlantic. However, different predictor variable–pCO_2(sw) relationships are likely to exist in these distinct ocean regions since each has a unique physical and biogeochemical setting.

Spatial and seasonal patterns of pCO_2(sw) revealed by RFR-CCS reflect interactions of physical and biological processes that differ substantially with latitude, season, and distance from shore (Figs. 6 and 7). For example, high annual mean pCO_2(sw) in a narrow band of the central coastal CCS reflects spring and summer upwelling; low annual mean pCO_2(sw) and CO₂ uptake in the northern coastal CCS and the offshore CCS in general reflects CO₂ drawdown by primary production, largely stimulated by nutrients delivered by coastal upwelling. Generally, across the study region, interpretations of pCO_2(sw) variability and the processes that drive it coincide with local-scale explanations in the coastal environment, suggesting high heterogeneity in coastal carbon cycling.

Finally, in the context of sea surface pCO₂ gap-filling strategies, this study highlights important factors that should be considered when working in coastal areas or regions that span the coastal to open-ocean continuum. For one, although a global gap-filled product may demonstrate mean annual values and average seasonal amplitudes of pCO_2(sw) that represent a broad region effectively, this does not mean that local-scale variability within that region has been captured just as well. Data-rich regions like the CCS confirm this notion, especially when variability at fixed time series sites like moored autonomous platforms is considered. Misrepresentation errors of this nature are especially concerning in dynamic nearshore environments, where local-scale processes can result in surface biogeochemical characteristics that change rapidly over short timescales. These rapid changes can have direct consequences for local biological responses and for CO₂ flux, both of which operate on relatively short timescales. To address potential errors associated with misrepresentation of pCO_2(sw) variability, the spatiotemporal coverage of carbon observing systems must be improved, especially at ocean margins. Further, innovative implementation and assessment of machine-learning approaches (Gregor et al., 2017; Gloege et al., 2021), biogeochemical models (DeVries et al., 2019; Friedlingstein et al., 2020), and ensemble approaches (Lebehot et al., 2019; Fay et al., 2021) should continue to be explored to best leverage the existing data.

Appendix A: Processing of predictor variables

SST (Fig. A1) was obtained from the NOAA daily Optimum Interpolation Sea Surface Temperature (OISST) analysis product (Reynolds et al., 2007; Huang et al., 2021). This data product combines satellite and in situ observations of SST using an optimum interpolation (OI) technique, providing daily SST values at 0.25^∘ resolution. We averaged daily gridded SST values from OISSTv2.1 for each month from 1998 to 2020 to obtain the required monthly 0.25^∘ resolution datasets to match with our gridded pCO₂.

Sea surface salinity (SSS) was obtained from the NASA Estimating the Circulation and Climate of the Ocean (ECCO) project. The ECCO2 state estimate (Menemenlis et al., 2008) uses a Green's function approach (Menemenlis et al., 2005) to make optimal adjustments to parameters, initial conditions, and boundary conditions of a general circulation model to produce a daily ocean state estimate. We averaged daily gridded SSS values from the ECCO2 state estimate for each month from 1998 to 2020 to obtain the required monthly 0.25^∘ resolution datasets to match with our gridded pCO₂.

Sea surface chlorophyll a concentration estimates (Chl; Fig. A2), based on Sea-Viewing Wide Field-of-View Sensor (SeaWiFS) and Moderate Resolution Imaging Spectroradiometer (MODIS) satellite data, were obtained from the Oregon State University (OSU) Ocean Productivity website (http://www.science.oregonstate.edu/ocean.productivity, last access: 16 June 2021). The OSU Ocean Productivity website provides both monthly and 8 d Chl files at either $\mathrm{1}/\mathrm{6}$^∘ or $\mathrm{1}/\mathrm{12}$^∘ resolution. We obtained monthly $\mathrm{1}/\mathrm{6}$^∘ resolution files for 1998–2002 (SeaWiFS-based) and 2003–2020 (MODIS-based) and interpolated each to a 0.25^∘ resolution grid using a standard two-dimensional linear interpolation for each monthly file. For high-latitude wintertime gaps in the Chl datasets, we interpolated Chl for each grid cell through time using one-dimensional linear interpolation when observations in the previous and subsequent month were available. To avoid anomalous values at the beginning and end of the time series, empty grid cells were filled with nearest-neighbor interpolation when a previous or subsequent observation was not available (Fig. A3). Chl was log₁₀-transformed to produce a distribution of values that was closer to normal before constructing the regression model.

Wind speed data (Fig. A4) were obtained from the ERA5 reanalysis product (Hersbach et al., 2020), produced by the European Centre for Medium-Range Weather Forecasts (ECMWF). The ERA5 atmospheric reanalysis provides a detailed record of atmospheric parameters from 1950 to the present day. We obtained monthly, 0.25^∘ resolution wind speed data at 10 m above the surface from the Copernicus Climate Change Service (C3S) Climate Data Store (CDS). Wind speed (U) was calculated from its vector components (north–south wind, v_w, and east–west wind, u_w).

$$\begin{array}{}\text{(A1)}& U=\sqrt{v_{\mathrm{w}}^{\mathrm{2}}+u_{\mathrm{w}}^{\mathrm{2}}}\end{array}$$

Atmospheric CO₂ partial pressure (pCO_2(atm)) was obtained from the NOAA marine boundary layer (MBL) reference (Dlugokencky et al., 2020). This data product is derived from weekly air samples of atmospheric CO₂ mole fraction (xCO₂) at a subset of sites from the NOAA Cooperative Global Air Sampling Network. The product is provided as weekly latitudinal averages with a resolution of sin(lat) = 0.5. We interpolated weekly xCO₂ values to monthly xCO₂ values relative to the middle of each month. To convert xCO₂ to pCO_2(atm), xCO₂ was multiplied by monthly sea level pressure (P) from NCEP reanalysis, which was corrected for water vapor pressure (VP${}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$) as described by Dickson et al. (2007).

$$\begin{array}{}\text{(A2)}& p{\mathrm{CO}}_{\mathrm{2}\left(\mathrm{sw}\right)}=x{\mathrm{CO}}_{\mathrm{2}}\left[P-{\mathrm{VP}}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}\right]\end{array}$$

Mixed layer depths (MLDs), based on output from the Hybrid Coordinate Ocean Model (HYCOM) (Chassignet et al., 2007), were obtained from the OSU Ocean Productivity website. We obtained monthly $\mathrm{1}/\mathrm{6}$^∘ resolution MLD files and interpolated each to a 0.25^∘ resolution grid using a standard two-dimensional linear interpolation for each monthly file. MLD was log₁₀-transformed to produce a distribution of values that was closer to normal before constructing the regression model.

Figure A1 Gridded means of SST from satellite observations from 1998–2020.

Figure A2 Gridded means of chlorophyll a concentration from satellite observations from 1998–2020.

Distance from shore (Dist.) for each grid cell was calculated using the dist2coast.m function from the Climate Data Toolbox for MATLAB (Greene et al., 2019), applied to each latitude–longitude grid cell. That function accepts input of latitude and longitude coordinates and returns the great circle distance to the nearest coastline.

Year (yr) was normalized to an epoch of 1997 (i.e., yr_norm = yr – 1997). Month of year (mn) was transformed into two separate predictor variables (mn_sin and mn_cos) using sine and cosine functions to maintain its cyclical nature (after Gregor et al., 2018).

$$\begin{array}{}\text{(A3)}& {\displaystyle }{\mathrm{mn}}_{\sin }& {\displaystyle }=\sin (\mathrm{2}\mathit{\pi }\cdot \mathrm{mn}/\mathrm{12})\text{(A4)}& {\displaystyle }{\mathrm{mn}}_{\cos }& {\displaystyle }=\cos (\mathrm{2}\mathit{\pi }\cdot \mathrm{mn}/\mathrm{12})\end{array}$$

Figure A3 Sea surface chlorophyll concentration at 55^∘ N, 135^∘ W. Data from satellite observations are in orange and interpolated data are in blue.

Figure A4 Gridded means of wind speed from ERA5 reanalysis from 1998–2020.

Appendix B: Supplementary figures and tables

Figure B1 The number of years containing a pCO_2(sw) observation within each month over the 23 years of our gridded pCO_2(sw) data product from 1998–2020.

Figure B2 Predictor variable feature importances calculated for the random forest regression model fit used to produce RFR-CCS (Sharp et al., 2022; https://doi.org/10.5281/zenodo.5523389).

Figure B3 Like Fig. 2 in the main text, showing monthly values of pCO_2(sw) from mooring observations (black), RFR-CCS (blue), the mooring-excluded RFR-CCS-Eval model (orange), and L17 (green). The envelopes around the black lines equal the standard deviations of all mooring observations within each month, representing the natural variability of the 3-hourly mooring measurements; the envelopes around the blue and orange lines represent the RFR-CCS and RFR-CCS-Eval results plus 1 standard uncertainty (43.6 µatm; Sect. 3.5); the envelopes around the green lines represents the L17 data product plus the RMSE of an independent data evaluation in the province most closely associated with the mooring locations (52.5 µatm; Table 3 of Laruelle et al., 2017; Province P7).

Figure B4 Monthly mean differences in pCO_2(sw) values between RFR-CCS-clim and L20 (top) and between RFR-CCS-clim and L17 (bottom).

Figure B5 Correlations (a, b) and the p values of those correlations (c, d) in each grid cell of RFR-CCS between pCO_2(sw) and SST (a, c) and between pCO_2(sw) and wind speed (b, d).

Table B1 Mean biases (MBs), root mean squared errors (RMSEs), and coefficients of determination (R²) for comparisons of RFR-CCS, the mooring-excluded RFR-CCS-Eval, and L17 to mooring observations.