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

A high-quality dataset of surface water partial pressure/fugacity of CO2 (pCO2w/fCO2w)1, comprised of over a million observations, and derived products are presented for the Northern Sea covering the timespan from 2002 through 2018. Prior to installation of automated pCO2 systems on cruise ships of the Royal Caribbean Cruise Lines and subsidiaries, very limited 25 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's ocean. The dataset, and derived quantities are provided on a 1-degree monthly grid at http://accession.nodc.noaa.gov/0207749, DOI:10.25921/2swk-9w56 (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 ̊ N to 28 ̊ N and 88 ̊ W to 62 ̊ W. The data and products are used for determination of status and trends of 30 ocean acidification, for quantifying air-sea CO2 fluxes, and for ground truthing models. Methodologies to derive the inorganic carbon system parameters and to calculate the fluxes from fCO2w are described.

. It's unfortunate that there's little data where the salinity variability is presumably the highest (i.e. in the Southern Caribbean where the large South American Rivers affect the region). What is the effect of this on the gridded data?  -I'm not sure of the value of section 1.3.3 unless these data are used in the presented dataset or data products (which is unclear) 185 This is to indicate that this UWpCO2 system is part of a larger effort. This also points readers to a possible opportunity to utilize these observations in conjunction with UWpCO2 data and data products presented -Section 2 could be a subsection under section 1 -Much of the information on lines 65-74 190 would be more appropriate in the methods (much of it is also repeated in the different methods sections) -The information on lines 270-294 would be better suited in section 3.5 -In section 4.1 you give much information which is suitable, and partly repeated, in section 5

As shown in
We've made significant edits along the lines suggested but retained the overall structure. E.g. section 1. Is about instrumentation, section 2 is about data from the 195 instrumentation. For readers who are more interested in the data it avoids needing to closely read the detail of the instruments themselves.

Minor comments:
Line 64: I'd prefer the term "raw data processing" over "data reduction". While the former is commonly used in the community, it is not intuitive to those outside what it actually involves 200 We are confused by this comment as it suggests that the recommended change to "raw data processing" is not intuitive. We agree with that as "raw data" is a bit ambiguous and we kept nomenclature as is.
In the introduction it is stated that the Explorer of the Seas changed her home port to Bayonne, NJ in 2008 while in section 1 it is stated that the new home port is Cape Liberty 205 Cruise Port. I realize these may be in the same place but it is nevertheless confusing.

Thanks for pointing this out it was changed
Please revise Line 171: Explain what flag questionable is (presumably WOCE 3)

Done 210
Line 242: I do not understand the method. Please explain. Line 340-341: While this is correct I find it helpful to instead state that when omega_Ar<1 dissolution is thermodynamically favored, and vice versa when omega_Ar>1. In living organisms both dissolution and precipitation of calcium carbonate is biologically mediated, and shells have been shown to survive well in water with omega_Ar 0.9.

215
Based on this reviewers comments, and that of reviewer one, we have changed the text to "When ΩAr is less than 1 dissolution is thermodynamically favored and when ΩAr is greater than 1 aragonite would have a tendency to precipitate." In section 4 you should define the difference between a dataset and data products. My 220 experience is that surprisingly many do not know the difference. It is unclear whether you consider the gridded data part of the dataset or a data product.
We've adopted following nomenclature [following Wanninkhof  The manuscript describes a 17 year dataset for surface water marine carbonate data collected using multiple ships within the Caribbean and a substantial set of derived data.
The manuscript appears to have been a little rushed. There are instances of unclear 240 statements, inconsistent naming, repetition, use of non-SI units, formatting errors and some structural issues. I have listed all comments referring to these issues under the section 'Minor comments' (See below).
These issues were also brought up by the other reviewers and we have addressed these and/or provided detail about changes (or why not) below. The issue of non-SI 245 units is noted but we use the community accepted terms of expressing fCO2 in µatm and pressure in millibar. This is done almost exclusively in ocean carbon cycle research.
I would suggest that the manuscript is re-considered after revision and my reasoning is 250 explained below within the Major comments.
Major comments: 1. The uncertainty information within the manuscript is inconsistent and/or incomplete. Some information is given for the fCO2 data but nothing is given for the temperature or salinity. No uncertainty information or statements are given for the derived datasets eg pH or aragonite saturation state or gas fluxes. This limited information will limit the 255 use of these data, or could result in users making incorrect assumptions about the uncertainties. It would be good if the authors could follow a standard framework or phrasing for presenting the uncertainty information e.g. BIPM 2008 framework and identifying if uncertainties are Type A or Type B and also identify which components of the uncertainty budget have been considered and which have been ignored. Its clear that the derived datasets 260 are unlikely to be considered to be 'truth' measurements, so the authors need to write some text to explain this, so that users of the dataset don't make the mistake of assuming that these data are truth. It may make sense, and/or make it easier for the reader, if all of the uncertainty information was grouped together into a common location (eg one table?) which can then be referred to within the different sections of the manuscript. ignores all vertical temperature gradients. However the dataset includes OISST data which could be used to perform a more accurate gas flux calculation (e.g. re-calculate pCO2 to a common depth, then perform a more accurate calculation). The authors could either provide the results using a more accurate gas flux calculation or highlight this issue to the user/reader and then refer to them to an example analysis that shows the impact of a lower accuracy gas 285 flux calculation and This is a description of the mapped product using the conventional bulk flux parameterization as described. This is the same approach as in most climatologies. We 295 are aware of the developments and claims of superior ways of determining the fluxes by normalizing to a common depth and using skin temperatures. We have added a short description in the uncertainty section and referenced two key papers addressing the cool skin effect controversy. This is mentioned in the Wanninkhof 2019 paper as well . 300 its not clear why the multi-linear regressions are performed and/or why anyone would want this output. These results and methods should introduced giving an explanation as to why they are useful. I'm not sure that this part of the dataset is needed though.
We have clarified that this is the means of mapping for gap filling.

305
As mentioned:" The gridded observations (1˚ by 1˚ by mo) represent about 10 % of the area of investigation from 15-28 ˚N and 88-62 ˚W over the period of investigation" so there has to be a means to fill in the missing cells. values ie +-1uatm? how has this value of 2uatm been estimated? This is provided in the reference that is now added , "estimated at better than 2 µatm (Pierrot et al., 2009)," The intercomparison showing 1 uatm agreement used common standards, pressure and temperature measurements. The 2 uatm incorporates uncertainties in Teq and P as described in the provided reference. 6. section 1.3.2 the precision, accuracy and sensitivity of these instruments are missing.
We have now stated that these can be found in Pierrot et al (2009 9. line 160. processing routines are mentioned but no detail is given. could an overview of these processing routines be provided in the appendices? this information would appear fairly important should anyone want to use these data and/or try and follow the same methods for a similar effort somewhere else in the world.

390
The processing routines are described in Pierrot et al. (2006) and as listed MATLAB routines are available from Pierrot on request.
10. line 266, I'd suggest 'While both differences include zero within their uncertainty..' Done 11. Suggest that section 3.4 (binning procedure) comes before the sections on the calculations 395 (as surely the binning is done first, then the calculations are performed). We've rearranged the sections and clarified procedures. The fCO2 and pCO2 was calculated and then the results were binned. For the fluxes we used the binned product.
So placing it before the fluxes would make sense but it does not flow as well. Added a reference to the following section "a bulk formulation is applied to the data from the 400 gridded mapped product (see 3.4):" 12. line 303, I think that CCMP data are available for 2018 eg http://data.remss.com/ccmp/v02.0/ Yes but this was not available when the calculations were performed. These will be 405 annual updates and when the 2019 data is included the winds will be updated

Introduction
Over the past 20 years a rapidly expanding program of measurements of surface water partial pressure of carbon dioxide (pCO2w) 1 , or the 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., 2012;Takahashi 495 et al., 2014;Wanninkhof et al. 2019 a, b). Marginal seas, that 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 Cruise Lines (RCCL), 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 RCCL ship Explorer of the Seas (EoS) was equipped with an underway pCO2 setup providing observations on 500 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. The surface water pCO2 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 carbon atlas (SOCAT) (Bakker et al., 2016) and global CO2 climatology (Takahashi et al., 2009;. 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 530 changing fCO2w. This effort has made the Northern Caribbean one of the few places in the world's 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 Caribbean ocean acidification product suite that maps ocean acidification conditions in the Caribbean (Gledhill et al., 2008; https://www.coral.noaa.gov/accrete/oaps.html; 535 https://cwcgom.aoml.noaa.gov/erddap/griddap/miamiacidification.graph). 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 the findable, accessible, interoperable and reusable (FAIR) principles (Wilkinson et al., 2016). The documentation of the cruises, the sampling methodology, and data reduction techniques are presented in brief. This is 540 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 on monthly scales, referred to as gridded observations. Then the so-called second inorganic carbon system parameter is calculated that is used in conjunction with fCO2w. 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 a software program CO2SYS (Pierrot et al., 2006) to calculate the other inorganic carbon 545 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 multi-linear regressions (MLR) are developed between the gridded fCO2w data, and sea surface temperature (SST), sea surface salinity (SSS), location (latitude (Lat )and longitude (Lon)), and mixed layer depth as independent variables. These regressions are applied at monthly and 1˚ by 1˚ spatial resolution to the region between 15˚ N and 28˚ N, and 88˚ W and 62˚ W using remotely sensed or modeled independent 550 variables, and used to calculated air-sea CO2 fluxes. These calculated parameters are called gridded mapped products. The procedures and datasets, including an uncertainty analyses, and tables of column headers of the product created are provided below.

Observations
The observational program is described in terms of the ships, voyages, and instrumentation. The ships predominantly sailed 555 in the Caribbean Sea, but also had tracks outside the region, including in the Northeast of the USA, to Bermuda and in the Mediterranean. The description of operations, data, and products presented here cover the Northern Caribbean Sea and  (Fig. 2). The cruises lasted between 7 and 14 days, and made about a half a dozen ports of call. The ships generally were 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 585 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 non-optimal 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 forward of the bow thruster and had fewer issues with bubble entrainment. On the EoS, the 590 underway pCO2 instrument was initially in a dedicated science laboratory built for purpose 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 of 2008. The Eqnx had an air intake near the bow one level below the main deck since the initial installation in 595 March 2015. The underway pCO2 systems on these ships made marine boundary layer (MBL) air observations (xCO2a) as described in Wanninkhof 2019b but these are not used in the data products presented. Typical cruise speeds were 22 knots which, with sampling every 2.5 minutes, yielded a fCO2w sample approximately every 1.7 km except for 35 minutes every 4.5 hours when four calibration gases, and a CO2-free reference gas were analyzed, followed by 5 atmospheric CO2 measurements for the ships with air intakes. This created a gap of 24 km (≈1/4˚) without fCO2w measurements 600

pCO2 system
The instrumentation is based on a community design described in Pierrot et al. (2009). The instruments were manufactured by General Oceanics Inc. and have performed to high accuracy specifications (Wanninkhof, 2013). The surface water drawn from the intake at about 4 l min -1 went through a 1.1 l sprayhead equilibrator with a water volume of 0.5 l and a headspace 630 of 0.6 l. The spray and agitation caused the CO2 in the headspace to equilibrate with the CO2 in water with a response time of about 2 minutes (Pierrot et al., 2009;Web et al., 2016). Thus, the air in the headspace reached 99.8 % equilibration in 12 minutes. As the air is recirculated and surface waters were relatively homogeneous on hourly timescales, 100. The gas entering the analyzer was dried by passing it through a thermo-electric 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 640 humidity. Every 27 hours the CO2 of LICOR model 6262 infrared analyzer were zeroed with the dry CO2-free air and spanned with the highest standard. The water vapor channel was zeroed with the dry CO2-free air but not spanned. The dry mole fraction of CO2 (xCO2 in parts per million, ppm), as calculated and output by the analyzer based on measured CO2 and water vapor levels, were recorded. The systems were automatically turned off when the ships entered port, and upon shutdown they were back-flushed with fresh water removing particles from the inline water filter thereby alleviating 645 clogging issues in the filter and reducing biofouling of water lines, filter and equilibrator.
Data from the ALos and Eqnx were transmitted to shore updated daily and displayed online in graphical format at https://www.aoml.noaa.gov/ocd/ocdweb/allure/allure_realtime.html and https://www.aoml.noaa.gov/ocd/ocdweb/equinox/equinox_realtime.html, respectively. These updates provided a near-real 650 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 when the ships returned to port. The instruments have shown agreement to within 1 µatm with other state-of-art systems in intercomparison studies (Nojiri et al., Feb. 2009 pers. com.). 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. 655 Deleted: Typical cruise speeds are 22 knots which, with sampling every 2.5 minutes, yields a sample approximately every 1.7 km except for 35 minutes every 4.5 hours when four calibration gases and a CO2-free reference gas are analyzed, followed by 5 670 atmospheric CO2 measurements for the ships with air intakes. This creates a gap of 24 km (≈1/4˚) without fCO2w measurements. …

Other instrumentation
The underway effort is part of a larger scientific operation lead by RSMAS called Oceanscope (https://oceanscope.rsmas.miami.edu/). Additional instrumentation onboard the ships include 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-AERI's are Fourier-Transform Infrared interferometers situated on 705 the deck viewing the sea surface away from the wake. Acoustic Doppler Current Profilers (ADCP) mounted on the hull of the ships are used to measure ocean currents.

pCO2 data 710
Full details on data acquisition with these systems and calculation of pCO2 and fCO2 can be found in Pierrot et al. (2009)..
Post-cruise, the xCO2 data were processed by first linearly interpolating each standard measured every ~4 hours 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 3 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 715 routines were developed to correct for these situations (available on request from D. Pierrot). The post-cruise corrected xCO2 values were used for calculation of pCO2 and fCO2 as described in the calculation section.

Thermosalinograph, sea surface temperature and salinity data
The SST data were obtained from a temperature probe (Seabird, model SBE38) near the intake. The salinity was determined with a thermosalinograph (Seabird, model SBE45), from the measured conductivity and temperature in the unit using the 720 internal software of the SBE45. The SST and SSS data were appended to the pCO2 data records in real time and also logged analyses typically had about 10 % or less humidity. Every 27 hours the LICOR model 6262 infrared analyzer was zeroed with the CO2free air and spanned with the highest standard. The dry mole fraction of CO2 (xCO2 in parts per million, ppm), as calculated and output by the analyzer based on measured CO2 and water vapor levels, were via another shipboard computer at more frequent intervals. The SSS data was not quality controlled and no corrections to the SSS data were made, other than removal of spikes and values that were out of range (<5 and > 40). As salinity has minimal effect on the calculated fCO2, bad or missing salinities were removed and substituted by linearly interpolated values to 760 eliminate gaps. When SST were not recorded in the CO2 files, or in error, 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 was 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.28˚C lower than SST for all the cruises. Per cruise the standard deviation (stdev) of the difference of Teq and SST was on the order of 0.04 ˚C. The in situ pCO2 data calculated 765 with Teq was flagged with a WOCE quality control flag of 3 which refers to data that is deemed questionable and has a larger uncertainty.
For the regionally mapped products on a 1-degree grid and monthly timescale (1˚ by 1˚ by mo), SST and SSS were obtained from the following sources: The SSS data were from a numerical model, the Hybrid Coordinate Ocean Model (HYCOM) (https://HYCOM.org/) and referred to as SSSHYCOM. The SST product for the region is the Optimum Interpolated SST, 770 OISST from https://www.esrl.noaa.gov/psd/., (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).

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 775 requirement for the regional flux maps. Instead, wind speeds were obtained from the updated cross-calibrated multi-platform wind product (CCMP-2) (Atlas et al., 2011). The mean scalar neutral wind at 10 m height <u10>, and its second moment <u10 2 > were used to calculate the fluxes. They were determined from the ¼ degree, 6-hourly product that was obtained from Remote Sensing Systems (RSS) (www.remss.com). This product relies heavily on the European Center for Median Weather Forecasting (ECMWF) assimilation scheme that uses in situ and remotely sensed assets, particularly (passive) radiometers 780 on satellites. The directional component uses scatterometer data. The ¼ degree, 6-hourly CCMP-2 product was binned and averaged in 1-degree grid boxes on monthly scales (1˚ by 1˚ by mo). In absence of CCMP-2 data for 2018, the wind product from European Reanalysis, ERA5, was used (https://www.ecmwf.int/en/forecasts/datasets/reanalysis-datasets/era5, 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 785 products in the Caribbean.

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 Deleted: The SST and salinity data, obtained from temperature 790 and conductivity measured in the thermosalinograph, were appended to the pCO2 data records in real-time and also logged via another shipboard computer at a more frequent interval. The thermosalinographs on the ships were factory calibrated on an annual basis. Post-calibrations showed no drift in the temperature 795 sensor but occasionally some drift in the conductivity. No corrections to the sea surface salinity (SSS) were made and salinity did not undergo further quality control other than removal of spikes and values that were out of range.

Moved up [1]:
The thermosalinographs on the ships were factory 800 calibrated on an annual basis. Post-calibrations showed no drift in the temperature sensor but occasionally some drift in the conductivity. No http://www.science.oregonstate.edu/ocean.productivity/index.php. The MLDs are based on a density contrast of 0.03 between surface and subsurface. Mixed layer depths (MLD) 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 (2019b, Figure 9), MLDs are 820 negatively correlated with fCO2w. They are provided in the gridded mapped product.

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. 825

Calculation of pCO2
The starting point in the calculations, which were aided by use of MATLAB routines following the procedures as in Pierrot et al. (2009), were calibrated (dry) xCO2.
The xCO2 values were converted to pCO2eq (µatm) values: where eq refers to equilibrator conditions. The 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: pCO2w = pCO2eq e (0.0423(SST-Teq)) (2) 835 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 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 840 the correction from Teq to SST using the coefficient of 0.0423 in Eq. 2 is 1.9 µatm under average conditions of SST= 27 ˚C and fCO2w =374 µatm. Using a temperature coefficient of 0.036 the temperature correction would be 1.6 µatm, or a 0.3 µatm difference.

Gridding procedure
Gridding of the observations of fCO2w, SST, and SSS was performed by binning and averaging the data in (1˚ by 1˚ by mo) where ∆fCO2 is (fCO2w-fCO2a), s is the seawater CO2 solubility (Weiss and Price, 1980), and k is the gas transfer velocity parameterized as a function of wind speed (Wanninkhof, 2014): where <u10 2 > is the monthly 2 nd moment of the wind speeds reported in CCMP-2. The 2 nd 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-hours and ¼˚ grid resolution.
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 fCO2MLR is 1035 determined for each grid cell.. There were significant cross-correlations between independent variables such that effect of the, often significant, year-to-year differences in coefficients were difficult to interpret. However, the fCO2wMLR from the annual MLRs faithfully reproduced the trends and variability.
The MLRs were produced for each year such that the mapped products could be extended for future years in a 1040 straightforward fashion. To determine if there were anomalous discontinuities between subsequent annual MLRs that could impact the timeseries, the difference in between fCO2w for subsequent months were plotted versus time in  (Table 1) are plotted in Figure 4a along with the gridded observations for the grid cells that span 1050 the longitude range from 88˚W to 62˚ W between 23˚ N and 24˚ N (see Fig.1). Figure 4a shows that the mapped product using annual MLRs show the increases in fCO2w over time in the region, and consistent differences in patterns between the East and the West for the three 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 (Figs. 4b and 4c). In particular, the strong minima in OISST at 1055 79˚ W is 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 Figure 4a. It illustrates that the mapped fCO2wMLR product is both influenced by the annual MLR, and the gridded MLDHycom, OISST, and SSSHycom.  Table 1. The error in the calculated fCO2w ranges from 5 to 9 Deleted: The strongest dependence was with SST.  (Table 1) were applied to the January data in 2004, 2011, and 2017 and are plotted in Figure 3a along with the observations (if available) for the grid boxes that span the longitude range from 88˚W to 62˚ W between 23˚ N and 24˚ N (see Figure 1). Figure 3a 1075 shows that the mapped product using annual MLRs reflected the increases in fCO2w in the region, and consistent differences in patterns between the East and the West for the three years. The mapped product showed a reasonable correspondence with the gridded observations. Some of the deviations between the gridded 1080 observations and mapped product were caused by the mismatch between SST and SSS in situ with the OISST and SSSHYCOM (Figures 3b and 3c). In particular, the strong minima in OISST at 79˚ W were likely due to the OISST capturing the lower SST near the coast of Cuba. This caused the mapped fCO2 product to be lower 1085 as well, as shown in Figure 3a. It illustrates that the gridded fCO2 product is both influenced by the annual MLR, and the MLD, SST, and SSS products used for mapping.…

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 1090 product: 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): where <u10 2 > is the monthly 2 nd moment of the wind speeds reported in CCMP-2. The 2 nd 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-hours and ¼˚ grid resolution. The number of wind speed observation in a (1˚ by 1˚ by mo) 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 1100 product. The 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 ˚C 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 a SST of 20 ˚C that would correspond to a Sc of 660. 1105 The annual fCO2MLR, and modeled and remote sensed products, OISST, MLDHYCOM, SSSHYCOM, were determined for all (1˚ by 1˚ by mo) grid cells. The monthly fCO2aMBL values were derived from the weekly average xCO2aMBL of the stations on https://www.esrl.noaa.gov/gmd/ccgg/flask.php). The second moments of the scalar winds, <u 2 > 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. The second row (in italics) for each annual entry is the error of the coefficient. 1135

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 andCai et al., 2010). The relationship of Cai et al. 2010 , TA= 57.3SSS+ 296.4, stdev = 5.5 was used as this was determined from observations that are similar to the 1140 conditions in the Caribbean.
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 SST and SSS were used to determined TA. For the pHT the gridded TA 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, K2 from Lueker et 1145Lueker et al. (2000; the KSO4dissociation constants from Dickson (1990); thr KF dissociation constants from Perez and Fraga (1987); and the total boron relationship with salinity from Uppström (1974).

Aragonite saturation state (ΩAr) 1150
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. That is, when ΩAr is less than 1 aragonite dissolution is thermodynamically favored and when ΩAr is greater than 1 it has a tendency to precipitate. It 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 one but no active precipitation takes place except under 1155 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: where [Ca 2+ ] is the total calcium concentration and is derived from salinity: [Ca 2+ ] = 0.02128 / 40.087 * (SSS/ 1.80655) = 293.84 S in mol/kg-SW (Riley and Tongudai, 1967).. [CO3 2-] is the total carbonate ion concentration determined from two of the inorganic carbon system parameters, and KAr'sp is the apparent solubility product of aragonite in seawater at a specified  (Millero et al., 1998;Lee et al., 2006;Takahashi et al., 2014 andCai et al., 2010) that show close agreement for the conditions in the Caribbean. Figure 4 shows several TA-salinity relationships over the salinity 1175 range encountered in the region along with error bars for the Cai et al. (2010) relationship depicting the uncertainty of 5.5 µmol kg -1 . The relationship of Cai et al. (2010), specifically developed for the Caribbean Sea (see insert of Figure 10 in Cai et al., 2010), TA= 57.3SSS+ 296.4, stdev = 5.5 was used here. … 1180 Deleted: As shown in Figure 4 the agreement between relationships was good, and choice of TA relationship did not have a 1230 determining influence on results. The calculated TA was used to determine pH with the gridded fCO2w product. The program CO2SYS for excel V2.2 (Pierrot et al., 2006) was used with the apparent CO2 dissociation constants, K1, K2 from Lueker et al. (2000); KSO4dissociation constants from Dickson (1990); KF 1235 dissociation constants from Perez and Fraga (1987); and Total Boron salinity relationship from Uppström (1974). The pH product was determined at the OISST and on the total scale (pHT) on the 1˚ by 1˚ by mo grid. ¶ This approach of calculating the pH from the monthly product 1240 differs from other created pH products (Lauvset et al., 2016;Jiang et al., 2015) where the pH is calculated from the in situ measurements of DIC and TA and then regressed and gridded. To examine the differences derived from using one or the other approach, both were ... Where pKAr'sp = -log KAr'sp , 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 mo) values of SSS, SST, fCO2w and TA, while the mapped product uses SSSHycom, OISST, fCO2wMLR and TA.

Uncertainty of observations, gridded and mapped products 1255
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 1260 quantities and nomenclature we follow the approach in Orr et al. (2018). The standard uncertainty is characterized by its standard deviation/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 for addition is (eA 2 +eB 2 +…..) 0.5 and for multiplication it is ((eA/A) 2 + (eB/B) 2 +…. ) 0.5 , where A and B are the variables, and eA is the standard error in variable A. 1265 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 at manufacturer, in intercomparison exercises, and at sea. SST measurements, at point of measurement, are accurate to 0.02 ˚C 1270 and SSS to 0.1 based on instrument specifications and annual calibrations.
The combined standard uncertainty in 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 stdev of the fCO2w in each cell was determined and then the average of the stdev for the 9924 cells with observations was taken. The average stdev was 3.4 ± 2.6 µatm 1275 (n=9224). The same procedure was followed for SST and SSS and yielded values of 0.22± 0.19 ˚C 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; ≈1 ˚C for  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 TA-SSS relationships was good, and choice of TA relationship did 1290 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 Figure 10 in Cai et al., 2010). The standard error for the relationship was 5.5 µmol kg -1 . The average stdev for salinity in the grid cells is 0.1 which translates in a uncertainty gridded TA of 5.8 µmol kg -1 . Thus combined standard uncertainty for TA is ((5.5/2375) 2 +(5.8/2375) 2 ) 0.5 2375 = 8 µmol kg -1 , where 2375 µmol kg -1 is the average TA. 1295 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, Figure12), the uncertainty in the constants dominates the calculated pH T. Using the uncertainty in the constants as presented in the program, an uncertainty of gridded 1300 TA of 8 µmol kg -1 and 3.4 µatm for fCO2 yields an 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 ΩAr using the same approach and uncertainties in independent parameters is ±0.20. Again it is the uncertainty in the dissociation constants that dominate the overall uncertainty.

1305
For the mapped products there is the additional uncertainty through 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 1310 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 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 (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 1315 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 stdev provided in the gridded products (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, 1320 both were compared for 2017 data. The pH and TA were calculated for every fCO2w observation and binned into a (1˚ x 1˚ x mo) grid. A MLR was from the calculated pH created for 2017: pHT(±0.005) = 0.0003194 Lon -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 1325 the approach used here of calculating the pH using the mapped fCO2wMLR and TA-SSS relationships on (1˚ x 1˚ x 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(MLR) yielded precise and consistent gridded pHT. 1330 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 and an uncertainty 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 1335 temperatures will have an impact on the fCO2w and K0. The magnitude, along with its applicability to the bulk flux formulation 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 implicitly a common reference depth is 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 1340 instead of the OISST would change the flux from -0.87 to -0.63 mol m -2 yr -1 a 27 % decrease. If OISST is used for bulk temperature and using a canonical value for difference in bulk and skin temperature of 0.17 ˚C, it would increase the uptake of CO2 by 18 % from -0.87 to -1.04 mol m -2 yr -1 or 18 %. Based on current knowledge, we believe that using the OISST and the resulting calculated fCO2wMLR as was done here yields the appropriate fluxes. FluxSST-Flux OISST =0.24 mol m -2 yr -1 1. Method refers to either the individual data point (measured); the girded value averaged on a 1˚ by1˚ by mo (gridded); or the interpolated products (mapped) 2. The systematic uncertainty is based on the difference between the different methods and products a indicated. 1350 Note, the combined standard uncertainties and systematic uncertainties are based on average conditions

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 here the file format and column headers is presented with a brief description when warranted.

Underway pCO2 data 1355
The quality controlled cruise data are posted at different locations. The individual cruise files with metadata can be found at https://www.aoml.noaa.gov/ocd/ocdweb/occ.html. Data can be found as part of the SOCAT holdings (Bakker et al. 2016) using an interactive graphical user interface https://ferret.pmel.noaa.gov/socat/las/. In addition, cruise files of the three ships are provided in annual directories at the National Center for Environmental Information (NCEI) (https://www.nodc.noaa.gov/ocads/oceans/VOS_Program/explorer.html). The data file structures are from the MATLAB 1360 data reduction program of Pierrot (pers. comm.).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 March 6, 2002, the EXPO code is 33KF20020330; for a ALos cruise starting

Deleted: AoS
November 25, 2018, the EXPO code is BHAF20181125, and for the Eqnx cruise departing her homeport on February 10, 2018, it is MLCE20180210. The individual cruise files at the sites above sometimes include data outside the study region.

Gridded data
The gridded datasets are the binned and averaged fCO2w, SST, and SSS observations on a (1˚ by 1˚ by mo) grid. The files include the auxiliary data obtained from remote sensing and interpolated data (OISST), data assimilation of remotely sensed 1375 winds (CCMP-2) ,and from the HYCOM model (SSSHYCOM). Calculated TA, pHT and ΩAr using procedures outlined above are provided in the file. The calculated fCO2w using the annual MLRs (Table 1) are provided as well. This gridded data set has spatial and temporal gaps as the ships did not transit through each pixel, and coverage is uneven. The number of observations differ for each grid cell and are listed in the gridded data set. For the auxiliary data the number of data points are fixed by the resolution of the data products except where part of the grid includes land which is masked. The column 1380 headers are provided in Table 3 and include units and descriptions when warranted.  (1990); KF from Perez and Fraga (1987) and Total Boron from Uppström (1974) ΩAr Aragonite saturation state calculated using CO2SYS with fCO2w_OBS and TA as input parameters and the same dissociation constants as used for pHT OISST ˚C Optimal interpolated sea surface temperature (Reynolds et al., 2007)

Mapped product 1395
The mapped product provides the data on a homogeneous (1˚ by 1˚ by mo) grid boxes utilizing the annual MLRs of fCO2w as a function of Lat, Lon, OISST, SSSHYCOM, and MLDHYCOM for the region from 15˚ N to 28˚ N and -62˚ to -88˚ (= 62˚ W to 88˚ W). The modeled and remotely sensed products OISST, SSSHYCOM, and MLDHYCOM and 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 1400 warranted. from Lueker et al. (2000); KSO4for Dickson (1990); KF from Perez and Fraga (1987) and Total Boron from Uppström (1974) ΩAr Aragonite saturation state calculated using CO2SYS with fCO2wMLR, TA, OISST, Deleted: The comparison between observed SST and OISST

1405
showed that the OISST was on average 0.25 ˚C lower than SST. The SSSHYCOM was 0.10 higher than SSS. The average difference between fCO2w_OBS and fCO2wMLR was 1.5 µatm with the fCO2w_OBS being higher. While the differences are within the standard deviation, they are possibly real due to known near-surface 1410 gradients in the region. No attempt was made to normalize the SST and SSS to OISST and SSSHYCOM. ¶

Monthly and Annual estimates for the Caribbean 2002-2018
Summary files of monthly and annual data and products covering the whole region from 15º N to 28˚ N and 88º W and 62˚ W are provided based on averaging or summing 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 are area weighted based on the area of the cell as provided in the mapped product (Table 4) according to: area cell/(total area/#of 1435 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 carbon (= 10 12 g C) per month or per year.  (2000); KSO4for Dickson (1990); KF from Perez and Fraga (1987) and Total Boron from Uppström (1974) ΩAr Aragonite saturation state calculated using CO2SYS with fCO2wMLR, TA, OISST, and SSSHYCOM as input parameters with same dissociation constants as used for pHT <u 2 > m 2 s -2 Second moment of the wind based on ¼˚ 6-h CCMP-2 product (Atlas et al., 2011) CO2_Flux mol m -2 mo -1 Monthly air-sea CO2 flux calculated according to Eqs. 5 and 6 CO2_FluxTotal Tg C mo -1 Total monthly air-sea CO2 flux calculated according to Eqs. 5 and 6 in Teragram carbon  (Wanninkhof et al., 2019a). The products cover the years 2002 through 2018 and will be updated annually.

Conclusions
The datasets from the cruise ships sailing the Caribbean Sea are a rich resource for studying the trends and patterns of 1490 inorganic carbon cycling and ocean acidification in the region. The scales of variability and data density are such that the (1˚ by 1˚ by mo) monthly gridding captures the magnitudes and trends of fCO2w and derived inorganic carbon products on 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 1495 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 near-shore fine scale models.

Team list 1500
This work was done on Royal Caribbean Cruise Lines (RCCL) ships who provided access, and personnel and infrastructure resources for the measurement campaign coordinated through the Rosenstiel School of Marine and Atmospheric Sciences      Lon is longitude; Lat is latitude, and SST is sea surface temperature. The MLR with the full number of independent 1780 parameters are used in the data products described the main text and the resulting data products are in Wanninkhof et al.  Tables A1 and A2. For the entire record, from 2002-2018, single regressions of fCO2w with position, SST, SSS, and MLD showed larger standard errors and coefficients of determination (r 2 ) as they do not capture the increase in fCO2w over time due to fCO2w 1790 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 multi-year changes in trends of fCO2w (Wanninkhof et al., 2019b). However, using the year as a variable in the regression provides a reasonable  Bates et al, 2014). Including SSS and MLD improves the multi-linear regression used for spatial gap filling to some extent (Tables A1-A2). The MLD is obtained from a numerical model, as no full 1815 observational record of mixed layer depths is available for the Caribbean. The residuals of fCO2w,MLR versus binned fCO2w for grid boxes with fCO2w observations show no spatial or seasonal pattern with any of the combinations of independent variables used in Eq. A1. ¶ 1820 ¶ 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).

1835
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 a MLR with either a subset or substitution of variables. Physical parameters were correlated with location (Lat, Lon) in the region. In particular, salinity showed broad correspondence with position. Therefore, substituting SSS for Lat and Lon provided similar magnitude and standard error in the coefficients of the independent variable, but a 10 % greater RMS in the estimated fCO2w: Eq. A4 which used 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), 1850 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 was the best option.