the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Order of magnitude wall time improvement of variational methane inversions by physical parallelization: a demonstration using TM54DVAR
Sander Houweling
Arjo Segers
Atmospheric inversions are used to constrain emissions of trace gases using atmospheric molefraction measurements. The fourdimensional variational (4DVAR) inversion approach allows optimization of emissions at a higher temporal and spatial resolution than ensemble or analytical approaches but provides limited opportunities for scalable parallelization because it is an iterative optimization method. Multidecadal variational inversions are needed to optimally extract information from the long measurement records of longlived atmospheric trace gases like carbon dioxide and methane. However, the wall time needed – up to months – complicates these multidecadal inversions. The physical parallelization (PP) method introduced by Chevallier (2013) addresses this problem for carbon dioxide inversions by splitting the period of the chemical transport model into blocks and running them in parallel. Here, we present a new implementation of the PP method which is suitable for methane inversions accounting for the chemical sink of methane. The performance of the PP method is tested in an 11year inversion using a TM54DVAR inversion setup that assimilates surface observations to optimize methane emissions at grid scale. Our PP implementation improves the wall time performance by a factor of 5 and shows excellent agreement with a full serial inversion in an identical configuration (global mean emissions difference =0.06 % with an interannual variation correlation R=0.99; regional mean emission difference <5 % and interannual variation R>0.94). The wall time improvement of the PP method increases with the size of the inversion period. The PP method is planned to be used in future releases of the Copernicus Atmosphere Monitoring Service (CAMS) multidecadal methane reanalysis.
Methane (CH_{4}) is the secondmost important greenhouse gas after carbon dioxide (CO_{2}), and its atmospheric abundance has increased by more than a factor of 2.5 since preindustrial times. Methane is responsible for 25 % of the anthropogenic radiative forcing despite its 200 times lower abundance than CO_{2}, due to its strong global warming potential (Myhre et al., 2013). Atmospheric inversions provide methane emission estimates by optimally combining the information in atmospheric observations and bottomup emissions (estimates using processbased models and inventories) along with corresponding error characteristics. Chemical transport models (CTMs) simulate the spatiotemporal distribution of the methane mole fractions in the atmosphere for a given set of emissions while also accounting for its atmospheric sink. Inversions use CTMs to disentangle the influences of atmospheric transport from the influences of emissions and sinks on the observed mole fractions (Naus et al., 2019; Pandey et al., 2019). A few studies have performed inversions on multidecadal scales to constrain emissions using the long measurement record of methane mole fractions. For example, “The Global Methane Budget 2000–2017” (Saunois et al., 2020) presents regional emission estimates from nine different inversion setups. The methane emissions reanalysis project under the Copernicus Atmosphere Monitoring Service (CAMS) performs multidecadal inversions using the TM54DVAR variational approach to provide regularly updated gridded methane emissions (Segers and Houweling, 2020).
Trace gas inversions adjust a state vector, which includes gridded emissions (and sometimes initial molefraction field and other parameters), to improve the agreement between model simulations and observations. The inversions minimize a Bayesian cost function that is defined based on the difference between the modeled and observed mole fractions as well as the magnitude of the emission adjustments, weighted with respective error covariances. There are three main approaches used in atmospheric inverse modeling: analytical, ensemble, and variational. The analytical approach is based on a closedform solution of Bayes' theorem (Gurney et al., 2002). It requires the calculation of observation sensitivities of each of the state vector elements separately. The large computational cost involved restricts its application to smallsize state vectors. The ensemble approach improves the computational performance by parameterizing the state vector sensitivities using a statistical ensemble (Peters et al., 2005). Still only a relatively smallsize state vector can be afforded using this approach.
The variational inversion approach was introduced to lift the state vector size restriction (Chevallier et al., 2005). In this approach, the cost function minimum is computed using an iterative procedure, with each iteration comprising a forward and an adjoint CTM run. The variational approach can be applied to weakly nonlinear inverse problems using a suitable steepest decent numerical minimizer (Naus et al., 2021; Pandey et al., 2016). Truncated posterior uncertainties can be obtained from variational inversions using the conjugate gradient minimizer for linear inverse problems (Meirink et al., 2008). A more robust but computationally expensive estimate of posterior uncertainties can be obtained using a Monte Carlo method (Chevallier et al., 2007; Pandey et al., 2016). However, as each iteration of a variational inversion uses the output of the previous iteration, the calculations can take months, depending on the spatial and temporal resolution. This long wall time limits the resolution, duration, and number of iterations that can be used in multidecadal variational inversions. To reduce this long wall time for CO_{2} inversions, Chevallier (2013) introduced the physical parallelization (PP) method. In this method, the full inversion period is split into a number of blocks, and the CTM runs for the blocks are performed in parallel within each iteration. Corrections are added to the simulated CO_{2} mole fractions in a block to account for emission adjustments (iteration minus prior emissions) in earlier blocks. This method reduced the wall time by an order of magnitude (sevenfold improvement for a 32year inversion), while keeping the inversionderived emission adjustments statistically consistent with a serial inversion. However, the original implementation of the PP method cannot be used for a reactive trace gas like methane as the method does not account for the atmospheric chemical sink.
Here, we present an improved PP method that accounts for the limited atmospheric lifetime of reactive trace gases such as methane, which has an atmospheric lifetime of about 9 years (mainly due to oxidation by OH radicals). The intention is to use this new PP implementation for the CAMS methane flux reanalysis, which aims to provide annually updated multidecadal emission estimates, within a production cycle of only a few months. In Sect. 2, we present our PP method. The method's performance is tested using an 11year inversion setup presented in Sect. 3. The wall time and optimized emissions of a PP inversion are compared to a serial inversion in an identical configuration. In Sect. 4, we discuss possible future improvements and applications of the PP method. Our conclusions are summarized in Sect. 5.
An inversion of an atmospheric trace gas minimizes a Bayesian cost function of the state vector x:
Here, y is the observation vector and x^{b} is the prior state vector. B and R are the error covariance matrices of the prior emissions and the observations, respectively. The vector m constitutes the modeled mole fractions corresponding to y; m is computed using a CTM operator H, which simulates the mole fractions given the emissions in the state x and the initial mole fraction field c_{0}:
In a variational inversion setup, the posterior solution of Eq. (1) is obtained by minimizing J using an iterative procedure that computes a new emission update x^{i+1} in each iteration i using the gradient:
H^{∗} represents the adjoint CTM operator, which is implemented using the adjoint code of the CTM. The inversion finishes when a predefined convergence criterion is met, such as a desired gradient norm reduction or simply a maximum number of iterations.
In the PP method presented in Chevallier (2013), the full period of the inversion is split into r overlapping time blocks, which can be run in parallel. Figure 1 schematically represents the main steps in the PP method used in the forward mode to calculate m^{i}. At the start of the inversion, a serial CTM run (without segmentation) is performed to calculate initial molefraction fields ${\mathit{c}}_{k}^{\mathrm{b}}$ for each block k using the prior emissions x^{b}. In an iteration, the block mole fractions for the iteration ${\mathit{m}}_{k}^{i}$ are computed using the block CTM operator H_{k}, the iteration emissions for the block ${\mathit{x}}_{k}^{i}$, the initial mole fraction for this block ${\mathit{c}}_{k}^{\mathrm{b}}$, and a mole fraction correction ${n}_{k}^{i}$:
Here, the scalar ${n}_{k}^{i}$ accounts for the global mean molefraction changes due to emission differences (x^{i}−x^{b}) during the inversion period that precedes the block k. The error due to this simplification is further reduced by using an overlap period between consecutive blocks, where modeled mole fractions from the succeeding block are discarded. The overlap period CTM run distributes the emission differences uniformly through the earth's atmosphere. The PP method by Chevallier (2013) was applied to CO_{2} inversions, where the scalar molefraction correction ${n}_{k}^{i}$ for block k was simply calculated as the sum of the emission differences from each preceding block (i.e., block 1, 2, 3 … k−1):
Here, E denotes a summation matrix used to compute the global sum of the elements of x_{l} and f is a scalar used to convert emissions to mole fractions assuming a uniform distribution of the emitted trace gas throughout the earth's atmosphere; f is calculated simply as the ratio between the number of moles in a unit emission and the number of moles of air in the atmosphere.
Methane has an atmospheric lifetime of about 9 years. Unlike CO_{2}, the molefraction impact of methane emission differences will be reduced due to atmospheric chemistry within the duration of a multidecadal inversion as well as within a PP inversion block. Therefore, in our new implementation of the PP, we use a molefraction correction vector ${\mathit{n}}_{k}^{i}$ (with size of ${\mathit{m}}_{k}^{i}$) instead of the scalar ${n}_{k}^{i}$ to apply separate corrections to each observation. We account for the limited lifetime of methane by implementing an atmospheric sink operator S. In addition, we use a CTM block sensitivity vector h_{k} to distribute global emission changes more precisely, taking into account the full 3D atmospheric transport and the sink rather than assuming a globally uniform distribution. Herein, h_{k} is computed at the start of the inversion by running H_{k} in forward mode with an uniform initial mole fraction field and zero emissions, i.e., ${\mathit{h}}_{k}=\phantom{\rule{0.125em}{0ex}}{H}_{k}\left({\mathit{c}}_{k}=\mathrm{1},\phantom{\rule{0.125em}{0ex}}{\mathit{x}}_{k}=\mathrm{0}\right)$; ${\mathit{n}}_{k}^{i}$ is computed as
Here, the scalar s_{k,l} accounts for the impact of atmospheric sinks on the global uniform molefraction change during the time period between the blocks k and l; s_{k,l} is generated using a sink operator S. We describe a formulation of S in the next section. Within block k itself, the impact of atmospheric sinks is accounted for by h_{k}.
Each iteration of a variational inversion computes a departures vector δm:
The adjoint CTM H^{∗} is run with δm to compute the local gradient of the cost function (Eq. 2). In the PP method, H^{∗} is split into blocks covering the same periods as used for the forward CTM simulation. In an iteration, each adjoint block is first run with the respective departures. Then, the modeled adjoint sensitivities of a block $\mathit{\delta}{\mathit{x}}_{k}^{i}$ are adjusted for the effects of departures of succeeding blocks by adding the adjoint molefraction correction scalar ${g}_{k}^{i}$:
Here ${\mathit{h}}_{l}^{T}\phantom{\rule{0.125em}{0ex}}\mathit{\delta}{\mathit{m}}_{l}^{i}$ is the matrix dot product of the two vectors, both of which have the same size. The correctness of the adjoint implementation of the PP method can be verified using the adjoint test (Meirink et al., 2008). The test checks whether the equality
is satisfied to an accuracy near the computing precision. M and M^{∗} denote the forward and adjoint model operators, 〈〉 denotes the inner product, and p and q are the arbitrary forward and adjoint model states.
In a PP inversion, the initial molefraction field c_{0} needs to be consistent with the observations, as a discrepancy between the two leads to large emission differences in the early months of the inversion period. This issue can easily be dealt with in a serial inversion using a spinup period and rejecting this period from the posterior solution. However, in a PP inversion, the large emission differences may result in large molefraction corrections, which increases the error in the PP approximation (see Eqs. 4 and 5). This can be avoided by taking a realistic c_{0} from the posterior molefraction simulations of another inversion covering the period before the PP inversion. If such an inversion is not available, c_{0} can be computed by performing an inversion for the 1year period preceding the PP inversion.
In summary, the main steps of the PP methane inversion are as follows:

Construct an initial molefraction field c_{0} consistent with observations at the start of the inversion.

Split the full period of the inversion into r overlapping time blocks.

To calculate ${\mathit{c}}_{k}^{\mathrm{b}}$, run the forward CTM serially with prior emissions x^{b} and save the simulated molefraction fields at the start time of each block.

Calculate the CTM block sensitivity vector h_{k} by running the CTM over each block with a uniform initial molefraction field of 1 and zero emissions, and sample the model output at the observation time and locations.

Prepare a sink operator S which accounts for the impact of atmospheric sinks on methane mole fractions during a period.

Perform the inversion by iteratively minimizing the cost function until the convergence condition is met using a forward and an adjoint run in each iteration:
 a.
Forward run:
 i.
Run all forward CTM blocks in parallel with the initial molefraction fields from step 3.
 ii.
Account for the emission changes relative to the prior emissions in preceding blocks by adding the corrections ${\mathit{n}}_{k}^{i}$ (Eq. 5).
 i.
 b.
Adjoint run:
 i.
Run all adjoint CTM blocks in parallel to calculate the adjoint emission sensitivities.
 ii.
Add the adjoint correction ${g}_{k}^{i}$ to account for the effect of departures in successive blocks (Eq. 8).
 i.
 a.
The CTM runs in steps 4, 6.a.i, and 6.b.i are performed in parallel. The steps without a CTM run (1, 2, 5, 6.a.ii, and 6.b.ii) require very little wall time. Step 3 is the most time consuming because a full serial CTM run is performed in the step.
We evaluate the performance of the PP method by comparing a PP inversion with a serial methane inversion. Both inversions are performed for an 11year period (1999–2010) with identical observations and prior emissions. We use the TM54DVAR inversion system (Bergamaschi et al., 2010; Meirink et al., 2008; Krol et al., 2005) with the settings used in Pandey et al. (2016). The TM5 CTM is run at $\mathrm{6}{}^{\circ}\times \mathrm{4}{}^{\circ}$ horizontal resolution and 25 vertical hybrid sigmapressure levels from the surface to the top of the atmosphere. The meteorological fields for this offline model are taken from the European Centre for MediumRange Weather Forecasts (ECMWF) ERAInterim reanalysis (Dee et al., 2011). We optimize a single category (“total”) of methane emissions at $\mathrm{6}{}^{\circ}\times \mathrm{4}{}^{\circ}$ spatial resolution and monthly temporal resolution. The posterior emissions of the two inversions are compared after integrating over the TRANSCOM regions shown in Fig. 2a.
The inversion assimilates surface observations from the NOAA Earth System Research Laboratory (ESRL) global cooperative air sampling network at on and offshore sites (Dlugokencky et al., 2011, 2020). The locations of the observation sites are shown in Fig. 2b. The prior covariance matrix B is constructed as follows: the diagonal elements of B are constructed assuming ±1σ uncertainties of 50 % of the emissions per grid cell per month. The offdiagonal elements are constructed by assuming the emissions to be correlated temporally using an exponential correlation function with an efolding timescale of 3 months, and spatially with a Gaussian correlation function using a length scale of 500 km (Houweling et al., 2014). Uncertainties of 1.4 ppb are assigned to methane observations. Our system also assigns a modeling representation error based on simulated local molefraction gradients (Basu et al., 2013). The prior emissions are taken from the same sources as in Pandey et al. (2016). There is no interannual variability in the prior emissions because the 2008 emissions are used for every year. The cost function J is minimized using the conjugate gradient minimizer, which is based on the Lanczos algorithm (Fisher and Courtier, 1995). The inversions use the convergence criterion of gradient norm reduction by a factor 1000, which is achieved after 19 iterations in both inversions.
In the PP inversion, we split the inversion period of 1999–2010 into 11 blocks of 21 months. The first 9 months of each block is the overlap period used for uniformly mixing the emission changes within the atmosphere, while the last 12 months provide modeled mole fractions for assimilating the observations. We parameterize the sink operator S, which computes the sink scaling factor s_{k,l} (Eq. 6), with an efolding decay function and a constant atmospheric lifetime of methane (τ) of 9 years.
Here, t_{l} and t_{k} are the start times of the blocks l and k, respectively. We found this simple parameterization with a constant lifetime is sufficient for our test inversion.
The input emissions of TM5 are mass fluxes (Tg yr^{−1}) and the output is in mole fractions (ppb). The methane emission changes are converted into mole fractions using an $f=\mathrm{0.361}\frac{\mathrm{ppb}}{\mathrm{Tg}}$. Successful implementation of the PP method in the adjoint mode was verified using the adjoint test (Eq. 9).
3.1 PP inversion errors
Herein, we evaluate the difference in modeled mole fractions and posterior emissions between the PP and serial inversions.
3.1.1 Molefraction errors
We first examine the quality of the inversionoptimized fit to the observation. Figure 3 shows the time series of the prior and posterior model simulations and the observations for two background sites representing each hemisphere: Barrow (Alaska) and the South Pole. The root mean square difference (RMSD) of the prior model with observations for Barrow (78 ppb) is three times higher than for the South Pole (28 pbb). Barrow observations show more highfrequency variations than the South Pole, as the Northern Hemisphere station is influenced by methane emissions from wetlands. The mole fractions simulated by the PP inversion are in good agreement with the results obtained from the serial inversion: RMSDs of 2 and 1 pbb for Barrow and the South Pole, respectively, which are only 2.5 % and 3.2 % of the prior RMSD. This shows that the PP inversion, starting from identical prior emissions, is able to match the observations at these sites as well as the serial inversion.
Figure 4 shows the average molefraction differences at all observation sites. The observation–prior RMSD for all observations combined is 67 ppb. The mean mismatch is −58 ppb because the 2008 bottomup emissions used as the prior emissions are larger than the mean posterior emission over 1999–2010. The average data uncertainty (mean of the square root of the diagonal elements of R) is 19 ppb (not shown). For both inversions, a good model fit (90 % reduction in mean of observation–model mole fraction mismatch) to the observations is achieved with a gradient norm reduction of 1000. The posterior simulation of both the serial and PP inversions reduce the RMSD to 20 ppb (mean $=\mathrm{2}$ ppb). The RMSD between PP and serial is 1.9 ppb (mean $=\mathrm{0.1}$ ppb), which is an order of magnitude smaller than the posterior–prior RMSD of 62 ppb (mean $=\mathrm{55}$ ppb). This shows that the implementation of the PP method has little impact on the inversion's ability to fit the observations.
3.1.2 Posterior emission errors
A good agreement between observations and posterior models does not guarantee that the inversions have produced similar posterior emissions. The physically parallelized CTM used in the PP inversion has lost some of the consistency of the full CTM used in the serial inversion and the PP emission errors will depend on the impact of this CTM simplification. Figure 5 shows mean emissions (averaged over 1999–2010) from the inversions integrated over the globe and TRANSCOM regions. We do not have a good estimate of the posterior uncertainties because a large number of variational inversion iterations are needed for the second derivative of the cost function to converge. Therefore, we evaluate PP performance by comparing the PP–serial emission differences against the emission adjustments performed by the serial inversion (serial–prior differences) and prior emission uncertainties. The serial inversion adjusts the global mean prior emissions of 544±11 Tg yr^{−1} by −22 Tg yr^{−1}. The PP inversion is in excellent agreement with the serial inversion in this respect. The two differ by 0.3 Tg yr^{−1} (0.06 %), which is 1 % of the difference between the prior and serial emissions. The global methane emissions are in general well constrained by the NOAA observations in a serial inversion, and the additional error introduced by the PP method only causes a 1 % error relative to the serial–prior emission mismatch. At regional scales, the serial inversion adjustment is the smallest for Australia: +0.4 Tg yr^{−1} from the prior emissions of 6.6±0.4 Tg yr^{−1}. The PP inversion adjusts the prior emissions here by +0.5 Tg yr^{−1}. The serial inversion changes the Eurasian temperate emissions the most, by −58 Tg yr^{−1}, where prior emissions are 135±8 Tg yr^{−1}. The PP inversion changes these emissions by −60 Tg yr^{−1}, i. e., a difference of 2 Tg yr^{−1}. The South American temperate region has the largest PP error relative to the serial–prior difference, of 2 Tg yr^{−1}. The serial emissions for this region are 6.5 Tg yr^{−1} higher than the prior emissions of 36±2.4 Tg yr^{−1}. In summary, mean PP emission estimates for the TRANSCOM regions deviate within <5 % from the prior emissions, while the serial–prior differences are up to 50 % of the prior emissions.
Figure 6 shows the interannual variability of the emission estimates. The global emissions time series of the PP and serial inversions show a very good agreement with a correlation coefficient R=0.99, explained by the large observational constraint. Over the TRANSCOM regions, the North American temperate region has the best agreement (R=1.0). All other regions have R higher than 0.98 except for Australia (0.96) and Europe (0.94). Figure 7 shows the intraannual variations of the emissions. At the global scale, the PP and serial time series match very well with R=1.00, whereas R between prior and serial emission variations is 0.93. The agreement between PP and serial time series is also very good for all TRANSCOM regions (R>0.98), despite low correlations between prior and serial emissions for some regions, e.g., R=0.13 for the South American temperate region. This shows that the PP inversion is able to reproduce the seasonal cycle of the emissions very well. Figure 8 shows the spatial distribution of the emission differences at grid scale. The mean (±1σ spread) of the differences between the serial and prior emissions is $\mathrm{8}\phantom{\rule{0.125em}{0ex}}\times \phantom{\rule{0.125em}{0ex}}{\mathrm{10}}^{\mathrm{3}}\phantom{\rule{0.125em}{0ex}}(\pm \mathrm{0.5})$ Tg gridbox^{−1} yr^{−1}, and is $\mathrm{9}\phantom{\rule{0.125em}{0ex}}\times \phantom{\rule{0.125em}{0ex}}{\mathrm{10}}^{\mathrm{5}}\phantom{\rule{0.125em}{0ex}}(\pm \mathrm{0.04})$ Tg gridbox^{−1} yr^{−1} for serial and PP inversions. Emission differences between the PP and serial inversions are visible over India and the South American temperate region. These differences are likely due to the lack of observational constraint in these regions (see Fig. 2). In summary, the combination of small differences in the mean emissions and the high correlations between intra and interannual time series shows that PP inversion can effectively reproduce results of serial inversion at regional scales.
3.2 Wall time
Table 1 compares the wall times used by the PP and serial inversions. The TM5 model in our inversion setup uses OpenMP parallelization and gives the best wall time performance on four central processing units (CPUs) on a single node (12core 2.6 GHz Intel Xeon E52690 v3). Using more CPUs reduces the performance, as the communication overhead within the CPUs becomes the bottleneck (note that the TM5MP version described in Williams et al., 2017, with improved parallel scaling, was not used in this study). In this configuration, a forward or adjoint TM5 CTM run of 1 year took about 15 min. Hence, an iteration of the serial inversion, consisting of 11 years of forward and adjoint runs, required 5 h. The PP inversion iterations were performed in 11 parallel blocks of 21 months each on four CPUs. A single PP iteration took 55 min, which is >5 times faster than the serial inversion. The main steps of PP implementation are listed in Sect. 2. In our inversion test, the initial molefraction fields c_{0} (step 1) were taken from an inversion using surface measurements that was not performed in this study. Steps 1, 2, 5, 6.a.ii, and 6.b.ii took negligible time. Step 3 took 2.5 h because it consists of a full serial TM5 forward run. Steps 4, 6.a.i, and 6.b.i consist of 11 TM5 run over blocks of 21 months which were run in parallel and took 25 min each. Note that an iteration took longer than the sum of the forward and adjoint block runs because of a few minutes waiting time for the computer cores to become available again. In total, the PP inversion took 20 h, 5 times less than the serial inversion which took an inversion's wall time of 101 h. Note that although the PP inversion took a shorter wall time, it needed extra CPU core hours for the additional 9month overlap, CTM block sensitivity, and initial molefraction computation runs. The PP inversion used a total of 700 CPU core hours, whereas the serial inversion used about 400 CPU core hours. Table 1 also provides a projection of the wall time improvement of a hypothetical 35year inversion (not performed in this study) based on the TM54DVAR inversion setup used in this study. A PP inversion would be 15 times faster for such a long period. Overall, we find that the PP method, which accounts for the atmospheric lifetime of methane, is able to effectively reproduce the posterior emissions of a 11year conventional serial inversion 5 times faster.
4.1 PP method applications
The utility of the PP method for inversion of a trace gas depends on the timescale of the influence of emissions on observations within the spatial domain of the CTM. Therefore, PP is mainly useful in global inversions of trace gases that have an atmospheric lifetime of a year or longer in the atmosphere. For a trace gas with a shorter lifetime, such as carbon monoxide with a 2 month lifetime, emission perturbations last for a short duration. A multidecadal inversion of such a trace gas can be broken down into many short inversions. These short inversions can be performed in parallel, and the posterior emission can be combined thereafter. A similar approach can be used for regional inversions of short and longlived trace gases because emission perturbations are quickly advected out of the regional CTM domain and hence do not influence observations for a long period.
The hydroxyl radical OH is the main sink of methane in the atmosphere. Zhang et al. (2018) showed that the satelliteobserved atmospheric signature of the methane sink is sufficiently distinct from that of methane emissions; hence, OH mole fractions can be optimized using synthetic shortwave infrared (SWIR) and thermal infrared (TIR) satellite observations. Following up on this, Maasakkers et al. (2019) and Zhang et al. (2021) used methane observations from the GOSAT satellite to optimize atmospheric OH fields along with methane emissions. The simultaneous optimization of OH with methane emissions introduces a nonlinearity into the inversion because methane loss rate depends on the product of methane and OH mole fractions. However, the changes to the methane mole fractions are expected to remain small during the inversions. Hence, the nonlinear effect is small and a quasilinearity is assumed to solve the inversion analytically using computation of the full Jacobian matrix of the CTM. Under a quasilinearity assumption, OH can be optimized in a PP methane inversion by introducing annual OH scaling factors into the state vector and the methane lifetimes in the sink operator can be scaled in each iteration to reflect the corresponding OH adjustments. Such an implementation can also be used in inversions optimizing OH using methyl chloroform (Naus et al., 2021).
4.2 Possible further improvements
The PP method accounts for changes in the background mole fractions due to emission changes in preceding blocks using a sink operator S, a CTM block sensitivity h, and an overlap between the consecutive blocks. In our test experiment, S is assumed to be an efolding decay function with an atmospheric lifetime of methane of 9 years, which we found to be sufficient for the annually repeating OH field used in our 11year CTM runs. Methane lifetime within the duration of a longer multidecadal inversion will vary due to climatological influences as well as possible trends and interannual variations in hydroxyl radical abundance. In such cases, S can be defined as a function of an annual lifetime vector for the specific CTM run. The lifetime vector can be calculated as the ratio of the annual sink and global methane burden simulated by the serial CTM run in step 3 of the PP method.
The overlap period between consecutive blocks in the PP method allows methane emission perturbations to mix within the CTM domain according to atmospheric transport. We used a 9month overlap in our test inversion setup. It was sufficient to estimate the total emissions from TRANSCOM regions using the surface observations. The 6month overlap used by Chevalier (2013) for CO_{2} inversion was found to be insufficient for a PP methane inversion, likely because of the differences between the source and sink distributions of methane and CO_{2}. Increasing the overlap period to 9 months and using a CTM block sensitivity vector solved this issue. We expect that a 1year overlap, equal to the interhemispheric mixing time, would be more than sufficient for all tracers, irrespective of their source–sink distribution and lifetime. A shorter overlap would improve the computational efficiency and wall time but reduce the accuracy of the physical parallelization of the CTM. The PP accuracy could be maintained with shorter overlap periods by using a molefraction correction vector per hemisphere rather than the single global vector used in this study. However, the computational resources and wall time saved by this would be partially spent on the additional block sensitivity runs. Our test inversions are performed at a relatively coarse horizontal resolution of $\mathrm{6}{}^{\circ}\times \mathrm{4}{}^{\circ}$ with 25 vertical hybrid sigmapressure levels. We do not expect the performance of the PP method to degrade significantly for higherresolution inversions if there is sufficient overlap between the blocks and the molefraction corrections are parameterized correctly. Furthermore, the performance gained by performing the inversions at higher resolution will likely outweigh the accuracy loss due to the assumptions made in the PP method because of the improved computational performance.
The PP method reduces the wall time of the CTM simulations in a variational inversion but introduces additional model errors because of the simplifications made. For our test inversion setup, these PPCTM model errors are minor, as the posterior PP emission estimates are in good agreement with the serial estimates. In future PP implementations, these PP–CTM errors can be accounted for in the observation error matrix R. The PP–CTM error can be calculated as the difference between the model output of a PP and a serial forward CTM run with randomly perturbed prior emissions.
4.3 Current CAMS inversion setup
In the future, the PP method will be implemented in the CAMS multidecadal methane emissions reanalysis setup. The European Commission has anticipated the need for reliable information about atmospheric composition of greenhouse gases through development of numerical systems that combine sophisticated physical models with measurements from a wide range of observing systems for an operational service, which is being implemented. The current CAMS methane flux reanalysis product (Segers and Houweling, 2020) uses the TM54DVAR inverse modeling system and provides measurementinformed monthly methane emission estimates. The latest release has two sets of methane emissions: (1) release v19r1 for 1990–2019 using surface observation; (2) release v19r1s for 2010–2019 using surface and also GOSAT satellite observations. The surface observations are mainly from the NOAA network (Dlugokencky et al., 2011). Methane emissions are optimized at $\mathrm{3}{}^{\circ}\times \mathrm{2}{}^{\circ}$ spatial resolution and monthly temporal resolution using TM5 with 34 vertical layers. If performed in serial mode, each iteration of the 1990–2019 inversion would take about 5–10 d, and the full inversion will require multiple months to finish. Segers and Houweling (2020) circumvent this issue by breaking down the full inversion into smaller inversions of 3year time windows that are performed in parallel. The target inversion at high resolution ($\mathrm{3}{}^{\circ}\times \mathrm{2}{}^{\circ}$, 34 layers) is preceded by a coarseresolution inversion ($\mathrm{6}{}^{\circ}\times \mathrm{4}{}^{\circ}$, 25 layers) that provides the initial molefraction fields and is processed serially. The highresolution inversion optimizes only the emissions and uses initial mole fractions for each 3year block obtained from molefraction fields of a coarseresolution inversion, which optimizes both emission and initial mole fractions. The 1990–2019 inversion using this approach still takes 3–4 months to finish and requires about 40 smaller inversions to provide the end result. These numbers of course depend on the parallel efficiency of the model and the computing server. The need for a coarseresolution serial sequence of inversions to provide initial molefraction fields limits the inversion period for which this method can be used. With the implementation of the PP method presented in this study, the wall time performance of the CAMS reanalysis inversions will improve in the future.
Regular surface observations of methane mole fractions started in early 1984, and by now the measurement record spans more than 35 years (Dlugokencky et al., 2011). An atmospheric inversion with a very large state vector is needed optimize emissions using such long measurement records at a grid scale. The variational inversion approach allows for optimization of a much larger state vector than the ensemble or analytical approaches. However, each iteration of a variational inversion uses the output of the previous iteration, limiting the opportunity for scalable parallelization. At the same time, the increase in the spatiotemporal resolution of CTMs needed to take full advantage of the rapidly improving precision and coverage of surface and satellite measurements results in a rapid increase in wall time.
We have developed the PP method for methane inversions, which improves the wall time of variational methane inversions by physical CTM parallelization while accounting for the atmospheric lifetime in forward and adjoint variational modes. We have tested the performance of this method using an 11year TM54DVAR inversion setup that consists of a conventional serial inversion and a PP inversion in an identical configuration. The PP method reduced the wall time by a factor of 5, while still showing excellent agreement with the posterior emissions from the serial inversion. The wall time improvement of using PP will be even larger for longer inversions, e.g., by a factor of 15 for a 35year inversion. The PP method makes multidecadal global inversions of longlived atmospheric trace gases more feasible. It will be implemented in the CAMS reanalysis setup which provides regular updates of multidecadal emission estimates by assimilating surface and satellite observations.
The TM54DVARPP version 1.0beta1 code used in this study for the simulations can be downloaded from Zenodo (https://doi.org/10.5281/zenodo.6326373, Pandey et al., 2022). The TM5 model is described in detail on http://tm5.sourceforge.net/ (Krol et al., 2022).
NOAA ESRL methane observations used in this study are available on Zenodo in the input folder of the TM54DVARPP code (https://doi.org/10.5281/zenodo.6326373, Pandey et al., 2022).
The study was designed by SH and AS. SP and AS developed the PP. SP implemented the PP method on TM54DVAR and did the performance test simulations. SP wrote the paper using contributions from all the coauthors.
The contact author has declared that neither they nor their coauthors have any competing interests.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
We thank the efforts of NOAA and other surfaceobservations networks for producing and maintaining the vital long record of global methane observations. The computations for this study were carried out on the Dutch national supercomputer Cartesius (https://userinfo.surfsara.nl/systems/cartesius; last access: 1 June 2022) maintained by SURFSara. CAMS is implemented by the European Centre for MediumRange Weather Forecasts on behalf of the European Commission.
This research has been supported by the Copernicus Atmosphere Monitoring Service (grant no: CAMS73).
This paper was edited by Slimane Bekki and reviewed by three anonymous referees.
Basu, S., Guerlet, S., Butz, A., Houweling, S., Hasekamp, O., Aben, I., Krummel, P., Steele, P., Langenfelds, R., Torn, M., Biraud, S., Stephens, B., Andrews, A., and Worthy, D.: Global CO_{2} fluxes estimated from GOSAT retrievals of total column CO_{2}, Atmos. Chem. Phys., 13, 8695–8717, https://doi.org/10.5194/acp1386952013, 2013.
Bergamaschi, P., Krol, M., Meirink, J. F., Dentener, F., Segers, A., van Aardenne, J., Monni, S., Vermeulen, A. T., Schmidt, M., Ramonet, M., Yver, C., Meinhardt, F., Nisbet, E. G., Fisher, R. E., O'Doherty, S., and Dlugokencky, E. J.: Inverse modeling of European CH_{4} emissions 2001–2006, J. Geophys. Res., 115, D22309, https://doi.org/10.1029/2010JD014180, 2010.
Chevallier, F.: On the parallelization of atmospheric inversions of CO_{2} surface fluxes within a variational framework, Geosci. Model Dev., 6, 783–790, https://doi.org/10.5194/gmd67832013, 2013.
Chevallier, F., Fisher, M., Peylin, P., Serrar, S., Bousquet, P., Bréon, F. M., Chédin, A., and Ciais, P.: Inferring CO_{2} sources and sinks from satellite observations: Method and application to TOVS data, J. Geophys. Res.Atmos., 110, 1–13, https://doi.org/10.1029/2005JD006390, 2005.
Chevallier, F., Breon, F.M., and Rayner, P. J.: Contribution of the Orbiting Carbon Observatory to the estimation of CO_{2} sources and sinks: Theoretical study in a variational data assimilation framework, J. Geophys. Res.Atmos., 112, D09307, https://doi.org/10.1029/2006JD007375, 2007.
Dee, D. P., Uppala, S. M., Simmons, A. J., Berrisford, P., Poli, P., Kobayashi, S., Andrae, U., Balmaseda, M. A., Balsamo, G., Bauer, P., Bechtold, P., Beljaars, A. C. M., van de Berg, L., Bidlot, J., Bormann, N., Delsol, C., Dragani, R., Fuentes, M., Geer, A. J., Haimberger, L., Healy, S. B., Hersbach, H., Hólm, E. V., Isaksen, L., Kållberg, P., Köhler, M., Matricardi, M., McNally, A. P., MongeSanz, B. M., Morcrette, J.J., Park, B.K., Peubey, C., de Rosnay, P., Tavolato, C., Thépaut, J.N., and Vitart, F.: The ERAInterim reanalysis: Configuration and performance of the data assimilation system, Q. J. Roy. Meteor. Soc., 137, 553–597, 2011.
Dlugokencky, E. J., Nisbet, E. G., Fisher, R. and Lowry, D.: Global atmospheric methane: Budget, changes and dangers, Philos. T. Roy. Soc. A, 369, 2058–2072, https://doi.org/10.1098/rsta.2010.0341, 2011.
Dlugokencky, E. J., Crotwell, A. M., Mund, J. W., Crotwell, M. J., and Thoning, K. W.: Atmospheric Methane Dry Air Mole Fractions from the NOAA GML Carbon Cycle Cooperative Global Air Sampling Network, Version: 202007, NOAA Global Monitoring Laboratory Data Repository [data set], https://doi.org/10.15138/VNCZM766, 2020.
Fisher, M. and Courtier, P.: Estimating the covariance matrices of analysis and forecast error in variational data assimilation, in: ECMWF Technical Memorandum 220, ECMWF, Reading, UK, https://doi.org/10.21957/1dxrasjit, 1995.
Gurney, K., Law, R., Denning, A., Rayner, P., Baker, D., Bousquet, P., Bruhwiler, L., Chen, H., Ciais, P., Fan, S., Fung, I., Gloor, M., Heimann, M., Higuchi, K., John, J., Maki, T., Maksyutov, S., Ken, M., Peylin, P., Prather, M., Pak, B. C., Randerson, J., Sarmiento, J., Taguchi, S., Takahashi, T., and Yuen, C.W.: Towards robust regional estimates of CO_{2} sources and sinks using atmospheric transport models, Nature, 415, 626–630, https://doi.org/10.1038/415626a, 2002.
Houweling, S., Krol, M., Bergamaschi, P., Frankenberg, C., Dlugokencky, E. J., Morino, I., Notholt, J., Sherlock, V., Wunch, D., Beck, V., Gerbig, C., Chen, H., Kort, E. A., Röckmann, T., and Aben, I.: A multiyear methane inversion using SCIAMACHY, accounting for systematic errors using TCCON measurements, Atmos. Chem. Phys., 14, 3991–4012, https://doi.org/10.5194/acp1439912014, 2014.
Krol, M., Houweling, S., Bregman, B., van den Broek, M., Segers, A., van Velthoven, P., Peters, W., Dentener, F., and Bergamaschi, P.: The twoway nested global chemistrytransport zoom model TM5: algorithm and applications, Atmos. Chem. Phys., 5, 417–432, https://doi.org/10.5194/acp54172005, 2005.
Krol, M., Sager, P., and Segers, A.: 4DVAR data assimilation system using TM5, https://sourceforge.net/projects/tm5/ (last access: 1 June 2022), 2013.
Maasakkers, J. D., Jacob, D. J., Sulprizio, M. P., Scarpelli, T. R., Nesser, H., Sheng, J.X., Zhang, Y., Hersher, M., Bloom, A. A., Bowman, K. W., Worden, J. R., JanssensMaenhout, G., and Parker, R. J.: Global distribution of methane emissions, emission trends, and OH concentrations and trends inferred from an inversion of GOSAT satellite data for 2010–2015, Atmos. Chem. Phys., 19, 7859–7881, https://doi.org/10.5194/acp1978592019, 2019.
Meirink, J. F., Bergamaschi, P., and Krol, M. C.: Fourdimensional variational data assimilation for inverse modelling of atmospheric methane emissions: method and comparison with synthesis inversion, Atmos. Chem. Phys., 8, 6341–6353, https://doi.org/10.5194/acp863412008, 2008.
Myhre, G., Shindell, D., Bréon, F.M., Collins, W., Fuglestvedt, J., Huang, J., Koch, D., Lamarque, J.F., Lee, D., Mendoza, B., Nakajima, T., Robock, A., Stephens, G., Takemura, T., and Zhang, H.: Anthropogenic and Natural Radiative Forcing, in: Climate Change 2013: The Physical Science Basis, Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, Cambridge University Press, Cambridge, UK and New York, NY, USA, https://www.ipcc.ch/report/ar5/wg1/ last access: 1 June 2022), 2013.
Naus, S., Montzka, S. A., Pandey, S., Basu, S., Dlugokencky, E. J., and Krol, M.: Constraints and biases in a tropospheric twobox model of OH, Atmos. Chem. Phys., 19, 407–424, https://doi.org/10.5194/acp194072019, 2019.
Naus, S., Montzka, S. A., Patra, P. K., and Krol, M. C.: A threedimensionalmodel inversion of methyl chloroform to constrain the atmospheric oxidative capacity, Atmos. Chem. Phys., 21, 4809–4824, https://doi.org/10.5194/acp2148092021, 2021.
Pandey, S., Houweling, S., Krol, M., Aben, I., Chevallier, F., Dlugokencky, E. J., Gatti, L. V., Gloor, E., Miller, J. B., Detmers, R., Machida, T., and Röckmann, T.: Inverse modeling of GOSATretrieved ratios of total column CH_{4} and CO_{2} for 2009 and 2010, Atmos. Chem. Phys., 16, 5043–5062, https://doi.org/10.5194/acp1650432016, 2016.
Pandey, S., Houweling, S., Krol, M., and Aben, I.: Influence of Atmospheric Transport on Estimates of Variability in the Global Methane Burden Geophysical Research Letters, Geophys. Res. Lett., 46, 2302–2311, https://doi.org/10.1029/2018GL081092, 2019.
Pandey, S., Houweling, S., and Segers, A.: TM54DVARPP Version 1.0beta.1, Zenodo [code], https://doi.org/10.5281/zenodo.6326373, 2022.
Peters, W., Miller, J. B., Whitaker, J., Denning, A. S., Hirsch, A., Krol, M. C., Zupanski, D., Bruhwiler, L., and Tans, P. P.: An ensemble data assimilation system to estimate CO_{2} surface fluxes from atmospheric trace gas observations, J. Geophys. Res., 110, D24304, https://doi.org/10.1029/2005JD006157, 2005.
Saunois, M., Stavert, A. R., Poulter, B., Bousquet, P., Canadell, J. G., Jackson, R. B., Raymond, P. A., Dlugokencky, E. J., Houweling, S., Patra, P. K., Ciais, P., Arora, V. K., Bastviken, D., Bergamaschi, P., Blake, D. R., Brailsford, G., Bruhwiler, L., Carlson, K. M., Carrol, M., Castaldi, S., Chandra, N., Crevoisier, C., Crill, P. M., Covey, K., Curry, C. L., Etiope, G., Frankenberg, C., Gedney, N., Hegglin, M. I., HöglundIsaksson, L., Hugelius, G., Ishizawa, M., Ito, A., JanssensMaenhout, G., Jensen, K. M., Joos, F., Kleinen, T., Krummel, P. B., Langenfelds, R. L., Laruelle, G. G., Liu, L., Machida, T., Maksyutov, S., McDonald, K. C., McNorton, J., Miller, P. A., Melton, J. R., Morino, I., Müller, J., MurguiaFlores, F., Naik, V., Niwa, Y., Noce, S., O'Doherty, S., Parker, R. J., Peng, C., Peng, S., Peters, G. P., Prigent, C., Prinn, R., Ramonet, M., Regnier, P., Riley, W. J., Rosentreter, J. A., Segers, A., Simpson, I. J., Shi, H., Smith, S. J., Steele, L. P., Thornton, B. F., Tian, H., Tohjima, Y., Tubiello, F. N., Tsuruta, A., Viovy, N., Voulgarakis, A., Weber, T. S., van Weele, M., van der Werf, G. R., Weiss, R. F., Worthy, D., Wunch, D., Yin, Y., Yoshida, Y., Zhang, W., Zhang, Z., Zhao, Y., Zheng, B., Zhu, Q., Zhu, Q., and Zhuang, Q.: The Global Methane Budget 2000–2017, Earth Syst. Sci. Data, 12, 1561–1623, https://doi.org/10.5194/essd1215612020, 2020.
Segers, A. and Houweling, S.: Validation of the CH_{4} surface flux inversion – reanalysis 1990–2017, 1–40, https://atmosphere.copernicus.eu/sites/default/files/202102/CAMS73_2018SC2_D73.2.4.12020_202012_validation_CH4_19902019_v2.pdf (last access: 27 March 2021), 2020.
Williams, J. E., Boersma, K. F., Le Sager, P., and Verstraeten, W. W.: The highresolution version of TM5MP for optimized satellite retrievals: description and validation, Geosci. Model Dev., 10, 721–750, https://doi.org/10.5194/gmd107212017, 2017.
Zhang, Y., Jacob, D. J., Maasakkers, J. D., Sulprizio, M. P., Sheng, J.X., Gautam, R., and Worden, J.: Monitoring global tropospheric OH concentrations using satellite observations of atmospheric methane, Atmos. Chem. Phys., 18, 15959–15973, https://doi.org/10.5194/acp18159592018, 2018.
Zhang, Y., Jacob, D. J., Lu, X., Maasakkers, J. D., Scarpelli, T. R., Sheng, J.X., Shen, L., Qu, Z., Sulprizio, M. P., Chang, J., Bloom, A. A., Ma, S., Worden, J., Parker, R. J., and Boesch, H.: Attribution of the accelerating increase in atmospheric methane during 2010–2018 by inverse analysis of GOSAT observations, Atmos. Chem. Phys., 21, 3643–3666, https://doi.org/10.5194/acp2136432021, 2021.