Remote Sensing of Sea Surface p CO 2 in the Bering Sea in Summer Based on a Mechanistic Semi-Analytical Algorithm ( MeSAA )

The Bering Sea, one of the largest and most productive marginal seas, is a crucial carbon sink for the marine carbonate system. However, restricted by the tough observation conditions, few underway datasets of sea surface partial pressure of carbon dioxide (pCO2) have been obtained, with most of them in the eastern areas. Satellite remote sensing data can provide valuable information covered by a large area synchronously with high temporal resolution for assessments of pCO2 that subsequently allow quantification of air-sea carbon dioxide 2 flux. However, pCO2 in the Bering Sea is controlled by multiple factors and thus it is hard to develop a remote sensing algorithm with empirical regression methods. In this paper pCO2 in the Bering Sea from July to September was derived based on a mechanistic semi-analytical algorithm (MeSAA). It was assumed that the observed pCO2 can be analytically expressed as the sum of individual components controlled by major factors. First, a reference water mass that was minimally influenced by biology and mixing was identified in the central basin, and then thermodynamic and biological effects were parameterized for the entire area. Finally, we estimated pCO2 with satellite temperature and chlorophyll data. Satellite results agreed well with the underway observations. Our study suggested that throughout the Bering Sea the biological effect on pCO2 was more than twice as important as temperature, and contributions of other effects were relatively small. Furthermore, satellite observations demonstrate that the spring phytoplankton bloom had a delayed effect on summer pCO2 but that the influence of this biological event varied regionally; it was more significant on the continental slope, with a later bloom, than that on the shelf with an early bloom. Overall, the MeSAA algorithm was not only able to estimate pCO2 in the Bering Sea for the first time, but also provided a quantitative analysis of the contribution of various processes that influence pCO2.


Introduction
The Bering Sea is one of the largest and most productive marginal seas in the world [1], with a sea surface area of 2.3 ˆ10 6 km 2 [2] and primary production up to 85-400 gCm ´2¨yr ´1 [3,4].Moreover, the Bering Strait is the only channel of Pacific water flowing into the Arctic Ocean, which transports heat, freshwater, nutrients, and carbon to the Arctic Ocean [5][6][7].The material transport has a great influence on the carbon source/sink pattern in the Arctic region [8][9][10][11][12].Therefore, the Bering Sea, a marginal sea connecting the Pacific Ocean and the Arctic Ocean, is a very important region for the global marine carbon cycle research [1,[13][14][15][16][17][18][19][20].
Air-sea carbon dioxide (CO 2 ) flux in the Bering Sea has been estimated based on measurements of sea surface partial pressure of carbon dioxide (pCO 2 ) and/or extrapolation of simple correlations of pCO 2 and environmental parameters.Takahashi et al. [21] estimated that the climatological annual uptake of CO 2 on the shelf of the Bering Sea was ´3.65 TgCyr ´1 (over an area of 6.94 ˆ10 5 km 2 ).Bates et al. [1] reported that the CO 2 sink on Bering Sea shelf was ´157 ˘35 TgCyr ´1 (area of 5 ˆ10 5 km 2 ) in 2008 using a multiple linear regression method.Based on the underway data obtained from July to September in 2008, Chen et al. [22] estimated the air-sea CO 2 fluxes on the shelf, slope, and central basin of the Bering Sea were ´9.4,´16.3, and ´5.1 mmolm -2 ¨d-1 , respectively, and the total annual uptake of CO 2 in the Bering Sea was ´34 TgCyr ´1 (area of 23.4 ˆ10 5 km 2 ).Cross et al. [23] reported an annual CO 2 sink of ´6.8 TgCyr ´1 by using in situ data from the eastern Bering shelf from April to December during 2008-2012.Fransson et al. [24] found that the pCO 2 exhibited an increasing trend from 1995 to 2001 in the upwelling region of the southern Bering Sea (53-55 ˝N, 174-194 ˝E), and that this area was a carbon source with average annual CO 2 outgassing of 4.7 ˘3.1 TgCyr ´1 (area of 2.9 ˆ10 5 km 2 ).
The aforementioned investigations have analyzed the carbon source/sink pattern in the Bering Sea, with a focus on the eastern part of the Bering Sea.However, the estimated air-sea CO 2 fluxes differed greatly due to limited temporal and spatial coverage of various cruise data; and interpolation results in large uncertainty.Satellite remote sensing has the advantage of observing a large area over time, and has been increasingly applied to the research on air-sea CO 2 flux (e.g., [25][26][27][28]).Many efforts have been made to estimate sea surface pCO 2 by developing linear or nonlinear regressions between pCO 2 and satellite-derived parameters, such as sea surface temperature (SST) (e.g., [25]), chlorophyll a concentration (chla) (e.g., [26]), salinity (e.g., [29]), colored dissolved organic matter (CDOM) (e.g., [30]), and mixed layer depth (e.g., [27]).However, in some complex waters pCO 2 variability was the result of a combination of multiple factors.It was generally difficult to find a straightforward, significant regression between pCO 2 and salinity, temperature, chla, and other parameters in a large area.Parard et al. [31,32] used a method combining self-organizing maps and multiple linear regression to estimate the ocean surface pCO 2 in the Baltic Sea from the remotely sensed sea surface temperature, chlorophyll, colored dissolved organic matter, net primary production, and mixed-layer depth; and the reconstructed satellite-derived pCO 2 had a good correlation with the in situ measurement.The method of self-organizing maps is a subfamily of neural network algorithms used to perform multidimensional classification, and it was a good way to deal with the complex coastal water.To a certain degree, the more relative the parameters used as the inputs, the more accuracy was obtained in the neural network; on the other hand, it is not easy to provide a clear mechanistic explanation of the structure of neurons.In addition to the direct multiple regression and neural network, some studies began to look into the carbonate system which controls the variation of pCO 2 [28,33].Hales et al. [28] proposed a mechanistic method to estimate pCO 2 in the upwelling-dominated North American Pacific coast, which determined empirical relationships between total dissolved inorganic carbon (DIC, a.k.aTCO 2 ), total alkalinity (TA), and satellite data (SST and chla) within various hydro-biological subregions, and then pCO 2 was calculated from the empirically predicted DIC and TA therein.Recently, Bai et al. [33] proposed a mechanistic semi-analytical algorithm (MeSAA) to estimate pCO 2 by which the variation of pCO 2 can be analytically expressed as the sum of individual components controlled by major factors.More specifically, the MeSAA algorithm first analyzes the main mechanisms causing pCO 2 variability, e.g., temperature/thermodynamics, biology, and mixing; and then parameterizes and quantifies the contributions of each factor to pCO 2 , thus estimating overall pCO 2 variability with satellite parameters.The first implementation of the MeSAA algorithm was in the East China Sea, a river-dominated ocean margin, and showed a good performance [33].
The majority of the Bering Sea has predominately oceanic characteristics, except for some northeast coastal areas influenced by river plumes, which are different from the river-dominated East China Sea; specific methods are needed to parameterize controlling factors on pCO 2 .To our best knowledge, there is still no literature referring to the remote sensing of pCO 2 in the Bering Sea.Therefore, in this paper, we applied the MeSAA method for estimating sea surface pCO 2 in the Bering Sea in summertime using satellite data.In Section 2, we give a general overview of the study area and data.In Section 3, we analyze pCO 2 distribution and major controlling factors, parameterize the thermodynamics and biological effects, and then present and validate satellite-derived pCO 2 .Finally, in Section 4, we discuss the MeSAA parameterization method and quantify the contribution of different controlling factors on pCO 2 variability.

Study Area
The Bering Sea is located in the northern part of the Pacific Ocean (51-66 ˝N and 160 ˝E-158 ˝W) (Figure 1a).It is bordered by Siberia in the west, Alaska in the east, and the Aleutian archipelago and the North Pacific Ocean in the south.The southwestern part of the Bering Sea is a basin (including, from east to west, the Kamchatka Basin, the Bowers Basin, and the Aleutian Basin), with depths over 3500 m in 40% of the area.The northeastern part is a shelf area with a depth of less than 200 m, whose area is about half of the Bering Sea.The shelf area has high productivity [3,4], while the basin has relatively low productivity [5].
Remote Sens. 2016, 8, 558 3 of 25 knowledge, there is still no literature referring to the remote sensing of pCO2 in the Bering Sea.Therefore, in this paper, we applied the MeSAA method for estimating sea surface pCO2 in the Bering Sea in summertime using satellite data.In Section 2, we give a general overview of the study area and data.In Section 3, we analyze pCO2 distribution and major controlling factors, parameterize the thermodynamics and biological effects, and then present and validate satellite-derived pCO2.Finally, in Section 4, we discuss the MeSAA parameterization method and quantify the contribution of different controlling factors on pCO2 variability.

Study Area
The Bering Sea is located in the northern part of the Pacific Ocean (51-66°N and 160°E-158°W) (Figure 1a).It is bordered by Siberia in the west, Alaska in the east, and the Aleutian archipelago and the North Pacific Ocean in the south.The southwestern part of the Bering Sea is a basin (including, from east to west, the Kamchatka Basin, the Bowers Basin, and the Aleutian Basin), with depths over 3500 m in 40% of the area.The northeastern part is a shelf area with a depth of less than 200 m, whose area is about half of the Bering Sea.The shelf area has high productivity [3,4], while the basin has relatively low productivity [5].Currents in the Bering Sea are complex (Figure 1a).Pacific water flows from east to west along the southern side of the Aleutian archipelago and merges into the Alaskan Stream, where it then enters the Bering Sea through the channels of the Aleutian archipelago.After the Alaskan Stream enters the Bering Sea, a portion of it flows northwestwards and, together with the Kamchatka Current (KC), flows southwards out of the Bering Sea along the west coast.Another portion of this water mass flows eastwards, together with the Aleutian North Slope Current (ANSC), and joins the Bering Slope Current (BSC) at ~167°W and flows northwards [34].A portion of the BSC flows northwards to form the Anadyr Current (AC), and another portion joins the KC, flows along the west coast and then flows out of the Bering Sea through the Kamchatka Strait [35,36].The ANSC, BSC, and KC constitute a gyre in the basin area [37].Therefore, the Alaskan Stream has great influence on the water properties in the ANSC, southern KC, and the basins [37].The water mass is relatively stable in the central basin Currents in the Bering Sea are complex (Figure 1a).Pacific water flows from east to west along the southern side of the Aleutian archipelago and merges into the Alaskan Stream, where it then enters the Bering Sea through the channels of the Aleutian archipelago.After the Alaskan Stream enters the Bering Sea, a portion of it flows northwestwards and, together with the Kamchatka Current (KC), flows southwards out of the Bering Sea along the west coast.Another portion of this water mass flows eastwards, together with the Aleutian North Slope Current (ANSC), and joins the Bering Slope Current (BSC) at ~167 ˝W and flows northwards [34].A portion of the BSC flows northwards to form the Anadyr Current (AC), and another portion joins the KC, flows along the west coast and then flows out of the Bering Sea through the Kamchatka Strait [35,36].The ANSC, BSC, and KC constitute a gyre in the basin area [37].Therefore, the Alaskan Stream has great influence on the water properties in the ANSC, southern KC, and the basins [37].The water mass is relatively stable in the central basin due to little influence from the currents and water exchange [35].The eastern shelf is influenced by the inputs from the Yukon River and the Kush Yanukovich River at 63 ˝N and 60 ˝N, respectively.Thus, the low-salinity Alaskan Coastal Current (ACC) is primarily from the mixed water mass of northeastern Pacific water and the river plumes [37][38][39].
Sea ice is a non-negligible factor for pCO 2 variability in the shelf regions of the Bering Sea, and its inter-annual variability causes different seawater properties between the south and north parts of the shelf [40].The ice reaches the maximum in February to March, and melts from July to October [41].The thickness, areal coverage, and southernmost extent of sea ice vary on multiple timescales [42][43][44].
Cross et al. [23] reported that the biological contribution to the carbon sink on the eastern shelf of the Bering Sea from April to December was 3.36 TgCyr ´1.They also found this area was a weak carbon source in autumn, but it was a sink from April to May, accounting for approximately one third of the annual carbon sink.Thus, the remaining summer flux was the greatest portion of the annual flux (~55%).In summer, there is generally an intense seasonal thermocline existing within a depth range of approximately 50-100 m (e.g., [2,45,46]); and the surface water is less influenced by vertical mixing in the basin region due to stratification, resulting in temperature and biology as the two main processes influencing pCO 2 .Thus, as a first step of pCO 2 -MeSAA algorithm development for the Bering Sea, we focus only on the summer season (July through September).

Data and Methods
In situ data included underway temperature, salinity, and pCO 2 in the Bering Sea during 2000-2013, downloaded from the Carbon Dioxide Information Analysis Center (CDIAC) [47].Data distribution is most abundant from May through September (Figure 1b), and were the richest in 2010 including the data from the basin, slope and shelf areas (Figure 1c).In addition, dissolved inorganic carbon (DIC), total alkalinity (TA), and apparent oxygen utilization (AOU) in August 2004 were also obtained from CDIAC [48].Individual cruises and data contributors are listed in Table 1.The satellite data were eight-day and monthly composite MODIS/Aqua products of SST and chla, with a spatial resolution of 4 km, and were obtained from the NASA ocean color website [49].Since the CDIAC dataset did not contain chla data, we matched the eight-day composite satellite chla with the CDIAC underway data in the same time and space windows.
Partial pressure of atmospheric CO 2 (pCO 2(air) ) was calculated from the CO 2 mole fraction of dry air (xCO 2 ), temperature, sea surface vapor pressure (pH 2 O), and sea-level total pressure (SLP) using Equations ( 1) and ( 2 The air-sea CO 2 flux (F) was calculated using the following equations: where K is the gas exchange velocity, and U 10 is monthly average wind speed at 10 m; Sc is the Schmidt number [50]; B is the solubility of CO 2 , a function of temperature and salinity [51]; and pCO 2(sea) is the sea surface pCO 2 .Since we did not aim to show the amount and distribution of air-sea CO 2 flux in this paper, we only used the traditional equation for K calculation (Equation ( 3)), and took the flux as one of the controlling factors on pCO 2 variation.We also calculate DIC accumulation (∆DIC air-sea ) from winter to the summer resulting from the invasion of air-sea CO 2 flux, by using the iteration method [12].First, we obtain the initial TA, DIC, salinity, and temperature in winter from cruise measurement and calculated the initial pCO 2 using CO2SYS (MatLab version [52,53]).Next, a new pCO 2 at the temperature in the next month is calculated with the temperature rise from March to July obtained from the climatological monthly satellite SST data.For the respective ∆pCO 2 (air-sea) , the air-sea CO 2 flux was calculated with Equations (3) and (4), and this amount of CO 2 for a monthly time step is then added to the DIC inventory in the mixed layer.Using the new DIC and constant TA, a new pCO 2 is calculated for the next time step, and so on.The detail value in the calculation will be discussed in Section 4.1.
The re-equilibrium time (τ) triggered by the air-sea CO 2 flux can be calculated with the following equations [54][55][56]: where H is the mixed layer depth; r ř CO 2 s{rCO 2paqq s denotes the ratio of DIC to free CO 2 in seawater, which is generally set as 100 [57]; and R is the Revelle Factor, characterizing the change in DIC caused by the change in (CO 2 ) in seawater and generally ranges from 8 to 15 [58,59].
For the above calculations, SLP data were obtained from the NOAA/ESRL Physical Sciences Division NCEP/NCAR Reanalysis data, with a resolution of 2.5 ˝ˆ2.5 ˝ [60].xCO 2 (ppm) data were obtained from the NOAA/ESRL Global Monitoring Division Carbontracker Project [61], with a spatial resolution of 3 ˝ˆ2 ˝.Climatological monthly average mixed layer depth data are from the website [62], with a spatial resolution of 2 ˝ˆ2 ˝.Monthly average wind speed data (WindSat data) are from the website [63], with a spatial resolution of 0.25 ˝ˆ0.25 ˝.

Correlations between Underway pCO 2 and Related Parameters
First, we wanted to examine the validity of simple regression methods in the Bering Sea.We analyzed the correlations between pCO 2 and temperature, salinity, and chla using all underway data from July to September.We found that pCO 2 did not show a relationship with temperature or salinity (Figure 2).In the basin pCO 2 roughly exhibited a decreasing relationship with increasing temperature (R 2 = 0.46) in the basin (Figure 2a), while the root-mean-square error (RMSE) in the basin was as high as 41.9 µatm.In the shelf area, relationship between pCO 2 and temperature was messy (R 2 = 0.16) (Figure 2c).On the whole, there was no apparent relationship between pCO 2 and temperature.pCO 2 had a weak positive relationship with salinity in the basin area (R 2 = 0.51, RMSE = 36.7 µatm, p = 0.0016) (Figure 2b), while there was no apparent relationship on the shelf (R 2 = 0.03) (Figure 2d).On the shelf, data with salinity ranging from 25 to 30 and the corresponding high chla (>2.5 mg/m 3 ) were located in the area influenced by the Yukon River plume (167.3 ˝W, 62.8 ˝N).On the whole, pCO 2 did not exhibit a decreasing relationship with increasing chla as expected, with R 2 lower than 0.2.In summary, the relationships between these parameters were not apparent and the influences of various mechanisms on pCO 2 were variable in different areas.Thus, it was difficult to establish an empirical algorithm to retrieve pCO 2 from the temperature, salinity, and chla using regression relationships.

Correlations between Underway pCO2 and Related Parameters
First, we wanted to examine the validity of simple regression methods in the Bering Sea.We analyzed the correlations between pCO2 and temperature, salinity, and chla using all underway data from July to September.We found that pCO2 did not show a relationship with temperature or salinity (Figure 2).In the basin pCO2 roughly exhibited a decreasing relationship with increasing temperature (R 2 = 0.46) in the basin (Figure 2a), while the root-mean-square error (RMSE) in the basin was as high as 41.9 μatm.In the shelf area, relationship between pCO2 and temperature was messy (R 2 = 0.16) (Figure 2c).On the whole, there was no apparent relationship between pCO2 and temperature.pCO2 had a weak positive relationship with salinity in the basin area (R 2 = 0.51, RMSE = 36.7 μatm, p = 0.0016) (Figure 2b), while there was no apparent relationship on the shelf (R 2 = 0.03) (Figure 2d).On the shelf, data with salinity ranging from 25 to 30 and the corresponding high chla (>2.5 mg/m 3 ) were located in the area influenced by the Yukon River plume (167.3°W,62.8°N).On the whole, pCO2 did not exhibit a decreasing relationship with increasing chla as expected, with R 2 lower than 0.2.In summary, the relationships between these parameters were not apparent and the influences of various mechanisms on pCO2 were variable in different areas.Thus, it was difficult to establish an empirical algorithm to retrieve pCO2 from the temperature, salinity, and chla using regression relationships.

The MeSAA Algorithm for pCO2
Based on the principle of the MeSAA method proposed by Bai et al. [33], sea surface pCO2 variation (△pCO2) can be expressed as the sum of individual pCO2 contributions associated with each process or controlling factor ( ∂ pCO2(factor-n)):

The MeSAA Algorithm for pCO 2
Based on the principle of the MeSAA method proposed by Bai et al. [33], sea surface pCO 2 variation (∆pCO 2 ) can be expressed as the sum of individual pCO 2 contributions associated with each process or controlling factor (BpCO 2(factor-n) ): where ∆V factor-n represents an individual process that varies independently of the other processes and ε is the residual error.In general, the major factors controlling pCO 2 variability are identified as: (1) the temperature-dependent thermodynamic effect (therm); (2) mixing between water masses with different values of marine carbonate parameters (mix); (3) the biological effect (bio); and (4) air-sea exchange (air-sea).
The surface water pCO 2 is linked directly to DIC, TA, and/or pH in the marine carbonate system in nonlinear relationships.The relationships among pCO 2 , DIC, TA, and pH are fully mathematically defined.As DIC and TA are measured in total concentration, they are mass-conservative during mixing processes; and DIC is directly linked by stoichiometry to nutrient removal/release and O 2 production or consumption during biological processes.Thus, many studies ( [33], and references therein), especially marine carbonate models, calculate pCO 2 from DIC and TA variations by estimating contributions to DIC by specific processes of mixing, biological effect and air-sea CO 2 exchange.Because our desired final product is satellite-derived DIC and TA, we chose to estimate pCO 2 directly to avoid error propagation from estimations of satellite-derived DIC and TA, even though our MeSAA approach was built upon the principle of DIC and TA changes and we also used carbonate system calculations as part of our approach.Thus, we assume ∆pCO 2 is analytically expressed as the sum of contribution from several independent processes, and the critical issue of MeSAA is how to derive the analytical expressions of each individual pCO 2 components using satellite-derived variables.
As described earlier, thermodynamics and biology are the two main factors influencing sea surface pCO 2 in the Bering Sea in summertime.Therefore, we implement the MeSAA algorithm in this paper as follows: first, we need to identify a water mass where temperature is relatively stable and pCO 2 is minimally influenced by mixing and biology; this water mass is then used as a reference (pCO 2(o) ).Second, we develop parameterization methods to quantify the individual thermodynamic (∆pCO 2(therm) ) and biological effects (∆pCO 2(bio) ), and, finally, add them to pCO 2(o) as in Equation ( 8 Since the July 2010 cruise was the only cruise that passed through the basin, slope, and shelf areas (Figure 3a), we used this dataset to analyze the properties of the water masses in different regions, in order to determine the reference water mass and the pCO 2(o) term.We divided the different regions into four zones along the cruise track according to latitudinal distribution of underway salinity and pCO 2 features (Figure 3b).Region 1: in the southern part of the basin (54-57 ˝N), the water masses were relatively stable, and pCO 2 (367.5 ˘18.0 µatm), salinity (32.98 ˘0.06), and temperature (7.57˘0.16 ˝C) varied very little with latitude.Region 2: in the northern part of the basin (57-59.5 ˝N), pCO 2 decreased with decreasing salinity and also was negatively related to temperature.The water in the northern basin might be influenced by mixing with the northwestwards flowing BSC with relatively low salinity and low pCO 2 .Region 3: in the slope area (59.5-61 ˝N), salinity, and pCO 2 reached their minima due to high biological production.According to Springer et al. [3], the area near 61 ˝N was located in the "Bering Sea Green Belt".This area was at the shelf-slope front, and the shelf-basin exchange processes, such as tidal mixing and eddies, provided high nutrient concentrations to the surface layer, resulting in primary production here up to 175-275 gCm ´2¨y ´1, 60% higher than that in the adjacent shelf area.Region 4: in the shelf area (above 61 ˝N), salinity and pCO 2 increased slightly compared with the slope area, however, pCO 2 was still relatively low (below 280 µatm) compared with the other regions.In this region we also found that there was a weak positive relationship between pCO 2 and salinity, which might be influenced by mixing with the northwards flowing BSC.The water mass in the southern basin was relatively stable as suggested by in situ pCO2, salinity, and temperature variability; chla was at a relative low level (below 0.5 mg/m 3 ) compared with other areas, seen from the climatological average (2002-2014) satellite chla image in July (Figure 3d in box C2).Therefore, we selected the water mass in the southern basin as the reference water mass (55.7-57°N, 171.4-174.2°W),shown in box A of Figure 3a.In a total of 238 sets of data, temperature (7.7 ± 0.1 °C) and salinity (32.98 ± 0.02) in the selected reference water changed very little (Figure 3c), indicating the limited influence of water mixing.Furthermore, the average temperature of the reference water mass (7.7 ± 0.1 °C) was within the temperature range (6 to 12 °C) of the basin and slope areas in summer, thus errors resulting from the wide temperature range of pCO2 normalization by the thermodynamic equation should be minimal [33,64].We also calculated the monthly composite MODIS/Aqua SST during 2002-2014 within this area, and its average standard deviation (SD) for individual months during these 13 years was within ± 0.9 °C, indicating a small inter-annual change in temperature.Overall, the selected water mass was relatively stable and could be used as a reference water mass for the remote sensing algorithm of pCO2.
The temperature-dependent thermodynamic effect on sea surface pCO2 has an exponential relationship of approximately 4.23%/°C [21, 59,65].Note that the slope of 0.0423 was derived for the global ocean, and this can be confirmed with carbonate system calculations, although the slope varied slightly with TA, DIC, salinity and temperature [33].Thus, we can analytically calculate the pCO2 at a certain temperature (T) by the exponential relationship as in Equation (9) [21, 59,65].Note that pCO2(T) in Equation (9) represents the first two items in the Equation (8), which combined the value of reference water (pCO2(o)) and the temperature-dependent thermodynamic effect.
Briefly, we used the average values in the reference water mass as pCO2(o) = 381.8± 5.1 and temperature T(o) = 7.7 ± 0.1 (Table 2), as well as the temperature effect on pCO2 under a constant TA, DIC, and salinity [65].pCO2(o) was an empirical value in practice, and we will prove its stability through further marine carbonate system analysis in Section 4.1.The water mass in the southern basin was relatively stable as suggested by in situ pCO 2 , salinity, and temperature variability; chla was at a relative low level (below 0.5 mg/m 3 ) compared with other areas, seen from the climatological average (2002-2014) satellite chla image in July (Figure 3d in box C2).Therefore, we selected the water mass in the southern basin as the reference water mass (55.7-57 ˝N, 171.4-174.2˝W), shown in box A of Figure 3a.In a total of 238 sets of data, temperature (7.7 ˘0.1 ˝C) and salinity (32.98 ˘0.02) in the selected reference water changed very little (Figure 3c), indicating the limited influence of water mixing.Furthermore, the average temperature of the reference water mass (7.7 ˘0.1 ˝C) was within the temperature range (6 to 12 ˝C) of the basin and slope areas in summer, thus errors resulting from the wide temperature range of pCO 2 normalization by the thermodynamic equation should be minimal [33,64].We also calculated the monthly composite MODIS/Aqua SST during 2002-2014 within this area, and its average standard deviation (SD) for individual months during these 13 years was within ˘0.9 ˝C, indicating a small inter-annual change in temperature.Overall, the selected water mass was relatively stable and could be used as a reference water mass for the remote sensing algorithm of pCO 2 .
The temperature-dependent thermodynamic effect on sea surface pCO 2 has an exponential relationship of approximately 4.23%/ ˝C [21, 59,65].Note that the slope of 0.0423 was derived for the global ocean, and this can be confirmed with carbonate system calculations, although the slope varied slightly with TA, DIC, salinity and temperature [33].Thus, we can analytically calculate the pCO 2 at a certain temperature (T) by the exponential relationship as in Equation (9) [21, 59,65].Note that pCO 2(T) in Equation (9) represents the first two items in the Equation (8), which combined the value of reference water (pCO 2(o) ) and the temperature-dependent thermodynamic effect.pCO 2pTq " pCO 2pOq ˆexpp0.0423ˆpT ´TpOq qq (9) Briefly, we used the average values in the reference water mass as pCO 2(o) = 381.8˘5.1 and temperature T (o) = 7.7 ˘0.1 (Table 2), as well as the temperature effect on pCO 2 under a constant TA, DIC, and salinity [65].pCO 2(o) was an empirical value in practice, and we will prove its stability through further marine carbonate system analysis in Section 4.1.

Quantification of Biological Effect
The biological effect on pCO 2 variability is an indirect and complex process that changes the DIC concentration through various processes, such as photosynthesis/respiration and/or CaCO 3 production/dissolution, and thus changes pCO 2 .Direct and accurate analytic expressions of these processes are still an open question [27,64].For practical satellite data processing we used a similar method to that in Bai et al. [33], which uses an integral to express pCO 2 drawdown with an increase of chla.In an area mainly controlled by biological activities, it can be assumed that pCO 2 and chla have the following relationship [33]: where ε represents the contribution from other factors independent of biological effects.The influence of the biological effect on pCO 2 exhibits an accumulation effect with time, so the ∆pCO 2(bio) can be expressed with the following equation [33]: where Chla n is the sea surface chlorophyll concentration, and Chla 0 is the background value.
A subset of data controlled mainly by the biological effect, with little influence by terrestrial input and mixing effect, needed to be determined to regress the coefficient B in Equation (10).Due of not having in situ chla data and few available chla data in the daily satellite products due to cloud cover, we matched the eight-day composite satellite chla data with underway data.Through analyzing all available match-up data in the spring and summer, we found that pCO 2 and chla had significant negative relationships during May to June 2010, which were located on slope and shelf areas (Figure 4a,b); therefore, we used this subset of data for further analysis.The thermodynamic effect should first be eliminated in order to remove this apparent effect when analyzing the other processes.We normalized pCO 2 (NpCO 2 ) to the average temperature (2.3 ˝C) of this subset (Figure 4c).There was a significant negative relationship between NpCO 2 and chla with R 2 = 0.757 (Figure 4c); furthermore, the consistent relationship in May and June indicated that pCO 2 drawdown caused by per chla production was relatively consistent in this area.Thus, we obtained the coefficient B in Equation (10) as: NpCO 2 " ´217.62 ˆlogpchlaq `300.13.
Consequently, pCO 2 variability caused by the biological effect can be expressed as Equation ( 13).The chla 0 in Equation ( 11) was taken empirically as 0.1 mg/m 3 , according to the generally low chla level in the basin area.
dpchlaq " ´217.62 ˆrlogpchla n q ´logp0.1qs Remote Sens. 2016, 8, 558 10 of 25 Consequently, pCO2 variability caused by the biological effect can be expressed as Equation ( 13).The chla0 in Equation ( 11) was taken empirically as 0.1 mg/m 3 , according to the generally low chla level in the basin area.

Results and Validation
We generated monthly composite satellite SST and chla data for July, August, and September of 2010 and 2011 based on the inversion models (Equations ( 8), (9), and ( 13)).The accuracy of satellitederived chla might be not good at the Bering Sea, and showed systematic overestimation for the SeaWiFS OC4V4 algorithm (i.e., [66,67]).The ocean optical community must improve the chla algorithm in some special sea, like the Bering Sea.However, such systematic deviation of satellitederived chla had been eliminated from the pCO2-MeSAA algorithm development.Because of not having in situ chla data, we matched the MODIS-derived chla data with underway pCO2 data.Through the regression, we get the key coefficient of ∂(pCO2@bio)/∂(log(chla)) and finally, calculate the pCO2 drawdown caused by phytoplankton carbon uptake with MODIS-derived chla (Equation ( 13)).In these procedures, the systematic deviation has less influence on the estimated pCO2 in MeSAA algorithm.
From the images of satellite-derived pCO2 (Figure 5), pCO2 in the basin was approximately 250-360 μatm, and lower than atmospheric pCO2 (387 μatm) in most of the Bering Sea.pCO2 on the slope was lower than that in the basin, ranging from 200 to 300 μatm; and pCO2 in the shelf area ranged from 100 to 250 μatm.pCO2 varied both spatially and temporally over the Bering Sea during the years studied.Due to high SST in August (10.38 ± 0.29 °C in 2010 and 10.48 ± 0.36 °C in 2011) pCO2 was generally higher in August (241.3± 13.1 μatm in 2010 and 252.8 ± 14.1 μatm in 2011) than in July and September, especially in the basin area with low chla.High chla concentration in the slope and shelf

Results and Validation
We generated monthly composite satellite SST and chla data for July, August, and September of 2010 and 2011 based on the inversion models (Equations ( 8), (9), and ( 13)).The accuracy of satellite-derived chla might be not good at the Bering Sea, and showed systematic overestimation for the SeaWiFS OC4V4 algorithm (i.e., [66,67]).The ocean optical community must improve the chla algorithm in some special sea, like the Bering Sea.However, such systematic deviation of satellite-derived chla had been eliminated from the pCO 2 -MeSAA algorithm development.Because of not having in situ chla data, we matched the MODIS-derived chla data with underway pCO 2 data.Through the regression, we get the key coefficient of B(pCO2@bio)/B(log(chla)) and finally, calculate the pCO 2 drawdown caused by phytoplankton carbon uptake with MODIS-derived chla (Equation ( 13)).In these procedures, the systematic deviation has less influence on the estimated pCO 2 in MeSAA algorithm.
From the images of satellite-derived pCO 2 (Figure 5), pCO 2 in the basin was approximately 250-360 µatm, and lower than atmospheric pCO 2 (387 µatm) in most of the Bering Sea.pCO 2 on the slope was lower than that in the basin, ranging from 200 to 300 µatm; and pCO 2 in the shelf area ranged from 100 to 250 µatm.areas, as well as some patchiness of algal blooms in the basin, resulted in under-saturated pCO2 with respective to the atmosphere, with pCO2 as low as ~100 μatm.For validation of satellite-derived pCO2, we matched up satellite data with underway measurements, and narrowed the temporal and spatial gap between these two datasets as much as possible.However, the MODIS/Aqua data (a spatial resolution of 4 km) had few available data in daily products due to the heavy cloud coverage.Therefore, we used the eight-day composite satellite data.There were only three cruises from July to September in our datasets-July 2010, August 2011, and September 2010-and the available data matched up primarily in the central basin and near the slope region (Figure 6a).Since the data in the central basin from July 2010 had been used to calculate the reference values (Box A in Figures 3a and 6a), we exclude this data in the validation in Figure 6b.In general, satellite-derived pCO2 was consistent with the underway pCO2, especially in the basin (Figure 6b).The RMSE for the entire validation data set in eight-day composite MODIS was 9.87% (N = 1274, SD = 35.11μatm).For the central basin, the average RMSE was 4.37%, ranging from 0.04% to 14.95% (N = 1052, SD = 17.67 μatm).For the slope, the average RMSE was 35.9%, ranging from 22.81% to 55.61% (N = 222, SD = 74.8μatm), and the satellite-derived pCO2 was about 50 μatm higher than the underway pCO2 (data in box D in Figure 6b, in which salinity was below 32.7).This is an interesting phenomenon and will be discussed in Section 4.2.
In order to check the validation on a larger scale, we also used the monthly-composite satellite data, and compared the results with underway data in general trends (Figure 7).The MeSAA prediction yielded the best consistency with data in the basin area (Figure 7d), as well as the offshore area near the west coast of the Bering Sea (Figure 7b).The data used for validation were independent of those used to develop the algorithm, thus the current comparisons have proved good reliability of MeSAA algorithm in the majority of the Bering Sea.In the shelf area, the satellite-derived pCO2 near the 50 m isobaths was, however, significantly underestimated (Figure 7c).Un-predicted high pCO2 in the shelf area might be a result of two possibilities: (1) the influence of the Yukon and Kush Yanukovich River plumes, where in situ pCO2 was up to ~390 μatm due to the re-mineralization of terrestrial organic matter; and (2) the effect of coccolithophorid blooms.We will discuss these in Section 4.2.For validation of satellite-derived pCO 2 , we matched up satellite data with underway measurements, and narrowed the temporal and spatial gap between these two datasets as much as possible.However, the MODIS/Aqua data (a spatial resolution of 4 km) had few available data in daily products due to the heavy cloud coverage.Therefore, we used the eight-day composite satellite data.There were only three cruises from July to September in our datasets-July 2010, August 2011, and September 2010-and the available data matched up primarily in the central basin and near the slope region (Figure 6a).Since the data in the central basin from July 2010 had been used to calculate the reference values (Box A in Figures 3a and 6a), we exclude this data in the validation in Figure 6b.In general, satellite-derived pCO 2 was consistent with the underway pCO 2 , especially in the basin (Figure 6b).The RMSE for the entire validation data set in eight-day composite MODIS was 9.87% (N = 1274, SD = 35.11µatm).For the central basin, the average RMSE was 4.37%, ranging from 0.04% to 14.95% (N = 1052, SD = 17.67 µatm).For the slope, the average RMSE was 35.9%, ranging from 22.81% to 55.61% (N = 222, SD = 74.8µatm), and the satellite-derived pCO 2 was about 50 µatm higher than the underway pCO 2 (data in box D in Figure 6b, in which salinity was below 32.7).This is an interesting phenomenon and will be discussed in Section 4.2.
In order to check the validation on a larger scale, we also used the monthly-composite satellite data, and compared the results with underway data in general trends (Figure 7).The MeSAA prediction yielded the best consistency with data in the basin area (Figure 7d), as well as the offshore area near the west coast of the Bering Sea (Figure 7b).The data used for validation were independent of those used to develop the algorithm, thus the current comparisons have proved good reliability of MeSAA algorithm in the majority of the Bering Sea.In the shelf area, the satellite-derived pCO 2 near the 50 m isobaths was, however, significantly underestimated (Figure 7c).Un-predicted high pCO 2 in the shelf area might be a result of two possibilities: (1) the influence of the Yukon and Kush Yanukovich River plumes, where in situ pCO 2 was up to ~390 µatm due to the re-mineralization of terrestrial organic matter; and (2) the effect of coccolithophorid blooms.We will discuss these in Section 4.2.In summary, the data we used to validate the satellite-derived pCO2 were totally independent data that had not been used in the algorithm development.Comparison results showed good agreement between satellite-derived pCO2 and underway data, especially in the central basin.With the mechanistic principle of MeSAA algorithm, the deviation between the satellite results and in situ data can highlight unconsidered processes; such deviations help to further understand the   In summary, the data we used to validate the satellite-derived pCO2 were totally independent data that had not been used in the algorithm development.Comparison results showed good agreement between satellite-derived pCO2 and underway data, especially in the central basin.With the mechanistic principle of MeSAA algorithm, the deviation between the satellite results and in situ data can highlight unconsidered processes; such deviations help to further understand the In summary, the data we used to validate the satellite-derived pCO 2 were totally independent data that had not been used in the algorithm development.Comparison results showed good agreement between satellite-derived pCO 2 and underway data, especially in the central basin.With the mechanistic principle of MeSAA algorithm, the deviation between the satellite results and in situ data can highlight unconsidered processes; such deviations help to further understand the mechanism and parameterization of new components that influence pCO 2 variability.Analysis of these issues will be presented in Sections 4.2 and 4.3.

Contribution of Different Controlling Factors to the Variation of pCO 2
Since the principle of MeSAA algorithm was the sum of individual pCO 2 variations associated with each process, we could identify the contribution of different controlling factors: NpCO 2pseaq " pCO 2pseaq@T obs ˚ep0.0423˚pTpoq´Tobs qq (14) ∆pCO 2pthermq " ∆pCO 2psea´airq ´N∆pCO 2psea´airq " ppCO 2pseaq ´pCO 2pairq q ´pNpCO 2pseaq ´pCO 2pairq q (15) B ˚" N∆pCO 2psea´airq " N pCO 2pseaq ´pCO 2pairq ( 16) To distinguish the temperature effect from the biological effect, Cross et al. [23] normalized pCO 2 (NpCO 2 ) obtained in the eastern Bering Sea shelf from April to December during 2008-2012 to the annual average temperature (Equation ( 14)).This calculation removes the localized effects of warming and cooling on pCO 2(sea) , such that the temperature effect at any one time relative to the annual average can be estimated as the difference between pCO 2(sea) and NpCO 2(sea) (Equation ( 15)) [23].Accordingly, NpCO 2(sea) variability should be the result of any additional process that modifies pCO 2(sea) other than thermodynamic effects, such as biological CO 2 uptake, air-sea exchange, and vertical/lateral transport of carbon [23].To simplify the calculation, Takahashi et al. [58], Bates et al. [1], and Cross et al. [23] all regarded the variations in NpCO 2(sea) as the result of biogeochemical effect (defined as B* in Equation ( 16)), which includes the other processes.For our MeSAA algorithms, however, we can quantify the contribution of biological uptake separately (Equation ( 13)), so the contribution of other processes (∆pCO 2(other) ) can be further separated from B * (Equation ( 17)).
Thus, using Equations ( 13)-( 17), we obtained influences of three types of processes that affected pCO 2 : ∆pCO 2(therm), ∆pCO 2(bio) , and ∆pCO 2(other) , in the Bering Sea during the summer with satellite data (shown in Figure 8).Here, we used the reference temperature T (0) = 7.7 ˝C in Equation (14).From July to September, ∆pCO 2(bio) was the primary influence on pCO 2 , and decreased it 100-200 µatm; while the contributions of ∆pCO 2(therm) and ∆pCO 2(other) were relatively low, each about ˘50 µatm.As the sea surface pCO 2 was generally lower than the atmospheric pCO 2 in summer, there was an atmospheric CO 2 influx into the surface ocean that increased pCO 2 .In addition, vertical mixing also increased the ocean pCO 2 due to high DIC water from the bottom layer; therefore, ∆pCO 2(other) were generally positive, with a value of 0 to 30 µatm in most of the Bering Sea.
The ∆pCO 2(bio) pattern was neither spatially nor temporally uniform in the Bering Sea (Figure 8).Relatively low ∆pCO 2(bio) (´50-´100 µatm) in the central basin corresponded to relatively low chla (0.6 mg/m 3 ) (Figures 3d and 8).∆pCO 2(therm) declined gradually from west to east (Figure 8), which followed the same pattern as SST.∆pCO 2(therm) and ∆pCO 2(other) in the majority of the Bering Sea was close to equilibrium, with a slight decrease from the southwest to the northeast, especially on the shelf, except in the Kamchatka Basin.
∆pCO 2(therm) was remarkably high in the Kamchatka Basin in the west of the Bering Sea (on the west of 170 ˝E), reaching 60 µatm in August 2011; both chla (~0.9 mg/m 3 ) and SST (~11.5 ˝C) in the Kamchatka Basin were higher than the central basin, therefore the contributions of ∆pCO 2(bio) and ∆pCO 2(therm) here were larger than in the other sub-regions; however the influence of these two effects was opposite and thus compensated for each other.Thus, the contribution of ∆pCO 2(other) was positive, but close to neutral.∆pCO 2(other) was greater in the Kamchatka Basin than in the other sub-regions.The complex currents in the Kamchatka Basin, with the Kamchatka Current flowing out through the Kamchatka Strait, Pacific water flowing in through the Near Strait, and anticyclone eddies [35,37], cause high sea surface pCO 2 variability.In general, the contribution of the biological effect in the summer on sea surface pCO 2 variability was roughly more than twice as important as the temperature effect.Other processes increased the sea surface pCO 2 but their influences were similar to, or slightly less than, the influence of the temperature effect.

Stability of Reference Water Mass
The reference water mass for the MeSAA model in the Bering Sea was chosen in the central basin, where currents and mixing had a minimal influence on sea surface pCO2.The location of the reference water mass was close to the entrance of the Alaskan Stream flowing into the Bering Sea from the Near Strait.Thus, the stability of the reference water mass would be influenced by the Alaskan Stream water properties and flow rate.Favorite [68] suggested that the flow rate of the Alaskan Stream from the Near Strait was constant (10 × 10 6 m 3 •s −1 ).Onishi et al. [69] also found that the Alaskan Stream was stable, with relatively high homogeneity in the vertical direction compared with the surrounding water masses [70].Therefore, we assume that the variation of the Alaskan Stream had little influence on the properties of the reference water mass.
Next we discuss carbonate system variability of the reference water mass for its source and seasonal evolution (Equations ( 18) and ( 19)).In MeSAA algorithms, our defined reference pCO2(o) in summer could be regarded as the result of well-mixed winter water (pCO2(winter)) going through a series of processes, such as warming, biological alteration (∆pCO2(bio)), and air-sea exchange (∆pCO2(airsea)), to evolve into the summer status (pCO2(summer)) (Equation ( 18)).We now try to quantify this seasonal evolution using marine carbonate system calculations; the comparison between selected pCO2(o) and calculated pCO2(summer) can check the stability of reference water mass.TA is assumed to be constant in this area based on the stable Alaskan Stream and limited influence from currents mixing in the central basin, thus variations in pCO2 will depend on DIC change (Equation ( 19)).

Stability of Reference Water Mass
The reference water mass for the MeSAA model in the Bering Sea was chosen in the central basin, where currents and mixing had a minimal influence on sea surface pCO 2 .The location of the reference water mass was close to the entrance of the Alaskan Stream flowing into the Bering Sea from the Near Strait.Thus, the stability of the reference water mass would be influenced by the Alaskan Stream water properties and flow rate.Favorite [68] suggested that the flow rate of the Alaskan Stream from the Near Strait was constant (10 ˆ10 6 m 3 ¨s´1 ).Onishi et al. [69] also found that the Alaskan Stream was stable, with relatively high homogeneity in the vertical direction compared with the surrounding water masses [70].Therefore, we assume that the variation of the Alaskan Stream had little influence on the properties of the reference water mass.
Next we discuss carbonate system variability of the reference water mass for its source and seasonal evolution (Equations ( 18) and ( 19)).In MeSAA algorithms, our defined reference pCO 2(o) in summer could be regarded as the result of well-mixed winter water (pCO 2(winter) ) going through a series of processes, such as warming, biological alteration (∆pCO 2(bio) ), and air-sea exchange (∆pCO 2(air-sea) ), to evolve into the summer status (pCO 2(summer) ) (Equation ( 18)).We now try to quantify this seasonal evolution using marine carbonate system calculations; the comparison between selected pCO 2(o) and calculated pCO 2(summer) can check the stability of reference water mass.TA is assumed to be constant in this area based on the stable Alaskan Stream and limited influence from currents mixing in the central basin, thus variations in pCO 2 will depend on DIC change (Equation ( 19)).pCO 2psummerq " pCO 2pwinterq ´∆pCO 2pbioq ´∆pCO 2pair´seaq (18) DIC psummer@sur f aceq " DIC pwinter@sur f aceq ´∆DIC air´sea " rDIC psummer@temp´min´layerq ´∆DIC bio s ´∆DIC air´sea First we determined DIC and TA values in winter surface water based on field measurements.In order to avoid differences between cruises, we estimated winter values from summer observations.We were able to do this as there was a temperature minimum layer (below 3 ˝C) at 75-150 m in summer under the thermocline in the Bering Sea basin; this layer was regarded as a remnant of the winter surface layer [45,46,71,72].Therefore, we used observations from the temperature minimum layer in summer as winter values; however, it should be pointed out that these values still include biologically produced DIC (∆DIC bio ), which should be removed to obtain a final winter value (Equation ( 19)).In the spring and summer, the organic carbon produced by the phytoplankton photosynthesis in the upper layer sinks to the lower layer, and is decomposed into inorganic carbon through biological respiration (∆DIC bio ).This process was recorded by the apparent oxygen utilization (AOU), and could be converted to carbon concentration based on the Redfield Ratio [73] (∆DIC bio = 106*AOU/138).There was only one set of inorganic carbon observations obtained during the summer in the PACIFICA database: data from an east-to-west section at 57 ˝N in August 2004 (Figure 9a).The minimum temperature (2-3 ˝C) occurred at the depth of 100 m at all stations during this cruise (Figure 9b), indicating that the properties of the water mass in this area were constant and stable.We used data from Station 38 (Table 2), which was closest to the position of our selected reference water mass.As AOU = 35.98µmol/kg in the temperature minimum layer, we obtained ∆DIC bio = 106*AOU/138 = 27.7 µmol/kg, thus the well-mixed winter DIC value was 2126.76 µmol/kg (Station 38 summer value minus ∆DIC bio ).With the unchanged TA, pCO 2(winter) was calculated to be 415.1 µatm.Since sea surface pCO 2 in the winter was higher than atmospheric pCO 2 (387.2µatm), CO 2 would have outgassed to the atmosphere.Thus, ∆DIC air-sea from winter to summer (set as March to July) could be calculated with the iteration method (see the description in Section 2.2), and the relative values in the iteration calculation were as follows.The initial conditions including the DIC, TA, temperature, and salinity were set as those in March (winter value) (Table 2); we assumed that the wind speed was constant as the average value during March to July (7.28 m/s from climatologic satellite data) and the temperature raising rate in surface seawater was constant at 0.035 ˝C/d [(T 0 -T winter )/150].The climatological monthly average mixing layer depth was 50.7 m during March to July (Section 2.2).We assume an atmospheric pCO 2 of 387.2 µatm from the average xCO 2 data during March to July (Equations ( 1) and ( 2)).After the iteration calculation from March to July (150 days), we obtained the final ∆DIC air-sea as 20.36 µmol/kg.If there were only air-sea exchange and warming effects on the surface water, the surface DIC in summer decreased to 2106.4 µmol/kg according the T (o) = 7.7 ˝C in July (Equation ( 19)), and then the calculated pCO 2(summer) was 450.9 µatm.Considering that chla in the reference water mass was approximately 0.2 mg/m 3 in July 2010, ∆pCO 2bio was estimated with Equation ( 13) to be 65 µatm, then the final pCO 2(summer) was 385.9 µatm.
In summary, starting with a pCO 2 in winter of 415.1 µatm, and after a series of processes as air-sea exchange and warming (increase by 35.8 µatm)) and biological effect (decrease by 65 µatm), pCO 2(summer) was about 385.9 µatm.This calculated summer value was close to the reference pCO 2(o) (381.8 µatm).
The above calculation used a traditional Redfield C:O ratio of 106:139 [73].Considering the Redfield ratio is various over different biogeochemical conditions (e.g., [74][75][76][77]), we used C:O ratios ranging between 94:131.8 and 134:162 [77] to test its sensitivity.Then, we obtained the calculated well-mixed winter DIC value as 2124.64-2128.74µmol/kg.With the unchanged TA, pCO 2(winter) was calculated to be 408.7-421.3µatm.The surface DIC in summer decreased to 2105.7-2106.9µmol/kg according the T (o) = 7.7 ˝C in July (Equation ( 19)), and finally the calculated pCO 2(summer) was 448.7-452.4µatm.The result (448.7-452.4µatm) considering the Redfield ratio variability was very close to that with the traditional Redfield ratio (450.9 µatm).In the central basin, since the biological level was very slow, ∆DIC bio was also small, thus the variability in Redfield ratio had little influence on the results.This agreement further suggests that our selected reference water mass is derived from well-mixed winter water and can be traced with quantitative processes, and thus the reference values were applicable to estimated satellite-derived sea surface pCO 2 in a long time series.The above estimation is based on simplified processes and is an approximate estimation of changes in carbonate system biogeochemistry.Future work is needed on larger scale satellite data to use the MeSAA algorithm during an entire calendar year.
This agreement further suggests that our selected reference water mass is derived from well-mixed winter water and can be traced with quantitative processes, and thus the reference values were applicable to estimated satellite-derived sea surface pCO2 in a long time series.The above estimation is based on simplified processes and is an approximate estimation of changes in carbonate system biogeochemistry.Future work is needed on larger scale satellite data to use the MeSAA algorithm during an entire calendar year.Since sea surface pCO2 had been drawn down to a very low level at the end of algal blooms, we could estimate the re-equilibrium time (τ ) of CO2 (Equations ( 5) and ( 6) in Section 2.2), which was triggered by the large gradient of air-sea pCO2.In the basin the biological effect was at a relatively low level and the algal bloom from April to June was not significant, therefore having little influence on pCO2 in the summer.On the shelf the calculated re-equilibrium time (τ ) was roughly 40-60 days, calculated with the climatological monthly mixed layer depth from March to June of 15-20 m.The spring bloom occurred in April on the shelf but in June the chla decreased to a relatively small value; then, after a period of 1-2 months, the air-sea exchange compensated for the biological signal and reestablished the air-sea equilibrium in July.Therefore, the influence of the spring bloom on the satellite-derived pCO2 on the shelf was minimal.The adequate re-equilibrium time by air-sea exchange in the central basin and on the shelf could explain the good performance of MeSAA algorithms in these regions (Figure 7).
In contrast to the other regions, slope chla reached its peak in May and still maintained a high level (2.3 mg/m 3 ) in June, and then decreased to 0.35 mg/m 3 in July (Figure 10a).Therefore, only one month was not enough for the re-equilibrium by air-sea exchange, and pCO2 in July still maintained a relatively low value, leading to a typical low chla and low pCO2 phenomenon.Such a delayed influence of the previous spring's biological effect could be a reason for the ~50 μatm overestimation of satellite-derived pCO2 on the slope in July 2010 (Figure 6b).Another possible reason for the overestimation could be due to tidal mixing and upwelling on the slope [80][81][82]; however, these overestimated data had salinity lower than 32.7 (Figure 6b), which would be less possible from the influence of the sub-surface high-pCO2, high salinity (>33 psu) water [46,81,83].
Specific phytoplankton functional groups that affect carbonate chemistry would also influence the pCO2.For example, coccolithophorid blooms have been found to frequently occur on the eastern shelf of the Bering Sea from August to October [84-89]; during the blooms, the cells release CO2 during calcification, consequently elevating the surface seawater pCO2 [88,89].This process offsets the DIC consumption by photosynthesis.On the shelf, we found NpCO2 (normalized to the reference Since sea surface pCO 2 had been drawn down to a very low level at the end of algal blooms, we could estimate the re-equilibrium time (τ) of CO 2 (Equations ( 5) and ( 6) in Section 2.2), which was triggered by the large gradient of air-sea pCO 2 .In the basin the biological effect was at a relatively low level and the algal bloom from April to June was not significant, therefore having little influence on pCO 2 in the summer.On the shelf the calculated re-equilibrium time (τ) was roughly 40-60 days, calculated with the climatological monthly mixed layer depth from March to June of 15-20 m.The spring bloom occurred in April on the shelf but in June the chla decreased to a relatively small value; then, after a period of 1-2 months, the air-sea exchange compensated for the biological signal and re-established the air-sea equilibrium in July.Therefore, the influence of the spring bloom on the satellite-derived pCO 2 on the shelf was minimal.The adequate re-equilibrium time by air-sea exchange in the central basin and on the shelf could explain the good performance of MeSAA algorithms in these regions (Figure 7).
In contrast to the other regions, slope chla reached its peak in May and still maintained a high level (2.3 mg/m 3 ) in June, and then decreased to 0.35 mg/m 3 in July (Figure 10a).Therefore, only one month was not enough for the re-equilibrium by air-sea exchange, and pCO 2 in July still maintained a relatively low value, leading to a typical low chla and low pCO 2 phenomenon.Such a delayed influence of the previous spring's biological effect could be a reason for the ~50 µatm overestimation of satellite-derived pCO 2 on the slope in July 2010 (Figure 6b).Another possible reason for the overestimation could be due to tidal mixing and upwelling on the slope [80][81][82]; however, these overestimated data had salinity lower than 32.7 (Figure 6b), which would be less possible from the influence of the sub-surface high-pCO 2 , high salinity (>33 psu) water [46,81,83].
Specific phytoplankton functional groups that affect carbonate chemistry would also influence the pCO 2 .For example, coccolithophorid blooms have been found to frequently occur on the eastern shelf of the Bering Sea from August to October [84-89]; during the blooms, the cells release CO 2 during calcification, consequently elevating the surface seawater pCO 2 [88,89].This process offsets the DIC consumption by photosynthesis.On the shelf, we found NpCO 2 (normalized to the reference temperature of 7.7 ˝C) and chla had a good positive correlation (Figure 11a) Meanwhile, the deviation between the satellite-derived pCO 2 and in situ data also had a significant negative correlation with chla; it also indicated that the contributions of other processes (not considered here) on pCO 2 would be relative to the biological effects (Figure 11b) Thus, we hypothesize that the other processes not considered might be influenced by coccolithophorid blooms, which could explain the underestimation in MeSAA validation in Figure 7c.Another possibility is the influence of the re-mineralization of terrestrial organic matter, which could result in a pCO 2 increase.We infer this influence from terrestrial sources because (1) the underestimated data corresponded with low salinity (below 31 salinity units in the sub-regions in Figure 7a); and (2) the deviation between the satellite-derived pCO 2 and in situ data also showed a positive relationship with salinity below 31 (Figure 11c).Unfortunately, we do not have the necessary data to clarify these two speculations at this time.Further efforts should be made to understand the complex biological and terrestrial effect on pCO 2 , and to obtain more information from the satellite data beyond chla, such as the identification of special functional groups and the timing of algal blooms (onset, duration, etc.).
Remote Sens. 2016, 8, 558 18 of 25 temperature of 7.7 °C) and chla had a good positive correlation (Figure 11a) Meanwhile, the deviation between the satellite-derived pCO2 and in situ data also had a significant negative correlation with chla; it also indicated that the contributions of other processes (not considered here) on pCO2 would be relative to the biological effects (Figure 11b) Thus, we hypothesize that the other processes not considered might be influenced by coccolithophorid blooms, which could explain the underestimation in MeSAA validation in Figure 7c.Another possibility is the influence of the remineralization of terrestrial organic matter, which could result in a pCO2 increase.We infer this influence from terrestrial sources because (1) the underestimated data corresponded with low salinity (below 31 salinity units in the sub-regions in Figure 7a); and (2) the deviation between the satellitederived pCO2 and in situ data also showed a positive relationship with salinity below 31 (Figure 11c).Unfortunately, we do not have the necessary data to clarify these two speculations at this time.
Further efforts should be made to understand the complex biological and terrestrial effect on pCO2, and to obtain more information from the satellite data beyond chla, such as the identification of special functional groups and the timing of algal blooms (onset, duration, etc.).

The Influence of Mixing
As a first step toward implementing pCO2-MeSAA for the Bering Sea, we chose the summer season to eliminate the influence of vertical mixing caused by the seasonal cooling; however, mixing still occurs during the summer in some areas of the Bering Sea.The main currents in the Bering Sea include the westwards flowing ANSC along with the Aleutian archipelago and the anticlockwise flowing BSC in the slope [34,35].Seasonal upwelling, caused by wind and eddies (due to the difference of currents velocity), on the southeastern slope and southern Bering Sea close to the Aleutian archipelago, would also occur during the summer (e.g., [2,46,80,90,91]) and bring high DIC and nutrient concentrations from subsurface water to the surface layer [24].The slope in the northeast was influenced by the broad, northwestwards flowing BSC and appeared to be a system of eddies rather than a continuous current [80,82], therefore the shelf-basin exchange processes of intensive tidal mixing and eddies likely play an important role in pCO2 variability.Mixing processes supply nutrients to the regions adjacent to the slope and enhance primary production, eventually leading to a biological drawdown of pCO2 [22].The shelf was affected by mixing with the Yukon River and Kush Yanukovich River plumes increasing pCO2 from organic carbon re-mineralization.
In this preliminary study we analyzed underway data in the central basin to access the influence of mixing on pCO2.There were only two cruises in July 2010 and August 2011 that measured pCO2 in the central basin during summer months.Since most of the corresponding satellite chla data in July Figure 11.The relationship between pCO 2 , chla, temperature, and salinity in the shelf in July 2010 (locations were marked as M2 in Figure 7).(a) The relationship between pCO 2 normalized to reference temperature 7.7 ˝C (NpCO 2 ) and chla; (b,c) are the relationships between δpCO 2 (the difference between satellite-derived pCO 2 and in situ pCO 2 ) and chla or salinity, respectively.The high value in salinity and temperature (red color in all figures) were located in the east-most leg in the cruise track of M2 in Figure 7.

The Influence of Mixing
As a first step toward implementing pCO 2 -MeSAA for the Bering Sea, we chose the summer season to eliminate the influence of vertical mixing caused by the seasonal cooling; however, mixing still occurs during the summer in some areas of the Bering Sea.The main currents in the Bering Sea include the westwards flowing ANSC along with the Aleutian archipelago and the anticlockwise flowing BSC in the slope [34,35].Seasonal upwelling, caused by wind and eddies (due to the difference of currents velocity), on the southeastern slope and southern Bering Sea close to the Aleutian archipelago, would also occur during the summer (e.g., [2,46,80,90,91]) and bring high DIC and nutrient concentrations from subsurface water to the surface layer [24].The slope in the northeast was influenced by the broad, northwestwards flowing BSC and appeared to be a system of eddies rather than a continuous current [80,82], therefore the shelf-basin exchange processes of intensive tidal mixing and eddies likely play an important role in pCO 2 variability.Mixing processes supply nutrients to the regions adjacent to the slope and enhance primary production, eventually leading to a biological drawdown of pCO 2 [22].The shelf was affected by mixing with the Yukon River and Kush Yanukovich River plumes increasing pCO 2 from organic carbon re-mineralization.
In this preliminary study we analyzed underway data in the central basin to access the influence of mixing on pCO 2 .There were only two cruises in July 2010 and August 2011 that measured pCO 2 in the central basin during summer months.Since most of the corresponding satellite chla data in July 2010 were not available due to the cloud cover, we were limited to using only the data from August 2011 for analysis.
We divided the central basin into three sub-areas, i.e., E1, E2, and E3 (Figure 12), because the underway data covered a vast spatial scale.The three different regions were defined using scatter plots of pCO 2 , temperature, and salinity.The relationship between pCO 2 and salinity showed that pCO 2 was influenced by water mixing in all three sub-areas (Figure 12b-d), in particular the sub-area E2 near the slope where temperature and salinity both showed a good positive relationship with pCO 2 ; this may be due to mixing of the westwards flowing branch of the BSC with the relatively low-temperature, low-salinity, and low-pCO 2 water from the branch of the westwards-flowing BSC [81].The difference between satellite-derived pCO 2 and the underway pCO 2 had some negative relationship with the in situ salinity in several of the sub-regions (Figure 12e-g), indicating the possibility to parameterize mixing using salinity data and adding this term in the MeSAA algorithm.However, the data from August 2011 in Figure 12 were also used to validate the satellite-derived pCO 2 (Figures 6 and 7); and the good agreement between the satellite results and the in situ data showed that horizontal mixing had a limited influence on sea surface pCO 2 in summer.
Remote Sens. 2016, 8, 558 19 of 25 2010 were not available due to the cloud cover, we were limited to using only the data from August 2011 for analysis.We divided the central basin into three sub-areas, i.e., E1, E2, and E3 (Figure 12), because the underway data covered a vast spatial scale.The three different regions were defined using scatter plots of pCO2, temperature, and salinity.The relationship between pCO2 and salinity showed that pCO2 was influenced by water mixing in all three sub-areas (Figure 12b-d), in particular the sub-area E2 near the slope where temperature and salinity both showed a good positive relationship with pCO2; this may be due to mixing of the westwards flowing branch of the BSC with the relatively lowtemperature, low-salinity, and low-pCO2 water from the branch of the westwards-flowing BSC [81].The difference between satellite-derived pCO2 and the underway pCO2 had some negative relationship with the in situ salinity in several of the sub-regions (Figure 12e-g), indicating the possibility to parameterize mixing using salinity data and adding this term in the MeSAA algorithm.However, the data from August 2011 in Figure 12 were also used to validate the satellite-derived pCO2 (Figures 6 and 7); and the good agreement between the satellite results and the in situ data showed that horizontal mixing had a limited influence on sea surface pCO2 in summer.Developing an adequate method for parameterizing the mixing effect is an ongoing process.We suggest that future work could focus on satellite salinity data, such as the Soil Moisture and Ocean Salinity (SMOS) and Aquarius/SAC-D satellites, as a proxy for mixing, or wind-driven upwelling could be parameterized by the satellite wind data.In the coastal area, the satellite-derived CDOM is a good proxy to give additional information of shelf water and the open sea (e.g., [92][93][94]).The CDOM behaves conservatively and has good negative correlations with salinity in many large river estuaries and plume systems where physical mixing prevails, thus CDOM can be a mixing indictor in river-plume dominated coastal oceans (e.g., [93,94]).We have used the satellite-derived CDOM as one of the inputs for pCO 2 algorithm in the Changjiang River-dominated East China Sea in our previous manuscript [33].The whole Bering Sea is not significantly influenced by large rivers (except in some northeastern inshoreareas), and the majority is basin.Actually, the controlling factors on pCO 2 variation are very complicated and the available data of multiple parameters are lacking in coastal areas, such as the localized relationship between CDOM and salinity, and the mixing behaviors of carbonate parameters, e.g., DIC and TA.Refining the methods for quantifying mixing processes in different sub-regions is needed as more data become available.

Conclusions
In this paper remote sensing inversion of the sea surface pCO 2 in the Bering Sea in summer was achieved with the satellite temperature and chlorophyll data using the MeSAA algorithm.The MeSAA algorithm was used to parameterize the individual effect of various factors that influence the variability of sea surface pCO 2 using a mechanistic analysis.First, a reference water mass with pCO 2 variability minimally influenced by mixing and biology was identified in the central basin, and then the thermodynamic equation and the parameterization of biological effect were added to the reference pCO 2 for applying to the whole area.We used the totally independent data (not used in the algorithm development) to validate the satellite products.Comparison results showed good agreement between satellite-derived pCO 2 and underway data, especially in the central basin.The mean relative error for the entire validation data set in eight-day composite MODIS was 9.87% (N = 1274, SD = 35.11µatm).
We had advanced the use of the pCO 2 -MeSAA algorithm, which was formerly applied in the East China Sea [33].The Bering Sea (excluding the coastal ocean) has predominately oceanic characteristics, while the East China Sea is a large river-dominated marginal sea [78,95,96], thus there was a need to modify the method to fit the present study region.In the East China Sea we used river and marine end-member values to determine the lateral mixing of inorganic carbon [33], whereas the reference water mass in the present study was used.Thereby, we were able to use carbonate system calculation to show the evolution of sea surface pCO 2 from winter to summer, therefore showing the usefulness of the reference water mass.
A significant feature of MeSAA is its mechanistic interpretation [33].Thus, we also discussed pCO 2 variability quantitatively by separating the contribution of different processes that influence pCO 2 .Moreover, the deviation between the satellite results and in situ data highlighted unconsidered processes; such deviations help to further understand the mechanism and parameterization of new components that influence pCO 2 variability.For the biological effect, we found that spring blooms had a delayed effect on pCO 2 in the Bering Sea in summer, but the degree of influence varied among regions.In general, the earlier the occurrence of the spring blooms, the weaker the influence of the previous biological effect on summer pCO 2 due to the sufficient re-equilibrium time by air-sea CO 2 exchange.The mixing effect was more difficult to discern using the present data.data [49], NOAA for providing the NCEP/NCAR Reanalysis data [60] and xCO2 data [61], Ifremer for providing climatological monthly average mixed layer depth data [62], and REMSS for providing monthly average wind speed data [63].Please contact the corresponding author (baiyan@sio.org.cn)regarding methods and data used in this paper.

Figure 1 .
Figure 1.(a) The isobaths and currents in the Bering Sea (black arrows denotes current, modified from Stabeno et al. [35]; Chen et al. [4]); Data distribution in month (b) and year (c) of underway pCO2 in the Bering Sea in 2000-2013 with the grey background of water depth.BSC: Bering Slope current; KC: Kamchatka Current; ANSC: Aleutian North Slope Current; ACC: Alaskan Coastal Current; AC: Anadyr Current.

Figure 1 .
Figure 1.(a) The isobaths and currents in the Bering Sea (black arrows denotes current, modified from Stabeno et al. [35]; Chen et al. [4]); Data distribution in month (b) and year (c) of underway pCO 2 in the Bering Sea in 2000-2013 with the grey background of water depth.BSC: Bering Slope current; KC: Kamchatka Current; ANSC: Aleutian North Slope Current; ACC: Alaskan Coastal Current; AC: Anadyr Current.

Figure 2 .
Figure 2. The relationship between pCO2, temperature, and salinity in the basin (a,b) and shelf (c,d) from July to September during 2000 to 2013.NpCO2 denoted the in situ pCO2 normalized to an average temperature of subset as 9.5 °C and 7.5 °C in basin (b) and shelf (d), respectively.

Figure 2 .
Figure 2. The relationship between pCO 2 , temperature, and salinity in the basin (a,b) and shelf (c,d) from July to September during 2000 to 2013.NpCO 2 denoted the in situ pCO 2 normalized to an average temperature of subset as 9.5 ˝C and 7.5 ˝C in basin (b) and shelf (d), respectively.

Figure 3 .
Figure 3. (a) Cruise track in July 2010 with temperature marked; (b) variation of pCO2 and salinity with the latitude; (c) relationship between temperature and salinity; (d) satellite climatology average (2002-2014) chla in July with the coastal areas shallower than 50 m depth masked.The three black lines in (a,b) divide the Bering Sea into four sub-regions according to latitude; Box A in (a) is the area of reference water mass (55.7-57°N, 171.4-174.2°E); the boxes of C1, C2, C3, and C4 in (d) are the typical area in the Kamchatka basin, the eastern basin, slope, and shelf, which are used in the discussion of Section 4.2.

Figure 3 .
Figure 3. (a) Cruise track in July 2010 with temperature marked; (b) variation of pCO 2 and salinity with the latitude; (c) relationship between temperature and salinity; (d) satellite climatology average (2002-2014) chla in July with the coastal areas shallower than 50 m depth masked.The three black lines in (a,b) divide the Bering Sea into four sub-regions according to latitude; Box A in (a) is the area of reference water mass (55.7-57 ˝N, 171.4-174.2˝E); the boxes of C1, C2, C3, and C4 in (d) are the typical area in the Kamchatka basin, the eastern basin, slope, and shelf, which are used in the discussion of Section 4.2.

Figure 4 .
Figure 4. Biological effect on pCO2.(a) The location of match-up data between eight-day composite satellite chla and the underway data in May and June in 2010; (b) the relationship between pCO2 and temperature; (c) the relationship between pCO2 normalized to the average temperature of 2.3 °C (NpCO2) and chla.

Figure 4 . 3 ˝C(
Figure 4. Biological effect on pCO 2 .(a) The location of match-up data between eight-day composite satellite chla and the underway data in May and June in 2010; (b) the relationship between pCO 2 and temperature; (c) the relationship between pCO 2 normalized to the average temperature of 2.3˝C pCO 2 varied both spatially and temporally over the Bering Sea during the years studied.Due to high SST in August (10.38˘0.29 ˝C in 2010 and 10.48 ˘0.36 ˝C in 2011) pCO 2 was generally higher in August (241.3˘13.1 µatm in 2010 and 252.8 ˘14.1 µatm in 2011) than in July and September, especially in the basin area with low chla.High chla concentration in the slopeand shelf areas, as well as some patchiness of algal blooms in the basin, resulted in under-saturated pCO 2 with respective to the atmosphere, with pCO 2 as low as ~100 µatm.

Figure 5 .
Figure 5. Monthly average satellite-derived pCO2 from July to September in 2010 and 2011.The black line represents isobaths of the 200 m depth, and the coastal areas shallower than 50 m depth are masked.

Figure 5 .
Figure 5. Monthly average satellite-derived pCO 2 from July to September in 2010 and 2011.The black line represents isobaths of the 200 m depth, and the coastal areas shallower than 50 m depth are masked.

Figure 6 .
Figure 6.The comparison between eight-day composite satellite-derived pCO2 and underway data.(a) Location of the match-up data; (b) scattering plot between satellite-derived pCO2 and underway data.Box D is highlighted for the derivation in the slope.

Figure 7 .
Figure 7.The comparison between monthly satellite-derived pCO2 with underway data.(a) The location of the match-up data in September 2010 (M1), July 2010 (M2), and August 2011 (M3) with the salinity marked; the comparison between satellite-derived pCO2 and underway pCO2 data in the western coastal area (M1) (b); the shelf (M2) (c); and the basin (M3) (d).

Figure 6 .
Figure 6.The comparison between eight-day composite satellite-derived pCO 2 and underway data.(a) Location of the match-up data; (b) scattering plot between satellite-derived pCO 2 and underway data.Box D is highlighted for the derivation in the slope.

Figure 6 .
Figure 6.The comparison between eight-day composite satellite-derived pCO2 and underway data.(a) Location of the match-up data; (b) scattering plot between satellite-derived pCO2 and underway data.Box D is highlighted for the derivation in the slope.

Figure 7 .
Figure 7.The comparison between monthly satellite-derived pCO2 with underway data.(a) The location of the match-up data in September 2010 (M1), July 2010 (M2), and August 2011 (M3) with the salinity marked; the comparison between satellite-derived pCO2 and underway pCO2 data in the western coastal area (M1) (b); the shelf (M2) (c); and the basin (M3) (d).

Figure 7 .
Figure 7.The comparison between monthly satellite-derived pCO 2 with underway data.(a) The location of the match-up data in September 2010 (M1), July 2010 (M2), and August 2011 (M3) with the salinity marked; the comparison between satellite-derived pCO 2 and underway pCO 2 data in the western coastal area (M1) (b); the shelf (M2) (c); and the basin (M3) (d).

Figure 9 .
Figure 9. (a) The stations from an east-to-west section at 57°N in the Bering Sea basin in August 2004; (b) Ttemperature profiles of all stations.

Figure 11 .
Figure11.The relationship between pCO2, chla, temperature, and salinity in the shelf in July 2010 (locations were marked as M2 in Figure7).(a) The relationship between pCO2 normalized to reference temperature 7.7 °C (NpCO2) and chla; (b,c) are the relationships between δpCO2 (the difference between satellite-derived pCO2 and in situ pCO2) and chla or salinity, respectively.The high value in salinity and temperature (red color in all figures) were located in the east-most leg in the cruise track of M2 in Figure7.

Figure 12 .
Figure 12.The relationship between pCO2, temperature, and salinity in the basin in August 2011.(a) Cruise tracks with data divided into three groups marked by E1, E2, and E3 (b-d) are the relationship between pCO2 normalized to the average temperature 10.8 °C (NpCO2) and salinity in E1, E2, and E3, respectively; (e-g) are the relationships between δpCO2 (the difference between satellite-derived pCO2 and in situ pCO2) and the underway salinity data.

Figure 12 .
Figure 12.The relationship between pCO 2 , temperature, and salinity in the basin in August 2011.(a) Cruise tracks with data divided into three groups marked by E1, E2, and E3 (b-d) are the relationship between pCO 2 normalized to the average temperature 10.8 ˝C (NpCO 2 ) and salinity in E1, E2, and E3, respectively; (e-g) are the relationships between δpCO 2 (the difference between satellite-derived pCO 2 and in situ pCO 2 ) and the underway salinity data.

Table 1 .
Basic information for underway data and discrete carbonate data from CDIAC used in this research.

Table 2 .
Values for the different parameters in the reference water mass and water in the minimum temperature layer in the summer.Data in the parentheses are the standard deviation from the average values of N = 238 data points.The values of the reference water in the basin from July 2010 (Box A in Figure 3).b The value of the temperature minimum layer at a depth of 100 m from August 2004.c The values of well-mixed winter water, calculated from the summer data in the temperature minimum layer (Station 38); see the text for detail. a