In the current era of increased CO₂ emissions and ocean acidification, it is important to understand the temporal changes in the ocean carbonate system variables, i.e. dissolved inorganic carbon (DIC), total alkalinity (TA), pH, and the fugacity of carbon dioxide (fCO₂). Opportunistic sampling and sparse observations through space and time make it difficult to estimate the trends and variations in the ocean carbonate system, especially in the ocean interior. As a result, the literature has focused mainly on temporally averaged spatial climatologies (Bennington et al., 2022; Broullón et al., 2020; Lauvset et al., 2016; Takahashi et al., 2014; Sabine et al., 2004). Part of the reason for this emphasis is that carbonate variables are not consistently measured through time. Large data gaps in time induce large uncertainties, which favor estimation of decadal-scale trends, but not on finer time scales. We seek to improve the time resolution of the carbonate variables to a monthly scale by using observation-based methods. Here, we refer to observation-based, or observation-centric, methods as statistical approaches using direct observations of the ocean carbonate system from moorings, cruises, floats or underway sensors, as opposed to numerical model estimates. By amalgamating multiple observation sources and by using the chemical equilibrium relationships between variables, we can exploit the observation-based information to its fullest. We develop a methodology that produces a monthly time-depth gridded product of the carbonate system, with uncertainties that vary over time. Our application is focused on carbonate system variables in the Labrador Sea and how they have changed over the last 24 years.

The ocean carbonate variables (DIC, TA, pH, and fCO₂) have often been examined individually. Surface ocean fCO₂ is used to calculate the CO₂ air-sea flux, thus quantifying the magnitudes of the ocean carbon sink (Landschützer et al., 2016). Moving into the water column, DIC has been used to investigate the spatial distribution and storage rate of anthropogenic CO₂ (Müller et al., 2023; Gruber et al., 2019). With the increase of anthropogenic CO₂ in the ocean, pH measurements quantify ocean acidification and the impact on marine life (Rastrick et al., 2018). Finally, total alkalinity has been used to understand the ocean’s buffering capacity in response to changes in pH (Middelburg et al., 2020). TA is also often used with one of the other carbonate variables to then estimate the rest of the carbonate system (Carter et al., 2017).

A few studies have examined the carbonate variables from a multivariate perspective (Gregor and Gruber, 2021; Carter et al., 2021; Bittig et al., 2018; Sauzède et al., 2017). While these multivariate studies estimated the four carbonate variables, the variables themselves were analyzed independently, as in without correlation or connection between the carbonate variables. This is important since the carbonate variables are related (and constrained) through chemical equilibrium equations (Sarmiento and Gruber, 2006). CO2SYS is a family of algorithms with implementations in MATLAB (Sharp et al., 2020), Python (Humphreys et al., 2023) and R (Gattuso et al., 2022), that solves the chemical equilibrium equations to calculate the state of the carbonate system that consists of the four measurable variables (DIC, TA, pH, and fCO₂), but also includes the aragonite saturation state (Ω), DIC speciation and the Revelle factor. Few observation-based methods have incorporated the chemical connection between carbonate variables within the analysis method, but it is often done as a data processing step to fill in gaps in a set of observations prior to being analyzed (Lauvset et al., 2021; Mackay and Watson, 2021) to transform the results to another variable (Jiang et al., 2019), or as a consistency check (i.e. CONTENT, Bittig et al., 2018).

We used two observation-based data sources for carbonate variables: GLODAP with subsurface bottle sampled data (Lauvset et al., 2022) and SOCAT with surface fCO₂ data from underway systems on ships (Bakker et al., 2016). There is also BGC-ARGO profiling floats with subsurface pH data (Bittig et al., 2019). It is most common to use one data source (e.g. GLODAP) as input data for studies, and the others (e.g. SOCAT and/or BGC-ARGO) as validation or comparison data (Bittig et al., 2018). Few papers have incorporated multiple data sources as inputs into their analysis, with the exception of Feely et al. (2023); Iida et al. (2020); Gregor and Gruber (2021) and ECCO-Darwin, which is a data-assimilative biogeochemical model that used SOCAT, GLODAP, SOCCOM and BGC-ARGO data as observational constraints in their model (Carroll et al., 2022). Combining multiple data sources is challenging, as different observation types will have different measurement uncertainties, as well as spatial and temporal sampling schemes. Research on these carbonate variables often reflects the spatial distribution of their observations, with relatively little analysis on fCO₂ in the interior ocean (Fiedler et al., 2013), but significant work analyzing fCO₂ in the surface waters (Chau et al., 2024; Bennington et al., 2022; Gloege et al., 2022; Yasunaka et al., 2019; Du et al., 2015; Zeng et al., 2014; Landschützer et al., 2013).

In this study, we estimate monthly time series of the four carbonate system variables in the Labrador Sea from 1993-2016 using a fusion of GLODAP and SOCAT data. We used a state space model approach that produced time-varying uncertainties based on data availability and data type. We also improved our effective observation availability by combining direct observations with pseudo-observations derived from observation data with CO2SYS. This methodology could be generalized to different regions, temporal or spatial scales, with the caveat that data availability and its temporal resolution is a limiting factor for obtaining meaningful results.

Our region of interest is the Labrador Sea, which we define as the region south of 65° N and west of 40° W within the Atlantic Arctic (ARCT) Longhurst biogeochemical province (Longhurst, 2007) ( Figure 1 ). Within the Longhurst province, biogeochemical properties are assumed to be homogeneous. The Labrador Sea in the northwest Atlantic is a key area for anthropogenic carbon storage due to its association with deep convection (Raimondi et al., 2021).

Two datasets were used in this study: GLODAPv2.2022 (Lauvset et al., 2022) and SOCAT2022 (Bakker et al., 2022, 2016). In our study region, these datasets provide complementary observations: GLODAP has direct observations of DIC, TA, pH, fCO₂, temperature, salinity and oxygen in the interior ocean for depths 0-4500 m from 1981-2021; and SOCAT has surface fCO₂ from 1990-2021. This version of GLODAP in the North Atlantic has large data gaps of multiple years before 1991 and after 2016. Our analysis will span 1993-2016 to minimize gaps larger than a year where no data for any carbonate variables exist. However, 1993 was chosen as a start date to mimic other oceanographic reanalysis data products with the same start date. Large gaps at the beginning and end of the time series will cause the overall decadal trend to be underestimated. We included SOCAT data from all cruises (flags A-D), but only GLODAP and SOCAT data with a WOCE flag 2 (good), thus removing the interpolated data for the GLODAP carbonate variables. The GLODAP and SOCAT observations within the Labrador Sea were separated into 20 depth layers and then averaged monthly. Layer divisions were chosen to have an approximately equal number (∼500) of individual DIC observations from GLODAP in each layer. This is related to near-surface DIC observations occurring in 40 months (14% of time series). The top of the depth layers in this study are: 0 m, 20 m, 53 m, 104 m, 209 m, 357 m, 505 m, 683 m, 890 m, 1100 m, 1360 m, 1580 m, 1830 m, 2090 m, 2370 m, 2590 m, 2810 m, 3000 m, 3200 m, and 3380 m.

Variables are sometimes reported for different ocean conditions, i.e. in situ ocean vs. laboratory temperature and pressure. Ocean conditions are defined by potential temperature (T), practical salinity (S) and pressure (P). We converted all variables to in situ ocean conditions to ensure compatibility across datasets. Subsurface ocean observations of fCO₂ are often standardized to laboratory conditions (P=0 dbar, T=20°C, and S=35 PSU), as is the case in the GLODAP dataset, since fCO₂ is not conservative to changes in T, S or P (Sarmiento and Gruber, 2006). Meanwhile, SOCAT reports surface fCO₂ at in situ conditions (P=0 dbar, T between −1.8 and 21.8°C, and S between 17.65 and 35.47 PSU). In order to make these two datasets comparable for use in the same analysis, GLODAP fCO₂ was converted to in situ conditions using chemical equilibrium equations (CO2SYS) and its corresponding input parameters: DIC, TA, T, S and P. The equilibrium constants used in CO2SYS are as described later in Section 3.1. pH was reported on total scale and in situ conditions by GLODAP.

GLODAP observations of the ocean interior are obtained from research cruises, and in our study region, observations are seasonally biased towards months May-July (63% of observations in the top 5 m are from May-July). The most sampled variable in GLODAP is DIC, followed by TA. Within GLODAPv2.2022, fCO₂ through the water column is rarely measured, but fCO₂ is more frequently available at the surface (notably between 2003-2012) owing to the ships of opportunity that collect the SOCAT data. The SOCAT data does not have a seasonal sampling bias (i.e. in the top 5 m there is a similar number of observations in each month). The sampling density of pH in GLODAP has been low at 0.08% of months, but with BGC-ARGO, pH availability since 2018 in the Labrador Sea has increased to a full 100% monthly coverage. Extra pH data from BGC-ARGO was considered to be included in the main results, but further investigation and quality control are needed, as we found the BGC-ARGO pH data has unexpected increasing trends in the Labrador Sea from 2018 to 2022. Instead, pH results with BGC-ARGO data are in the Supplementary Material ( Supplementary Figure S19 ) and serve to demonstrate how the state space model is able to adapt to changes in time series characteristics present in observations.

We also carried out data augmentation, that when DIC and TA were available, fCO₂ and pH were calculated by CO2SYS, as denoted by CO2SYS(DIC, TA). CO2SYS solves the chemical equilibrium equations for the carbonate system (Lewis and Wallace, 1998) (see Section 3.1 for more details). The pseudo-observations fill in gaps when direct, monthly observations are not available, and this greatly increases the effective availability of fCO₂ information. Note, that the data augmentation was performed after monthly averaging the data, as further described in Section 3.2 Implementation.

Since the datasets differ in time range, depth range, and spatial and temporal sampling schemes ( Figure 1 ), we checked how their T, S and fCO₂ data compared. The following comparisons were done as time series plots: (i) T and S of GLODAP and SOCAT, (ii) fCO₂ direct from SOCAT vs fCO₂ calculated via CO2SYS(DIC,TA) and GLODAP data. In the top 20 m, for temperature and salinity, SOCAT has a larger standard deviation than the others, but after being monthly averaged the three data sets have equivalent overall mean and standard deviations across our time period of interest ( Supplementary Figure S20 ). The fCO₂ calculated from GLODAP DIC and TA were compared to fCO₂ observations from SOCAT, and the calculated fCO₂ had poor precision but were sufficient for our analysis. 37 observations were able to be co-located within a month, in the top 5 m and the closest points within 5° of one another. The differences ( Supplementary Figure S16 ) showed an average bias of 1.56 µatm, with GLODAP calculated values slightly overestimating the fCO₂ on average and there is also a large RMSE of 50.26 µatm representing a large error between the calculated fCO₂ and the closest direct measurement within the month. Part of this bias would be due to measurement uncertainty of the input variables (DIC and TA) and the uncertainty in the equilibrium constants used in CO2SYS (Álvarez et al., 2020; Millero, 2007). The equilibrium constants are valid only for a set range of temperature and salinity values, however, in the subpolar region of the Labrador Sea, our temperatures do sometimes go below the recommended range.

The goal of the analysis is to estimate the state of ocean carbonate variables through time and over depth for the domain of interest of the Labrador Sea. x t is the state of the carbonate system at a given time t, and is as a column vector of length n_x = 4L, where L is the number of depths-layers (see Section 2). The state vector contains all the variables DIC, TA, pH and fCO₂ at each depth-layer through the water column (from surface to deep). The analysis times range from t = 1 , … , n t on a regular interval (in our case, monthly). Estimation of the state vector over time makes use of all available observations, and the framework used for the analysis is a state space model. This is comprised of two equations: an observation Equation 1 and a persistence Equation 2, which are defined below. Once these are defined and their parameters specified, the carbonate state is then estimated through time using the Kalman filter/smoother algorithm (Anderson and Moore, 1979).

The observation equation relates the state ( x t ) to available observations of the carbonate system through the following equation:

Here, y t   is a column vector that contains the observations of the carbonate system at time t. This vector is of length n y t and the length varies at each time step with the number of observations available at time t. The y t vector includes observations of the different carbonate variables from multiple datasets including the pseudo-observations (see Supplementary Material ). A t is an n y t × n x observation indicator matrix comprising 1’s and 0’s to match each observation in y t with its corresponding term in the state vector, x t . A t allows the simultaneous use of observations from multiple datasets, even observations at the same time t. It is also notable that A t varies over time, reflecting the fact the variables that are sampled will change over time. To complete the relationship between the mean state ( x t ) and the direct observations ( y t ), we include observation error, v t ∼ N ( 0 , Σ v , t ) , that follows a multivariate normal distribution with variance-covariance matrix Σ _v,t .

The persistence equation allows the carbonate state to have a memory (or be auto-correlated) in time. It is represented as a random walk process Σ v , t

Specifically, this means that the new carbonate state ( x t ) follows the previous state ( x t − 1 ) but with some random perturbation or persistence error added. The persistence error ( ( w t ∼ N ( 0 , Σ w ) ) follows a multivariate normal distribution with variance-covariance matrix Σ w . The random walk, a widely used formulation in other statistical models (e.g. Markov Chain Monte Carlo), was chosen as the simplest model possible to achieve a non-parametric temporal persistence of the state that adapts to available observations, without the mean-reverting quantities of an auto-regressive model. The degree to which the state can change from one time to the next is controlled by the magnitude of the diagonal elements of Σ w . This covariance matrix was structured to also allow for correlations between adjacent depth-layers, but not between the different carbonate variables (so it has a block-diagonal structure). Unlike the observation error covariance matrix Σ v , t which changes over time based on which observations are available, the persistence error covariance matrix Σ w is fixed in time.

The state space model is defined by Equations 1 and 2, both of which provide information for estimating the carbonate state through time. This linear, Gaussian framework is robust and reliable, and its properties are well understood and interpretable. The Kalman filter/smoother algorithm provides the optimal solution for the state space model and yields the state of the carbonate system ( x t ). Specifically, it provides estimates for the monthly mean values of DIC, TA, pH, and fCO₂ at all depths, and each of their time-varying variances. The latter can then be used to produce error bars or confidence intervals. The algorithm uses Equations 1 and 2 in a sequential, or recursive, estimation procedure. This relies on the Kalman filter, which operates as a forward-in-time recursion such that at each time step t, an estimate of the new state vector (the four variables at the L depths) is available via a one-step ahead prediction using (2). This persistence estimate is then updated to be closer to any available observations. If there are no observations at time t, the observations Equation 1 is not used, only the persistence Equation 2. The Kalman smoother then further refines these estimates through a backwards-in-time recursion. The ratio of the variances for the observation error (diagonal elements of Σ v , t ) to the persistence error (diagonal elements of Σ w ) is a key quantity that dictates how closely the state estimates follow the observations. The smaller the Σ v , t : Σ w ratio the closer our estimates are drawn to the observations otherwise the more the estimate will ignore the observations and predict a constant unchanging line. This is discussed in more detail in the Supplementary Material , along with details of the Kalman filter/smoother algorithm.

The Σ v , t matrix was estimated based on the variance of observations, and the Σ w matrix used maximum likelihood estimation. For details of how the matrix terms were estimated, the reader is referred to the Supplementary Material .

The state of the carbonate system (i.e. the four variables DIC, TA, pH and fCO₂) was estimated monthly from 1993-2016 for the 20 depth-layers from 0 m - 3500 m using the state space model. These results are presented as time-depth plots for each variable ( Figures 2 - 5 ). The estimated time series for the surface layer and an intermediate depth-layer were plotted in full to contrast their differences due to oceanic variability and observation sparsity, while all time series are presented Supplementary Figures S10-S13 . Depth profiles of the calculated annual rates of increase ( Supplementary Figure S18 ) visualizes the numbers that are listed on the right side of Figures 2A - 5A . The results for each variable are presented in turn.

DIC increases with depth, as expected, and there is a general increase of DIC over time as shown by the depth-weighted column mean ( Figure 2 ). The time series of the near-surface waters ( Figure 2B ) have larger variability than the sub-surface waters ( Figure 2C ). The waters deeper than 100 m show linearly increasing trends of DIC at rates of 0.57 to 0.76 µmol kg ^{− 1} year ^{− 1}. However, examining the individual time series for each depth ( Supplementary Figure S10 ), for depths 200-800 m and for most below 1800 m the estimated DIC is unchanging in the 1990s and increasing after 2000. These deep water results are consistent with the anthropogenic increase presented in Boteler et al. (2023) and later shown by Olivarez et al. (2024) to be a delayed effect of the Pinatubo 1991 eruption. The intermediate waters however were unexpected, though a similar change in trend is seen in oxygen ( Supplementary Figure S9a ) which goes from a lack of trend in the 1990s to a decline in oxygen in the 2000s. Note also that the state space model is resistant to outliers with only a small divergent of the mean state towards the observation, as seen at the end of 2010 ( Figure 2C ). A list of cruises that contributed to outliers in monthly averaged DIC observations is provided in the Supplementary Material .

At the near surface (top 50 m), there is a small seasonal cycle in DIC with lower values in the summer due to biological production and a shallow mixed layer. There is also a pronounced 4-year cycle, which could be an artifact of the data, be actual inter-annual variability, or be a combination. The 4-year cycle is also seen in the near-surface temperature and oxygen time series ( Supplementary Figures S7b, S9b ). It could be related to biological activity as well as related to the North Atlantic Oscillation (NAO) and Atlantic Multidecadal Oscillation (AMO). The NAO and AMO have 4 year and 4-5 year periodic signals, respectively, as identified in a spectral analysis and F-test of the frequencies (Thomson, 1982). By visually comparing them, we see that high DIC values (years 2000 and 2004) coincide with moments of high NAO index. The other years of high DIC (years 1996, 2009, and 2014) are just a year out of phase with both the NAO and a high AMO index.

fCO₂ shows that it is steadily increasing over time at rates between 0.70 and 2.45 µatm year ^{− 1} ( Figure 3 ), which we generally attribute to the atmospheric CO₂ increase of ∼1.93 µatm year ^{− 1} (Lebehot et al., 2019), with the near-surface fCO₂ increase being quite close to that of atmospheric CO₂. The time-depth plot also shows that fCO₂ increases with depth as a consequence of the increase in pressure, as expected. At the near-surface, fCO₂ is increasing at a rate of 1.48±0.04 µatm year ^{− 1}, which agrees with other ocean observation-based estimates (1.47 ± 0.06 µatm year ^{− 1}) and is close to the model-based estimates (1.9 ± 0.09 µatm year ^{− 1}) (Lebehot et al., 2019). It is however interesting to note that at depths of 1000-2000 m, the rate of fCO₂ increase is faster than the atmosphere, between 1.72 and 2.2 µatm year ^{− 1}. It is likely the deep convection in the Labrador Sea is taking the winterfCO₂ down to the intermediate waters.

The fCO₂ near-surface time series in Figure 3B shows the presence of direct observations from both GLODAP and SOCAT, as well as pseudo-observations. Between 2005 and 2012, the high frequency of SOCAT observations reveals a seasonal cycle peaking around February ( Supplementary Figure S15 ), which is similar to other pCO₂ products in the Labrador Sea (Takahashi et al., 2002; Rödenbeck et al., 2013; Zeng et al., 2014; Arruda et al., 2024). The increasing trend in the near-surface is more obvious for fCO₂ than for DIC as the fCO₂ seasonal cycle presented here has an amplitude smaller than the increasing trend. This amplitude is however smaller than expected but it is based on the current data and the uncertainty parameters chosen. The monthly fCO₂ data between 2005 and 2012 have such a large scatter that a larger observation error would be preferred, but in that case the mean state estimate resulted in a linear line passing through them. We reduced the observation error to show that the seasonal cycle can be shown to be at the correct phase, though with an underestimated amplitude.

The width of the confidence interval shows the uncertainty of the mean state, but note that this uncertainty is not necessarily equal to the spread of the data. Figure 3B and Figure 4B show confidence intervals that are narrower than the spread of the data. The width of the confidence interval is based on three factors: frequency of data, observation error, and persistence error. The persistence error was estimated using maximum likelihood. The first two are based on the data and are illustrated in Figure 3B for the confidence interval narrow between 2005 and 2012 due to the increased frequency of SOCAT data and the small observation error for SOCAT within Σ v , t (see Supplementary Material ).

Of the four variables analyzed, fCO₂ sub-surface benefited the most from the use of pseudoobservations, and without them, the fCO₂ time series estimates would not be possible. In the deeper waters ( Figure 3C ) GLODAPv2 observations are available only one time, and the rest are pseudo-observations that were calculated using CO2SYS.

In the Labrador Sea, ocean acidification (i.e. declining pH) is occurring at rates between −0.0018 and −0.0007 year⁻¹ ( Figure 4 ). In the near-surface (top 20 m) we see higher mean pH values (pH > 8.14), but of the lowest rate of ocean acidification. This rate is smaller than other estimates for the region and time frame (−0.0014 year⁻¹ Chau et al. (2024) and −0.002 year⁻¹ Lauvset et al. (2015)) and may be due to our data which has a larger variance than change over time. Through the water column, the pH value lowers with depth, but the rate of ocean acidification does not change linearly with depth ( Supplementary Figure S18 ). The rate of ocean acidification is stronger for depths 50 to 2000 m. Below 2000 m around the aragonite saturation horizon (Ω = 1; solid line), the rate of ocean acidification is slower than the waters above or below it.

Total Alkalinity deeper than 50 m has minimal variability, but above 50 m lowers after 2006 ( Figure 5A ), which was found to be correlated with a decrease in salinity. This strong correlation between TA and salinity allowed us to calculate pseudo-observations of TA, but it also makes interpretation of TA difficult since its changes are related to salinity. Hence we use salinity-corrected alkalinity to aid in interpretation. Salinity normalized total alkalinity (nTA) (Friis et al., 2003) removes the expected variability of salinity, allowing us to distinguish changes in alkalinity that are not related to changes in salinity. After salinity-normalization, the trends in the top 50 m have been reduced by more than half, and in the top 20 m the trend is no longer statistically significant, indicating the trend was due to the changes in salinity ( Figure 5B ). Overall, nTA fluctuates less than ±10 µmol kg⁻¹ about the mean, which typically is considered a negligible variation. There are a few high nTA values (nTA >2315 µmol kg⁻¹) with most occurring around 2010 at 46°W and 52°N. Specifically, there is a single TA observation in May 9, 2010 (https://doi.org/10.25921/zp3g-cm29; doi: 10.5285/cf2d9ba9-d51d-3b7c-e053-8486abc0f5fd) that did not appear as an outlier when looking at the original GLODAP data as a time series (i.e. it was within the acceptable range of all the observations, Supplementary Figures S2–S5 ). However, by being a TA = 2315.8 µmol kg⁻¹ and no other observations in that month to balance it out, after monthly averaging the data it stood out as a potential outlier. There was enough variability in the surrounding years that we left this observation in, so as to not force the data into perfection. However, we note that this high TA value has cascading effects on other calculated quantities, as will be seen in the Revelle factor.

Using the estimated state for DIC and TA with CO2SYS, we calculated further derived quantities of the carbonate system: the aragonite and calcite saturation state, the Revelle factor and the speciation of DIC ( Figure 6 ). We see in Figure 6E that for every depth-layer the aragonite saturation state is steadily decreasing. At the near-surface there is also some seasonality that is correlated with that of DIC. The same time-depth characteristics are also seen for calcite saturation. The calcite saturation has a range of 1.18 to 3.81, thus the saturation horizon (when Ω = 1) is not shown. The aragonite saturation horizon (Ω = 1), as shown on previous plots ( Figures 2 - 5 ), is seen to go from our 16th depth-layer (2590-2810 m) to our 15th (2370-2590 m), with this shoaling between depth-layers occurring around 2003. Considering these depth-layers span 220 m each, the shoaling could be up to a max width of 440 m over 24 years, relating to an increase of 18 m year ^{− 1}.

The Revelle factor was calculated from the TA and DIC results, γ D I C = 3 T A · D I C − 2 D I C 2 ( 2 D I C − T A ) ( T A − D I C ) (Sarmiento and Gruber, 2006). For all depths after 2005, the Revelle factor increases, indicating that fCO₂ becomes more sensitive over time to changes in DIC. This increase is paused in 2010, but this is related to the anomalous high TA value. Time series of the top few layers ( Supplementary Figure S14 ) show in more detail these changes in the Revelle factor. For the near-surface this implies a lowering of the ocean’s capacity to take up more atmospheric carbon. Since the Revelle factor is calculated from TA, this increased sensitivity is a reflection of the freshening in the area ( Supplementary Figure S8 ). For the interior ocean, the Revelle factor has been used to determine the concentration of anthropogenic DIC (Terhaar et al., 2022; Sabine et al., 2004).

Our state space model also calculates 95% confidence intervals for all the carbonate variables. Hereafter, confidence intervals are referred to as uncertainties, though sometimes the related 95% margin of errors are reported instead (distance from mean estimate to upper uncertainty bound) as in Figure 7 . The uncertainties produced from our time-correlated method are lower than the uncertainties from using CO2SYS on DIC and TA values which treat the carbonate variables as being independent through time, rather than auto-correlated (see Table 1 for a detailed comparison). Overall, our time-correlated method has lowered the uncertainty by more than 64%. Our uncertainty estimates are driven by the data distribution over time, as dictated by the opportunistic sampling protocol. In other words, the uncertainties varied over time based on the presence/absence of observations ( Figure 7 ) and the average variability of observations within that depth-layer. In the near-surface waters of DIC, TA and pH the uncertainty is higher due to larger variability of the data. For fCO₂ at the surface, between 2005 and 2012 the uncertainty is reduced due to the presence of SOCAT data. The fCO₂ between 700-1600 m also has higher uncertainty than the waters above or below due to the pseudo-observations having larger variability.

We developed and implemented a multivariate, time series method to estimate the depth-time variability of the carbonate system in the Labrador Sea. This incorporated the four carbonate variables: DIC, TA, pH and fCO₂. Our method takes into account the temporal memory of the carbonate system, and as a consequence is able to considerably lower the uncertainty of pH by 71% and by 64% for fCO₂, compared to methods that assume independence in time. The estimated time series of DIC and TA were also used to calculate other aspects of the carbonate system, such as the aragonite and calcite saturation state, the Revelle factor and the speciation of DIC (i.e. carbonic acid, bicarbonate and carbonate). Our results for the carbonate system quantify their variations with depth and highlight temporal trends that are assumed to be due to anthropogenic sources: increases in DIC (0.57 to 1.16 µmol kg ^{− 1} year ^{− 1}) and fCO₂ (0.70 to 2.45 µatm year ^{− 1}), as well as the acidification of pH (-0.0018 to -0.0009 year ^{− 1}).

We used a statistical time series method that amalgamates multiple data sources and takes account of their unique features. As well as direct observations from GLODAP and SOCAT, CO2SYS was used with DIC and TA observations to calculate pseudo-observations. The pseudo-observations more than doubled the amount of information on pH and fCO₂, which allowed for more detailed and realistic estimates of their depth-time variability, including the estimation of anthropogenic trends of pH and fCO₂. Using pseudo-observations also made pH and fCO₂ estimates chemically balanced with the rest of the carbonate system. DIC and TA were chosen to be in our CO2SYS calculation as we prioritized variables that were abundant over time, but as ARGO data continues to increase, pH could be considered instead of TA as the second input variable for CO2SYS. Using an input combination of one T, P-dependent (pH or fCO₂) and one non-dependent (DIC or TA) variable provides lower bias uncertainty for CO2SYS outputs (Raimondi et al., 2019). The distribution of data through time is also considered in the estimation of the uncertainties. Our time-varying uncertainties are smaller in the presence of observations, and larger in the data gaps. This is seen in Figure 7 where there are vertical bands of lower uncertainty at times when observations are present. Our results were on monthly time steps, but any interval could be used that is deemed reasonable given the observation distribution. If the gaps between observations are large (i.e. more than a year) then the state estimate will have flat unchanging sections between observations due to the Kalman filter based estimation procedure.

The mean and variance of the carbonate state are not constant over time due to anthropogenic sources of CO₂. DIC observations analyzed using a 1-year sliding window showed that not only is the mean state increasing over time and has inter-annual variability, but the standard deviation (or standard error) of those observations is also increasing over time in the near-surface and shows inter-annual variability for multiple depths ( Supplementary Figure S17 ). This behavior could provide insight into the changes in the DIC seasonal cycle reported by Landschützer et al. (2018). The change of DIC standard error over time could be modelled by having our observation error ( Σ v , t ) change over time. Improvements can be made in the estimation of the observation error parameter ( Σ v , t ), a key quantity in the analysis.

While the Kalman Smoother is resistant to outliers (i.e. the state estimate does only a blimp towards them and does not go directly to them, Figure 2C ), we found it was necessary to remove some extreme outliers. The motivation was because the small blip in our DIC estimates due to an outlier in 2015 around 700 m had a larger impact on our calculated quantities of the Revelle factor, bicarbonate and carbonate, which showed a local maximum at this depth. While no other outliers stood out after the monthly averaging, other suspected outliers in the original GLODAP data were removed as described in the Supplementary Material along with a list of the data and its cruise information.

Spatio-temporal estimates of the carbonate system have been produced by other methods, ranging from machine learning (Bittig et al., 2018) to Earth System Models (Carroll et al., 2020). The leading observation-based method for the estimation of the carbonate system is CANYON-B (Bittig et al., 2018). It uses a Bayesian neural network approach that employs non-linear relationships to estimate monthly spatial climatologies of the four carbonate variables and nutrients. Machine Learning methods are dependent on the data they are trained on. CANYON-B may find limitations in the training of their model in regions like the Labrador Sea, where the data has been mostly collected at the same time each year (e.g. AR7W from Labrador to Greenland each May). The CANYON-B framework employs year-to-year correlation whereas our method uses more localized time correlation (month-to-month). For validation of results, it is common to train the model on one data source (like GLODAP) and test on another (i.e. SOCAT or ARGO) (Bittig et al., 2018). Meanwhile, ECCO-Darwin (Carroll et al., 2020) used multiple data sources (i.e. GLODAP, SOCAT, ARGO) as observational constraints for their Earth System Model. Iida et al. (2020) also used multiple data sources (GLODAP, SOCAT, Satelite Data) with multiple linear regression and chemical relationships to estimate the carbonate system. Using multiple data inputs was a feature we also implemented in our observation-centric framework.

Our results yield time-varying uncertainties that are on average 64-71% lower compared to time independent methods such as calculations with CO2SYS on independent observations (see Table 1 ). There are multiple metrics for uncertainties. The time-varying uncertainties are defined as 95% confidence intervals that are calculated as part of the statistical model using the Kalman smoother algorithm. The time-independent CO2SYS calculates uncertainties using Gaussian propagation (Orr et al., 2018). From our method for the years 1993-2016, the carbonate system variables have the following average uncertainties: 5.14 µmol kg ^{− 1} for DIC, 4.51 µmol kg ^{− 1} for TA, 0.005 for pH and 8.36 µatm for fCO₂. When CO2SYS is used on time-independent observations, the estimated pH has an average uncertainty of 0.03 and fCO₂ uncertainty of 24 µatm. Carter et al. (2021) estimated spatial climatologies for the carbonate system using a mixed estimation of linear regression method (i.e. LIRv3) and neural networks, and in the North Atlantic their uncertainties averaged: 7.7 µmol kg ^{− 1} for DIC, 5.0 µmol kg ^{− 1} for TA, and 0.009 for pH. The neural network based Canyon-B uncertainty estimates averaged at: 8.6 µmol kg ^{− 1} for DIC, 5.7 µmol kg ^{− 1} for TA, and 0.008 for pH. Compared to our uncertainties, Carter et al. (2021) and Bittig et al. (2018) are slightly larger for DIC and lower for TA, but overall they are comparable values.

Our pH in the top 20 m shows the acidification in the Labrador Sea at a rate of 0.0009 year ^{− 1}, which is almost half the rate of other studies of the Labrador Sea of -0.0014 to -0.0016 year ^{− 1} (Chau et al., 2024), meanwhile our sub-surface waters do have more comparable rates of -0.0015 to -0.0018 year ^{− 1}. With the large variance of surface pH observations, the signal-to-noise ratio is much higher in the surface waters than interior ocean ( Figures 4B, C ). This large observation variance paired with sporadic pH GLODAP observations, means our results are showing a steady decline as opposed to a seasonal cycle, as reported in the neighboring Irminger Sea, Bates et al. (2014) that peaks in the summer and spans 8.07 to 8.15 for years 1982-2006. The pseudo-observations we calculate for in the surface layer have larger uncertainty than the direct observations, which already have large variability themselves. This is an issue for pH and the fCO₂ seasonal cycle has a smaller amplitude than expected. A limitation of using pseudo-observations in the surface is that they add more uncertainty to a system, making it harder to distinguish process-driven variability like seasonal cycles from natural noise of observation; a limitation that other methods overcome by using global climatologies within their analysis. While the pseudo-observations have their surface limitations, their use in sub-surface hints at inter-annual variability, with a slow acidification between 2004 and 2010 before resuming ( Supplementary Figure S12 ), though further analysis with more data should be performed.

Our estimated time series for the carbonate system in the Labrador Sea produces comparable results to others (Chau et al., 2024; Feely et al., 2023; Gregor and Gruber, 2021). The average across time of our near-surface layer is summarized in Table 2 . The other studies ran global analyses but by extracting the Labrador Sea from their results, we find our DIC, TA, fCO₂, aragonite saturation and calcite saturation are close or overlapped with their values. Our average surface pH of 8.16 is on the high end of the other studies. These studies we compare with overlap in time span, and overlap with data inputs. Gregor and Gruber (2021) uses GLODAP TA, SOCAT pCO₂ and global climatologies to calculate the other variables using CO2SYS and Geospatial Random Cluster Ensemble Regression. Chau et al. (2024) uses SOCAT fCO₂, pCO₂ climatology and CMEMS data in CO2SYS and a neural network ensemble approach (CMEMS-LSCE-FFNN). Feely et al. (2023) uses SOCAT and DIC, TA and pH from GLODAP and the Common Online Data Analysis Platform in North America (CODAP-NA) and pools all into one observational database before it was interpolated using Data Interpolating Variational Analysis.

Our average TA results below 100 m were 2302 µmol kg ^{− 1}, which is the same as Broullón et al. (2019) reported for an annual mean TA in the Labrador Sea for X based on GLODAPv2.2019 data. In the surface layer, Broullón et al. (2019) reports a TA of 2300 µmol kg ^{− 1}, but our estimated TA is lower at 2278 µmol kg ^{− 1}. The variability seen in TA is correlated with salinity, as confirmed by the salinity-normalized TA ( Figure 5 ). The freshening of the Labrador Sea in this timeframe has not been shown before, actually it has been described as having no change while the eastern North Atlantic in the 2010s had a notable freshening and cooling event (Holliday et al., 2020). We see cooler surface waters in the Labrador Sea in the 2010s compared to 2000s and these cooler waters are suggested to play a role in this eastern freshening Fox et al. (2022). But returning to the Labrador Sea, the freshening we see as of 2006 comes earlier than expected, as an ocean simulation has estimated around 13 years for the subpolar North Atlantic to adjust to the increase of annual Greenland freshwater discharge which increased as of 2000 (Dukhovskoy et al., 2021). Ocean simulations have also shown that freshwater from the Arctic could cause freshening in the Labrador Sea (Zhang et al., 2021).

Using an Earth System model, Jiang et al. (2019) predicts surface pH changing from 8.1 to 7.75 for 2000 to 2100 Jiang et al. (2019), see their Figure 6 ), implying a declining pH of 0.0035 year ^{− 1}, 4 times larger than our estimate. For a surface Revelle factor in the year 2000, our results agree with (Jiang et al., 2019) that the Revelle factor is around 15 in the Labrador Sea. Our aragonite saturation horizon in 2004 is found at 2370 m, which is very close to the 2300 m reported by Azetsu-Scott et al. (2010) using 2003-2005 Labrador Sea data. Our max increase of aragonite saturation horizon of 440 m from 1993-2016 agrees with Tanhua et al. (2007) who reported a 400 m increase in the North Atlantic from pre-industrial to 2004 and projected a 700 m increase from pre-industrial to 2050. The aragonite saturation horizon is increasing at a faster rate in the North Atlantic than in the Pacific, which was reported to have increased at 1-2 m year ^{− 1} (Feely et al., 2012).

Carbonate (CO₃ ^–2) is rarely discussed in observation-based work. Using CO2SYS to calculate carbonate from our DIC and TA, our carbonate estimate for the year 1994 is 143 µmol kg ^{− 1} in the top 20 m, which is similar to Orr et al. (2005), whose biogeochemical model estimated an average surface carbonate to be around 160 µmol kg ^{− 1} at 60°N for the year 1994, and our 2010 carbonate estimate of 134 µmol kg ^{− 1} is similar to the 120-140 µmol kg ^{− 1} Jiang et al. (2023) estimated for the Labrador Sea. Our near-surface carbonate has a seasonal cycle but also a declining trend of -1.06 µmol kg ^{− 1} year ^{− 1} (intercept of 151.1), projected forward to the year 2100 would yield a carbonate value of 38 µmol kg ^{− 1}, which is close to the SSP3-7.0 projections of Jiang et al. (2023) (see their Figure 8).

Examining the individual time series of DIC for each depth ( Supplementary Figure S10 ) we see the effects of the Pinatubo 1991 eruption. Our estimated DIC below 1800 m is unchanging in the 1990s and increasing after 2000. These deep-water results are consistent with anthropogenic carbon results in Boteler et al. (2023) and has since been shown by Olivarez et al. (2024) that the Pinatubo 1991 eruption had a delayed effect of not showing an increase in ocean carbon until after the year 2000 (using pCFC-12). This temporal pattern of an unchanging 1990s and increase after 2000 was also reported for our DIC between 200-800 m, which may have a separate biological reason for our oxygen results at those depths are also showing a change in trend around the year 2000. A longer time series would help in de-tangling its processes for the various depths, as this multi-decadal variability might also be in the surface layers but hidden within the seasonal cycle. The near-surface DIC as well as oxygen also have a 4-year cycle that seems related to the North Atlantic Oscillation (NAO) and Atlantic Multidecadal Oscillation (AMO).

The fCO₂ in the near-surface was estimated to be increasing at a rate of 1.48 µatm year ^{− 1}. This is lower than the atmospheric CO₂ increase of ∼1.9 µatm year ^{− 1} but the same as other observation based estimates for the surface North Atlantic of 1.47 ±0.06 µatm year ^{− 1} for years 1992-2014 (Lebehot et al., 2019) and Atlantic increase of 1.67 ±0.02 µatm year ^{− 1} for years 1990-2018 (Gregor and Gruber, 2021). The near-surface fCO₂ values always remain below atmospheric partial pressure and with the ocean increasing at a slower rate, the gradient between atmospheric and ocean CO₂ is increasing. Moving to the deeper waters, sub-surface measurements of fCO₂ are rare, but they could prove useful for tracking the anthropogenic increase of CO₂ in the ocean. While DIC is a popular choice for estimating the anthropogenic change in the ocean interior, fCO₂ has a lower range of values and variability through the water column than DIC, making it more sensitive to anthropogenic changes through time and easier to identify variations within the ocean interior. Similar temporal variability to fCO₂ was also seen in the Revelle factor, which is described as the sensitivity of fCO₂ to the changes in DIC, and has been used to infer anthropogenic carbon in the ocean interior. The Revelle factor showed a steady rate of increase in the waters deeper than 100 m ( Figure 6D ). More sub-surface measurements of fCO₂ would be useful to further investigate the Revelle factor and fCO₂ and how their roles in quantifying anthropogenic change are distinctive or overlapping.

In summary, we present an observation-centric method to estimate the mean state carbonate system and its uncertainties for the Labrador Sea for the period of 1993-2016. The state space model framework allowed for the multivariate estimation of the carbonate system, and for the fusion of multiple data sources: GLODAP cruise data, SOCAT ships of opportunity, and calculated pseudo-observations via CO2SYS. Our analysis focused on estimating time-depth distributions and temporal trends at multiple depth layers through the water column in the Labrador Sea for carbonate system variables. Our analysis was transparent, making no assumptions or model constraints on the variables to let the data drive the estimation, the changes over time, the uncertainties, and thus the story of the carbonate system.