Articles | Volume 12, issue 3
Data description paper
06 Jul 2020
Data description paper |  | 06 Jul 2020

A 17-year dataset of surface water fugacity of CO2 along with calculated pH, aragonite saturation state and air–sea CO2 fluxes in the northern Caribbean Sea

Rik Wanninkhof, Denis Pierrot, Kevin Sullivan, Leticia Barbero, and Joaquin Triñanes

A high-quality dataset of surface water fugacity of CO2 (fCO2w)1, consisting of over a million observations, and derived products are presented for the northern Caribbean Sea, covering the time span from 2002 through 2018. Prior to installation of automated pCO2 systems on cruise ships of Royal Caribbean International and subsidiaries, very limited surface water carbon data were available in this region. With this observational program, the northern Caribbean Sea has now become one of the best-sampled regions for pCO2 of the world ocean. The dataset and derived quantities are binned and averaged on a 1 monthly grid and are available at (last access: 30 June 2020) (; Wanninkhof et al., 2019a). The derived quantities include total alkalinity (TA), acidity (pH), aragonite saturation state (ΩAr) and air–sea CO2 flux and cover the region from 15 to 28 N and 88 to 62 W. The gridded data and products are used for determination of status and trends of ocean acidification, for quantifying air–sea CO2 fluxes and for ground-truthing models. Methodologies to derive the TA, pH and ΩAr and to calculate the fluxes from fCO2w temperature and salinity are described.

1 Introduction

Over the past 20 years a rapidly expanding program of measurements of surface water partial pressure of carbon dioxide (pCO2w)1 or fugacity of CO2 (fCO2w)1 has provided data to determine air–sea CO2 fluxes and rates of ocean acidification on local to global scales (e.g., Boutin et al., 2008; Degrandpre et al., 2002; Evans et al., 2015; Schuster et al., 2013; Takahashi et al., 2014; Wanninkhof et al., 2019a, b). Marginal seas, which have historically had a dearth of measurements, have been targeted for increased observations. Through an industry, academic and federal partnership between the U.S. National Oceanic and Atmospheric Administration (NOAA), Royal Caribbean International (RCI) and the University of Miami, cruise ships were outfitted with automated surface water pCO2 systems, also called underway pCO2 systems (Pierrot et al., 2009). In 2002, the RCI ship Explorer of the Seas (EoS) was equipped with an underway pCO2 setup providing observations on alternating weekly transects from Miami, FL, to the northeastern and northwestern Caribbean. In 2015, the EoS was repositioned out of the Atlantic, and the Celebrity Equinox (Eqnx) was outfitted with an underway pCO2 system. Additionally, an underway pCO2 system was installed on the Allure of the Seas (ALoS) in 2016. The Eqnx and ALoS covered similar transects as the EoS but on a more irregular and seasonal basis. A total of 582 cruises covered the region from 2002 through 2018. A map of all cruise tracks is shown in Fig. 1. The number of cruises per year covering the Caribbean Sea and adjacent western Atlantic is provided in Fig. 2. There are fewer cruises in the middle part of the record, when the EoS was diverted to other routes outside the area and eventually repositioned.

Figure 1Map with the cruise tracks of the EoS, Eqnx and ALoS from 2002 though 2018. The gray scales show bathymetry, with white being less than 200 m. The green rectangle depicts the region where the data and products are compared in Fig. 4.

Figure 2Histogram of number of cruises per year used in this work.


The surface water fCO2 observational dataset and derived products including total alkalinity (TA), acidity (pH), aragonite saturation state (ΩAr) and air–sea CO2 flux are of importance for determining the anthropogenic carbon uptake and to assess trends and impacts of ocean acidification. The observational data are provided to the global Surface Ocean CO2 Atlas (SOCAT; Bakker et al., 2016) and global CO2 climatology (Takahashi et al., 2009, 2018). These data are the main source of fCO2 observations available in the region, and the high frequency of measurements provides a seasonally resolved picture of changing fCO2w. This effort has made the northern Caribbean one of the few places in the world ocean where such regional observational density has been established. The data and mapped products are interpreted in Wanninkhof et al. (2019b), who show large decadal changes in trends of surface water fCO2 and associated changes in air–sea CO2 fluxes.

The data from the first part of this record form the basis of the ocean acidification product suite that maps ocean acidification conditions in the Caribbean (Gledhill et al., 2008;, last access: 12 April 2020;, last access: 12 April 2020). The large, high-quality and well-resolved dataset is also used to validate models (Gomez et al., 2020).

For optimal application, datasets and associated data products are fully documented here and are readily accessible according to findable, accessible, interoperable and reusable (FAIR) principles (Wilkinson et al., 2016). In addition, the data underwent substantial quality control by our group and through the SOCAT quality control and check procedures. The documentation of the cruises, the sampling methodology and data reduction techniques are presented in brief. This is followed by a description of the approaches to calculate the different inorganic carbon system parameters. The procedure is to bin and average the fCO2w on a 1 by 1 grid at monthly scales, referred to as gridded observations. Then the so-called second inorganic carbon system parameter is calculated. We estimate the total alkalinity (TA) based on robust relationships of TA with salinity (Cai et al., 2010; Takahashi et al., 2014; Lee et al., 2006; Millero et al., 1998) and use the software program CO2SYS (Pierrot et al., 2006) to calculate the other inorganic carbon system parameters of interest, in this case pH and ΩAr. These data products are presented at monthly scales binned and averaged on a 1 by 1 grid, referred to as gridded products. Annual multilinear regressions (MLRs) are developed between the gridded fCO2w data and sea surface temperature (SST), sea surface salinity (SSS), location (latitude: lat; longitude: long) and mixed-layer depth (MLD) as independent variables. These regressions are applied at monthly and 1 by 1 spatial resolution to the region between 15 and 28 N and between 88 and 62 W using remotely sensed or modeled independent variables, and they are used to calculate air–sea CO2 fluxes. These calculated parameters are called gridded mapped products. The procedures and datasets, including uncertainty analyses and tables of column headers of the product created, are provided below.

2 Observations

The observational program is described in terms of the ships, voyages and instrumentation. The ships predominantly sailed in the Caribbean Sea but also had tracks outside the region, including in the northeast of the USA as well as in Bermuda and in the Mediterranean. The description of operations, data and products presented here covers the northern Caribbean Sea and western Atlantic, north of the Caribbean island chain covering the region from 15 to 28 N and −62 to −88 (=62 to 88 W; Fig. 1). The home ports of the ships, where passengers embark and disembark, are Miami and Fort Lauderdale, FL. The Explorer of the Seas (EoS) changed home port from Miami, FL, to Cape Liberty Cruise Port, NJ, in 2008 and changed its routes at that point to include cruises with Bermuda as a port of call. In 2015 the EoS was repositioned to the Pacific, and the underway pCO2 system was removed. From 2015 onward the Celebrity Equinox (Eqnx) and from 2016 onward the Allure of the Seas (ALoS) covered the area. The Eqnx spent the summers of 2015 and 2016 in the Mediterranean, causing seasonal data gaps in the Caribbean.

2.1 Cruises

The EoS had 331 cruises from 2002 to 2015, the Eqnx completed 135 cruises from 2015 through 2018, and the ALoS performed 116 cruises in the study area from 2016 through 2018. Temporal coverage over the 17 years shows occupations every two weeks at the beginning of the record, from 2002 to 2007, and at the end, from 2014 through 2018, with fewer occupations in the years in between (Fig. 2). The cruises lasted between 7 and 14 d and made about a half a dozen ports of call. The ships were generally in port from early morning to late afternoon and transited between ports at night except for long runs (e.g., from Miami to San Juan), when the ship sailed continuously for several days. Ports are listed in the metadata accompanying the original data.

The systems were installed in different locations for the three ships, but each had a dedicated seawater intake near the bow. The EoS had an intake in the bow thruster tube (≈3 m depth) that was nonoptimal due to bubble entrainment during bow thruster operations and heavy seas. These observations have been culled from the datasets. The ALoS and Eqnx had their intake at ≈5 m depth, but it was forward from the bow thruster and had fewer issues with bubble entrainment. On the EoS, the underway pCO2 instrument was initially in a dedicated science laboratory built for use on the ship and located amidships about 100 m from the intake. In 2008, a new system was placed in the engineering space closer to the bow with no apparent change in performance. For the ALoS and Eqnx, the instruments were near the bow intake in the engineering space, about 5 m from the intake. The EoS had an air intake mounted on a mast at the forward-most point of the main deck in August 2008. The Eqnx had an air intake near the bow, one level below the main deck, since the initial installation in March 2015. The underway pCO2 systems on these ships made marine boundary layer (MBL) air observations (xCO2a) as described in Wanninkhof et al. (2019b), but these are not used in the data products presented. Typical cruise speeds were 22 kn, which, with sampling every 2.5 min, yielded a fCO2w sample approximately every 1.7 km except for 35 min every 4.5 h, when four calibration gases and a CO2-free reference gas were analyzed, followed by five atmospheric CO2 measurements for the ships with air intakes. This created a gap of 24 km (1/4) without fCO2w measurements.

2.2 Instrumentation

2.2.1pCO2 system

The instrumentation is based on a community design described in Pierrot et al. (2009). The instruments were manufactured by General Oceanics Inc., and for the cruises described here they have performed to high-accuracy specifications. The surface water drawn from the intake at about 4 L min−1 went through a 1.1 L spray head equilibrator with a water volume of 0.5 L and a headspace of 0.6 L. The spray and agitation caused the CO2 in the headspace to equilibrate with the CO2 in the water with a response time of about 2 min (Pierrot et al., 2009; Webb et al., 2016). Thus, the air in the headspace reached 99.8 % equilibration in 12 min. As the air is recirculated, and surface waters were relatively homogeneous on hourly timescales, 100.0 % equilibration was assumed. The four calibration gases supplied by the Global Monitoring Division (GMD) of the Environmental Science Research Laboratory (ESRL) of the NOAA were traceable to the WMO CO2 mole fraction scale. The CO2 concentrations of the standards spanned the range of surface water values encountered along the ship tracks (≈280–480 ppm).

The gas entering the analyzer was dried by passing it through a thermoelectric cooler at 5 C and a PermaPure drier. The standards did not contain water vapor. The air and equilibrator headspace analyses typically had about 10 % or less humidity. Every 27 h the CO2 signal of the LI-COR model 6262 infrared analyzer was zeroed with the dry CO2-free air and spanned with the highest standard. The water vapor channel was zeroed if it showed a reading of greater than 0.5 mmol mol−1 for the dry CO2-free air. The water vapor channel values are used for a minor correction to the xCO2, and any inaccuracy in the H2O results in <0.5 ppm change in the analyzer output of dry mole fraction of CO2 (xCO2 in parts per million, ppm). xCO2 as calculated and output by the analyzer based on measured CO2 and water vapor levels were recorded. The seawater circulation in the systems was automatically turned off when the ships entered port and back-flushed with fresh water, removing particles from the inline water filter and thereby alleviating clogging issues in the filter and reducing biofouling of water lines, the filter and the equilibrator.

Data from the ALoS and Eqnx were transmitted to shore daily and displayed online in graphical format at (last access: 30 June 2020) and (last access: 30 June 2020), respectively. These updates provided a near-real-time opportunity to look at the response of surface water CO2 to episodic events, such as passage of hurricanes, and yielded timely indications of instrument malfunction that could be remedied by company engineers on board or when the ships returned to port. The instruments have shown agreement to within 1 µatm with other state-of-the-art systems in intercomparison studies (Yukihiro Nojiri, Ocean pCO2 System Intercomparison at Hasaki, Japan, personal communication, February 2009). Several different versions of the instrument have been deployed on the ships over the years, but overall measurement principles and accuracies, estimated at better than 2 µatm (Pierrot et al., 2009; Wanninkhof et al., 2013), were maintained.

2.2.2 Thermosalinograph

Temperature and salinity were measured with a flow-through Sea-Bird SBE45 thermosalinograph (TSG) that was in a seawater flow line parallel to the pCO2 equilibrator. A SBE38 remote temperature probe was situated near the inlet before the pump and was used as the SST measurement. The TSGs and temperature probes were maintained by collaborators from the Marine Technical group at the Rosenstiel School of Marine and Atmospheric Sciences (RSMAS) at the University of Miami (U. Miami). The TSGs on the ships were factory-calibrated on an annual basis. Postcalibrations showed no drift in the temperature sensor but occasionally some drift in the conductivity.

2.2.3 Other instrumentation

The underway effort is part of a larger scientific operation lead by RSMAS called Oceanscope (, last access: 12 April 2020). Additional instrumentation on board the ships includes the Marine-Atmospheric Emitted Radiance Interferometer (M-AERI) to assess surface skin temperature retrievals from a number of radiometers on earth-observation satellites (Minnet et al., 2001). The M-AERIs are Fourier-transform infrared interferometers situated on the deck viewing the sea surface away from the wake. Acoustic Doppler current profilers (ADCPs) mounted on the hull of the ships are used to measure ocean currents.

3 Datasets

3.1fCO2 data

Full details on data acquisition with these systems and calculation of pCO2 and fCO2 can be found in Pierrot et al. (2009).

Postcruise, the xCO2 data were processed by first linearly interpolating each standard measured every ∼4 h to the time of a sample measurement and then recalculating the air and water xCO2 values based on the linear regression of the interpolated standard values at the time of sample measurement. For three cruises, the analyzer output showed negative water vapor values due to the condition of desiccant chemicals and thus yielded erroneous dry xCO2 values. Separate processing routines were developed to correct for these situations (available on request from Denis Pierrot). The postcruise corrected xCO2values were used for calculation of pCO2 and fCO2 as described in the calculation section.

3.2 Thermosalinograph, sea surface temperature and salinity data

The SST data were obtained from a temperature probe (Sea-Bird, model SBE38) near the intake. The salinity was determined with a thermosalinograph (Sea-Bird, model SBE45) from the measured conductivity and temperature in the unit using the internal software of the SBE45. The SST and SSS data were appended to the pCO2 data records in real time and also logged via another shipboard computer at more frequent intervals. The SSS data were not quality-controlled, and no corrections to the SSS data were made other than the removal of spikes and values that were out of range (<5 and >40). As salinity has a minimal effect on the calculated fCO2, bad or missing salinities were removed and substituted by linearly interpolated values to eliminate gaps. When SST data were not recorded in the CO2 files and in the event of errors, the SST gaps were filled from the high-resolution SST data files maintained by the RSMAS Oceanscope project. On the rare occasion (<0.1 %) that SST was not recorded at all, SST data were estimated from the equilibrator temperature data (Teq) after applying a constant offset between Teq and SST using SST data before and after the gap. On average the Teq was 0.12±0.28C lower than SST for all the cruises. Per cruise the standard deviation (SD) of the difference of Teq and SST was on the order of 0.04 C. The in situ fCO2 data calculated with Teq were flagged with a WOCE quality control flag of 3, which refers to data that are deemed questionable and have a larger uncertainty (, last access: 12 April 2020).

For the regionally mapped products on a 1 grid and monthly timescale (1 by 1 by month), SST and SSS were obtained from the following sources: the SSS data were from a numerical model, the Hybrid Coordinate Ocean Model (HYCOM;, last access: 12 April 2020), and are referred to as SSSHYCOM. The SST product for the region was the optimum interpolated SST (OISST) from (last access: 12 April 2020) (Reynolds et al., 2007). It uses data from ships, buoys and satellites to generate the fields. For the OISST the reference SST is from buoys, and the other SST data obtained from ships and other platforms were adjusted to the buoy data by subtracting 0.14 C in the OISST product (Reynolds et al., 2007).

3.3 Wind speed data

Winds were measured on the ships, but these data are not used as they are not synoptic for the whole region, which is a requirement for the regional flux maps. Instead, wind speeds were obtained from the updated cross-calibrated multiplatform wind product (CCMP-2; Atlas et al., 2011). The mean scalar neutral wind at 10 m height <u10> and its second moment <u102> were used to calculate the fluxes. They were determined from the 1∕4, 6-hourly product that was obtained from remote-sensing systems (RSSs;, last access: 6 June 2020). This product relies heavily on the European Centre for Medium-Range Weather Forecasts (ECMWF) assimilation scheme that uses in situ and remotely sensed assets, particularly (passive) radiometers on satellites. The directional component uses scatterometer data. The 1∕4, 6-hourly CCMP-2 product was binned and averaged in 1 grid boxes at monthly scales (1 by 1 by month). In the absence of CCMP-2 data for 2018, the wind product from European Reanalysis, ERA5, was used (, last access: 6 June 2020. Copernicus Climate Change Service, 2017). The ERA5 wind data are at 31 km and 3-hourly resolution but were binned and averaged in the same manner as the CCMP-2 winds. There were no apparent biases between the scalar winds for the two products in the Caribbean.

3.4 Mixed-layer depth data

No MLD determinations were made on the cruises, and limited observational estimates from other sources are available. The MLDs provided here are from the same numerical model (HYCOM) as used for the mapped SSS and obtained from (last access: 12 April 2020). The MLDs are based on a density contrast of 0.03 between surface and subsurface. MLDs are used as an independent variable in the MLRs to map the fCO2w values. They are also needed if mixed-layer dissolved inorganic carbon (DIC) inventories are desired and to determine the effect of mixed-layer depth on the changes in fCO2w. As shown in Wanninkhof et al. (2019b, Fig. 9), MLDs are negatively correlated with fCO2w. They are provided in the gridded mapped product.

4 Calculations

The calculations of the concentrations and fluxes follow standard procedures as described below. The mapping procedures are detailed. The section on uncertainty includes an example of different possible means of mapping and the effect on the final product.

4.1 Calculation of pCO2

The starting point in the calculations, which were aided by the use of MATLAB routines following the procedures as in Pierrot et al. (2009), was calibrated (dry) xCO2 values.

The xCO2 values were converted to pCO2eq (µatm) values:

(1) p CO 2 eq = xCO 2 eq ( P eq - pH 2 O ) ,

where eq refers to equilibrator conditions. Peq is the pressure in the equilibrator headspace, and pH2O is the water vapor pressure calculated according to Eq. (10) in Weiss and Price (1980). The pCO2eq was corrected to surface water values using the intake temperature (SST) and the temperature of water in the equilibrator (Teq) according to the empirical relationship that Takahashi et al. (1993) developed for North Atlantic surface waters:

(2) p CO 2 w = p CO 2 eq e ( 0.0423 ( SST - T eq ) ) .

This empirical correction for temperature is widely used, but it is of note that applying the thermodynamic relationships for carbonate dissociation constants yields different temperature dependencies that are a function of temperature. For the average SST in the Caribbean of 27.0 C the coefficient of temperature dependence varies from 0.036 to 0.040 using commonly used constants as provided in inorganic carbon system programs such as CO2SYS (Pierrot et al., 2006) compared to the coefficient of 0.0423 (or 4.23 % C−1) used above. On average the difference between SST and Teq is 0.12 C for all the cruises such that the correction from Teq to SST using the coefficient of 0.0423 in Eq. (2) is 1.9 µatm under average conditions of SST =27C and pCO2w =374µatm. Using a temperature coefficient of 0.036 the temperature correction would be 1.6 µatm, or a 0.3 µatm difference.

4.2 Calculation of fCO2 in air and water

The fCO2 is the pCO2 corrected for nonideality of CO2 solubility in water using the virial equation of state (Weiss, 1974). The correction can be expressed as fCO2a,w = eg(SST,P) pCO2a,w and

(3) g ( T , P ) = [ ( - 1636.75 + 12.0408 T - 0.0327957 T 2 + 0.0000316528 T 3 ) + 2 ( 1 - xCO 2 10 - 6 ) 2 ( 57.7 - 0.118 T ) ( P / 1013.25 ) ] / ( 82.0575 T ) ,

where T is in kelvin, xCO2 is in parts per million, and P is in millibars.

Under average conditions in the Caribbean, the function eg(SST,P)0.997 and fCO2w will be ≈1.2µatm less than pCO2w. As the corrections from partial pressure to fugacity in air and water are approximately the same, the difference between ΔpCO2 (= pCO2w-pCO2a) and ΔfCO2 (= fCO2w-fCO2a), which were used to determine the fluxes (Eq. 5), is negligible (<0.1µatm).

4.3 Gridding procedure

Gridding of the observations of fCO2w, SST and SSS was performed by binning and averaging the data in (1 by 1 by month) cells. At typical ship speeds of 22 kn, the ship covered 1 in about 2.5 h, taking 60 measurements. This would yield about 250 measurements per month assuming weekly cruises through the area. The actual number of measurements per grid cell ranged from 8 to 500. The higher number of observations per cell was mostly in the latter part of the record, when the Eqnx and ALoS operated in the same area. The total number of observations from March 2002 through December 2018 was 1.13 million, and the total number of (1 by 1 by month) grid cells with observations was 9224.

The gridding facilitated the colocation and merging of products such as MLDHYCOM, SSSHYCOM, OISST, <u2>CCMP and marine boundary layer xCO2 (xCO2aMBL) into the gridded observational dataset for further interpretation. The gridded data aided comparison of in situ SST with OISST and SSS with SSSHYCOM. The average difference between SSS and SSSHYCOM for the 2002–2018 data was -0.1±0.28 (n=9224), and for SST and OISST the difference was 0.25±0.40C (n=9224), with the in situ SSS being lower and SST being higher. While both differences include zero within their standard deviation, the temperature difference is in agreement with expected cooler near-surface temperatures that could lead to lower fCO2w values, which in turn have a large impact on the calculated air–sea fluxes.

Table 1Coefficients for the MLR for each year*.

* These annual regressions were used to create the mapped fCO2wMLR fields using the 1 by 1 by month gridded data product. fCO2wMLR = a longitude + b latitude + c OISST + d MLDHYCOM + e SSSHYCOM+f. The second row (in italics) for each annual entry is the error of the coefficient.

Download Print Version | Download XLSX

4.4 Mapping procedures for fCO2w and fluxes

4.4.1 Mapping fCO2w using a multilinear regression

The gridded observations (1 by 1 by month) represented about 10 % of the area and months of investigation from 15 to 28 N and 88 to 62 W over the period of investigation. To interpolate (“map”) the data in space (334 grid cells) and time (192 months) a multilinear regression (MLR) approach was used to determine fCO2w in each grid cell. For each year from 2002 through 2018 the gridded fCO2w observations were regressed against the 1 monthly gridded values of position (lat and long), SST, MLD and SSS. Other permutations of independent parameters were tested but yielded less robust fits (see Appendix A). Annual MLRs were created as fCO2w levels change over time in response to increasing atmospheric CO2 levels. If fCO2w kept pace with atmospheric CO2 increase, this would translate to a linear trend of fCO2w of 2.13 µatm yr−1 over the time period. As shown in Wanninkhof et al. (2019b), multiyear (>4-year intervals) trends varied from −4 to 4 µatm yr−1, with an average linear trend in fCO2w of 1.4 µatm yr−1 from 2002 to 2018.

The annual MLRs were created of the form

(4) f CO 2 w = a lon + b lat + c OISST + d MLD HYCOM + e SSS HYCOM + f ,

with the coefficients for each year along with the standard error in fCO2w and standard error in each of the coefficients provided in Table 1. The standard error in the fCO2wMLR ranged from 5 to 9 µatm for each year. Using the locations, OISST (the optimal interpolated SST product), MLDHYCOM and the SSSHYCOM (output of the HYCOM model) the fCO2wMLR are determined for each grid cell. There were significant cross-correlations between independent variables such that the effect of the often significant year-to-year differences in coefficients was difficult to interpret. However, the fCO2wMLR from the annual MLRs faithfully reproduced the trends and variability.

Figure 3The difference between fCO2wMLR for subsequent months plotted versus time. The solid squares are the differences between December and January where different MLRs are used.


The MLRs were produced for each year such that the mapped products could be extended for future years in a straightforward fashion. To determine if there were anomalous discontinuities between subsequent annual MLRs that could impact the time series, the difference between fCO2w for subsequent months was plotted versus time in Fig. 3. No significant discontinuities were observed between December and January. Only for December and January 2008/2009, 2009/2010 and 2016/2017 do there appear to be slight differences in the pattern of monthly progressions, but such anomalies are observed during other times of year as well. Using an MLR that includes year as one of the coefficients (see Appendix A, Eqs. A3 and A4) provides a slightly worse fit. Moreover, using such a fit would necessitate recalculating the mapped products every time a new year is added.

To illustrate the differences between gridded and mapped products, the results of the fCO2wMLR calculated with the MLRs for 2004, 2011 and 2017 (Table 1) are plotted in Fig. 4a along with the gridded observations for the grid cells that span the longitude range from 88 to 62 W between 23 and 24 N (see Fig. 1). Figure 4a shows that the mapped product using annual MLRs shows the increases in fCO2w over time in the region and consistent differences in patterns between the east and the west for the 3 years. The mapped product showed a reasonable correspondence with the gridded observations. Some of the differences between the gridded observations and mapped product were caused by the mismatch between SST and SSS in situ with the OISST and SSSHYCOM (Fig. 4b and c). In particular, the strong minima in OISST at 79 W are not seen in the SST. This is likely because the OISST captured the lower SST near the coast of Cuba. This caused the fCO2wMLR product to be lower as well, as shown in Fig. 4a. It illustrates that the mapped fCO2wMLR product is influenced by both the annual MLR and the gridded MLDHYCOM, OISST and SSSHYCOM.

Figure 4Zonal section of gridded and mapped products between 23 and 24 N for January 2004 (black), January 2011 (blue) and January 2017 (red). (a) fCO2w; (b) SST; (c) SSS. The lines with small solid circles are the mapped product; the larger solid circles are the gridded data, with standard deviation of data in the box shown as error bars.


4.4.2 Determining the air–sea CO2 fluxes for the region

For the determination of the air–sea CO2 flux (FCO2, mol m−2 yr−1), a bulk formulation was applied using the gridded mapped product

(5) F CO 2 = k K 0 Δ f CO 2 ,

where ΔfCO2 is (fCO2w-fCO2a), K0 is the seawater CO2 solubility that is a function of temperature and salinity (Weiss and Price, 1980), and k is the gas transfer velocity parameterized as a function of wind speed (Wanninkhof, 2014).

(6) k = 0.251 < u 10 2 > ( S c / 660 ) - 1 / 2 ,

where <u102> is the monthly second moment of the wind speeds reported in CCMP-2. The second moment accounts for the impact of variability of the wind speed on k. It is determined by taking the monthly average of the sum of squares of the wind speed in CCMP-2 provided at 6 h and 1∕4 grid resolution. The number of wind speed observations in a (1 by 1 by month) grid is 1920. This sample density captures the frequency spectrum of winds in the grid boxes except that extreme wind events such as hurricanes are not fully represented due to the local nature of the extremes and inherent smoothing in the CCMP-2 product. Sc is the Schmidt number of CO2 in seawater, defined as the kinematic viscosity of seawater divided by the molecular diffusion coefficient of CO2. It is determined as a function of temperature from Wanninkhof (2014). At the average temperature of 27 C the ScCO2 equals 475. Over the typical range of SST in the Caribbean from 24 to 30 C, the (Sc∕660)-1/2 will vary from 1.1 to 1.27, indicating that the gas transfer velocity will be 27 % higher at an SST of 30 C compared to an SST of 20 C that would correspond to an Sc of 660.

The annual fCO2wMLR values as well as the modeled and remotely sensed products OISST, MLDHYCOM and SSSHYCOM were determined for all (1 by 1 by month) grid cells. The monthly fCO2aMBL values were derived from the weekly average xCO2aMBL of the stations on Key Biscayne (KEY) and Ragged Point Barbados (RPB; CarbonTracker Team, 2019;, last access: 20 June 2020). The second moments of the scalar winds <u2> from CCMP-2 were averaged on the same grids. As no CCMP-2 product was available in 2018, the ERA5 wind product was used for the last year of the record.

4.5 Gridded and mapped products for alkalinity and pHT

Total alkalinity (TA) was determined from salinity. For estimation of TA, several algorithms have been developed with salinity (Fig. 5; Millero et al., 1998; Lee et al., 2006; Takahashi et al., 2014 and Cai et al., 2010). The relationship of Cai et al. (2010), TA =57.3 SSS +296.4, SD =5.5, was used as this was determined from observations that are similar to the conditions in the Caribbean.

Figure 5Relationships of surface water TA with salinity. The relationship of Cai et al. (2010), with standard error shown as error bars, is used to calculate pH and ΩAr from fCO2w.


The pHT, the pH on the total scale at SST, was subsequently determined from the calculated TA and fCO2w. For the gridded products the gridded SSS was used to determine TA. For the pHT the gridded TA, SST and fCO2w were used. For the mapped products the gridded OISST and SSSHycom were applied in the calculations. The pHT was calculated using the program CO2SYS for Excel V2.2 (Pierrot et al., 2006) with the apparent CO2 dissociation constants K1 and K2 from Lueker et al. (2000), the KSO4- dissociation constant from Dickson (1990), the KF dissociation constant from Perez and Fraga (1987), and the total boron relationship with salinity from Uppström (1974).

Table 2Estimated uncertainties in measured, gridded and mapped variables.

a Method refers to either the individual data point (measured), the values averaged on a 1 by 1 by month grid, or the interpolated products (mapped). b The systematic uncertainty is based on the difference between the different methods and products as indicated. Note that the combined standard uncertainties and systematic uncertainties are based on average conditions.

Download Print Version | Download XLSX

Table 3Column headers, units and description for the monthly 1 gridded observational product (1 by 1 by month).

Download Print Version | Download XLSX

4.6 Aragonite saturation state (ΩAr)

The aragonite saturation state (ΩAr) indicates the level of supersaturation or undersaturation of seawater with respect to the mineral aragonite, a polymorph of calcium carbonate, and part of the skeletal structure of many marine calcifiers. When ΩAr is less than 1, aragonite dissolution is thermodynamically favored, and when ΩAr is greater than 1, it has a tendency to precipitate. The ΩAr is used as an indicator of ecosystem health with regards to ocean acidification (Mollica et al., 2018). In warm tropical regions surface water saturation states are well above 1, but no active precipitation takes place except under unusual circumstances in shallow waters, in a precipitation process called whitings (Purkis et al., 2017). The ΩAr is not measured directly and is defined as the product of calcium and carbonate ion concentrations divided by the solubility product of aragonite:

(7) Ω Ar = Ca 2 + CO 3 2 - K Ar sp - 1 ,

where Ca2+=0.02128/40.087×(SSS/1.80655)=293.84 SSS (in moles per kilogram of seawater) is the total calcium concentration, derived from salinity (Riley and Tongudai, 1967). CO32- is the total carbonate ion concentration determined from two of the inorganic carbon system parameters, and KArsp is the apparent solubility product of aragonite in seawater at a specified salinity, temperature and pressure. In this work CO32- is determined from the gridded fCO2w and calculated TA using the CO2SYS program. For surface waters, KArsp is

(8) p K Ar sp = - [ - 171.945 - 0.077993 T + 2903.293 / T + 71.595 log ( T ) + ( - 0.068393 + 0.0017276 T + 88.135 / T ) S 0.5 - 0.10018 S + 0.0059415 S 1.5 ] ,

where pKArsp=-logKArsp, T is temperature in kelvin (K) and S is salinity (Mucci, 1983). As with pHT, the ΩAr-gridded product was determined from the gridded (1 by 1 by month) values of SSS, SST, fCO2w and TA, while the mapped product uses SSSHYCOM, OISST, fCO2wMLR and TA.

4.7 Uncertainty of observations and the gridded and mapped products

The uncertainty of the products is difficult to quantify due to the many factors, calculations and interpolations influencing the overall uncertainty. Moreover, the uncertainty estimate includes a random and a systematic component. The latter can have a large influence on interpretations, particularly on the calculated air–sea CO2 fluxes. Below we address the uncertainty in terms of standard errors in the observations, in the gridded products and in the mapped products. For the calculated quantities and nomenclature we follow the approach in Orr et al. (2018). The standard uncertainty is characterized by its standard deviation or error of the measured quantities and the standard uncertainties of the input variables. The propagated uncertainty of a calculated variable is called the combined standard uncertainty. Error propagation is (eA2+eB2+)0.5 for addition and ((eA/A)2+(eB/B)2+)0.5 for multiplication, where A and B are the variables, and eA is the standard error in variable A.

The individual measurements of fCO2w have a combined standard uncertainty of less than 2 µatm based on a propagation of error of instrument response, equilibrator efficiency standardization, and temperatures and pressures at equilibration and at the sea surface (Pierrot et al., 2006). The performance and output data from the UWpCO2 systems have been checked by the manufacturer, in intercomparison exercises and at sea. At the point of measurement, SST measurements are accurate to 0.02 C and SSS measurements to 0.1 based on instrument specifications and annual calibrations.

The combined standard uncertainty in the gridded product will vary based on the number of measurements. It includes the actual variability in the 1 by 1 cells. To estimate the uncertainty per grid cell, the standard deviation of the fCO2w in each cell was determined, and then the average of the standard deviations for the 9924 cells with observations was taken. The average standard deviation was 3.4±2.6µatm (n=9224). The same procedure was followed for SST and SSS and yielded values of 0.22±0.19C for SST and 0.10±0.10 for SSS. These were relatively small uncertainties compared to the monthly spatial range of ≈20µatm for fCO2w, ≈1C for SST and ≈1 for SSS. The amplitudes of the seasonal cycle of ≈40µatm for fCO2w and ≈4C for SST were significantly greater than the standard uncertainties as well.

The calculated parameters for the gridded products, TA, pHT and ΩAr have an added uncertainty due to the uncertainty in the constants and parameterizations. The agreement between the TA–SSS relationships was good, and choice of TA relationship did not have a determining influence on results. We used TA =57.3 SSS_OBS +296.4 specifically developed for the subtropical western Atlantic (see insert of Fig. 10 in Cai et al., 2010). The standard error for the relationship was 5.5 µmol kg−1. The average standard deviation for salinity in the grid cells is 0.1, which translates to an uncertainty in gridded TA of 5.7 µmol kg−1. Thus the combined standard uncertainty for TA is 5.52+5.720.5=7.9µmol kg−1.

pHT is determined from fCO2w and calculated TA with associated uncertainties. An added uncertainty for the calculated pHT is the uncertainty in the dissociation constants that are used to calculate pHT. These uncertainties can be calculated using a modified version of CO2SYS (Orr et al., 2018). As shown in Orr et al. (2018, Fig. 12), the uncertainty in the constants dominates the calculated pHT. Using the uncertainty in the constants as presented in the program, an uncertainty of gridded TA of 8 µmol kg−1 and 3.4 µatm for fCO2 yields a combined standard uncertainty in pH of 0.0075. For comparison, the uncertainty in pH would be 0.0070 if state-of-the-art measurement uncertainties in TA of ±3µmol kg−1 and ±2µatm for fCO2 are used. Similarly, the uncertainty of ΩArusing the same approach and uncertainties in independent parameters are ±0.20. Again it is the uncertainty in the dissociation constants that dominate the overall uncertainty.

For the mapped products there is additional uncertainty due to the use of regressions. The standard errors of the annual regressions in fCO2wMLR and the coefficients of the independent parameters are provided in Table 1. The average standard error for the fCO2wMLR for the 17 years is 6.4±1.2µatm. Following the approach above for the other independent parameters and propagating these uncertainties yields a combined standard uncertainty of 0.0090 in mapped pH and 0.21 in mapped ΩAr.

Systematic errors are introduced by the different SST and SSS data that are used for the gridded and mapped products, with the former using the gridded measured SST and SSS and the latter using the SSSHYCOM and OISST products. The magnitude of the systematic uncertainty between the mapped and gridded product is estimated from the difference in the gridded and mapped parameters for cells that have both gridded and mapped products (for columns headers, see Table 2). Another possible source of systematic uncertainty is using the gridded and mapped fCO2sw, SSS and SST to calculate the TA, pHT and ΩAr rather than using the in situ values to calculate the parameters and then gridding them. This uncertainty is small based on the low uncertainty of variables in each cell as shown from their standard deviations provided in the gridded products (for column headers, see Table 3).

This is confirmed by comparing this method with the approach of calculating the pH and then gridding and mapping as done in Lauvset et al. (2016) and Jiang et al. (2015). To examine the differences derived from using one or the other approach, both were compared for 2017 data. The pH and TA were calculated for every fCO2w observation and binned into a (1×1× mo) grid. An MLR from the calculated pH was created for 2017:

(9) pH T ( ± 0.005 ) = 0.0003194 long - 0.00046744 lat - 0.00965183 SST + 0.00019602 MLD + 0.00069378 SSS + 8.3240 , r 2 = 0.89 ( n = 1244 ) .

This MLR was then applied to the independent variables for each grid box to determine pHT(MLR). This was compared to the approach used here of calculating the pH using the mapped fCO2wMLRand TA–SSS relationships on (1×1× mo) grids, called pHT(fCO2w,TA). The two approaches provided similar results, with pHT(fCO2w,TA) – pHT(MLR) =-0.0001±0.005 for 2017. The small difference showed a pattern with SST (Fig. 6) but not with the other independent variables. The differences using either pHT(fCO2w,TA) or pHT,MLR are an order of magnitude smaller than the combined standard uncertainties such that the approach of using mapped fCO2w and TA to determine pHT yielded precise and consistent gridded pHT.

Figure 6pHT calculated from gridded fCO2w and TA estimated from SSS – pHT(fCO2,TA) – as done in this work minus pHT calculated from observed fCO2w and TA from SSS – pHT(MLR) – plotted against temperature for the 2017 data.


Table 4Column headers, units and description for the monthly 1 mapped product (1 by 1 by month) for the whole region.

Download Print Version | Download XLSX

Table 5Column headers, units and description for the monthly averaged mapped product.

Download Print Version | Download XLSX

Table 6Column headers, units and description for the annual averaged mapped product.

Download Print Version | Download XLSX

Calculated air–sea CO2 fluxes have a significant uncertainty as they are driven by relatively small air–water concentration differences. Using the uncertainty in fCO2wMLR of 6.4 µatm, an uncertainty in fCO2a of 1 µatm, and uncertainties in k of 20 % (Wanninkhof, 2014) and K0 of 0.002 (Weiss, 1974), the corresponding combined standard uncertainty is 21 % for the flux. For flux calculations the systematic error, or bias, is a big issue. The near-surface temperature gradient and skin temperatures will have an impact on the fCO2w and K0. The magnitude, along with its applicability to the bulk flux formulation, is under debate (McGillis and Wanninkhof, 2006; Woolf et al., 2016). For the mapped product, the OISST is used. The OISST uses a variety of temperature data, including remote sensing of the skin temperature, and the product is adjusted to buoy temperatures nominally at 1 m (Reynolds et al., 2007) such that a common reference depth is implicitly used for fCO2wMLR. As shown in Wanninkhof et al. (2019b) the OISST was on average 0.25 C lower than SST, and using the SST instead of the OISST would change the flux from −0.87 to −0.63 mol m−2 yr−1, a 27 % decrease. Using OISST for bulk temperature and a canonical value of 0.17 C for the difference in bulk and skin temperature would increase the uptake of CO2 by 18 %, from −0.87 to −1.04 mol m−2 yr−1. Based on current knowledge, we believe that using the OISST and the resulting calculated fCO2wMLR as was done here yields the appropriate fluxes.

5 Datasets and data products

Several different data products are provided in conjunction with this paper. The methodology to create the products is described above, and the file format and column headers are presented here with a brief description when warranted.

5.1 Underway fCO2 data

The quality-controlled cruise data are posted at different locations. The individual cruise files with metadata can be found at (last access: 30 June 2020). Data can be found as part of the SOCAT holdings (Bakker et al., 2016) using an interactive graphical user interface (, last access: 2 July 2020). In addition, cruise files of the three ships are provided in annual directories at the National Centers for Environmental Information (NCEI;, last access: 12 April 2020). The data file structures are from the MATLAB data reduction program of Pierrot (Denis Pierrot, personal communication, January 2020). The primary identifier for the cruises is the EXPO code, which is the International Council for the Exploration of the Sea (ICES) ship code and the day the ship starts the cruise. Examples are as follows: for a cruise of the EoS starting on 6 March 2002, the EXPO code is 33KF20020306; for an ALoS cruise starting on 25 November 2018, the EXPO code is BHAF20181125; and for the Eqnx cruise departing her home port on 10 February 2018, it is MLCE20180210. The individual cruise files at the sites above sometimes include data outside the study region.

5.2 Gridded data

The gridded datasets are the binned and averaged fCO2w, SST and SSS observations on a (1 by 1 by month) grid. The files include the auxiliary data obtained from remote sensing and interpolated data (OISST), data assimilation of remotely sensed winds (CCMP-2), and from the HYCOM model (SSSHYCOM). TA, pHT and ΩAr calculated using procedures outlined above are provided in the file. The fCO2w values calculated using the annual MLRs (Table 1) are provided as well. This gridded dataset has spatial and temporal gaps as the ships did not transit through each pixel, and coverage is uneven. The number of observations differs for each grid cell and is listed in the tables of the gridded dataset. For the auxiliary data the number of data points is fixed by the resolution of the data products except where part of the grid includes land, which is masked. The column headers are provided in Table 3 and include units and descriptions when warranted.

5.3 Mapped product

The mapped product provides the data in homogeneous (1 by 1 by month) grid boxes utilizing the annual MLRs of fCO2w as a function of lat, long, OISST, SSSHYCOM and MLDHYCOM for the region from 15 to 28 N and -62 to -88 (= 62 to 88 W). The modeled and remotely sensed products OISST, SSSHYCOM and MLDHYCOM as well as position were used in the MLR (Eq. 4) as the independent parameters. The mapped product includes the air–sea CO2 fluxes in the region as a specific flux (mol m−2 yr−1) for each grid box. The column headers are provided in Table 4, including units and descriptions when warranted.

5.4 Monthly and annual estimates for the Caribbean, 2002–2018

Summary files of monthly and annual data and products covering the whole region from 15 to 28 N and from 88 to 62 W are provided based on the average or the sum of the data in the mapped products. The column headers for the monthly and annual products are similar (Tables 5 and 6). For the monthly files the average of each parameter for each cell is area-weighted based on the area of the cell as provided in the mapped product (Table 4) according to cell area / (total area / number of cells). The total CO2 mass flux (CO2_FluxTotal) is the integral of the monthly area-weighted CO2 fluxes (mol m−2 yr−1) expressed in teragrams of carbon (=1012 g C) per month or per year.

6 Data availability

The observations are available at two locations in slightly different formats, but all files are stored by ship and cruise. The data are submitted to SOCAT at least once a year such that they can be posted in the annual updates of SOCAT (, last access: 2 July 2020, Bakker et al., 2016). The permanent depository of the data is at the NCEI, where the data are stored per cruise in directories listed per year (, last access: 2 July 2020, Wanninkhof et al., 2020). The gridded observations and mapped products described herein are posted in directories at the AOML and the NCEI. The dataset and derived quantities are provided on a 1 monthly grid at (last access: 30 June 2020) (; Wanninkhof et al., 2019a). The products cover the years 2002 through 2018 and will be updated annually.

7 Conclusions

The datasets from the cruise ships sailing the Caribbean Sea are a rich resource for studying the trends and patterns of inorganic carbon cycling and ocean acidification in the region. The scales of variability and data density are such that the (1 by 1 by month) monthly gridding captures the magnitudes and trends of fCO2w and derived inorganic carbon products at seasonal to interannual scales. Using annual MLRs to interpolate fCO2w, with position, SST, SSS and MLD as independent variables, yielded accurate monthly products (Wanninkhof et al., 2019a). A comprehensive investigation of the changes in decadal trends based on the dataset and products was presented in Wanninkhof et al. (2019b). The combined standard uncertainties and systematic offsets of gridding and mapping were estimated from comparing fCO2w observations with gridded and mapped products including TA, pHT and ΩAr. The MLRs capture the spatial and temporal variability in fCO2w and calculated pHT and ΩAr well in the region. The datasets and products are invaluable for model initiation and validation and serve as boundary conditions for nearshore fine-scale models.

Appendix A

Multilinear regressions (MLRs) for fCO2w were applied to the data on a (1 by 1 by month) grid using sequentially fewer independent parameters, and the increase in residual was determined. This analysis was performed as two of the independent parameters used – sea surface salinity (SSS) and mixed-layer depth (MLD) – are modeled, and their accuracy is not readily known. The functional form for the multilinear regression (MLR) fit is

(A1) f CO 2 w , MLR = a long + b lat + c SST + [ d MLD ] + [ e SSS ] + f ,

SST is sea surface temperature, long is longitude and lat is latitude. The MLR with the full number of independent parameters is used in the data products described in the main text, and the resulting data products are in Wanninkhof et al. (2019a). The terms in brackets indicate the parameters omitted in the estimates here. The coefficients for the MLRs and their standard errors for each year without MLD as well as the MLRs without MLD and SSS are provided in Tables A1 and A2. SST and location are the strongest predictors of fCO2w levels in the region. The increase in the average root mean square of the residual or error (RMSE) excluding MLD and SSS in the annual MLRs is shown in Tables A1 and A2. The average RMSE of the calculated fCO2w,MLR increases by 8±5.5 % by excluding MLD and increases by 11±5.4 % when MLD and SSS are omitted with the annual differences provided in the last column of Tables A1 and A2.

For the entire record, from 2002 to 2018, single regressions of fCO2w with position, SST, SSS and MLD showed larger standard errors and coefficients of determination (r2) as they do not capture the increase in fCO2w over time due to fCO2w increases in response to increasing atmospheric CO2 levels. Regressing against ΔfCO2, which in principle should not have a trend over time if surface water CO2 levels keep up with atmospheric CO2 increases, did not improve the correlation with the independent parameters. This was attributed to the relatively small magnitude of ΔfCO2 and the observed multiyear changes in trends of fCO2w (Wanninkhof et al., 2019b). However, using the year as a variable in the regression provides a reasonable means of extrapolating over the entire time–space domain with a single regression from March 2002 through February 2018:

(A2) f CO 2 w , MLR ( ± 7.8 ) = 23.3 ( ± 2.0 ) + 1.45 ( ± 0.02 ) ( YR - 2002 ) + 10.23 ( ± 0.06 ) SST + 1.19 ( ± 0.03 ) lat - 0.50 ( ± 0.01 ) long , r 2 = 0.84 ,

where YR is the integer calendar year. Thus (YR-2002) is the year since the start of the record. The coefficient for the (YR-2002) term of 1.45 reflects the annual increase in surface water fCO2w due to atmospheric fCO2a increase, which averages 2.1 µatm yr−1 over the 2002–2018 time period.

The standard error in the fCO2w,MLR of 7.8 µatm in Eq. (A2) is generally greater than the standard error in the annual algorithms used to fill the gaps for the annual estimates (Eq. A1) that range from 4.7 to 9.9 µatm, depending on the year, with an average of 6.4 µatm (Table 1 in main text).

The importance of the different independent variables for the fCO2w,MLR can, in part, be discerned from the standard error in the coefficients, but since variables are cross-correlated, other means are investigated, such as creating an MLR with either a subset or substitution of variables. Physical parameters were correlated with location (lat, long) in the region. In particular, salinity showed broad correspondence with position. Therefore, substituting SSS for lat and long provided a similar magnitude and standard error in the coefficients of the independent variable but a 10 % greater root mean square (rms) in the estimated fCO2w :

(A3) f CO 2 w , MLR ( ± 8.8 ) = - 238 ( ± 9.8 ) + 1.36 ( ± 0.02 ) ( YR - 2002 ) + 10.11 ( ± 0.06 ) SST + 9.07 ( ± 0.25 ) SSS , r 2 = 0.79 .

Finally, a simple two-parameter linear fit with YR and SST had reasonable predictive capability:

(A4) f CO 2 w , MLR ( ± 9.4 ) = 107.2 ( ± 1.7 ) + 1.38 ( ± 0.02 ) ( YR - 2002 ) + 9.50 ( ± 0.06 ) SST , r 2 = 0.77 .

Equation (A4) uses only SST, and time showed an increase in the standard error of the derived independent variable fCO2w compared to the other permutations of the MLR. This simple relationship did show some biases with location (not shown), and for gap filling to create uniform monthly fields of fCO2w, the full annual regressions (Table 1) using all independent parameters – SST, SSS, MLD and location – were the best option.

Table A1Coefficients for the MLR, fCO2w,MLR = a long + b lat + c SST + e SSS + f.

a The second row (in italics) for each annual entry is the error of the coefficient. b The increase in root mean square error (RMSE) of fCO2w,MLR compared to Table 1, which includes MLD as an input.

Download Print Version | Download XLSX

Table A2Coefficients for the MLR, fCO2w,MLR = a long + b lat + c SST + f.

a This table is the same as Table A1 except that MLD in addition to SSS is omitted as an independent variable. b The increase in root mean square error (RMSE) of fCO2w,MLR compared to Table 1, which included MLD and SSS as independent variables. c The second row (in italics) for each annual entry is the error of the coefficient.

Download Print Version | Download XLSX

Author contributions

RW, JT, LB, DP and KS contributed to writing and editing the documents. JT performed most of the gridding and binning and provided the model and remotely sensed data from the sources listed in the text. LB, KS and DP performed maintenance and data reduction and liaised with all parties involved in the operations.

Competing interests

The authors declare that they have no conflict of interest.


This work would not have been possible without support from Royal Caribbean International, who have provided access to their ships and significant financial, personnel and infrastructure resources for the measurement campaign coordinated through the Rosenstiel School of Marine and Atmospheric Sciences (RSMAS) of the University of Miami. Peter Ortner, Elizabeth Williams, Don Cucchiara and Chip Maxwell of the Marine Technical group at RSMAS have been instrumental in maintaining the science operations. David Munro (INSTAR, ESRL/GMD) provided the KEY and RPB CO2 data. NOAA optimal interpolated SST data were provided by the NOAA/OAR/ESRL/PSD. The NOAA Oceanic and Atmospheric Research (OAR) office is acknowledged for financial support, in particular the Ocean Observations and Monitoring Division (OOMD; fund reference no. 100007298) and the NOAA/OAR Ocean Acidification Program.

Financial support

This research has been supported by the NOAA/OAR/OOMD (fund reference no. 100007298).

Review statement

This paper was edited by David Carlson and reviewed by three anonymous referees.


Atlas, R., Hoffman, R. N., Ardizzone, J., Leidner, S. M., Jusem, J. C., Smith, D. K., and Gombos, D.: A cross-calibrated multiplatform ocean surface wind velocity product for meteorological and oceanographic applications, B. Am. Meteorol. Soc., 92, 157–174,, 2011. 

Bakker, D. C. E., Pfeil, B., Landa, C. S., Metzl, N., O'Brien, K. M., Olsen, A., Smith, K., Cosca, C., Harasawa, S., Jones, S. D., Nakaoka, S., Nojiri, Y., Schuster, U., Steinhoff, T., Sweeney, C., Takahashi, T., Tilbrook, B., Wada, C., Wanninkhof, R., Alin, S. R., Balestrini, C. F., Barbero, L., Bates, N. R., Bianchi, A. A., Bonou, F., Boutin, J., Bozec, Y., Burger, E. F., Cai, W.-J., Castle, R. D., Chen, L., Chierici, M., Currie, K., Evans, W., Featherstone, C., Feely, R. A., Fransson, A., Goyet, C., Greenwood, N., Gregor, L., Hankin, S., Hardman-Mountford, N. J., Harlay, J., Hauck, J., Hoppema, M., Humphreys, M. P., Hunt, C. W., Huss, B., Ibánhez, J. S. P., Johannessen, T., Keeling, R., Kitidis, V., Körtzinger, A., Kozyr, A., Krasakopoulou, E., Kuwata, A., Landschützer, P., Lauvset, S. K., Lefèvre, N., Lo Monaco, C., Manke, A., Mathis, J. T., Merlivat, L., Millero, F. J., Monteiro, P. M. S., Munro, D. R., Murata, A., Newberger, T., Omar, A. M., Ono, T., Paterson, K., Pearce, D., Pierrot, D., Robbins, L. L., Saito, S., Salisbury, J., Schlitzer, R., Schneider, B., Schweitzer, R., Sieger, R., Skjelvan, I., Sullivan, K. F., Sutherland, S. C., Sutton, A. J., Tadokoro, K., Telszewski, M., Tuma, M., van Heuven, S. M. A. C., Vandemark, D., Ward, B., Watson, A. J., and Xu, S.: A multi-decade record of high-quality fCO2 data in version 3 of the Surface Ocean CO2 Atlas (SOCAT), Earth Syst. Sci. Data, 8, 383–413,, 2016. 

Boutin, J., Merlivat, L., Henocq, C., Martin, N., and Sallee, J. B.: Air–sea CO2 flux variability in frontal regions of the Southern Ocean from CARbon Interface OCean Atmosphere drifters, Limnol. Oceanogr., 53, 2062–2079, 2008. 

Cai, W.-J., Hu, X., Huang, W.-J., Wang, Y., Peng, T.-H., and Zhang, X.: Alkalinity distribution in the Western North Atlantic Ocean margins, J Geophys. Res., 115, C08014,, 2010. 

CarbonTracker Team: Compilation of near real time atmospheric carbon dioxide data; obspack_co2_1_NRT_v4.4.2_2019-06-10; NOAA Earth System Research Laboratory, Global Monitoring Division,, 2019. 

Copernicus Climate Change Service (C3S): ERA5: Fifth generation of ECMWF atmospheric reanalyses of the global climate. Copernicus Climate Change Service Climate Data Store (CDS), available at: (last access: 9 July 2019), 2017. 

DeGrandpre, M. D., Olbu, G. J., Beatty, C. M., and Hammar, T. R.: Air–sea CO2 fluxes on the US Middle Atlantic Bight, Deep-Sea Res. Pt. II, 49, 4355–4367, 2002. 

Dickson, A. G.: Standard potential of the reaction: AgCl(s) + 1/2 H2(g) = Ag(s) + HCl(aq), and the standard acidity constant of the ion HSO4- in synthetic seawater from 273.15 to 318.15 K, J. Chem. Thermodyn., 22, 113–127, 1990. 

Evans, W., Mathis, J. T., Cross, J. N., Bates, N. R., Frey, K. E., Else, B. G. T., Papkyriakou, T. N., DeGrandpre, M. D., Islam, F., Cai, W.-J., Chen, B., Yamamoto-Kawai, M., Carmack, E., Williams, W. J., and Takahashi, T.: Sea-air CO2 exchange in the western Arctic coastal ocean, Global Biogeochem. Cy., 29, 1190–1209,, 2015. 

Gledhill, D. K., Wanninkhof, R., Millero, F. J., and Eakin, M.: Ocean acidification of the greater Caribbean region 1996–2006, J Geophys. Res., 113, C10031,, 2008. 

Gomez, F. A., Wanninkhof, R., Barbero, L., Lee, S.-K., and Hernandez Jr., F. J.: Seasonal patterns of surface inorganic carbon system variables in the Gulf of Mexico inferred from a regional high-resolution ocean biogeochemical model, Biogeosciences, 17, 1685–1700,, 2020. 

Jiang, L.-Q., R. A. Feely, B. R. Carter,, D. J. Greeley, D. K. G., and, and Arzayus, K. M.: Climatological distribution of aragonite saturation state in the global oceans, Global Biogeochem. Cy., 29, 1656–1673,, 2015. 

Lauvset, S. K., Key, R. M., Olsen, A., van Heuven, S., Velo, A., Lin, X., Schirnick, C., Kozyr, A., Tanhua, T., Hoppema, M., Jutterström, S., Steinfeldt, R., Jeansson, E., Ishii, M., Perez, F. F., Suzuki, T., and Watelet, S.: A new global interior ocean mapped climatology: the 1×1 GLODAP version 2, Earth Syst. Sci. Data, 8, 325–340,, 2016. 

Lee, K., Tong, L. T., Millero, F. J., Sabine, C. L., Dickson, A. G., Goyet, C., Park, G.-H., Wanninkhof, R., Feely, R. A., and Key, R. M.: Global relationships of total alkalinity with salinity and temperature in surface waters of the world's oceans Geophys. Res. Lett., 33, L19605,, 2006. 

Lueker, T. J., Dickson, A. G., and Keeling, C. D.: Ocean pCO2 calculated from dissolved inorganic carbon, alkalinity, and equations for K1 and K2; validation based on laboratory measurements of CO2 in gas and seawater at equilibrium, Mar. Chem., 70, 105–119, 2000. 

McGillis, W. and Wanninkhof, R.: Aqueous CO2 gradients for air-sea flux estimates, Mar. Chem., 98, 100–108, 2006. 

Millero, F. J., Lee, K., and Roche, M.: Distribution of alkalinity in the surface waters of the major oceans, Mar. Chem., 60, 111–130, 1998. 

Minnett, P. J., Knuteson, R. O., Best, F. A., Osborne, B. J., Hanafin, J. A., and Brown, O. B.: The Marine-Atmospheric Emitted Radiance Interferometer: A high-accuracy, seagoing infrared spectroradiometer, J. Atmos. Ocean. Tech., 18, 994–1013,<0994:TMAERI>2.0.CO;2, 2001. 

Mollica, N. R., Guo, W., Cohen, A. L., Huang, K.-F., Foster, G. L., Donald, H. K., and Solow, A. R.: Ocean acidification affects coral growth by reducing skeletal density, P. Natl. Acad. Sci. USA, 115, 1754,, 2018. 

Mucci, A.: The solubility of calcite and aragonite in seawater at various salinities, temperatures, and one atmosphere total pressure, Am. J. Sci., 283, 780–799,, 1983. 

Orr, J. C., Epitalon, J.-M., Dickson, A. G., and Gattuso, J.-P.: Routine uncertainty propagation for the marine carbon dioxide system, Mar. Chem., 207, 84–107,, 2018. 

Perez, F. F. and Fraga, F.: Association constant of fluoride and hydrogen ions in seawater, Mar. Chem., 21, 161–168,, 1987. 

Pierrot, D., Lewis, E., and Wallace, D. W. R.: MS Excel program developed for CO2 system calculations, ORNL/CDIAC-105a, Carbon Dioxide Information Analysis Center, Oak Ridge National Laboratory, U.S. Department of Energy, Oak Ridge, Tennessee,, 2006. 

Pierrot, D., Neil, C., Sullivan, K., Castle, R., Wanninkhof, R., Lueger, H., Johannson, T., Olsen, A., Feely, R. A., and Cosca, C. E.: Recommendations for autonomous underway pCO2 measuring systems and data reduction routines, Deep-Sea Res. Pt. II, 56, 512–522, 2009. 

Purkis, S., Cavalcante, G., Rohtla, L., Oehlert, A., Harris, P. M., and Swart, P.: Hydrodynamic control of whitings on Great Bahama Bank, Geology, 45, 939–942,, 2017. 

Reynolds, R. W., Smith, T. M., Liu, C., Chelton, D. B., Casey, K. S., and Schlax, M. G.: Daily high-resolution blended analyses for sea surface temperature, J. Climate, 20, 5473–5496, 2007. 

Riley, J. P. and Tongudai, M.: The major cation/chlorinity ratios in sea water, Chem. Geol., 2, 263–269,, 1967. 

Schuster, U., McKinley, G. A., Bates, N., Chevallier, F., Doney, S. C., Fay, A. R., González-Dávila, M., Gruber, N., Jones, S., Krijnen, J., Landschützer, P., Lefèvre, N., Manizza, M., Mathis, J., Metzl, N., Olsen, A., Rios, A. F., Rödenbeck, C., Santana-Casiano, J. M., Takahashi, T., Wanninkhof, R., and Watson, A. J.: An assessment of the Atlantic and Arctic sea–air CO2 fluxes, 1990–2009, Biogeosciences, 10, 607–627,, 2013. 

Takahashi, T., Olafsson, J., Goddard, J. G., Chipman, D. W., and Sutherland, S. C.: Seasonal variation of CO2 and nutrients in the high-latitude surface oceans: a comparative study, Global Biogeochem. Cy., 7, 843–878,, 1993. 

Takahashi, T., Sutherland, S. C., Wanninkhof, R., Sweeney, C., Feely, R. A., Chipman, D. W., Hales, B., Friederich, G., Chavez, F., Sabine, C., Watson, A., Bakker, D. C. E., Schuster, U., Metzl, N., Inoue, H. Y., Ishii, M., Midorikawa, T., Nojiri, Y., Koertzinger, A., Steinhoff, T., Hoppema, M., Olafsson, J., Arnarson, T. S., Tilbrook, B., Johannessen, T., Olsen, A., Bellerby, R., Wong, C. S., Delille, B., Bates, N. R., and de Baar, H. J. W.: Climatological mean and decadal change in surface ocean pCO2, and net sea-air CO2 flux over the global oceans, Deep-Sea Res. Pt. II, 2009, 554–577,, 2009. 

Takahashi, T., Sutherland, S. C., Chipman, D. W., Goddard, J. C., Ho, C., Newberger, T., Sweeney, C., and Munro, D. W.: Climatological distributions of pH, pCO2, total CO2, alkalinity, and CaCO3 saturation in the global surface ocean, and temporal changes at selected locations, Mar. Chem., 164, 95-125,, 2014. 

Takahashi, T., Sutherland, S. C., and Kozyr, A.: Global Ocean Surface Water Partial Pressure of CO2 Database: Measurements Performed During 1957–2017 (LDEO Database Version 2017; NCEI Accession 0160492), Version 4.4, NOAA National Centers for Environmental Information, LDEOv2017, 2018. 

Uppström, L. R.: The boron/chloronity ratio of deep-sea water from the Pacific Ocean, Deep-Sea Res., 21, 161–162, 1974. 

Wanninkhof, R.: Relationship between wind speed and gas exchange over the ocean revisited, Limnol. Oceanogr.-Meth., 12, 351–362,, 2014. 

Wanninkhof, R., Bakker, D., Bates, N., Olsen, A., and Steinhoff, T.: Incorporation of alternative sensors in the SOCAT database and adjustments to dataset quality control flags, Carbon Dioxide Information Analysis Center, Oak Ridge National Laboratory, US Department of Energy, Oak Ridge, Tennessee, 25,, 2013. 

Wanninkhof, R., Pierrot, D., Sullivan, K. F., Barbero, L., Triñanes, J. A.: A 17-year dataset of surface water fugacity of carbon dioxide (fCO2), along with calculated pH, aragonite saturation state, and air-sea CO2 fluxes in the Caribbean Sea, Gulf of Mexico and North-East Atlantic Ocean covering the timespan from 2003-03-01 to 2018-12-31, (NCEI Accession 0207749), NOAA National Centers for Environmental Information, Dataset, available at: (last access: 30 June 2020),, 2019a. 

Wanninkhof, R., Triñanes, J., Park, G.-H., Gledhill, D., and Olsen, A.: Large decadal changes in air-sea CO2 fluxes in the Caribbean Sea, J. Geophys. Res.-Oceans, 124, 6960–6982,, 2019b. 

Wanninkhof, R., Sullivan, K. F., and Pierrot, D.: Partial pressure of carbon dioxide, temperature, salinity and other variables collected from surface underway observations during Ships of Opportunity (SOOP) C/S Allure of the Seas lines in the Caribbean Sea, Gulf of Mexico and North Atlantic Ocean in 2019 (NCEI Accession 0184598), NOAA National Centers for Environmental Information, Dataset,, 2020. 

Webb, J. R., Maher, D. T., and Santos, I. R.: Automated, in situ measurements of dissolved CO2, CH4, and δ13C values using cavity enhanced laser absorption spectrometry: Comparing response times of air-water equilibrators, Limnol. Oceanogr.-Meth., 14, 323–337,, 2016. 

Weiss, R. F.: Carbon dioxide in water and seawater: the solubility of a non-ideal gas, Mar. Chem., 2, 203–215, 1974. 

Weiss, R. F. and Price, B. A.: Nitrous oxide solubility in water and seawater, Mar. Chem., 8, 347–359, 1980.  

Wilkinson, M. D., Dumontier, M., Aalbersberg, I. J., Appleton, G., Axton, M., Baak, A., Blomberg, N., Boiten, J.-W., da Silva Santos, L. B., Bourne, P. E., Bouwman, J., Brookes, A. J., Clark, T., Crosas, M., Dillo, I., Dumon, O., Edmunds, S., Evelo, C. T., Finkers, R., Gonzalez-Beltran, A., Gray, A. J. G., Groth, P., Goble, C., Grethe, J. S., Heringa, J., 't Hoen, P. A. C., Hooft, R., Kuhn, T., Kok, R., Kok, J., Lusher, S. J., Martone, M. E., Mons, A., Packer, A. L., Persson, B., Rocca-Serra, P., Roos, M., van Schaik, R., Sansone, S.-A., Schultes, E., Sengstag, T., Slater, T., Strawn, G., Swertz, M. A., Thompson, M., van der Lei, J., van Mulligen, E., Velterop, J., Waagmeester, A., Wittenburg, P., Wolstencroft, K., Zhao, J., and Mons, B.: The FAIR Guiding Principles for scientific data management and stewardship, Scientific Data, 3, 160018,, 2016. 

Woolf, D. K., Land, P. E., Shutler, J. D., Goddijn-Murphy, L. M., and Donlon, C. J.: On the calculation of air-sea fluxes of CO2 in the presence of temperature and salinity gradients, J. Geophys. Res.Oceans, 121, 1229–1248,, 2016. 


fCO2 is the pCO2 corrected for the nonideal behavior of CO2 (Weiss, 1974). In surface water fCO2w≈0.997 pCO2w

Short summary
This paper describes a 17-year dataset of over a million data points of automated partial pressure of CO2 (pCO2) measurements on large luxury cruise ships of Royal Caribbean Cruise Lines (RCCL). These data are used to provide trends of ocean acidification and air–sea CO2 fluxes. The effort was possible through a unique continuing industry (RCCL), academic (University of Miami) and governmental (NOAA) partnership.