A machine-learning reconstruction of sea surface pCO₂ in the North American Atlantic Coastal Ocean Margin from 1993 to 2021

Insufficient spatiotemporal coverage of observations of the surface partial pressure of CO₂ (pCO₂) has hindered precise carbon cycle studies in coastal oceans and justifies the development of spatially and temporally continuous pCO₂ data products. Earlier pCO₂ products have difficulties in capturing the heterogeneity of regional variations and decadal trends of pCO₂ in the North American Atlantic Coastal Ocean Margin (NAACOM). This study developed a regional reconstructed pCO₂ product for the NAACOM (Reconstructed Coastal Acidification Database-pCO₂, or ReCAD-NAACOM-pCO₂) using a two-step approach combining random forest regression and linear regression. The product provides monthly pCO₂ data at 0.25° spatial resolution from 1993 to 2021, enabling investigation of regional spatial differences, seasonal cycles, and decadal changes in pCO₂. The observation-based reconstruction was trained using Surface Ocean CO₂ Atlas (SOCAT) observations as observational values, with various satellite-derived and reanalysis environmental variables known to control sea surface pCO₂ as model inputs. The product shows high accuracy during the model training, validation, and independent test phases, demonstrating robustness and a capability to accurately reconstruct pCO₂ in regions or periods lacking direct observational data. Compared with all the observation samples from SOCAT, the pCO₂ product yields a determination coefficient of 0.92, a root-mean-square error of 12.70 µatm, and an accumulative uncertainty of 23.25 µatm. The ReCAD-NAACOM-pCO₂ product demonstrates its capability to resolve seasonal cycles, regional-scale variations, and decadal trends of pCO₂ along the NAACOM. This new product provides reliable pCO₂ data for more precise studies of coastal carbon dynamics in the NAACOM region. The dataset is publicly accessible at https://doi.org/10.5281/zenodo.14038561 (Wu et al., 2024a) and will be updated regularly.

1 Introduction

Accurate and comprehensive datasets of the sea surface partial pressure of CO₂ (pCO₂) are necessary for quantifying coastal CO₂ uptake and assessing the impact of climate change on coastal ocean ecosystems. On a global scale, the coastal ocean, covering 8.4 % (30.4 × 10⁶ km²) of the global ocean surface area (Chen et al., 2013; Dai et al., 2022), plays a significant role in the global carbon budget, accounting for approximately 10.9 % of the global ocean CO₂ uptake from the atmosphere (0.25 of 2.3 Pg C yr⁻¹) on the global average (Dai et al., 2022; Friedlingstein et al., 2023). However, on regional scales, areal-based CO₂ uptake in specific coastal regions are often much greater than those in open oceans despite their less distinguishable global means (Dai et al., 2022). This is because sea surface pCO₂ is highly variable due to the influence of various physical and biogeochemical processes in coastal oceans, such as riverine input, upwelling, tidal mixing, and large-scale circulations (Laruelle et al., 2018; Roobaert et al., 2024b). Thus, accurately quantifying the CO₂ uptake in specific coastal regions becomes particularly challenging when only using observations due to the incomplete coverage of pCO₂ data in space and time.

This study focuses on the North American Atlantic Coastal Ocean Margin (NAACOM; Fig. 1). The entire region is defined as the area within 400 km of the coastline and is divided into six subregions following Fennel et al. (2019) based on their geographic locations: the Gulf of Mexico (GoMx), South Atlantic Bight (SAB), Mid-Atlantic Bight (MAB), Gulf of Maine (GoMe), Scotian Shelf (SS), and Gulf of St. Lawrence and Grand Banks (GStL&GB). The carbonate system in the NAACOM is influenced by large-scale circulations (Fig. 1), including the Gulf Stream and Labrador Current, as well as local processes like river discharge, export from marshes, and upwelling dynamics (Cai et al., 2020; Fennel et al., 2019; Wang et al., 2013). These complex physical and biogeochemical processes contribute to substantial spatial and temporal heterogeneity in sea surface pCO₂ across the NAACOM (Cai et al., 2020). Elucidating the driving mechanisms of these spatiotemporal pCO₂ variations necessitates extensive data coverage in time and space in this region. Over the past 2 decades, coastal field investigation efforts in this region have substantially increased through programs like the East Coast Ocean Acidification (ECOA) and Gulf of Mexico Ecosystems and Carbon Cruise (GOMECC) (Cai et al., 2020; Wang et al., 2013; Wanninkhof et al., 2015). Ongoing measurements from these cruises, combined with ongoing measurements from volunteer observing ships and buoys, are quality-controlled and compiled in the Surface Ocean CO₂ Atlas (SOCAT) database (Bakker et al., 2016), substantially advancing our understanding of coastal inorganic carbon chemistry along the NAACOM (Cai et al., 2020).

Despite significant progress in observational efforts, the spatial and temporal coverage of pCO₂ data remains limited in the NAACOM, with observations encompassing only 2.9 % of grid cells during the period 1993–2021 (Fig. 2). Observations are concentrated in the southern region, with fewer samples available during winter. This data scarcity introduces substantial uncertainty into air–sea CO₂ exchange quantification and hinders comprehensive understanding of coastal inorganic carbon dynamics, particularly in areas north of Cape Cod where measurements are very sparse (Fig. 2). For example, reported air–sea CO₂ fluxes for the GoMe exhibit a wide range spanning from −0.50 to +2.50 mol C m⁻² yr⁻¹, with conflicting reports characterizing it as a CO₂ source (Fennel and Wilkin, 2009; Vandemark et al., 2011), CO₂-neutral (Signorini et al., 2013), or a CO₂ sink (Cahill et al., 2016; Rutherford et al., 2021), underscoring the need for improved pCO₂ data coverage.

Figure 1 Topography (m) and large-scale circulation along the North American Atlantic Coastal Ocean Margin (NAACOM). The study region, defined as coastal areas extending 400 km offshore, is indicated by blue shading. The thick black line is the 200 m isobath, which roughly marks the shelf break and typically defines the continental shelf boundary. The Gulf Stream (thick red dashed line with an arrow) flows northward along the eastern coast of the United States before veering eastward into the open Atlantic Ocean around Cape Hatteras. The Labrador Current (thick light-blue dashed line with an arrow) flows southward along the eastern coast of Canada before meeting the Gulf Stream. Following Fennel et al. (2019), the study region is divided into six subregions by straight orange lines: the Gulf of Mexico (GoMx), South Atlantic Bight (SAB), Mid-Atlantic Bight (MAB), Gulf of Maine (GoMe), Scotian Shelf (SS), and Gulf of St. Lawrence and Grand Banks (GStL&GB). Dashed contour lines indicate bathymetric depths of 50 and 100 m on the shelf (from the coastline to the 200 m isobath) and 1000, 2000, 3000, and 4000 m from the shelf break to the open ocean.

Recently, various pCO₂ products, global or regional, with full coverage in time and space were developed as essential supplements to observations. These products usually employed diverse algorithms, environmental proxies from satellites and reanalysis products as model inputs, and SOCAT observations as constraints to reconstruct the pCO₂ field with full temporal and spatial coverage. The development of those products has significantly advanced our understanding of inorganic carbon chemistry and the ocean carbon cycle. For example, seven global pCO₂ products were used to evaluate the ocean CO₂ uptake in the Global Carbon Budget 2023 edition (Friedlingstein et al., 2023). However, most of these products reconstruct pCO₂ in the open ocean, with coastal regions often being extrapolated or excluded. Currently, only one pCO₂ product has been developed specifically for the coastal ocean on a global scale (Laruelle et al., 2017; Roobaert et al., 2024a). This product was recently combined with an open-ocean product to create a global reconstruction of the ocean CO₂ sink (Landschützer et al., 2020) and has since been utilized to narrow the variability in global reconstructions (Fay et al., 2021). However, global products primarily aim to ensure high accuracy of parameters on a global average; they may not guarantee equivalent accuracy for spatiotemporal variations on the regional scale. In comparison, regional pCO₂ products have demonstrated superior capability in resolving detailed small-scale variations.

Within the NAACOM region, several area-specific pCO₂ products have been reconstructed, focusing on specific regions such as the GoMx (e.g., Chen and Hu, 2019; Fu et al., 2020; Lohrenz and Cai, 2006) and the SAB and MAB (e.g., Wang et al., 2024; Xu et al., 2020). These regional and global pCO₂ products are valuable for validating model estimations (Roobaert et al., 2022; Ross et al., 2023). However, existing products often have limitations in spatial coverage, temporal resolution, or trend analysis capabilities. For instance, Chen and Hu (2019) provided a high-resolution (4 km) pCO₂ product for the GoMx, but this product faces challenges in capturing decadal changes in pCO₂ (Wu et al., 2024b). Conversely, Xu et al. (2020) successfully captured decadal trends of pCO₂ but only as area-averaged pCO₂ time series for the SAB and MAB, lacking comprehensive spatial coverage. Signorini et al. (2013) reconstructed a product using multiple linear regression (MLR) covering the areas from the SAB to SS, but this only spans 8 years (2003–2010). Despite these valuable efforts, there remains a lack of comprehensive data products that adequately capture regional variations, seasonal cycles, and decadal changes in pCO₂ simultaneously for the entire NAACOM.

This study aims to develop a regional pCO₂ product specifically designed for the NAACOM, encompassing coastal regions extending 400 km offshore from the GoMx to the GB (Fig. 1). We integrated random forest and linear regression methods with hydrological parameters from satellite observations and reanalysis data to generate a monthly reconstructed pCO₂ product at 0.25° spatial resolution spanning the period from 1993 to 2021. The pCO₂ product, termed the Reconstructed Coastal Acidification Database or ReCAD-NAACOM-pCO₂, is specifically designed to resolve the spatial variations, seasonal cycles, and decadal changes of pCO₂ along the NAACOM.

The structure of this paper is as follows: Sect. 2 details the methodology used to reconstruct ReCAD-NAACOM-pCO₂ and describes the datasets employed. Section 3 evaluates the accuracy of the reconstructed product, performance, and applicability in resolving seasonal cycles, regional variations, and decadal trends of pCO₂. Sections 4 and 5 provide links to access the dataset and codes used to generate the dataset and figures presented in this study. The final section summarizes the conclusions. ReCAD-NAACOM-pCO₂ demonstrates enhanced capability in resolving spatial variations and capturing the seasonal cycle and decadal trends of pCO₂ compared to the global products across different subregions in the NAACOM. This product offers improved insights into coastal carbon dynamics in this complex region, addressing the need for a comprehensive pCO₂ dataset in the NAACOM. Applications of this data product to examine the processes controlling the spatial variability, seasonal cycle, and decadal trends of pCO₂ and air–sea CO₂ flux will be published separately.

2 Data and methods

2.1 Observational data from SOCAT

The observational data for training the regression model were measurements of the seawater fugacity of CO₂ (fCO₂) extracted from the SOCAT database (2023 edition). fCO₂ represents the pCO₂ corrected for the nonideal behavior of the gas in seawater, and both are commonly used in oceanographic studies. SOCAT compiles quality-controlled fCO₂ measurements from various platforms, including research vessels, commercial ships, and moorings (Bakker et al., 2016). This study used the monthly gridded SOCAT coastal product with a spatial resolution of 0.25° × 0.25° (but with data gaps). The gridded product incorporated measurements with quality flags A and B (uncertainty of 2 µatm) and C and D (uncertainty of 5 µatm) (Bakker et al., 2016). Over the period of 1993–2021, the SOCAT product encompassed 55 347 grid cells within our study area (Fig. 2), accounting for approximately 2.9 % of the total number of grid cells in the NAACOM. The observational data show a lower sampling density in the areas north of Cape Cod and the western and southern GoMx (blue box in Fig. 2). The temporal distribution of the samples exhibits a notable bias, with reduced collection during winter (Fig. 2d). Despite these spatial and temporal heterogeneities, the SOCAT observations provide coverage across all subregions and seasons of the NAACOM (Fig. 2). This comprehensive, albeit sparse, coverage facilitates the reconstruction of the fCO₂ and pCO₂ field through interpolation and regression techniques.

Figure 2 Spatial distribution of sea surface pCO₂ observations from the SOCAT database (2023 edition) in the NAACOM across the four seasons from 1993 to 2021. Grid samples with data were counted by season: (a) spring (March to May), (b) summer (June to August), (c) fall (September to November), and (d) winter (December to February). The study region is divided into northern (blue box) and southern (red box) areas at approximately 41.5° N (Cape Cod). The number and percentage of the grid samples are indicated for each region by season. The color scale represents pCO₂ values (µatm). A higher sampling density is evident in the southern area. Winter shows the lowest overall sampling coverage. Note that the SOCAT database provides quality-controlled fCO₂ measurements as the default parameters, which are subsequently converted into pCO₂ using Eq. (2).

2.2 Model design

The procedures for developing and reconstructing the pCO₂ product are illustrated in Fig. 3. Initially, the input variables and sea surface fCO₂ data were matched to create a comprehensive dataset. To maintain consistency with the SOCAT database, which reports seawater CO₂ concentrations as fCO₂, we adopted fCO₂ as the model training label and first-step output variable in our model. During the model development phase, fCO₂ measurements served as training labels for the machine-learning algorithm. The matched dataset was then divided into two sets: X1, encompassing the periods 1993–2003 and 2006–2021, and X2, covering 2004–2005. Set X1 was randomly subdivided further, with 80 % allocated for model training and the remaining 20 % for the validation test. Set X2 served as an independent test set. The model training set (80 % of X1) was used to develop a two-step RFR–LR (random forest regression–linear regression) model. The RFR model is designed to capture complex, nonlinear relationships between the input variables and the target variable (i.e., fCO₂), while the LR model is subsequently applied to mitigate potential systematic biases in RFR-derived fCO₂ values arising from spatiotemporal heterogeneities in the SOCAT observational dataset (Fig. 2). RFR, an ensemble learning technique, combines multiple decision trees to produce more accurate and stable predictions (Breiman, 2001; Lu et al., 2019). Each decision tree in the RFR model is trained on a randomly selected subset of the input data, with the final prediction derived from the average output of all the trees. This approach mitigates overfitting and enhances the generalization performance of the model, making it particularly suitable for large datasets with complex, nonlinear variable relationships. The RFR model was trained using 10-fold cross-validation with optimized hyperparameters, including a minimum leaf size of 1, a bagging method for ensemble aggregation, and 300 learning cycles after tuning. After RFR model training, the LR model was applied to the RFR-estimated fCO₂ (fCO_2est) output to make sure that the RFR model was not systematically biased:

$$\begin{array}{}\text{(1)}& f{\mathrm{CO}}_{\mathrm{2}\mathrm{obs}}=a\times f{\mathrm{CO}}_{\mathrm{2}\mathrm{est}}+b+\mathit{\epsilon },\end{array}$$

where fCO_2obs is the observed fCO₂ from SOCAT, a is the linear regression coefficient, b is the intercept, and ε is the residual that the linear model cannot resolve. This additional step was implemented to mitigate potential systematic bias in the RFR model that could arise from areas with a higher sampling density, thereby ensuring a more balanced representation across the entire study region. The comparison between the estimated pCO₂ before and after LR calibration is presented in Appendix A. The calibration was applied to each grid cell individually. To increase the data pool for linear regression, samples within a 5 × 5 grid window in space (i.e., 1.25° × 1.25°) were aggregated for LR model development. As the available measurements could not cover every grid cell and were insufficient to produce continuous spatial maps of the calibration coefficients (i.e., a and b in Eq. 1), we employed a locally interpolated regression strategy similar to that of Carter et al. (2018). Mathematically, given the spatial and temporal continuity of fCO_2est and fCO_2obs, the coefficients a and b must also be continuous in space and time. Therefore, we linearly interpolated the coefficients a and b across the NAACOM. The interpolated coefficients were subsequently used to adjust the RFR-derived fCO_2est.

Figure 3 A flowchart of the two-step machine-learning regression model for generating the reconstructed pCO₂ product. The grey boxes represent the input and output datasets. The blue boxes illustrate the model training, validation testing, and independent test processes. The orange boxes represent the final trained model for predicting the reconstructed product. The two models in the orange boxes are identical. The training data, consisting of paired input variables (longitude, latitude, month, sea surface temperature (SST), sea surface salinity (SSS), sea surface height (SSH), atmospheric pCO₂ (pCO_2air), and the corresponding sea surface fCO₂ (fCO_2sea) labels), are divided into two sets: X1 (1993–2003 and 2006–2021) and X2 (2004–2005). X1 is randomly divided further into subsets for the model training set (80 %) and the validation set (20 %). The predictive model combines a random forest regression (RFR) algorithm and a linear regression (LR) algorithm. The trained and validated regression model is then applied to all satellite and reanalysis data (without gaps) to generate the 3D reconstructed fCO_2sea product, which is then converted into pCO_2sea with satellite SST data.

The validation set, comprising 20 % of X1 randomly sampled from 1993–2003 and 2006–2021, serves as a critical monitoring step for model evaluation. This subset plays two key roles: first, it tests hyperparameter tuning by providing independent performance metrics on unseen data, and second, it helps detect potential overfitting by monitoring the divergence between training and validation performance. While the validation set itself cannot prevent overfitting, it enables the detection of overfitting patterns when the performance of the model improves on training data but deteriorates on validation data. Through this continuous evaluation process, the validation set ensures more robust model development and helps achieve better generalization capabilities.

The independent test set (X2), covering the years 2004–2005, serves as a critical evaluation period specifically designed to assess the reliability of the model in predicting values for years that were completely excluded from both the training and validation phases. Because we intentionally withhold these 2 years from model development, this approach directly tests the capability of the model to generate reliable predictions and fill temporal data gaps for periods without observational data.

Finally, the trained model is applied to all the satellite and reanalysis data to generate the final gap-free reconstructed fCO₂ data. As most products reported seawater CO₂ concentrations as pCO₂, we subsequently converted the reconstructed fCO₂ values into pCO₂ using the following equation (Takahashi et al., 2020), and our final product reports both fCO₂ and pCO₂:

$$\begin{array}{}\text{(2)}& p{\mathrm{CO}}_{\mathrm{2}}=f{\mathrm{CO}}_{\mathrm{2}}\times (\mathrm{1.00436}-\mathrm{4.669}\times {\mathrm{10}}^{-\mathrm{5}}\times \mathrm{SST}).\end{array}$$

2.3 Regression model input variables from satellite and reanalysis

The input variables for training the regression model are longitude (long), latitude (lat), month, sea surface temperature (SST), sea surface salinity (SSS), sea surface height (SSH), and atmospheric pCO₂ (pCO_2air). Longitude, latitude, and month serve as spatiotemporal predictors, enabling the algorithm to identify and capture regional and seasonal variability in fCO₂ within the study area (Su et al., 2020; Yang et al., 2024). SST, SSS, and SSH are critical variables that characterize the physical and biogeochemical ocean settings, which play a crucial role in determining the spatial and temporal variability of fCO₂. pCO_2air represents the atmospheric forcing in the air–sea CO₂ exchange. Including pCO_2air is essential for accurately assessing the decadal pCO₂ trend.

SST data were obtained from the National Oceanic and Atmospheric Administration (NOAA) Optimum Interpolation Sea Surface Temperature (OISST) v2.1 product (Huang et al., 2021). The OISST dataset is a global gridded SST analysis that blends observations from various sources, including satellites, ships, and buoys. The dataset employs an optimum interpolation technique to combine these observations and generate a daily SST field at a spatial resolution of 0.25° × 0.25°. For this study, the daily SST data were averaged to create a monthly product.

SSS data were obtained from the Simple Ocean Data Assimilation (SODA) v3.15.2 product (Carton et al., 2018). SODA is a comprehensive reanalysis dataset that integrates a global ocean model with observational data to estimate ocean state variables consistently. The SODA system assimilates observations from multiple sources, including floats, moorings, and ship-based measurements, thereby constraining the model output and enhancing the accuracy of the represented ocean physical properties, including SSS. The SODA v3.15.2 product offers monthly SSS data with a temporal resolution of 1 month and a spatial resolution of 0.5° × 0.5°, which were linearly interpolated to a 0.25° × 0.25° grid resolution to maintain consistency with other input variables and the gridded SOCAT fCO₂ data. Note that such interpolation could potentially introduce additional errors. We doubled the SSS uncertainty in the region, assuming that this would encompass its true uncertainty (see Appendix B).

SSH data were extracted from the Global Ocean Gridded L4 Sea Surface Heights (CMEMS, 2021) created by the Copernicus Marine Environment Monitoring Service (CMEMS). Since 1993 (ongoing) this product has provided daily SSH data derived from altimeters with a spatial resolution of 0.25° × 0.25°. Daily SSH data were averaged to monthly means.

pCO_2air data converted from the mole fraction of CO₂ in dry air (xCO_2air) were downloaded from the NOAA Marine Boundary Layer (MBL) reference product (Lan et al., 2023). The MBL reference product provides weekly zonal average xCO_2air measurements from a global observation network. The xCO_2air data were linearly interpolated to the same spatial and temporal resolution as the other input variables (0.25° × 0.25°, monthly). xCO_2air was converted into pCO_2air with the equation

$$\begin{array}{}\text{(3)}& p{\mathrm{CO}}_{\mathrm{2}\mathrm{air}}={\mathrm{xCO}}_{\mathrm{2}\mathrm{air}}\times (P-p_{\mathrm{w}}),\end{array}$$

where P is the total atmospheric pressure on the sea surface, which was downloaded from the fifth-generation reanalysis (ERA5) of the European Centre for Medium-Range Weather Forecasts (ECMWF) (Hersbach et al., 2019), and p_w is the water vapor pressure, which was calculated using the formula of Weiss and Price (1980), SST from OISST, and SSS from SODA.

2.4 Evaluation of the models

The accuracy of the model outputs was assessed using several statistical metrics, including the coefficient of determination (R²), root-mean-square error (RMSE), mean absolute error (MAE), and mean bias error (MBE). These metrics were calculated for the training and validation set phases as well as for the independent validation set:

$$\begin{array}{}\text{(4)}& {\displaystyle }{R}^{\mathrm{2}}=\mathrm{1}-{\sum }_i^N{\left(y_{\mathrm{obs},i}-y_{\mathrm{est},i}\right)}^{\mathrm{2}}/{\sum }_i^N{\left(y_{\mathrm{obs},i}-\stackrel{\mathrm{‾}}{y_{\mathrm{obs}}}\right)}^{\mathrm{2}},\text{(5)}& {\displaystyle }\mathrm{RMSE}=\sqrt{{\displaystyle \frac{\mathrm{1}}{N}}{\sum }_i^N{\left(y_{\mathrm{obs},i}-y_{\mathrm{est},i}\right)}^{\mathrm{2}}},\text{(6)}& {\displaystyle }\mathrm{MAE}={\displaystyle \frac{\mathrm{1}}{N}}{\sum }_i^N\left|y_{\mathrm{obs},i}-y_{\mathrm{est},i}\right|,\text{(7)}& {\displaystyle }\mathrm{MBE}={\displaystyle \frac{\mathrm{1}}{N}}{\sum }_i^N\left(y_{\mathrm{obs},i}-y_{\mathrm{est},i}\right),\end{array}$$

where i denotes the ith sample, y_obs and $\stackrel{\mathrm{‾}}{y_{\mathrm{obs}}}$ are the observed pCO₂ values from SOCAT and their average, y_est represents the predicted pCO₂ values from the final model, and N is the total number of matched samples.

2.5 Uncertainty of the reconstructed pCO₂

The uncertainty of the estimated pCO₂ in our product for each grid cell was accumulated from four sources of uncertainties: the direct pCO₂ measurement uncertainty from SOCAT (u_obs), gridding uncertainty (u_grid), mapping uncertainty (u_map), and the uncertainty accumulated from the input variables (u_inputs). The first three sources of uncertainty were calculated according to the approach used by earlier reconstructed pCO₂ products (Landschützer et al., 2014; Roobaert et al., 2024a; Sharp et al., 2022). u_obs is inherited from the SOCAT observations. The SOCAT database uses discrete samples with quality flags A and B (accuracy <2 µatm) and C and D (accuracy <5 µatm) to create the gridded file. Adopting a conservative approach, we used the maximum u_obs of 5 µatm. u_grid was calculated as the standard deviation of the samples used to calculate the gridded fCO₂ in each grid cell. u_map is introduced by reconstructing the pCO₂ using the RFR–LR model. It was evaluated as the RMSE between the reconstructed pCO₂ and the observed pCO₂ values following Roobaert et al. (2024a) and Sharp et al. (2022). Given that the derivation of u_obs, u_grid, and u_map is contingent upon SOCAT observations, these three uncertainties and the total uncertainty $u_{p{\mathrm{CO}}_{\mathrm{2}}}$ are reported on a subregional basis.

In addition to these three sources of uncertainty, this study incorporated cumulative uncertainties from input variables (u_inputs), including SST, SSS, SSH, and pCO_2air. These satellite-derived or reanalysis-based variables inherently possess uncertainties that propagate nonlinearly through the regression model, ultimately affecting the estimated pCO₂ values (Wang et al., 2021, 2023). We employed a Monte Carlo simulation to calculate u_inputs. For each input variable (SST, SSS, SSH, and pCO_2air), we added white noise following a normal distribution $N(\mathrm{0},u_{x_i})$, where $u_{x_i}$ is the uncertainty of the respective input variable x_i. We then recalculated pCO₂ using these noise-added inputs and determined the resulting changes in pCO₂. This process was repeated 100 times for each input variable, and the resulting uncertainty in pCO₂ from each variable was calculated as the standard deviation of the differences between the original reconstructed pCO₂ and the pCO₂ values after adding noise to each grid cell. The final u_inputs was computed as the square root of the quadratic sum of these individual uncertainties from the four input variables. Detailed procedures for determining u_inputs are described in Appendix B.

Assuming that these sources are independent, the uncertainty of the estimated gridded pCO₂ in our product, $u_{p{\mathrm{CO}}_{\mathrm{2}}}$, was calculated using the error propagation (Hughes and Hase, 2010; Taylor, 1997):

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

2.6 Comparison with the global reconstructed pCO₂ product

The ReCAD-NAACOM-pCO₂ product was evaluated through comparisons with seven reconstructed pCO₂ products developed for the global ocean and used in the Global Carbon Budget 2023 edition (Friedlingstein et al., 2023) and one reconstructed pCO₂ product specifically developed for the global coastal ocean (ULB_SOMFNN_coastal_v2; Roobaert et al., 2024a). These data products reconstructed pCO_2sea data using different machine-learning algorithms. Detailed information on the products is summarized in Table 1.

Table 1 References for the global pCO₂ products used for comparison with ReCAD-NAACOM-pCO₂ in this study. The abbreviations in the Methods column are RFRE for random-forest-based regression ensemble, SOM-FFN for self-organizing map–feed-forward network, MLR for multiple linear regression, FFNN for feed-forward neural network, XGB for the eXtreme Gradient Boosting algorithm, and GRaCER for geospatial random-cluster ensemble regression.

3 Results and discussion

3.1 Evaluating the regression model performance

Our product employs a two-step RFR–LR algorithm to retrieve pCO₂. The initial RFR step accurately captures most seasonal and decadal pCO₂ variations across all six subregions (Appendix A). When comparing only at matching grid cells where SOCAT measurements are available, the differences (N=12) in the monthly mean climatology between the SOCAT- and RFR-derived pCO₂ are less than 2 µatm on average, with standard deviations below 5 µatm across all the subregions (Fig. A1). However, the RFR-derived pCO₂ shows lower accuracy in capturing long-term pCO₂ changes in the GoMe and SAB. The subsequent LR calibration improves the performance significantly: R² values increase from 0.69 to 0.81 in the GoMe and from 0.83 to 0.93 in the SAB, while the RMSE decreases from 12.43 to 10.51 µatm in the GoMe and from 10.83 to 8.12 µatm in the SAB (Fig. A2).

The ReCAD-NAACOM-pCO₂ product demonstrated robust performance and high accuracy in capturing pCO₂ variability across the NAACOM (Fig. 4). During the model training phase, the product achieved an R² of 0.96, an RMSE of 9.1 µatm, an MAE of 5.92 µatm, and an MBE of 0.05 µatm (Fig. 4a). The model demonstrated comparable performance metrics during the validation phase (Fig. 4b). To further evaluate the generalizability and robustness of the model, we also conducted an independent test using data from 2004 to 2005 in which not all the data samples were included in the model training and validation sets. During this independent test phase, the pCO₂ product maintained high accuracy, with R²=0.64, RMSE = 27.2 µatm, MAE = 18.86 µatm, and MBE = 0.07 µatm (Fig. 4c). Additionally, most independent validation samples were distributed around the 1:1 correspondence line, proving the ability of the models to predict pCO₂ across unsampled spatial and temporal domains without overfitting. The model consistently demonstrated strong performance during the training, validation, and independent test phases across all the subregions (Table 2). Overall, compared with all the available samples in SOCAT, it achieved an R² of 0.92, an RMSE of 12.70 µatm, an MAE of 7.55 µatm, and an MBE of 0.13 µatm for the entire NAACOM (Table 2), highlighting the generalizability of the ReCAD-NAACOM-pCO₂ product and robustness in effectively capturing the variability in pCO₂ and providing reliable predictions of pCO₂ across the studied regions.

Figure 4 Evaluation of the regression model for reconstructing the ReCAD-NAACOM-pCO₂ product. The density scatterplots compare the product-estimated pCO₂ (pCO_2est) with the in situ SOCAT observations (pCO_2obs) during the (a) model training phase (80 % of the samples during the periods 1993–2003 and 2006–2021), (b) validation phase (20 % of the samples during the periods 1993–2003 and 2006–2021), and (c) independent test phase (samples during the period 2004–2005). The statistical metrics include the coefficient of determination (R²), root-mean-square error (RMSE), mean absolute error (MAE), mean bias error (MBE), and number of samples (N). The color bar represents the number of data points in each bin.

Table 2 Performance of the regression model during the model training, validation, and independent test phases across the different subregions. The metrics include the coefficient of determination (R²), root-mean-square error (RMSE), mean absolute error (MAE), and mean bias error (MBE). The subregions are the Gulf of Mexico (GoMx), South Atlantic Bight (SAB), Mid-Atlantic Bight (MAB), Gulf of Maine (GoMe), Scotian Shelf (SS), and Gulf of St. Lawrence and Grand Banks (GStL&GB).

3.2 Spatial distribution of the product bias

The ReCAD-NAACOM-pCO₂ product exhibited a negligible area-mean bias of +0.13 µatm with a standard deviation of 12.70 µatm when compared to all SOCAT observation grid cells across the entire NAACOM (Fig. 5 and Table 2). This small average difference suggests no consistent overestimation or underestimation by the regression model, indicating the reliability of the product in estimating the monthly and annual mean climatology of pCO₂ across the entire NAACOM region.

Figure 5 Spatial distribution of the MBE between the ReCAD-NAACOM-pCO₂ product and SOCAT observations across the NAACOM. The MBE is calculated for each grid cell as the average difference between product estimates and SOCAT observations. Positive values (red) indicate product overestimation, while negative values (blue) indicate underestimation relative to SOCAT. Regional MBE values with 1 standard deviation are shown for each subregion, corresponding to the values in the last column of Table 2. The overall bias error for the NAACOM is +0.13 ± 12.97 µatm. Following Fennel et al. (2019), the study region is divided into six subregions using straight orange lines: the GoMx, SAB, MAB, GoMe, SS, and GStL&GB. The thick black line is the 200 m isobath, which roughly marks the shelf break and typically defines the continental shelf boundary.

While the area-averaged difference is small, the differences are distributed heterogeneously in space. Larger differences (absolute difference >10 µatm) tend to occur in nearshore regions, particularly along the coastlines of the GoMx and SAB, as well as in northern areas such as the GoMe, SS, and GStL&GB (Fig. 5). These regional variations can be attributed to complex coastal processes such as terrestrial inputs, sparse observations in the northern areas (Lavoie et al., 2021; Rutherford et al., 2021; Salisbury and Jönsson, 2018), and less accurate satellite observations in the nearshore regions (Song et al., 2023). Conversely, smaller differences (absolute difference <2.5 µatm) are observed in the central parts of the GoMx, offshore regions of the SAB and MAB, and some nearshore regions of the SS and GB, which is likely due to more stable oceanic conditions in those regions. The regional MBEs for different machine-learning development phases (training, validation, and test sets) are detailed in Table 2. Despite these regional differences, the MBEs of both the validation set (−1.0 to 1.0 µatm) and the independent test set (−4.5 to 7.5 µatm) demonstrate minimal values across the subregions (Table 2), underscoring the effectiveness of the product in capturing the broader pCO₂ patterns across the NAACOM.

3.3 Evaluating the capacity of the product to capture pCO₂ seasonality

One of the primary objectives of this product is to capture the seasonal cycle of pCO₂ across the NAACOM region. Figure 6 showcases the applicability of the product in capturing the pCO₂ seasonal cycles across the southern and northern areas of the NAACOM (red and blue boxes in Fig. 2). The comparison of monthly climatologies between the gap-filled product and SOCAT observations reveals strong agreement in the southern region despite the coverage difference, with the product-estimated monthly means being only 3.05 ± 5.60 µatm higher than those of SOCAT (Fig. 6a) and suggesting that our product effectively captures the seasonal cycle where data are abundant.

Figure 6 Monthly mean climatology of pCO₂ in the southern and northern areas of the NAACOM from 1993 to 2021. The subregions are (a) the southern areas with the red box in Fig. 2 and (b) the northern areas with the blue box in Fig. 2. Two data representations are shown: (1) SOCAT observations (black curves), which may be influenced by missing data, and (2) the complete gap-filled product output (red curves). The error bars denote 1 standard deviation of the monthly mean climatology of pCO₂. The numbers indicate the mean difference (± 1 standard deviation) in the monthly climatological pCO₂ calculated from the two sources, with positive values indicating higher product estimates compared to SOCAT observations. The x axis shows the months (1–12, where 1 represents January), and the y axis shows pCO₂ (µatm).

In the northern region where SOCAT data are sparse, the gap-filling ability of the product is also demonstrated well. In the northern region, the area-averaged monthly pCO₂ climatology calculated from the continuous reconstructed product is 22 ± 11.12 µatm lower than the SOCAT observations, which can be attributed to the limited observational coverage in this area. This area is characterized by sparse sampling, with the observational density approximately 50 % lower than in the southern region (Fig. 2) due to the smaller area and limited cruise coverage. For instance, the GStL region only has one summer cruise in the SOCAT database (Fig. 2b), and the SS and GoMe have particularly sparse winter observations (Fig. 2d). The higher latitudes typically exhibit larger seasonal amplitudes in pCO₂, making the limited sampling from SOCAT particularly problematic for accurate characterization. Our gap-free product provides comprehensive spatial and temporal coverage, enabling more robust analysis of pCO₂ patterns and variability in these historically undersampled regions.

Over the 29-year period, the product predicts smaller monthly standard deviations in the southern region (less than 40 µatm; error bars in Fig. 6a), suggesting higher model accuracy and less interannual variability in these areas. Conversely, larger monthly standard deviations are observed in the northern areas, suggesting potentially lower accuracy and remarkable interannual variability. However, the larger interannual variability in these areas may be an artifact of the limited observational data available for regression model training, resulting in greater uncertainty in the predictions. Despite differences in the mean monthly climatology, the similar seasonal pCO₂ cycles calculated from SOCAT and the reconstructed product demonstrate the ability of the ReCAD-NAACOM-pCO₂ product to represent seasonal pCO₂ variability across diverse coastal environments. Nevertheless, there exist larger differences between the observations and reconstructed pCO₂ in some months and regions (Fig. 6b), highlighting the importance of the gap-free product in an unbiased understanding of regional carbon cycles (Ren et al., 2024). Detailed sea surface pCO₂ seasonal cycles and their controlling mechanisms across different subregions of the NAACOM will be presented in our subsequent work.

3.4 Evaluating the ability of the products to capture regional variation by comparing them to global products

The ReCAD-NAACOM-pCO₂ product demonstrates the capability to resolve fine-scale regional spatial distributions of pCO₂. Figure 7 illustrates the spatial distribution of the annual mean climatology of pCO₂ across the NAACOM as observed by SOCAT and predicted by different global open- and coastal-ocean pCO₂ products. Despite being affected by missing data, SOCAT observations (Fig. 7a) reveal significant regional variations in pCO₂. In the Louisiana Shelf (LAS) estuary plume region (box 1 in Fig. 7), pCO₂ values consistently remain below 340 µatm, while the West Florida Shelf (WFS; box 2 in Fig. 7) exhibits elevated values exceeding 400 µatm. These contrasting patterns have been reported in previous regional studies (Kealoha et al., 2020; Robbins et al., 2018; Wu et al., 2024b).

Figure 7 Spatial distribution of the annual mean pCO₂ climatology in the NAACOM from different sources. (a) SOCAT observations, (b) the ReCAD-NAACOM-pCO₂ product, (c) the ensemble mean of the seven global open-ocean pCO₂ products listed in Table 1, and (d) the coastal pCO₂ product ULB_SOMFFN_coastal_v2 (Roobaert et al., 2024a). The black contour delineates the coastal-ocean margin. The three boxes represent subregions in the NAACOM: box 1 for the Louisiana Shelf (LAS), box 2 for the West Florida Shelf (WFS), box 3 for the entire northern region, and box 4 for the southern GStL (S.GStL). Mean pCO₂ ± standard deviations of all the grid cells are provided for each dataset. The color scale represents pCO₂ (µatm).

The ReCAD-NAACOM-pCO₂ product demonstrates superior alignment with SOCAT observations in capturing those regional features that have been reported in previous observation-based studies (Fig. 7b), accurately representing the low pCO₂ values in the LAS Mississippi River plume (box 1) and the elevated pCO₂ levels in the WFS (box 2). In contrast, the global reconstructions of pCO₂, represented by the ensemble of seven open-ocean pCO₂ products (Fig. 7c), face challenges in resolving these regional pCO₂ variations, as previously discussed in Wu et al. (2024b). The coastal pCO₂ product from Roobaert et al. (2024a; ULB_SOMFFN_coastal_v2) also captures some small-scale structures like the low pCO₂ in the LAS (Fig. 7d), but the ReCAD-NAACOM-pCO₂ product exhibits values that are closer to the observations. In the northern region (box 3), the ReCAD-NAACOM-pCO₂ product predicts higher pCO₂ levels that are closer to observations in the nearshore region (Fig. 7b). This is not surprising, as ULB_SOMFNN_coastal_v2 is a global product known for its high accuracy on the global average.

In addition to these previously documented regional variations, our product reveals several notable features not previously captured by observations or other existing products. For instance, the GoMe displays intermediate pCO₂ levels of around 380 µatm, which is distinctly higher than surrounding waters at comparable latitudes, a feature previously documented by a pCO₂ product reconstructed using multiple linear regression (Signorini et al., 2013) and 5-year (2004–2009) mooring and cruise data (Vandemark et al., 2011). However, this contradicts two other studies based on numerical models (Cahill et al., 2016; Rutherford et al., 2021). In the southern GStL (S.GStL; box 4 in Fig. 7), pCO₂ values are slightly higher compared to adjacent waters at similar latitudes, aligning with high nutrient concentrations typically observed in these river-influenced waters (Lavoie et al., 2021). These regional patterns could not be captured completely by the global products (Fig. 7c and d). The ability of the ReCAD-NAACOM-pCO₂ product to resolve such regional features demonstrates its potential value for investigating coastal carbon dynamics and their responses to local and regional forcing factors in the NAACOM.

3.5 Evaluating the capacity of the product to detect decadal linear trends of pCO₂

Using pCO₂ products to accurately reconstruct pCO₂ linear trends in coastal regions presents significant challenges due to the high spatial heterogeneity of coastal pCO₂ dynamics. This heterogeneity often leads to sea surface pCO₂ changes that deviate from atmospheric trends (Laruelle et al., 2018). Even when utilizing similar observational datasets, derived products may not consistently reflect the underlying trends. For instance, Wu et al. (2024b) examined the ability of various products to reflect pCO₂ changes in the GoMx, a region where pCO₂ trends exhibit significant spatial variability. Despite this heterogeneity, seven global open-ocean products (listed in Table 1) indicate trends similar to atmospheric pCO₂ across the entire GoMx without regional differences. In contrast, the GoMx-specific regional product developed by Chen and Hu (2019) demonstrates no significant overall trend. The discrepancy in trend detection stems primarily from the design of the regression model and the selection of the input variables. These factors are critical in capturing the complex spatiotemporal variability of coastal pCO₂ and its long-term evolution.

To assess the capability of the product in resolving decadal pCO₂ trends, we conducted an analysis of the pCO₂ evolution using three distinct regions within the NAACOM (three boxes in Fig. 7) as representative examples (Fig. 8). Decadal trends of deseasonalized time series were calculated following the protocol established by Sutton et al. (2022). The LAS (box 1 in Fig. 7) has been identified as an increasing CO₂ sink characterized by a negative pCO₂ rate increase from 2002 to 2021 (Wu et al., 2024b). Our product results for the extended period of 1993–2021 indicate that pCO₂ increased at a rate of +0.44 ± 0.11 µatm yr⁻¹ (Fig. 8a). This rate is significantly lower than the observed atmospheric pCO₂ increase in this region during 2002–2021, which is approximately +1.8 µatm yr⁻¹. These findings corroborate our previous conclusion that the LAS is an increasing CO₂ sink, demonstrating the capability of our product in revealing long-term pCO₂ trends in this dynamic river plume region, extending the analysis period by nearly a decade compared to previous studies. In contrast, the WFS (box 2 in Fig. 7) exhibits an accelerated pCO₂ increase that is faster than the atmospheric pCO₂ of around +2.0 µatm yr⁻¹ (Fig. 8b), aligning with observations reported by Robbins et al. (2018), who found a transition from a CO₂ sink to a source in this region during the 1990s.

Both ReCAD-NAACOM-pCO₂ and SOCAT consistently report a pCO₂ trend of around +2.3 to +2.5 µatm yr⁻¹ in the northern region (box 3 in Fig. 7) over 1993–2021 (Fig. 8c), which is faster than the atmospheric pCO₂ increase (around +2.0 µatm yr⁻¹), suggesting that these areas are becoming a decreasing CO₂ sink. However, limited observational data in this area necessitate cautious interpretation and warrant further validation in future research. Overall, the spatiotemporal heterogeneity in surface-ocean pCO₂ trends across the NAACOM underscores the importance of long-term monitoring to elucidate the drivers of these trends, particularly in regions influenced by major current systems and in areas with limited observational data.

Figure 8 Decadal linear trends of sea surface pCO₂ in three regions of the NAACOM from 1993 to 2021. The blue and red dots are monthly average pCO₂ values (deseasonalized) calculated from SOCAT observations and the reconstructed ReCAD-NAACOM-pCO₂, respectively. The thick lines are linear fitted regression lines. The three regions are the boxes in Fig. 7: the (a) LAS (northern Gulf of Mexico shelf river plume region), (b) WFS, and (c) northern areas. Linear trends are calculated following the protocol established by Sutton et al. (2022). The numbers in parentheses are the number of months with data and p values.

3.6 Evaluating the uncertainty of the product

The uncertainty of the reconstructed pCO₂ values in each grid cell was estimated by accumulating uncertainties from mapping (u_map), gridding (u_grid), measurement (u_obs), and input variables (u_inputs; see Sect. 2.5 for further details on the calculation). To maintain a conservative estimate, we adopted the larger value of 5 µatm as u_obs for all the data points. The gridded fCO₂ values from SOCAT are reported as the averages of all samples collected within each grid cell. Accordingly, u_grid was quantified as the standard deviation of samples within each grid cell, calculated across six subregions. Following the previous literature (Roobaert et al., 2024a; Sharp et al., 2022), u_map was calculated using the RMSE values of the model validation phase reported in Table 2. The uncertainty from the validation set (20 % of X1) was chosen for its sample size that was larger than the independent test set (X2) and for consistency with the 10-fold cross-validation results while avoiding potential underestimation from the training set. u_inputs was calculated using a Monte Carlo simulation (Appendix B). These four sources of uncertainty were evaluated across different subregions of the NAACOM, as shown in Table 3. u_map contributes the largest portion to the total number of uncertainties across all the sub-subregions, with a maximum value of up to 20.12 µatm in the GoMe. Overall, the ReCAD-NAACOM-pCO₂ product demonstrates uncertainties ranging from 16 to 28 µatm across the six subregions and an average uncertainty of 23.25 µatm for the entire NAACOM.

Table 3 Uncertainty estimates for the ReCAD-NAACOM-pCO₂ product across the different subregions of the NAACOM. u_obs, u_grid, u_map, and u_inputs represent the measurement uncertainty, gridding uncertainty, mapping uncertainty, and uncertainty accumulated from input variables, respectively (see Sect. 2.5 for further details). $u_{p_{\mathrm{CO}\mathrm{2}}}$ is the total combined uncertainty (µatm). The subregions are the GoMx, SAB, MAB, GoMe, SS, and GStL&GB.

Our uncertainty estimation employs a conservative estimation using maximum values at the calculation step. This approach likely overestimates the true uncertainty. Despite this conservative method, our calculated uncertainty for the Atlantic margins is comparable to the 43.4 µatm reported by Sharp et al. (2022) for areas within 100 km of the North American Pacific margins, suggesting good performance of our product. It is important to note that our uncertainty calculation assumed independence among all the sources, which is a simplification. Recent research has highlighted that these uncertainties are often correlated (Ford et al., 2024). Future studies should consider these inter-variable correlations to refine uncertainty estimates. In addition, the uncertainties reported in this section and provided in the NetCDF file represent the propagated errors for individual pCO₂ values in each grid cell. Methods to calculate uncertainties in regional averages of pCO₂ or air–sea CO₂ fluxes over specific spatial and temporal domains are detailed in Roobaert et al. (2024a) and Landschützer et al. (2014).

3.7 Challenges and limitations

Even though ReCAD-NAACOM-pCO₂ resolves regional pCO₂ variability with high accuracy in the NAACOM, this product still has room for improvement in the future. Potential areas of improvement include the 0.25° spatial resolution, which is inadequate for resolving submesoscale variability at the scale of 0.1–10 km (McWilliams, 1985). Furthermore, during the independent validation phase, the accuracy of the model-predicted values decreased in the GoMe (R²=0.49) and GoMx (R²=0.46) (Table 2), which may be due to the complex biological and physical conditions in the estuary plume regions in these two gulfs. In this study, we opted not to include chlorophyll-a (Chl-a) concentrations and wind speeds as input variables for model training and prediction. This decision was primarily due to the limited temporal coverage of satellite-derived Chl-a data, which only extends back to 1997 with the launch of the Sea-viewing Wide Field-of-view Sensor (SeaWiFS) satellite (O'Reilly et al., 1998). The inclusion of Chl-a would have restricted the temporal range of our model, potentially limiting its ability to capture long-term trends and variability in pCO₂. Future versions of our model will aim to address this limitation. One potential approach is to develop a two-phase model: one phase for the period before 1997 without Chl-a data and another for the post-1997 period incorporating Chl-a information. Alternatively, we may explore methods to reconstruct historical Chl-a data or use proxy variables that correlate with biological productivity and are available for the entire study period.

In our previous work, we demonstrated that incorporating wind speeds and sea surface roughness data derived from synthetic aperture radar (SAR) could enhance model performance in predicting pCO₂ at submesoscale resolutions (Wang et al., 2024). In this work, we evaluated the inclusion of wind speed as an input variable in our model. However, at the 0.25° resolution employed here, the addition of wind speed data did not significantly improve the model performance (it only increased the R² by 0.1). Moreover, when using the same Monte Carlo simulation approach applied to other variables, incorporating wind speeds would introduce an additional 6 µatm uncertainty to pCO₂ estimates, doubling the input-related uncertainties. Consequently, we excluded wind speeds from our regression model to reduce input-related uncertainties. Despite this omission, our product demonstrates robust capability in resolving regional variations, seasonal cycles, and decadal trends in pCO₂, making it valuable for future studies.

4 Data availability

The reconstructed fCO₂ and pCO₂ and the uncertainty in ReCAD (v1.1) are available as a NetCDF file at https://doi.org/10.5281/zenodo.14038561 (Wu et al., 2024a) and will be updated regularly.

5 Code availability

The Python and MATLAB code used to process the data and create the figures included in this paper is provided at https://github.com/zelunwu/ReCAD (Wu, 2024).

6 Conclusions

The ReCAD-NAACOM-pCO₂ product developed in this study represents a significant advancement in our ability to detect the spatial variations, seasonal cycles, and decadal changes of surface-ocean pCO₂ dynamics in the NAACOM. By leveraging a two-step approach combining random forest and linear regression and a set of environmental predictors, we have created a high-resolution, long-term dataset (1993–2021 period) that captures the complex spatial and temporal variability of pCO₂ across the region. On average, compared with all available samples from the SOCAT observations in our study region, the product has an R² of 0.92, an RMSE of 12.70 µatm, an MAE of 7.55 µatm, and an MBE of 0.13 µatm for the entire NAACOM, with an average uncertainty of 23.25 µatm. The key findings from this study are the following:

1. The product demonstrates high accuracy and reliability, as evidenced by strong performance metrics during the training, validation, and independent test phases across the six subregions.

2. Distinct seasonal cycles are observed between the southern and northern subregions, with the product capturing nuanced features such as elevated pCO₂ levels during fall and winter in the northern areas.

3. Comparison with global products highlights the superior ability of the ReCAD-NAACOM-pCO₂ product to resolve small-scale coastal features and variability.

4. The pCO₂ product successfully reconstructed decadal linear trends that were consistent with previous studies while also revealing a rapid increase in pCO₂ in the northern region of the NAACOM.

While areas of future improvement exist, such as increasing spatial resolution and enhancing accuracy in estuary-plume-influenced regions, the ReCAD-NAACOM-pCO₂ product provides a robust foundation for studying coastal carbon dynamics. This dataset will be valuable for investigating air–sea CO₂ fluxes, assessing ocean acidification impacts, and understanding the role of coastal systems in the NAACOM.

Future research should validate the reconstructed trends, particularly in areas with limited observational data, and explore the mechanisms driving the spatiotemporal variability in pCO₂ across the NAACOM region. Additionally, the methodologies developed here can contribute to a more comprehensive understanding of coastal-ocean carbon dynamics in the face of climate change and have the potential to be applied globally.

Appendix A: Before and after LR calibration

Figure A1 Comparisons of the monthly pCO₂ climatology with SOCAT observations across the six subregions: evaluations before and after LR calibration. Values indicate the mean difference (± 1 standard deviation; blue: before LR and red: after LR) between model-estimated pCO₂ and SOCAT observations over the 12-month period, computed only at grid points and times where SOCAT measurements are available.

Figure A2 Comparisons of deseasonalized monthly pCO₂ anomalies with SOCAT observations across the six subregions: evaluations before and after LR calibration. Values indicate the R² and RMSE between model-estimated pCO₂ and SOCAT observations (blue: before LR; red: after LR), computed only at grid points and times where SOCAT measurements are available.

Appendix B: Monte Carlo simulation in calculating u_inputs

A crucial step in calculating u_inputs is determining the uncertainties of the input variables. In our reconstructed model, there were four variables that needed to be evaluated: SST, SSS, SSH, and pCO_2air. Our general principle was to adopt conservative estimates, using the largest reported uncertainty for each product when available.

SST errors are provided within the OISST product at the grid level. On the global average, OISST reports a mean bias and RMSE of −0.04 and 0.24 °C when compared with the observations on the global average (Huang et al., 2021). For our study region, we calculated the mean SST error across all the grid cells, yielding a value of 0.23 °C.

The SODA database assimilates observational data but does not directly provide SSS error estimates. Given this limitation in uncertainty reporting, we derived an estimate based on the RMSE between the model SSS and observations near our study region, as reported by Carton et al. (2018). Their analysis (their Fig. 8) indicates an RMSE exceeding 0.3 psu in the vicinity of our area of interest. In addition, interpolating the 0.5° SSS data to 0.25° resolution could potentially introduce more errors. To maintain a conservative approach in our uncertainty quantification, we doubled the uncertainty and adopted a value of 0.6 psu as the SSS uncertainty for our calculations.

SSH errors are directly provided in the dataset, which has a mean uncertainty of 1.8 cm in our study region.

pCO_2air, calculated from xCO_2air (MBL references), has a global mean uncertainty of 0.22 ppm.

To propagate these input uncertainties to the final pCO₂ estimate, a Monte Carlo simulation approach was implemented:

1. For each input variable x_i, random perturbations ε_i were generated following a normal distribution N(0,u_i), where u_i represents the uncertainty of the respective variable listed above.

2. Perturbed inputs (x_i+ε_i) were used to calculate pCO₂ with the established model.

3. The difference (Δ_i) between the reconstructed pCO₂ before and after adding the perturbation was computed.

4. Steps 1, 2, and 3 were iterated 100 times for each input variable.

5. The uncertainty contribution from each variable was quantified as the standard deviation of the 100 Δ_i values in each grid cell.

The total uncertainty attributed to the input variables (u_inputs) was then calculated as the square root of the quadratic sum of individual uncertainties:

$$\begin{array}{}\text{(B1)}& u_{\mathrm{inputs}}=\sqrt{u_{\mathrm{SST}}^{\mathrm{2}}+u_{\mathrm{SSS}}^{\mathrm{2}}+u_{\mathrm{SSH}}^{\mathrm{2}}+u_{p{\mathrm{CO}}_{\mathrm{2}\mathrm{air}}}^{\mathrm{2}}}.\end{array}$$

The largest uncertainties propagated from these variables are sourced from SSS and SSH (Fig. B1a and c). Simulating salinity in coastal regions is still challenging due to complex land–ocean interaction. For the SSH, the greatest uncertainties were observed in the GoMe and GStL. Overall, u_inputs is largest in the West Florida Shelf and nearshore waters around the GoMe, with a mean u_inputs uncertainty of 5.9 ± 4.7 µatm for the entire NAACOM.

Figure B1 Uncertainties of pCO₂ accumulated from the different input variables of the model.