A feed-forward neural network is used to create a monthly climatology of the sea surface fugacity of CO2 (fCO2) on a 1° × 1° spatial resolution. Using 127 880 data points from 1990 to 2011 in the track-gridded database of the Surface Ocean CO2 Atlas version 2.0 (Bakker et al.), the model yields a global mean fCO2 increase rate of 1.50 μatm yr−1. The rate was used to normalize multiple years’ fCO2 observations to the reference year of 2000. A total of 73 265 data points from the normalized data were used to model the global fCO2 climatology. The model simulates monthly fCO2 distributions that agree well with observations and yields an anthropogenic CO2 update of −1.9 to −2.3 PgC yr−1. The range reflects the uncertainty related to using different wind products for the flux calculation. This estimate is in good agreement with the recently derived best estimate by Wanninkhof et al. The model product benefits from a finer spatial resolution compared to the product of Lamont–Doherty Earth Observatory (Takahashi et al.), which is currently the most frequently used product. It therefore has the potential to improve estimates of the global ocean CO2 uptake. The method’s benefits include but are not limited to the following: (i) a fixed structure is not required to model fCO2 as a nonlinear function of biogeochemical variables, (ii) only one neural network configuration is sufficient to model global fCO2 in all seasons, and (iii) the model can be extended to produce global fCO2 maps at a higher resolution in time and space as long as the required data for input variables are available.
It has been estimated that the global ocean absorbs nearly half the total emissions of anthropogenic carbon dioxide (CO2) from the atmosphere (Sabine et al. 2004; Jacobson et al. 2007; Gruber et al. 2009; Takahashi et al. 2009). The surface ocean shows much more CO2 variability spatially than the atmosphere, as is shown in Komhyr et al. (1985) and Takahashi et al. (2009); hence, it is import to accurately quantify the oceanic CO2 distribution in order to improve our understanding of the ocean carbon sink. Meanwhile, in spite of many years’ efforts by researchers to measure the sea surface CO2 partial pressure (e.g., Feely et al. 1987; Inoue et al. 1995; Cooper et al. 1998; Bates et al. 1996; Wong et al. 1999, 2010; Murphy et al. 2001; González-Davila 2003; Keeling et al. 2004; Midorikawa et al. 2006; Luger et al. 2006; Schuster and Watson 2007; Borges et al. 2008), in situ measurements are still insufficient in most parts of the oceans because of technical and financial restrictions. This circumstance has motivated the search for robust and reliable methods to interpolate available CO2 data spatially and temporarily to obtain basinwide (e.g., Zeng et al. 2002; Lefèvre et al. 2005; Chierici et al. 2006; Sarma et al. 2006; Jamet et al. 2007; Friedrich and Oschlies 2009; Telszewski et al. 2009; Takamura et al. 2010; Landschützer et al. 2013; Nakaoka et al. 2013) and global ocean CO2 maps (Takahashi et al. 2002, 2009). Conventional interpolation methods model CO2 as a function of a number of biogeochemical variables. Recently, self-organizing maps (SOM), a type of artificial neural network, gained popularity in mapping the surface ocean CO2 (Lefèvre et al. 2005; Friedrich and Oschlies 2009; Telszewski et al. 2009; Nakaoka et al. 2013; Landschützer et al. 2013). Using an artificial neural network differs distinctively from conventional statistical methods in that it is not required to formulate the model structure a priori.
In this work we introduce a new concept for mapping the surface ocean fCO2 (fugacity of CO2) using a feed-forward neural network (FNN). As an artificial neural network, an FNN models and interpolates data nonlinearly, whereas a SOM establishes discrete estimates using a category map of limited size that associates a target variable with dependent variables (e.g., Telszewski et al. 2009). Previous studies applied a SOM to confined physical areas where biogeochemical variables were assumed to be able to interpret the fCO2 distribution well in a season. The concept is to exclude the temporal and spatial variables from the dependent variables of fCO2 so that scattered, multiple years’ observations can be pooled together to train a neural network. The new concept of our model includes the dependence of fCO2 on time and space in the first place and considers biogeochemical variables as proxy variables for complementing insufficient observations. We gain two benefits from the new concept: (i) the model simulates observations faithfully where and when observations are sufficient and (ii) the monthly global climatology can be modeled by just one neural network.
a. Model equations
Our analysis is based on a well-understood relationship stating that the fCO2 is a nonlinear function of time and space, which in our model are represented by month (MON), latitude (LAT), and longitude (LON), that is,
Available observations are scarce with respect to the biogeochemical properties of the surface ocean; hence, a direct interpolation would include abnormally large uncertainties in many areas. One way to improve a direct interpolation is to include proxy variables that connect fCO2 variability with time and space to extend the model to unobserved areas. Following the works of Lefèvre et al. (2005), Friedrich and Oschlies (2009), Telszewski et al. (2009), Landschützer et al. (2013), and Nakaoka et al. (2013), we choose sea surface temperature (SST), sea surface salinity (SSS), and chlorophyll-a concentration (CHL) as proxy variables. Therefore, our basic model has the form of
We exclude the mixed layer depth (MLD) used by Telszewski et al. (2009), Nakaoka et al. (2013), and Landschützer et al. (2013), because openly available MLD data are scarce and our preliminary analysis using World Ocean Atlas 1994 (WOA94) data (http://www.nodc.noaa.gov/OC5/WOA94/mix.html; Monterey and Levitus 1997) indicates that MLD (not shown here) has the least effect on fCO2 compared to SST, SSS, and CHL.
The MON and LON variables cannot be used directly because of their periodical properties. Therefore, we replaced each with two of its transformations, resulting in an updated model equation of the form
b. Neural network
Figure 1 illustrates the neural network used to model Eq. (3). It is a typical FNN having a layered structure and including three layers known as input, hidden, and output layers (Svozil et al. 1997). In our study, the input layer consists of eight neurons with one neuron corresponding to one independent variable in Eq. (3); the output layer has only one neuron whose output gives the scaled fCO2; and the number of neurons in the hidden layer is empirically determined, which will be discussed later.
A neuron in the input layer takes the value of an independent variable and passes it to all neurons in the hidden layer without any transformation. A neuron in the hidden layer adds a bias to the weighted sum of all inputs (outputs from its upstream neurons in the input layer) and uses the transformed sum as an input for the output neuron, that is,
in which M is the total number of input neurons. The neuron in the output layer transforms its inputs in the same way.
The sigmoid function is chosen because it increases monotonically from 0 to 1 for x from −∞ to +∞ and the derivative of y with respect to x has a simple form that affiliates the so-called error backpropagation algorithm (Rumelhart et al. 1986) for neural network training. The training algorithm adjusts the bias and weighting factors according to the negative gradient of the error cost function:
in which P is the total number of training patterns (a pattern comprises a set of observations for fCO2 and the proxy variables), y is the output of the FNN, and d is the target variable (i.e., the scaled fCO2). We used the Levenberg–Marquardt method (Wilamowski and Yu 2010) to minimize the error cost function.
c. Neural network training
It is well known that, given a sufficiently large number of hidden neurons, an FNN can approximate any finite function (Hornik 1991; Blum and Li 1991). That indicates that if the fCO2 values are unique in all training patterns and correspond to unique values of MON, LAT, LON, SST, SSS, and CHL, Eq. (3) can be fitted accurately with a larger number of hidden neurons. An FNN configured as such remembers all training patterns and thus loses the capability of filtering out erroneous information in the source data. Unfortunately, there is no theoretically approved method for choosing the right number of hidden neurons (Svozil et al. 1997). Our criterion for selecting the number is that the standard deviation of the difference between model outputs and observations (std1) is no larger than the standard deviations of repeated fCO2 observations averaged over all months and grid boxes (std2) by 20%, that is, |std1 − std2|/std2 < 0.2. At the experimental stage, we tried 4, 8, 16, 32, 64, and 128 hidden neurons. Eventually we settled on 32 according to the criterion.
d. Flux calculation
We further calculate the sea–air CO2 flux as a function of the partial pressure of CO2 (pCO2) using
where kp is the gas transfer coefficient, kw is the gas transfer velocity, and k0 is the solubility of CO2 in seawater. We use the formulation of Wanninkhof (1992) for the coefficient kw (cm h−1) with the scaling factor of Takahashi et al. (2009); that is,
where u (m s−1) is the wind speed 10 m above the sea level and Sc is the Schmidt number, which varies with SST (°C; see Wanninkhof 1992):
The solubility k0 (mol L−1 atm−1) is a function of temperature T (K) = SST + 273.15 and SSS (Weiss 1974):
We multiply kw by 7.03 to convert its unit to (m month−1) and k0 by 0.001 to convert its unit to (mol m−3 μatm−1), so that Eq. (6) gives flux in (mol m−2 month−1).
The conversion between fCO2 and pCO2 is done by the following formulation (see Weiss 1974):
where R = 82.057 46 (cm3 atm mol−1 K−1), Ps is the sea level air pressure used to substitute the equalibrator pressure, and B and δ are the virial coefficients for CO2. Coefficient B (cm3 mol−1) is given by
and δ (cm3 mol−1) by
a. CO2 data
We extracted monthly fCO2 data from the track-gridded database of the Surface Ocean CO2 Atlas (SOCAT) version 2.0 (http://www.socat.info/; Pfeil et al. 2013; Sabine et al. 2013; Bakker et al. 2013). The database has a 1° × 1° spatial resolution and includes global measurements from 1970 to 2011. We excluded some data points using four criteria: (i) year before 1990, (ii) fCO2 values smaller than 200 μatm or larger than 600 μatm, (iii) ocean depth smaller than 1000 m, and (iv) salinity smaller than 30.0. The first criterion is for reducing the uncertainty of normalizing fCO2 to the reference year of 2000 for comparison with the global CO2 maps of Takahashi et al. (2009). The second criterion is for reducing the potential effect of extreme data points on the extrapolation of the FNN. The last two criteria are for identifying seawater in open oceans.
We obtained 127 880 data points or patterns from multiple years’ data for training the FNN to estimate the global mean increase rate of fCO2. The rate was then used to normalize multiple years’ fCO2 to obtain the monthly fCO2 climatology, which comprises the means of repeated observations in the same months and same grid boxes. From the climatology, we obtained 73 265 data points for training the FNN to model the global fCO2 climatology. Figure 2 shows the composite sampling map of SOCAT in 1990–2011.
We prepared more CO2 data for training the FNN to evaluate the difference between the results of the FNN and Takahashi et al. (2009). In addition to the four criteria, we excluded data points in the equatorial Pacific, 10°N–10°S, in the El Niño periods, which were selected according to the National Oceanic and Atmospheric Administration (NOAA)’s analysis (http://www.cpc.ncep.noaa.gov/products/analysis_monitoring/ensostuff/ensoyears.shtml), resulting in 71 287 data points. Unless specifically stated in latter discussions, our FNN model results refer to those obtained by using the first CO2 dataset.
b. Proxy variable data
The monthly SST data of 1990–2011 were extracted from the NOAA optimum interpolation (OI) version 2 product (http://www.esrl.noaa.gov/psd/data/gridded/data.noaa.oisst.v2.html; Reynolds et al. 2002). We calculated the SST climatology using SST in that period. The SSS climatology was extracted from the World Ocean Atlas 2009 (WOA09) product (http://www.nodc.noaa.gov/OC5/WOA09/netcdf_data.html; Antonov et al. 2010). The WOA09 SSS covers the period of June 1890–December 2008. The CHL climatology was calculated using the Moderate Resolution Imaging Spectroradiometer (MODIS) Aqua and Sea-viewing Wide Field-of-view Sensor (SeaWiFS) climatology (http://oceancolor.gsfc.nasa.gov/cgi/l3). The 9-km gridded satellite data were binned into 1° × 1° grid boxes. The mean of the two CHLs was used as the CHL climatology for our model. At the time we obtained the data, Aqua CHL covers the period of July 2002–July 2013 and SeaWiFS CHL from September 1997 to December 2010.
c. Data for flux calculation
The required monthly wind speed and air pressure were calculated using the National Centers for Environmental Prediction Final Operational Global Analysis (NCEP-FNL) data (http://rda.ucar.edu/datasets/ds083.2/) of the years 2000–11, as the data before July 1999 are unavailable. The product has a spatial resolution of 1° × 1° and a temporal resolution of 6 h. We used the atmospheric pCO2 in the database (http://www.ldeo.columbia.edu/res/pi/CO2/) of Takahashi et al. (2009) in order to compare with their flux results. For estimating the global flux, we filled our unmodeled grid boxes, in which CHL is unavailable, using their seawater pCO2 values. The unmodeled area size is about 6% of our modeled area size.
d. Data scaling
Equation (4) indicates that data used as input can take any real values. However, an FNN usually initializes bias and weight parameters randomly in a predetermined range for generalizing its use with all kinds of inputs; therefore, scaling input data speeds up the convergence within the training process. We scale the values υ of an input variable by their mean and standard deviation σ as
so that they have a zero mean and a standard deviation of one. The observed fCO2 is scaled using
This confines the scaled fCO2 to between 0.1 and 0.9. The reason for the confinement is that Eq. (4) results in y values within the range of 0–1 and that a small change in y near zero or one requires a very large change in x, which will slow the convergence of training.
Before being scaled by Eq. (13), chlorophyll is transform by
The transformation helps to reduce the strongly skewed CHL values; that is, a majority number of CHL values are very small. The skewed distribution potentially leads to a trained FNN to have neuron parameters more optimized toward fitting small CHL well, thus weakening its prediction power for large CHL.
4. Results and discussions
a. CO2 increase rate
As the surface ocean CO2 concentration varies in responding to changes in biogeochemical conditions and to the increase of atmospheric CO2 concentration, the increase rate of the surface ocean CO2 showed geographical differences (see the summary in Takahashi et al. 2009). A constant CO2 increase rate was usually used to normalize multiple years’ CO2 data to a reference year (e.g., Takahashi et al. 2009; Nakaoka et al. 2013). Instead of adopting an available rate for our fCO2 data normalization, we used the FNN to estimate the global mean rate by repeating the following procedures until the error in Eq. (5) became stable: (i) modeling Eq. (3) with the fCO2 and SST of 1990–2011 and the climatology of SSS and CHL; (ii) linearly fitting the difference between the modeled and observed fCO2 against time (year) to estimate the annual increase rate; and (iii) adjusting the observed fCO2 by the rate and then going to procedure (i) with the adjusted fCO2 as the target.
The procedures are equivalent to treating the fCO2 increase rate as an error term in the FNN modeling and treating the seasonal and spatial variations of fCO2 as an error term in linear fitting. By recursively trying to guess the variations, removing the variations to guess the trend, and detrending the observed fCO2 to guess the variations again, the procedures effectively reduced the seasonal and spatial variations of CO2 in the residue. As a result, the trend signal became stronger with the increasing number of iterations. At first, a direct linear fit of the observed fCO2 resulted in a rate of 0.855 μatm yr−1 and a correlation coefficient of 0.144. After the first recursion, the FNN improved the correlation to 0.332 and gave a global mean rate of 1.085 μatm yr−1. The magnitude of improvement decreased with the number of iterations. The mean FNN error became stable after repeating procedures (i) to (iii) 5 times. Finally, the FNN gave a rate of 1.502 μatm yr−1 and a correlation coefficient of 0.443 (Fig. 3). The stabled increase rate is consistent with previously used values (e.g., Takahashi et al. 2009).
b. Nonlinear fitting
The annual increase rate of 1.502 μatm yr−1 was used to normalize monthly fCO2 of 1990–2011 to the reference year 2000. The FNN was further trained by using the normalized fCO2 and the climatology of SST, SSS, and CHL. To evaluate the predictive performance, we randomly took 90% of the 73 265 patterns for training use and the remaining 10% for validation. By using a different random seed, we could test cases with different subsamples. The results are all similar to those shown in Fig. 4. The linear regression of modeled fCO2 versus observed fCO2 in the training samples gives a slope of 0.784 and a correlation coefficient of R = 0.886. The mean and standard deviation of the bias (FNN – SOCAT) are −0.01 and 15.8 μatm, respectively. Considering the mean standard deviation of repeated observations of 13.2 μatm, these values indicate that the bias is negligibly small and that the FNN models the nonlinear relationship of Eq. (3) well. For the validation subsamples, the slope and correlation coefficient are 0.767 and 0.871, respectively; and the mean and standard deviation of the bias are −0.16 and 16.0 μatm, respectively. The validation results further indicate that the FNN model is acceptable for interpolating fCO2 to unsampled areas. Finally, we used all samples to train the FNN in order to produce global fCO2 maps. It yields a slope of 0.778, a correlation coefficient of 0.882, and a mean and standard deviation of the bias of −0.01 and 15.4 μatm, respectively.
Including space and time variables in the model yields a good representation of the seasonal cycle of fCO2. Figure 5 shows four examples for the western Pacific near the meridian line of 145°E, where the volunteer cargo ship observation of the National Institute for Environmental Studies supplies every monthly observation data (Nakaoka et al. 2013). They reveal three seasonal patterns shared by other oceans: (i) temperature-driven seasonality of low fCO2 in winter and high fCO2 in summer (Fig. 5a), (ii) fCO2 drawdown by biological activities from spring to summer (Fig. 5b), and (iii) second fCO2 drawdown by biological activities in autumn (Fig. 5c). In many areas in the tropical zone from about 20°S to 20°N, persistently high temperature, and low biological activity lead to less pronounced seasonal changes with high fluctuations (Fig. 5d). The coupling of CO2 variability with temperature and biological activities is in good agreement with previous findings (e.g., Bates 2001; Takahashi et al. 2002; Zeng et al. 2002).
For the convenience of discussion, we use FNN1 and FNN2 to represent the FNN trained with and without data points in the El Niño periods (see section 3b), respectively. Similarly, we use OBS1 for SOCAT observations normalized to the year 2000 and OBS2 for the same data but without data points in the El Niño periods. And we use the name T09 for the Lamont–Doherty Earth Observatory (LDEO) climatology of Takahashi et al. (2009), who interpolated the monthly pCO2 to obtain global distributions based on a lateral, two-dimensional advection–diffusion transport for surface ocean. The distributions of the annual mean of fCO2 difference among models and observations are shown in Fig. 6. Monthly statistics of the fCO2 differences are listed in Table 1.
Our model results show small monthly-mean biases (within 1.5 μatm) compared to the observations (Table 1). The overall bias of both FNN1 and FNN2 is zero. The T09 results also show small monthly-mean biases (within 1.1 μatm), but a small negative overall bias (0.4 μatm) and a larger overall uncertainty (17.5 μatm) compared with FNN1 and FNN2 (15.4 μatm). Between FNN2 and T09, small monthly-mean differences (within 0.5 μatm) occur in April, November, and December (Table 1). Relatively large monthly-mean differences (up to 4.2 μatm) occur in other months. Annually, the mean difference is about 0.9 ± 15.8 μatm.
Large differences between FNN2 and T09 occur near the equator and in areas where observations are relatively scarce (Fig. 6a). Compared to the observations, large differences seem to be associated with undersampling and large CHL uncertainty for the FNN and the T09 models (Figs. 6b,c). As shown in Fig. 5d, the observed fCO2 in the tropical area includes large variability throughout the year, while the amplitude of the seasonal variation is relatively small. As a result, both the FNN2 and T09 models show relatively large differences against OBS in the equatorial areas (Figs. 6b,c).
We put both modeled and observed fCO2 for September, when the T09−FNN2 difference is the largest, in Fig. 7. The T09−FNN2 difference fluctuates at the basin scale in most areas (Fig. 7a). T09 shows particularly high fCO2 unlike FNN2 in the zone near 30°S and near the Antarctic. By comparing both models’ results (Figs. 7b,c) with the observations (Fig. 7d), it becomes clear that a large T09−FNN2 difference is likely caused by T09’s overestimation and FNN2’s underestimation. In contrast, the T09 fCO2 is generally smaller than the FNN2 fCO2 in the equatorial areas. We attribute the discrepancy to the model uncertainty caused by relatively large fCO2 fluctuation in the same month.
Nevertheless, both the FNN and T09 models yield results that agree generally well with observations throughout the year (Fig. 8). Our monthly climatology maps, which include El Niño data points, illustrate clear seasonal fCO2 variations. The FNN product has the benefit of a finer spatial resolution compared to the T09 product, but it has the disadvantage of larger areas of missing data because of missing CHL inputs.
We evaluated the uncertainty of the FNN caused by using different proxy variable data. As the evaluation on a second dataset did not meaningfully change the conclusions, the results are not shown here.
d. CO2 flux
The gas transfer coefficient kp in Eq. (6) is primarily a function of wind speed (Takahashi et al. 2009); hence, the flux calculation can differ significantly depending on the wind product used. By applying the wind speed from a different wind product to the pCO2 of T09, Wanninkhof et al. (2013) yielded a global sea–air CO2 flux of −1.18 PgC yr−1, which is 15% less negative than that of T09 (−1.38 PgC yr−1). Using the wind and pCO2 data of T09, our formulation also yields a global air–sea CO2 flux of −1.39 PgC yr−1 even though different SST and SSS data were used (Table 2). However, when the T09 wind was replaced by the NCEP–FNL wind (see section 3c), we yielded a global flux of −1.04 PgC yr−1 for the T09 pCO2, which is about 25% less negative. Being aware of this uncertainty, we present our flux estimate as the range resulting from both wind data. For comparison reasons, we filled 6% of the areas, where there were no direct estimates from the FNN, by T09 CO2 data to estimate the CO2 flux in this study. Fluxes thus calculated are summarized in Table 2.
Overall, our results suggest a global sea–air CO2 flux ranging from −1.31 to −1.66 PgC yr−1 by including the El Niño periods (FNN1) and from −1.22 to −1.57 PgC yr−1 by excluding those periods (FNN2). They are about 15%–25% more negative than the estimate in T09 when the same wind product is used. Averaging zonally, the Southern Ocean south of 20°S contributes the most (about −1.17 PgC yr−1) to the total flux, while the northern oceans north of 20°N contribute roughly two-thirds the amount (about −0.75 PgC yr−1). In contrast, the equatorial zone from 20°S to 20°N appears to be a prominent source with a net flux of 0.61 PgC yr−1. The relative zonal contributions agree well with the numbers obtained by Takahashi et al. (2009).
As the FNN uncertainty caused by data points in the El Niño periods is small, we report our model flux based on the FNN1 results. By adding a correction of −0.2 PgC yr−1 for undersampling (Wanninkhof et al. 2013) to the direct estimate, we obtain a net annual CO2 uptake of −1.5 to −1.9 PgC yr−1 for the global ocean. Furthermore, by including an estimate for the riverine CO2 flux of 0.4 PgC yr−1 to the oceans (e.g., Jacobson et al. 2007), we yield an anthropogenic CO2 uptake of −1.9 to −2.3 PgC yr−1. These values are in good agreement with the best estimate for the anthropogenic CO2 uptake by Wanninkhof et al. (2013).
Our global fCO2 climatology derived from a feed-forward neural network, including the time and space variables, reveals a great potential to produce an fCO2 product for the modeling community. The model generally simulates the seasonal variation of fCO2 well and yields monthly distributions that are in good agreement with observations. Our product also agrees well with the LDEO product (Takahashi et al. 2009), but it has a much higher spatial resolution. Our method has the advantage of not having to define a fixed model relationship between fCO2 and its dependent variables. Its advantage over the SOM method (e.g., Telszewski et al. 2009) is that the interpolation of fCO2 is continuous, whereas an SOM establishes discrete estimates using a category map. Regarding the method of Landschützer et al. (2013) that uses both SOM and FNN, our method is much simpler in that only one neural network configuration is sufficient to model the global fCO2 in all seasons, whereas Landschützer et al. (2013) focused their method on the interannual variability of the Atlantic Ocean carbon sink. We achieved this by including transformed time and space as input parameters. The model is robust and can be extended to produce global fCO2 with a higher resolution in time and space as long as required data for input variables are available.
Our anthropogenic CO2 uptake estimate agrees well with the best estimate of Wanninkhof et al. (2013), but it appears to be about 15%–25% larger than the estimate of Takahashi et al. (2009). The uncertainty caused by including or excluding data in the El Niño periods is small compared to the difference between models. The difference caused by using different wind products can be as large as the difference between models and provides the largest CO2 flux uncertainty of our method.
We thank all the scientists and institutes who contributed to the SOCAT database.