Variability of North Atlantic CO 2 fluxes for the 2000–2017 period estimated from atmospheric inverse analyses

. We present new estimates of the regional North Atlantic (15˚N – 80˚N) CO 2 flux for the 2000 – 2017 period using atmospheric CO 2 measurements from the NOAA long term surface site network in combination with an atmospheric carbon cycle data assimilation system (GEOSChem – LETKF). We assess the sensitivity of flux estimates to alternative ocean CO 2 prior flux distribution and to the specification of uncertainties associated with ocean fluxes., We present a new scheme to 15 characterize uncertainty in ocean prior fluxes, derived from a set of eight surface pCO 2 –based ocean flux products, which reflects uncertainties associated with measurement density and pCO 2 –interpolation methods. This scheme provides improved model performance, in comparison to fixed prior uncertainty schemes, based on metrics of model – observation differences at the network of surface sites. Long term average posterior flux estimates for the 2000 – 2017 period from our GEOSChem– LETKF analyses are -0.255±0.037 PgC y -1 for the subtropical basin (15˚N–50˚N), and -0.203±0.037 PgC y -1 for the subpolar 20 region (50˚N–80˚N, eastern boundary at 20˚E). Our basin – scale estimates of interannual variability (IAV) are 0.036±0.006 PgC y -1 and 0.034±0.009 PgC y -1 for subtropical and subpolar regions respectively. We find statistically significant trend in carbon uptake for the subtropical and subpolar North Atlantic of -0.064±0.007 and -0.063±0.008 PgC y -1 decade -1 ; these trends are of comparable magnitude to estimates from surface ocean pCO 2 –based flux products, but larger, by a factor of 3-4, than trends estimated from global ocean biogeochemistry models.


Introduction
The ocean plays a key role in the global carbon budget, accounting for 2.5±0.6 PgC y -1 of net CO2 uptake from the atmosphere during the last decade (period 2009-2019), a level equivalent to ~26% of global fossil CO2 emissions (Friedlingstein et al., 2019). The North Atlantic ocean has been identified a region of significant net oceanic CO2 uptake in a range of recent analyses , Landschützer et al., 2013, Lebehot et al., 2019, and also the location of the largest Northern Hemisphere 30 uptake of anthropogenic CO2 in recent decades (Gruber et al., 2019, Khatiwala et al., 2013, Sabine et al., 2004. Recent estimates of net air-sea CO2 fluxes derived from sea surface partial pressure CO2 measurements (pCO2) indicate net annual uptake for the North Atlantic over the past decade (2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018) with a range of 0.35-0.55 PgC y -1 ( Landschutzer et al., 2016;Rodenbeck et al., 2013;Zeng et al., 2015;Watson et al., 2020), and equivalent to about 14%-22% of the global net ocean carbon ocean sink reported for this period. Regionally aggregated air-sea CO2 fluxes over the North Atlantic basin also display 35 significant variability on interannual (Watson et al., 2009) and decadal timescales (Landschützer et al., 2016. Based on analyses of surface pCO2 measurements, variations in regional pCO2 trends were observed in the subtropical and subpolar regions, potentially associated with large-scale climate oscillations such as the North Atlantic Oscillation and the Atlantic Multi-decadal Variation (McKinley et al., 2011, Macovei et al., 2020. Devries et al. (2019) estimated a negative trend (i.e., a strengthening ocean sink) in North Atlantic CO2 uptake for the 2000-2009 period based on 40 analysis of pCO2-based estimates and ocean models. Lebehot et al. (2019) find statistically significant trends in surface ocean CO2 fugacity (fCO2) for the 1992-2014 period from both observation-based surface mapping methods and from the CMIP5 Earth System models.
Recent analyses of North Atlantic air-sea CO2 fluxes have primarily been based on 'bottom-up' methods of varying complexity which use interpolated surface ocean pCO2 distributions (derived from in situ pCO2 measurements) in combination 45 with parameterizations of air-sea gas exchange (e.g., Landschützer et al., 2013;Rödenbeck et al., 2015;Takahashi et al., 2002Takahashi et al., , 2009). Estimates of air-sea CO2 fluxes have also been derived by alternative methods such as global ocean biogeochemical models (e.g., Buitenhuis et al.,2013;Law et al., 2017), and 'top-down' methods which involve the application of inverse analyses or data assimilation methods to atmospheric and oceanic CO2 measurements (e.g., Gruber et al., 2009, Mikaloff-Fletcher et al., 2006, Gurney et al., 2003, Peylin et al., 2013. 'Top-down' analyses estimate surface CO2 fluxes by using 50 information on observed gradients in atmospheric CO2 together with atmospheric transport constraints (typically from 3 D atmospheric models) and prior information on the magnitude and associated uncertainties of surface CO2 flux distributions (Rödenbeck et al., 2003;van der Laan-Luijkx et al., 2017;Peters et al., 2005;Peylin et al., 2013;Chevallier et al., 2014;Gaubert et al., 2019).
Previous studies also note that estimates of carbon fluxes from the atmospheric inverse method are sensitive to the specification 55 of the prior flux distribution and its associated uncertainty distribution (Carouge et al., 2010;Chatterjee et al., 2013;Peylin et al., 2013); While there have been recent studies evaluating the sensitivity of land-based carbon flux estimates to specification of the prior flux and its uncertainty, there has been far less examination of ocean flux estimates from inverse studies. Several global inverse model assessments of the past decade have relied on the climatological ocean-atmosphere CO2 flux database of Takahashi et al. (2009) to specify prior ocean fluxes. In view of the limited information available on the temporal and spatial 60 variability of ocean carbon fluxes from this climatological ocean database, these inverse analyses have adopted different approaches to the specification of prior uncertainty for ocean fluxes, ranging from uncertainties derived from a separate ocean model inversion (in the case of Nassar et al., 2011), to a specified percentage of the prior flux magnitude (Feng et al., 2016, Liu et al., 2016. In this study we present a new long term estimate of North Atlantic air-sea CO2 fluxes for recent decades (period 2000-2017) 65 using atmospheric inverse methods. We focus in particular on the specification of prior ocean fluxes (including sensitivity of flux estimates to alternative prior flux distributions) and on their associated flux uncertainties. To our knowledge these influences on inverse estimates of North Atlantic CO2 flux have not been assessed previously. We use the carbon cycle data assimilation system GEOSChem-LETKF (GCL, described further in Sect. 2) which combines the global atmospheric CO2 transport model GEOS-Chem (Nassar et al., 2010) with the Localized Ensemble Transform Kalman Filter (LETKF) data 70 assimilation system (Hunt et al., 2007;Miyoshi et al., 2007;Liu et al., 2019). In recent years several new global air-sea CO2 flux products have been developed based on mappings of ocean surface pCO2 measurements (e.g., Landschutzer et al., 2016, Rodenbeck et al., 2014, Watson et al., 2020 reported in the intercomparison of Roedenbeck et al., 2015). These ocean flux distributions are frequently derived from interpolations of surface ocean pCO2 measurements from the SOCAT database (Bakker et al., 2016) together with parameterizations of air-sea gas exchange. Following recent updates, the surface 75 ocean pCO2 database SOCATv2020 (https://www.socat.info/index.php/data-access/), now includes over 28 million surface ocean carbon measurements. The SOCAT database provides a valuable resource towards the development of bottom-up estimates of ocean-atmosphere CO2 fluxes, and a compilation of these flux products is reported in the recent Global Carbon Budget (Friedlingstein et al., 2020). The increased range of global air-sea CO2 flux products available (beyond the Takahashi et al., 2009 climatology) provides a valuable opportunity to develop an improved representation of air-sea CO2 flux variability 80 and a more robust characterization of the uncertainties associated with ocean carbon fluxes. In this study we employ some of the recently developed ocean CO2 flux products to provide a new method of characterizing the prior ocean flux uncertainty used for atmospheric inverse analyses. The methodology is based on the ensemble spread of the multiple ocean flux products, and reflects underlying uncertainties in these products, such as those associated with sampling density of the surface measurements and interpolation method employed. It provides a spatially and temporally variable specification of prior flux 85 uncertainty that will be of value to the inverse modeling community.
The remainder of the paper is organized as follows: Section 2 covers the methodology of the atmospheric inverse analysis, outlining the carbon cycle data assimilation system (GEOSChem-LETKF), the atmospheric CO2 observations, and specifications of prior fluxes and uncertainties. Further details of the methodology are presented in the Appendix. In Sect. 3 we present GEOSChem-LETKF assessments of alternative specifications of ocean prior fluxes and flux uncertainties, and 90 then use these results to derive long term estimates of North Atlantic CO2 fluxes for the 2000-2017 period. We also summarize specific characteristics of North Atlantic CO2 fluxes derived from these analyses, namely, long term means, trends, and interannual variability of fluxes, and compare our results with other recent relevant studies.

Overview 95
Our analysis employs the global GEOS-Chem atmospheric chemistry transport model together with the Localised Ensemble Kalman Filter (LETKF) (described in Sect. 2.2) and atmospheric CO2 observations from the NOAA-ESRL network of surface sites (Sect. 2.3). Section 2.4 describes the compilation of the set of air-sea CO2 flux products and the derivation of the prior flux uncertainty specification for the North Atlantic based on the ensemble spread of these products. Section 3 presents model results, including sensitivity analyses assessing different prior flux representations and flux uncertainty schemes (Sect. 3.1), 100 and regional CO2 flux estimates for the 2000-2017 period from the GEOSChem-LETKF system (Sect. 3.2). Further details on the methods, model, observations and uncertainty calculations are presented in the sections below and in the Appendix.

The GEOSChem-Localized Ensemble Kalman Filter (GCL) system
The GEOS-Chem atmospheric chemistry transport model has been used in a range of previous investigations into atmospheric CO2 and applied in conjunction with inverse analyses to estimate surface carbon fluxes (Nassar et al., 2010(Nassar et al., , 2011105 Suntharalingam et al., 2005;Liu et al., 2016). In this analysis we employ GEOSChem v11-01 at a horizontal resolution of 2° latitude by 2.5° longitude, with 47 levels in the vertical. Model transport fields are provided by GEOS-5 assimilated meteorological data from the NASA Global Modeling and Assimilation Office (GMAO, Rienecker et al., 2008). The GEOSChem configuration employed here primarily follows that of Nassar et al. (2011), but with updated representation of prior fluxes; more detail on the prior CO2 fluxes and uncertainties implemented in this study is given in Sect. 2.4. 110 The Localized Ensemble Transform Kalman Filter (LETKF) is a data assimilation system which provides an estimate given a prior (or "background") estimate of the current state based on past and current data (in this case, the atmospheric CO2 mole fraction observations). The general framework of the LETKF is described in Hunt et al. (2007); it has been adapted by Miyoshi et al. (2007) to provide gridscale localized analysis of flux estimates. The LETKF system has been used to estimate CO2 fluxes in a range of previous studies (e.g, Kang et al., 2012;Liu et al., 2016Liu et al., , 2019. The LETKF provides iterative estimates of the 115 time evolution of the system state, x, (here representing the grid-scale surface carbon fluxes). Each step involves a forecast stage (based on a physical model of the system evolution) and a state estimation stage (the 'analysis' step), which combines system observations, y, together with the background forecast, , to derive the improved state estimate. The observation operator H provides the mapping from the state space to the observation space; in this study H is provided by the GEOS-Chem atmospheric model. 120 In this analysis we employ the complete GEOSChem-LETKF (GCL) data assimilation system to conduct sensitivity analyses on the ocean prior fluxes, and to provide a long term flux estimate of surface CO2 fluxes for the North Atlantic for the period 2000-2017. We report a posteriori fluxes on monthly timescales for the 2000-2017 period; the optimized monthly fluxes are derived from four sequential weeks of assimilation cycles, as further described below. Our methods follow the implementation of the LETKF system by Liu et al. (2019), who have extended the previous carbon data assimilation system of Kang et al. (2011Kang et al. ( , 2012. The study of Kang et al. (2011) assimilated meteorological data and atmospheric CO2 concentrations to provide estimated atmospheric CO2 concentrations as part of the state estimate. Kang et al. (2012) extended this method to also provide estimates of surface carbon fluxes. Both these LETKF studies assimilated meteorological data and atmospheric CO2 concentrations and employed a short assimilation window of 6 hours in order to maintain linear behaviour of the ensemble perturbations (Kang et al., 2011(Kang et al., , 2012. In addition, Kang et al. (2012) also tested longer assimilation windows (up to 3 weeks) 130 for LETKF formulations that assimilated atmospheric CO2 concentrations alone (eliminating the assimilation of the meteorological data). The LETKF system of Liu et al. (2019) extended the Kang et al. (2011Kang et al. ( , 2012 analyses by incorporating the GEOSChem atmospheric model as the forecast model, along with its representation of surface CO2 fluxes which provide the prior flux specification for the forecast step. However, Liu et al. (2019) assimilate only atmospheric CO2 measurements (i.e., no assimilation of meteorological measurements), and use an assimilation window of 7 days; the duration of the 135 assimilation window was selected to maximize the correlation between observations and surface fluxes. The GEOSChem-LETKF system employed in our study follows the Liu et al. (2019) formulation; atmospheric CO2 measurements are assimilated at 7 day timescales, with the LETKF analysis step providing updates of the surface fluxes and associated uncertainties required as initial conditions for the next weekly forecast step. We report monthly flux estimates following four assimilation cycles. Further details on the LETKF and the governing equations for flux estimation are provided in Appendix 140 A.

Atmospheric CO2 Observations
Atmospheric CO2 observations used for this study are taken from the NOAA-ESRL GLOBALVIEWplus Observation  Table A1 of the Appendix. The specification of observational uncertainty associated with incorporation of the atmospheric CO2 measurements into the LETKF is derived using the methods of Chevallier et al. (2010); we use the standard deviation of measurement variability from detrended and deseasonalized CO2 time series at each measurement site. The resulting specification of observational uncertainty varies between 0.16 ppm (for 150 stations in and around the Southern Ocean) to over 5 ppm (for stations in continental interiors) (see Appendix Table A1 for more details).

Specification of Prior CO2 Fluxes and Associated Flux Uncertainties
The GEOSChem model CO2 simulation employed in this study includes representation of fossil fuel emissions, air-sea fluxes  Quéré et al., 2018), and land biosphere fluxes from the Joint UK Land Environment Simulator (JULES, Clark et al., 2011).
The focus of our study is on North Atlantic Ocean CO2 fluxes, and we investigate the representation of ocean prior fluxes and prior flux uncertainty in more detail. Firstly, in Sect. 3.1, in a set of sensitivity analyses, we compare the implementation of 160 three different representations of ocean CO2 fluxes that have been used to specify prior fluxes in recent inverse analyses: (i) the widely used Takahashi et al. (2009)  employed relatively simple characterizations of the prior ocean flux uncertainty, for example, based on a fixed proportion of the grid-scale or regional prior flux (Nassar et al., 2011, Liu et al., 2016, Feng et al., 2016. In Sect. 3.1, we employ both fixed flux uncertainties, and also present an alternative scheme derived from the ensemble spread of ocean CO2 flux products, as described below. The prior ocean flux distributions employed in atmospheric inversions are frequently derived from interpolations of the surface 170 ocean pCO2 database (e.g., SOCAT, Bakker et al., 2016) in combination with ocean-atmosphere gas exchange parameterizations. Uncertainties in the derived products stem from uncertainties in the input data (e.g., density of measurements), interpolation methods, and gas-transfer parameterizations (Landschutzer et al., 2013). However, some ocean regions, the North Atlantic in particular, have a higher density of pCO2 measurements and more consistent flux estimates from pCO2-based products , Landschutzer et al., 2013. Here we exploit the recent expansion of pCO2-based 175 ocean flux products to outline a new specification of ocean prior flux uncertainty based on the ensemble-spread of the different flux products (the "spread-based" uncertainty scheme). Towards the development of the spread-based scheme, we have compiled a set of eight global gridded interannually varying ocean-atmosphere CO2 flux products. These are Landschutzer et al., 2016, Rodenbeck et al., 2014, Denvil-Sommer et al., 2019, Zeng et al., 2015, Gregor et al., 2019, and Watson et al., 2020 The spread-based prior flux uncertainty scheme uses a diagnostic derived from the variation among the set of ocean atmosphere carbon flux products (see Eq. (1)). This scheme specifies lower uncertainty levels where alternative prior flux representations are in accord (e.g., when well-constrained by availability of surface pCO2 measurements), and higher uncertainty levels where the prior flux distributions differ significantly (typically in under-sampled regions or those of significant flux variability). This specification follows previously used methods to characterize uncertainties in ocean flux 185 distributions (e.g., Bopp et al., 2013). For this spread-based uncertainty specification, the gridded prior flux uncertainty, U(i,j) (for a gridcell with coordinates (i,j)) is specified as the standard deviation of the spread of the different prior flux products.
Thus, the uncertainty U(i,j) is calculated as: Here K is the total number of the prior ocean flux products considered, and subscript k refers to an individual flux product.
( , ) represents the gridded monthly flux for each prior ocean flux and ( , ) is the gridded monthly mean across all prior 190 ocean flux products. These prior flux uncertainties are estimated on monthly timescales and also account for interannual variations. The uncertainty statistics of the prior ocean flux distributions will be dependent on the uncertainties associated with the respective inputs and methods of constructing the flux products. Ocean-atmosphere carbon flux products derived from surface ocean pCO2 measurements are generally subject to two main sources of uncertainty: (i) in the specification of the surface CO2 partial pressure difference across the air-sea interface, and (ii) in the specification of the gas-exchange coefficient 195 used to derive fluxes (e.g., see discussion of Landschutzer et al., 2013;Watson et al., 2020). In the extended database of 8 pCO2-based flux products that we present above, the majority of the flux products (seven of the eight) rely on the surface ocean pCO2 data of the SOCAT database (Bakker et al., 2016). These flux products will be subject to similar uncertainties associated with data coverage in different ocean regions, although the uncertainties due to differences among surface interpolation methods may vary. 200 In this study we account for spatial correlations in the prior ocean fluxes, by inclusion of off diagonal elements in the background error covariance matrix P b (Appendix A1 Eq. A3). We follow the recommendations of Jones et al. (2012) on autocorrelation length scales in the surface ocean. That study derived spatial autocorrelation functions for air-sea fluxes from an analysis of the surface ocean pCO2 database reported in Takahashi et al. (2009), combined with a gas exchange parameterization. We currently do not account for spatial correlation in land fluxes, but will investigate this in future analyses. 205

Results and Discussion
We first present in Sect. 3.1 results of short-term sensitivity tests that compare the influence of different prior ocean flux distributions and prior ocean flux uncertainty schemes on GCL estimates of North Atlantic (NA) CO2 fluxes. Using these analyses as a basis, in Sect. 3.2 we conduct a multi-year GCL analysis of North Atlantic CO2 ocean fluxes for the 2000-2017 period. We also report on derived characteristics of regionally aggregated North Atlantic subtropical and subpolar fluxes (long 210 term means, trends and interannual variability) and compare these GCL results with recent estimates from other methodologies, including global ocean biogeochemical models (GOBMs), other atmospheric inverse studies, and surface pCO2-based data products.

Sensitivity tests on specification of prior flux uncertainty
In this section we investigate, via sensitivity analyses, the application of the spread-based prior flux uncertainty scheme 215 CO2 fluxes. The alternative specifications of prior flux uncertainty for ocean fluxes employed include (a) fixed percentagebased levels (U1:60% of prior flux, and U2:120% of prior flux), and (b) gridded flux uncertainties representing the variation or 'spread' of the different ocean flux data products at each location, and based on the standard deviation of the variation among the prior fluxes (U3: spread-based uncertainty; see Eq. (1)). The selection of the fixed percentage prior uncertainty 220 levels used in the sensitivity analyses was based on the range of variability seen for the individual prior flux distributions ( Fig.   1) for the sub-regions of the North Atlantic. These ranged from average levels of ~60% for the subtropical North Atlantic to levels greater than 120% for the subpolar North Atlantic, hence we have selected a level of U1:60% to characterize the lower sensitivity case, and U2:120% for the higher case. We apply the alternative flux uncertainty specifications to the three different ocean prior flux distributions discussed in Sect. 2.4, namely: (i) the Takahashi  2000; the length of spin-up was determined by recommendations on the duration required for stabilization of tropospheric CO2 gradients (e.g., Gurney et al., 2002), and following methods used for previous GEOSChem CO2 analyses (e.g., Nassar et al., 2010). The year 2003 was selected for sensitivity tests as the first viable year following spin-up. Analyses of inter-annual 230 variability in Atlantic CO2 (e.g., Landschutzer et al., 2013;Schuster et al., 2013) do not find 2003 to be an anomalous year for regional ocean fluxes. We evaluate the sensitivity of posterior ocean flux estimates with three different prior ocean uncertainty schemes U1, U2, and U3, described above; these are applied in turn for each of the three prior ocean flux distributions (Ta, La and Ro). Figure 1 presents the seasonal variation of the spatial distribution of the spread-based prior ocean flux uncertainty U3 (3 month averages for the year of 2003). Figure 1 demonstrates that over the course of the year, and particularly in the 235 Northern Hemisphere winter months, the spread-based uncertainty scheme (U3) provides a looser constraint on prior fluxes (i.e., levels of prior flux uncertainty > 120%) than the U1 and U2 schemes in the subpolar region, and a tighter constraint in the subtropical region (levels < 60%).  Table 1 summarizes the prior and posterior ocean flux estimates for the global and North Atlantic region (sub-divided into subpolar and subtropical regions) from the respective sensitivity tests. The distribution of prior flux for the subtropical North Atlantic shows closer agreement among the three source representations (Ta, La and Ro), with regional variation of 0.05 PgC 250 y -1 , in comparison to a regional variation of ~0.1 PgC y -1 for the subpolar region. Under the constraints provided by the atmospheric CO2 observations all posterior flux estimates for the North Atlantic show increased uptake (Table 1), indicating that all three representations of ocean prior flux underestimate the regional net atmosphere-ocean flux for the 2003 period.
Largest changes in the regional posterior fluxes are estimated under the U3 specification of prior flux uncertainty. In addition, our estimates indicate a larger increase in CO2 uptake in the subpolar basin (~0.05 PgC y -1 , changing from a prior flux range 255 of -0.13 to -0.23 PgC y -1 to posterior flux range of -0.18 to -0.27 PgC y -1 , for the U3 scenarios), in comparison to the smaller magnitude change for the subtropical North Atlantic basin (of ~ 0.04 PgC y -1 from -0.17 to -0.22 PgC y -1 to -0.22 to -0.26 PgC y -1 for the U3 scenarios) We note that the increases in estimated uptake for the North Atlantic basins are relatively smaller (on average in the range 10-20%) than the increased uptake estimated on the global scale (~30-50% changes, see Table 1), indicating that prior flux 265 representations of North Atlantic carbon uptake are more consistent with the constraints from atmospheric CO2 measurements than for other regions of the global ocean.
The U3 flux uncertainty specification is derived from the variation among a set of ocean-atmosphere carbon flux products (Eq. (1)). This scheme specifies lower uncertainty levels where alternative prior flux representations are in accord (e.g., when well constrained by availability of surface pCO2 measurements, as in the subtropical North Atlantic), and higher uncertainty 270 levels where the prior flux distributions differ significantly (typically in under-sampled regions or those of significant flux variability, such as the subpolar North Atlantic). We further assess the value of the U3 scheme using a metric of GCL modeled atmospheric CO2 concentration; specifically, estimates of the model-observation mismatch for the year 2003 at the NOAA network station sites in the North Atlantic using the a posteriori fluxes associated with the sensitivity analyses of this section (Appendix Table A2). The results summarized in Table A2 indicates that scheme U3 provides the smallest magnitude model-275 observation mismatch for the individual North Atlantic sites and for the global network average. For the long term analyses of the remainder of this study, therefore, we use the U3 spread-based flux uncertainty scheme in preference to the fixed level flux uncertainty schemes used in many previous inverse analyses.

Multi-year analyses of North Atlantic CO2 fluxes
In this section we present results of a multi-year GCL analysis (for the period 2000-2017), calculating regional estimates of   Atlantic subpolar region are in closer accord (Fig. 2b) with differences of less than 0.2 PgC y -1 (the CT estimate is an exception indicating variations of greater than 0.3 PgC y -1 from the other estimates). A potential reason for the anomalous behaviour of the CT estimate in the North Atlantic is the underlying prior flux uncertainties used in the analysis which give a loose constraint on the prior ocean fluxes and allow the ocean fluxes deviate far from the prior because of the impaction from atmospheric CO2

signals. 310
We also note that Peylin et al. (2013) Table 2); (ii) the estimated inter-annual variability (IAV) of fluxes (Table 3); and (iii) the long term trends ( Table 4). The IAV is calculated following methods of Rödenbeck et al. (2015) (i.e., derived from the standard deviation 325 of the residuals of a 12 month running mean over the CO2 flux time series). trend (panels (e) and (f)). The GCL estimates (red stars) are shown in comparison to other atmospheric inverse analyses (red symbols), surface ocean pCO2 products (blue) and global ocean biogeochemistry models (GOBMs, purple). Also shown are the estimated mean values from each sub-group of analyses (filled cross symbols) with their minimum-maximum range. Circled symbols in panel (e) and (f) indicate a statistically significant trend. Table 2. Summary metrics of GEOSChem-LETKF North Atlantic (NA) CO2 flux estimates, and comparison with independent estimates (from atmospheric inverse analyses, surface pCO2 mappings, and Global Ocean Biogeochemistry models (GOBMs)) for the period 2000-345 2017. Listed are estimates for the long term mean. The metrics listed in this table are plotted in Fig. 3a and 3b.
Long term mean (PgC y -1 ) NA Subtropics We also present in Fig. 3 the equivalent estimates from other independent assessments, including (i) other atmospheric inverse 350 analyses, (ii) surface ocean pCO2-based analyses, and (iii) analyses from global ocean biogeochemistry models (GOBMs).
For the North Atlantic subtropical region, the long term mean of the GCL posterior flux estimate is -0.255±0.037 PgC y -1 ( Figure 3a and Table 2). It lies in the range spanned by the other inverse analyses (-0.31 to -0.20 PgC y -1 , of slightly larger magnitude than CAMS, but smaller than CT and CTE), and is in good agreement with the mean-value estimates from surface pCO2-based methods and GOBMs. For the North Atlantic subpolar region, the GCL estimate of the long term mean uptake 355 is -0.203±0.037 PgC y -1 (Table 2), close to the inverse estimate of the CAMS analysis, and of smaller magnitude (by ~0.1 PgC y -1 ) than the inverse estimates of CT and CTE. The GCL estimate is consistent with the mean estimate from pCO2-based products, and within the range of flux estimates from GOBMs (from -0.341 to -0.197 PgC y -1 ).

Interannual variability 360
The interannual variability (IAV) of CO2 flux estimates derived from the GCL is 0.036±0.006 PgC y -1 for the North Atlantic subtropics and 0.034±0.009 PgC y -1 for the North Atlantic subpolar region ( Fig. 3c and 3d, Table 3). The IAV estimates from the different inverse analyses for both the North Atlantic subtropics and subpolar regions display a large range (0.032 to 0.084 PgC y -1 and 0.023 to 0.114 PgC y -1 respectively), than the ranges displayed by GOBMs (0.014 to 0.027 PgC y -1 and 0.015 to 0.024 PgC y -1 respectively) and pCO2-based fluxes (0.029 to 0.050 PgC y -1 and 0.009 to 0.037 PgC y -1 respectively). The 365 larger range of IAV from atmospheric inverse analyses is influenced especially by high magnitude IAV estimates from the CarbonTracker (CT) inverse analysis. Potential causes of the differences among atmospheric inversion between the GCL and CAMS IAV estimates and those of the CarbonTracker estimates are the different prior ocean fluxes employed by the inverse analyses, and the relative weighting assigned to the influence of atmospheric CO2 observations (Jacobson et al., 2019). The  Hauck et al. 2020). An important driver of IAV is the variability in biological carbon export; the lower variability observed in the GOBMs could result from opposing changes in biological vs. circulation impacts on carbon export, which potentially reduces the sensitivity of the 375 GOBM air-sea carbon fluxes to climate variability (Landschutzer et al., 2013, DeVries et al., 2019. Table 3. Summary metrics of GEOSChem-LETKF North Atlantic (NA) CO2 flux estimates, and comparison with independent estimates (from atmospheric inverse analyses, surface pCO2 mappings, and Global Ocean Biogeochemistry models (GOBMs)) for the period 2000-2017. Listed are estimates for the interannual variability of the regional fluxes over the period. The metrics listed in this table are plotted in Fig. 3 c, d.

385
Interannual Variability (IAV) (PgC y -1 ) NA Subtropics  Table 4 and Fig. 3e and 3f. We also highlight 390 in the table and figure panels, the trend estimates that are statistically significant (significant at the 95% confidence level, Montgomery et al., 2012). Our GCL analyses indicate statistically significant trends for the 2000-2017 period of -0.064±0.007 PgC y -1 decade -1 in the North Atlantic subtropical basin, and 0.063±0.008 PgC y -1 decade -1 in the subpolar region. These estimated trends are of similar magnitude to those estimated from surface ocean pCO2 products, but of much larger magnitude (by a factor of 3-4) than decadal trends estimated from the GOBMs (Fig. 3e, Table 4). Our findings are similar to those of 395 Devries et al. (2019), who noted that decadal trend estimates of North Atlantic CO2 uptake for the 2000s from the SOCOM inter-comparison of pCO2-based flux products were larger than those from the GOBMS in their analysis (see Fig. 3 of their study). Table 4. Summary metrics of GEOSChem-LETKF North Atlantic (NA) CO2 flux estimates, and comparison with independent estimates 420 (from atmospheric inverse analyses, surface pCO2 mappings, and Global Ocean Biogeochemistry models (GOBMs)) for the period 2000-2017. Listed are estimates for the trend of the regional fluxes over the period. The metrics listed in this table are plotted in Fig. 3 e, f of the main study.
Trend (PgC y -1 decade -1 ) NA Subtropics In this study we present a new long term estimate of North Atlantic air-sea CO2 fluxes for recent decades (period 2000-2017) using the atmospheric carbon cycle data assimilation system GEOSChem-LETKF. We focus, in particular, on the specification 430 of prior ocean fluxes, including sensitivity of flux estimates to alternative prior flux distributions, and on the specification of uncertainties associated with ocean fluxes. Towards this we have developed the 'spread-based' flux uncertainty scheme which represents the variability among a set of different prior ocean CO2 flux representations. The scheme ascribes higher levels of uncertainty to regions with larger discrepancies among ocean CO2 prior flux representation that arise from uncertainties associated with measurement density and pCO2-interpolation methods (Sect. 2.4). The spread-based flux uncertainty scheme 435 provides improved performance in comparison to schemes with fixed prior flux uncertainty levels, based on an assessment metric of differences in model-observation values for atmospheric CO2 at North Atlantic measurement sites of the NOAA-GLOBALVIEWCO2 network (Sect. 3.1). It provides a valuable new means of specifying prior flux uncertainties for atmospheric inverse analyses of ocean CO2 fluxes.
We have used the spread-based flux uncertainty scheme in the GEOSChem-LETKF to derive estimates of CO2 fluxes in the 440 North Atlantic for the 2000-2017 period. Long term mean estimates of the regional ocean CO2 uptake are -0.255±0.037 PgC y -1 for the North Atlantic subtropics and -0.203±0.037 PgC y -1 for North Atlantic subpolar region., and are consistent with recent regional flux estimates from surface pCO2-based methods and global ocean biogeochemistry models (GOBMs). The GEOSChem-LETKF estimates of interannual variability in air-sea CO2 fluxes are 0.036±0.006 PgC y -1 (North Atlantic subtropics) and 0.034±0.009 PgC y -1 (North Atlantic subpolar). In common with estimates from other atmospheric CO2 inverse 445 studies, the magnitude of IAV derived from the GEOSChem-LETKF is larger than corresponding estimates from GOBMs.
Our GEOSChem-LETKF estimates also indicate statistically significant trends of increasing CO2 uptake for the North Atlantic subtropical and subpolar region (estimated trend of -0.064±0.007 and -0.063±0.008 PgC y -1 decade -1 respectively). These trends are of comparable magnitude to those estimated from surface pCO2-based flux products, but much larger than those derived from global ocean biogeochemistry models for the 2000-2017 period. Estimates of inter-annual variability and long 450 term trends derived from our GEOSChem-LETKF analyses are generally more robust for the North Atlantic subtropics than for the NA subpolar region, and characterized by smaller uncertainty bounds. Limiting factors affecting estimates for the North Atlantic subpolar region include higher levels of uncertainty associated with specification of prior fluxes (Fig. 1), and the observational uncertainty at the atmospheric measurement CO2 sites in these high northern latitudes (Table A1). The number of regional atmospheric CO2 measurement sites available to constrain North Atlantic subpolar fluxes are also relatively few in 455 comparison to the subtropical region. Improved ocean CO2 flux estimates and associated metrics for this North Atlantic region will be obtained by provision of additional high accuracy marine boundary layer CO2 measurements for the region from fixed surface sites and from ships and buoys (Wanninkhof et al., 2019).
Appendix A: The Local Ensemble Transform Kalman Filter (LETKF) system 460 Here we briefly describe the LETKF system used for estimation of surface CO2 fluxes. The methodology follows that of Hunt et al. (2007) and Miyoshi et al. (2007), and additional detail is provided in these publications. The LETKF has been previously used in meteorological forecasting, and more recently in atmospheric CO2 data assimilation (e.g., Liu et al. 2019Liu et al. , 2016Kang et al. 2012). The LETKF provides iterative estimates of the time evolution of the system state, x, (here representing the gridscale surface carbon fluxes, of dimension m). Each step involves a forecast stage (based on a physical model of the system 465 evolution) and a state estimation stage (the 'analysis' step), which combines system observations, y (of dimension n), together with the background forecast, , to derive the improved state estimate. The observation operator H provides the mapping from the state space to the observation space; in this study H is provided by the GEOSChem atmospheric model. In the analysis step, the surface carbon flux estimates are obtained by minimization of a cost function (Eq. A1) which accounts for deviations of the system state x, from the background forecast, , and for the mismatch between observations (y) and their modeled 470 representations (Hx) : B represents the background flux covariance matrix, and R represents the observation covariance matrix.
In the LETKF system, an ensemble of model simulations is used to calculate the sample mean and covariance of the system state; thus, the background state is given by ( ( ) : = 1.2, … ) for k ensemble members. The sample mean ̅ and covariance of the background state vector given by : 475 is an m×k matrix whose ith column is ( ) − ̅ . is the background state covariance matrix (m×m).
Similarly the analysis state is represented by the ensemble ( ( ) : = 1,2, … ) with its sample mean and covariance given by: where is the m×k matrix whose ith column is ( ) − ̅ .
The analysis state and covariance ̅ and are updated based on the background information ̅ and observations y through the following equations: 480 The ensemble ( ) of background observation vectors is defined by: where is the n×k matrix whose ith column is ( ( ) − ), and w is a Gaussian random vector with mean = 0 and covariance = ( − 1) .Then the analogues of analysis equations (6) and (7) are: Following Hunt et al. (2007) and Miyoshi et al. (2007) (refer to these publications for the complete LETKF derivation) the overall analysis equation is: 485 The LETKF allows flexibility in the choice of observations to be assimilated at each grid point, based on the distance (r) of the observations from the gridpoint. The localization weighting function f(r) is given by: where L is an observation localization length which can be predefined to determine the outer boundary of the influence of the observations; i.e., the localization weighting function drops to zero at a value of = 2. 10 3

(A14)
The observation localization is realized by multiplying the inverse of the localization function f(r) with the observational error 490 covariance R.
The parameter L represents the horizontal localization radius, and is set to 2000 km for this study, following Liu et al. (2016).
The localization radius is used in the LETKF in a latitude-dependent weighting function which characterizes the spatial scale of the region within which atmospheric CO2 observations are assimilated at each gridpoint (Miyoshi et al. 2007).   Author contributions. ZC and PS designed the study. ZC, PS, JZ and NZ developed the model. ZC, PS, AW, and US discussed 530 the design of simulations. ZC performed the simulations and analysis and wrote the initial manuscript. All authors contributed to the writing of the paper.
Competing interests. The authors declare that they have no conflict of interest.
Disclaimer. The work reflects only the author's view,the Europen Commission and their executive agency are not responsible for any use that may be made of the information the work contains. 535