Progress in Earth andPlanetary Science

Niwa et al.
Progress in Earth and Planetary Science(2022) 9:42
https://doi.org/10.1186/s40645-022-00502-6

Open Access

METHODOLOGY

Toward a long‑term atmospheric ­CO2
inversion for elucidating natural carbon fluxes:
technical notes of NISMON‑CO2 v2021.1
Yosuke Niwa1,2* , Kentaro Ishijima2, Akihiko Ito1,3,4 and Yosuke Iida2,5

Abstract
Accurate estimates of the carbon dioxide ­(CO2) fluxes at the earth’s surface are imperative for comprehending thecarbon cycle mechanisms and providing reliable global warming predictions. Furthermore, they can also provide
valuable science-based information that will be helpful in reducing human-induced C­ O2 emissions. Inverse analysis
is a prominent method of quantitatively estimating spatiotemporal variations in C
­ O2 fluxes; however, it involves acertain level of uncertainty and requires technical refinement, specifically to improve the horizontal resolution sothat local fluxes can be compared with other estimates made at the regional or national level. In this study, a novel
set of inversion schemes was incorporated into a state-of-the-art inverse analysis system named NISMON-CO2. The
introduced schemes include a grid conversion, observational weighting, and anisotropic prior error covariance, thedetails of which are described. Moreover, pseudo-observation experiments were performed to examine the effect of
the new schemes and to assess the reliability of NISMON-CO2 for long-term analysis with practical inhomogeneous
observations. The experiment results evidently demonstrate the advantages of the grid conversion scheme for highresolution flux estimates (1° × 1°), with notable improvements being achieved through the observational weighting
and anisotropic prior error covariance. Furthermore, the estimated seasonal and interannual variations in regional C­ O2
fluxes were confirmed to be reliable, although some potential bias in terms of global land–ocean partitioning wasobserved. Thus, these results are useful for interpreting the flux variations that result from real-observation inverse
analysis by NISMON-CO2 ver. 2021.1.
Keywords: Inverse analysis, Carbon dioxide, Carbon cycle, Atmospheric transport model, Data assimilation
1 Introduction
Assessing the carbon dioxide ­(CO2) fluxes at the earth’s
surface, between the atmosphere–land and atmosphere–
oceans, is essential for studying the global carbon cycle.
Within the atmosphere, the mole fraction of C
­ O2 hasdistinct seasonal variations caused by natural ­CO2 fluxesthat are superimposed with the human-induced increasing trend. These natural fluxes are induced by dynamic
photosynthesis and respiratory activity in the terrestrial*Correspondence: niwa.yosuke@nies.go.jp
1
National Institute for Environmental Studies, 16‑2 Onogawa, Tsukuba, Ibaraki
305‑8506, Japan
Full list of author information is available at the end of the article

biosphere and by the sea–air C­ O2 exchange driven byocean thermodynamic activity, oceanic surface mixing, and biological activity. In addition to seasonal variations, these natural fluxes exhibit notable interannual
variations in response to major climate events, such as
El Niño (e.g., Keeling et al. 1995). Changes in physical or
biological land and ocean conditions can further influence the ­CO2 flux trends on the decadal scales. However,
the mechanisms underlying the natural C­ O2 flux variations are still unclear, resulting in considerable uncertainties in global warming predictions made by earth systemmodels equipped with climate–carbon cycle feedbacks(Friedlingstein et al. 2006).

© The Author(s) 2022. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which
permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to theoriginal author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images orother third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line
to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory
regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of thislicence, visit http://​creat​iveco​mmons.​org/​licen​ses/​by/4.​0/.

Niwa et al. Progress in Earth and Planetary Science(2022) 9:42

An inverse analysis of atmospheric ­CO2, where an
atmospheric transport model is used to connect the surface fluxes to the observed mole fraction of atmospheric­CO2, is a prominent method that provides quantitative observation-based estimates of spatiotemporal surface flux variations. In fact, a number of ground-based
stations and other observational platforms have beenused to measure atmospheric ­CO2 for decades; throughinverse analyses, they have revealed considerable interannual variations in region-by-region ­CO2 fluxes (Bousquet et al. 2000; Patra et al. 2005a,b; Baker et al. 2006;
Chevallier et al. 2010; Peylin et al. 2013). Over the last
decade, the four-dimensional variational method (4DVar) and ensemble Kalman smoother (or filter), which
are well-known state-of-the-art data assimilation methods utilized in meteorological and oceanic research,
have become common inverse analysis methods andhave made use of a wide variety of observations (Peters
et al. 2007; Chevallier et al. 2010; van der Laan-Luijkx
et al. 2017; Niwa et al. 2017b; Liu et al. 2021). However,Friedlingstein et al. (2020) demonstrated that inverse
analyses have significantly large ranges of interannual
flux variation estimates, even for zonal means of certain
latitudinal ranges, further verifying the limited reliabilityof such inverse analyses.
Recently, societal momentum to reduce C­ O2 emissions has gained traction globally in response to the
1.5 °C/2.0 °C target of the Paris agreement. Consequently,­CO2 flux data from inverse analyses are attracting global
attention because of their potential to provide useful
information for confirming national inventories. Recent
studies by Chevallier (2021) and Deng et al. (2022) compared inversion fluxes with national inventories basedon agriculture, forestry, and other land-use changes(LUCs). However, they used data from global-scale inversions with horizontal resolutions around hundreds ofkilometers. Therefore, further development is required
to increase the horizontal resolution so that fluxes can becompared more consistently with national inventories.
This is a cumbersome task, because a several-decadeslong analysis is required to investigate long-term trends,
which is computationally demanding.
Niwa et al. (2017a, b) developed a sophisticated andcomputationally efficient inversion system based on
the Nonhydrostatic Icosahedral Atmospheric Model(NICAM: Tomita and Satoh 2004; Satoh et al. 2008,2014). The ­CO2 inverse analysis using this model istermed the NICAM-based Inverse Simulation for Monitoring ­CO2 (NISMON-CO2) (Niwa et al. 2021), and the­CO2 flux estimates made using this system were included
in a multi-inversion analysis of the Global Carbon Project (GCP) (Friedlingstein et al. 2020, 2022). The generalfeatures of NISMON-CO2 were described by Niwa et al.

Page 2 of 19(2017b), who performed a pseudo-observation experiment. However, the experiment was performed for only
a year under ideal conditions; therefore, the practical performance of NISMON-CO2 for long-term analysis is yet
to be demonstrated.
In this study, we describe the latest inversion settingsof NISMON-CO2, specifically those intended to improve
the inversion performance. Furthermore, we performed
a long-term inverse analysis using pseudo-observations
that emulate actual observation networks. We focusedon the use of NISMON-CO2 in investigating interannual variations in regional ­CO2 fluxes and on the effectiveness of the newly introduced techniques. Apart fromthis study, an inverse analysis was performed using NISMON-CO2 with the same inversion settings, but instead
using actual observations for 1990–2020, whose data are
publicly available as NISMON-CO2 ver. 2021.1 (Niwa
2020). This study also aims to present a detailed description of NISMON-CO2 ver. 2021.1.

2 Methods
2.1 NICAM‑TM and NISMON‑CO2

Initially, we describe the fundamental features of the
NISMON-CO2 inversion system. One distinct featureof NICAM, on which NISMON-CO2 is based, is that it
adopts an icosahedral grid system and calculates the continuity equation using the finite volume method (Tomita
and Satoh 2004). Thanks to the finite volume method, theNICAM-based transport model (NICAM-TM) can completely conserve tracer mass without any numerical mass
fixer, which motivates the usage of NICAM for simulating atmospheric transport of long-lived species such as­CO2 (Niwa et al. 2011a, b). Simulations with NICAMTM are performed iteratively in an inverse analysis of
NISMON-CO2. Further details about NICAM-TM aredescribed by Niwa et al. (2017a).
The inverse analysis system NISMON-CO2 uses
4D-Var to optimize surface ­CO2 fluxes consistently withatmospheric observations (Niwa et al. 2017b). In 4D-Var,the cost function J, whose minimum state is sought, isdefined as:
J (δx) =1 T −11
δx B δx + (M(δx) − d)T R−1 (M(δx) − d), (1)22where δx represents the model parameter to be optimized, which comprises deviations from its prior estimate xpri of surface fluxes and the initial global mole
fraction offset. Meanwhile, M(.) is an operator representing the model calculation, and d is a vector definedby the observational vector y as d = y − M x pri . The
matrices B and R are the error covariance matrices ofprior estimates and the model–observation mismatches,

Niwa et al. Progress in Earth and Planetary Science(2022) 9:42

Page 3 of 19

respectively. The iterative calculation of 4D-Var uses thegradient of the cost function:∇δx J = B−1 δx + MT R−1 (M(δx) − d)(2)

to find the optimal model parameters (i.e., posterior estimates) that minimize the cost function. MT
in Eq. (2) represents the transpose of M, where M isthe tangent linear operator of M(.). In fact, the relationM�x ≈ M(�x) is satisfied because the model calculates
atmospheric transport linearly and does not include any
chemical reactions. The calculations of M(.) and MT are
performed by the forward and adjoint simulations using
NICAM-TM, which only simulates ­CO2 transport in theso-called offline mode with prescribed meteorological
data. The meteorological data are prepared in the “online
mode,” in which meteorological fields are simulated bynudging horizontal winds toward those of the reanalysis
data. In this study, we used the same model settings andprescribed meteorological data as those of Niwa et al.
(2017a). The forward and adjoint transport simulations
were performed with the horizontal grid level “glevel-5”
of NICAM, whose mean grid interval was 223 km. The
model contains 40 layers in vertical up to a height ofapproximately 45 km above the sea level.
The model operator M(.) includes spatiotemporalinterpolation from the model grids and temporal steps
to observational locations and times and the conversionof flux data from latitude–longitude grids to NICAM’s
icosahedral grid, as well as atmospheric transport. Therefore, the matrix M and its transpose can be expressed as:

M = HM′ G,(3)

MT = GT M′T HT ,(4)studies (Niwa et al. 2011a, 2012, 2017a, b, 2021) adopted
simple area averaging, in which data from within a certain distance in the surrounding latitude–longitude grids
were averaged. Subsequently, flux values were slightly
modified using scaling to ensure mass conservation.
In this study, we have updated to a more sophisticatedmethod of grid conversion, which no longer requires
mass-conserving scaling, thus retaining its linearity.
The new method initially divides the latitude–longitude
grid data into a finer-scale latitude–longitude grid, as
depicted in Fig. 1 (the bold dark-gray mesh becomes thethin light-gray mesh), i.e., rasterizing. Next, each finerscale grid is allocated to a hexagonal grid whose borders
are trigonometrically defined in latitude–longitude coordinates. In Fig. 1, the colored grid areas are allocated tothe centered icosahedral grid, which derives contributions from the surrounding latitude–longitude grids.
Consequently, the flux value converted into the ith icosahedral grid f′i can be expressed as:

fi′ =

Ni

j=1

 
ai,j f Gj /Ai ,where Gj denotes the original latitude–longitude grid thatoverlaps the ith icosahedral grid and Ni is the number of
those grids (Ni = 9 in the case of Fig. 1). The partial areaof each original latitude–longitude grid overlapped by
the ith icosahedral grid is derived from the finer latitude–
longitude grids, which is represented as ai,j (individual

where H and G represent the spatiotemporal interpolation and the flux data conversion, respectively. Matrix M′
is the integral calculation of atmospheric transport. This
study newly introduces G, which allows the flux data contained in δx to be arrayed in latitude–longitude grids, so
the optimization can be performed in that space, but notin the model grid space (i.e., the icosahedral grids). The
detailed conversion scheme is described in the next section. In NISMON-CO2, the optimization is performed
using a quasi-Newton method-based scheme termed thePreconditioned Optimizing Utility for Large-dimensional
analyses (POpULar: Fujii and Kamachi 2003; Fujii 2005).
2.2 Grid conversion scheme

The icosahedral NICAM grid shapes hexagons or pentagons, which requires an elaborate conversion fromlatitude–longitude grids, because flux input datasets are
originally prepared in latitude–longitude grids. Previous(5)

Fig. 1 Illustration of flux conversion from latitude–longitude grids
to icosahedral grids. The black line indicates the boundaries of the
icosahedral grids, while the dark-gray and light-gray lines represent
original and finer latitude–longitude grids, respectively. In theconversion, fluxes in colored areas are allocated to the centeredicosahedral gridNiwa et al. Progress in Earth and Planetary Science(2022) 9:42

colored areas in Fig. 1). The area of the ith icosahedral
grid is denoted as Ai. Because every finer latitude–longitude grid is allocated to one of the surrounding icosahedral grids, the total amount of global data is conserved.
Because of the linearity of Eq. (5) and the absence of any
mass-conserving scaling, this grid conversion scheme isentirely linear and thus represented by the matrix G in
Eq. (3). This property enables us to easily develop theadjoint code represented by GT in Eq. (4). By introducing
G and GT into the model operator and its adjoint, respectively, the control variables can be arrayed in latitude–
longitude grids.
Optimizing fluxes in the latitude–longitude grid has
two technical advantages. First, conversion of flux datafrom icosahedral grids to latitude–longitude grids after
inversion for comparison with other datasets is not
required. Such conversion often causes deterioration ofmass conservation and requires additional modification.
Second, the flux data can be provided at a higher resolution than the icosahedral grid data. In this study, we optimized the flux data at a 1° × 1° resolution, which is finer
than the resolution at the icosahedral grid (approximately
223 km). This allows the optimized fluxes to retain thehigh-resolution information originally contained in theprior flux data. During the grid conversion in this study,
we divided 1° × 1° latitude–longitude grids by 10 × 10 toobtain finer grids.
2.3 Flux model

The surface ­CO2 fluxes input to NICAM-TM can be
described as:Page 4 of 19

55-year Reanalysis (JRA-55: Kobayashi et al. 2015;
Harada et al. 2016). This is similar to the method of Olsen
and Randerson (2004). The nature of the terrestrial biosphere causes βGPP and βRE to show different diurnal variations (e.g., βGPP is zero at night because photosynthesis
is inactive, which is not the case for βRE).
In this study, fossil fuel emissions were not optimized,
i.e., they were fixed as ffos, whereas the flux deviationsof fGPP, fRE, and focn were optimized. Furthermore,
the LUC and biomass burning fluxes vary with scaling factors, whose deviations from 1, �αLUC and �αfire
were optimized. Therefore, the LUC and biomass burning fluxes were not allowed to change where their basic
fluxes were set to zero in the optimization.
A monthly temporal resolution was selected for theflux parameters f and �α, except for the ocean flux
focn , whose temporal resolution is set annually. This isbecause ocean flux variations are smaller than terrestrial
flux variations, and some erroneous values could arisein the optimized ocean fluxes due to “leakage” of terrestrial signals to the oceans. An annual temporal resolution
may help mitigate such flux leakage; however, this does
not allow the inversion to optimize seasonal variations inocean fluxes.
In NISMON-CO2 ver. 2021.1, where real observations
were used, we used fossil fuel emission data from the
GCP-Gridded Fossil Emission Dataset (GCP-GridFED:
Jones et al. 2021) ver. 2021.2 for ffos, and the GPP, RE, and
LUC data were obtained from the Vegetation Integrative
SImulator for Trace gases (VISIT: Ito and Inatomi 2012;
Ito 2019, 2021) for fGPP, fRE, and fLUC, respectively. The
fCO2 (x, t) = ffos (x, t) − βGPP fGPP (x, t) + �fGPP (x, t) + βRE fRE (x, t) + �fRE (x, t)+ (1 + �αLUC (x, t))fLUC (x, t) + (1 + �αfire (x, t))ffire (x, t)+ focn (x, t) + �focn (x, t),where x and t represent the flux location and time,
respectively. Fluxes from fossil fuel use and cement production, gross primary production (GPP) and respiration(RE) of the terrestrial biosphere, LUC, biomass burning,and oceans are denoted as ffos, fGPP, fRE, fLUC, ffire, and focn,
respectively; they are prescribed using flux datasets andhere have a monthly temporal resolution.
In Eq. (6), we separate the terrestrial biosphere fluxes
into GPP and RE components following the method of
Lokupitiya et al. (2008) because GPP and RE vary differently. The coefficients, βGPP and βRE, are scaling factorsthat produce diurnal variations. βGPP and βRE distribute
fluxes at a three-hourly resolution from the monthly
fGPP and fRE, respectively, whose values are derived fromthe downward shortwave radiation at the earth’s surface
and 2-m height air temperatures data from the Japanese(6)biomass burning emission data were obtained from the
Global Fire Emissions Database (GFED) ver. 4.1s (van derWerf et al. 2017) for ffire, and the air–sea ­CO2 exchange
data were taken from the Japan Meteorological Agency(JMA) (Iida et al. 2015, 2021) for focn. Except for focn,these data were originally provided at a higher resolution
than 1° × 1°; however, for ease of analysis, we uniformly
aggregated the data at 1° × 1°. Some of these flux datasets
were replaced in the pseudo-observation experiment inthis study, which is described in the next section.
Note that although the terrestrial flux is separated
into GPP, RE, LUC, and biomass burning, they are notintended to be fully optimized independently. As demonstrated by Niwa et al. (2021), inversion can identify
drastic flux signals, such as a large-scale biomass burning
event; however, estimating fluxes independently for eachNiwa et al. Progress in Earth and Planetary Science(2022) 9:42

Page 5 of 19

component is difficult because the inversion only uses
observations of atmospheric ­CO2 mole fractions. The
separation of those fluxes is intended to represent the different diurnal variations in GPP and RE and to consider
future possible constraints such as those from carbonylsulfide (Campbell et al. 2015).
2.4 Prior error covariance

A 4D-Var calculation requires not only a prior estimate
but also its prior error covariance matrix, as denotedby B in Eqs. (1) and (2). In the inverse experiment by
Niwa et al. (2017b), the prior error covariances (i.e., theoff-diagonal elements of B) were simply designed witha Gaussian function commonly used for the globe, inwhich a correlation scale length was arbitrarily determined. However, that assumption is too simple because
the error covariances should dynamically vary in space
and time according to the flux mechanisms.
In this study, we define the error covariance matrix
using an ensemble set of prior fluxes as follows. First, we
prepare a monthly long-term flux dataset and consider
each year of data as a member of the ensemble. Therefore, the number of years equals the number of ensemble
members, and each member contains 12 temporal components. For a prior flux component n, we approximate
the covariance matrix using m members asmB̃n, cyc =1  k  k T
δxn δxn ,
m−1
k=1

(7)

where the subscript “cyc” represents cyclostationary,
whose matrix is applied repeatedly for each year. The
superscript k is an ensemble index, and δxnk describes thedeviation of the flux of the kth member from the ensemble mean (i.e., the long-term mean). However, due to thelimited number of ensemble members (years), erroneous
correlations occur in remote areas, i.e., so-called sampling errors. Therefore, to localize covariances, we apply
the Schur product of B̃n, cyc with a correlation matrix C as

Bn,cyc = γn B̃n,cyc ◦ C.

(8)

Here, we use the Gaussian function for each element of
C as
 2
lij
Cij = exp − 2 ,(9)Lwhere the indices i and j represent different flux locations, whose distance is denoted by lij. L is a constantand set to 2000 km, which is long enough to only damp
erroneous correlations in remote areas. In practice, we let
Cij = 0 where lij > 4000 km to minimize the computationalburden. In Eq. (8), the scaling factor γn is introduced asa tuning parameter to adjust the degree to which fluxes
are subjected to change by observational constraints.
Furthermore, we replaced the resulting negative covariances with zeros, as we consider negative correlations tobe primarily generated by sampling errors. The techniquedescribed in Eq. (8) has been used as “the covariancelocalization” in an ensemble Kalman filter assimilation(Lorenc 2003); however, we have applied it here for the
4D-Var inverse analysis.
For the prior covariances of the terrestrial biosphere
GPP and RE, we used a long-term VISIT simulation for
1901–2020, generating 120 members for the covariancecalculation (Eq. 7). After applying the correlation matrix,
the resulting covariance was inflated by γn = 4 so that the
globally integrated variance was approximately (2 Pg C
­yr−1)2. The defined uncertainty of 2 Pg C ­yr−1 is derived
from the difference between the global total net flux ofthe prior data and that estimated from the growth rate of
the atmospheric observations; this is under the assumption that terrestrial flux is the dominant driver of global
flux variations. To determine the error covariance of theprior ocean flux, we used the JMA ocean flux data for
1990–2019. Because we optimized the annual meansof ocean fluxes, each annual mean was considered one
ensemble member (with 30 members in total). The covariance for the ocean fluxes was derived similarly to thatof the terrestrial biosphere using Eqs. (8) and (9) with
γn = 2. The resulting global variance was approximately(0.17 Pg C ­yr−1)2; this number is rather arbitrary compared with that for the terrestrial flux, but it was derived
after several trial experiments. When γn was increased
from 2, we obtained unrealistic variations in posterior
ocean fluxes. Although 0.17 Pg C ­yr−1 is too small as a
global ocean flux uncertainty, such small prior errors areneeded to stabilize ocean flux estimates.
For the other flux components such as LUC and biomass burning, we do not consider error correlations anddefine each prior flux error covariance matrix Bn,cyc as adiagonal matrix, whose diagonal elements (i.e., variance)
are all set to 4.0; this yields errors of 200%. We combined
all Bn,cyc and applied them to every year, whose values arearrayed in the same grid and resolution as those of theprior flux, i.e., 1° × 1°. Finally, we obtained the full matrixof B by including the variance of the global initial mole
fraction offset of (0.5 ppm)2.
Examples of the prior flux error covarianceinare given 
1/2 1/2
Fig. 2, which shows error correlations = Bij / Bii Bjj
of a given prior GPP for April and July. The original error
correlations were significant in both months, even at alarge distance. That feature was more pronounced in July,whereas larger correlations were present at a short distance in April (gray crosses in Fig. 2a and b). After performing the calculation in Eq. (8), those correlations

Niwa et al. Progress in Earth and Planetary Science(2022) 9:42

Page 6 of 19

Fig. 2 Error correlations of a prior gross primary production (GPP) flux at near 60.5° N, 100.5° E along with distance for April (a) and July (b). Grayand red crosses represent correlations before and after localization (Eq. 8), respectively. Spatial structures of the localized error correlations of a andb are shown in c and d, respectively. Furthermore, e and f are the same as c and d but for a different prior GPP flux at near 30.5° N, 110.5° E

became localized, whereas the above-mentioned feature
was maintained (red crosses in Fig. 2a and b). The error
correlations generated using this method have anisotropic spatial structures, as shown in Fig. 2c and d. Furthermore, they differ by location as shown in Fig. 2e andf. In fact, these anisotropic correlation features arise fromyear-to-year variations in the long-term VISIT data,
which are driven by changes in meteorological conditions(shortwave radiation, temperature, moisture, etc.). We
expect them to reflect actual biogeochemical mechanisms better than the previously used simple isotropic
correlations.
2.5 Observations

The observational data input to NISMON-CO2 wasderived from the dataset named Observation Package(ObsPack)-GLOBALVIEWplus and the related nearreal-time version ObsPack-NRT, provided through the
National Oceanic and Atmospheric Administration(Masarie et al. 2014). Additionally, we used other independently provided data. In NISMON-CO2 ver. 2021.1,
versions 6.1_2021_03-01 (Schuldt et al. 2021a) and6.1.1_2021-05-17 (Schuldt et al. 2021b) of ObsPackGLOBALVIEWplus and ObsPack-NRT were used,
respectively; all data are listed in Niwa (2020). To demonstrate the potential performance of NISMON-CO2 ver.
2021.1, we performed pseudo-observation inversions forthis study by emulating observational data at actual sites
and times. The locations of the observations are shownin Fig. 3. For this dataset, we used all available observations provided from selected institutes (Niwa 2020), of
which some sites have a constant data record over theentire period, and some do not. Therefore, several temporal observational gaps were present; we can investigate
to determine the effect of these data gaps on flux estimations using the pseudo-observation experiment. As
shown in Fig. 3, we only used surface observations (flask
samplings and in situ measurements at ground-based

Niwa et al. Progress in Earth and Planetary Science(2022) 9:42

Page 7 of 19

Fig. 3 Locations of observations used in the pseudo-observation inversion experiment. The observations are the same as those of NISMON-CO2
v2021.1 (Niwa 2020) but for those available during 2007–2018. Flask air sampling and in situ observations at ground-based stations are depicted byfilled magenta and unfilled purple triangles, respectively. Light blue circles denote shipboard observational locations. Tower observation sites aredenoted by cyan crosses

stations and tower sites, and shipboard observations).
This is because NISMON-CO2 ver. 2021.1 was intended
to aid in the inversion comparison of GCP (Friedlingstein, et al. 2022), which is designed to use only surface
observations so that aircraft data could be used as independent data for evaluating the inversions.
2.6 Observation–model mismatch errors

We defined the observation–model mismatch error
covariance matrix R as:

Rii =

1
,
(βri )2 Ni(10)

where ri denotes the standard deviation of the mole fractions around the ith observation and Ni represents thenumber of observations within a certain spatiotemporal
range of the ith observation. In this study, we apply a spatiotemporal range of one week, a 1000 km horizontal
diameter circle, and a 1 km vertical depth. Furthermore,min
we introduce a scaling factor β so that χ 2 := 2Jm,

where m is the number of observations and Jmin is theminimum of the cost function, should be < 1 (Tarantola
2005). We calculated ri from the simulated daily mole
fraction variations of ­CO2 at the site of observation i
ranging from one week prior and one week after theobservational timing; the simulation was performed
using NICAM-TM with prior fluxes. As denoted by

Eq. (10), R is a diagonal matrix, which assumes that allobservations are independent from each other. However,this is a simplistic assumption, specifically in cases whereobservations were obtained with high density. Therefore,
we inflate the variances for such areas by introducing Ni,
which reduces the weights of such dense observations inspace and time. This variance inflation prevents continuous in situ observations from imposing exceedingly
strong constraints compared with discrete flask sampling
observations. For NISMON-CO2 ver. 2021.1, we adopted
a strategy of using as many available observations as possible. Consequently, this resulted in duplicated observations because several institutes share the sameobservational platforms; Ni places a single-site constraint
on such collocated observations.
2.7 Pseudo‑observation inversion experiment

To investigate the general performance of NISMON-CO2
and effectiveness of the proposed inversion settings, we
performed an inverse analysis under ideal conditions;that is, we know the “true” state of fluxes. We prepared
two datasets, one of which was regarded as “true” andthe other one was used as a prior estimate. First, we performed a forward transport simulation of atmospheric­CO2 using the true flux dataset and extracted the simulated ­CO2 mole fractions at the above-mentioned observational points and times. Next, we performed inverse
experiments incorporating those pseudo-observations.

Niwa et al. Progress in Earth and Planetary Science(2022) 9:42

Page 8 of 19

Niwa et al. (2017b) conducted a similar experiment; however, the experiment in this study was designed morepractically, using a combination of continuous in situ
and discrete flask measurements (Fig. 3) and covering
a longer analysis period. Furthermore, by virtue of theimplemented grid conversion scheme, we optimized
fluxes in a latitude–longitude grid of 1° × 1° rather than in
the icosahedral grid.
The forward simulation used to generate the pseudoobservations (true run) employed the same NICAM-TM
settings used for the inversion, except for the horizontal
resolution. The horizontal resolution is “glevel-6,” whosemean grid interval is 112 km. This is one level higher
than that used for the inverse analysis (glevel-5, 223 km).
Therefore, the true run can simulate much smaller mole
fraction variations. This can provide insights into the
influence of errors in model representativeness, although
the mean grid interval of 112 km is not sufficiently small
to fully consider this aspect.
The true fluxes used in this study were obtained fromthe posterior fluxes of NISMON-CO2 ver. 2021.1 (Niwa
2020), which were generated from real observations;
therefore, the generated pseudo-observations of atmospheric ­CO2 mole fractions have realistic seasonal andinterannual variations, and the true fluxes are assumed
to reflect this. Thus, except for the fossil fuel emissions,which were taken from the same GCP-GridFED dataset,
we made moderate changes to the prior fluxes from thosein NISMON-CO2 ver. 2021.1 (described in Sect. 2.3). For
the terrestrial fluxes (fGPP, fRE, fLUC, and ffire), we employed
the same VISIT and GFEDv4.1s datasets as in NISMONCO2 ver. 2021.1, but they were climatologically averaged for 2006 to 2018. Meanwhile, the ocean fluxes (focn)
were based on different ocean flux data from Takahashi
et al. (2009), which provide climatological fluxes for thereference year 2000. Therefore, excluding the fossil fuelemissions, the prior fluxes do not contain interannual
variations.
The inversion calculation was performed from January 2007 to March 2018, with the first 12 months andlast three months disregarded as the spin-up and spindown times of the analysis, respectively. Therefore, theinversion fluxes were evaluated over a ten-year period.
Note that this analysis period is nearly one-third as
short as that of NISMON ver. 2021.1 (1990–2020) forthe ease of performing the multiple inverse experiments
described in the next section.
2.8 Sensitivity tests

In addition to the control experiment performed with thesettings mentioned above (referred to as CTL), we performed four inversions with different settings (Table 1).
The first (referred to as NO_ERR) was conducted using
the pseudo-observations produced by the model withthe same horizontal resolution as the inversion (i.e.,
glevel-5), ensuring that the model’s representation error
is no longer present. For the other three tests, we set each
updated method back to those used in previous studies(Niwa et al. 2017b, 2021). The second inversion (referred
to as ICO) was performed by optimizing fluxes on the
icosahedral grids. The third experiment (referred to as
NO_WTS) did not consider the observational weighting applied to the observation–model mismatch error,i.e., Ni = 1 in Eq. (10). In this study, for a fair comparison
with the CTL case, we changed β in Eq. (10) so that thecost function (Eq. 1) is nearly equivalent to that of CTLat the beginning of the iterative calculation. The fourth
inversion (referred to as ISO_ERR) employed the isotropic prior error covariance that was simply definedusing a Gaussian function with correlation scale lengthsof 500 km and 1000 km for land and ocean fluxes, respectively. In this experiment, also for a fair comparison, each
γn in Eq. (8) was tuned so that the global integrated valueof the resulting global prior error covariance was nearly
equivalent to that of CTL.

3 Results
3.1 Cost function and number of iterations

In the inverse analysis by Niwa et al. (2017b), the optimization for estimating 12 monthly fluxes required
approximately 50 iterations. The experiment assumed
weekly observation sampling at 65 sites. The optimization was expected to require over 50 iterations whenusing an expanded observation network with more than

Table 1 List of the sensitivity tests
Experiment nameCTL
NO_ERR

Optimization grid
Latitude–longitude (1° × 1°)

Latitude–longitude (1° × 1°)

Observation generator
resolution

Observational
weightingPrior error covarianceglevel-6

Yes

Anisotropic

glevel-5

Yes

Anisotropic

ICO

Icosahedral (~ 223 km)

glevel-6

Yes

Anisotropic

NO_WTS

Latitude–longitude (1° × 1°)

glevel-6No

Anisotropic

glevel-6

Yes

Gaussian

ISO_ERR

Latitude–longitude (1° × 1°)

Niwa et al. Progress in Earth and Planetary Science(2022) 9:42

a hundred sites consisting of continuous measurements
and flask air sampling (Fig. 3). In fact, it required a few
hundred iterations. Moreover, Fig. 4 depicts the changes
in cost function from the beginning to the 300th iteration for the CTL experiment and the other four sensitivity cases. As shown in Fig. 4, all cost functions decreasesmoothly but do not fully converge at the 300th iteration. However, we did confirm that seasonal and interannual variations in the optimized fluxes, which will be
discussed later, mostly converged at the 300th iteration.
Although they did not converge perfectly, 300 iterations
seemed to be the limit in practical use due to the computational demand. Hereafter, we present results fromfluxes optimized by 300 iterations for all experiments.
3.2 Comparison of ­CO2 mole fractions
with the pseudo‑observations

The comparison of C
­ O2 mole fractions with the pseudoobservations at Minamitorishima (24.3° N, 154.0° E) andSyowa Station (69.0° S, 39.6° E) is depicted in Fig. 5, which
are well-established background stations in the northernand southern hemispheres, respectively. At both sites,
the model with the posterior flux of CTL reproduced
atmospheric ­CO2 mole fractions without bias, indicatingthat the inversion reasonably incorporated the pseudoobservations. Furthermore, the root-mean-square (RMS)
differences denoted in the lower panels of Fig. 5 indicate
that, except for NO_ERR, the compatibility with thepseudo-observations is similar across all sensitivity tests.
The significantly smaller RMS differences in NO_ERR are
attributed to the absence of the transport model error.
The inset of Fig. 5a shows the ­CO2 mole fraction timeseries zoomed for January 2017, including the lowerresolution pseudo-observation used for NO_ERR. The

Page 9 of 19

CTL model accurately reproduced elevated events. However, their magnitudes are underestimated even after the
inversion and are rather comparable to those of the lowresolution pseudo-observation. This indicates that theinversion did not fix the surface fluxes and consequentatmospheric mole fractions beyond the model’s representation error (i.e., not overfitting to the observations).
3.3 CO2 flux maps and comparison with icosahedral grid
optimization

An example of the optimized fluxes is depicted in Fig. 6,
which shows the distributions of C
­ O2 fluxes averaged
for July 2011. (Fossil fuel emissions are excluded.) Here,we compare the optimized flux of CTL with that of ICO,
as well as the true (i.e., NISMON-CO2 ver. 2021.1) andprior fluxes. Both the CTL and ICO cases successfully
retrieved the true flux pattern from the prior fluxes, asevidenced by the continental sinks in Eurasia and theeast side of North America, sources in central Africa andSouth Asia, and the southwest side of North America.
However, the south of the Amazon shows notable discrepancies from the true fluxes. This is probably becauseof the lack of inland observations and convective verticalmixing, leading to poorly captured flux signals by surface
stations.
Although the ocean fluxes were optimized annually,
the flux pattern was retrieved to some extent. The tongue
shape of the sources in the eastern equatorial Pacific,
which persisted throughout the entire year, was suitably retrieved, and the strength of the sources increased
toward the true fluxes.
Notably, CTL, where fluxes were optimized in the latitude–longitude grid, showed a similar global distribution
to the ICO case while resolving smaller flux distributions.
For instance, over the Himalayas, CTL shows a flux pattern with finer-scale variations (not necessarily closer tothe true), reflecting the complexity of the topography.
These results indicate that the newly introduced latitude–longitude grid optimization in this study functions
reasonably well.
3.4 Seasonal and interannual variations in regional C­ O2
fluxesFig. 4 Cost function changes with iterations for the CTL case (blue)
and the other four sensitivity tests of NO_ERR (light blue), ICO (lightgreen), NO_WTS (orange), and ISO_ERR (pink)

Here, we present seasonal and interannual variations in­CO2 fluxes via various iterations to demonstrate how
those temporal variations converge as the iterative calculation progresses. The seasonal variations in nonfossil
fuel ­CO2 fluxes in the 11 TransCom terrestrial regions(Gurney et al. 2002; Baker et al. 2006) averaged for 2008–2017 are illustrated in Fig. 7. In most areas, the optimized
seasonal variation is highly consistent with that of thetrue flux. Furthermore, the optimization converges wellwithin 50 iterations for northern high-latitude areas,

Niwa et al. Progress in Earth and Planetary Science(2022) 9:42

Page 10 of 19

Fig. 5 Upper panels indicate the time series of ­CO2 mole fractions at Minamitorishima (24.3° N, 154.0° E) (a) and Syowa Station (69.0° S, 39.6° E),
Antarctica (b). Black and red lines represent the pseudo-observation and the model with the CTL posterior flux, respectively. The inset of panel(a) depicts the ­CO2 mole fraction time series zoomed for January 2017, including the lower-resolution pseudo-observation (gray line). (Thex-axis covers the whole period of January 2017, and the y-axis ranges from 406 to 420 ppm.) Lower panels indicate the differences between the
CTL model and the pseudo-observation (model–observation) at Minamitorishima (c) and Syowa Station (d), where the number indicates the
root-mean-square difference for each sensitivity experiment. Note that the NO_ERR case is compared with the pseudo-observations generated
from the lower-resolution model

which contain dense observations (Europe, Boreal Asia,
and Boreal North America). The flux variations converged within 100 iterations, even in other regions.
However, for interannual variations, convergencerequired more iterations. For example, in Fig. 8, interannual variations in regional C­ O2 fluxes had not converged
after 100 iterations, even in northern high-latitude areas.
Moreover, interannual variations converged much more
slowly than seasonal variations toward the 300the iteration. Also, significant differences still exist between 200
and 300 iterations (e.g., Boreal N. America), indicatingthat 300 iterations may not be fully sufficient for convergence. This is consistent with the decreasing cost function even at the 300th iteration (Fig. 4).
Nevertheless, after 300 iterations, the posterior flux
revealed interannual variations consistent with those of
the true flux, specifically in northern high-latitude areas.
In other areas, notable peaks and troughs of the interannual variations were more or less well retrieved, albeit

with some noticeable errors (e.g., consistently smaller
than the true flux in South America).
The slower convergence speed of the interannual variations than that of the averaged seasonal variations is
attributed to their small degree of variation. Fluxes innorthern high-latitude areas have a seasonal amplitude ofapproximately 10 Pg C y­ r−1, whereas their annual averages vary by approximately 1 Pg C y­r−1. These results
imply that optimizing such small flux variations requiresa large number of iterations. This also applies to ocean
fluxes. The yearly variations in ocean fluxes for the 11
TransCom ocean regions are depicted in Fig. 9. Although
the posterior ocean fluxes changed gradually, they
approached the true fluxes from the prior data. In fact,
the ocean flux variations were much smaller than thoseof terrestrial flux variations; the range of interannual
variations was at most 0.4 Pg ­yr−1. The inversion mostly
reproduced those interannual variations and 10-year
trends (generally decreasing trends), although some areas

Niwa et al. Progress in Earth and Planetary Science(2022) 9:42

Page 11 of 19

Fig. 6 Distributions of ­CO2 fluxes averaged for July 2011 from the true (a) and prior (b) flux datasets and the posterior fluxes of the CTL (c) and ICO(d) cases. White color areas depict fluxes close to zero. Fossil fuel emissions, which were not optimized, were excluded

had persistent biases (approximately 0.1–0.2 Pg ­Cyr−1).
However, the posterior fluxes do not fully represent the
recent uptake increases in the North Pacific and North
Atlantic (Fig. 9b and d), despite those regions having relatively dense observation networks of marine stations andships (Fig. 3). This is probably because these areas arelocated on the downwind sides of the continents. Terrestrial flux signals from the continents are transported over
the oceans with higher magnitudes than those of ocean
fluxes, which might make it difficult to clearly distinguishocean fluxes in the inversion.
Figures 7, 8 and 9 depict seasonal and interannual
variations that were derived under the assumption of nomodel error (NO_ERR, blue dotted line). The model error
can be grasped from the inset of Fig. 5a and the RMS differences in Fig. 5c and d. As shown in Fig. 7, the posterior
flux of NO_ERR showed approximately the same patterns in terrestrial seasonal variations as CTL. Although
the NO_ERR case was better than the CTL case for theinterannual variations of the terrestrial fluxes, it did not
perfectly follow the true flux (Fig. 8). Conversely, NO_
ERR showed little improvement in interannual variations

over the oceans (Fig. 9). Figures 10 and 11 confirm those
NO_ERR features, in which the RMS errors of the yearly
posterior flux variations for NO_ERR are better over terrestrial areas and similar to those of CTL over the oceans(see blue and cyan bars).
3.5 Global land–ocean partitioning

Table 2 summarizes the global, land, ocean, and total
carbon budgets, and their interannual variations are
depicted in Fig. 12. For 2008 to 2017, the global mean
true sinks over the land and the ocean were 2.62 and 2.04
Pg C ­yr−1, respectively (estimated using NISMON-CO2
ver. 2021.1 with real observations), and the prior estimates were 4.61 and 1.41 Pg C y­ r−1, respectively. Meanwhile, the inversion estimates by CTL were 3.39 and 1.29
Pg C ­yr−1, respectively; thus, the estimated global land–
ocean partitioning had a bias of approximately 0.8 Pg C
­yr−1. However, the global total sink was nearly identicalto that of the true flux, indicating that the mass conservation was thus ensured. In fact, the bias was coherently
present over the entire analysis period (Fig. 12). Nevertheless, Fig. 12 depicts that the inversion successfully

Niwa et al. Progress in Earth and Planetary Science(2022) 9:42

Page 12 of 19

Fig. 7 Monthly variations of nonfossil fuel ­CO2 fluxes averaged for 2008–2017 for the 11 TransCom terrestrial regions (Gurney et al. 2002; Baker et al.
2006). Light-gray and dark-gray solid lines represent the prior and true fluxes, respectively. The other colors represent posterior fluxes from different
iterations (10–300) of CTL. The blue dotted line denotes the 300th iteration result in the case where model errors do not exist (NO_ERR)

reproduced the interannual variation patterns from theconstant prior estimates, even for the ocean, which moderately changes. Figure 12 also indicates that 100–150
iterations are required for convergence of the land–ocean
partitioning, which is less than those required for estimates of regional interannual variations (Figs. 8 and 9).
Compared with the CTL case, the NO_ERR case did notnecessarily provide better estimates of land–ocean partitioning and displayed the same coherent bias and similar interannual variations (Table 2; Fig. 12). Furthermore,
Fig. 12 shows that the differences between the posterior
and true fluxes of the land and the ocean are gradually
increasing over time. From the former five-year period(2008–2012) to the latter five-year period (2013–2017),
the mean differences between the posterior and truefluxes of the global land and ocean changed from 0.66 to
0.89 Pg C y­ r−1 and from 0.69 to 0.82 Pg C y­ r−1, respectively. This is attributable to the fact that prior ocean

fluxes do not have interannual variations but are accompanied by small prior uncertainties, which makes itslightly difficult to keep up with the increasing trend of
the ocean sink and, consequently, the decreasing trend of
the land sink.
3.6 Sensitivity to inversion settings

As shown in Table 2, NO_ERR (2.74 × ­10−1 gC­m−2 ­day−1) had the smallest RMS error for the overall
fluxes during 2008–2017, which were calculated against
the true fluxes, while the second smallest value wasobtained in CTL (3.49 × ­10−1 gC m
­ −2 ­day−1). The smaller
RMS error of NO_ERR compared with that of CTL is
due to fluxes over land. This could be attributed to more
inhomogeneous and larger land fluxes than the ocean,and mole fractions simulated from terrestrial fluxes are
more sensitive to the model transport resolution. The
difference between CTL and NO_ERR is the modelNiwa et al. Progress in Earth and Planetary Science(2022) 9:42

Page 13 of 19

Fig. 8 Same as Fig. 7, but for yearly variations from 2008 to 2017

reproducibility of atmospheric mole fractions, which isderived from the horizontal resolution of the model that
produced the pseudo-observations. Therefore, the result
indicates that a change in the model horizontal resolution could highly impact the terrestrial flux estimates
than the ocean flux estimates.
The other sensitivity tests, ICO, NO_WTS, and ISO_
ERR, produced similar values for global land and ocean
uptakes, i.e., similar biases; however, their RMS errors forthe overall fluxes were larger than those of CTL (Table 2).
In the ICO case, the larger RMS error came from the terrestrial fluxes, a feature that can also be found in interannual flux variations (Fig. 10). This result implies thatthe latitude–longitude optimization introduced in thisstudy not only provided higher resolution fluxes but also
improved the inversion estimates. The NO_WTS caseshowed larger RMS errors over the land and ocean, which
is true for the overall fluxes (Table 2) and regional interannual variations (Figs. 10 and 11). This impact is highest
in the sensitivity tests, indicating that the spatiotemporalweighting of the observation–model mismatch errors
plays a crucial role. The ISO_ERR case also had a larger
RMS error than the CTL case for the overall terrestrial
fluxes (Table 2). However, this was not necessarily thecase for the regional interannual variations; ISO_ERR
had smaller errors in some regions than CTL (Fig. 10).

4 Discussion
The temporal variations in ­CO2 fluxes were retrieved
well in this study, demonstrating the reliability of theflux estimates by NISMON-CO2. It is also promising thata weakly decreasing trend was retrieved in some ocean
regions, because detecting such small flux changes is critical for monitoring carbon sink capacities of natural reservoirs and gaining insights into the future global carbon
cycle considering the global warming conditions.
However, areas in the southern hemisphere and tropics still have substantial errors in flux estimates, possibly
due to a sparse observational network there. The inversion reproducing atmospheric mole fractions at the two

Niwa et al. Progress in Earth and Planetary Science(2022) 9:42

Page 14 of 19

Fig. 9 Same as Fig. 8, but for the 11 TransCom ocean regions

Fig. 10 Root-mean-square errors of the yearly fluxes against the truefluxes for the TransCom terrestrial regions. Results of the CTL and theother four sensitivity tests are illustrated here

sites without any bias is depicted in Fig. 5; however, thosesites are in regions with persistent flux biases (i.e., North
Pacific and Southern Ocean, Fig. 9). This suggests that

Fig. 11 Same as Fig. 10, but for the TransCom ocean regions

fluxes in other areas compensate to make the modeledmole fractions consistent with the observations. Suchflux compensations could occur where observations aresparse, which may lead to further errors in flux estimates.

Niwa et al. Progress in Earth and Planetary Science(2022) 9:42

Page 15 of 19

Table 2 Global-scale net sinks and errors for 2008–2017. The sink values do not include fossil fuel emissions. These errors are derived
using root-mean-square errors against the true fluxes over the land, the ocean, and the global total carbon fluxes
Flux data

Sinks (Pg C ­yr−1)
Land

Root-mean-square errors (× ­10−1 g C ­m−2 ­day−1)
Ocean

Total

Land

Ocean

Total

True

2.62

2.04

4.67

n/a

n/a

n/a

Prior

4.61

1.41

6.02

5.92

0.33

5.93

CTL

3.39

1.29

4.68

3.48

0.29

3.49

NO_ERR

3.38

1.29

4.67

2.73

0.29

2.74

ICO

3.38

1.31

4.69

4.19

0.26

4.19

NO_WTS

3.39

1.26

4.65

4.61

0.32

4.62

ISO_ERR

3.38

1.31

4.69

4.52

0.29

4.53

Fig. 12 Yearly variations of the global land (a), ocean (b), and total (c) ­CO2 fluxes. (Fossil fuel emissions are excluded.) Light-gray and dark-gray solidlines represent the prior and true fluxes, respectively. The other colors indicate posterior fluxes from different iterations (10–300) of CTL. The blue
dotted line denotes the 300th iteration result where model errors are nonexistent (NO_ERR)

This is a typical inversion problem for long-lived species,
such as ­CO2.
The persistent bias and long-term trend error in the
global carbon budgets indicates difficulty in separating
the land and ocean budgets with the current observation
network. This land–ocean partitioning error could be
mitigated by using the best ocean flux dataset available

for the prior estimate; however, assessing which dataset
is the best is difficult. Therefore, for the prior ocean fluxin the real-observation inversion of NISMON-CO2 ver.
2021.1, we used the JMA’s yearly varying ocean flux data
(Iida et al. 2021), which is likely representing the actual
ocean sink trend better than the climatological data
(Takahashi et al. 2009) used in the pseudo-observation

Niwa et al. Progress in Earth and Planetary Science(2022) 9:42

experiment. Performing multiple inversions with different ocean flux datasets could also be useful to account for
ocean flux uncertainties. Furthermore, using additional­O2/N2 data, as Rödenbeck et al. (2008) did, might beanother approach to improve the land–ocean partitioning. However, they are left for future studies.
It should be noted that we did not select observationalsites based on their temporal coverage. Therefore, wellreproduced interannual variations, as shown in Figs. 8
and 9, indicate that such temporal gaps in the observational data do not significantly affect the estimation ofinterannual variations. Nonetheless, some interannual
variations remain poorly estimated and might be affectedby the observational gaps. Such observational gap effects
could be mitigated by selecting only observations with
sufficient temporal coverage, as some other inversion
studies do (e.g., Chevallier et al. 2010). However, reducing observations would increase overall flux errors. Furthermore, due to changes in meteorological conditions,misestimation of interannual variations could occur evenwith temporally sufficient observations. When winddirection at a site changes from year to year, its observational constraints may also change. Therefore, we
adopted a strategy to use as many observations as possible without considering any site-specific selection interms of temporal data coverage. This research demonstrates that such a strategy can still reproduce reasonableinterannual variations. It should also be noted that thisinversion only used surface observations. However, various aircraft data are available, which could provide effective constraints in flux estimates (Niwa et al. 2012, 2021).
Furthermore, satellite observations (Crisp et al. 2008;
Yokota et al. 2009) may provide different insights fromsurface and aircraft observations, although assessing
biases and errors of satellite observations can be difficult.
Incorporating such observations is one way to improveNISMON-CO2.
In contrast to a previous study by Niwa et al. (2017b),
we used a higher resolution model to prepare pseudoobservations (~ 112 km). This provided insight into therepresentation error of the model, although the actual
degree of the model’s representation error is probably
much greater. The inversion performance in retrieving seasonal and interannual variations in surface fluxes
depicted in this study is possibly overestimated due to thesmall model errors. Nevertheless, the difference between
CTL and NO_ERR provides valuable insight into theeffect of model resolution, from which we can predicta degree of improvement in a future inversion with ahigher resolution model.
One major update to the inversion settings of Niwa
et al. (2017b) is the implementation of a grid conversion in the inverse calculation, which enables a

Page 16 of 19

high-resolution latitude–longitude grid optimization.
The computational time required for the grid conversion
is mere minutes, and the convergence speed is almost thesame as that of ICO (Fig. 4); therefore, this implementation does not induce any substantial increase in computational costs. Because it uses the same transport model
resolution as before, the amount of information that canbe extracted from the observations does not necessarily
increase. Nevertheless, it is valuable because the original
information in the prior flux datasets does not deteriorate. Furthermore, such highly resolved fluxes could be
more useful for regional- or national-scale budget analyses, demand for which has recently increased for use inmonitoring and reducing national carbon emissions. This
scheme can be widely utilized because it is also applicable to inversion systems that use a conventional latitude–
longitude grid model. The grid conversion scheme can beapplied to convert from higher- to lower-resolution latitude–longitude grids.
Other updates to the inversion settings include theobservational weighting and anisotropic prior error
covariance. Although both of these updates improved theinversion results, the impact of observational weighting
was considerably greater. Weighting is performed using
Eq. (10) with three parameters (horizontal lengths, vertical lengths, and temporal range for counting nearby
observations), which is much simpler than changing
the model resolution or the prior flux error covariancedesign. In this study, these parameters were determined
after a few preliminary sensitivity tests but were not thoroughly investigated. Therefore, equation modification orparameter sweeping could efficiently improve the inversion. The impact of the anisotropic prior error covarianceis comparable to that of observational weighting on the
RMS error of the overall fluxes (Table 2), whereas theRMS errors of the regional and temporal variations are
marginal (Figs. 10 and 11). In fact, the improved result
from CTL compared to ISO_ERR is partially attributed
to the fact that the true flux was derived from NISMONCO2 ver. 2021.1, which used the same inversion settings,
including the prior error covariance, as those of CTL.
This might have given CTL some advantages. Nonetheless, the difference between CTL and ISO_ERR indicatespotential improvements made by the anisotropic error
covariance and inversion estimates could be improved byshaping the error covariance according to real flux variations. Because a long-term flux dataset is used to derive
not only the prior estimate but also its error covariance,
preparing such a flux dataset has become more vital. Furthermore, the anisotropic error covariance is also valuable for providing high-resolution fluxes with the gridconversion scheme. If isotropic error covariance is used,
the patterns of flux increments from prior fluxes would

Niwa et al. Progress in Earth and Planetary Science(2022) 9:42

be blurred even when the optimization is performed in
the latitude–longitude grid as the isotropic error covariance works on spatial smoothing. However, as depicted
in Fig. 2, the anisotropic error covariance retains the finescale patterns derived from flux dynamics.

5 Conclusions
In this study, we described the updated inversion settings for NISMON-CO2. Furthermore, to demonstrate
the practical application of long-term inverse analysis, we conducted 10-year inversion experiments using
pseudo-observations that emulate an inhomogeneous
observation network consisting of continuous measurements and flask air samplings with actual data gaps.
Despite a persistent global bias in the land–ocean partitioning, the seasonal and interannual variations were
sufficiently retrieved in the experiment, demonstrating
the reliability of NISMON-CO2 and its real-observation
inversion product (NISMON-CO2 ver. 2021.1). Moreover, sensitivity experiments revealed that the flux estimation improvements were attributed to grid conversion,
observational weighting, and anisotropic prior error
covariance, all of which were newly introduced to NISMON-CO2. Hence, these schemes will be useful as thedefault settings of NISMON-CO2 and can also be applied
to inversion analyses of other atmospheric constituents(e.g., methane).
Since the adoption of the Paris Agreement in 2015, everyone globally is working toward reducing greenhousegas emissions with the common goal to limit global
warming well below 2 °C, preferably below 1.5 °C, compared with pre-industrial levels. To verify the achievements in ­CO2 emission reductions, regular evaluationsof surface ­CO2 fluxes are imperative. To this end, we will
operationally conduct inverse analysis using NISMONCO2 and provide science-based estimates of C
­ O2 fluxes,
as well as developing the system further.
Abbreviations
NICAM: Nonhydrostatic Icosahedral Atmospheric Model; NISMON-CO2:
NICAM-based Inverse Simulation for Monitoring ­CO2; NICAM-TM: NICAMbased Transport Model; 4D-Var: Four-dimensional variational method; GCP:
Global Carbon Project; GPP: Gross primary production; RE: Respiration; LUC:
Land-use change; GridFED: Gridded Fossil Emission Dataset; VISIT: Vegetation
Integrative SImulator for Trace gases; GFED: Global Fire Emissions Database;
JMA: Japan Meteorological Agency; ObsPack: Observation Package.
Acknowledgements
Comments from two anonymous reviewers, which helped us improve themanuscript, were greatly appreciated. The inverse analyses and systemdevelopment were performed on the supercomputer systems at the National
Institute for Environmental Studies and Meteorological Research Institutes(NEC SX-Aurora TSUBASA and FUJITSU PRIMERGY CX2550M5). YN is gratefulto all the staff working on those supercomputer systems. YN also thanks Prof.
Masaki Satoh and other members at JAMSTEC, RIKEN, the University of Tokyo,
and NIES for developing NICAM. Furthermore, we acknowledge Ray Langenfelds, Paul Krummel, Zoe Loh, Doug Worthy, Martin Steinbacher, Juha Hatakka,

Page 17 of 19

Tuula Aalto, Luciana V. Gatti, Kazuyuki Saito, Michel Ramonet, Marc Delmotte,
Gilles Bentz, Morgan Lopez, Francois Gheusi, N. Mihalopoulos, Josep-Anton
Morgui, Olivier Laurent, Britton Stephens, Hitoshi Mukai, Toshinobu Machida,
Motoki Sasakawa, Shohei Nomura, Yukio Terao, Shin-Ichiro Nakaoka, Yasunori
Tohjima, Cathrine Lund Myhre, Ove Hermanssen, Gordon Brailsford, Sylvia
Nichol, Arlyn Andrews, Ed Dlugokencky, John Lee, Colm Sweeney, Kirk
Thoning, Pieter Tans, Stephan De Wekker, Marc L. Fischer, Dan Jaffe, Kathryn
McKain, Brian Viner, John B. Miller, Anna Karion, Charles Miller, Christopher D.
Sloop, Peter Bakwin, Ralph Keeling, Shane Clark, Bill Paplawsky, Adam Cox,
Stephen Walker, Eric Morgan, Eric Hintsa, Shuji Aoki, Shinji Morimoto, Daisuke
Goto, and Kenneth Schuldt for providing their observations through the
ObsPack framework or directly.
Author contributions
YN conceived and designed the study and conducted the inversion experiments. KI developed the grid conversion scheme. AI and YI made the important prior flux data from VISIT and the JMA ocean flux, respectively, and helped
interpret the land and ocean flux variations in the analysis. All authors readand approved the final manuscript.
Funding
This study was supported by the Environment Research and Technology
Development Fund of the Environmental Restoration and Conservation
Agency of Japan (JPMEERF21S20810). The work was also supported by the
JSPS KAKENHI, Grant Number JP19K03976 and FLAGSHIP2020, MEXT, withinpriority study 4 (advancement of meteorological and global environmental
predictions utilizing observational “Big Data”).
Availability of data and material
The following flux datasets used in this study are available at: GridFED v2021.2:
https://​doi.​org/​10.​5281/​zenodo.​55651​99. VISIT: https://​doi.​org/​10.​17595/​
20210​521.​001. JMA ocean flux: https://​www.​data.​jma.​go.​jp/​gmd/​kaiyou/​engli​
sh/​co2_​flux/​co2_​flux_​data_​en.​html. GFED v4.1s: https://​www.​globa​lfire​data.​org/. NISMON-CO2 v2021.1: https://​doi.​org/​10.​17595/​20201​127.​001.

Declarations
Competing interests
The authors declare that they have no competing interest.
Author details
1
National Institute for Environmental Studies, 16‑2 Onogawa, Tsukuba,
Ibaraki 305‑8506, Japan. 2 Meteorological Research Institute, 1‑1 Nagamine,
Tsukuba, Ibaraki 305‑0052, Japan. 3 Japan Agency for Marine-Earth Scienceand Technology, 3173‑25 Showa‑machi, Kanawaza‑ku, Yokohama, Kanagawa
236‑0001, Japan. 4 Nagoya University, Furo‑cho, Chikusa‑ku, Nagoya, Aichi
464‑8601, Japan. 5 Japan Meteorological Agency, 3‑6‑9 Toranomon, Minato
City, Tokyo 105‑8431, Japan.
Received: 19 January 2022 Accepted: 2 August 2022

References
Baker DF, Law RM, Gurney KR, Rayner P, Peylin P, Denning AS, Bousquet P, Bruhwiler L, Chen Y, Ciais P, Fung IY, Heimann M, John J, Maki T, Maksyutov S,
Masarie K, Prather M, Pak B, Taguchi S, Zhu Z (2006) TransCom 3 inversion
intercomparison: impact of transport model errors on the interannual
variability of regional ­CO2 fluxes, 1988–2003. Glob Biogeochem Cycles.
https://​doi.​org/​10.​1029/​2004G​B0024​39
Bousquet P, Peylin P, Ciais P, Quéré CL, Friedlingstein P, Tans PP (2000) Regionalchanges in carbon dioxide fluxes of land and oceans since 1980. Science
290(5495):1342–1346. https://​doi.​org/​10.​1126/​scien​ce.​290.​5495.​1342
Campbell JE, Whelan ME, Seibt U, Smith SJ, Berry JA, Hilton TW (2015) Atmospheric carbonyl sulfide sources from anthropogenic activity: implications
for carbon cycle constraints. Geophys Res Lett 42(8):3004–3010. https://​doi.​org/​10.​1002/​2015G​L0634​45
Chevallier F, Ciais P, Conway TJ, Aalto T, Anderson BE, Bousquet P, Brunke EG,
Ciattaglia L, Esaki Y, Fröhlich M, Gomez A, Gomez-Pelaez AJ, Haszpra L,

Niwa et al. Progress in Earth and Planetary Science(2022) 9:42

Krummel PB, Langenfelds RL, Leuenberger M, Machida T, Maignan F, Matsueda H, Morguí JA, Mukai H, Nakazawa T, Peylin P, Ramonet M, Rivier L,
Sawa Y, Schmidt M, Steele LP, Vay SA et al (2010) ­CO2 surface fluxes at grid
point scale estimated from a global 21 year reanalysis of atmosphericmeasurements. J Geophys Res. https://​doi.​org/​10.​1029/​2010J​D0138​87
Chevallier F (2021) Fluxes of carbon dioxide from managed ecosystems
estimated by national inventories compared to atmospheric inverse
modeling. Geophys Res Lett 48(15):2021–093565. https://​doi.​org/​10.​1029/​2021G​L0935​65
Crisp D, Miller CE, DeCola PL (2008) NASA orbiting carbon observatory: measuring the column averaged carbon dioxide mole fraction from space. J
Appl Remote Sens 2:023508. https://​doi.​org/​10.​1117/1.​28984​57
Deng Z, Ciais P, Tzompa-Sosa ZA, Saunois M, Qiu C, Tan C, Sun T, Ke P, Cui Y,
Tanaka K, Lin X, Thompson RL, Tian H, Yao Y, Huang Y, Lauerwald R, Jain
AK, Xu X, Bastos A, Sitch S, Palmer PI, Lauvaux T, d’Aspremont A, Giron
C, Benoit A, Poulter B, Chang J, Petrescu AMR, Davis SJ, Liu Z et al (2022)
Comparing national greenhouse gas budgets reported in UNFCCC inventories against atmospheric inversions. Earth Syst Sci Data 14:1639–1675.
https://​doi.​org/​10.​5194/​essd-​14-​1639-​2022
Friedlingstein P, Cox P, Betts R, Bopp L, von Bloh W, Brovkin V, Cadule P, Doney
S, Eby M, Fung I, Bala G, John J, Jones C, Joos F, Kato T, Kawamiya M, Knorr
W, Lindsay K, Matthews HD, Raddatz T, Rayner P, Reick C, Roeckner E,
Schnitzler K-G, Schnur R, Strassmann K, Weaver AJ, Yoshikawa C, Zeng N(2006) Climate–carbon cycle feedback analysis: Results from the C4MIP
model intercomparison. J Clim 19(14):3337–3353. https://​doi.​org/​10.​
1175/​JCLI3​800.1
Friedlingstein P, O’Sullivan M, Jones MW, Andrew RM, Hauck J, Olsen A,
Peters GP, Peters W, Pongratz J, Sitch S, Le Quéré C, Canadell JG, Ciais P,
Jackson RB, Alin S, Aragão LEOC, Arneth A, Arora V, Bates NR, Becker M,
Benoit-Cattin A, Bittig HC, Bopp L, Bultan S, Chandra N, Chevallier F, Chini
LP, Evans W, Florentie L, Forster PM et al (2020) Global carbon budget
2020. Earth Syst Sci Data 12(4):3269–3340. https://​doi.​org/​10.​5194/​
essd-​12-​3269-​2020
Friedlingstein P, Jones MW, O’Sullivan M, Andrew RM, Bakker DCE, Hauck J, Le
Quéré C, Peters GP, Peters W, Pongratz J, Sitch S, Canadell JG, Ciais P, Jackson RB, Alin SR, Anthoni P, Bates NR, Becker M, Bellouin N, Bopp L, Chau
TTT, Chevallier F, Chini LP, Cronin M, Currie KI, Decharme B, Djeutchouang
L, Dou X, Evans W, Feely RA et al (2022) Global carbon budget 2021. Earth
Syst Sci Data 14:1917–2005. https://​doi.​org/​10.​5194/​essd-​14-​1917-​2022
Fujii Y, Kamachi M (2003) A nonlinear preconditioned quasi-newton method
without inversion of a first-guess covariance matrix in variational analyses. Tellus A 55(5):450–454. https://​doi.​org/​10.​1034/j.​1600-​0870.​2003.​
00030.x
Fujii Y (2005) Preconditioned optimizing utility for large-dimensional
analyses (POpULar). J Oceanogr 61(1):167–181. https://​doi.​org/​10.​1007/​
s10872-​005-​0029-z
Gurney KR, Law RM, Denning AS, Rayner PJ, Baker D, Bousquet P, Bruhwiler L,
Chen Y-H, Ciais P, Fan S, Fung IY, Gloor M, Heimann M, Higuchi K, John
J, Maki T, Maksyutov S, Masarie K, Peylin P, Prather M, Pak BC, Randerson
J, Sarmiento J, Taguchi S, Takahashi T, Yuen C-W (2002) Towards robust
regional estimates of C
­ O2 sources and sinks using atmospheric transportmodels. Nature 415(6872):626–630. https://​doi.​org/​10.​1038/​41562​6aHarada Y, Kamahori H, Kobayashi C, Endo H, Kobayashi S, Ota Y, Onoda H,
Onogi K, Miyaoka K, Takahashi K (2016) The JRA-55 reanalysis: representation of atmospheric circulation and climate variability. J Meteor Soc Jpn
94(3):269–302. https://​doi.​org/​10.​2151/​jmsj.​2016-​015
Iida Y, Kojima A, Takatani Y, Nakano T, Sugimoto H, Midorikawa T, Ishii M (2015)
Trends in ­pCO2 and sea–air C­ O2 flux over the global open oceans for thelast two decades. J. Oceanogr. 71(6):637–661. https://​doi.​org/​10.​1007/​
s10872-​015-​0306-4
Iida Y, Takatani Y, Kojima A, Ishii M (2021) Global trends of ocean C­ O2 sink andocean acidification: an observation-based reconstruction of surface
ocean inorganic carbon variables. J Oceanogr 77(2):323–358. https://​doi.​
org/​10.​1007/​s10872-​020-​00571-5
Ito A, Inatomi M (2012) Use of a process-based model for assessing themethane budgets of global terrestrial ecosystems and evaluationof uncertainty. Biogeosciences 9:759–773. https://​doi.​org/​10.​5194/​
bg-9-​759-​2012
Ito A (2019) Disequilibrium of terrestrial ecosystem ­CO2 budget caused by
disturbance-induced emissions and non-CO2 carbon export flows: aPage 18 of 19

global model assessment. Earth Syst Dyn 10(4):685–709. https://​doi.​
org/​10.​5194/​esd-​10-​685-​2019
Ito A (2021) Output data of greenhouse gas budget and carbon cycle simulated by the VISIT terrestrial ecosystem model, ver.2021.1. NIES, https://​doi.​org/​10.​17595/​20210​521.​001 (Reference date: 2021/05/21)
Jones MW, Andrew RM, Peters GP, Janssens-Maenhout G, De-Gol AJ, Ciais P,
Patra PK, Chevallier F, Le Quéré C (2021) Gridded fossil ­CO2 emissions
and related ­O2 combustion consistent with national inventories 1959–
2018. Sci. Data 8(1):2. https://​doi.​org/​10.​1038/​s41597-​020-​00779-6
Keeling CD, Whorf TP, Wahlen M, van der Plichtt J (1995) Interannual
extremes in the rate of rise of atmospheric carbon dioxide since 1980.
Nature 375(6533):666–670. https://​doi.​org/​10.​1038/​37566​6a0
Kobayashi S, Ota Y, Harada Y, Ebita A, Moriya M, Onoda H, Onogi K, Kamahori
H, Kobayashi C, Endo H, Miyaoka K, Takahashi K (2015) The JRA-55
reanalysis: general specifications and basic characteristics. J Meteor
Soc Jpn 93(1):5–48. https://​doi.​org/​10.​2151/​jmsj.​2015-​001
Liu J, Baskaran L, Bowman K, Schimel D, Bloom AA, Parazoo NC, Oda T, Carroll D, Menemenlis D, Joiner J, Commane R, Daube B, Gatti LV, McKain
K, Miller J, Stephens BB, Sweeney C, Wofsy S (2021) Carbon monitoring
system flux net biosphere exchange 2020 (CMS-Flux NBE 2020). Earth
Syst Sci Data 13(2):299–330. https://​doi.​org/​10.​5194/​essd-​13-​299-​2021
Lokupitiya RS, Zupanski D, Denning AS, Kawa SR, Gurney KR, Zupanski
M (2008) Estimation of global C­ O2 fluxes at regional scale using the
maximum likelihood ensemble filter. J Geophys Res. https://​doi.​org/​10.​1029/​2007J​D0096​79
Lorenc AC (2003) The potential of the ensemble Kalman filter for NWP-a
comparison with 4D-Var. Q J Roy Meteor Soc 129(595):3183–3203.
https://​doi.​org/​10.​1256/​qj.​02.​132
Masarie KA, Peters W, Jacobson AR, Tans PP (2014) ObsPack: a framework forthe preparation, delivery, and attribution of atmospheric greenhousegas measurements. Earth Syst Sci Data 6(2):375–384. https://​doi.​org/​
10.​5194/​essd-6-​375-​2014
Niwa Y (2020) Long-term global ­CO2 fluxes estimated by NICAM-based
inverse simulation for monitoring C­ O2 (NISMON-CO2), ver.2021.1.
Center for Global Environmental Research, NIES. https://​doi.​org/​10.​
17595/​20201​127.​001 (Reference date: 2021/12/24)
Niwa Y, Patra PK, Sawa Y, Machida T, Matsueda H, Belikov D, Maki T, Ikegami
M, Imasu R, Maksyutov S, Oda T, Satoh M, Takigawa M (2011a) Threedimensional variations of atmospheric C­ O2: aircraft measurements andmulti-transport model simulations. Atoms Chem Phys 11(24):13359–
13375. https://​doi.​org/​10.​5194/​acp-​11-​13359-​2011
Niwa Y, Tomita H, Satoh M, Imasu R (2011b) A three-dimensional icosahedral grid advection scheme preserving monotonicity and consistency
with continuity for atmospheric tracer transport. J Meteor Soc Japan
89(3):255–268. https://​doi.​org/​10.​2151/​jmsj.​2011-​306
Niwa Y, Machida T, Sawa Y, Matsueda H, Schuck TJ, Brenninkmeijer CAM,
Imasu R, Satoh M (2012) Imposing strong constraints on tropicalterrestrial ­CO2 fluxes using passenger aircraft based measurements. J
Geophys Res. https://​doi.​org/​10.​1029/​2012J​D0174​74
Niwa Y, Tomita H, Satoh M, Imasu R, Sawa Y, Tsuboi K, Matsueda H, Machida
T, Sasakawa M, Belan B, Saigusa N (2017a) A 4D-Var inversion system
based on the icosahedral grid model (NICAM-TM 4D-Var v1.0)—Part
1: Offline forward and adjoint transport models. Geosci Model Dev
10(3):1157–1174. https://​doi.​org/​10.​5194/​gmd-​10-​1157-​2017a
Niwa Y, Fujii Y, Sawa Y, Iida Y, Ito A, Satoh M, Imasu R, Tsuboi K, Matsueda H,
Saigusa N (2017b) A 4D-Var inversion system based on the icosahedral
grid model (NICAM-TM 4D-Var v1.0)—Part 2: Optimization scheme andidentical twin experiment of atmospheric ­CO2 inversion. Geosci Model
Dev 10(6):2201–2219. https://​doi.​org/​10.​5194/​gmd-​10-​2201-​2017b
Niwa Y, Sawa Y, Nara H, Machida T, Matsueda H, Umezawa T, Ito A, Nakaoka
S-I, Tanimoto H, Tohjima Y (2021) Estimation of fire-induced carbon
emissions from equatorial asia in 2015 using in situ aircraft and shipobservations. Atoms Chem Phys 21(12):9455–9473. https://​doi.​org/​10.​
5194/​acp-​21-​9455-​2021
Olsen SC, Randerson JT (2004) Differences between surface and column
atmospheric ­CO2 and implications for carbon cycle research. J Geophys
Res. https://​doi.​org/​10.​1029/​2003J​D0039​68
Patra PK, Ishizawa M, Maksyutov S, Nakazawa T, Inoue G (2005a) Role ofbiomass burning and climate anomalies for land-atmosphere carbon
fluxes based on inverse modeling of atmospheric C­ O2. Glob Biogeochem
Cycles. https://​doi.​org/​10.​1029/​2004G​B0022​58

Niwa et al. Progress in Earth and Planetary Science(2022) 9:42

Patra PK, Maksyutov S, Ishizawa M, Nakazawa T, Takahashi T, Ukita J (2005b)
Interannual and decadal changes in the sea-air C­ O2 flux from atmospheric ­CO2 inverse modeling. Glob Biogeochem Cycles. https://​doi.​org/​
10.​1029/​2004G​B0022​57
Peters W, Jacobson AR, Sweeney C, Andrews AE, Conway TJ, Masarie K, Miller
JB, Bruhwiler LMP, Pétron G, Hirsch AI, Worthy DEJ, van der Werf GR,
Randerson JT, Wennberg PO, Krol MC, Tans PP (2007) An atmosphericperspective on north american carbon dioxide exchange: CarbonTracker.
Proc Natl Acad Sci (USA) 104(48):18925–18930. https://​doi.​org/​10.​1073/​
pnas.​07089​86104
Peylin P, Law RM, Gurney KR, Chevallier F, Jacobson AR, Maki T, Niwa Y, Patra
PK, Peters W, Rayner PJ, Rödenbeck C, van der Laan-Luijkx IT, Zhang X(2013) Global atmospheric carbon budget: results from an ensemble ofatmospheric ­CO2 inversions. Biogeosciences 10(10):6699–6720. https://​doi.​org/​10.​5194/​bg-​10-​6699-​2013
Rödenbeck C, Quérér CL, Heimann M, Keeling RF (2008) Interannual variability
in oceanic biogeochemical processes inferred by inversion of atmospheric ­O2/N2 and ­CO2 data. Tellus B 60(5):685–705. https://​doi.​org/​10.​1111/j.​1600-​0889.​2008.​00375.x
Satoh M, Matsuno T, Tomita H, Miura H, Nasuno T, Iga S (2008) Nonhydrostatic
icosahedral atmospheric model (NICAM) for global cloud resolving
simulations. J Comput Phys 227:3486–3514. https://​doi.​org/​10.​1016/j.​jcp.​
2007.​02.​006
Satoh M, Tomita H, Yashiro H, Miura H, Kodama C, Seiki T, Noda AT, Yamada Y,
Goto D, Sawada M, Miyoshi T, Niwa Y, Hara M, Ohno T, Iga S, Arakawa T,
Inoue T, Kubokawa H (2014) The non-hydrostatic icosahedral atmospheric model: description and development. Prog Earth Planet Sci
1(1):1–32. https://​doi.​org/​10.​1186/​s40645-​014-​0018-1
Schuldt KN, Mund J, Luijkx IT, Aalto T, Abshire JB, Aikin K, Andrews A, Aoki S,
Apadula F, Baier B, Bakwin P, Bartyzel J, Bentz G, Bergamaschi P, Beyersdorf
A, Biermann T, Biraud SC, Boenisch H, Bowling D, Brailsford G, Chen G,
Chen H, Chmura L, Clark S, Climadat S, Colomb A, Commane R, ConilS, Cox A, Cristofanelli P et al (2021a) Multi-laboratory compilation of
atmospheric carbon dioxide data for the period 1957–2019; obspack_co2_1_GLOBALVIEWplus_v6.1_2021a–03–01. NOAA Global Monitoring
Laboratory, https://​doi.​org/​10.​25925/​20201​204
Schuldt KN, Jacobson AR, Aalto T, Andrews A, Bakwin P, Bergamaschi P, Biermann T, Biraud SC, Chen H, Colomb A, Conil S, Cristofanelli P, De Mazière
M, De Wekker S, Delmotte M, Dlugokencky E, Emmenegger L, Fischer
ML, Hatakka J, Heliasz M, Hermanssen O, Holst J, Jaffe D, Karion A, Kazan
V, Keronen P, Kominkova K, Kubistin D, Laurent O, Laurila T, et al (2021b)
Multi-laboratory compilation of atmospheric carbon dioxide data for theperiod 2020–2021b; obspack_co2_1_NRT_v6.1.1_2021b–05–17. NOAA
Earth System Research Laboratory, Global Monitoring Laboratory, https://​doi.​org/​10.​25925/​20210​517
Takahashi T, Sutherland SC, Wanninkhof R, Sweeney C, Feely RA, Chipman DW,
Hales B, Friederich G, Chavez F, Sabine C, Watson A, Bakker DCE, Schuster
U, Metzl N, Yoshikawa-Inoue H, Ishii M, Midorikawa T, Nojiri Y, Körtzinger A,
Steinhoff T, Hoppema M, Olafsson J, Arnarson TS, Tilbrook B, Johannessen
T, Olsen A, Bellerby R, Wong CS, Delille B, Bates NR, de Baar HJW (2009)
Climatological mean and decadal change in surface ocean ­pCO2, and net
sea-air ­CO2 flux over the global oceans. Deep-Sea Res Part II 56:554–577
Tarantola A (2005) Inverse problem theory and methods for model parameter
estimation. Society for Industrial and Applied Mathematics, Philadelphia,
p 342
Tomita H, Satoh M (2004) A new dynamical framework of nonhydrostatic
global model using the icosahedral grid. Fluid Dyn Res 34:357–400.
https://​doi.​org/​10.​1016/j.​fluid​dyn.​2004.​03.​003
van der Laan-Luijkx IT, van der Velde IR, van der Veen E, Tsuruta A, Stanislawska
K, Babenhauserheide A, Zhang HF, Liu Y, He W, Chen H, Masarie KA, Krol
MC, Peters W (2017) The CarbonTracker data assimilation shell (CTDAS)
v1.0: implementation and global carbon balance 2001–2015. Geosci
Model Dev 10(7):2785–2800. https://​doi.​org/​10.​5194/​gmd-​10-​2785-​2017
van der Werf GR, Randerson JT, Giglio L, van Leeuwen TT, Chen Y, Rogers BM,
Mu M, van Marle MJE, Morton DC, Collatz GJ, Yokelson RJ, Kasibhatla PS(2017) Global fire emissions estimates during 1997–2016. Earth Syst Sci
Data 9(2):697–720. https://​doi.​org/​10.​5194/​essd-9-​697-​2017
Yokota T, Yoshida Y, Eguchi N, Ota Y, Tanaka T, Watanabe H, Maksyutov S (2009)
Global concentrations of C
­ O2 and C­ H4 retrieved from GOSAT: first preliminary results. SOLA 5:160–163. https://​doi.​org/​10.​2151/​sola.​2009-​041

Page 19 of 19

Publisher’s NoteSpringer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.