Quantifying the Atmospheric CO2 Forcing Effect on Surface Ocean pCO2 in the North Pacific Subtropical Gyre in the Past Two Decades

Despite the well-recognized importance in understanding the long term impact of anthropogenic release of atmospheric CO2 (its partial pressure named as pCO2air) on surface seawater pCO2 (pCO2sw), it has been difficult to quantify the trends or changing rates of pCO2sw driven by increasing atmospheric CO2 forcing (pCO2swatm_forced) due to its combination with the natural variability of pCO2sw (pCO2swnat_forced) and the requirement of long time series data records. Here, using a novel satellite-based pCO2sw model with inputs of ocean color and other ancillary data between 2002 and 2019, we address this challenge for a mooring station at the Hawaii Ocean Time-series Station in the North Pacific subtropical gyre. Specifically, using the developed pCO2sw model, we differentiated and separately quantified the interannual-decadal trends of pCO2swnat_forced and pCO2swatm_forced. Between 2002 and 2019, both pCO2sw and pCO2air show significant increases at rates of 1.7 ± 0.1 μatm yr–1 and 2.2 ± 0.1 μatm yr–1, respectively. Correspondingly, the changing rate in pCO2swnat_forced is mainly driven by large scale forcing such as Pacific Decadal Oscillation, with a negative rate (-0.5 ± 0.2 μatm yr–1) and a positive rate (0.6 ± 0.3 μatm yr–1) before and after 2013. The pCO2swatm_forced shows a smaller increasing rate of 1.4 ± 0.1 μatm yr–1 than that of the modeled pCO2sw, varying in different time intervals in response to the variations in atmospheric pCO2. The findings of decoupled trends in pCO2swatm_forced and pCO2swnat_forced highlight the necessity to differentiate the two toward a better understanding of the long term oceanic absorption of anthropogenic CO2 and the anthropogenic impact on the changing surface ocean carbonic chemistry.

Despite the well-recognized importance in understanding the long term impact of anthropogenic release of atmospheric CO 2 (its partial pressure named as pCO 2 air) on surface seawater pCO 2 (pCO 2 sw), it has been difficult to quantify the trends or changing rates of pCO 2 sw driven by increasing atmospheric CO 2 forcing (pCO 2 sw atm_forced ) due to its combination with the natural variability of pCO 2 sw (pCO 2 sw nat_forced ) and the requirement of long time series data records. Here, using a novel satellite-based pCO 2 sw model with inputs of ocean color and other ancillary data between 2002 and 2019, we address this challenge for a mooring station at the Hawaii Ocean Time-series Station in the North Pacific subtropical gyre. Specifically, using the developed pCO 2 sw model, we differentiated and separately quantified the interannual-decadal trends of pCO 2 sw nat_forced and pCO 2 sw atm_forced . Between 2002 and 2019, both pCO 2 sw and pCO 2 air show significant increases at rates of 1.7 ± 0.1 µatm yr −1 and 2.2 ± 0.1 µatm yr −1 , respectively. Correspondingly, the changing rate in pCO 2 sw nat_forced is mainly driven by large scale forcing such as Pacific Decadal Oscillation, with a negative rate (-0.5 ± 0.2 µatm yr −1 ) and a positive rate (0.6 ± 0.3 µatm yr −1 ) before and after 2013. The pCO 2 sw atm_forced shows a smaller increasing rate of 1.4 ± 0.1 µatm yr −1 than that of the modeled pCO 2 sw, varying in different time intervals in response to the variations in atmospheric pCO 2 . The findings of decoupled trends in pCO 2 sw atm_forced and pCO 2 sw nat_forced highlight the necessity to differentiate the two toward a better understanding of the long term oceanic absorption of anthropogenic CO 2 and the anthropogenic impact on the changing surface ocean carbonic chemistry.

INTRODUCTION
Since industrialization, the global ocean has been a major sink of the increasing atmospheric CO 2 , absorbing ∼25% of anthropogenic CO 2 in recent years (Sabine et al., 2004a;Friedlingstein et al., 2019;Gruber et al., 2019). On one hand, the continuous ocean sink of atmospheric CO 2 (its partial pressure is named as pCO 2 air) is changing ocean carbonic chemistry and the ocean carbon cycle (Borges et al., 2005;Cai et al., 2006;Fujii et al., 2009;Landshützer et al., 2013;Wanninkhof et al., 2013;Xiu and Chai, 2014). On the other hand, the resulting ocean acidification has great potential to degrade marine ecosystems and marine biota, particularly the calcifying organisms such as shellfish and corals (Widdicombe and Spicer, 2008;Doney, 2010;Fabricius et al., 2011;Dickinson et al., 2012;Chan and Connolly, 2013;Davis et al., 2017). Both impacts are closely related to the sustainable development of the marine biota and ecology. Therefore, the anthropogenic effect on surface seawater carbonic chemistry and the potential of the ocean in absorbing anthropogenic CO 2 in the changing world are pressing concerns of the environmental research community.
At present, the study on anthropogenic CO 2 at the sea surface is quite limited. Instead, there are many studies on anthropogenic CO 2 in the ocean interior. The anthropogenic CO 2 stored in the ocean exists in various forms of carbon, originating from the cumulative CO 2 emissions from human activities (e.g., fossil fuel combustion, cement production, etc.) since the beginning of the Industrial Revolution. Several chemical and isotopic tracer approaches have been attempted to estimate the size of this pool of anthropogenic CO 2 (e.g., Sabine et al., 2002Sabine et al., , 2004bLee et al., 2003;Quay et al., 2017;Gruber et al., 2019). However, due to the sparse measurements of chemical tracers in space and time, there is still significant uncertainty in the long-term accumulation rates of anthropogenic CO 2 and the potential of the ocean in continued absorption of anthropogenic CO 2 , making it important to investigate the anthropogenic CO 2 variabilities at the sea surface.
One approach for tracking changes in surface pCO 2 sw is through the collection of autonomous underway and mooring observations over the global ocean (e.g., Surface Ocean CO 2 Atlas (SOCAT, Bakker et al., 2016;Sutton et al., 2019). Many studies focused on the overall variabilities in surface pCO 2 sw and CO 2 flux (e.g., Rödenbeck et al., 2015;Landshützer et al., 2016Landshützer et al., , 2019Gregor et al., 2019;Denvil-Sommer et al., 2019;Iida et al., 2020 among others). However, because of the absence of isotope tracers in the autonomous observing systems of surface pCO 2 sw, it is very challenging to estimate how much anthropogenic CO 2 emissions is driving measured pCO 2 sw. Alternatively, it is known that surface pCO 2 sw is affected by both increasing atmospheric CO 2 forcing (mainly caused by anthropogenic CO 2 emissions) and natural oceanic forcing (e.g., driven by oceanic physical and biogeochemical dynamics) (Fennel et al., 2008;Ikawa et al., 2013;Xue et al., 2016). The effect of atmospheric CO 2 forcing on surface pCO 2 sw (named as pCO 2 sw atm_forced hereafter) actually refers to the changes of surface pCO 2 sw driven by the increase of atmospheric CO 2 . Since the increase of atmospheric CO 2 is due to anthropogenic emissions, the changing rates of pCO 2 sw atm_forced in the past decades can be used to infer the interannual-decadal variations of the anthropogenic signals in surface pCO 2 sw. Yet the pCO 2 sw atm_forced should be differentiated from the total observed pCO 2 sw because of the combination of the natural variability in surface pCO 2 sw (pCO 2 sw nat_forced ). Here pCO 2 sw nat_forced refers to the remaining pCO 2 sw component without atmospheric CO 2 forcing effect, which could be influenced by different physical and biogeochemical processes in the ocean, including the biological activities (i.e., photosynthesis and respiration), ocean warming driven by climate change and anthropogenic CO 2 emissions, and ocean mixing, etc. The effect of all these different processes was regarded as the overall natural oceanic forcing effect on surface pCO 2 sw, It should be clarified that, although we regard all these different oceanic processes to be natural, their changes can still not be completely due to "natural" forcing because these changes in 2002-2019 inherently and implicitly contain atmospheric forcing.
Long time data records are needed to quantify the interannualdecadal trends of pCO 2 sw atm_forced and pCO 2 sw nat_forced in the ocean. Indeed, Sutton et al. (2019) analyzed the time scale of trend detection using 40 autonomous mooring time series of total observed surface pCO 2 sw over the globe, and found that anthropogenic trend detection requires a minimum 8 and 16 years of data records for the sites studies in open ocean and coastal regions, respectively. However, the current global time series observation network of surface pCO 2 sw just starts to approach these time scales, which has made it difficult to track the atmospheric forcing effect for most oceanic environments where the moorings are deployed. Several recent studies attempted to examine the anthropogenic trend in pCO 2 sw based on underway measurements in the past decades (Takahashi et al., 2009(Takahashi et al., , 2014McKinley et al., 2011). For example, Takahashi et al. (2009Takahashi et al. ( , 2014 found that pCO 2 sw is increasing at varying rates of 1.2 ± 0.5∼2.1 ± 0.5 µatm yr −1 in different ocean basins. However, the ship-based measurements are quite limited in both spatial and temporal coverage, leading to many uncertainties in the derived rates. More importantly, these rates are not exactly referring to the atmospheric forcing rates of surface pCO 2 sw, because of the combination of natural variability (i.e., pCO 2 sw nat_forced ) as mentioned above and the difficulty to differentiate and quantify both pCO 2 sw atm_forced and pCO 2 sw nat_forced using in situ observations of surface pCO 2 sw alone.
When combined with in situ surface pCO 2 sw observations, satellite remote sensing has become an important tool for synoptic estimation of surface pCO 2 sw (e.g., Lohrenz et al., 2010Lohrenz et al., , 2018Hales et al., 2012;Signorini et al., 2013;Bai et al., 2015;. Without a spectroscopic method for direct measurements of surface pCO 2 sw from space, it is possible to develop satellite-based pCO 2 sw models through correlations with other related environmental variables. A satellite-based surface pCO 2 sw model also makes it possible to differentiate pCO 2 sw atm_forced from pCO 2 sw nat_forced . Indeed, satellite data accumulated in the past 20 years show great potential to quantify the interannual-decadal trends of the atmospheric forcing effect on pCO 2 sw. However, the past remote sensing studies mainly focused on the retrieval of seasonal surface pCO 2 sw (e.g., Lefèvre et al., 2005;Chierici et al., 2009;Zhu et al., 2009;Borges et al., 2010;Jo et al., 2012;Tao et al., 2012;Marrec et al., 2015;Parard et al., 2015;Le et al., 2019), and are quite limited in predicting interannual variability because of their insufficient parameterization of increasing atmospheric CO 2 forcing (Shadwick et al., 2010;. Therefore, the satellite-based pCO 2 sw algorithms need to be refined to FIGURE 1 | The geolocation of the study site WHOTS (annotated in black star), and the general climate mode of the North Pacific in terms of PDO based on the HadISST data set (Rayner et al., 2003) for the period 1870-2019. enable their capabilities in assessing the interannual trends of pCO 2 sw atm_forced and pCO 2 sw nat_forced .
The Woods Hole Oceanographic Institution Hawaii Ocean Time-series Station (WHOTS) near Hawaii in the North Pacific Subtropical Gyre (NPSG) maintains high resolution surface pCO 2 sw observations. It provides an important open ocean reference for Hawaiian coral reefs (Dore et al., 2003;Sutton et al., 2017;Terlouw et al., 2019), thus is important to know the interannual-decadal trends of the atmospheric forcing effect on surface pCO 2 sw for a better understanding of the long-term ocean acidification and oceanic absorption of anthropogenic CO 2 . The WHOTS station was selected mainly because it has sufficient field data records for anthropogenic trend detection as mentioned above. WHOTS is located at station ALOHA (A Long-term Oligotrophic Habitat Assessment) (Karl and Church, 2018) in the NPSG (Figure 1), under the large-scale climate forcing of Pacific Decadal Oscillation (PDO). Several published studies investigated the interannual variability of the upper ocean carbon cycle at this station (Dore et al., 2003(Dore et al., , 2009Keeling et al., 2004;Palevsky and Quay, 2017). For example, based on a 14-year time series (1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002) at ALOHA, Brix et al. (2004) found that surface pCO 2 sw and isotopic 13 C/ 12 C showed long-term increase and decrease (yet no rates were provided), respectively, and they attributed it to the uptake of isotopically light anthropogenic CO 2 from the atmosphere. Using the same data time series of pCO 2 sw at ALOHA, Dore et al. (2003) found that the significant decrease in CO 2 sink in 1989-2001 was driven by the climate variability in salinity (Lukas and Santiago-Mandujano, 2008). Later based on a longer data record of 19 years (1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007) at ALOHA, Dore et al. (2009) presented a pCO 2 sw increasing rate of 1.88 µatm yr −1 . In contrast, with a synthesis of 35 years of observations in the North Pacific, Takahashi et al. (2006) found the interannualdecadal change in surface pCO 2 sw is mostly correlated with the increases of sea surface temperature (SST) and anthropogenic CO 2 . Therefore, it is necessary to further investigate the effects of both anthropogenic CO 2 emissions and the climate-driven natural variability in the ocean on surface pCO 2 sw. However, to date, no studies have differentiated these two forcing effects.
Considering the importance of addressing this knowledge gap to promote our understanding of the ocean capability in absorbing anthropogenic CO 2 in the long run, here we for the first time differentiate the atmospheric forcing and natural forcing effects on surface pCO 2 sw, that's, pCO 2 sw atm_forced and pCO 2 sw nat_forced , based on a novel satellite-based pCO 2 sw model developed in this study. Specifically, the objectives of this study include: (1) develop a satellite-based surface pCO 2 sw model at WHOTS, which should be able to capture the interannual-decadal variabilities in pCO 2 sw and differentiate pCO 2 sw atm_forced and pCO 2 sw nat_forced , and (2) quantify the interannual-decadal trends of both terms in the past 2 decades, and understand its implications for ocean acidification and long term oceanic uptake of anthropogenic CO 2 . Although the study was conducted at the WHOTS station, the findings in this study may provide insight on the interannual-decadal trends of pCO 2 sw driven by atmospheric and natural forcing effects, respectively, in other global subtropical open ocean regions. More importantly, the approach developed in this study can be extended to other regions with sufficient data available.

Data
The WHOTS station (22.7 • N, 158 • W) is located in the subtropical oligotrophic region of the North Pacific and is operated by the Woods Hole Oceanographic Institution (WHOI). Field data time series [including surface pCO 2 sw, and pCO 2 air, SST, and sea surface salinity (SSS)] at this station collected between 2004 and 2018 at led by NOAA's Pacific Marine Environmental Laboratory and were obtained from the National Centers for Environmental information (NCEI) 1 (Sutton et al., 2012). Specifically, the pCO 2 data were measured by a nondispersive infrared gas analyzer (LI-COR TM , model LI-820), which has a sampling frequency of every 3 h, with an accuracy of 2 µatm (or better) (Sutton et al., 2014;Sabine et al., 2020). Surface pCO 2 sw data were collected at a water depth of <0.5 m, and the pCO 2 air data were collected at 1.2 m above the sea surface. SST and SSS were obtained from a CTD (SBE16) integrated in the autonomous CO 2 mooring system. The details of data collection, processing, and quality control can be found in Sutton et al. (2014). These data were binned to daily time series to remove the diurnal variations (i.e., 0.4∼3.4 µatm), which are not considered in this study. The data time series were then averaged at monthly scales as presented in in Figures 2A,B,D. The Hawaii Ocean Time-series (HOT) program also maintains ship-based monthly sampling of surface pCO 2 sw calculated from dissolved inorganic carbon (DIC) and total alkalinity (TA) at this location ( Figure 2D). We chose to use the high-frequency data from the WHOTS buoy mainly to assure that there are sufficient data available to develop the machine learning pCO 2 sw model and the monthly averages of the modeled pCO 2 sw should have lower bias than the monthly observed pCO 2 sw at HOT. NASA standard daily SST (Figure 2A) and 8 day Chlorophylla (Chla, mg m −3 ) ( Figure 2C) Level-3 data products (version R2018.0) covering the study region for the period of July 2002-December 2019 with a spatial resolution of ∼4 km were downloaded from the NASA Goddard Space Flight Center (GSFC) 2 . These Level-3 data products were derived from measurements by the Moderate Resolution Imaging Spectroradiometer (MODIS) on the Aqua satellite.
Clearly there are lots of data gaps in the field measurements (e.g., SST, pCO 2 air, pCO 2 sw, Figures 2A,C). Full record of SST is obtained from MODIS. For a full data record of pCO 2 air at WHOTS between 2002 and 2019, daily time series of atmospheric xCO 2 (in unit of ppm) at Mauna Loa Observatory (MLO) in Hawaii between 2002 and 2019 were obtained from the NOAA ESRL Global Monitoring Laboratory (2019), and this data was regarded as the atmospheric xCO 2 at WHOTS over the study period considering the close distance between Mauna Loa and WHOTS. To calculate the corresponding pCO 2 air at WHOTS from the atmospheric xCO 2 following the standard operating procedures (Weiss, 1974;Dickson et al., 2007), ancillary daily data of sea surface air pressure (in unit of atm) and specific humidity (in unit of%) were obtained from the National Centers for Environmental Prediction (NCEP), with a spatial resolution of 2.5 • . The derived pCO 2 air ( Figure 2D) together with the MODIS data (Figures 2A,C) were used to estimate pCO 2 sw between 2002 and 2019 based on the developed pCO 2 sw model. It should be clarified that, for broader impact, one main reason in choosing MODIS SST and NCEP ancillary data instead of other in situ data at the WHOTS mooring was to demonstrate our model capability in dealing with the uncertainties in each parameter, particularly when extending our method to other locations or regions where field measurements could be limited.

Methods
Surface pCO 2 sw is mainly controlled by four oceanic processesthe thermodynamic effect, biological activity, physical mixing, and air-sea CO 2 exchange (Fennel et al., 2008;Ikawa et al., 2013;Xue et al., 2016). Accordingly, satellite-derived variables of SST, SSS, and Chla are commonly used to estimate surface pCO 2 sw from remote sensing in past studies (Olsen et al., 2004;Ono et al., 2004;Lohrenz and Cai, 2006;Sarma et al., 2006;Lohrenz et al., 2010Lohrenz et al., , 2018Nakaoka et al., 2013;Chen et al., 2016Chen et al., , 2017. However, these algorithms are quite limited in capturing the long-term trend in pCO 2 sw, mainly because of the insufficient parameterization of the anthropogenic or atmospheric CO 2 forcing effect on pCO 2 sw. Feely et al. (2006), and Landshützer et al. (2013Landshützer et al. ( , 2016 have investigated the interannual and decadal variations of pCO 2 sw and CO 2 flux under the anthropogenic CO 2 forcing, yet to better quantify this effect, further studies are needed to differentiate the warming effect of SST from the atmospheric effect on surface pCO 2 sw and quantify both effects separately. Dore et al. (2003) found that the significant increase of pCO 2 sw at ALOHA in 1989-2001 was mainly caused by the increase of SSS due to excess evaporation over this period, 2 https://oceancolor.gsfc.nasa.gov/ suggesting that the physical changes in the subtropical North Pacific may affect the ocean biogeochemistry including surface pCO 2 sw. Yet in this study, SSS was found to have little effect on pCO 2 sw (R = 0.102 at p > 0.05, which explains 1% nges in pCO 2 sw) at the WHOTS station over the period of 2004∼2018, as also found by Sutton et al. (2017) which shows a small effect (<5%) of salinity changes on pCO 2 sw increase. The SMOS satellite maintains the longest SSS data record since 2009 (Font et al., 2009(Font et al., , 2013, however, a comparison between the field SSS and SMOS-derived SSS shows a very large uncertainty of 1.1 for SSS ranging between 34.5 and 35.5 at WHOTS. As such, SSS was not used in the model. The mixed layer depth (MLD) could drive the interannual dynamics of surface pH at ALOHA (Dore et al., 2009), yet considering the lack of MLD data from remote sensing and the covariations of SST and MLD dynamics, we chose to use SST alone to indicate the effect of warming and mixing on surface pCO 2 sw. Therefore, the inputs of the satellite pCO 2 sw algorithm included observed SST and pCO 2 air, and concurrent MODIS-derived Chla, as well as Julian day (Jday) normalized sinusoidally to "tune" the seasonal cycles of pCO 2 sw (Friedrich and Oschlies, 2009;Signorini et al., 2013;Chen et al., 2016Chen et al., , 2017, and the output was modeled pCO 2 sw (Eq. 1). In total, there were 3074 matched data samples between 2004 and 2017. Within this dataset, data samples collected in 2016 (N = 311) were kept for independent validation considering its near full coverage in each month (other years do not); the remaining were randomly divided into two groups: one for model training (N = 1,934), and the other for model validation (N = 829).
Various approaches have been used to model pCO 2 sw from remote sensing, such as polynomial regression, mechanistic semianalytical approach, machine-learning approaches (Friedrich and Oschlies, 2009;Jo et al., 2012;Landshützer et al., 2013;Bai et al., 2015;Moussa et al., 2016;Lohrenz et al., 2018).  did extensive comparisons of these approaches and found that, the Random Forest based Regression Ensemble (RFRE) was the most robust one in modeling pCO 2 sw. Therefore, this approach was used in this study with model parameters locally tuned for the WHOTS station (Eq. 1). RFRE is one type of machine learning technique, which ensembles many weighted regression trees to implement the random forest algorithm (Breiman, 1996(Breiman, , 2001James et al., 2013) in Matlab (R2017a). For better model generalization, the RFRE takes advantage of each regression tree via bootstrap aggregation (or bagging) (Breiman, 1996;James et al., 2013) in model parameterization. In the model training phase, the ensemble regression trees grow independently on a drawn bootstrap replica of the training dataset. That's, each regression tree can randomly select a subset of predictors at each split and can involve many splits in the algorithm. This manipulation greatly reduces the correlations among the developed regression trees, resulting in improved independency among the regression trees. The mean square error was used as loss function to adjust the model performance in each iteration. Briefly, there are two important parameters to define the RFRE model structure: the minimum leaf size and the number of regression trees. Leaf size refers to the number of data samples used in each node of a regression tree, and its minimum thus determines the splits and depth of a regression tree. By trial and error, these two parameters were optimized to be 8 and 28, respectively. With these settings, the RFRE model became stable and had the best model statistics, thus it was used to predict pCO 2 sw. See  for more details of the RFRE approach.
We varied SST, Chla, and pCO 2 air by ± 1 • C, and ± 20%, and ± 5 µatm, respectively, to examine the sensitivity of the model to changes in each variable. The changes are based on the uncertainties in the MODIS-derived SST and Chla (Gregg and Casey, 2004;Mélin et al., 2007;Hu et al., 2009) as well as on the seasonal variations in pCO 2 air.
The modeled pCO 2 sw is the sum of pCO 2 sw atm_forced and pCO 2 sw nat_forced . Just as its name implies, the pCO 2 sw nat_forced refers to the pCO 2 sw without atmospheric CO 2 forcing, thus based on the model developed following Eq. 1, the pCO 2 sw nat_forced was calculated by assuming that the pCO 2 air remained at the same level as in in the start year (i.e., 2002) of the study period (Eq. 2). The pCO 2 sw atm_forced was defined as the difference between the modeled pCO 2 sw (Eq. 1) and pCO 2 sw nat_forced (Eq. 3). To quantify the natural forcing effect, the net atmospheric CO 2 forcing effect over the study period (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017) remained at exactly zero by keeping the pCO 2 air values in the model at the same level as in 2002. By doing so, both the derived pCO 2 sw nat_forced and pCO 2 sw atm_forced are relative quantities to the year of 2002, which should be higher than those derived by referring to pre-industrialization. However, either referring to 2002 or other years only affects the absolute values of these quantities, and they would affect the changing rates of trends in both pCO 2 sw nat_forced and pCO 2 sw atm_forced in the past two decades that we are interested in.
pCO 2 sw nat_forced = f RFRE [SST, log 10 (Chla), pCO 2 air @2002 , where the pCO 2 air @2002 means the pCO 2 air data in 2002-2019 remained at the same level as in 2002 by assuming that there is no additional atmospheric effect referred to 2002.
Trends in pCO 2 sw, pCO 2 sw atm_forced , pCO 2 sw nat_forced , pCO 2 air, SST, and Chla were quantified based on their monthly anomalies, which were derived by subtracting the monthly climatologies from the monthly averages between 2002 and 2019 using least-square technique. Figure 3 shows the performance of the RFRE-based pCO 2 sw algorithm in both model training and validation. Clearly, most of the data pairs of the observed and modeled pCO 2 sw followed closely along the 1:1 line, with a RMSD of 2.2 µatm (0.6%) and R 2 of 0.98. The additional independent validation (Figure 4) using the data time series in 2016 also shows good consistency between the observed pCO 2 sw and modeled pCO 2 sw, with a RMSD of 4.3 µatm (1.1%) and R 2 of 0.87.

RESULTS
The RFRE model is more sensitive to changes in SST and pCO 2 air than to changes in Chla (Figure 5). Statistically, with + 1 • C (-1 • C) added to SST, the modeled pCO 2 sw was higher (lower) than the original pCO 2 sw, with RMSD of 9.7 µatm (2.6%) [8.0 µatm (2.1%)], R 2 of 0.89 (0.93), and MB of 8.5 µatm (-6.8 µatm). The resulting pCO 2 sw shows slight underestimation and overestimation in cases of 20% increase and 20% decrease in Chla, with MB of 1.3 and −1.3 µatm, respectively. With + 5 µatm in pCO 2 air, the new pCO 2 sw was estimated higher than the original pCO 2 sw, with RMSD of 5.7 µatm (1.6%), R 2 of 0.90, and MB of 3.7 µatm. With −5 µatm in pCO 2 air, the new pCO 2 sw was  In the North Pacific subtropical gyre at the WHOTS station, time series of pCO 2 sw between 2002 and 2019 was obtained using this RRFE-based pCO 2 sw algorithm, with good consistency to the observed pCO 2 sw in the overlapped time periods (Figure 6). Overall, the pCO 2 sw follows the same seasonal pattern as SST from high values in summer to low values in winter, with a seasonal magnitude of ∼50 µatm, in the opposite phase of pCO 2 air (Figure 6). In addition, the pCO 2 sw was lower than the pCO 2 air most of the time over the years, suggesting a continuous CO 2 flux from the atmosphere to the ocean.
Both pCO 2 sw and pCO 2 air show significant increase between 2002 and 2019 (Figure 6). After removing the seasonality signals, statistically, the pCO 2 sw had a mean rate of 1.7 ± 0.1 µatm yr −1 (R 2 = 0.80, at p < 0.05), lower than the rate of pCO 2 air (2.2 ± 0.1 µatm yr −1 , R 2 = 0.99, at p < 0.05), as shown in Figure 7. The pCO 2 sw nat_forced shows a significant increasing rate of 0.2 ± 0.1 µatm yr −1 (R 2 = 0.07, at p < 0.05) on average in the study period. In contrast, the pCO 2 sw atm_forced , which is just driven by the atmospheric CO 2 forcing, had a mean rate of 1.4 ± 0.1 µatm yr −1 (R 2 = 0.84, at p < 0.05), but tended to plateau since 2016. Indeed, the pCO 2 sw without the thermodynamic effect (i.e., pCO 2 nonT, Chen and Hu, 2019) had similar interannual patterns as pCO 2ant at a mean rate of 1.2 ± 0.1 µatm yr −1 . Correspondingly, the Chla time series did not show any trends over the years while the SST was increasing at an overall rate of 0.03 ± 0.01 • C yr −1 (R 2 = 0.07, at p < 0.05). This warming trend could be influencing the pCO 2natural trend.
Clearly, there are some visible trends (e.g., <10 years) particularly in SST and pCO 2 sw nat_forced different from those over the 20-year time frame (Figure 7). To further investigate the trends in each variable, we quantified the rates of each for a variety of periods starting between 2002 and 2015, ending between 2006 and 2019, with durations ranging from 5 to 18 years (Figure 8). It is found that, at confidence level of >95%, the SST had a negative and positive rate of −0.1 ± 0.02 • C yr −1 and 0.1 ± 0.05 • C yr −1 for periods ending in ≤2013 and >2013, respectively ( Figure 8A). Again, the Chla did not show any trend over the years. Correspondingly, pCO 2 sw nat_forced shows a very similar pattern as the rates in SST, with a negative rate of −0.5 ± 0.2 µatm yr −1 for periods ending in ≤2013, and a positive rate of 0.6 ± 0.3 µatm yr −1 for periods ending in >2013. The anthropogenic forcing on atmospheric pCO 2 tends to accelerate over the study period consistent with the published studies (Canadell et al., 2007), with a rate of 1.7 ± 0.1 µatm yr −1 for periods ending in ≤2011, and a rate of 2.3 ± 0.2 µatm yr −1 for periods ending in beyond 2011, and the acceleration is getting even stronger (2.4 ± 0.1 µatm yr −1 ) after 2016. As a result, the pCO 2 sw shows a lower rate ( Correspondingly, the pCO 2 sw atm_forced shows similar but significantly weaken signals (at p < 0.05) in these three time frames, with rates of 1.6 ± 0.3 µatm yr −1 , 1.8 ± 0.5µatm yr −1 , and 0.9 ± 0.5 µatm yr −1 , respectively.
The sensitivity of the pCO 2 sw model to each input variable indicates not only the model's capacity in tolerating the uncertainty of each variable, but also the model's response to real changes in each variable. Specifically, the positive feedback of modeled pCO 2 sw to changes in SST are consistent with the thermodynamic effect on pCO 2 sw (increased SST leads to an increase in pCO 2 sw and vice versa). The negative response of the pCO 2 sw model to Chla suggests that the increase (decrease) in Chla indicates stronger (weaker) biological uptake of oceanic CO 2 , therefore, the resulting modeled pCO 2 sw was lower (higher) than without the Chla perturbation. Although the Chla level at the WHOTS station is consistently low (Figure 2C), the sensitivity analysis here suggests the necessity of including Chla in the model to better modulate the seasonal variations of surface pCO 2 sw. Yet it should be noted that, Chla is only a proxy to indicate the overall biological activities that could affect surface pCO 2 sw. Although there is no visible change in surface Chla, still there could be possible changes in the phytoplankton community and net community production. The insignificant responses of the pCO 2 sw model to the 20% change in Chla suggest the model is insensitive to uncertainties in the satellite Chla. For the same reason, the biological uptake of CO 2 tends to have a quite limited effect on pCO 2 sw in the oligotrophic ocean, consistent with previous studies . For regions where satellite Chla is not available due to severe cloud coverage (e.g., some tropical and high latitude zones), a first examination of the Chla effect on surface pCO 2 sw using field observations (if there are) is suggested to determine the potential bias that would be resulted in pCO 2 sw if Chla is not included in the model. The changes of pCO 2 air directly affect the gradient between pCO 2 air and pCO 2 sw, which drives the air-sea CO 2 exchange, thus, it is reasonable to see a positive response of the pCO 2 sw model to changes in pCO 2 air. The resulting increase (MB = 3.7 µatm) in pCO 2 sw was slightly weaker than the assigned increase of 5 µatm in pCO 2 air, which may be due to the ocean's increasing Revelle Factor and reduced buffering capacity of seawater (Fassbender et al., 2017).

Interannual Changes of pCO 2 sw Driven by Natural and Atmospheric Forcing
In response to the accelerating rates of pCO 2 air, the modeled surface pCO 2 sw shows different rates at various time intervals. Specifically, the 5 year pCO 2 sw trends we derived for the periods of 2007-2011, 2008-2012, and 2009-2013 are high at rates of 2.5, 2.1, and 2.5 µatm yr −1 , respectively, which are higher than the relatively low rates in period of 2003-2007 visually interpreted from Figure 1 in Dore et al. (2009). To further examine the trends in pCO 2 sw, we analyzed the ship-based monthly pCO 2 sw datasets at ALOHA from HOT program (used in Dore et al., 2009). Indeed, the 5 year HOT-based pCO 2 sw trends starting in 2007-2008 did show low values, but these low values are insignificant at p > 0.05, yet no such statistics was available in Dore et al. (2009). For the 5 year pCO 2 sw trend starting in 2009, the HOT-based pCO 2 sw and our modeled pCO 2 sw show close trends of 2.5 and 2.2 µatm yr −1 , respectively, at p < 0.05. Meanwhile, the overall trend we detected in surface pCO 2 sw (i.e., 1.7 ± 0.1µatm yr −1 ) in period of 2002-2019 was a bit smaller than that (i.e., 1.88 µatm yr −1 ) in period of 1988-2007 found in Dore et al. (2009) and that (i.e., 2.4 µatm yr −1 ) in period of 2003-2014 presented in Sutton et al. (2017). This could be reasonable considering the different physical and biogeochemical dynamics on decadal time scales and the acceleration of ocean acidification in the western North Pacific (Ono et al., 2019). Besides, it should be noted that the ship-based monthly pCO 2 sw dataset is derived from measurements of DIC and TA collected approximately once a month to compose this monthly dataset. In contrast, our monthly pCO 2 sw is based on the daily modeled pCO 2 sw and is validated thoroughly with daily-averaged in situ measurements at WHOTS. Therefore, the trends in the modeled pCO 2 sw we derived here should be reliable with high confidence. Also, the mooring measures pCO 2 sw at surface of <0.5 m, while the ship-based HOT data were based on the mean measurements within 0-30 m, which could be another potential source for the discrepancy. In the North Pacific subtropical gyre (represented by the WHOTS station), the interannual changes of surface pCO 2 sw is mainly driven by both SST and pCO 2 air (Figures 7, 8 and Table 1), consistent with the published studies . Despite the little impact of SSS on pCO 2 sw shown in our study period (2002-2019), a further experiment with SSS added into our model was conducted. It shows that the inclusion of SSS did not result in any significant difference in the modeled pCO 2 sw and pCO 2 sw nat_forced . Considering the important impact of SSS on pCO 2 sw in 1989-2007 presented in Dore et al. (2003), it seems that the effect of SSS depends on the specific study periods. Here we prefer to exclude SSS from our model mainly considering the large error (i.e., 1.1) in the SMOS SSS at present. With more accurate SSS data available from satellites in the future, it could be possible to include SSS to better model the variations of pCO 2 sw, particularly the effect of rainfall minus precipitation on pCO 2 sw in any time periods. However, most of the published studies directly regarded the interannual trend of pCO 2 sw as the trend of anthropogenic pCO 2 sw. It should be noted that the anthropogenic pCO 2 sw refers to the pCO 2 sw impacted by atmospheric CO 2 increases, thus most of the reported anthropogenic trend of pCO 2 sw actually refers to  the total rate of pCO 2 sw (Takahashi et al., 2009(Takahashi et al., , 2014McKinley et al., 2011;Sutton et al., 2019), which also includes the natural variability of pCO 2 sw driven by the general oceanic processes (e.g., thermodynamics, ocean mixing, biological activities).
In this study, both the natural and atmospheric CO 2 forcing effects on pCO 2 sw were separately quantified. The rates in pCO 2 sw nat_forced over the study period follow a similar pattern as those in SST with a correlation coefficient (R) of 0.82, indicating that the interannual trend signals in pCO 2 sw nat_forced are mainly driven by SST, at least over the study period of 2002-2019. The cooling characteristics in SST between 2002 and 2012 resulted in a significant negative rate in pCO 2 sw nat_forced , and the warming effect since 2013, which were also reported in previous studies (Sutton et al., 2017;Terlouw et al., 2019), leads to a significant positive rate in pCO 2 sw nat_forced . In addition to the global warming effect on SST, the interannual SST dynamics could also be attributed to the changes in MLD because of the ocean mixing effect on SST. As such, the interannual variations in pCO 2 sw nat_forced could also be driven by the MLD changes, and more DIC enriched waters would be entrained into the surface when MLD deepens and SST decreases (Dore et al., 2009). Overall, it seems that the rate of pCO 2 sw nat_forced tends to correspond to decadal oscillations in SST between cooling and warming periods associated with PDO (Yasunaka et al., 2014;Newman et al., 2016;Landshützer et al., 2019). Indeed, the interannual PDO ( Figure 7E) shows very similar variation patterns to the SST ( Figure 7A) with a significant correlation of R = 0.53 (Table 1). Specifically, the PDO decreased progressively from 2004 to 2012, was low in 2011-2012, reached a maximum in 2015, and then decreased from 2015 to 2019. As a result, the pCO 2 sw nat_forced also shows a significant correlation (R = 0.41, see Table 1) with PDO, suggesting the large scale climate forcing also contribute to the natural oceanic forcing effect on surface pCO 2 sw.
With the exclusion of pCO 2 sw nat_forced , the pCO 2 sw atm_forced rates were significantly smaller than the corresponding pCO 2 sw rates in various time intervals (Figure 8). Although pCO 2 sw atm_forced is mainly driven by the oceanic uptake of increasing atmospheric CO 2 (R = 0.91), it shows distinctively different patterns in changing rates from that of the pCO 2 air over various time intervals in 2002-2019. This different response of pCO 2 sw atm_forced toward pCO 2 air seems mainly caused by the buffering effect of dissolved CO 2 in seawater (Egleston et al., 2010). However, for the tendency of pCO 2 sw atm_forced to plateau after 2016, there could be several potential explanations depending on the condition of air-sea CO 2 fluxes. Specifically, it would be reasonable to observe a plateau signal in pCO 2 sw atm_forced if there is little change in air-sea CO 2 fluxes after 2016; yet if the dissolved CO 2 keeps increasing after 2016, the little response in pCO 2 sw atm_forced would tend to suggest that a larger fraction of dissolved CO 2 stays in forms of other carbonate species (i.e., HCO 3− , CO 32− ), significantly lowering the Revelle factor and enhancing the ocean's buffering capacity in recent years; and if there is a decrease in air-sea CO 2 fluxes after 2016, it would be likely that a fration of bicarbonate and carbonate species are converted to dissolved CO 2 , which would lower the ocean's buffering capacity and promote ocean acidification. Xue and Cai (2020) found that TA minus DIC can be used as a proxy for deciphering ocean acidification. Here using the ship-based monthly TA and DIC data in the study period, we found a significant decreasing trend in TA minus DIC over the years (Figure 9A), which suggests a strong ocean acidification in the study period. However, the changing rates of TA minus DIC is distinctively higher in recent years since 2014 (Figure 9B), suggesting a stronger ocean acidification and weaker buffering capacity in the past few years. Indeed, ocean acidification has shifted the carbonate chemistry speciation and lowered the CaCO 3 saturation state (Orr et al., 2005;Doney et al., 2009;Krug et al., 2011), yet further studies are needed to investigate and quantify the changing patterns of the air-sea CO 2 flux and the carbonate species over the past decades. In general, the oceanic uptake of anthropogenic CO 2 is resulting in more rapid changes in carbonic chemistry in the surface ocean and accelerating ocean acidification Ono et al., 2019), yet a revisit of such phenomenon is needed when more satellite/field data are available in the coming years.

Implications
Long time series data are required to investigate the anthropogenic effect on surface pCO 2 sw. However, the field data are always limited in both spatial and temporal coverage. For example, few of the 40 global pCO 2 sw mooring stations have data coverage of >10 years (Sutton et al., 2019), and the global field pCO 2 sw database (i.e., SOCAT or LDEO, Bakker et al., 2016;Takahashi et al., 2019), although greatly accumulated in recent years, still has data gaps in some regions and at some time intervals. More importantly, it is impossible or difficult to separate the pCO 2 sw atm_forced and pCO 2 sw nat_forced signals apart based on purely field measurements to better quantify the anthropogenic forcing impact on surface pCO 2 sw. Instead, with the related environmental variables observed from satellites, surface pCO 2 sw models using satellite data and other ancillary data can be developed and applied to the full satellite data record over the past ∼20 years. Besides, SSS measurements from SMOS and SMAP satellite have been available since 2009 and 2015, respectively, with longer and accurate data records available, the interannual and decadal trends in surface pCO 2 sw as well as the natural forcing and atmospheric CO 2 forcing components can be further studied. The recovered long time series of pCO 2 sw can be used to quantify both pCO 2 sw atm_forced and pCO 2 sw nat_forced accordingly. The findings of decoupled changing rates in pCO 2 sw atm_forced and pCO 2 sw nat_forced in this study highlight the necessity of differentiating the two, in order to have a better understanding of the long term oceanic absorption of anthropogenic CO 2 and its buffering capacity in the long term. Therefore, this study sets a template for future study to examine both natural and anthropogenic or atmospheric CO 2 forcing effects on pCO 2 sw in various oceanic systems over the past   In (B), the X-axis and Y-axis represent the start year and end year, respectively, of each rate analyzed, the diagonal lines (i.e., 5, 10, and 15 years) indicate the length of trend periods, and a white cross is superposed on the plot when the p-value was >0.05.
decades, toward an improved understanding of anthropogenic forcing on surface pCO 2 sw. Specifically, the pCO 2 sw in the North Pacific subtropical gyre shows various increase rates in response to the increasing pCO 2 air between 2002 and 2019. The accelerating increase rates in pCO 2 air and the weaker rates in pCO 2 sw indicate stronger gradients between pCO 2 air and pCO 2 sw, which implies an accelerated oceanic CO 2 uptake and ocean acidification. If the warming effect continues following the decadal pattern in SST in recent years since 2010, a steady rate of ∼0.8 ± 0.1 µatm yr −1 in pCO 2 sw nat_forced (see Figure 8E) would be expected in the coming few years. The weaker rate in pCO 2 sw atm_forced in recent years in response to the accelerating rate in pCO 2 air implies a lower ocean buffering capacity leading to more rapidly changing oceanic carbon chemistry and ocean acidification, yet further study in this field is needed to promote our knowledge and understanding.
Based on observations at WHOTS, the present work demonstrated the necessity in differentiating the atmospheric forcing and natural forcing effects on surface pCO 2 sw, and show unprecedented information on their interannual-decadal trends over both short and long time scales. The WHOTS station is located in the North Pacific Subtropical Gyre, therefore, the results and findings should be referential to understand the overall surface pCO 2 sw dynamics for a broader impact of the ocean in absorbing anthropogenic CO 2 , particularly under both anthropogenic CO 2 forcing and natural oceanic forcing (Henson et al., 2016).
More importantly, the pCO 2 sw model was developed using satellite-derived environmental data and other ancillary data, thus the model is capable to tolerate the uncertainties involved in each variable as demonstrated in the sensitivity analysis. This is of great importance and significance to locations or areas where very limited data are available. Specifically, with these limited field observations of surface pCO 2 sw, it would be possible to develop a surface pCO 2 sw model with related environmental variables from satellite and ancillary data from NCEP to differentiate the two forcing effects following Eqs 1 and 2. With nearly 20 years of satellite data records, it would be straightforward to extend the current study to other oceanic regions to investigate the interannual-decadal surface pCO 2 sw dynamics by differentiating the atmospheric forcing and natural forcing effects toward a better understanding of the ocean in absorbing anthropogenic CO 2 and its impact on the surface ocean carbonate chemistry.

CONCLUSION
The rate of anthropogenic or atmospheric CO 2 forcing pCO 2 sw in surface seawater has been difficult to characterize because of the interaction of natural variability in pCO 2 sw and the requirement of long time series data records. In this study, we show that a remote sensing algorithm applied to the WHOTS station in the North Pacific subtropical gyre can reveal the interannual-decadal variability of surface pCO 2 sw between 2002 and 2019. Such an ability enables the separation of atmospheric CO 2 forced pCO 2 sw (pCO 2 sw atm_forced ) from natural variability in pCO 2 sw (pCO 2 sw nat_forced ). We believe that this is the first time such atmospheric CO 2 forced pCO 2 sw and natural oceanic processes driven pCO 2 sw are mathematically differentiated and their interannual-decadal changing rates are statistically quantified. Results show unprecedented information on their interannual-decadal rates over both short and long time scales at the WHOTS site. With the availability of ocean color data and other ancillary data globally, it is straightforward to extend the current study to other oceanic regions.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. This data can be found here: The WHOTS mooring dataset analyzed in this study is available at National Centers for Environmental information (NCEI) under https://www.nodc.noaa.gov/ocads/ oceans/Moorings/. The MODIS ocean color data are available at the NASA Goddard Space Flight Center (GSFC) https: //oceancolor.gsfc.nasa.gov/. The atmospheric xCO2 data at Mauna Loa is available at the NOAA ESRL Global Monitoring Laboratory under https://www.esrl.noaa.gov/gmd/dv/data/ index.php.

AUTHOR CONTRIBUTIONS
SC designed the study, processed, analyzed the data, and wrote the manuscript. AS contributed the mooring data and main concept definition. CH and FC contributed data analysis and manuscript writing. All authors commented on the manuscript.

FUNDING
This work was supported by the National Natural Science Foundation of China (NSFC) projects (42030708, 41906159, and 41730536). The pCO 2 observations at WHOTS were supported by the Office of Oceanic and Atmospheric Research of the NOAA, U.S. Department of Commerce, including resources from the Global Ocean Monitoring and Observation program.