Toward a long-term atmospheric CO2 inversion for elucidating natural carbon fluxes: technical notes of NISMON-CO2 v2021.1 – Progress in Earth and Planetary Science

NICAM-TM and NISMON-CO₂

Initially, we describe the fundamental features of the NISMON-CO₂ inversion system. One distinct feature of NICAM, on which NISMON-CO₂ 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, the NICAM-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 CO₂ (Niwa et al. 2011a, b). Simulations with NICAM-TM are performed iteratively in an inverse analysis of NISMON-CO₂. Further details about NICAM-TM are described by Niwa et al. (2017a).

The inverse analysis system NISMON-CO₂ uses 4D-Var to optimize surface CO₂ fluxes consistently with atmospheric observations (Niwa et al. 2017b). In 4D-Var, the cost function J, whose minimum state is sought, is defined as:

where δx represents the model parameter to be optimized, which comprises deviations from its prior estimate x^pri of surface fluxes and the initial global mole fraction offset. Meanwhile, M(.) is an operator representing the model calculation, and d is a vector defined by the observational vector y as ({mathbf{d}} = {mathbf{y}} – Mleft( {{mathbf{x}}^{{{text{pri}}}} } right)). The matrices B and R are the error covariance matrices of prior estimates and the model–observation mismatches, respectively. The iterative calculation of 4D-Var uses the gradient of the cost function:

to find the optimal model parameters (i.e., posterior estimates) that minimize the cost function. M^T in Eq. (2) represents the transpose of M, where M is the tangent linear operator of M(.). In fact, the relation ({mathbf{M}}Delta {mathbf{x}} approx Mleft( {Delta {mathbf{x}}} right)) is satisfied because the model calculates atmospheric transport linearly and does not include any chemical reactions. The calculations of M(.) and M^T are performed by the forward and adjoint simulations using NICAM-TM, which only simulates CO₂ transport in the so-called offline mode with prescribed meteorological data. The meteorological data are prepared in the “online mode,” in which meteorological fields are simulated by nudging horizontal winds toward those of the reanalysis data. In this study, we used the same model settings and prescribed 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 of approximately 45 km above the sea level.

The model operator M(.) includes spatiotemporal interpolation from the model grids and temporal steps to observational locations and times and the conversion of 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:

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 not in the model grid space (i.e., the icosahedral grids). The detailed conversion scheme is described in the next section. In NISMON-CO₂, the optimization is performed using a quasi-Newton method-based scheme termed the Preconditioned Optimizing Utility for Large-dimensional analyses (POpULar: Fujii and Kamachi 2003; Fujii 2005).

Grid conversion scheme

The icosahedral NICAM grid shapes hexagons or pentagons, which requires an elaborate conversion from latitude–longitude grids, because flux input datasets are originally prepared in latitude–longitude grids. Previous 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 sophisticated method 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 the thin light-gray mesh), i.e., rasterizing. Next, each finer-scale 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 to the 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:

where G_j denotes the original latitude–longitude grid that overlaps the ith icosahedral grid and N_i is the number of those grids (N_i = 9 in the case of Fig. 1). The partial area of each original latitude–longitude grid overlapped by the ith icosahedral grid is derived from the finer latitude–longitude grids, which is represented as a_i,j (individual colored areas in Fig. 1). The area of the ith icosahedral grid is denoted as A_i. 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 is entirely linear and thus represented by the matrix G in Eq. (3). This property enables us to easily develop the adjoint code represented by G^T in Eq. (4). By introducing G and G^T into the model operator and its adjoint, respectively, the control variables can be arrayed in latitude–longitude grids.

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 the conversion, fluxes in colored areas are allocated to the centered icosahedral grid

Optimizing fluxes in the latitude–longitude grid has two technical advantages. First, conversion of flux data from icosahedral grids to latitude–longitude grids after inversion for comparison with other datasets is not required. Such conversion often causes deterioration of mass 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 the high-resolution information originally contained in the prior flux data. During the grid conversion in this study, we divided 1° × 1° latitude–longitude grids by 10 × 10 to obtain finer grids.

Flux model

The surface CO₂ fluxes input to NICAM-TM can be described as:

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 f_fos, f_GPP, f_RE, f_LUC, f_fire, and f_ocn, respectively; they are prescribed using flux datasets and here 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, (beta_{{{text{GPP}}}}) and (beta_{{{text{RE}}}}), are scaling factors that produce diurnal variations. (beta_{{{text{GPP}}}}) and (beta_{{{text{RE}}}}) distribute fluxes at a three-hourly resolution from the monthly f_GPP and f_RE, respectively, whose values are derived from the downward shortwave radiation at the earth’s surface and 2-m height air temperatures data from the Japanese 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 (beta_{{{text{GPP}}}}) and (beta_{{{text{RE}}}}) to show different diurnal variations (e.g., (beta_{{{text{GPP}}}}) is zero at night because photosynthesis is inactive, which is not the case for (beta_{{{text{RE}}}})).

In this study, fossil fuel emissions were not optimized, i.e., they were fixed as f_fos, whereas the flux deviations of ({Delta }f_{{{text{GPP}}}}), ({Delta }f_{{{text{RE}}}}), and ({Delta }f_{{{text{ocn}}}}) were optimized. Furthermore, the LUC and biomass burning fluxes vary with scaling factors, whose deviations from 1, ({Delta }alpha_{{{text{LUC}}}}) and ({Delta }alpha_{{{text{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 the flux parameters ({Delta }f) and ({Delta }alpha), except for the ocean flux ({Delta }f_{{{text{ocn}}}}), whose temporal resolution is set annually. This is because ocean flux variations are smaller than terrestrial flux variations, and some erroneous values could arise in 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 in ocean fluxes.

In NISMON-CO₂ 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 f_fos, 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 f_GPP, f_RE, and f_LUC, respectively. The biomass burning emission data were obtained from the Global Fire Emissions Database (GFED) ver. 4.1s (van der Werf et al. 2017) for f_fire, and the air–sea CO₂ exchange data were taken from the Japan Meteorological Agency (JMA) (Iida et al. 2015, 2021) for f_ocn. Except for f_ocn, 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 in this 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 not intended 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 component is difficult because the inversion only uses observations of atmospheric CO₂ 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 carbonyl sulfide (Campbell et al. 2015).

Prior error covariance

A 4D-Var calculation requires not only a prior estimate but also its prior error covariance matrix, as denoted by B in Eqs. (1) and (2). In the inverse experiment by Niwa et al. (2017b), the prior error covariances (i.e., the off-diagonal elements of B) were simply designed with a Gaussian function commonly used for the globe, in which 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 as

where the subscript “cyc” represents cyclostationary, whose matrix is applied repeatedly for each year. The superscript k is an ensemble index, and ({updelta }{mathbf{x}}_{n}^{k}) describes the deviation of the flux of the kth member from the ensemble mean (i.e., the long-term mean). However, due to the limited 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 ({tilde{mathbf{B}}}_{{n,{text{ cyc}}}}) with a correlation matrix C as

Here, we use the Gaussian function for each element of C as

where the indices i and j represent different flux locations, whose distance is denoted by l_ij. L is a constant and set to 2000 km, which is long enough to only damp erroneous correlations in remote areas. In practice, we let C_ij = 0 where l_ij > 4000 km to minimize the computational burden. In Eq. (8), the scaling factor (gamma_{n}) is introduced as a 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 to be primarily generated by sampling errors. The technique described in Eq. (8) has been used as “the covariance localization” 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 covariance calculation (Eq. 7). After applying the correlation matrix, the resulting covariance was inflated by (gamma_{n} = 4) so that the globally integrated variance was approximately (2 Pg C yr⁻¹)². The defined uncertainty of 2 Pg C yr⁻¹ is derived from the difference between the global total net flux of the 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 the prior ocean flux, we used the JMA ocean flux data for 1990–2019. Because we optimized the annual means of 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 that of the terrestrial biosphere using Eqs. (8) and (9) with (gamma_{n} = 2). The resulting global variance was approximately (0.17 Pg C yr⁻¹)²; this number is rather arbitrary compared with that for the terrestrial flux, but it was derived after several trial experiments. When (gamma_{n}) was increased from 2, we obtained unrealistic variations in posterior ocean fluxes. Although 0.17 Pg C yr⁻¹ is too small as a global ocean flux uncertainty, such small prior errors are needed to stabilize ocean flux estimates.

For the other flux components such as LUC and biomass burning, we do not consider error correlations and define each prior flux error covariance matrix B_n,cyc as a diagonal matrix, whose diagonal elements (i.e., variance) are all set to 4.0; this yields errors of 200%. We combined all B_n,cyc and applied them to every year, whose values are arrayed in the same grid and resolution as those of the prior flux, i.e., 1° × 1°. Finally, we obtained the full matrix of B by including the variance of the global initial mole fraction offset of (0.5 ppm)².

Examples of the prior flux error covariance are given in Fig. 2, which shows error correlations (left( { = B_{ij} /left( {B_{ii}^{1/2} B_{jj}^{1/2} } right)} right)) of a given prior GPP for April and July. The original error correlations were significant in both months, even at a large 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 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 and f. In fact, these anisotropic correlation features arise from year-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.

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). Gray and red crosses represent correlations before and after localization (Eq. 8), respectively. Spatial structures of the localized error correlations of a and b 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

Observations

The observational data input to NISMON-CO₂ was derived from the dataset named Observation Package (ObsPack)-GLOBALVIEWplus and the related near-real-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-CO₂ ver. 2021.1, versions 6.1_2021_03-01 (Schuldt et al. 2021a) and 6.1.1_2021-05-17 (Schuldt et al. 2021b) of ObsPack-GLOBALVIEWplus and ObsPack-NRT were used, respectively; all data are listed in Niwa (2020). To demonstrate the potential performance of NISMON-CO₂ ver. 2021.1, we performed pseudo-observation inversions for this study by emulating observational data at actual sites and times. The locations of the observations are shown in 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 the entire 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 stations and tower sites, and shipboard observations). This is because NISMON-CO₂ 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.

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

Observation–model mismatch errors

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

where r_i denotes the standard deviation of the mole fractions around the ith observation and N_i represents the number 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, we introduce a scaling factor (beta) so that (chi^{2} left( {: = frac{{2J_{{{text{min}}}} }}{m}} right)), where m is the number of observations and J_min is the minimum of the cost function, should be < 1 (Tarantola 2005). We calculated r_i from the simulated daily mole fraction variations of CO₂ at the site of observation i ranging from one week prior and one week after the observational timing; the simulation was performed using NICAM-TM with prior fluxes. As denoted by Eq. (10), R is a diagonal matrix, which assumes that all observations are independent from each other. However, this is a simplistic assumption, specifically in cases where observations were obtained with high density. Therefore, we inflate the variances for such areas by introducing N_i, which reduces the weights of such dense observations in space and time. This variance inflation prevents continuous in situ observations from imposing exceedingly strong constraints compared with discrete flask sampling observations. For NISMON-CO₂ 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 same observational platforms; N_i places a single-site constraint on such collocated observations.

Pseudo-observation inversion experiment

To investigate the general performance of NISMON-CO₂ 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” and the other one was used as a prior estimate. First, we performed a forward transport simulation of atmospheric CO₂ using the true flux dataset and extracted the simulated CO₂ mole fractions at the above-mentioned observational points and times. Next, we performed inverse experiments incorporating those pseudo-observations. Niwa et al. (2017b) conducted a similar experiment; however, the experiment in this study was designed more practically, using a combination of continuous in situ and discrete flask measurements (Fig. 3) and covering a longer analysis period. Furthermore, by virtue of the implemented 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 pseudo-observations (true run) employed the same NICAM-TM settings used for the inversion, except for the horizontal resolution. The horizontal resolution is “glevel-6,” whose mean 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 from the posterior fluxes of NISMON-CO₂ ver. 2021.1 (Niwa 2020), which were generated from real observations; therefore, the generated pseudo-observations of atmospheric CO₂ mole fractions have realistic seasonal and interannual 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 those in NISMON-CO₂ ver. 2021.1 (described in Sect. 2.3). For the terrestrial fluxes (f_GPP, f_RE, f_LUC, and f_fire), we employed the same VISIT and GFEDv4.1s datasets as in NISMON-CO₂ ver. 2021.1, but they were climatologically averaged for 2006 to 2018. Meanwhile, the ocean fluxes (f_ocn) were based on different ocean flux data from Takahashi et al. (2009), which provide climatological fluxes for the reference year 2000. Therefore, excluding the fossil fuel emissions, the prior fluxes do not contain interannual variations.

The inversion calculation was performed from January 2007 to March 2018, with the first 12 months and last three months disregarded as the spin-up and spin-down times of the analysis, respectively. Therefore, the inversion 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) for the ease of performing the multiple inverse experiments described in the next section.

Sensitivity tests

In addition to the control experiment performed with the settings 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 with the 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., N_i = 1 in Eq. (10). In this study, for a fair comparison with the CTL case, we changed β in Eq. (10) so that the cost function (Eq. 1) is nearly equivalent to that of CTL at the beginning of the iterative calculation. The fourth inversion (referred to as ISO_ERR) employed the isotropic prior error covariance that was simply defined using a Gaussian function with correlation scale lengths of 500 km and 1000 km for land and ocean fluxes, respectively. In this experiment, also for a fair comparison, each (gamma_{n}) in Eq. (8) was tuned so that the global integrated value of the resulting global prior error covariance was nearly equivalent to that of CTL.