Our study focuses on a limited region of the Central Labrador Sea, spanning from 55.5°N to 57.5°N and from 51.5°W to 53.5°W (orange box, Figure 1 ). The specific area selected was based on mapping of past deep convection events (Marshall and Schott, 1999), which has led to deployments of several moorings in this location, and is believed to be representative of a significant area of the deeper mixing within the Central Labrador Sea. However, the deep convection activity in the Central Labrador is also known to be dynamic, and there is inter-annual variability of the area of deep water formation (Rühs et al., 2021). The red box in Figure 1 shows, for comparison, the resolution of early global-based pCO₂ products produced by Takahashi et al., 2002 and Takahashi et al., 2009 (4°x5° grid). In contrast, a smaller grid cell (2°x2°) was chosen for analysis in this study to exclude data-points collected over the continental shelves and slopes that surround the Central Labrador Sea.

The “non-SOCAT” underway measurements presented and utilized here (see Figure 1 ) are from the Atlantic Zone Off-Shore Monitoring Program (AZOMP) and Atlantic Zone Monitoring Program (AZMP) of Canada’s Department of Fisheries and Oceans (DFO), with both programs taking place on-board CCGS Hudson, between 2016 and 2019 (CCGS Hudson data submitted to SOCAT versions 2023 and 2024). Further, underway measurements collected from the CCGS Amundsen are also included in this study, with data from 2017 to 2020 (pending submission to SOCAT). Finally, observations over the Atlantic continental shelves (data: Cyr et al., 2022; Gibb et al., 2023) are also included in Figure 1 (Labrador Shelf area) to show possible opportunities for future expansion of the Canadian pCO₂ observation network. All of these “non-SOCAT” data were quality controlled and are either submitted or pending submission to newer versions of the SOCAT database ( Supplementary Table 4 in Supplementary Materials ). In the meantime they can be requested from the authors listed in the data availability section of this manuscript.

Within the grid cell chosen in this study (approximately 50,000 km², orange box, Figure 1 ), the “non-SOCAT” underway data includes three crossings (13 days of observations), from the years 2015, 2016 and 2018. In comparison, the SOCATv2021 dataset includes only two crossings of the area (four days of observations), from 2008 and 2016.

It is important to point out the importance of the mooring datasets presented here (details in Table 1 ) that made this study possible due to the temporal/seasonal coverage that they provide in comparison to ship-based studies. When ship-based underway-observations are so sporadic and limited, observations from mooring deployments become essential. The discussion in the remainder of this paper focuses on the data collected within the 2°x2° orange box in the Central Labrador Sea ( Figure 1 ).

The global pCO₂ products compared here are the climatologies of Takahashi et al. (2002); Takahashi et al. (2009) and Fay et al. (2023) (discussed in Fay et al., 2024); these climatologies are referred to here as T2002, T2009 and T2023, respectively. The T2023 climatology is provided as ΔfCO2, we therefore recalculated fCO2 by adding the atmospheric fCO2 (for reference year 2010 as in Fay et al., 2024), and converted to pCO2 using surface temperature (Weiss, 1974) for the reference year 2010 (Multi Observation Global Ocean ARMOR3D L4 – Copernicus - Guinehut et al., 2012).

We also compare six other observation-based products from the harmonization of pCO₂ products provided by Gregor and Fay (2021) and discussed in Fay et al. (2021), which includes multiple linear regression models (MLR), machine learning ensemble (ML6), mixed layer scheme (MLS) and three neural network-based models (references and details of each product are given in Table 2 ). All these products used SOCAT observations for reconstruction of pCO₂ on a model grid.

We calculated a regional pCO₂ climatology directly from pCO₂ observations and compared it with the climatologies calculated from the pCO₂ global products. We also used the pCO₂ observations and pCO₂ from each of the global products to calculate air-sea CO₂ fluxes, and consequently the climatologies of the fluxes.

The air-sea CO₂ fluxes (FCO₂ ^(air-sea)) were calculated using the following equation (Equation 1):

where kCO₂ is the transfer velocity according to Wanninkhof (1992), K₀ is the solubility constant following Weiss (1974), with pCO₂ ^(air) and pCO₂ ^(sea) as the atmospheric and ocean surface pCO₂. From here on we will use ΔpCO₂ = pCO₂ ^(air) – pCO₂ ^(sea). The air-sea flux is dependent on the transfer velocity (kCO₂), which is strongly dependent on wind speed. Since we are dealing with a data-poor region, we fixed the choice of wind product and wind parameterization in order to focus on the effect of pCO₂ data-coverage. Specifically, we used wind speed estimates from the reanalysis product ERA5 (Hersbach et al., 2020), which has been widely used in global products, including those compared here, due to its high spatio-temporal resolution for the flux calculation (see also Atamanchuk et al., 2020).

We calculated the monthly uncertainty of air-sea CO₂ fluxes (∂FCO₂) by propagating the monthly errors of wind (∂U) and ΔpCO₂ (∂ΔpCO₂), which we believe to be the largest sources of uncertainties for this climatological approach, using the following equation (Equation 2):

The monthly climatology, i.e. the seasonality of sea-surface pCO₂, was calculated as a monthly average and monthly standard deviation to create our climatological reference for the Central Labrador Sea using pCO₂ observations from 2000 to 2020. The climatology was calculated with the pCO₂ and fluxes data compiled for this study and for each of the global products of monthly time series. We then analyze how well the global products compare with the directly observed seasonal variability of pCO₂ and air-sea CO₂ fluxes in the Central Labrador Sea.

To compile pCO₂ observations collected over 21 years for a climatological monthly averaging approach, it is necessary to correct for the increase in atmospheric (and surface ocean) pCO₂ over time. For that, we used the Icelandic atmospheric time series between 1992 and 2020 (Dlugokencky et al., 2021), as well as observations from Sable Island between 1993 and 2019 (Worthy, 2023). For the Iceland station, an increase of 2.16 µatm/year was found, and for the Sable Island station there was an increase of 2.08 µatm/year. Both time series showed a similar rate of increase, with the slopes of the two linear least squares regression being statistically indistinguishable (p-value>0.05). This rate of atmospheric increase used here is consistent with the 2.2 µatm/year rate reported in Raimondi et al., 2021, for the period 1996-2016 in the same region. Therefore, for simplicity, we used 2.1 µatm/year for adjusting surface water pCO₂ to the common reference year of 2020.

Figure 1 shows that there are large data gaps throughout most of the Labrador Sea domain, especially in the Central, Northern, and Labrador Shelf regions. We decided to show all observations in and around the study area (even though they are not all used in our analysis) to emphasize the major observational gap that exists in this region, despite the region’s potential significance for exchange of gases and carbon between the atmosphere and the deep ocean. Even with the addition of the “non-SOCAT” data presented here, we are still far from having anything close to representative observational coverage for the Northwestern Atlantic Ocean, including the Labrador Sea (Central and Northern) and Canadian shelves (see Duke et al., 2023b).

The pCO₂ data from the Central Labrador Sea moorings ( Figure 1 ; Table 1 ) show a strong seasonal cycle ( Figure 2 ), with relatively high, near-equilibrium pCO₂ in winter (JFM) weakening the uptake of atmospheric CO₂, and low pCO₂ values in spring (AMJ) and summer (JAS), increasing the difference with the atmospheric pCO₂ and thus driving a strong CO₂ sink. The timing of the decline in pCO₂ in mid-spring (referred to here as the “spring-decline”) varies from year to year, and a second less pronounced drop in pCO₂ may occur in the fall (OND) as well. This overall seasonality is driven partly by biological activity, including a strong decrease of pCO₂ coinciding with the start of the spring bloom, and partly by abiotic controls (i.e. changes in temperature and vertical mixing) as pCO₂ increases steadily after summer until the end of winter. As the wintertime cooling sets in, increased solubility would drive CO₂ fluxes into the ocean, however, the deepening of the mixed layer carrying a high pCO₂ signal from respiration are mixed into the surface layer, thus driving CO₂ fluxes out of the ocean (outgassing) (Körtzinger et al., 2008a; Martz et al., 2009). This will lead to a maximum winter-time pCO₂ as observed in other high latitude regions (Iida et al., 2015).

The variability of our pCO₂ climatology increases after including the underway observations with the mooring data in the analysis, however the overall seasonality remains consistent and well represented ( Figure 2 ), with the expected high pCO₂ in winter (409 ± 7 µatm), followed by a steep decline through spring until mid-summer, when it reaches the minimum in July (down to 250 µatm). After averaging all available data, both from moorings and underway systems, we can confirm that over the entire seasonal cycle the region acts as a year-around sink for atmospheric pCO₂ (i.e. pCO₂ ^(sea) < pCO₂ ^(air)). There are only a few days when ocean surface pCO₂ may exceed atmospheric pCO₂, and this may happen right before the “spring-decline”, when observations show increased variability (during May and June).

As seen in Figure 2 , even after the inclusion of the underway observations, the majority of the observations discussed here come from the three mooring deployments, showing their major importance for this otherwise under-sampled region. Using the monthly average of the pCO₂ observations shown in Figure 2 , we produced an estimate for the seasonal climatology for the Central Labrador Sea, which is used as the reference “observation-based” climatology.

When comparing pCO₂ seasonal climatologies from the different global products ( Figures 3 , 4 ), they also characterize the Central Labrador Sea as a region of atmospheric CO₂ uptake (sink). Most products follow the overall pattern seen in the observations ( Figure 2 ). However, the timing and amplitude of the seasonal cycle of pCO₂ is not consistent between the products and the observational data. For example, most products indicate an earlier “spring-decline” of pCO₂ compared to the observation-based estimate (March/April vs May, Figure 2 ). There is also a shift of timing of the summer minimum (earlier summer minimum, except in Takahashi et al. (2009), that shows a minimum in August). Most of the products underestimate pCO₂ in winter, spring and summer when compared to the observation-based estimate (black line), with a bias (product - observations) over time ranging from -80 to +40 µatm ( Figure 4 ). To a lesser degree, there is also an overestimation of pCO₂ by most products in late-summer and fall. The products MPI, JENA, NIES and JMA showed the highest seasonal amplitudes (winter maximum – summer minimum) of around 150 µatm, the observation-based estimate however showed an even higher amplitude of almost 200 µatm.

We note an especially strong difference between the two early Takahashi climatologies (Takahashi et al., 2002, 2009), with the seasonality from Takahashi et al. (2009) being the least consistent with observations and with the other global products. These differences may be due to the different extrapolation techniques combined with the fast increase in observations of pCO₂, with the new observations being mostly on the border of the 4°x5° grid (see Figure 1 , red box) around the Greenland shelf and slope by the Nuka Arctica underway system (Olsen et al., 2008), therefore potentially skewing the expected seasonality of the Central Labrador Sea for this climatology.

We keep these earlier Takahashi climatologies in the discussion since they have been used as benchmarks for comparisons in earlier studies (e.g. Lüger et al., 2004; Körtzinger et al., 2008b; Landschützer et al., 2014; Lauderdale et al., 2016). The new Takahashi climatology (T2023 - Fay et al., 2023) shows a better seasonal cycle (based on SOCATv2022 with more observations available), agreeing with the other products, although showing an early increase in pCO₂ in the summer, and an overestimation of pCO₂ from summer through fall. Other more recent gap-filling methods discussed in the study may, however, be more appropriate for this regional scale analysis due to their finer resolution, although they also have their strengths and weaknesses based on statistical metrics.

Of the pCO₂ products with a 1°x1° degree resolution, the CSIR product (Gregor et al., 2019) has the lowest mean absolute error (MAE) and root mean square error (RMSE) related to the observation-based climatology of this study (MAE=16.5 µatm; RMSE=30.7 µatm). This is influenced by its summer-fall values being almost equal to our observation-based product (MAE=3.56 µatm). The CSIR winter-spring values are underestimated compared to the observations, however their bias is on the low-end compared to the other products. The annual average bias was slightly more negative for CSIR (-14.3 µatm) than JMA (-13.4 µatm) (Iida et al., 2021), however further examination of the monthly bias ( Figure 4 ) shows multiple months with positive bias in the fall that will cancel out some of the negative bias in the winter. The higher spread of bias values for the JMA product is reflected in its higher MAE and RMSE (MAE=21.2 µatm; RMSE=27.3 µatm). Recent products with finer resolution (0.25°x0.25°,e.g. Chau et al., 2024; Gregor et al., 2024) could lead to improvements in the results for this region.

Although the data products discussed here are all intended as global-scale products, these should be tested to assess their skill in different basins or even at regional levels, such as in this study. Rödenbeck et al. (2015), for example, recommended checks on the consistency between such products and the use of multiple products in such comparisons.

Figure 5 shows the seasonal variability of the calculated air-sea CO₂ fluxes and the comparison of these observation-based estimates with the same global products discussed above. Similar to pCO₂, the majority of estimates of air-sea CO₂ fluxes from the global products indicate a pattern of overestimation from winter to spring and underestimation from summer to fall when compared to the observation-based estimate in this study. The monthly average bias (product - observations) for each global product is shown in Figure 6 , ranging from -0.42 to +0.5 molC m^-2month^-1, with most products showing stronger fluxes than the observations in the first half of the year (winter: -0.42 to +0.38 molC m^-2month^-1; spring: -0.35 to +0.16 molC m^-2month^-1), and weaker fluxes in the second half of the year (summer: +0.12 to +0.28 molC m^-2month^-1; fall: -0.09 to +0.52 molC m^-2month^-1).

Bias values are negative (stronger fluxes than observation-based estimate) for winter-spring and positive (weaker fluxes than observation-based estimate) for summer-fall. The seasonal positive and negative biases tend to cancel out, thus leading many products to have an annual average bias that is low. Hence the mean absolute error (MAE) and RMSE metrics (see Table 3 ) are more appropriate for discussion here. Of the 1°x1° degree products, CMEMS, JENA and JMA show the lowest MAE and RMSE. The climatology from Takahashi et al., 2002 also shows good metrics when compared to the observations, which does not hold for the newer version of Takahashi et al., 2009 (T2009), as shown in Figures 6 , 7 . Notably, T2009 stands out from the other products in showing a strong overestimation of fluxes during winter. The climatology from Fay et al., 2023 (T2023) shows a slight improvement in the metrics when compared to T2009 (T2023 includes additional observations from SOCATv2022). Overall, the JMA product has the lowest bias, however the MAE and RSME metrics in both pCO₂ and flux suggest the CMEMS, JENA and JMA as the best options when compared to the reference observation-based climatology.

When averaged seasonally, summer and fall are the seasons with the highest fluxes based on observations, and are also the seasons with larger inconsistencies between the global products and observation-based estimate ( Figure 7 ; Table 3 ). In the fall, only the T2009 product corresponds closely to the observations (with overlapping error-bars) and in summer, only T2002 and NIES. The NIES product is the only product that classifies the Central Labrador Sea as a source of CO₂ to the atmosphere during winter and shows small positive values (close to equilibrium or a weak source) within the seasonal variability during the fall. The MPI product overall classifies the winter as a sink, but its large uncertainty does not preclude that some winters may have out-gassing periods. T2023 and NIES are the products with higher amplitudes of 0.62 and 0.64 molC m^-2month^-1, respectively. Annually, all products show consistent representation of an ocean sink, with the observation-based estimate being -4.0 ± 2.2 molC m^-2year^-1, and most products showing a slight underestimation of the flux when compared to our observation-based estimate ( Figure 8 ).

High latitude regions such as the Labrador Sea are amongst the poorest in terms of pCO₂ data coverage, even while their significance for global air-sea fluxes and net carbon storage is high. Therefore, high latitude regions usually fall within the regions with highest uncertainty and errors in both regional and global gap-filling estimates (Gloege et al., 2021). For example, there remains controversy whether the Southern Ocean acts as a strong or weak sink (Sutton et al., 2021). However, the seasonality at the regional scales and the strength of the inter-annual variability is poorly characterized, due to large winter-gaps in observations increasing the uncertainty in the Southern Ocean (Mackay et al., 2022; Wu and Qi, 2022). Similarly, in the North Atlantic Ocean, the Labrador Sea is one of the regions with the fewest observations, also leading to high uncertainties.

Here we identify some key sources of uncertainties for the air-sea flux estimates. Firstly, the choice of wind products and wind parametrization for the bulk-formula calculation of CO₂ fluxes are among the most important sources of errors. Different parameterization choices for the gas transfer coefficient can alter the intensity of the CO₂ fluxes estimates in the region by an average of ± 20% or approximately 0.08 molC m^-2month^-1 (Atamanchuk et al., 2020). In this study, the parameterization of Wanninkhof (1992) was used for our observation-based reference climatology and for the global products compared in this study.

When using different wind products (e.g. NCEP and CCMP products) with the same parameterization, the observation-based estimates can vary by as much as -0.17 molC m^-2month^-1, notably during the period of strong summer uptake. This can lead to an almost 50% decrease of the intensity of the summertime carbon sink in the Central Labrador Sea, by switching from ERA5 to NCEP. Differences between CCMP and ERA5 are less pronounced, with a maximum difference around -0.05 molC m^-2month^-1 during the summer. Overall, differences between fluxes calculated using these three wind products are largest in summer and fall, and less pronounced in winter and spring, being consistent with the flux formulation and ΔpCO₂, that is, when air-sea pCO₂ gradient are larger, differences of using different wind products are also more pronounced. In this study, ERA5 was used to calculate CO2 fluxes, for the estimate presented in this study and for the global products.

The uncertainties of wind products and wind parametrization have been discussed previously (e.g. Moore et al., 2008; Koelling et al., 2017; Woolf et al., 2019; Atamanchuk et al., 2020), and clearly represent a major problem for global estimates of air-sea CO₂ fluxes. Another source of uncertainty is the measurement of the surface temperature, which in most cases occurs at the depth of a ship’s intake of water, which is located typically well below the air-sea interface (e.g. at 5–10 m). This implies a need for an adjustment of the temperature (and pCO₂) to reflect actual surface conditions (Watson et al., 2020). Finally, the cool and salty skin-temperature effect offers potential for major bias, with these two factors together having the potential to increase the global oceanic uptake by as much as factor of two (Watson et al., 2020). Other important sources of uncertainties include: (1) uncertainty related to the pCO₂ measurements (Bender et al., 2002; Wanninkhof et al., 2019; Dong et al., 2024); (2) gap-filling model uncertainties (e.g. data-coverage uncertainty; Duke et al., 2023a); and (3) uncertainty from the wind measurements that feed global wind-products (Roobaert et al., 2018; Chiodi et al., 2019; Wright et al., 2021; Fang and An, 2022).

Figure 9 shows the seasonal relationship of air-sea CO₂ fluxes, pCO₂ (or ΔpCO₂), surface temperature and wind. The fluxes are more intense and more variable starting in spring, through summer and fall. The contribution of each variable towards predicting the variability of CO₂ fluxes was explored using multiple linear least squares regression. We found ΔpCO₂ alone was able to describe 62% of the calculated flux variability (R² = 0.62 and p-value = 0.0023). Based on a regression of CO₂ fluxes with both ΔpCO₂ and wind, given the amplitude of variability, these two variables were able to describe 84% of the flux variability (R² = 0.84 and p-value = 0.0002), thus, the variable wind (U) improves the regression together with ΔpCO₂. However, we found that wind alone cannot explain the variability of fluxes (large p-value). Temperature alone also cannot explain the variability of fluxes in the Central Labrador Sea (large p-value), due to the opposite expected effect of the temperature changes in ΔpCO₂ and thus in the fluxes as well. Therefore, measurements of surface ocean pCO₂ remains the most important variable for constraining and improving the estimates of air-sea CO₂ fluxes in this region (consistent with this, Dong et al., 2024 found larger standard deviations in reconstructions due to a recent decline in SOCAT observations), followed by resolving/improving the wind (U) data products and parametrization. It is important to point out that these relationships are particular to our region of interest and also dependent on the type of gas exchange parameterization (in this case, using Wanninkhof (1992)).