The Bering Sea is located in the North Pacific Ocean, bordered by Russia to the west and Alaska to the east. The Bering Sea system consists of a deep offshore basin (with a maximum depth of 3,500 m) and a wide continental shelf (of<200 m) (Askren, 1972; Coachman, 1986; Stabeno et al., 1999). The central Bering Sea basin is characterized as an iron-limited, high-nutrient, low-chlorophyll region (Leblanc et al., 2005; Sugie et al., 2013). However, the eastern shelf of the Bering Sea is broad and shallow, being >500 km wide, with an average depth of only 70 m, resulting in a long growing season and high annual primary production (Rho and Whitledge, 2007). Tidal mixing along the shelf break also results in a highly productive region, often referred to as the green belt, where nitrate from the deep basin and iron from the shelf converge to mix into the transmissive zone (Springer et al., 1996). This high primary productivity across the shelf and slope in turn supports a wide variety of pelagic and benthic predators, as well as rich fishery resources (Fissel et al., 2016). For the coast of Alaska, streams flow into the Bering Sea, and eddy-induced downwelling transports relatively warm seawater from the surface layer (Miura et al., 2002; Mizobata et al., 2002). Substantial portions of the shelf are normally covered by sea ice during winter months, and sea ice retreat typically occurs in April–May (Stabeno and Bell, 2019).

The circulation in the Bering Sea basin is commonly described as a cyclonic gyre, with the Kamchatka Current (KC) flowing southward forming the western boundary current and the Bering Slope Current (BSC) flowing northward forming the eastern boundary current ( Figure 1 ). Circulation in the Bering Sea is strongly affected by the Alaskan Stream (AS), which flows through many channels in the Aleutian Islands into the Bering Sea. Inflow into the Bering Sea is balanced by outflow through the Kamchatka Strait, so circulation in the Bering Sea basin may be more properly described as a continuation of the North Pacific subarctic circulation (Stabeno and Reed, 1994; Stabeno et al., 1998). Circulation on the eastern Bering Sea shelf is broadly northwestward. The net northward transport through the Bering Strait, although important for the Arctic Ocean, has little effect on the circulation in the Bering Sea basin. However, it does play a dominant role in determining the circulation of the northern shelf (Stabeno et al., 1999).

A series of pCO₂ data was collected from the Surface Ocean CO₂ Atlas (SOCAT; available at http://www.socat.info/) (Pfeil et al., 2013; Sabine et al., 2013). The latest SOCAT version, SOCAT v2021, contains 30.6 million quality-controlled fCO₂ observations from 1957 to 2020 (Bakker et al., 2016). At the time of this analysis, SOCATv2022 was the most up to date version available. In this study, we selected the best fCO₂ by its quality flag (with an accuracy of better than 5 μatm); the total amount of the original observations was 411,497 during 2003 to 2019, and the range was from north of the Aleutian Islands to south of the Bering Strait (approximately from 52 N to 66 N).

To obtain collocated and contemporaneous field-measured pCO₂ and satellite data products, we gridded the observations into 4 km spatial resolution gridded sea surface pCO₂ data and did not directly use the gridded dataset listed in SOCATv2022, which has a spatial resolution of 1 degree × 1 degree. In addition to the consideration of sampling biases caused by uneven coverage, we averaged the underway data with a 1/24° (~4-km) monthly grid and obtained 24,695 gridded fCO₂ data. The gridded fCO₂ data density is highly uneven, with most of the surface fCO₂ data collected on the Bering continental shelf during the summer ( Figures 2A–D ). In addition, we converted fCO₂ to pCO₂ by using an empirical formula (Weiss, 1974a; Körtzinger, 1999). The measured years of SOCAT data in the Bering Sea from 2003 to 2019 in the shelf and basin are shown in Figures 2E, F . During the study period, there are enough data in 2009-2019 and data are unavailable in 2005-2008.

Following the 4-km grid of the observed data, we matched a variety of environmental parameters to construct matched datasets for training and testing of the retrieved model. Although the spatial and temporal resolutions of these parameters are different, we used the nearest neighbor interpolation method to adjust the spatial resolution to a standard monthly 4-km resolution grid in this study.

Sea surface temperatures were downloaded from the National Oceanic and Atmospheric Administration (NOAA). We used the 1/4° Daily Optimum Interpolation Sea Surface Temperature (OI-SST), which is a long-term climate data record that incorporates observations from different platforms (satellites, ships, buoys, and Argo floats) into a regular global grid, in this research (available at https://www.ncei.noaa.gov/products/optimum-interpolation-sst) (Reynolds et al., 2007; Banzon et al., 2016; Huang et al., 2021). The associated chlorophyll a (Chla) concentrations were taken from a global 4-km-resolution monthly data product compiled by NASA’s MODIS Ocean Science Team (Nasa Goddard Space Flight Cente, 2014), which was calculated using an empirical relationship derived from in situ measurements of Chla and blue-to-green band ratios of in situ remote-sensing reflectance (Carder et al., 2004). The sea surface height (SSH) and mixed-layer depth (MLD) were extracted monthly from the Global Ocean Physics Reanalysis Glorys12V1 (PHYS 001-030) products at a resolution of 1/12° (~8 km) from the Copernicus Marine Environment Monitoring Service (CMEMS) (https://doi.org/10.48670/moi-00021). We extracted the dry air mole fractions of CO₂ (xCO₂) from CarbonTracker (version 2019B; https://doi.org/10.25925/20201008), which is a modeling system developed by NOAA that provides global estimates of surface atmosphere fluxes and four-dimensional CO₂ data (Jacobson et al., 2020). The bathymetry was extracted from the global ETOPO1 database (from NOAA, 2022) (Amante and Eakins, 2009). The list of all products used in the calculations is summarized in Table 1 .

In the process of calculating the flux, we added extra parameters to estimate the air–sea CO₂ flux, including sea surface 10-m wind and atmospheric pressure. Sea surface 10-m wind data were from the cross-calibrated multiplatform (CCMP) vector wind analysis dataset Mears et al., 2022, which combines satellite measurements, in situ measurements, and a background wind field into complete maps of ocean winds every 6 h with biases lower (by almost 0.5 m/s) than those of buoys (Mears et al., 2019). We also compared many satellite and reanalysis wind speed products and finally chose the CCMP wind speed product due to its suitable temporal coverage and accurate characterization. Atmospheric pressure data were from the European Centre for Medium-Range Weather Forecasts Reanalysis v5 reanalysis (ERA5) dataset, which embodies a detailed record of the global atmosphere, land surface, and ocean waves from 1950 onward, providing an hourly output throughout with 37 pressure levels (Hersbach et al., 2020). We used the sea surface pressure from a 0.25° square latitude–longitude grid. The sea ice concentration (SIC) data were obtained from the National Snow and Ice Data Center (NSIDC) Climate Data Record at a spatial resolution of 25 km and a monthly temporal resolution (Meier et al., 2013).

The global air–sea CO₂ flux is often estimated by a bulk method in which in situ pCO₂ measurements in seawater are combined with a wind-speed-dependent gas transfer velocity (Weiss, 1974b; Wanninkhof, 2014):

where f l u x (mol m⁻² yr⁻¹) is the air–sea CO₂ flux; k w is the gas transfer velocity; s o l is the solubility of CO₂ in seawater (mol m⁻³ µatm⁻¹); p C O 2 w a t e r (µatm) and p C O 2 a i r (µatm) are the partial pressure of surface ocean CO₂ and atmospheric CO₂, respectively, in the marine boundary layer; SIC is the sea ice concentration of ocean area covered by sea ice; U 10 is the wind speed 10 m above sea level; and Sc denotes the Schmidt number calculated from SST and sea surface salinity (SSS). We used SIC in this study because sea ice, an imperfect barrier to gas exchange, has an important effect on air-sea CO₂ flux (Long et al., 2011; Loose et al., 2014; Butterworth and Miller, 2016). However, whether the effect of sea ice on gas exchange is nonlinear (Loose et al., 2009; Loose et al., 2014) or linear (Butterworth and Miller, 2016; Prytherch et al., 2017) is still under debate, and for simplicity, only linear ice corrections were used in this work. In addition, we used 99% when the SIC was greater than 99% to allow for air-sea CO₂ exchange through fracture, lead, and brine channels (Semiletov et al., 2004). All data used here can be sourced in Section 2.2. The sink–source pattern of CO₂ flux was determined by the difference between ocean pCO₂ and atmospheric pCO₂ (ΔpCO₂); that is, a positive value of ΔpCO₂ corresponds to a CO₂ source, while a negative value corresponds to a CO₂ sink.

The carbon sink in a total area is usually conducted by the area integral of air-sea CO₂ flux. In the Bering Sea, we summed the results of multiplying each grid flux by the area of each grid of three parts (the basin, the shelf and the area of the Aleutian Islands coast). The net CO₂ flux in a subregion was estimated by multiplying the mean CO₂ flux density among the available pixels and normalizing to the total area. In addition, the high-latitude pCO₂ is missing due to the unavailable Chla. In the processing of ocean color satellite data, atmospheric correction usually fails under high solar zenith angles with weak sunlight signals, which usually occurs in winter in high-latitude areas. Therefore, there are large blank areas of missing Chla and corresponding reconstructed pCO₂ and flux in Jan. and Dec. To reduce the effect of this vacancy on the area integral of air-sea CO₂ flux estimation, we used the seasonal average, which is calculated through adjacent months at the same location, to replace the vacancy.

The root-mean-square error (RMSE) and correlation coefficient (R²) were used as the standard statistical metrics to measure model performance in this study. The RMSE was calculated for the dataset as follows:

where n is the index of the samples, O b s are the observation measurements, and E s t are estimates from the model.

The seasonal difference calculation of pCO₂ was calculated by using the following four steps. First, a third-order polynomial (to account for trends) and fourth-order harmonic function (to reproduce the seasonality) were fitted to all grid data to reproduce the full seasonal cycle (Graven et al., 2013).

where t is time in years and T is the period, chosen here as 1 yr. Second, we recreate the seasonal cycle of a certain year by fitting equation (4) to every full analysis year, including the year before and after that, to create a three-year running time series. 2003 and 2019 are reconstructed using the two following or preceding years. Third, from the resulting harmonic function f(t) segments, the mean of the months January, February and March was used as the winter average and July, August and September as the summer average. Then, we calculated winter minus summer for a certain year as the seasonal difference. Finally, the trends in these seasonal differences are calculated from the slope of the linear regression line fit to the 17-year time series.

Based on the large match-up dataset, we tried various machine learning methods, including support vector machine, eXtreme gradient boosting, Gaussian-process-based kernel machine learning approach (GPR), and random forest. Among them, the GPR method was found to have optimal results in the retrieval of pCO₂ in the Bering Sea (the comparison of algorithm performance can be found in the Table 2 ).

The GPR approach allows handling the model selection issue within a Bayesian framework in a completely automatic way, thus offering the potential advantage of avoiding the traditional empirical and tricky tuning of the free parameters of the model (Pasolli et al., 2010). GPRs are usually highly flexible and accurate for prediction over new inputs closer to training data points, and recent studies have demonstrated the effectiveness of GPR for time-series gap-filling applications (Mateo-Sanchis et al., 2018; Pipia et al., 2019; Belda et al., 2020). Furthermore, GPR has exhibited robust performance in some sea surface biophysical retrieval approaches and therefore has outperformed other machine learning regression algorithms, such as random forest and artificial neural networks, especially in retrieving biophysical parameters from satellite data (Verrelst et al., 2012; Verrelst et al., 2015; Rivera-Caicedo et al., 2017). Pipia et al. (2021) proposed a GPR-based gap-filling strategy and generation of multi-orbit cloud-free green leaf area index maps with an unprecedented level of detail. Svendsen et al. (2020) used a Gaussian process to estimate the chlorophyll content, colored dissolved matter, and inorganic suspended matter from multispectral data acquired by the Sentinel-3 OLCI sensor. Blix et al. (2018) used the GPR machine learning model to estimate the complex waters against in situ measurements over Lake Balaton in 2017 based on the Ocean and Land Color Instrument (OLCI) onboard the Sentinel 3A satellite. The effectiveness of this approach in the context of estimation of biophysical parameters offers a new alternative to state-of-the-art regression methods such as those based on artificial neural networks and support vector machines (Pasolli et al., 2010). Based on the characteristics of the Bering Sea and the previous research, we inferred that the strengths of GPR on biophysical environmental parameters (Verrelst et al., 2012; Verrelst et al., 2015; Rivera-Caicedo et al., 2017) and the advantages in learning from small samples (Ly et al., 2021) can help the algorithm to reconstruct sea surface pCO₂ well in the Bering Sea. Thus, we chose to use the GPR method to retrieve pCO₂ in the Bering Sea.

We used time and location data from in situ pCO₂ measurements to identify a set of pCO₂ values with matched environmental parameters. These gridded monthly standard products were used in the pCO₂ model development. The predictor variables of the model were a set of environmental parameters associated with the Bering Sea carbonate system based on some previous studies (Song et al., 2016; Chau et al., 2022), including SST, Chla, MLD, SSH, xCO₂, and bathymetry.

In many other studies, these parameters represent physical, chemical, and biological proxies determining the distribution of sea surface pCO₂. Here, both SST and MLD were used in the model to account for the influence of thermodynamic and DIC variation in pCO₂ (Landschützer et al., 2013). Meanwhile, Chla is commonly used as a “tuning” parameter to emphasize biological activity and primary productivity (Chen and Hu, 2017). SSH, however, was used to characterize the horizontal current flows and to capture the mixing and convergence of water masses at the ocean surface (Chau et al., 2022). In addition, bathymetry was collected to identify the two different carbonate systems: the ocean basin and the shelf.

We tried to incorporate sea surface salinity into the input data. Due to the large uncertainty of various salinity datasets in the Bering Sea, the input data produced low-quality results, so we removed the salinity parameter. The SSH usually contains contributions from high- and low-frequency changes in steric (temperature and salinity) signals and ocean circulation (Meyssignac and Cazenave, 2012; Zhou et al., 2012; Volkov et al., 2022). Here, we used SSH to relatively optimize the performance of our algorithm in the lack of high-quality salinity data. Remote sensing reflectance has been used as input in other studies; however, in the Bering Sea, the model performance with remote-sensing reflectance as inputs was not satisfactory and may form a redundancy of input quantities similar to chlorophyll, as chlorophyll products are usually retrieved from remotely sensed reflectance data. Thus, after trying various combinations, we finally chose six independent environmental predictors: SST, SSH, MLD, Chla, xCO₂, and bathymetry in the pCO₂ model development; details of the datasets are given in Section 2.2.

We divided the underway pCO₂ match-up dataset into three parts. First, to ensure the final independent validation, we selected data from the two longest cruises: One was from a cruise in the Bering Basin in October 2009, and the other was from a cruise in the continental shelf sea in July and August 2019 ( Figure 3A ). These are the longest cruises in each of the two regions to ensure the enough data number and spatial coverage, and we used them as a third completely independent data to validate the algorithm performance, beyond the training dataset and testing data. These cruise-based observations were used for an extra-validation, which put a higher demand on the robustness of the algorithm. To verify that the two cruises are not special, we randomly chose random sets of cruises covering the basin and shelf for independent validation and developed a model with the same parameters and methods. The accuracy of the four models is very close to that of the original model, demonstrating that our modeling is stable with changes in the independent validation (four cases of the model validation results from different random selections of independent cruises can be found in the Supplemental Material ).

Then, we derived a relationship between the pCO₂ field and environmental predictors based on a GPR method, with 75% of the remaining gridded pCO₂ (produced from SOCAT v2022 cruise data, see more details in Section 2.2) used as the training dataset and 25% as the model testing dataset ( Figures 3B, C ). We therefore could establish nonlinear relationships between six independent environmental predictors (SST, SSH, MLD, Chla, xCO₂, and bathymetry) and pCO₂. The reconstructed pCO₂ field matches the SOCAT data well. In both the basin and continental shelf, a high agreement between the training pCO₂ and the gridded SOCAT v2022 data is demonstrated, with an overall mean R² of 0.93 and RMSE ~ 19.18 µatm (based on 16,985 observations) ( Figures 3D, E ). This good fit can also be found in the testing, with an RMSE = 12.56 μatm in the basin and 23.08 μatm in the shelf ( Figures 3F, G ).

In independent validation, we obtained a mean R² of 0.94 and RMSE ~ 23.49 µatm (based on 1424 observations) ( Figure 4 ). Although there are some errors in regions of extreme gradient variations, we were still able to capture the correct trends in these rapidly changing regions ( Figures 4A, B ). The validation indicated that the GPR method can accurately reconstruct pCO₂ in the Bering Sea ( Figure 4C ). However, because of the lack of in situ measurements, our estimates along the western coast of the Bering Sea (e.g., the Gulf of Anadyr) are tentative extrapolations without sufficient validation. We note that there is no in situ data in the self of the Bering Sea during 2005-2008 and no data available in the basin before 2008, and these data gap may induce the uncertainty of time series analysis. We try to develop satellite retrieval algorithm and to reconstruct the satellite-based pCO₂ from 2003 to 2019, but we will also emphasis such uncertainty and discuss the time series from 2009 to 2019.

Based on a novel GPR method with multisource remote-sensing data, we generated monthly 1/24° (~4 km) resolution pCO₂ and reconstructed the air–sea CO₂ flux during the period from January 2003 to December 2019. The 17-year climatological monthly pCO₂ data exhibit marked seasonal and spatial variation (ranging from 250–450 µatm) in the Bering Sea ( Figure 5 ). In the Bering Sea basin, sea surface pCO₂ decreases continuously from April to July and reaches its annual minimum (326.2 μatm) in August and rises from September to November, eventually reaching a peak (441.2 μatm) in December; a similar pattern can be found in the shelf area where sea surface pCO₂ reaches its annual minimum (291.6 μatm) in May with a peak (416.3 μatm) in January. Although both the continental shelf and the ocean basin exhibit typical subpolar pCO₂ seasonal patterns, with low pCO₂ in summer and high pCO₂ in winter and spring, the pCO₂ in the shelf region declines 2–3 months earlier than that in the ocean basin area and remains at lower values for longer periods. Our reconstructed dataset reconstructs the large summer gradient in pCO₂ between the basin and shelf well ( Figure 5 ), and this phenomenon is also confirmed by the measured data ( Figure 2 ). The coastal region of Alaska has experienced extremely high pCO₂ in the summer. The seasonal amplitude of seawater pCO₂ is almost a factor of 3–4 greater than that of the atmospheric pCO₂ ( Figure 6A ).

Based on the satellite-derived pCO₂, sea surface wind speed, xCO₂, and SST, we estimated the long-time-series air–sea CO₂ flux in the Bering Sea from 2003 to 2019. The seasonal variation in air–sea CO₂ fluxes in the Bering Sea is closely related to the seawater pCO₂ pattern, with a mean air-sea CO₂ flux of approximately ± 0.7 mol m⁻² mon⁻¹ ( Figure 7 ). The Bering Sea basin exhibits dramatic seasonal amplitudes, with sea water releasing CO₂ to the atmosphere from November to April and becoming a sink during June to October ( Figure 6B ). The continental shelf, however, continues to behave as a sink for the atmosphere and reaches its maximum absorption in spring ( Figure 6B ), except along the Alaskan coast, which exhibits a weak source in summer, especially from June to August. In addition, seawater in the southern Bering Sea (coastal Aleutian Islands) is a total source to the atmosphere, with the strongest emission occurring in spring. As mentioned in Section 2.3, remote-sensing reflectance data and Chla products under high solar zenith angles are usually invalid in winter in high-latitude areas. Therefore, there are large blank areas of missing Chla and corresponding reconstructed pCO₂ and flux in Jan. and Dec. in Figures 5 , 6 , which would lead to uncertainty in the estimation.

We averaged the valid pCO₂ data for each month and then calculated the annual average based on the monthly average results. From 2003 to 2019, atmospheric pCO₂ in the Bering Sea increased steadily (2.13 ± 0.09 μatm yr⁻¹, P< 0.01). Meanwhile, sea surface pCO₂ in the basin area increased more rapidly than that in the atmosphere (2.82 ± 0.63 μatm yr⁻¹ in 2003-2019 and 2.95 ± 0.57 μatm yr⁻¹ in 2009-2019, P< 0.01). However, there is no clear trend on the continental shelf, where sea surface pCO₂ fluctuates in the range of 347 ± 12 μatm, consistently maintaining a low value of<360 μatm with seasonal fluctuations of ±63 μatm ( Figure 8A ). As a result, the trend of the annual mean air–sea pCO₂ difference (ΔpCO₂) in the basin area rises at a slow rate, while it declines rapidly in the shelf area. Small variations occurred in ΔpCO₂ during 2003-2008, which might contain some uncertainty as the in-situ data are unavailable during this period. During 2009-2019, the decreasing rate of ΔpCO₂ (4.37 μatm yr^-1) in the shelf is approximately three times that of the increasing rate in the basin; against the background of continuously increasing atmospheric CO₂, sea surface pCO₂ exhibited insignificant changes in continental shelf, which also exacerbated the air-sea ΔpCO₂ (Figure 8C). Thus, the Bering Sea basin was a slowly developing CO₂ source, while the carbon sink in the Bering Sea shelf area has been increasing significantly. Meanwhile, the annual mean wind speed in the Bering Sea changed little during the study period, and the mean wind speed in the sea basin was always 1 m s⁻¹ greater than that on the shelf ( Figure 8B ). Note that in 2008-2019, the mean wind speed exhibited a slight increasing trend: 0.052 m s^-1 yr^-1 (shelf, P = 0.013) and 0.043 m s^-1 yr^-1 (basin, P = 0.038). Then, we analyzed the annual variation of sea surface pCO₂ and atmospheric pCO₂ in each month to check if the trends are various in season ( Figures 8D, E ). During 2003-2019, the atmospheric pCO₂ increased at a similar rate (2.19 ± 0.02 μatm yr^-1) in all months ( Figure 8B ), while sea surface pCO₂ exhibited different trends in different months. In general, except April-June, the growth rates of sea surface pCO₂ were around 1.63 ± 0.98 μatm yr^-1, p-value< 0.05). However, the growth rates were not significant in April-June (p-value > 0.05), indicating that sea surface pCO₂ fluctuated within a certain range without significant annual trends in these months. This may be due to the strong photosynthesis and high primary productivity during May to August and the biological effects facilitate the carbonate buffering system to some extent. When weakened productivity and deep DIC upwelling occur simultaneously, large growth rates occurred in October (2.6 μatm yr^-1) and November (2.5 μatm yr^-1), which even higher than the atmospheric growth rate ( Figure 8F ).

We estimated the annual absorption of the sea basin, continental shelf, and whole Bering Sea ( Figure 9 ) with the area integral of the air-sea CO₂ flux. Note here that the total CO₂ uptake does not exactly match the trend of the mean air–sea CO₂ flux because CO₂ uptake is the integral of the gridded flux with respect to the gridded area. Consequently, southern grids will produce greater CO₂ uptake or release than those in the north under the same carbon flux. As a result, the whole Bering Sea exhibited an increasing carbon sink with an integral increase in air-sea CO₂ fluxes from 6 to 19 TgC over 17 years. The variation of total air-sea carbon flux was insignificant in 2003-2008, and consistently increased from 2011 to 2015, maintaining at 15-19 TgC since 2015. However, due to the insufficient observations for the continental shelf in 2003-2009, we cannot confirm whether the GPR algorithm underestimates air-sea CO₂ sink, which may have contributed to the rapid increase in 2011-2015. Although the Bering Sea is generally characterized as an important carbon sink with an average air–sea CO₂ flux of 8.87 TgC yr⁻¹, this absorption is almost entirely driven by the on-shelf regions and the green belt. The deep basin, in contrast, is a high-nutrient, low-chlorophyll region with increasing sea surface pCO₂ and acts as a CO₂ source with increasing CO₂ emissions during the study time period. The upwelling along the coastal Aleutian Islands contributed to the regional high pCO₂ (Kelley and Hood, 1973; Sapozhnikov et al., 2011), whereas this upwelling-induced CO₂ emission did not exhibit a significant trend over the study period from 2003 to 2019. The increased wind speed from 2008-2019 also contributed to a larger air–sea CO₂ flux. In conclusion, the carbonate systems in the basin and continental shelf of the Bering Sea evolved in different directions under the background of the continuous rise in atmospheric pCO₂.

To explore the drivers of air–sea CO₂ fluxes in the Bering Sea, we used a method similar to that developed by Takahashi et al. (2002) and calculated the effects of temperature (T) and all other non-temperature processes (nonT) on the sea surface pCO₂ differences. We used temperature-normalized pCO₂ to remove the temperature control; this can be expressed as follows:

where the quantity n p CO 2 is the temperature-normalized pCO₂; S S T is the sea surface temperature; and S S T mean is the climatological annual mean calculated from the optimum interpolation sea surface temperature during 2003–2019, as (∂ln pCO₂/∂T) is +4.23% °C⁻¹ (Takahashi et al., 1993). Therefore, the effects of temperature (T) and all other non-temperature processes on the sea surface pCO₂ differences can be estimated as follows:

where T ( p CO 2 ) represents the impact of temperature (temperature effect) and nonT ( p CO 2 ) weights the impact of all other processes, such as the biological effect, upwelling and other mixing processes. p CO 2 ( atm ) is atmospheric pCO₂. The sum of temperature and non-temperature pCO₂ components is equal to the seasonal variability at any given time under the assumption that the annual mean pCO₂ should be well represented (Takahashi et al., 2002; Cross et al., 2014). Because we aim at the mechanisms affecting seasonal air-sea CO₂ fluxes, the atmospheric pCO₂ was selected as a reference when calculating nonT ( p CO 2 ) (Wang et al., 2022). Hence, the sum of T ( p CO 2 )  and  nonT ( p CO 2 ) is equal to ΔpCO₂(= p CO 2 − p CO 2 ( atm ) ) at any given time in this study.

Based on Eqs. (5) and (6), we estimated the effects of temperature and non-temperature effects in the Bering Sea from pCO₂ data estimated from satellite data. Figure 10 shows that the spatial distributions of T ( p CO 2 )  and  nonT ( p CO 2 ) differ. Bates et al. (2011) indicated that the biological effect, that is, the balance of photosynthesis and respiration, is the most dominant process driving variability in non-temperature processes on the continental shelf in the Bering Sea. Our findings also suggested that non-temperature effects, including biological contributions, upwelling and seasonal mixing, together dominate the surface pCO_2. As T ( p CO 2 )  and  nonT ( p CO 2 ) . have opposite effects, and the non-temperature effect is a factor of ~2–4 greater than the temperature effect. The mean temperature effect of pCO₂ is 12.7 μatm, whereas the mean non-temperature effect is −51.8 μatm. This difference is consistent with the previous conclusion that temperature is not the primary control mechanism in the Bering Sea (Takahashi et al., 2002; Song et al., 2016). In addition, the effect of temperature on pCO₂ is positive, and the strongest temperature effects occur within the Norton Sound (>40 μatm) ( Figure 10A ). In contrast, the effect of the non-temperature component is negative and more widespread, with thmost impacted area being in the Gulf of Anadyr (<-100 μatm) ( Figure 10B ).

To analyze the pCO₂ distributions in more detail, we classified water masses in the Bering Sea with respect to the subregions (see Figure 1A ). Based on the previous method of water mass classification in the Bering Sea, we discuss the characteristics of six water masses in the Bering Sea: basin winter water (BWW), basin summer water (BSW), shelf winter water (SWW), shelf summer water (SSW), melt water (MW), and Alaskan coastal water (ACW) (Schumacher et al., 1979; Coachman, 1986; Clement et al., 2005; Jinping et al., 2006; Danielson et al., 2017; Lin et al., 2019; Yamashita et al., 2019; Abe et al., 2021; Hirawake et al., 2021; Wang et al., 2022). The classified boundaries of temperature and salinity are given in Table 3 . The Bering Sea basin and shelf water bodies can be distinguished from summer and winter by a boundary of 4°C. Salinity can be used to distinguish among the sea basin, shelf, and coast. Low-salinity water near the coast has two sources: melt ice and riverine input, which can be distinguished by temperature. Due to the mixing of MW and ACW, we classified the mixing water masses as either MW or ACW according to their temperature. Due to the lack of cruise data in the nearshore area, we cannot distinguish the mixing of MW with river water. During thermodynamic dominance, pCO₂ increases with higher temperature. However, in subpolar regions, pCO₂ increases in winter when the seasonal mixed-layer depth is deepening, with DIC-rich water from deep upwelling to the sea surface (Sarmiento, 2013; Fay and Mckinley, 2017; Gallego et al., 2018). In addition, the vertical mixing of high-DIC subsurface waters results in sea surface pCO₂ increasing. We also found that the pCO₂ in the basin was significantly higher in winter (BWW) than in summer (BSW), suggesting that the Bering Sea basin is not a region of dominant temperature effect on pCO₂ ( Figure 11 ). Similar seasonal variations occurred in the shelf and basin with lower pCO₂ in summer water than in winter water ( Figures 11C, F ). Meanwhile, this seasonal difference is more significant in the shelf, which may be because phytoplankton production is higher in the shelf area than in the high-nutrient low-chlorophyll basin. The pCO₂ coastal region of Alaska is much more strongly influenced by the low-salinity MW and ACW; however, different types of changes are generated. The MW diluted the seawater, resulting in a lower pCO₂ ( Figure 11C ). In the river plume area, pCO₂ is often at a low level because of the photosynthesis of phytoplankton. However, the estuarine areas and nearshore areas are also affected by a large amount of riverine humic and other organic matter from the catchment basin, and strong respiration, mineralization and degradation processes lead to high pCO₂ in rivers and estuarine areas. Nishimura et al. (2012) observed high concentrations of dissolved Fe and humic-type fluorescence intensity in a nearly peak river water discharge period in the northern Bering Sea shelf (Yukon River estuary region and St. Lawrence Island polynya region). Strong mineralization and degradation of high concentrations of terrestrial organic matter will lead to high pCO₂ in estuaries and offshore areas. Murata (2006) also indicated that coccolithophorid Emiliania huxleyi blooms occurred in coastal Alaska simultaneously in summer with CO₂ changes (increase in CO₂). During the growth of coccoliths when the calcium carbonate shell is formed, CO₂ is released at the same time, resulting in an increase in sea water pCO₂. It is different from diatom blooms with CO₂ drawdown (Murata, 2006).

We also calculated the temperature and non-temperature effects on pCO₂ in each of the six typical water masses ( Figure 12 ). The opposite trend between temperature and non-temperature effects is evident in the Bering Sea basin and the continental shelf, with temperature providing positive values and negative values from non-temperature effects in summer and vice versa in winter. However, regardless of the seasonal change, the non-temperature effect always outweighs the temperature effect, becoming the most dominant control mechanism in the Bering Sea, which can usually be considered a result of seasonal variations in DIC (Takahashi et al., 2009; Gallego et al., 2018). In addition, this non-temperature effect reaches a maximum in the summer over the continental shelf, which may be due to the combined effect of low DIC resulting from the high net community productivity (primary productivity – community respiration) produced and low DIC diluted by high summer river runoff (Mathis et al., 2010; Mathis et al., 2011; Cross et al., 2012).

Many studies have found that the amplitude of pCO₂ seasonal differences have been increasing substantially in recent decades based on observed data or model-based predictions. Landschützer et al. (2018) indicated that the seasonal differences increase at a rate of 1.5 ± 1.1 μatm per decade observed at the Hydrostation ‘S’/Bermuda Atlantic time-series study site and at the Hawaiian Ocean time-series station, which was 3.8 ± 2.4 μatm per decade. Furthermore, model-based predictions also support the increasing seasonal difference in pCO₂ in the South Atlantic, Pacific, and North Atlantic (Doney et al., 2009; Hauck and Völker, 2015; Mcneil and Sasse, 2016).

In the Bering Sea basin, the satellite-derived pCO₂ reveals no remarkable change in the seasonal amplitude (1.3 μatm per decade); however, it has been rapidly increasing on the continental shelf from 2003 to 2019, approaching 14 μatm per decade (the calculation method can be found in the supplemental materials ). Although a considerable amount of winter data is missing in the shelf of the Bering Sea from our reconstructed pCO₂ dataset (which both resulted from the failure of satellite processing at the large zenith angle in high latitudes in winter and the inability to observe water signals due to ice cover), there is some uncertainty in the time-series trend estimation. We regarded that such missing data might be systematic, and the trend estimated from satellite products could be referable, as Landschützer et al. (2018) also showed the average increasing rate at 55°N can be ~8 μatm per decade. In addition, the seasonal cycle in high-latitude areas has a maximum in winter, leading to a positive winter-minus-summer difference in pCO₂. Whether in the sea basin or on the shelf, Chla is relatively stable in autumn and winter, while the peak caused by blooms in spring and summer fluctuates (1.5-2.3 μg/L in the continental shelf and 2.5-3.0 μg/L in the basin) with significant interannual variation ( Figure 13A ). The temperature shows a large seasonal amplitude ( Figure 13B ). The sea surface height of the basin increases throughout the study time period, while there is no significant trend in the shelf ( Figure 13C ). In addition, the mixed-layer depth shows the most complex interannual variability, with an abrupt decrease in the mixed-layer depth occurring in the basin (since 2013), while at the same time, the shelf area is increasing with a smaller trend ( Figure 13D ). Compared to the Bering Sea basin, the greater seasonal amplitude of pCO₂ in the continental Bering Sea is mainly the result of decreasing pCO₂ in summer ( Figure 13E ), combined with the stable growth of atmospheric xCO₂. The final interannual variation in ΔpCO₂ is shown in Figure 13F .

Previous studies have identified two possible mechanisms that may drive the increase in seasonal differences: the long-term increase in the mean concentration of CO₂ in the surface ocean and the reaction of the added CO₂ with the carbonate ions in seawater (Landschützer et al., 2018). Taking these hypotheses into consideration, we analyzed the mechanisms in the continental Bering Sea. First and foremost, because sea surface pCO₂ increases, the pCO₂ variation with temperature also increases; thus, even the same seasonal temperature difference can strengthen the seasonal pCO₂ amplitude. This effect is more applicable to explain the variation in seasonal amplitude in the basin. The continental Bering Sea remained at 348 ± 12 μatm during 2003–2019, exhibiting an unclear trend. Second, the reaction of the continuously added CO₂ with the carbonate ions in continental seawater results in a reduction in the capacity of the surface ocean CO₂ system. Hauck et al. (2015) showed that the seasonal drawdown of carbon by biological production contributes an even larger share to the total CO₂ uptake as the Revelle factor increases in this century in the Southern Ocean and other regions governed by strong seasonality (Hauck et al., 2015). In other words, the same amount of biological production will lead to stronger CO₂ uptake over the course of the century. The strong seasonal decrease in DIC might increase the uptake of CO₂ in the Bering Sea. Furthermore, this is particularly vital in systems with strong seasonality, such as at high latitudes, where darkness inhibits biological production in winter. Increased seasonal variability in sea surface pCO₂ largely enhances the influence of ocean acidification on marine organisms by exposing them earlier to higher levels of ocean acidification (Doney et al., 2009; Bates et al., 2014; Lauvset et al., 2015). This leads to transitions that cross critical thresholds harmful to marine ecosystems and fisheries, such as hypercapnia and low saturation regarding calcium carbonate (Gruber et al., 2012; Lauvset et al., 2015; Mcneil and Sasse, 2016).