Large Decadal Changes in Air‐Sea CO2 Fluxes in the Caribbean Sea

Sixteen years of surface water CO2 data from autonomous systems on cruise ships sailing in the Caribbean Sea and Western North Atlantic show marked changes on interannual timescales. The measured changes in fugacity (≈partial pressure) of CO2 in surface water, fCO2w, are based on over a million observations. Seasonally the patterns are similar to other oligotrophic subtropical regions with an amplitude of fCO2w of ≈40 μatm with low wintertime values, causing the area to be a sink, and high summertime values making it a source of CO2 to the atmosphere. On annual scales there was negligible increase of fCO2w from 2002 to 2010 and a rapid increase from 2010 to 2018. Correspondingly, the trend of air‐sea CO2 flux from 2002 to 2010 was strongly negative (increasing uptake or sink) at −0.05 ± 0.01 (mol m year) year and positive (decreasing uptake) at 0.02 ± 0.02 (mol m year) year from 2010‐2018. For the whole period from 2002 to 2018, the fCO2w lagged the atmospheric CO2 increase by 24 %, causing an increase in CO2 uptake. The average flux into the ocean for the 16 years is −0.20 ± 0.16 mol m year with the uncertainty reflecting the standard deviation in annual means. The change in multiannual trend in fCO2w is modulated by several factors, notably changes in sea surface temperature and ocean mixed layer depth that, in turn, affected the physical and biological processes controlling fCO2w. Plain Language Summary Through a unique collaboration with Royal Caribbean Cruise Lines several cruise ships were outfitted with automated surface water carbon dioxide (CO2) measurement systems, providing weekly observations in the Caribbean Sea over the past 16 years. From over a million measurements, the increase in surface water CO2 in response to rising atmospheric levels was accurately monitored. The region is, on average, a carbon dioxide sink with large multiyear differences. For the first 8 years, the surface water levels do not change appreciably causing an increase in difference between the air and sea concentrations. This increased differential drives in increase in the air to sea flux of CO2. Carbon dioxide levels in surface water in the following 8 years actually increased faster than the atmosphere, thereby decreasing the differential and subsequent flux rate. The cause of changes in trends appears associated with changes in the ocean biogeochemistry, linked to sea surface temperature and mixed layer depth, particularly in the middle part of the record.


Introduction
There has been a large increase in measurements of fugacity of carbon dioxide in surface water (fCO 2w ) over the past two decades after realization that the quality of automated measurements was high enough to quantify global ocean CO 2 uptake (Takahashi et al., 2002). Recommendations to provide sustained systematic fCO 2w measurements in the world's oceans to estimate global air-sea fluxes were implemented (Bender et al., 2002). The Caribbean Sea was targeted for observations with passenger cruise ships, providing a unique platform for the unattended systems (Figure 1). Through a federal, industry, academic partnership between NOAA, Royal Caribbean Cruise Lines (RCCL) and the University of Miami, the cruise ship Explorer of the Seas (EoS) was outfitted with an automated surface water CO 2 measurement system (Pierrot et al., 2009) providing observations on alternating weekly transects in the Eastern and Western Caribbean. The biweekly occupations of the same locations were optimal to determine regional changes in surface water CO 2 levels in The Caribbean Sea is largely oligotrophic with low nutrient concentrations in surface water. However, several local environmental factors can impact its nutrient and carbon dynamics. Episodic riverine influence from the Orinoco River increase biological productivity in the southwestern part of the region, based on observed lower salinity and remotely sensed color (Lopez et al., 2013). Runoff and other coastal influences from the islands have an impact on the local CO 2 dynamics and flux estimates (Melendez Oyola et al., 2018). Hurricanes leave an imprint on fCO 2w through large effluxes during the extreme wind events and subsequent decreases in fCO 2w caused by lower SST and enhancement of biological productivity due to enhanced nutrient exchange from below the mixed layer (Huang & Imberger, 2010;Wanninkhof et al., 2007). The riverine, island, and hurricane influences are significant at local and short time scales but do not show a large impact on seasonality or the regional scale patterns and multidecadal trends of air-sea CO 2 fluxes for the whole region.
The canonical view that the surface inorganic carbon dynamics in the region are largely driven by temperature and processes impacted by temperature holds largely true. However, like the ocean time series stations in the oligotrophic gyres, near Hawaii (HOT) and Bermuda (BATS), there are multiyear trends that cause the surface water CO 2 increases to deviate from the expected rate of increase (Bates et al., 2014;Dore et al., 2009). The divergent multiyear trends observed in the Caribbean are significantly greater than observed at BATS and HOT for the time period in the references listed. This is unexpected because the shallow MLD, relatively low productivity, and strong stratification should cause the fCO 2w increases to be largely controlled by atmospheric CO 2 increases, leading to fCO 2w increases of about 2 μatm/year.
The multiyear changes impact the rate of ocean acidification (OA), which is of concern in the region where coral reef cover has declined by over 50% in the last 5 decades due to multiple human stressors (Jackson et al., 2014). The current observation-based estimates of OA for the region rely on empirical correlations between inorganic carbon chemistry versus SST and SSS from the earlier years of this observation program from which an ocean acidification product suite was developed (Gledhill et al., 2008). In this model, annual increases in atmospheric fCO 2 are directly coupled to surface ocean acidification rates with subannual variability governed primarily by SST. The results of the current effort indicate that this approach can be finetuned with the 16-year record of observations showing the multiannual anomalies in trends in fCO 2w .
Seasonal changes in fCO 2w in the Caribbean were first detailed by Olsen et al. (2004) based on the first year (2002) of cruise data from the EoS. The fCO 2w levels were closely related to SST and, thus, had lower values in winter and higher in the summer. This annual cycle of fCO 2w is such that in the wintertime the region is a CO 2 sink and during the summer is a CO 2 source. Olsen et al. (2004) created a simple but robust algorithm relating fCO 2w to SST and location based on the 2002 data: fCO 2w ¼ 10:18×SST þ 0:5249×Lat−0:2921×Lon þ 52:19; r 2 ¼ 0: 87 (1) where SST is in°C; Lat is latitude in fractional degrees (North); and Lon is longitude as negative degrees (East). The root-mean-square (RMS) uncertainty of calculated fCO 2w was 5.7 μatm.
Equation (1) reflects the strong relationship of fCO 2w with SST. It also indicates that at constant SST, fCO 2w will be on average 5 μatm higher in the southern part (15°N) than in the northern part (25°N) and 8 μatm higher in the western part (88°W = −88°E) compared to the eastern boundary of the region (62°W). These spatial gradients are relatively small compared to the ≈40-μatm seasonal amplitude associated with a 4°C change in SST from winter to summer in the region. Park and Wanninkhof (2012) analyzed 8 years of fCO 2w data from 2002 to 2009, in the northeastern part of the region, north of Puerto Rico, predominantly from the EoS. During this time period the fCO 2w increased ≈9 μatm, compared to expectations that it should have changed by ≈18 μatm in response to atmospheric CO 2 increases. The stagnation was attributed primarily to wintertime fCO 2w increasing very slowly over this time period. This was caused by anomalously high wintertime SST in 2002-2003and lower wintertime SST in 2008. Winter MLD also increased toward the end of the time period. As a result, the air-water fugacity of CO 2 difference, ΔfCO 2 , became more negative in wintertime causing an increase of the CO 2 sink strength in the region for the 2000-2009 period.
Here we build on the previous work, and extend the time series and analyses from 2002 to the spring of 2018. Data products such as wind speed needed to determine the fluxes along with remotely sensed SST, modeled SSS, Normalized Fluorescence Line Height (NFLH), and modeled MLD are used to interpret to results. The sources of these data are detailed in the data section. Instrumental procedures, time series analyses, mapping procedures, and flux calculations are described in the methods section. The seasonal cycles and multiyear trends of the variables of interest are described, and air-sea CO 2 fluxes over the region are quantified. The marked changes from a stagnant annual average fCO 2w over the first part of the record from 2002 to 2010 to a larger than atmospheric rate of increase in fCO 2w from 2010 to 2018 are quantified. The causes thereof are discussed and deconvolved into temperature, gas exchange, and a residual attributed to biological and mixing (B&M) components.

Data and Data Products
In this paper, three different datasets are used and compared. The following nomenclature is used: observational datasets from the cruise ships using the time and position of measurements are referred to as observations; these data are binned and averaged on a 1°× 1°space scale and monthly timescale (1°× 1°× mo) and are called gridded data (product). The gridded data are useful to determine biases between observed SST and SSS, and remotely sensed SST and modeled SSS. They are also used to determine if gridding and extrapolation procedures impact the interpretation. For the whole region spanning 16°N-28°N and 88°W-62°W, the fCO 2w data are extrapolated using multilinear regressions (MLRs) with position, SST, SSS, and MLD as described below. They are used in conjunction with monthly remotely sensed SST, wind, and NFLH, and modeled SSS on the same grid and are referred to as the (gridded) mapped product. This product is the basis of the trend estimates and the explanation of trends. Using the gridded mapped product, the trends in fCO 2w and flux for the whole time period and area can be determined without artifacts caused by data gaps.
The field observations were obtained from the three cruise ships listed above totaling over 1 million data points over the 16-year period from March 2002 to 2018. They are posted on the website: www.aoml.noaa. gov/ocd/ocdweb/occ.html per cruise and ship, along with pertinent metadata. The data are archived in the NCEI Ocean Carbon Data System (www.nodc.noaa.gov/ocads/). Data are also collated annually in the Surface Ocean Carbon ATlas, SOCAT (Bakker et al., 2016) and LDEO (Takahashi et al., 2017) data products.
Here we use the original data as served from the AOML website. The fCO 2 data are the same for all sites but the quality control procedures and flags differ slightly. Only data with quality control flag 2 (= good data) are used, which are over 98 % of the observations.
For the determination of the fugacity of CO 2 in air, fCO 2a the average monthly atmospheric CO 2 mole fraction (XCO 2a ) product based on weekly measurements at the stations at Key Biscayne, FL, USA (KEY; 26°N, 80°W), and Ragged Point, Barbados (RPB; 13°N, 59°W) are used. They are obtained from the flask sampling network of the Global Monitoring Division of the Earth System Research Laboratory (GMD/ERSL/NOAA; Dlugokencky et al., 2017). These values are quality controlled at GMD and fit to the annual seasonal cycle (see www.esrl.noaa.gov/gmd/ccgg/mbl/data.php). The monthly product is considered representative for the whole region. The difference in monthly mean between the KEY and RPB products from 2002-2018 is 0.6 ± 0.6 ppm (n = 189) with KEY values being on average higher. Support data used for analyses are on the same 1°× 1°× mo grid as the mapped fCO 2w . Wind speeds are from the updated cross-calibrated multiplatform wind product 2 (Atlas et al., 2011). The 1°× 1°× mo mean neutral wind at 10-m height <u 10 >, and its second moment <u 10 2 > is determined from the ¼ degree, 6-hourly product from Remote Sensing Systems (www.remss.com). The SST product is the Optimum Interpolated SST, OISST (Reynolds et al., 2007). It uses data from ships, buoys, and satellites to generate the analyzed fields. For the OISST, the reference SST is from buoys, and the SST from ships in this product are adjusted to the buoy data by subtracting 0.14°C. The OISST is on average 0.25 ± 0.40°C (n = 8,191) lower than SST measurements from the cruise ships when comparing the 1°-gridded data product for the particular month with OISST product. The SST probes on the cruise ships have accuracies better than 0.01°C, and while the offset is within the standard deviation, it is assumed the difference is real and attributed to the near-surface cool skin (Soloviev & Schlüssel, 1996).
No MLD determinations were made on the cruises and limited observational estimates from other sources are available. The profiling float program, Argo, produces a global MLD product in near real time (Holte et al., 2017). However, Argo only has limited deployments in the Caribbean and no MLD product is available for the region. The MLDs used here are from a numerical model (HYCOM; https://hycom.org/regional). The MLDs are based on a density contrast of 0.03 between surface and subsurface (de Boyer Montégut et al., 2004).
The SSS data used were also extracted from this model, SSS hycom . The SSS hycom is on average 0.1 ± 0.28 (n = 8,191) greater than the ship-based 1°× 1°× mo data product. The difference does not impact the results or interpretation as the SSS has a weak influence on the fCO 2w calculated from the MLRs.
NFLH, which is used as a proxy for biological productivity, was retrieved from the moderate resolution imaging spectrophotometer, MODIS. MODIS flies on the satellites AQUA and TERRA. NFLH is a relative measure of water-leaving radiance associated with chlorophyll fluorescence. The level 3 data used are at 1°× 1°× mo resolution and are obtained from http://oceancolor.gsfc.nasa.gov/cgi/l3 website.

Instrumentation
The fCO 2w levels are measured by infrared (IR) analysis of headspace gas equilibrated with water in an equilibration chamber. The gas in the headspace is partially dried (≈>75%) before measurements. The procedures and instrumentation are based on a community design described in Pierrot et al. (2009). The instruments are manufactured by General Oceanics INC. and have performed to high accuracy specifications. The instruments have shown agreement to within 1 μatm with other state-of-art systems in intercomparison studies (Y. Nojiri et al., personal communication, February 2009). The IR sensor (LI-COR 6262) is calibrated every 4.5 hr against four standard gasses supplied by the GMD/ESRL/NOAA, traceable to the WMO CO 2 mole fraction scale. The CO 2 concentrations in the standards span the range of values encountered along the ship tracks. 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, have remained the same.
Pressure in the headspace and temperature of the water in the equilibrator and near the seawater intake, assumed equivalent to SST, are measured at high accuracy. Equilibrator temperature is obtained with a Hart RS thermistor accurate to ±0.002°C. The pressure is measured by a Setra model 270 barometer which is accurate to 0.1 mbar (hPa). SST and salinity are determined by a Seabird thermosalinograph accurate to 0.01°C and 0.01 based on factory specifications. As conductivity values used to calculate salinity change due to fouling, the in situ observations probably are less accurate than specification but the uncertainty in salinity has negligible impact on the calculation of fCO 2w .
The water intake for the underway pCO 2 instrument and thermosalinograph is near the bow thruster cavity on the ships at a depth of 1.5 m. The shallow depth of the intake and operation of the bow thruster can cause bubble entrainment near ports and in heavy seas. Data obtained under these conditions are flagged and not used in the analyses. Cruises follow a fixed route with occasional deviations due to weather, such as averting heavy seas during and immediately after hurricanes. Typical cruise speeds are 22 knots that with sampling every 2.5 min, yields a sample approximately every 1.7 km except for 35 min every 4.5 hr when four calibration gases and a CO 2 -free reference gas are analyzed followed by four atmospheric CO 2 measurements. This causes a distance of 24 km (≈1/4°) without fCO 2w measurements. The air for the atmospheric measurements is obtained from an intake mounted near the bow of the ship and piped to the infrared analyzer. Air measurements do not cover the entire record. They were started on the EoS in 2008 and on the Equinox in May 2015. Due to challenges installing the airline that needs to be strung up eight decks and forward in the ship's interior there is no air intake line on the Allure of the Seas. As there are extensive periods when no air measurements are available, the values used here were determined from two locations in the MBL CO 2 product as described in the data section.

Calculation of fCO 2
The calibrated values of XCO 2 from the IR analyzer are accurate to within 0.2 ppm based on the certified values of the calibration gases that are traceable to the WMO scale, and uncertainties in their interpolation. The XCO 2 values are converted to fCO 2eq (μatm) values following Pierrot et al. (2009): where eq refers to equilibrator conditions. P eq is the pressure in equilibrator headspace and the pH 2 O is the water vapor pressure calculated according to equation (10) in Weiss and Price (1980). The T eq is the temperature of water in the equilibrator. The function f(T,P) is the fugacity correction (Weiss, 1974). The fugacity is the partial pressure of CO 2 corrected for nonideality of the gas. The virial coefficients used in the correction as well as the appropriateness of the correction are under debate. The magnitude of the correction (f(T,P) ≈0.997) for converting from pCO 2 to fCO 2 is about ≈1 μatm at the fCO 2 ranges encountered, so the uncertainties in the correction have no effect on the conclusions. The fCO 2eq is corrected to surface water values using the intake temperature (SST) according to the empirical relationship that Takahashi et al. (1993) developed for North Atlantic surface waters: on average SST − T eq = 0.12 ± 0.28°C (n = 667,865). The SST is higher than T eq because of the warm SST in 10.1029/2019JC015366 Journal of Geophysical Research: Oceans the region and lower air temperatures in the air-conditioned ship, in contrast to most ships where water warms once inside.
The fCO 2a (μatm) are calculated from XCO 2a according to: here P is the ambient pressure at sea level, pH 2 O is determined at SST, and f (SST, P) is the fugacity correction.

Calculation of Air-Sea CO 2 Flux
For the determination of the air-sea CO 2 flux, F CO2 (mol m −2 year −1 ), a bulk formulation is applied to the data from the gridded mapped product: where ΔfCO 2 is (fCO 2w − fCO 2a ), s is the seawater CO 2 solubility (Weiss & Price, 1980), and k is the gas transfer velocity parameterized as a function of wind speed (Wanninkhof, 2014): where <u 10 2 > is the monthly second moment of the 6-hourly wind speeds reported in cross-calibrated multiplatform wind product 2, accounting for the impact of variability of the wind speed on k. The Sc is the Schmidt number of CO 2 in seawater determined as a function of temperature from Wanninkhof (2014).

Gridded Mapped fCO 2w Product Using Multilinear Regressions
Monthly fCO 2w fields for the region were created in the following fashion. Annual MLRs were determined from the fCO 2w data that were binned and averaged on a 1°× 1°× mo grid similar to the approach of Olsen et al. (2004). To take advantage of improved regional products of salinity and mixed layer depth these variables were included in the MLRs. For SST both a linear and logarithmic dependence were checked but little difference was observed in the fit, and a linear SST dependence was used. The functional form for the MLR fit is as follows: MLRs were determined for each year, so no assumptions about annual fCO 2w trends had to be made. The coefficients and their standard errors for each annual MLR are provided in Table 1. SST and location are the strongest predictors. Multicollinearity was checked for the MLRs created for each year. To determine if two or more of the independent parameters (predictors) were correlated, variance inflation factors (VIFS) were calculated and presented in Table 1. The VIFs were well below 5, considered the threshold for multicollinearity. The exception is for 2008 when the VIF for Latitude was 5.7. The residuals of mapped fCO 2w,mapped versus the gridded fCO 2w data product (1°× 1°× mo) for pixels with fCO 2w observations show no spatial or seasonal patterns.

Determination of controls of fCO 2w
For investigation of processes controlling fCO 2w , the monthly changes in fCO 2w are decomposed into a temperature term, an air-sea CO 2 flux term, and a residual that is controlled by biological processes and mixing: where δ is the monthly time derivative such that δfCO 2w (Total) is the monthly change in fCO 2w . The δfCO 2w (SST) is the monthly change in fCO 2w due to temperature change: This is referred to as the thermodynamic effect, where the subscripts i and i + 1 are the values for consecutive months. δfCO 2w (F CO2 ) is the monthly change in fCO 2w due to the air-sea CO 2 flux. This is determined in the  Journal of Geophysical Research: Oceans calculated. Then DIC (i+1) = DIC (i) + F CO2 /MLD. Accordingly, fCO 2w(i+1) is calculated with CO2SYS using DIC (i+1) , TAlk (i+1) , SST (i+1) , and SSS (i+1) . The change in fCO 2w due to F CO2 is as follows: The δfCO 2w (B&M) is determined as the residual: The δfCO 2w (B&M) is the monthly change caused by other processes that impact fCO 2w , including biological drawdown due to biological production, and remineralization and respiration in the mixed layer, the sum which is referred to as net community production. The term also includes effects of mixing and entrainment of waters with different DIC concentrations (and fCO 2w ) from below the mixed layer and from lateral transport. The δfCO 2w (B&M) also implicitly contains processes that change alkalinity in the mixed layer. Alkalinity increases will lower the fCO 2w . While some of the alkalinity effects could, in principle, be determined from changes in salinity due to the strong relationship between salinity and alkalinity, month to month salinity changes are small such that the alkalinity effects on fCO 2w are not studied separately but rather included in this residual term.
The fCO 2w is not a state variable and the deconvolution described in (8) is not strictly additive. For instance, the δfCO 2w (SST) does not change the state variable DIC, while the δfCO 2w (F CO2 ) and δfCO 2w (B&M) terms do impact DIC. However, on monthly basis where the changes are relatively small, a linear approximation as used here is justified. More formal deconvolutions to tease out the processes influencing fCO 2w have been performed (Metzl et al., 2010;Takahashi et al., 1993;Thomas et al., 2008) using the state variables DIC and TAlk, and the derivatives of fCO 2w with DIC, TAlk, temperature, and salinity. In our case where DIC and TAlk were not measured but rather estimated using salinity and fCO 2w , the uncertainty in the estimated of DIC and TAlk changes, along with use of a modeled mixed layer, warrant the simpler approach of using a residual B&M term that includes several processes.
The uncertainties in the deconvolution cannot be easily quantified. They include uncertainty in the flux of ≈20%-30% (Wanninkhof, 2014), uncertainty in MLD estimate, spatial and temporal variability in the grid cell, and nonlinearity of processes changing fCO 2w . Therefore, no uncertainty estimate is provided due to challenges in properly determining and propagating the errors in each term. Biases in δfCO 2w (Total), δfCO 2w (SST), and δfCO 2w (F CO2 ) will be reflected in the δfCO 2w (B&M) term. Accordingly, this is interpreted with caution.

Deseasonalization Using Harmonic Fits
The seasonal cycle is the primary mode of variability in the region and it exhibits a regular pattern. The fCO 2w and SST can be well represented with a harmonic fit. In turn, this fit can be subtracted from the monthly values to more clearly discern the multiannual trends. The seasonal cycles of fCO 2w and SST of the gridded data product (1°× 1°× mo) are fit to a harmonic function as detailed in Gruber et al. (1998) and Park and Wanninkhof (2012). For fCO 2w , fCO 2w;fit ¼ fCO 2w;mean þ a sin 2 π mo=12 ð Þþb cos 2 π mo=12 ð Þ where fCO 2w,mean is the average of all the fCO 2w values for all grid points over the 16-year period (373.6 μatm), mo is the month (1-12), and a and b are fitting parameters optimized through a cost function. Using a = −20.4 and b = −10.6 results in an average root-mean-square error, RMSE between the data and the fit of 1.3 μatm. The same procedure is used for SST, with an average SST of 26.95°C and a RMSE of 0.005°C. The monthly average values of fCO 2w and SST determined with harmonic fits covering the 16-year period are subtracted from the monthly fCO 2w and SST values for the region to obtain a seasonality detrended, or deseasonalized, product. Using the deseasonalized residuals of fCO 2w and SST improve the detection of longerterm trends and increases the confidence level of the trends. The detrending procedures do not bias the results and the same multiannual patterns are discerned with the original nondeseasonalized observations and mapped products.

Results and Discussion
The fCO 2w levels are influenced by processes operating over a range of temporal and spatial scales, from hours to decades, and from local to regional. The patterns of fCO 2w and air-sea CO 2 fluxes in time and space are discussed for the study region shown in Figure 1 nominally covering 16°N -28°N and 88°W-62°W. Unless noted, the mapped product is used which is on a (1°× 1°× mo) grid from the annually derived MLRs (Table 1).
The lack of observations in the Southeastern Gulf of Mexico north of 25°N can cause the interpolation with MLRs to yield erroneous results in this area as it is affected by other influences than typical for the Caribbean Sea. For example, shelf specific processes on the West Florida shelf, and Mississippi outflow can impact the fCO 2w in the Southeastern Gulf of Mexico (Robbins et al., 2018;Xue et al., 2016;Yang et al., 2015). The area is included in the calculation of the total air-sea CO 2 flux, but the results might not be fully representative.
The complete observational dataset of fCO 2w and SST from March 2002 to 2018 are plotted against time in Figure 2. The power spectrum of fCO 2w shows a singular strong peak at annual time scales. This corresponds with the power spectrum of temperature. This differs from the analysis of mooring based CO 2 data near Bermuda by Bates et al. (2001) that showed a strong daily peak and a broad peak centered on a week based on analysis of monthly data. In both the Bermuda case and this study, the fCO 2w and SST spectra match closely, reflecting the role of SST as a direct and indirect driver of fCO 2w . The lack of a peak in SST and fCO 2w on daily scale for the Caribbean data set is attributed to the signal not being resolved at regional scale of analyses and the ships covering the region at a set schedule. That is, the cruise ships cover the same locations at the same time of day for each cruise itinerary, and the ships spend more time at sea during nighttime such that the diurnal SST and fCO 2w cycles at are not well resolved.

Spatial Differences
Spatial differences are small compared to seasonal changes but they are observable. This can also be inferred from the MLRs of fCO 2w (Table 1) that show strong dependencies with position. Based on the MLRs, the average fCO 2w values decrease towards the north and toward the east. As temperature is the strongest predictor in the MLRs and there is a correlation between location and temperature, the spatial trends are primarily related to SST. Maps of fCO 2w , SST, and SSS are shown in Figures 3a-3f using February and August 2017 as representative examples. The largest spatial differences in fCO 2w are apparent in the wintertime (Figure 3a) with high values in the southwest and lower values to the northeast. In summertime (Figure 3b), the highest fCO 2w values are in the northwestern part of the region. The relationship of fCO 2w with SST is apparent in the patterns. The SSS shows spatial and seasonal variations (Figures 3e and  3f), but the correspondence with fCO 2w is not clearly observed, except in the southeastern part of the region where lower SSS waters typically also have lower fCO 2w in February (Figures 3a and 3e).
The differences in the annual cycle for four representative locations shown in Figure 1 are highlighted in Figure 4 where the gridded fCO 2w , SST, SSS, and MLD products are plotted versus month in 2017. Monthly mean gridded data products from the Florida Current, South of the Florida Keys (pixel centered at 24.5°N, 81.5°W); the Outer Bahamas Banks (Turks and Caicos, 21.5°N, 70.5°W); the Lesser Antilles (18.5°N, 64.5°W); and the western Caribbean (Cayman Islands, 19.5°N, 81.5°W) are shown. The spatial differences in fCO 2w and SST are smaller than the extremes in the seasonal cycle (Figures 4a and 4b). Local anomalies are apparent, particularly in SSS (Figure 4c) where a strong minimum is observed in the Lesser Antilles region in the summer and fall, attributed to a combination of Orinoco river outflow and regional precipitation anomalies (Ibánhez et al., 2017). Annual SSS anomalies are seen in other locations as well, such as the minimum observed in summer at the Florida current location. The fCO 2w can be directly affected by Journal of Geophysical Research: Oceans dilution with rainfall (Turk et al., 2010) or, in case of riverine runoff, indirectly as well though enhanced biological productivity from excess nutrients or alkalinity carried by rivers. The low fCO 2w levels near the Lesser Antilles appear largely caused by dilution as the observed decrease in fCO 2w corresponds to the thermodynamic effect that is 17.5 μatm per unit salinity change (Takahashi et al., 1993). For the low salinities in the Florida current, there is no corresponding reduction in fCO 2w , likely because of the increasing SST. The model derived MLD for four locations are given in Figure 4d. All show a seasonal cycle with late summer MLD minima. The Florida current site shows MLD as shallow as 15 m during the summer. The monthly MLD product shows significantly deeper MLD in the western Caribbean compared to the other locations throughout much of 2017.
Local deviations from basin-wide patterns are apparent for other years as well. For example, a maximum in fCO 2w of 445 μatm that was observed in the Florida current in July 2016 (not shown) was ≈20 μatm higher than values observed in the adjacent months of 2016 and the fCO 2w values observed in July 2017 (Figure 4a). The corresponding SST in July 2016, of 30.3°C, is only ≈0.4°C higher than July 2017. Thus, the fCO 2w anomaly in July 2016 is about twofold larger than expected from a thermodynamic temperature effect of 4.23%°C −1 . This high fCO 2w value is attributed to entrainment of West Florida Shelf water and is an example of local anomalies that are present throughout the record. However, while local monthly anomalies in SST, SSS, and fCO 2w are apparent for different locations, and for different years, the Journal of Geophysical Research: Oceans region as a whole behaves similarly and, as elaborated upon below, can be viewed as a single entity to assess seasonal to interannual variability and trends.

Subannual (Seasonal) Patterns
Different factors influence the seasonal fCO 2w cycle but seasonal changes in SST are the dominant influence on fCO 2w . The thermodynamic dependence of fCO 2w to SST is 0.0423°C −1 (4.23% per degree centigrade) due to decreasing solubility and shifting of carbonate equilibrium constants favoring fCO 2w increase. The SST changes by ≈4°C over the annual cycle in the Caribbean and thus has the potential to change fCO 2w by ≈60 μatm.
Seasonal changes in MLD can affect fCO 2w in differing and opposing ways. Deepening mixed layers can entrain colder water that will thermodynamically decrease the fCO 2w . However, the waters can contain higher DIC, which will increase fCO 2w . Entrainment and mixing can also supply nutrients, and the resulting primary production can decrease fCO 2w . Variability of fCO 2w due to mixed layer dynamics and nutrient supply are not well understood (Fawcett et al., 2014). The gradients in inorganic carbon and nutrients below the mixed layer up to the maximum wintertime MLD are small based on GO-SHIP line A22 cruise data (cchdo. The average seasonal cycle for MLD, NFLH, SST, and fCO 2w based on the monthly mapped products for the entire region are shown in Figure 5 with the sd of the monthly values over the 16 years depicted as error bars. Multiannual trends have not been subtracted. The MLD (Figure 5a) shows shallow depths of 21 ± 6 m from May to September and rapid deepening thereafter reaching a maximum average depth of 54 ± 14 m in January after which a shallowing occurs in early spring. The NFLH (Figure 5b) shows a minimum from March to May. The decrease corresponds roughly with the mixed layer shallowing and possibly caused by limited diffusion of nutrients from below during this time.
The seasonal signal in SST (Figure 5c) and fCO 2w (Figure 5d) can be well represented using a harmonic function (equation (12)). The seasonal range in fCO 2w is 40 μatm and 4°C in SST with minimum values in February and maximum values in August/September. There is no distinct annual cycle in salinity (not shown).
The close correspondence in phase and shape of the seasonal cycle of SST and fCO 2w indicates a strong correlation between them. Figure 6 shows the monthly averaged fCO 2w plotted against SST following Lefèvre  (2002). This diagram illustrates the factors influencing the seasonal cycle in fCO 2w . The average change in fCO 2w with respect to temperature based on a linear regression of monthly fCO 2w versus SST is 9.5 μatm°C −1 (2.5%°C −1 ). This is about 60% of the thermodynamic trend, which is shown as a dashed line in Figure 6 indicating that there are counteracting factors impacting the seasonal cycle of fCO 2w .
The progression of SST and fCO 2w through the year can be delineated by periods of warming from March to August, cooling from September to February, and transitions in February-March, and August-September, respectively (Figure 5c). The seasonal trends in fCO 2w are caused by interplay of temperature changes, gas transfer, mixing, and biologically mediated processes. The partitioning of the processes (equation (8)) for the seasonal cycle using gridded product are shown in Figure 7. The deconvolution of fCO 2w for the full 16-year record is presented in the attribution section below.
As depicted in Figures 6 and 7, during periods of warming (March-September) the fCO 2w increases by ≈8-10 μatm/mo or 2.5%-3.1%°C −1 . The thermodynamic effect, which is predominant in the annual cycle, increases fCO 2w . The net effect of B&M is a drawdown of fCO 2w during this warming period. The effect of air-sea gas exchange is of the same magnitude as the B&M impacts (≈1 to 3 μatm/mo) but changes sign during the warming period as it depends on the sign of the ΔfCO 2 gradient. It contributes to increasing fCO 2w when the area is a net sink from March to May and to decreasing fCO 2w from June to September. During the cooling period from September to February, fCO 2w decreases by ≈8-10 μatm/mo, 2.5%-3%°C −1 . Again, the thermodynamic effect dominates, with counteracting effects of B&M. The effect of air-sea fluxes reverses from a CO 2 sink in October to a small CO 2 source for the remainder of the cooling period. The B&M effect increases fCO 2w during the cooling period except in February. In the transition month from cooling to warming (February-March), the flux into the ocean and B&M are of equal but opposite magnitude. In September, when the transition from warming to cooling occurs, gas evasion counteracts the small net warming with a small contribution of B&M.
The B&M term is calculated as a residual (equation (11)) between the observed monthly fCO 2w change, the thermodynamic effect, and the effect of gas transfer. The observed monthly change in fCO 2w and the thermodynamic effect are well constrained but biases in the gas flux can be significant at 20%-30%. Generally, the air-sea CO 2 flux term is a smaller contributor. The δfCO 2w (B&M) is negative from February to August and positive from September to January. This does not match the NFLH annual progression that shows a minimum in spring and a maximum in December ( Figure 5b). However, δfCO 2w (B&M) does show broad correspondence with MLD ( Figure 5a). This suggests that during periods of shallow mixed layers the biological production is lower but still draws down the fCO 2w while during times of mixed layer deepening nutrients are entrained, enhancing productivity and elevating the NFLH. However, due to the deeper MLD the effect on fCO 2w is smaller. Also, with deeper MLD remineralized carbon in the form of DIC are brought to the surface that causes the B&M effect to increase fCO 2w .

Temporal Trend
There are clear trends and patterns on multiannual scales in the region. The fCO 2w observations versus time show a well-defined seasonal cycle and long-term trends (Figure 2). For the entire observational record from March 2002 to February 2018, the increase in fCO 2w is 1.30 ± 0.003 μatm/year but there is a large change in the middle part of the 16-year record, such that it can represented as two linear segments (Figure 2). From March 2002 to February 2010, the observed fCO 2w decreases by −1.37 ± 0.017 μatm/year with a reversal in trend to 3.69 ± 0.011 μatm/year from March 2010 to February 2018. The uncertainties are the standard errors in the slope and are small due to the large number of data points. As data coverage is sparser from 2007 to 2011 than at the beginning and end of record, the exact year of turnaround in trend is not well defined in the observations. A similar trend in the mapped product for the whole region is seen using the deseasonalized monthly residuals, albeit with lower, and only positive trends (Figure 8a). It shows an increase of 1.66 ± 0.09 μatm/year for the entire record. The trend is 0.60 ± 0.21 μatm/year from 2002 to 2010 and 2.55 ± 0.28 μatm/year from 2010 to 2018. While the trends for the observations and mapped values differ, the pattern of little change or deceases in fCO 2w for the first 8 years followed by a strong positive trend is the same. The differences between the observations and the gridded mapped product are attributed to a combination of stronger trends in fCO 2w and SST in the central region along the tracks where the observations are taken and the use of MLRs in the mapped product that causes features to be smoothed.
To determine the timeframe when the trend changes, 3-year increments are compared, using the mapped product (Table 2). This follows the approach of Fay and McKinley (2013) and allows trends for different increments and lengths to be compared. There is a broad minimum in fCO 2w trends centered on the 2008-2011 time frame that cannot be attributed to a single anomalous year. The residuals of the detrended/deseasonalized fCO 2w data (Figure 8a  The trends can be put in context of the rising atmospheric CO 2 levels. At steady state, the fCO 2w would keep up with the atmospheric CO 2 rise of, on average, 2.21 ppm/year from 2002 to 2018 based on the data from the KEY and RPB stations. This rise in mole fraction of CO 2 in the marine boundary layer, XCO 2a translates to 2.13 μatm/year at 100% humidity and 27°C for fCO 2w . Thus, for the first part of the record, the fCO 2w shows an appreciable lag while in the second part of the record the increase is greater than expected from invasion of CO 2 from the atmosphere. Over the entire 16 years, the fCO 2w has increased but it is appreciably less than expected if the Caribbean Sea had kept up with the atmospheric increase leading to a deficit of 13 μatm (24%) compared to the expected trend by the end of the record. The short-dashed line in Figure 8a is the expected trend in the residual if the fCO 2w followed atmospheric CO 2 increase and SST trends (Figure 8b).

Attribution of Trends in fCO 2w
For the attribution of multiyear trends, we first describe the statistical correlations between the variables in the gridded mapped product. This is followed by the description of how SST and MLD impact the trends. Finally, the effect of temperature, gas exchange, biology, and mixing (B&M) on multiyear trends fCO 2w are described. No single overriding cause for the large change in trend in fCO 2w is apparent but rather it is a combination of factors, including multiyear seasonal anomalies in the middle part of the record that cause the large changes in fCO 2w trends.
The correlation analysis of the gridded data provides a summary of the correlation of fCO 2w with the environmental parameters as well as correlations between the different parameters ( Figure 9). For the whole mapped product, the strongest dependency of fCO 2w is with temperature (positive) but there are also statistically significant correlations (at >90% confidence level) with wind (negative), position (negative), MLD (negative), and NFLH (positive). The fCO 2w shows insignificant statistical dependency on SSS. The sign of the correlations indicates that fCO 2w increases in response to the SST increase and decrease in response to increasing wind. The fCO 2w is negatively related with mixed layer depth indicating that deeper mixed layers lead to lower fCO 2w . The NFLH shows a positive correlation with fCO 2w , similar to what was shown in the subannual analysis, contrary to the expectation that greater fluorescence would indicate higher productivity and thus lower fCO 2w . The unexpected signs of the correlations of fCO 2w with MLD and NFLH are, in part, due to parameters being cross-correlated such that the statistical relations do not always imply causality. For example, SST and MLD are negatively correlated and the strong dependency of fCO 2w on SST will decrease the effect of entrainment of DIC due to mixed layer deepening. The NFLH and mixed layer depth are positively correlated, while fCO 2w and MLD are negatively correlated confounding simple causal relationships of the fCO 2w with individual parameters. The positive sign and weak correlations suggest that NFLH is not a strong predictor of the effect of productivity on fCO 2w in the region.   Of the parameters investigated, SST has the greatest influence on fCO 2w . Based on the relationship of ∂fCO 2w (∂SST) −1 of 0.042°C −1 , which under average conditions in the Caribbean Sea causes fCO 2w to change by ≈16 μatm per°C, the trend in observed fCO 2w due to surface ocean cooling/warming would be −2.9 μatm/year from 2002 to 2010 and 3.1 μatm year from 2010 to 2018 for the observed SST, and −1.0 and 1.2 μatm/year, using the mapped OISST product. Figure 8a includes trend lines depicting the changes in fCO 2w expected due to changes in atmospheric CO 2 and SST (thermodynamic effect) for the OISST mapped product in these periods.
The fCO 2w trends are affected by seasonal changes in MLD as shown in Park and Wanninkhof (2012). They showed that the wintertime fCO 2w decreased in the northeastern part of the study region from 2002 to 2009. This was attributed to longer periods with deep mixed layer depth (>30 m) providing more volume for CO 2 uptake and thus small changes in fCO 2w due to gas exchange. The summertime fCO 2w increased over the same time period. However, on annual scale the fCO 2w remained invariant thereby increasing the uptake significantly. Following this study, a more detailed investigation of the seasonal controls on fCO 2w trends Table 2 The 3 Table 3.
The wintertime and summertime data are also investigated to determine differences in trends between gridded data products and mapped products for SST and MLD. For the 1°× 1°gridded data, the trends in in situ measured SST are compared with the OISST, denoted in Table 3 as Obs and Mapped@Obs, respectively. This indicates if there are biases in seasonal trends using the MLRs to estimate fCO 2w , compared to measured fCO 2w .
Over the entire period 2002-2018, the trends during summer and winter for the mapped product and gridded data are consistent but with some differences (Table 3). Wintertime SST shows a slight warming for the mapped product and gridded data of ≈0.02 to 0.05°C/year. The summertime OISST for the whole region show an increase as well of 0.03 ± 0.01°C/year while observations along the cruise tracks are showing a nonsignificant cooling with a change of −0.02 ± 0.03°C/year. The OISST at the sampling locations show no trend but overall the observations, gridded data, and mapped products agree to within their uncertainties. The MLD shows an appreciable increase of ≈0.3 m/year over the 16-year time period with general correspondence between summer and winter.
When breaking down the summer/winter trends into two time periods, differences are more apparent between the gridded data and the gridded mapped products. From 2002 to 2010, the gridded SST data along the cruise tracks show a decrease of ≈−0.16°C/year in wintertime and no significant change in summer. The whole region, however, does not show a strong trend in OISST in winter. The differences are in the western region, which was not well sampled by the ships during this time period. The area shows a slight warming in the OISST product. The MLDs show a decrease of 0.3 m/year in summer and a large increase of 1 m/year in wintertime at the observation. Summertime shallowing is the same but there is no deepening in the mapped product. This wintertime MLD increase at the observational grid points correlates with the SST decrease during this period and similarly the lack of SST and MLD change in the mapped product. As a result, the fCO 2w shows a decrease of −0.06 (regional mapped product) to −2.64 (mapped) μatm/year during the wintertime and an increase of 1.37 (mapped) to 1.12 μatm/year (gridded) in the summer. Thus, during 2002-2010, the steady (mapped) or decreasing fCO 2w (gridded) are largely driven by the trends in wintertime similar to the results for one location as presented in Park and Wanninkhof (2012).
For the 2010-2018 time period, there are coherent positive trends for both winter and summer with SST increasing 0.2°C/year. The MLD increases, particularly strongly in the summertime, and the fCO 2w increases on the order of 2 to 4 μatm/year. Winter and summertime trends for fCO 2w are also more similar compared to the 2002-2010 time frame where the seasonal trends diverged. The increases in fCO 2w and SST are significantly higher for the gridded data than for the mapped product for the whole region, again indicating that the trends in fCO 2w and forcing functions are stronger in the areas along the cruises' tracks.

Journal of Geophysical Research: Oceans
To further elucidate the physical and biogeochemical factors contributing to the observed trends in fCO 2w over the entire period, the same deconvolution is used as for the seasonal attribution (equation (8)). The monthly values of δfCO 2w (SST), δfCO 2w (F CO2 ), and δfCO 2w (B&M) are determined from 2002 to 2018. Then the16-year average monthly values (Figure 7) are subtracted to deseasonalize the δfCO 2w data in order to determine the monthly anomalies over the 16-year record. To discern seasonal anomalous patterns, 3-month averages are used ( Figure 10). The deseasonalized record shows a general anticorrelation of δfCO 2w (SST) and δfCO 2w (B&M) for the times with large deviations. This suggests that during anomalous warming, there is increased biological drawdown, and during cooling higher fCO 2w likely caused by mixing. However, for small δfCO 2w , the δfCO 2w (SST), and the δfCO 2w (B&M) can act synergistically or antagonistically.
The δfCO 2w (F CO2 ) shows smaller seasonal anomalies compared to δfCO 2w (SST) and δfCO 2w (B&M). However, the δfCO 2w (F CO2 ) shows persistent trends and is negative at  Changes in MLD are an important contributor to the changes in δfCO 2w through its effect on SST, F CO2 , and B&M. The MLD influences ∂fCO 2w (SST) through cooling during MLD deepening and accelerated warming during shallowing. The ∂fCO 2w (F CO2 ) is impacted by MLD as deeper mixed layers cause smaller changes in fCO 2w due to the volume effect. The MLD deepening impacts ∂fCO 2w (B&M) through differing and opposing means. Increased nutrient supply can enhance productivity that lowers fCO 2w , but increased DIC supply will increase fCO 2w . As shown in Figure 11,   The air-sea CO 2 fluxes in the region are quantified by using the fCO 2w mapped product (fCO 2w,mapped ). It is merged with the monthly averaged fCO 2a determined from the MBL product to determine the monthly ΔfCO 2 . The OISST, SSS hycom , and the second moment of the winds, <u 2 > are used to determine the monthly fluxes according to equations (5) and (6).
Overall, the seasonal progression in fCO 2w controls an air-sea CO 2 flux pattern that is characteristic for the subtropical and tropical regions of CO 2 outgassing in the summer and invasion in the winter (Bates et al., 2014;Takahashi et al., 2002). The monthly averaged CO 2 flux, <u 2 >, and ΔfCO 2 over the 16-year period are shown in Figure 12  The differences are primarily caused by differences in uptake during the winter months. In summertime, the region is a CO 2 source and shows less year-to-year variability than wintertime.
Annual anomalies of note are the maxima in uptake in 2013 and 2017 that are caused predominantly by large wintertime uptake. In 2017, there also is reduced summertime evasion (Figure 12a). While high wintertime uptake for 2013 and 2017 are mostly caused by the lower ΔfCO 2 , in 2017 higher than average winter time winds contribute to the larger uptake as well. In addition to the effect high wintertime winds on the flux in 2017, low winds in the summer of 2015 significantly decreased the CO 2 evasion (Figures 12a and 12b).
The large changes in annual fluxes over preclude a clear indication that the CO 2 exchange in this marginal sea is appreciably different than the surrounding oceans as suggested in the ocean-dominated margin hypotheses of Dai et al. (2013). The hypothesis is that in waters the Caribbean Sea exchange with the open ocean through the passages and the subsequent vertical mixing causes slightly elevated fCO 2w and effluxes compared to the Atlantic Ocean basin. This was, in part-based interpretation of the results of Olsen et al. (2004), where Dai et al. (2013) suggest that fCO 2w are about 5 μatm higher in the Caribbean Sea than in the open ocean. As shown in Figure 12, 2003 had higher ΔfCO 2 and fluxes than seen in the following  (Takahashi et al., 2002) shows that on average the 14°N to 18°N latitude band in the open ocean is a weak source while the 18°N-26°N band is a net sink. These meridional trends and values are similar to observed in the region of the current study Caribbean Sea.
The means of determining the total CO 2 flux over the area has an impact on the results. All global and regional air-sea CO 2 flux estimates rely on interpolation of sparse data such that the large Caribbean dataset, and gridded data product can be used to provide a view of the differences. To estimate the effect of different temperatures, the average difference between the OISST product and measured temperatures on the ship of 0.25°C is added the OISST and the mapped fCO 2w was determined with the MLRs in Table 1. This leads to mapped fCO 2w values that are on average 2 μatm higher and the average specific uptake decreases by 25% from −0.20 to −0.15 mol m −2 year −1 . This illustrates that accurate temperature measurements and means to interpolate the fCO 2w and SST from the intake depth to the water surface are important (Woolf et al., 2016).
The fluxes determined using gridded ΔfCO 2 data, and the mapped ΔfCO 2 based on the MLRs yield differences as well. The former only covers part of the region and has temporal gaps. Qualitatively, the two products agree. Both show similar large seasonal variations with the gridded data showing slightly greater amplitudes (Figure 12a). However, the average fluxes for the entire period differ greatly. They are −0.19 ± 0.70 mol m −2 year −1 for the MLR based mapped product and −0.05 ± 0.70 mol m −2 year −1 for estimate based on the gridded data. The differences in magnitudes of fluxes are attributed to lack of observations in the winter months, along with the use of observed SST in the gridded observations and OISST in the gridded mapped product. The OISST yields more negative fluxes as described above.

Relationship With Climate Reorganizations
The changes in fCO 2w and CO 2 fluxes on seasonal to interannual timescales are closely related to SST and MLD. In turn, the regional scale SST anomalies are associated with changes in the physical environment influenced by large-scale climate reorganizations. The climate cycles in the Caribbean are linked, through teleconnections, to the El Niño-Southern Oscillation, ENSO; the Atlantic Multidecadal Oscillation, AMO; and North Atlantic Oscillation, NAO (Birkmark, 2014;Ibánhez et al., 2017;Lefèvre et al., 2013Lefèvre et al., , 2014Thomas et al., 2008). Since the meteorological and oceanic responses to the climate indices are complex and sometimes lagged, only a broad view of responses of air-sea CO 2 fluxes to these climate patterns are investigated using the seasonal indices that have the greatest impact, for the NAO this is the wintertime (December-February), for the AMO spring (March-May), and for the ENSO the early winter (November-January).
Time series of these seasonal climate indices and annual CO 2 fluxes are shown in Figure 13, each showing significant anomalies in the middle part of the record (2009)(2010)(2011) where the change in flux trends occurred. The annual fluxes in Figure 14 are plotted at year's end to reflect the expected lag between the indices and the annual cumulated flux. The overall trend in the NAO is positive for the time period, much like the preceding two decades (Thomas et al., 2008), but it shows a large minimum in 2010 that corresponds to the minimum in CO 2 flux (maximum uptake). The AMO shows a large negative followed by a positive excursion in 2009/2010 similar to the reversals in SST and fCO 2w anomalies, and the ∂pCO 2w (SST) (Figure 11) observed during this time. These changes in annual fluxes in response to the AMO are opposite to the changes observed in the Equatorial Atlantic (Ibánhez et al., 2017). The air-sea CO 2 fluxes show a positive relationship with ENSO with greater uptake during the negative phase of the ENSO (El Niño conditions). This is synergistic with the decreased evasion in the Equatorial Pacific albeit fivefold less in magnitude compared to the El Niño induced depressed outgassing in the Equatorial Pacific

10.1029/2019JC015366
Journal of Geophysical Research: Oceans upwelling region (Feely et al., 2006). The mechanisms for the lower fluxes in the Equatorial Pacific and Caribbean differ as well. The increased uptake in the Caribbean in response El Niño is due to prevailing lower SSTs. In the Equatorial Pacific, the decrease fCO 2w is due to decreased upwelling of waters with high DIC. In contrast, the greater uptake in the Caribbean is antagonistic with the Western Equatorial Atlantic where outgassing is enhanced during El Niño phase of ENSO (Lefèvre et al., 2013). The last year of the record (2017) shows an increase in CO 2 uptake along with a decrease in the ENSO and NAO indices. Although the magnitudes of the indices are not extreme, the change from 2016 to 2017 are large for the NAO and AMO. The increasing MLD and <u 2 >, and decreasing SST in response to the decreasing NAO and increasing AMO, leads to increased CO 2 uptake.
The correspondence of the fluxes with climate indices is apparent, in particular when considering year-toyear changes, in addition to absolute magnitude. The ENSO and NAO minima correspond to lower SST and deeper MLD contributing to lower fCO 2w and a minimum in flux, or greater uptake. The change from low trend in fCO 2w in 2010/2011 to a much greater rate of increase corresponds to extrema in all the climate indices. The winters of 2010 and 2011 showed anomalous low SSTs for much of the North Atlantic, including the Caribbean (Buchan et al., 2013). Rivière and Drouard (2015) point out the contrasting phase of the NAO during 2010 and 2014. The negative NAO index in 2010 corresponds to low SST and low fCO 2w while maximum in the NAO in 2014 is associated with higher SSTs and fCO 2w . This suggests that the 3-to 8-year trends are strongly influenced by large single-year anomalies in the climate indices.
The effect of the AMO in forcing of the fluxes is less apparent but the large swing in 2009/2010 corresponds to large changes in the trend of fCO 2w . Indeed, the midpart of the record show that all indices have large anomalies ( Figure 13) corresponding with a minimum in fCO 2w and largest uptake. However, the multifacetted impacts and monthly variation of the indices along with the large seasonal cycle of fCO 2w precludes a simple statistical or strong causal relationship for the entire period.

Conclusions
Through a corporate, academic and federal partnership the Caribbean Sea has changed from an undersampled area for surface water CO 2 and air-sea CO 2 fluxes to one of the best-observed regions in the world. The high sampling density and robust interpolation methods utilizing the strong correlation of fCO 2w with SST and other variables make it possible to discern strong decadal trends for the whole area in spite of seasonal variability that is 10 times greater. The region shows an average net uptake of CO 2 of −11.1 Tg C year −1 for 2002-2018 with annual fluxes ranging from +8.2 (evasion) in 2003 to −22.4 Tg C year −1 in 2017 in contrast to the annual climatological values in the region of 0.4 Tg C year −1 (Takahashi et al., 2009). Our results show that the general view of net outgassing in low latitude seas needs to be reassessed. Moreover, decadal variability of regional air-sea fluxes in oligotrophic marginal seas should be accounted for. In this work, much lower than expected fCO 2w increases prevailed for 8 years with major transition in patterns and trends occurring around 2010, closely associated with changes in SST trends and MLD that have a direct impact and indirect effect, though biological processes, on fCO 2w . The fCO 2w minima corresponds to anomalous low wintertime SSTs and extrema in several major climate indices, most notably the NAO. The fCO 2w hardly changed from 2002 to 2010 and strongly increased from 2010 to 2018. Over the entire time period the surface ocean fCO 2w lagged the atmospheric CO 2 increase by 24% causing the uptake to increase over time. Annual fluxes show increasing uptake from 2002 to 2010 and a small decrease in uptake thereafter with a reversal in the last year full of the record (2017) to larger uptakes. The changes in flux are less dramatic that the fCO 2w trend from 2010 to 2018 as the XCO 2a and fCO 2w are both increasing such that the ΔfCO 2 , and thus, the flux does not change as much as fCO 2w . The trends and pattern are similar for the observations and the interpolated products, but magnitudes vary indicating that quantitatively assessing changes in fCO 2w and air-sea CO 2 fluxes requires both high density and quality data without temporal gaps, and robust approaches for spatial and temporal interpolation.