Ocean carbon sink assessment via temperature and salinity data assimilation into a global ocean biogeochemistry model

Global ocean biogeochemistry models are frequently used to derive a comprehensive estimate of the global ocean carbon uptake. These models are designed to represent the most important processes of the ocean carbon cycle, but the idealized process representation and uncertainties in the initialization of model variables lead to errors in their predictions. Here, observations of ocean physics (temperature and salinity) are assimilated into the ocean biogeochemistry model FESOM2.1-REcoM3 over the period 2010–2020 to study the effect on the air–sea carbon dioxide (CO₂) flux and other biogeochemical (BGC) variables. The assimilation nearly halves the model–observation differences in sea surface temperature and salinity, with modest effects on the modeled ecosystem and CO₂ fluxes. The main effects of the assimilation on the air–sea CO₂ flux occur on small scales in highly dynamic regions, which pose challenges to ocean models. Its largest imprint is in the Southern Ocean during winter. South of 50° S, winter CO₂ outgassing is reduced; thus the regional CO₂ uptake increases by 0.18 Pg C yr⁻¹ through the assimilation. Other particularly strong regional effects on the air–sea CO₂ flux are located in the area of the North Atlantic Current (NAC). However, the effect on the global ocean carbon uptake is a comparatively small increase by 0.05 Pg C yr⁻¹ induced by the assimilation, yielding a global mean uptake of 2.78 Pg C yr⁻¹ for the period 2010–2020.

1 Introduction

The ocean plays a pivotal role in regulating the global carbon budget and thereby mitigating the impacts of anthropogenic carbon dioxide (CO₂) emissions on the Earth's climate. Since the 1960s, the ocean has consistently absorbed around 25 % of anthropogenic CO₂ emissions annually (Friedlingstein et al., 2023) and has cumulatively taken up 26 %–34 % of fossil and land-use-change CO₂ emissions since the onset of the Industrial Revolution (Crisp et al., 2022). However, quantification of air–sea CO₂ flux still remains challenging. Air–sea CO₂ flux is usually inferred from the gradient of partial pressure (pCO₂) or fugacity of CO₂ across the air–sea interface (Wanninkhof, 2014). However, during 2010–2020, which constitutes the best-sampled decade in terms of surface ocean pCO₂ observations so far, observations covered merely 3 % of the monthly global ocean (as calculated from the 1° × 1° SOCAT product; Bakker et al., 2023). While the North Atlantic and North Pacific are comparably well observed, data remain scarce in vast regions, such as the Indian Ocean, the South Pacific, and areas south of 30° S during austral winter, where less than 1 % of SOCAT grid cells has been sampled. Although these observations are thought to be representative of a larger area (Jones et al., 2012; Hauck et al., 2020), challenges in deriving a comprehensive global estimate of the global ocean CO₂ uptake arise due to substantial spatial and temporal pCO₂ variations and potential biases induced by the irregular sampling pattern (Denvil-Sommer et al., 2021; Gloege et al., 2021; Hauck et al., 2023b). Particularly in the Southern Ocean, the uncertainty is considerable (Gerber et al., 2009; Gloege et al., 2021), and estimates of the mean flux range from −0.37 to −1.25 Pg C yr⁻¹ for the period 2010–2018 (data provided by Hauck et al., 2023b).

In the Global Carbon Budget, estimates of the ocean carbon sink were initially derived from hindcast simulations of global ocean biogeochemistry models (GOBMs) (Le Quéré et al., 2009; Wanninkhof et al., 2013; Hauck et al., 2020). More recently, air–sea CO₂ flux estimates were added based on regression and machine learning techniques, interpolating pCO₂ observations to achieve global coverage through advanced statistical methods (referred to as pCO₂ products; Rödenbeck et al., 2015). Furthermore, atmospheric transport models that ingest atmospheric CO₂ measurements were employed to estimate the ocean carbon uptake (referred to as atmospheric inversions; Peylin et al., 2013). Although the different estimation methods have provided valuable and robust insights into large-scale patterns of oceanic carbon uptake (Gruber et al., 2009), discrepancies have emerged. Assessments based on pCO₂ products tend to yield larger estimates of the ocean carbon sink, with stronger trends towards more uptake, compared to estimates based on models (Friedlingstein et al., 2023; Terhaar et al., 2022). The larger estimates are supported by ocean interior observations (Müller et al., 2023), atmospheric oxygen data, and atmospheric inversions (Friedlingstein et al., 2023). For the years 2010–2020, pCO₂ products included in the Global Carbon Project suggest a mean oceanic sink of 3±0.4 Pg C yr⁻¹, while the mean of Global Carbon Project GOBMs is 2.5±0.4 Pg C yr⁻¹ (data provided by Friedlingstein et al., 2023). Trends over the same time period are 0.7 and 0.3 Pg C yr⁻¹ per decade, respectively.

Machine learning estimates perform well when trained with sufficient data (Gloege et al., 2021). Their performance is less reliable in data-sparse areas. Particularly in the Southern Ocean, many pCO₂ products show diverging results from one another and are likely biased towards more ocean uptake (Hauck et al., 2023b). However even in parts of the North Pacific, which is undersampled in the 2010s, some pCO₂ products show spurious decadal trends (Mayot et al., 2024). Models provide process-driven estimates of the CO₂ flux across the entire global ocean, drawing from the theory of ocean dynamics and biological and chemical processes (Hauck et al., 2020; Fennel et al., 2022). Despite the growing confidence in our mechanistic understanding of the ocean carbon cycle (Crisp et al., 2022), models are also subject to uncertainty. This uncertainty stems from uncertainties in model parametrization, model spin-up, and initial conditions; unresolved sub-grid-scale processes; and uncertainties in the atmospheric forcing (Hauck et al., 2020; Terhaar et al., 2024).

Data assimilation (DA) can be employed to address the emerging discrepancies between pCO₂ products and models (Carroll et al., 2020). Several studies assimilating ocean surface pCO₂ have focused on specific regions (e.g., a baseline state of air–sea CO₂ fluxes in the Southern Ocean; Verdy and Mazloff, 2017), short time periods (e.g., optimized biogeochemical (BGC) initial fields for the period 2009–2011 in Brix et al., 2015), or the climatological mean state (e.g., corrections of large-scale pCO₂ model biases in While et al., 2012). These studies capture the assimilated pCO₂ observations well, while obeying physical laws and BGC equations. Data assimilation can also be used to provide a better understanding of various components of the ocean carbon cycle, such as the transport of anthropogenic CO₂ in the ocean (e.g., a reconstruction of anthropogenic carbon storage since 1770 in Gerber et al., 2009), regional and interannual variability in the air–sea CO₂ flux (e.g., global reanalysis in Ford and Barciela, 2017; Carroll et al., 2020; Valsala and Maksyutov, 2010), the biological carbon pump (e.g., carbon export at a nutrient-rich and nutrient-poor site and estimation of BGC parameters related to air–sea CO₂ fluxes; Sursham, 2018; Hemmings et al., 2008), and specific ecosystems (e.g., the northwest European shelf ecosystem in Ciavatta et al., 2016, 2018). So far, however, there is no data assimilation product that provides a long-term, annually updated estimate of global ocean CO₂ uptake.

While previous studies indicate that the available BGC observations, when assimilated in isolation, are too sparse to constrain the modeled carbon cycle (Verdy and Mazloff, 2017; Spring et al., 2021), the assimilation of physical variables is expected to have a significant indirect effect on the modeled air–sea CO₂ fluxes (Bernardello et al., 2024). This is because the uptake of atmospheric CO₂ depends in large part on the physical carbon transport between the surface, the mixed layer, and the deep ocean in the form of dissolved inorganic carbon (DIC) through mixing, upwelling, and subduction (Doney et al., 2004). According to current knowledge, ocean physics is the dominant driver of interannual variability in the global air–sea CO₂ flux and is also responsible for stagnation and acceleration of the CO₂ uptake on decadal scales (Doney et al., 2009; Keppler and Landschützer, 2019; Mayot et al., 2023; Liao et al., 2020; DeVries et al., 2017). Related to the strong control that physics exerts on the interannual variability in air–sea CO₂ fluxes, it was shown in one idealized study that assimilating ocean physics at the initial state of a model simulation has a stronger and more positive impact on the modeled carbon cycle on interannual timescales than assimilating the BGC initial state (Fransner et al., 2020). However, the relative importance of uncertainties in physical and biogeochemical fields generally remains an open research question (e.g., Séférian et al., 2014; Li et al., 2016; Lebehot et al., 2019). Therefore, we use ensemble-based data assimilation of physical observations into a global ocean general circulation model coupled to a biogeochemistry model, aiming to improve the modeled air–sea CO₂ flux for the years 2010–2020. For this, we continuously assimilate temperature and salinity observations from remote sensing at the surface and from in situ profile measurements for 11 years and update the modeled temperature, salinity, horizontal velocities, and sea surface height using an ensemble Kalman filter variant (Nerger et al., 2012).

Several difficulties are associated with physics DA into GOBMs. A common issue is erroneous equatorial upwelling leading to unrealistically high biological productivity in the tropics (Park et al., 2018; Gasparin et al., 2021; Raghukumar et al., 2015). Furthermore, any coupled ecosystem model is adapted to its associated physical model with its strengths and weaknesses through carefully selected parameter values and a spin-up to near-equilibrium. Accordingly, the modeled carbon cycle may be sensitive to deviations from the physical state that is typical for this model (Kriest et al., 2020; Spring et al., 2021). Potentially, this leads to biases in the carbon cycle through physics DA. Such effects highlight where physical model errors are compensated for by BGC parameters, and DA may thereby reveal critical areas for potentially unrealistic BGC model behavior in projections in a changing climate (Löptien and Dietze, 2019). The question therefore arises to what extent an ecosystem model coupled to a data-assimilated physical model also represents a more realistic biogeochemistry and which mechanisms drive the response of the CO₂ flux in physics DA approaches. One possible driver is the physical transport of DIC and alkalinity because velocities and diffusion are changed by the DA, affecting in particular the upwelling of carbon-rich waters and subduction, which are important to capture the ocean storage of anthropogenic carbon (Davila et al., 2022). Furthermore, physics DA may change pCO₂ directly through its temperature dependence, an effect emphasized by Verdy and Mazloff (2017). Additionally, the modeled biological pump might be altered, for example, through the temperature dependency of phytoplankton growth or through effects of stratification on nutrient availability.

In this study, we describe the response of the model's air–sea CO₂ flux to physics DA and identify the underlying mechanistic drivers. To this end, we differentiate between the thermally, DIC-, and alkalinity-induced components and changes in lateral and vertical transport through mixing and advection. We focus, firstly, on the global air–sea CO₂ flux. Secondly, we investigate the Southern Ocean given the relevant impact of DA in Southern Ocean winter in our study. Thirdly, we present regions in the North Atlantic given observational coverage and relevant local processes there.

2 Methods

2.1 Model FESOM2.1-REcoM3

The oceanic model component, FESOM2.1, computes the advection and diffusion of passive biogeochemical tracers. The model is based on hydrostatic primitive equations under the Boussinesq approximation and utilizes a finite-volume discretization approach with surface triangles projected vertically to form prisms. Surface salinity (S), temperature (T), and biogeochemical tracers are located at the vertices of triangles (nodes), while the horizontal velocities are centered at the triangles (elements). The model allows a variable mesh resolution (see Sect. 2.2) and incorporates parametrizations for diffusion and eddy-stirring along isoneutral surfaces, for which parametrized mixing is scaled by mesh resolution (Danilov et al., 2017). Vertical mixing is parametrized through the K-profile parametrization (KPP) scheme, and the mixing depth is specified through a “boundary layer” (the layer of active mixing, which may have a vertical structure because the mixing of all properties across the layer is not instantaneous, as opposed to the mixed layer which is defined by already well mixed properties; Large et al., 1994), with an additional vertical-mixing scheme used in the Southern Ocean (Monin–Obukhov parametrization; Timmermann and Beckmann, 2004). The sea surface salinity (SSS) is restored towards the World Ocean Atlas climatology through a fictional surface flux with v_SSS=0.17 m d⁻¹ according to Eq. (1) and as in Gürses et al. (2023):

$$\begin{array}{}\text{(1)}& ({\mathrm{SSS}}_{\mathrm{clim}}-{\mathrm{SSS}}_{\mathrm{model}})\times v_{\mathrm{SSS}}\times (h_{\mathrm{surf}}{)}^{-\mathrm{1}},\end{array}$$

with surface-layer thickness h_surf. A detailed description of FESOM2.1 and a model assessment are provided by Danilov et al. (2017) and Scholz et al. (2019, 2022).

The ocean biogeochemistry component, the Regulated Ecosystem Model version 3 (REcoM3), describes processes in the ocean carbon cycle and represents oceanic carbon in the form of dissolved inorganic carbon, dissolved organic carbon, plankton, and detritus (Gürses et al., 2023). REcoM3 contains 28 BGC tracers listed in Table A1. There are two phytoplankton groups (diatoms and small phytoplankton with implicit representation of calcifier), two zooplankton groups (mixed and polar macro-zooplankton) (Karakuş et al., 2021), and two classes of detritus. REcoM3 includes variable intracellular stoichiometry with ratios of C : N : Chl : CaCO₃ for the small phytoplankton and C : N : Chl : Si for diatoms, which are propagated to zooplankton and detritus (Schartau et al., 2007; Hohn, 2008). The publicly available routines to model the ocean carbonate system (mocsy2.0; Orr and Epitalon, 2015) are used to compute pCO₂ and air–sea CO₂ flux, employing the gas-exchange parametrization of Wanninkhof (2014). Alkalinity is restored by a fictional surface flux of 10 m yr⁻¹ (as in Hauck et al., 2013; Schourup-Kristensen et al., 2014; Gürses et al., 2023). The current model version FESOM2.1-REcoM3 was assessed by Gürses et al. (2023), and previous versions were evaluated and applied in global and regional studies of the ocean carbon cycle and planktonic ecosystems (Hauck et al., 2013; Schourup-Kristensen et al., 2014; Hauck et al., 2020; Karakuş et al., 2021).

2.2 Simulation setup

The model setup for both simulations closely follows Gürses et al. (2023). The mesh resolution is nominally 1°, ranging between 120 and 20 km, with enhanced resolution in the equatorial belt and north of 50° N (126 858 surface nodes). It has 47 vertical layers with thicknesses ranging from 5 m at the surface to 250 m in the deep ocean, as described by Scholz et al. (2019, CORE mesh). The model time step is 45 min. For atmospheric forcing, JRA55‐do v.1.5.0 is used, a reanalysis product tailored for driving ocean–sea–ice models (Tsujino et al., 2018a). The atmospheric CO₂ mixing ratio values were taken from the Global Carbon Budget (Joos and Spahni, 2008; Ballantyne et al., 2012; Friedlingstein et al., 2023). We use model restart fields from Gürses et al. (2023), where the model was spun up by repeating the 1961 JRA forcing for 189 years with preindustrial atmospheric CO₂ conditions, followed by a period from 1800 to 1957 with increasing atmospheric CO₂. Subsequently, simulations were continued with historical JRA forcing from 1958 to 2009. During the assimilation window (2010–2020), we conduct two ensemble simulations to study the impact of data assimilation (DA): one without DA (referred to as FREE) and another identical setup applying DA (referred to as ASML). For each simulation, the ensemble mean for the following variables is written as output: temperature, salinity, velocity, boundary-layer depth, surface pCO₂, DIC, alkalinity, nutrients, chlorophyll, net primary production, and biological export through sinking of detritus at 190 m. For the year 2020, additional output is available for individual ensemble members, mixed-layer depth, physical sources or sinks of DIC and alkalinity through horizontal and vertical advection and diffusion, and biological net sources or sinks of DIC and alkalinity through combined processes. For DIC, the net biological term is the sum of photosynthesis, respiration, remineralization of dissolved organic carbon, and formation and dissolution of calcite (Gürses et al., 2023, Eq. A6). For alkalinity, the net biological term is the sum of nitrogen assimilation and remineralization and the formation and dissolution of calcite (Gürses et al., 2023, Eq. A7).

2.3 Data assimilation

2.3.1 Assimilated observations

The assimilated observations are sea surface temperature (SST), sea surface salinity, and profiles of temperature and salinity. SST observations are from the Operational Sea Surface Temperature and Ice Analysis (OSTIA) dataset (Good et al., 2020). OSTIA provides daily gap-free maps of SST at a horizontal resolution of 0.05° × 0.05°, compiled from in situ and satellite data from infrared and microwave radiometers. The OSTIA observations were averaged to the FESOM2.1 model grid because their spatial resolution is higher than the nominal resolution of the model grid. An observation error standard deviation of 0.8 °C is prescribed for the DA following Nerger et al. (2020). Observations are excluded in the DA process if the difference between the model and observation exceeds 3 times the observation error standard deviation, thus 2.4 °C, and at grid points with sea ice in the model, as in Tang et al. (2020) and Mu et al. (2022). This exclusion keeps the model stable despite large differences between the model and observations at these sites, particularly as water temperature and salinity develop differently under sea ice than under the influence of the atmosphere (Tang et al., 2020). Instead, a “gentler” correction is made by assimilating neighboring points. After the initial phase, about 7 % of SST observations are excluded because of the 2.4 °C threshold. Nevertheless, the data assimilation still has a strong effect in areas where these large model–observation discrepancies are typically found (the North Atlantic, Japan, and the Southern Ocean).

The assimilated SSS data are from the European Space Agency (ESA) Sea Surface Salinity Climate Change Initiative (CCI) v03.21 dataset (Boutin et al., 2021). ESA-CCI contains daily data at a spatial resolution of 50 km, albeit not capturing temporal variability below weekly. The ESA-CCI observations are averaged to the FESOM2.1 model grid. We prescribe a constant observation error standard deviation of 0.5 psu following Nerger et al. (2024). Like for the SST data, SSS observations are excluded at locations where sea ice is present in the model.

The assimilated temperature and salinity profiles are from the EN.4.2.2 dataset (Good et al., 2013). The EN4 dataset contains quality-controlled profiles from various in situ ocean profiling instruments. To assimilate the profiles, the observations are assigned to the respective model layers (depth range) in the vertical. In the horizontal, the model values are computed as the average of the grid points of the triangle enclosing the observation. The observation error standard deviation is set to 0.8 °C for temperature and to 0.5 psu for salinity without excluding observations, as in Tang et al. (2020).

2.3.2 Assimilation method and implementation

For the assimilation, we use the Localized Error-Subspace Transform Kalman Filter (LESTKF; Nerger et al., 2012). The LESTKF sequentially updates the model forecast, incorporating observations when and where available. The model state and error covariance are represented by an ensemble simulation. Thereby, the assimilation of temperature and salinity affects the state of the physical model in its whole, including the horizontal velocities and sea surface height. A review of the LESTKF and other filters frequently used in geophysics can be found in Vetra-Carvalho et al. (2018). The assimilation is implemented using the Parallel Data Assimilation Framework (PDAF version 2.1), a software environment for data assimilation. PDAF is an open-source project and provides fully implemented DA algorithms (Nerger et al., 2020; http://pdaf.awi.de, last access: 28 January 2025). The current implementation builds on the works of Mu et al. (2022), who used DA of ocean temperature and salinity for sea-ice forecasts with FESOM2.0 coupled to an atmospheric model, and Tang et al. (2020), who studied the dynamic impact of oceanic DA into FESOM1.4 onto a coupled atmospheric component.

With localization of the LESTKF, the observation error is increased for an increasing horizontal distance between an observation and a model grid point, which weighs down the influence of a more distant observation. This avoids the model being influenced by observations at distant locations through spurious ensemble estimated correlations. We use a localization radius of 200 km and choose a fifth-order polynomial weighting function that mimics a Gaussian function (Gaspari and Cohn, 1999). We apply daily analysis steps at 00:00 UTC model time, assimilating all available observations for the day. The DA process only directly updates the physical model variables temperature, salinity, horizontal velocities, and sea surface height. After each assimilation step, corrections are applied to the analysis state to ensure the consistency of model physics: salinity is set to a minimum value of zero, and temperature is set to a minimum value of −2 °C, if the value is otherwise below. The increment of sea surface height (SSH) is limited to 2 standard deviations of the ensemble. While in the simulation the correction was necessary for about 10 % of SSH updates and 10⁻⁵ % of temperature values, the correction of salinity was never required. The analysis step is followed by an ensemble forecast of 1 d.

The ensemble size is 40, a compromise to balance computational resources while ensuring a sufficiently large ensemble with enough variability even in the deep ocean. The ensemble is generated through an initial perturbation of sea surface height, horizontal and vertical velocities, temperature, salinity, and sea-ice concentration based on the implementation of Tang et al. (2020). This initial ensemble perturbation is generated by second-order exact sampling (Pham, 2001) from a model trajectory of FESOM2.1. With this method, the leading empirical orthogonal functions (EOFs) of a model trajectory are used to generate an ensemble perturbation that contains the leading patterns of model variability. A timescale must be chosen for the variability that is represented by the ensemble. Here, we chose variability on a weekly timescale (Tang et al., 2020).

To maintain ensemble spread, we apply a perturbed atmospheric forcing with an autoregressive perturbation (perturb_m,n) at every model time step (n) to each ensemble member (m):

$$\begin{array}{}\text{(2)}& {\mathrm{perturb}}_{m,n+\mathrm{1}}=(\mathrm{1}-\mathrm{arc})\cdot {\mathrm{perturb}}_{m,n}+\mathrm{arc}\cdot s\cdot {\mathrm{rand}}_m,\end{array}$$

where rand_m is a stochastic element, again generated by second-order exact sampling from a 72 d long trajectory of atmospheric forcing fields that captures patterns of day-to-day atmospheric variability. The autoregression coefficient (arc) can be used to tune how quickly the perturbation changes and is set to the inverse number of model steps per day. s is a scaling factor for each perturbed atmospheric forcing field. For specific humidity, downwelling longwave radiation and air temperature s=10 is used. The perturbation of winds is set to the smaller value s=2 because the air–sea CO₂ flux in the model is particularly sensitive to perturbations of the wind fields. Due to the functioning of the Kalman filter (which updates the model error covariance in each analysis step to reflect the new reduced uncertainty), the ensemble spread decays at each analysis step. As the method relies on a sufficiently large ensemble spread, an inflation of the ensemble covariance is applied (Pham et al., 1998). Thereby, the ensemble covariance matrix is amplified by a factor of $\mathrm{1}/\mathit{\rho }$ before entering the updating step. This so-called forgetting factor downweighs that past observations have reduced the model uncertainty (see, e.g., Nerger et al., 2005). The forgetting factor is tuned to maintain model uncertainty, where ρ=1 means no inflation and smaller values mean larger inflation. Here, we use a time-varying forgetting factor between ρ=0.95 and ρ=1. The strongest inflation (ρ=0.95) is applied during the first 2 weeks of the DA process. This is when the DA increments are largest because the model state estimates are furthest from the observations. During the following 75 d, ρ is increased to 0.99. From month 17 onward, the forgetting factor is set to either 0.99 or 1.0 depending on the ensemble standard deviation of temperature.

The ensemble standard deviation of the local instantaneous air–sea CO₂ fluxes that results from the perturbation of physical fields is larger than that of the global CO₂ flux, with a mean standard deviation of 0.32 mmol C m⁻² d⁻¹ for monthly means of local fluxes compared to a standard deviation of 0.0068 mmol C m⁻² d⁻¹ (0.01 Pg C yr⁻¹) for the annual global flux in FREE in the year 2020. The largest ensemble standard deviation (Fig. A1a) is generated in the Southern Ocean, the North Atlantic, and the North Pacific, which corresponds to regions of high uncertainty in existing CO₂ flux estimates (Pérez et al., 2024; Hauck et al., 2023a; Mayot et al., 2024). However, the modeled standard deviation should not be understood as the true uncertainty of the model but as a value dependent on tuning (Evensen, 2003).

2.4 Data analysis

We present CO₂ flux estimates for the period 2010–2020 that are compared to the Regional Carbon Cycle Assessment and Processes 2 (RECCAP2) global air–sea CO₂ flux estimates (DeVries et al., 2023). The RECCAP2 pCO₂ products account for oceanic outgassing of river carbon into the atmosphere. To make them comparable with our estimate stemming from a model without river carbon input, we apply a river flux adjustment (Friedlingstein et al., 2023; Regnier et al., 2022) to the RECCAP2 pCO₂ products. Thus, we quantify the anthropogenic perturbation of the ocean carbon sink (as S_OCEAN in the Global Carbon Budget; Friedlingstein et al., 2023; Hauck et al., 2020) and not the contemporary net air–sea CO₂ flux with outgassing of river carbon (as in the original RECCAP2 pCO₂ products).

To study the effect of DA on the CO₂ flux, we define regions where the time-mean air–sea CO₂ flux difference ASML − FREE ($\mathrm{\Delta }{F}_{{\mathrm{CO}}_{\mathrm{2}}}$) is pronounced based on the biome definition of Fay and McKinley (2014). Originally, these are, going polewards from the subtropics in each hemisphere, the subtropical seasonally stratified biome (STSS), the subpolar seasonally stratified biome (SPSS), and the sea-ice biome (ICE). In the Southern Ocean (denoted by subscript _SO) within the STSS_SO, we differentiate between the area where $\mathrm{\Delta }{F}_{{\mathrm{CO}}_{\mathrm{2}}}$ is positive (the assimilation leads to a flux change directed out of the ocean), referred to as region “STSS_SO+”, and the area where $\mathrm{\Delta }{F}_{{\mathrm{CO}}_{\mathrm{2}}}$ is negative, referred to as region “STSS_SO−”. All Southern Ocean regions are outlined in Fig. 5a. In the North Atlantic (denoted by subscript _NA), we consider four coherent regions within the STSS_NA and SPSS_NA outlined in Fig. 7a. The regions “Central STSS_NA−” and “Western STSS_NA+” are located in the North Atlantic STSS_NA biome and are defined by $\mathrm{\Delta }{F}_{{\mathrm{CO}}_{\mathrm{2}}}$ less than −1 mmol C m⁻² d⁻¹ and $\mathrm{\Delta }{F}_{{\mathrm{CO}}_{\mathrm{2}}}$ greater than 1 mmol C m⁻² d⁻¹, respectively. The regions “Newfoundland Basin_NA+” and “Coastal SPSS_NA−” are part of the SPSS_NA. The former is located east of Newfoundland and south of Greenland and is defined by $\mathrm{\Delta }{F}_{{\mathrm{CO}}_{\mathrm{2}}}$ greater than 3 mmol C m⁻² d⁻¹, and the latter is located off the North American coast and defined by $\mathrm{\Delta }{F}_{{\mathrm{CO}}_{\mathrm{2}}}$ less than −1 mmol C m⁻² d⁻¹. The Central STSS_NA− and Western STSS_NA+ lie on the warm side of the North Atlantic Current (NAC), and the Newfoundland Basin_NA+ and Coastal SPSS_NA− lie on the cold side of the NAC, which is evident from the modeled surface velocity field (Fig. A2a).

Within these regions, we identify the time of the year when the DA affects air–sea CO₂ flux and calculate the difference ASML − FREE for physical and biogeochemical fields. In order to assess the dynamic DA effects on surface pCO₂, it is useful to distinguish between different variables that constitute the change in pCO₂. Oceanic pCO₂ varies mainly with temperature, DIC, and alkalinity. Thus, we decompose changes in pCO₂ into their contributions from changes in SST, surface DIC, and surface alkalinity (Alk). For that, we apply the following approximations of Sarmiento and Gruber (2006) and Takahashi et al. (1993):

$$\begin{array}{}\text{(3)}& {\displaystyle }\mathrm{\Delta }p{\mathrm{CO}}_{\mathrm{2},\mathrm{DIC}}={\displaystyle \frac{p{\mathrm{CO}}_{\mathrm{2}}}{\mathrm{DIC}}}\cdot {\mathit{\gamma }}_{\mathrm{DIC}}\cdot \mathrm{\Delta }\mathrm{DIC},\text{(4)}& {\displaystyle }\mathrm{\Delta }p{\mathrm{CO}}_{\mathrm{2},\mathrm{Alk}}={\displaystyle \frac{p{\mathrm{CO}}_{\mathrm{2}}}{\mathrm{Alk}}}\cdot {\mathit{\gamma }}_{\mathrm{Alk}}\cdot \mathrm{\Delta }\mathrm{Alk},\text{(5)}& {\displaystyle }\mathrm{\Delta }p{\mathrm{CO}}_{\mathrm{2},\mathrm{SST}}=p{\mathrm{CO}}_{\mathrm{2}}\cdot \exp \left(\mathrm{0.0423}\mathit{°}{\mathrm{C}}^{-\mathrm{1}}\cdot \mathrm{\Delta }\mathrm{SST}\right).\end{array}$$

Here, differences between ASML and FREE are denoted by Δ; else, the average of ASML and FREE is used for the computation. The sensitivities γ_DIC and γ_Alk describe how pCO₂ varies with changes in one variable while keeping all other variables constant. For the sensitivities, we use an approximation derived from seawater carbonate chemistry following Sarmiento and Gruber (2006):

$$\begin{array}{}\text{(6)}& {\displaystyle }{\mathit{\gamma }}_{\mathrm{DIC}}={\displaystyle \frac{\mathrm{3}\cdot \mathrm{Alk}\cdot \mathrm{DIC}-\mathrm{2}\cdot {\mathrm{DIC}}^{\mathrm{2}}}{(\mathrm{2}\cdot \mathrm{DIC}-\mathrm{Alk})(\mathrm{Alk}-\mathrm{DIC})}},\text{(7)}& {\displaystyle }{\mathit{\gamma }}_{\mathrm{Alk}}={\displaystyle \frac{-{\mathrm{Alk}}^{\mathrm{2}}}{(\mathrm{2}\cdot \mathrm{DIC}-\mathrm{Alk})(\mathrm{Alk}-\mathrm{DIC})}}.\end{array}$$

Based on the range of valid values for γ_DIC and γ_Alk according to the explicit formulation by Egleston et al. (2010), values are excluded above 18 and below −19, respectively. This affects parts of the Southern Ocean SPSS_SO and ICE_SO biome (see white areas in Fig. 6b and c). Finally, the effect on the air–sea CO₂ flux relates directly to the pCO₂ difference at each grid point, as detailed in Orr et al. (2017, Eqs. 6–15):

$$\begin{array}{}\text{(8)}& \mathrm{\Delta }{F}_{{\mathrm{CO}}_{\mathrm{2}}}=\mathit{\alpha }\cdot k_w\cdot \mathrm{\Delta }p{\mathrm{CO}}_{\mathrm{2}},\end{array}$$

where α is the solubility of CO₂ in seawater and k_w is the gas-transfer velocity.

To evaluate the impact of DA on ocean physics, we compare the simulated SST and SSS to the assimilated observations (Sect. 2.3.1). For temperature and salinity at depth, we use the EN4-OA product (Good et al., 2013, updated to version 4.2.2). EN4-OA is an objective analysis ingesting the assimilated EN4 profile data, interpolated to global coverage on 42 depth levels. Furthermore, we compare the sea-ice concentration with remote-sensing observations from OSI-SAF 2010–2020 (EUMETSAT, 2022), the mixed-layer depth in the year 2020 with the profile-observation-based climatology of de Boyer Montégut et al. (2004, updated version 2023), and the horizontal near-surface velocities 2010–2020 with the drifter-based climatology of Laurindo et al. (2017a).

To evaluate the impact of the DA on biogeochemistry, we compare model outputs with observational datasets of surface pCO₂, DIC, alkalinity, and surface chlorophyll. To evaluate surface pCO₂, we use observations from the Surface Ocean CO₂ Atlas (SOCAT Version 2023; Bakker et al., 2023, 2016), which are provided as a monthly gridded and quality-controlled compilation. To assess DIC and alkalinity, we compare the modeled surface fields to the GLODAPv2.2023 bottle data (Lauvset et al., 2024). At depth, we compare the model output to the GLODAPv2 DIC and alkalinity climatology (Lauvset et al., 2016), which is based on observations from the period 1972–2013 and normalized to 2002. To evaluate global surface chlorophyll, we use observations from ESA-CCI, which is a multi-sensor satellite ocean-color chlorophyll a dataset with monthly global coverage (Sathyendranath et al., 2021). In addition, for the Southern Ocean, we use the mean of three satellite products (Johnson et al., 2013) that were processed with more suitable algorithms for southern high latitudes. For each observation type (OBS), we define the improvement as

$$\begin{array}{}\text{(9)}& {\mathrm{improvement}}_{\mathrm{OBS}}=|\mathrm{FREE}-\mathrm{OBS}|-|\mathrm{ASML}-\mathrm{OBS}|.\end{array}$$

3 Results

3.1 Effect of DA on ocean physics

Before we investigate the CO₂ flux, we first evaluate the effect of DA on the modeled physics. In particular, we compare the model output of both simulations with the assimilated observations to verify that the assimilation brings them into better agreement with the observations. The assimilation improves the agreement with the assimilated SST observations. On a global average, the SST in FREE is 0.14 °C colder than the observations, which is the result of an extensive cold bias in the tropics and subtropics and a warm bias in the Southern Ocean south of 40° S (Fig. 1a; mean state of SST in Fig. A3a). In addition, there are regional SST differences FREE − OBS near strong currents and in eddy-rich regions in particular, such as the NAC, the Kuroshio Current, and the Southern Subtropical Front. These SST differences are estimated to lead to a solubility-driven global air–sea flux difference of −0.06 Pg C yr⁻¹ (Eqs. 5 and 8). The assimilation increases SST in the tropics and subtropics and reduces SST south of 40° S, with particularly large effects in the Southern Ocean and in the North Atlantic (difference ASML − FREE in Fig. 1b). Thereby, the global mean model–observation difference is reduced from −0.14 to −0.12 °C and from 0.59 to 0.32 °C in absolute terms. This assimilation-induced change in SST is estimated to drive a direct solubility-driven effect on the global air–sea CO₂ flux of −0.14 Pg C yr⁻¹ (Eqs. 5 and 8). However, this global attribution is subject to high uncertainty due to the non-linear dependency of pCO₂ on temperature and because regionally large effects with opposite signs lead to uncertainty in the global mean.

Figure 1 Effect of data assimilation on sea surface temperature (SST) and sea surface salinity (SSS). All panels show the mean over the period 2010–2020. (a) The model–observation difference in SST (FREE − OSTIA). (b) The difference ASML − FREE. (c) The improvement of monthly averaged model SST relative to OSTIA, where positive denotes that the assimilation brings the model closer to observations (Eq. 9). (d–f) The same for SSS, computed with SSS from ESA-CCI.

The assimilation also improves the agreement with the assimilated SSS observations. Additional experiments with and without salinity restoring towards climatology show that the best agreement with the SSS-CCI observations is achieved by simultaneously using assimilation and restoring. A benefit of the additional use of restoring is the global coverage of the SSS climatology. FREE shows a global SSS bias (0.49 psu; Fig. 1d). The assimilation leads to a global surface freshening (Fig. 1e). There are only a few regions where SSS in FREE is fresher than the observations and where the DA consequently increases the salinity, for example, in parts of the North Atlantic. The assimilation improves the model–observation agreement in 91 % of the observed ocean area, particularly in the North Atlantic Central STSS_NA− and in the Southern Ocean STSS_SO (Fig. 1f). Tests with the assimilation of temperature alone show negative side effects of temperature assimilation on SSS in some locations (not shown). In the final setup with combined assimilation, negative effects on SSS are found in 9 % of the observed area. Globally, the mean absolute difference is reduced from 0.32 to 0.17 psu relative to the SSS observations. The direct solubility-driven effect of salinity differences on the global air–sea CO₂ flux is estimated to be negligible.

Figure 2 Effect of data assimilation on zonally averaged temperature and salinity in the upper 1000 m. All panels show the mean over the period 2010–2020. (a) The model–observation difference in temperature (FREE − EN4-OA). (b) The difference ASML − FREE. (c) The improvement of monthly averaged temperature relative to EN4-OA. (d–f) The same for salinity.

The assimilation leads to a better agreement with subsurface temperature and salinity data from the EN4-OA product in the upper 1000 m. In the upper 100–200 m of the ocean, the model–observation difference in temperature follows the surface signal (compare Figs. 1a and 2a), and the difference is reduced by the assimilation (Fig. 2b and c). At intermediate depth (roughly 200–500 m), a subsurface warm bias exists in FREE in the Southern Hemisphere at mid-latitudes (Fig. 2; mean state in Fig. A4a). This bias affects the South Pacific, South Atlantic, and southern Indian Ocean (not shown). The bias might be connected to the model's surface warm bias in the formation region of Antarctic intermediate water (Fig. 1a). Further model–observation differences exist at greater depth than 500 m, where the model's temperature is colder than the observations at almost all latitudes but warmer than the observations north of 60° N. At most latitudes and depths, the effect of the assimilation is to reduce the model–observation differences (Fig. 2c).

The model is more saline than the observations from the surface down to a depth of about 1000 m for most latitudes (Fig. 2d). This shows that the model–observation difference in this depth range follows the surface signal. The exceptions to this are at high latitudes below 200 m, where FREE is fresher than the observations. At all other latitudes, the assimilation acts towards a freshening, with the strongest effect near the surface (Fig. 2e). This improves the agreement with observations, particularly near the surface (Fig. 2d). However, the improvement is smaller at depth and even becomes negative for some latitudes at greater depth. This might be due to the limited amount of assimilated in situ salinity profiles.

The effect of the assimilation on temperature and salinity is most pronounced in the upper 1000 m and, below that, mostly decreases with depth (not shown). After the second year of assimilation, the mean absolute difference between ASML and FREE stabilizes in the range 0.35–0.36 °C for SST and 0.20–0.25 psu for SSS, while the effect of DA on subsurface temperature and salinity keeps increasing throughout the years 2010–2020.

Figure 3 Effect of data assimilation on Antarctic sea-ice concentration in September. All panels show differences in the sea-ice concentration averaged for the month September over the period 2010–2020. The 15 % line for FREE, ASML, and OSI-SAF observations is shown as a dashed, continuous, or dotted line in panel (a) or (b). (a) The difference between FREE and OSI-SAF observations. (b) The difference ASML − FREE. (c) The improvement of September mean sea-ice concentration.

Sea ice reacts dynamically to the changed ocean physical state. In the Southern Ocean, FREE is characterized by a lower sea-ice concentration compared to OSI-SAF observations. The sea-ice extent, here defined as the area where the sea-ice concentration is more than 15 %, reaches a maximum in September. The maximum extent is smaller in FREE than in OSI-SAF, which is demonstrated by the 15 % line surrounding that area for FREE and OSI-SAF (Fig. 3a; mean state in Fig. A5) and by the sea-ice concentration difference for the month September (Fig. 3b). Through DA, a higher Antarctic sea-ice concentration is obtained. This improves the agreement with OSI-SAF (Fig. 3c). During all other seasons, the assimilation leads to a higher sea-ice concentration in the Antarctic, a larger sea-ice extent, and a better agreement with OSI-SAF as well (only September is shown). In the Arctic, the differences between FREE, ASML, and OSI-SAF are regionally different (not shown).

The boundary-layer depth and mixed-layer depth are mostly reduced through DA. In particular, deep-water formation events characterized by a mixed-layer depth of more than 1000 m occur less frequently in ASML (not shown). This improves the agreement with the profile-observation-based mixed-layer climatology of de Boyer Montégut et al. (2004), reducing the mean absolute difference to the climatology from 27 to 19 m (comparison of mixer-layer depth in Fig. A6). In addition, the absolute difference of near-surface horizontal velocities to the drifter-observation-based climatology of Laurindo et al. (2017a) is reduced by about 10 % through DA (comparison of surface velocities in Fig. A7). The biological productivity near the Equator is stable in ASML and FREE, indicating that FESOM2.1-REcoM3 does not suffer from the erroneous upwelling known from previous DA studies (Park et al., 2018). The meridional overturning, however, shows spurious structures, which may point to hidden assimilation artifacts on vertical velocities. Throughout the assimilation period, spurious, spatially limited, and often deep overturning structures emerge, evolve through several months or years, and disappear in the tropical Indian, Pacific, and Atlantic basins (not shown). Thereby, the surface overturning cell sometimes breaks apart where it should extend over the Equator, exposing the bottom cell to the surface (Fig. A8b). Transport in the North Atlantic at 26.5° N, an indicator for the strength of the Atlantic Meridional Overturning Circulation, is between 8–9 Sv in FREE. In ASML, during the first 2 years of assimilation, transport at 26.5° N decreases to below 3 Sv and, during the following years, recovers to 7–8 Sv (2016–2020). One possible cause is the effect of data assimilation on the eddy parametrization (Gent and Mcwilliams, 1990). The parametrized eddy activity is relevant for the dynamics in the deep ocean, and corrupting it may have a negative impact on the large-scale oceanic circulation, as described in Sidorenko (2004, Chap. 5.5 onwards) for a previous version of the ocean model FESOM.

In summary, the ASML temperature and salinity fields from the surface to several hundred meters below and the mixed-layer depth are in good agreement with observations, and the agreement of horizontal near-surface velocities with observations is improved. This can be interpreted as an indication that the velocity field in the upper part of the ocean is also well represented. Although the spurious effects on deep-ocean circulation should be further addressed in future work, we are confident that the DA provides an improved physical state in the upper ocean, which serves as an improved basis to estimate the air–sea CO₂ flux.

3.2 Effect of DA on global CO₂ flux

The ocean absorbs 2.78 Pg C yr⁻¹ in ASML and 2.83 Pg C yr⁻¹ in FREE during 2010–2020 (Fig. 4b); thus the assimilation decreases the global mean oceanic CO₂ uptake by 0.05 Pg C yr⁻¹. The temporal evolution of the annual global CO₂ flux is similar in ASML and FREE (Fig. 4a). The first assimilation year, 2010, stands out because it is one of the very few years during which the assimilation increases the oceanic CO₂ uptake. This slightly reduces the trend in CO₂ uptake in 2010–2020 from $-\mathrm{0.40}\pm \mathrm{0.09}$ Pg C yr⁻¹ per decade in FREE to $-\mathrm{0.38}\pm \mathrm{0.11}$ Pg C yr⁻¹ per decade in ASML (negative: into the ocean). The trend thereby remains within its confidence interval. Furthermore, the assimilation slightly reduces the interannual variability in the global mean oceanic uptake, demonstrated by a standard deviation of detrended annual means of 0.11 Pg C yr⁻¹ in FREE and 0.08 Pg C yr⁻¹ in ASML (not significantly different according to an F-test). Through DA, the ensemble standard deviation of the global CO₂ flux is reduced from $\mathrm{1.0}\times {\mathrm{10}}^{-\mathrm{2}}$ Pg C yr⁻¹ in FREE to $\mathrm{0.7}\times {\mathrm{10}}^{-\mathrm{2}}$ Pg C yr⁻¹ in ASML in the year 2020.

Figure 4 Effect of data assimilation on the air–sea CO₂ flux (negative: into the ocean). (a) Annual time series of global flux in FESOM2.1-REcoM3 with ASML (black) and FREE (violet) and RECCAP2 estimates (DeVries et al., 2023) with pCO₂ products (orange) and GOBMs (blue) and their respective means (bold lines). Here, the river flux adjustment (−0.65 Pg C yr⁻¹) was applied to the pCO₂ products. (b) Time-mean global flux 2010–2018 in ASML (black) and FREE (violet) and RECCAP estimates grouped by method (DeVries et al., 2023). Crosses represent individual estimates (e.g., individual GOBMs), and bars represent the method mean (e.g., mean of 12 GOBMs). Here, the river flux term was applied to all estimates except the models following the Global Carbon Budget methodology (Friedlingstein et al., 2023). For FESOM2.1-REcoM3, the time-mean 2010–2020 is additionally shown (horizontal lines). (c) Spatial distribution of CO₂ flux averaged over the period 2010–2020 in ASML. (d) Spatial distribution of CO₂ flux difference ASML − FREE averaged over the period 2010–2020. (e) Zonal averages of CO₂ flux 2010–2020 in ASML and FREE and (f) their difference.

The strongest time-mean air–sea CO₂ flux is found at mid- and high latitudes (Fig. 4c). The large-scale pattern of the CO₂ flux is generally very similar in FREE and in ASML (FREE not shown). The largest local changes through DA, both towards stronger or weaker CO₂ fluxes, occur in the North Atlantic in the area of the NAC and in the coastal North Pacific (Fig. 4d). The most prominent large-scale effect, though, is in the Southern Ocean (Fig. 4e and f). South of 50° S, the area-integrated CO₂ uptake increases by 0.18 Pg C yr⁻¹ through the assimilation. In contrast, the uptake decreases by 0.07 Pg C yr⁻¹ between 40–50° S. With the exception of the Southern Ocean, CO₂ uptake decreases in all world oceans by a small amount (Fig. 4d).

3.3 Effect of DA on regional CO₂ fluxes and their drivers

3.3.1 Southern Ocean

In the Southern Ocean, the ocean takes up CO₂ in the annual average (Fig. 5a), with regionally heterogeneous effects of DA (Fig. 5b). While the effect of DA on surface pCO₂ and the air–sea CO₂ flux can almost entirely be explained by the combined variation in DIC and alkalinity at most latitudes north of 40° S, the thermal effect also needs to be considered in the Southern Ocean (zonal mean pCO₂ effects in Fig. A9a). In the following, we examine how the assimilation influences the air–sea CO₂ flux across individual regions in the Southern Ocean.

Figure 5 Effect of data assimilation on Southern Ocean CO₂ flux (negative: into the ocean) and its seasonality averaged over the period 2010–2020. Additionally, lines in panels (a) and (b) denote the regions, and the green hatching denotes the STSS_SO+. (a) Map of mean CO₂ flux in ASML. (b) Map of CO₂ flux difference ASML − FREE. (c–f) Seasonal cycle of air–sea CO₂ flux by region. Shading indicates the range of ensemble members in the year 2020. (g–j) Seasonal air–sea CO₂ flux difference ASML − FREE by region. Note the different scales.

STSS_SO. In the northernmost biome of the Southern Ocean, the subtropical seasonally stratified biome (STSS_SO), the mean oceanic CO₂ uptake is comparably high (Fig. 5a). The uptake is largest in austral winter and spring (June to November; Fig. 5c and d). The part of the STSS_SO characterized by a positive CO₂ flux difference ASML − FREE (positive difference: reduced uptake through assimilation), which we call the STSS_SO+, roughly forms an outer northerly ring around the STSS_SO biome (hatched area in Fig. 5a and b). The reduction in CO₂ uptake in the STSS_SO+ is greatest in winter and spring from July to October (Fig. 5g).

Figure 6 Drivers of the effects of data assimilation on pCO₂ in the Southern Ocean. Panels (a), (b), and (c) show the effects of SST, DIC, and alkalinity differences ASML − FREE simulations on surface pCO₂. Additionally, hatching inside the STSS_SO indicates where net pCO₂ is increased through the assimilation (STSS_SO+). (d, e) Seasonal cycle of SST averaged over the regions STSS_SO+ and STSS_SO− and (f, g) the difference ASML − FREE for each region. (h) Potential density profiles for the SPSS_SO, with FREE (violet line) and ASML (dashed black line) based on daily T and S and with EN4-OA (dotted green line) based on monthly T and S. (i) DIC profiles for the SPSS_SO, showing FREE (violet line), ASML (dashed black line) from 2010–2020, and climatological DIC from GLODAP. Shading in panels (d), (e), (h), and (i) indicates the range of ensemble members in the year 2020.

The increase in pCO₂ in the STSS_SO+ is partly driven by lowered alkalinity and partly by increased surface DIC (Fig. 6b and c). These, along with the colder SST and fresher SSS in the STSS_SO+ (Fig. 1b and e), are indications for a year-round stronger influence of subantarctic waters. This is evident from typical water properties in the subantarctic and subtropical Southern Ocean. In the subantarctic, surface DIC is higher, surface alkalinity is lower, temperature is colder, and salinity is lower (maps of SST, SSS, DIC, and alkalinity in Fig. A10). In the fragmented area of the STSS_SO+, different factors contribute to regional changes in the surface DIC and alkalinity budget in ASML (sources minus sinks of DIC and alkalinity in Fig. A11). Depending on location, an increased upward transport of DIC through mixing, an increase in DIC through a reduced biological sink of DIC in spring, or a decrease in alkalinity through changes in horizontal and vertical advection dominates. The seasonality of the effect of DA on the air–sea CO₂ flux in the STSS_SO+ (Fig. 5c and g) is determined by seasonal temperature differences between ASML and FREE (Fig. 6d and f). During summer, SST is slightly reduced (Fig. 6f), which lowers pCO₂ (Fig. 6a). This counteracts the effects of DIC and alkalinity on pCO₂ (Fig. 6b and c) and thus dampens the overall DA effect on the air–sea CO₂ flux during summer.

The part of the STSS_SO characterized by a negative CO₂ flux difference ASML − FREE, which we call the STSS_SO−, is a fragmented region and roughly consists of segments of an inner southerly ring (non-hatched area in Fig. 5a and b). Here, the increase in CO₂ uptake through DA is largest in summer and fall (November to April; Fig. 5h). The reduction in pCO₂ is driven by increased alkalinity and partly also by lower surface DIC (Fig. 6b and c, non-hatched area). These, together with higher SST in ASML than FREE in the STSS_SO− regions (Fig. 1b), indicate a higher presence of subtropical waters (see characteristics of subtropical waters in Fig. A10). Where there is lower DIC in the STSS_SO− in ASML (Fig. 6b), this can mostly be explained by an increased biological sink of DIC, with the addition of sharply defined local changes in horizontal advection of DIC and alkalinity (Fig. A11). Additionally, seasonal temperature effects occur. During winter, SST is higher in ASML than in FREE (Fig. 6e and g). This increases pCO₂ in the STSS_SO− (Fig. 6a), counteracting the effects of lower DIC and higher alkalinity on pCO₂ and dampening the overall DA effect during winter.

The contrasting effects in the STSS_SO indicate a horizontal shift in water masses within the STSS_SO biome. In the center of the STSS_SO, the Subantarctic Front is located, which is associated with the Antarctic Circumpolar Current (ACC) and characterized by a strong gradient in SST, SSS, and various other tracers (Chapman et al., 2020). Because SST and SSS are directly influenced and improved by the assimilation, the position of this front is also expected to change as a result of the assimilation, leading to a horizontal relocation of waters separated by the front. With the relocation of the front, dynamic shifts in regional characteristics occur, such as the amount of DIC and alkalinity transported vertically through mixing and the biological sources and sinks of DIC and alkalinity.

SPSS_SO. Further south, in the subpolar seasonally stratified biome (SPSS_SO), the ocean absorbs CO₂ all year round (Fig. 5a). The oceanic uptake is increased through the assimilation, shown by a negative flux difference ASML − FREE in Fig. 5b. The largest difference between ASML and FREE is seen in spring from September to October (Fig. 5i). Due to the seasonally varying effect of DA, the seasonal cycle of the CO₂ flux in the SPSS_SO is altered. In ASML, the CO₂ uptake is weakest in February, gets stronger in fall (MAM), stagnates in winter (JJA), and resumes growth in spring (SON), reaching peak uptake in November (Fig. 5e). In FREE, the CO₂ uptake weakens in winter, is weakest in September, and gets stronger afterwards, reaching peak uptake in December.

In the SPSS_SO, the increased CO₂ uptake and lower surface pCO₂ during winter and spring are driven by a combination of colder temperatures and lower DIC (Fig. 6a and b), which outweighs the opposite effect of a decrease in alkalinity on pCO₂ (Fig. 6c; relative importance of thermal effect in Fig. A12a). Surface DIC is generally high due to upward transport of carbon-rich deep water (e.g., Hauck et al., 2023a). The reason for lower surface DIC in ASML is that the upward transport through mixing is reduced (Fig. A11) through a more stable stratification, which is also evident from a reduced density in the upper 300 m and an increased density below that (Fig. 6h). Thereby, the densities in the SPSS_SO agree better with densities calculated from EN4-OA. The boundary layer and mixed layer in winter and spring are shallower and thereby in better agreement with the observation-based climatology (Fig. A6). Vertical mixing within the boundary layer affects the vertical profiles of DIC and alkalinity, towards lower DIC in ASML above 100 m and higher DIC below (Fig. 6i). The vertical profile of DIC in ASML is closer to GLODAP DIC observations, although some differences to GLODAP still exist. Beside the fact that the differences in stratification and boundary-layer depth affect the vertical DIC profile, they also imply fewer available surface nutrients in ASML. Probably due to a combination of lower nutrient availability and colder surface temperature, ASML features lower net primary production (NPP), lower chlorophyll concentrations, and a lower phytoplankton biomass in the SPSS_SO (not shown). Thereby, the modeled biogeochemical cycle adjusts to the lower transport of nutrients to the surface by transferring less organic material to depth, ultimately acting to compensate about 60 % of the difference in physical transport of DIC (Fig. A13a) and adding to the reduction in surface alkalinity (Fig. A13b). Within the SPSS_SO (roughly south of 50° S), differences between FREE and ASML in terms of the temperature effect on pCO₂, vertical transport of DIC and alkalinity, and biological sources and sinks are larger than at any other latitude (Fig. A13).

ICE_SO. In the seasonally ice-covered biome (ICE_SO) surrounding the Antarctic continent, the time-mean CO₂ flux is smaller than in other biomes (Fig. 5a). In this region, the ocean absorbs CO₂ during summer and there is a smaller outgassing during winter (Fig. 5f), as the region is mostly ice-covered in winter (see sea-ice concentration in September in Fig. 3). In the northern part of the ICE_SO biome, close to the SPSS_SO, the effect of the assimilation is similar to the effect within the SPSS_SO itself (Fig. 5b). Here, the assimilation acts to increase ocean CO₂ uptake or to weaken CO₂ outgassing during winter and spring (Fig. 5i and j). Thereby, interestingly, the assimilation hinders outgassing of CO₂ from May to November in ASML in the ICE_SO biome (Fig. 5f; comparison of winter outgassing with other estimates in Fig. A14). The reduced outgassing and decreased pCO₂ during winter and spring are driven by similar processes to within the SPSS_SO. Again, lower surface DIC and colder temperatures (Fig. 6a and b) outweigh the opposite effect of a decrease in alkalinity on pCO₂ (Fig. 6c). As in the SPSS_SO, the reason for the decrease in pCO₂ is reduced surface DIC and increased DIC below 100 m as a result of less upward transport of DIC through mixing (Fig. A11) in a more stable stratification due to surface freshening (Fig. 1e). In addition, as the surface temperature is lower in ASML (Fig. 1b), the winter sea-ice concentration is higher (Fig. 3b), which prevents winter outgassing of CO₂. In the southern part of the ICE_SO biome, near the Antarctic continent, the effect of the DA on the CO₂ flux is small.

In summary, in the Southern Ocean, the main effects of the DA on the CO₂ flux are, firstly, an increase in the uptake in the SPSS_SO caused by surface cooling and by a more stable stratification and thus less upward transport of naturally carbon-rich water through mixing and, secondly, an overall lower CO₂ uptake in the STSS_SO as a consequence of a spatial redistribution of fluxes near the Subantarctic Front.

3.3.2 North Atlantic

In the North Atlantic, the assimilation has noticeable effects on the CO₂ flux in the area of the North Atlantic Current, where the ocean absorbs CO₂ in the annual average (Fig. 7a). During summer, however, the ocean releases CO₂ while the sea surface warms (Fig. 7c–f). In the Central STSS_NA−, the effect of the DA is to prevent outgassing during summer (Fig. 7c and g). In the Western STSS_NA+ and in the Newfoundland Basin_NA+, the ocean CO₂ uptake is decreased during winter (Fig. 7d, e, h, and j). The regionally different dynamics of the effects of the assimilation that drive these differences in the air–sea CO₂ flux in the North Atlantic are investigated next.

Figure 7 Effect of data assimilation on North Atlantic CO₂ flux (negative: into the ocean) and its seasonality averaged over the period 2010–2020. (a) Map of mean CO₂ flux in ASML. (b) Map of CO₂ flux difference ASML − FREE. (c–f) Seasonal cycle of air–sea CO₂ flux by region. Shading indicates the range of ensemble members in the year 2020. (g–j) Seasonal air–sea CO₂ flux difference ASML − FREE by region. Note the different scales.

Figure 8 Drivers of the effects of data assimilation on pCO₂ in the North Atlantic. Panels (a), (b), and (c) show the effects of SST, surface DIC, and alkalinity differences ASML − FREE on surface pCO₂. (d) Difference in boundary-layer depth (ASML − FREE) for spring (MAM) 2010–2020, where positive denotes a shallower boundary layer in ASML. (e) Difference in surface chlorophyll (ASML − FREE) for spring (MAM) 2010–2020. (f) Potential density profiles for the Newfoundland Basin_NA+ region, with FREE (violet line) and ASML (dashed black line) based on daily T and S and with EN4-OA (dotted green line) based on monthly T and S. (g) DIC profiles for the Newfoundland Basin_NA+ region, showing FREE (violet line), ASML (dashed black line) from 2010–2020, and climatological DIC from GLODAP. The shading in panels (f) and (g) indicates the range of ensemble members in the year 2020.

Central STSS_NA−. In the Central STSS_NA−, the effect of the DA is overall towards a more negative flux in CO₂ from May to November (Fig. 7g). Thus, spring and fall CO₂ uptake are increased and summer outgassing is prevented in ASML (Fig. 7c). The reason for decreased surface pCO₂ is higher alkalinity in ASML (Fig. 8c). In this region, the alkalinity effect, which reduces pCO₂, outweighs the opposing effects of DIC and SST on pCO₂ (Fig. 8a and b). A higher alkalinity could point to the presence of waters of subtropical origin transported northward with the NAC (Völker et al., 2002). Other fingerprints of waters transported by the NAC are a warm SST particularly in winter, a higher salinity, and higher DIC than that of North Atlantic subpolar waters (maps of SST, SSS, DIC, and alkalinity in Fig. A15; Völker et al., 2002). The assimilation causes a change in these properties towards a higher SST, higher salinity, and higher DIC in the Central STSS_NA−. Simultaneously, ASML represents a deeper boundary layer in this region (Fig. 8d). While changes in the North Atlantic mixed-layer depth overall result in a spatial pattern in ASML that more closely aligns with the pattern in the observation-based mixed-layer climatology, the modeled mixed layer in the simulations is still overall deeper than in the climatology, leading to less agreement in the Central STSS_NA− (Fig. A6). Likely facilitated by higher SST and more available nutrients through deeper mixing in winter and spring, ASML features a higher biological sink of DIC above 190 m (Fig. A16d); more biological carbon export through sinking of detritus at 190 m; and more column-integrated phytoplankton biomass and surface chlorophyll in spring, which is illustrated by the example of surface chlorophyll difference between ASML and FREE in Fig. 8e. In combination, the higher alkalinity associated with NAC transport and the higher biological sink of DIC result in lowered surface pCO₂ and higher oceanic uptake.

Western STSS_NA+. In the Western STSS_NA+, the DA reduces the CO₂ uptake and increases pCO₂ mainly during winter as a direct effect of increased SST (Fig. 8a). The direct thermal effect is dominant over the combined effect of DIC and alkalinity (relative importance of thermal effect in Fig. A12b). The latter have effects comparable in magnitude to SST but mostly cancel each other out (Fig. 8b–c). The effect of DA on surface properties (SST, SSS, DIC, and alkalinity) in the Western STSS_NA+ is similar to the effect in the Central STSS_NA−, which indicates a higher influence of subtropical waters in both regions.

Newfoundland Basin_NA+. In the Newfoundland Basin_NA+, the dominant effect of DA is a reduction in the CO₂ uptake and an increase in pCO₂ mainly during winter as a direct effect of increased SST (Fig. 8a). In addition, ASML also features a more stable stratification due to lower density at the surface than FREE (Fig. 8f), which mostly affects DIC at 50-400 m depth through reduced subduction of DIC (Fig. 8g). Furthermore, ASML represents less surface chlorophyll in the Newfoundland Basin_NA+ (Fig. 8e) as a result of a redistribution of biomass from the surface to 50–400 m depth due to spring mixing (not shown). The downward mixing of biomass results in an increase in the biological sink of DIC above 50 m, likely due to more primary production near the surface, and a decrease in the biological sink at 50–400 m, likely due to more remineralization at this depth. However, the differences in the biological sink of DIC are compensated by the mixing of DIC (profiles not shown). Overall, differences in the regional DIC profile to the observational GLODAP climatology slightly increase (Fig. 8g).

Coastal SPSS_NA−. In the Coastal SPSS_NA−, pCO₂ is reduced and the ocean CO₂ uptake is increased in ASML during winter and spring (Fig. 7f and j). The reduction in pCO₂ is facilitated by colder SST (Fig. 8a). This might be due to subpolar water masses penetrating further south along the coast in ASML because the location where the current separates from the coast is further south in ASML (velocities in Fig. A2).

In summary, DA affects the CO₂ flux in the North Atlantic mainly through changes in SST, combined with changes in horizontal advection of DIC and alkalinity near the NAC. Changes in the vertical mixing of DIC and alkalinity are largely compensated by feedbacks in biogeochemical cycles. Which of these effects is dominant, however, varies from region to region.

3.4 Comparison with biogeochemical observations

3.4.1 pCO₂ (SOCAT)

To evaluate the modeled air–sea CO₂ flux based on observations, surface pCO₂ is the most informative variable, as it is closely related to the air–sea CO₂ flux. Effects of the DA on the modeled ecosystem and associated carbon fluxes, along with thermal and dynamical effects that affect the CO₂ flux, are all included in pCO₂. The global mean of absolute monthly model–observation differences to the available SOCAT pCO₂ observations is 27.26 µatm for FREE. For ASML, the difference is slightly larger at 27.60 µatm. On global average, pCO₂ is higher than in SOCAT by 3.70 µatm in FREE and 4.59 µatm in ASML, as regions with positive and negative differences to SOCAT compensate (Fig. 9a). As an illustration of the regional changes through DA, the absolute differences in pCO₂ amount to 8.08 µatm (absolute difference ASML − FREE calculated at every grid point then averaged globally), which is ±27 % of the mean absolute model–observation difference. A linear offline estimation demonstrates that this change in pCO₂ would lead to an absolute change in the air–sea CO₂ flux by 1.06 mmol C m⁻² d⁻¹ on average (Eq. 8).

Figure 9 Partial pressure of CO₂ (pCO₂) at the surface averaged over the years 2010–2020. Panels (a)–(c) show the difference between FREE and SOCAT observations in (a) the global ocean, (b) the Southern Ocean, and (c) the North Atlantic; panels (d)–(f) show the impact of the assimilation as “improvement” relative to SOCAT observations computed from monthly mean pCO₂ in the same regions. Positive values (green) denote a reduced difference to SOCAT.

Overall, FREE and ASML show very similar regional pCO₂ differences compared to SOCAT (difference of FREE and SOCAT in Fig. 9a; difference of ASML and SOCAT not shown). In the subtropical and tropical Atlantic and the subtropical Pacific, FREE and ASML have higher pCO₂ than SOCAT, while, in the equatorial Pacific, pCO₂ is lower. At high latitudes, FREE and ASML represent mostly lower pCO₂ than SOCAT.

In the Southern Ocean, the simulations represent lower pCO₂ than SOCAT in the SPSS_SO and ICE_SO biomes in the annual mean (Fig. 9c), which is dominated by summer differences to SOCAT (not shown) when most observations are available. Through the assimilation, pCO₂ is slightly increased in summer and mostly reduced in winter (not shown), leading to an overall better agreement with SOCAT (Fig. 9e). In contrast, in the STSS, FREE and ASML represent higher pCO₂ than SOCAT, and, through the assimilation, the agreement with SOCAT decreases.

In the North Atlantic, the simulations and SOCAT show a similar large-scale pattern, namely that pCO₂ is higher in the subtropics (ASML around 400 µatm) than in the subpolar regions (ASML around 280 µatm). However, this latitudinal difference of pCO₂ is stronger in the simulations compared to SOCAT, meaning that, in the subtropics, pCO₂ in the simulations is higher than in SOCAT (Fig. 9d), while it is lower in the subpolar regions. Furthermore, in both simulations, there is a pronounced pCO₂ surface gradient in the NAC and North Atlantic Subpolar Gyre region, whose position is changed by the assimilation and which appears to be further northward in SOCAT. Thereby, the assimilation overall leads to a better agreement with SOCAT, in particular through a decrease in pCO₂ in the Central STSS_NA−, where the average difference is reduced from 26 µatm (FREE − SOCAT) to 1 µatm (ASML − SOCAT). However, in the Newfoundland Basin_NA+, the average difference is reversed from −17 µatm (FREE − SOCAT) to 13 µatm (ASML − SOCAT), which is associated with a larger absolute discrepancy of ASML and SOCAT.

3.4.2 DIC and alkalinity (GLODAP)

DIC and alkalinity are two of the most important variables from which pCO₂ is derived (Sect. 3.3). Comparing them with observations provides more insights into the strengths and weaknesses of the modeled carbonate system than a comparison with pCO₂ observations alone. The FESOM2.1-REcoM3 simulations represent higher surface DIC than GLODAP bottle observations (Lauvset et al., 2024, gridded monthly means) on average (Fig. 10a), with a global mean surface difference FREE − GLODAP of 6.46 mmol C m⁻³ for DIC. Although fewer DIC observations are available than pCO₂ observations, similarities between the respective model–observation differences for DIC and pCO₂ can be recognized. For example, DIC in the model is lower in the tropical and subtropical Atlantic than GLODAP and is higher in the polar Atlantic. This is consistent with SOCAT pCO₂ observations in the same areas. The model–observation differences to GLODAP DIC and SOCAT pCO₂ are also consistent with each other in the North Pacific. The assimilation induces absolute changes in surface DIC of 6.33 mmol C m⁻³ on global average, with regional differences in sign. These changes slightly reduce the mean absolute difference to the surface observations from 32.78 to 32.15 mmol C m⁻³ and yield a mixed picture of the improvement (Fig. 10b).

Figure 10 Comparison of the model results with surface DIC and alkalinity bottle observations from GLODAP globally over the years 2010–2020. (a) Surface DIC differences FREE − GLODAP. (b) Improvement of monthly surface DIC relative to GLODAP. (c, d) For alkalinity.

While the trend in surface DIC due to anthropogenic input makes it necessary to compare the model with contemporaneous observations at the ocean surface, a comparison with climatological data is meaningful below a depth of approximately 200 m. In fact, the modeled global distribution of DIC at depth is overall similar to that in the GLODAP climatology for both simulations (zonal mean DIC surface to 1000 m depth in Fig. A17). For example, the model results and GLODAP datasets show that DIC is lowest in the isopycnals of the subtropical gyres (2050–2150 mmol C m⁻³; Fig. A17a) and that DIC mostly increases with depth and is higher in the Pacific (2420 mmol C m⁻³ at 1000 m in the North Pacific) than in the Atlantic (2320 mmol C m⁻³ below 3000 m in the South Atlantic). However, depending on the ocean basin and depth, there can be both negative and positive differences between the simulations and the GLODAP climatology, which are in the order of 20 mmol C m⁻³ (Fig. A17c). On a global average, the assimilation leads to an increase in DIC between 200–600 m depth and a reduction in DIC between the surface and 200 m, with the largest effect in the upper 400 m (Fig. A17b). This leads to an improved agreement with the GLODAP climatology, with the largest global mean improvement at a depth of 400 m (2.5 mmol C m⁻³; Fig. A17d). Below 1000 m depth, the global mean absolute difference FREE − ASML is only 1–2 mmol m⁻³ of DIC and alkalinity and is therefore substantially smaller than at the surface.

The comparison with GLODAP bottle alkalinity at the surface shows a similar spatial pattern to that of DIC (see Fig. 10a and c). The magnitude of the bias is also comparable (14 mmol Alk m⁻³). The global mean of the absolute difference ASML − FREE of surface alkalinity is 7.72 mmol Alk m⁻³. The assimilation leads to a reduction in the absolute difference of the model alkalinity to GLODAP from 34.34 to 32.60 mmol Alk m⁻³. Since the effects of physics assimilation on alkalinity and DIC are regionally consistent, regions of improved or deteriorated agreement with GLODAP often coincide for both variables (compare Fig. 10b and d). Because changes in DIC and alkalinity have an opposing effect on the CO₂ flux, it is likely that their correlation results in compensating effects. A linear estimate shows that the joint effect of DIC and alkalinity changes is responsible for a change in the CO₂ flux in the order of 1.22 mmol C m⁻² d⁻¹ on average, and, globally integrated, the assimilation-induced changes in DIC and alkalinity lead to an estimated net increase in the air–sea CO₂ flux in the order of 0.50 Pg C yr⁻¹ (Eqs. 4, 3, and 8). However, this linear offline estimate is subject to a large uncertainty because regionally large effects with opposite sign lead to uncertainty in the global mean.

3.4.3 Surface chlorophyll (OC-CCI)

The representation of chlorophyll by the model is of interest as a proxy for primary production. Surface chlorophyll reflects the phytoplankton state and biomass; therefore, effects of the DA on the biological model state can be seen in the total surface chlorophyll concentration. A comparison of the modeled surface chlorophyll with remotely sensed chlorophyll from OC-CCI reveals that both simulations feature a higher surface chlorophyll concentration than OC-CCI (FREE − OBS in Fig. 11a and c). In FREE, the difference to OC-CCI is 0.02 mg m⁻³ on global average, with low deviations in the tropics and an enhanced difference north of 30° N (0.12 mg m⁻³) and south of 30° S (0.24 mg m⁻³). Apart from this, both simulations capture the global distribution of chlorophyll well. The simulations show the seasonal maxima in each hemisphere around 1 month earlier in the year (not shown). South of 30° S, FREE is in better agreement with chlorophyll a from the Southern Ocean-specific chlorophyll product of Johnson et al. (2013) (Fig. 11b) than with OC-CCI data (Fig. 11a).

Figure 11 Surface chlorophyll for the years 2010–2020. (a–c) Difference between FREE and SOCAT observations in (a) the global ocean, (b) the Southern Ocean, and (c) the North Atlantic; (d–f) impact of the assimilation as “improvement” relative to the observations in the same regions. Panels (a), (c), (d), and (f) compare to monthly OC-CCI observations, and panels (b) and (e) refer to the climatology for 1998–2019 by Johnson et al. (2013).

On global average, the assimilation slightly reduces the differences between model and OC-CCI data, from a global mean absolute difference of 0.31 to 0.29 mg m⁻³. The assimilation changes the chlorophyll concentration by an absolute value of 0.05 mg m⁻³ on average, which is 15 % of the global mean absolute difference to OC-CCI. There are regions in which assimilation leads to a reduction in chlorophyll and thus to better agreement with the satellite products, for example, in the North Atlantic Subpolar Gyre and the Southern Ocean SPSS_SO (Fig. 11e and f). In contrast, the model reacts to the DA with an increase in chlorophyll in the North Atlantic Central STSS_NA− and the Argentine Basin, which leads to poorer agreement.

4 Discussion

The improvement in temperature and salinity overall leads to a heterogeneous picture in biogeochemistry. While near-surface temperature and salinity fields are improved through DA almost everywhere, the global mean absolute difference in modeled surface pCO₂ to SOCAT remains similar in ASML compared to FREE, and this also applies to the model–observation differences for surface chlorophyll, DIC, and alkalinity (Sect. 3.4). Where improvements in one BGC variable occur, these do not necessarily lead to consistent improvement in all BGC variables. For example, the representation of pCO₂ improves, while that of chlorophyll deteriorates in the North Atlantic Central STSS_NA− (Figs. 11f and 9f). In the Southern Ocean SPSS_SO, the reduction in modeled surface chlorophyll in spring and the increase in pCO₂ in summer lead to a better agreement with pCO₂ observations, yet the available observations of DIC and alkalinity do not resolve the regional scales to evaluate the corresponding changes in these variables (Figs. 9, 10, and 11f). The uncertainty represented by the ensemble is reduced by the DA, which has the most obvious effect on the directly assimilated fields (SST in Fig. 6d and e and density in Fig. 8f). The ensemble standard deviation of the CO₂ flux, where it is large in FREE, is constrained by the DA to globally more uniform and smaller values (Figs. 5c–f, 7c–f, and A1). Only in the North Pacific is the standard deviation of CO₂ fluxes equally high in ASML and FREE, precisely in a region that also presents a challenge for pCO₂ products (compare Fig. A1 and Mayot et al., 2024, Fig. 5a). In the rest of the ocean, the reduced uncertainty represented by the ensemble does not necessarily coincide with improved agreement with BGC observations. One possible reason for improvement of model–data mismatch in one variable with worsening in another may lie in inconsistencies between the observational datasets. Another reason may be missing processes in the model and the use of constant BGC model parameters. Those parameters are responsible for linking changes between ecosystem variables, and, in reality, they vary across space and time depending on species composition in the ecosystem (Mamnun et al., 2023, Chap. 3). Overly simplified links between ecosystem variables can lead to canceling errors, which means that the state of one variable may worsen as a result of improving the other through DA (as in Ford and Barciela, 2017). For example, surface chlorophyll (Fig. 11f) and pCO₂ (Fig. 9f) in the central Greenland Sea deteriorate in response to improvements in SST (Fig. 1c), SSS (Fig. 1f), and sea-ice concentration (not shown). This could indicate that the BGC parametrization compensates for flaws in the free-running physical model in this region. The parameter mismatch might cause difficulties in modeling the change in BGC variables under the ongoing loss of Arctic sea ice (Chen et al., 2016).

The major effects of physics DA on BGC variables seem to follow changes in SST and are largely uniform over the full period of DA (Sect. 3.4). Surface chlorophyll changes show a pattern similar to SST changes (Figs. 1 and 11). The modeled phytoplankton growth is temperature-dependent (Gürses et al., 2023). Furthermore, indirect temperature effects on plankton dynamics due to stratification and mixing changes contribute, although those can have heterogeneous effects and the correlation of chlorophyll and boundary-layer depth is less clear (not shown). The changes in surface DIC and alkalinity show similar spatial patterns with regional heterogeneity (Sect. 3.3), again with the major changes being coherent with the changes in SST (Fig. 1). Furthermore, the effects of the assimilation on DIC and on temperature in the upper 1000 m correlate regionally: cooling through DA at intermediate depth (Fig. 2b) is usually accompanied by higher DIC in ASML (Fig. A17b), while warming through DA near the surface occurs together with reduced DIC in ASML. An overall more stable ocean stratification in the upper hundreds of meters explains why. On global average, the assimilation leads to lower DIC above 200 m and higher DIC between 200–600 m depth. In regions of substantial DA effects on vertical transport of DIC, for example, in the Central STSS_NA− or in the SPSS_SO (Sect. 3.3), the modeled biogeochemical cycles adjust dynamically to the altered vertical transport. The resulting changes in biological sources and sinks of DIC compensate for 20 %–70 % of the changes in vertical transport of DIC (Fig. A13a). In addition to changes in stability and mixing, the assimilation affects the distribution of DIC and alkalinity through local changes in near-surface horizontal transport. As the horizontal distribution of surface DIC, alkalinity, and SST is governed by latitudinal gradients and common pathways of transport (Figs. A10 and A15), all of them undergo similar changes as the SST field is modified. An exception to this is in the STSS_SO, where regional shifts along contrasting surface gradients of DIC, alkalinity, and temperature affect the respective variables differently (Sect. 3.3). These shifts change the spatial pattern of air–sea CO₂ fluxes. With the exception of the Southern Ocean, zonally averaged changes in surface pCO₂ are dominated by the combined effects of surface alkalinity and DIC on pCO₂ (Fig. A9a). Because alkalinity and DIC are usually modified according to the same pattern through mechanisms acting on both, their effects on pCO₂ are anticorrelated (Fig. A9b). The direct thermal effect on pCO₂ can still be the largest locally, for example, in the North Atlantic Newfoundland Basin_NA+ (Fig. A12b). While the DA dynamically induces changes in surface pCO₂ everywhere, the strongest effects on the air–sea CO₂ flux are at high latitudes, where pCO₂ changes are amplified by high wind velocities.

The net effect of DA on the global air–sea CO₂ flux varies from year to year between −0.12 and 0.15 Pg C yr⁻¹, which is small compared to the changes in regional CO₂ fluxes. The global net effect of lateral redistribution of alkalinity and DIC at the ocean surface is a result of compensation between regions where alkalinity and DIC are added and removed. Similarly, regional SST effects on surface pCO₂ mostly balance out globally because DA primarily induces a correction of regional SST biases, reducing the mean absolute difference to the observations from 0.59 to 0.32 °C, rather than changing the global mean SST, which differs by only 0.02 °C between FREE and ASML. DA-induced differences in vertical transport of DIC are comparably large south of 50° S, but approximately 95 % of them are balanced globally by opposing changes in vertical transport further north (vertical transport of DIC in Fig. A13a). In particular, the effect of DA on subduction of DIC through vertical advection into the ocean's deeper layers (not shown), which is the rate-limiting step on oceanic uptake of anthropogenic CO₂ emissions (DeVries, 2022), appears small, which may be due to an insufficient number of deep observations. Besides, experiments on longer timescales might be necessary to generate a visible effect of deep-circulation changes on the ocean's carbon cycle (Cao et al., 2009), which could, however, lead to imbalances in the CO₂ flux (Lebehot et al., 2019; Kriest et al., 2020; Primeau and Deleersnijder, 2009). Another possible reason why the DA effect on the global CO₂ flux in our simulation is small is the variable stoichiometry in REcoM. The dynamic biological functioning reduces the sensitivity of critical fields, like DIC, to physical changes (Buchanan et al., 2018). Furthermore, negative feedback effects between surface alkalinity, DIC, atmospheric pCO₂, and air–sea fluxes might reduce the overall response (Bunsen et al., 2024).

The overall impact of the DA on the air–sea CO₂ flux on a global scale is modest (0.05 Pg C yr⁻¹) compared to the differences between other estimates (e.g., a standard deviation of 0.45 Pg C yr⁻¹ of GOBMs in DeVries et al., 2023). The global air–sea CO₂ flux estimates of FREE and ASML fall in the range of previous model estimates and in the range of previous pCO₂ products (Fig. 4a and b) for the period 2010–2018, during which comparable estimates are available (DeVries et al., 2023). We compare here to two other data-assimilating BGC model approaches, namely ECCO-Darwin (global; Carroll et al., 2020) and B-SOSE, which is restricted to the Southern Ocean (Verdy and Mazloff, 2017). Both approaches use linearized least-squares optimization data assimilation methods (4D-var/adjoint and Green's function; Wunsch, 1996; Menemenlis et al., 2005). However, the largest difference to our study is probably that they assimilate BGC observations in addition to physical data. Thus, as expected, the effect on pCO₂ in our study is smaller (3 %) than in ECCO-Darwin and B-SOSE, where a reduction in pCO₂ model–data misfit of 6 % and 64 % was reported, respectively (here given as quadratic misfit). The global CO₂ flux (2010–2018) is smaller in FESOM2.1-REcoM3 (−2.73 Pg C yr⁻¹ in FREE and −2.78 Pg C yr⁻¹ in ASML) than in ECCO-Darwin (−3.13 Pg C yr⁻¹). The discrepancy between the CO₂ flux estimates based on models and pCO₂ products is an area of active research and not fully resolved (Friedlingstein et al., 2023; DeVries et al., 2023). On the one hand, model biases in the Atlantic Meridional Overturning Circulation, in Southern Ocean ventilation, and possibly in the surface ocean carbonate chemistry were suggested as reasons why models might underestimate the global mean CO₂ uptake in recent decades (Friedlingstein et al., 2023; Terhaar et al., 2024, 2022). On the other hand, the sparsity of observations is a concern for the pCO₂ products. According to one testbed simulation, the pCO₂ products reflect the global mean and the seasonal cycle relatively well, while the decadal variability may be overestimated (Gloege et al., 2021). An overestimation of the decadal trend, as suggested by Hauck et al. (2023b), could explain the high estimates of the pCO₂ products for the present-day global mean CO₂ flux. In contrast, for the North Atlantic, it was argued that pCO₂ is comparatively well constrained by observations in the last decade but not in the 1980s, which has an erroneous influence on the long-term trend (Pérez et al., 2024).

The effects of data assimilation on the CO₂ flux are most pronounced in the Southern Ocean STSS_SO and SPSS_SO in winter. Verdy and Mazloff (2017) also found the largest effects of assimilation on the CO₂ flux in this region. Although the region is of crucial importance for the global ocean carbon sink, it also has the greatest uncertainty due to the lack of ship-based winter observations (Friedlingstein et al., 2023; Hauck et al., 2020). In the last decade, the number of winter observations has increased due to the introduction of biogeochemical Argo floats (Johnson et al., 2017; Williams et al., 2017), although the float-based pCO₂ derived from pH measurements and estimated alkalinity is subject to higher uncertainty compared to direct pCO₂ measurements (Williams et al., 2017; Bakker et al., 2016). Machine learning approaches incorporating BGC Argo float observations suggest a stronger winter outgassing around Antarctica, particularly south of 50° S in the SPSS_SO and ICE_SO biomes, for 2015–2017 (Bushinsky et al., 2019; Gray et al., 2018). This results in a lower estimate of annual Southern Ocean CO₂ uptake in the float products. One suggestion in the literature is that model inadequacies in the representation of mixing and upwelling in the Southern Ocean might cause the discrepancy between float products and models (Gray et al., 2018). However, improvements in the modeled ocean physics and changes in mixing through data assimilation do not lead to closer agreement between the FESOM2.1-REcoM3 estimate and the float products (comparison of FESOM2.1-REcoM3, float products, and B-SOSE in Fig. A14). In contrast, ASML shows even weaker winter outgassing and stronger summer uptake south of 50° S than FREE, which brings the FESOM2.1-REcoM3 estimate further away from the float products. However, ASML is brought close to B-SOSE in terms of winter outgassing in the Antarctic polar ocean south of 60° S and winter uptake in the STSS_SO around 40° S. Additionally, airborne CO₂ flux estimates and direct pCO₂ measurements stemming from a sail drone have questioned the estimates of winter outgassing based on the BGC floats, attributing the high pCO₂ values either to possible biases in the floats' measuring devices or to anomalously high pCO₂ in the years 2015–2016 (Long et al., 2021; Sutton et al., 2021).

5 Conclusion

We apply data assimilation of temperature and salinity into a global ocean–biogeochemical model to improve the physical state for the years 2010–2020. The simulation is then assessed with regard to the effects on the biogeochemical variables. The experiments show that the effect of data assimilation (DA) on biogeochemical variables is mostly related to temperature changes. While the air–sea CO₂ flux and pCO₂ are directly affected by sea surface temperature, the DA also induces indirect changes to pCO₂ through dissolved inorganic carbon (DIC) and alkalinity. Globally integrated, these are more relevant for pCO₂ than the direct temperature effect. However, which of these factors has a dominant effect on pCO₂ varies locally. The assimilation leads to regional shifts in areas of CO₂ outgassing and uptake. Local effects on the air–sea CO₂ flux are particularly large in dynamic regions such as the North Atlantic Current and near the Subantarctic Front, whose pathways are challenging for the model to resolve without DA. The largest effect on the air–sea CO₂ flux occurs in the Southern Ocean during winter. In the simulation with assimilation, the uptake south of 50° S is increased due to shallower mixing and surface cooling, and the uptake northward of that (40–50° S) is weakened. In this area of the ocean, the uncertainty in current estimates of CO₂ fluxes is particularly high. Overall, the uncertainty inherent to the biogeochemical model appears to be larger than the uncertainties induced through physical biases in the free-running model. Locally, the changes in surface pCO₂, chlorophyll, alkalinity, and DIC caused by the assimilation range between about 15 % and 30 % of the mean absolute model–observation difference. However, local improvements in one variable do not necessarily come along with improvements across other observed biogeochemical variables. Therefore, globally, physics DA does not generally improve the difference between the model and the observations. In total, the effect of physics DA on the global ocean carbon uptake at 0.05 Pg C yr⁻¹ is small compared to the spread between previous estimates of models, pCO₂ products, and other DA estimates. While the assimilation of temperature and salinity improves the representation of these two and also of mixed-layer depth, sea-ice concentration, and horizontal near-surface velocities, possible errors in the vertical velocities and overturning circulation are not eliminated. Further biogeochemical variables are only indirectly affected. To this end, the additional assimilation of biogeochemical observations is an obvious next step to reduce the uncertainty stemming from the ecosystem model and to improve the model–observation differences for biogeochemical variables.

Appendix A

Table A1 List of tracers in REcoM3

Figure A1 Ensemble standard deviation of monthly air–sea CO₂ flux in the year 2020 in (a) FREE and (b) ASML.

Figure A2 North Atlantic surface velocities: (a) time-mean in ASML and (b) difference ASML − FREE.

Figure A3 Time-mean sea surface (a) temperature and (b) salinity in ASML.

Figure A4 Zonally averaged time-mean (a) temperature and (b) salinity in ASML.

Figure A5 September mean Antarctic sea-ice concentration in ASML.

Figure A6 Mixed-layer depth (a, b) in FREE and (c, d) ASML in the year 2020, (e, f) the profile-based climatology v2023 of de Boyer Montégut et al. (2004), and (g, h) the improvement through DA relative to the climatology. On the left: time-mean mixed layer. On the right: maximum of monthly mean mixed layer. For FREE and ASML (a–d), the mean absolute difference to the climatology is given in the bottom-left corner.

Figure A7 Near-surface velocities (a, b) in FREE and (c, d) ASML for the period 2010–2020, (e, f) the climatology from drifter observations of Laurindo et al. (2017a), and (g, h) the improvement through DA relative to the climatology. On the left: zonal velocities. On the right: meridional velocities. For FREE and ASML (a–d), the mean absolute difference to the climatology is given in the bottom-left corner.

Figure A8 Global meridional overturning in (a) FREE and (b) ASML and (c) the difference ASML − FREE.

Figure A9 The net difference ASML − FREE of surface pCO₂ by latitude (panel (a), blue line) and the offline-approximated effects causing that pCO₂ difference for the period 2010–2020. Thermal effect (panel (a), red line) and non-thermal effect calculated, firstly as the residual, i.e., net-minus-thermal (panels (a) and (b), dotted light-green lines), and secondly as the sum of alkalinity and DIC effects (panels (a) and (b), solid light-green lines), and effects of alkalinity and DIC individually (panel (b), orange and dark-green lines). The shaded areas in the background indicate the zonal extent of defined biomes in the Southern Ocean: ICE_SO in light blue, SPSS_SO in blue, and STSS_SO in pink. Colors blend where the regions overlap.

Figure A10 Southern Ocean time-mean sea surface (a) temperature, (b) salinity, (c) DIC, and (d) alkalinity in ASML.

Figure A11 The difference ASML − FREE of source and sink terms for the ocean's DIC and alkalinity content integrated over 0–190 m in the Southern Ocean in the year 2020. Transport terms include advection and diffusion of DIC and alkalinity. Biological terms for DIC are the sum of photosynthesis, respiration, remineralization of dissolved organic carbon, and formation and dissolution of calcite. Biological terms for alkalinity are the sum of nitrogen assimilation and remineralization and the formation and dissolution of calcite.

Figure A12 Linear offline estimate of the dominance of the thermal versus the non-thermal effect through the assimilation of pCO₂ in the Southern Ocean and North Atlantic for the period 2010–2020.

Figure A13 The difference ASML − FREE of source and sink terms for the ocean's (a) DIC and (b) alkalinity content integrated over 0–190 m per 1° latitude in the year 2020. Transport terms include advection and diffusion of DIC and alkalinity. Meridional transport is averaged across bins of 5° latitude. In panel (b), vertical and horizontal transport are summed up for readability. Biological terms for DIC are the sum of photosynthesis, respiration, remineralization of dissolved organic carbon, and the formation and dissolution of calcite. Biological terms for alkalinity are the sum of nitrogen assimilation and remineralization and the formation and dissolution of calcite.

Figure A14 Zonally averaged winter (JJA) air–sea CO₂ flux (negative: into the ocean) in FREE, ASML, and previous estimates (Hauck et al., 2023a; Verdy and Mazloff, 2017).

Figure A15 North Atlantic time-mean sea surface (a) temperature, (b) salinity, (c) DIC, and (d) alkalinity in ASML.

Figure A16 The difference ASML − FREE of source and sink terms for the ocean's DIC and alkalinity content integrated over 0–190 m in the North Atlantic in the year 2020. Transport terms include advection and diffusion of DIC and alkalinity. Biological terms for DIC are the sum of photosynthesis, respiration, remineralization of dissolved organic carbon, and the formation and dissolution of calcite. Biological terms for alkalinity are the sum of nitrogen assimilation and remineralization and the formation and dissolution of calcite.

Figure A17 Zonally averaged DIC: (a) time-mean in ASML, (b) difference ASML − FREE, (c) difference FREE − OBS based on the GLODAP climatology (Lauvset et al., 2016), and (d) improvement respective to GLODAP.

Figure A18 Surface velocities in the Southern Ocean, (a) time-mean in ASML, and (b) difference ASML − FREE.