FootNet v1.0: development of a machine learning emulator of atmospheric transport

There has been a proliferation of dense observing systems to monitor greenhouse gas (GHG) concentrations over the past decade. Estimating emissions with these observations is often done using an atmospheric transport model to characterize the source–receptor relationship, which is commonly termed the measurement “footprint”. Computing and storing footprints using full-physics models is becoming expensive due to the requirement to simulate atmospheric transport at high resolution. We present the development of FootNet, a deep-learning emulator of footprints at the kilometer scale. We train and evaluate the emulator using footprints simulated with a Lagrangian particle dispersion model (LPDM). FootNet predicts the magnitudes and extents of footprints in near real time with high fidelity. We identify the relative importance of input variables of FootNet for improving the interpretability of the model. Surface winds and a precomputed Gaussian plume from the receptor are identified as the most important variables for footprint emulation. The FootNet emulator developed here may help address the computational bottleneck of flux inversions using dense observations.

1 Introduction

Monitoring anthropogenic greenhouse gas (GHG) emissions is important for ensuring the success of the Paris Agreement's long-term goal of mitigating climate change (IPCC, 2022). To that end, there has been a proliferation of dense observing systems over the past decade to better track GHG emissions. Substantial efforts have been made to expand observation networks to better quantify urban GHG emissions, as the majority of the world's population lives in urban areas and the degree of urbanization is projected to increase in the future (United Nations Publications, 2019). For example, the Northeast Corridor GHG observation network was established to quantify emissions of carbon dioxide and methane using tower-based in situ measurements in urban regions in the northeastern United States (Karion et al., 2020). The BErkeley Atmospheric CO₂ Observation Network (BEACO₂N; Shusterman et al., 2016) utilizes low-cost sensors to increase the spatial density of measurements, which could be used to estimate urban emissions on intra-city scales in the San Francisco (SF) Bay Area. The proliferation of urban GHG observation networks allows for decadal analyses of GHG emissions and provides information for improving the efficiency of GHG reduction policies (Mitchell et al., 2018; Lauvaux et al., 2020). There has been a coincident expansion in spaceborne GHG monitoring instruments, which provide similarly dense observations, such as NASA's Orbiting Carbon Observatory-2 (OCO-2) and OCO-3, the TROPOspheric Monitoring Instrument (TROPOMI) on board the Copernicus Sentinel-5 Precursor (S5P) satellite (Veefkind et al., 2012), MethaneSat for methane, and a planned constellation of GHG monitoring satellites (e.g., GOSAT-GW).

The increased number of observational datasets places more constraints on estimating GHG emissions. However, current methods do not scale well with the increasing number of observations. Inferring GHG emissions using atmospheric observations is conventionally done via atmospheric flux inversions (e.g., Jiang et al., 2017; White et al., 2019; Turner et al., 2020). The state of the art in atmospheric flux inversions relies on Eulerian models or Lagrangian particle dispersion models (LPDMs) to simulate atmospheric transport, which provides a means of relating observations to surface fluxes. For example, the four-dimensional variational (4DVar) method uses the adjoint of Eulerian models to calculate sensitivities of GHG concentrations to surface fluxes (Baker et al., 2006; Henze et al., 2007; Jiang et al., 2017; Qu et al., 2022). Kalman filters are also widely used in flux inversions, which calculate covariance matrices between prior fluxes and GHG concentrations simulated by Eulerian models in order to estimate posterior fluxes (Feng et al., 2009; Kang et al., 2011; Miyazaki et al., 2017, 2020). Alternatively, LPDMs can be used to calculate the sensitivity of each observation to its upwind sources by simulating the trajectories of an ensemble of particles advected backward in time (Lin et al., 2004; Fasoli et al., 2018; Jones et al., 2007b; Pisso et al., 2019). The sensitivity of each receptor to its upwind sources, termed the receptor's “footprint”, can then be used to estimate fluxes inversely (e.g., Stohl et al., 2003, 2009; Jones et al., 2007a; Lin et al., 2004, 2021; Stein et al., 2015; Turner et al., 2020). These methods based on full-physics models are becoming prohibitively expensive due to the large computational burden of running high-resolution atmospheric transport models for dense observing systems. The 4DVar method runs the forward and adjoint models iteratively in order to optimize the a posteriori emission, which is hard to parallelize. Kalman filters could benefit from parallelism. However, they still require the forward model, and the computational cost scales up with the number of processors used (e.g., Houtekamer and Mitchell, 2001).

Here we present a machine-learning-based emulator, FootNet, to efficiently calculate footprints of ground-based receptors with high fidelity at kilometer-scale spatial resolution. The footprint emulator reduces the computational and storage cost of Lagrangian model flux inversion systems by 2–3 orders of magnitude, which will better accommodate the increased number of GHG observations. We show the evaluation of the performance of FootNet using independent datasets. Finally, we assess the relative importance of the input variables of FootNet using the permute-and-predict (PaP) method.

2 Construction of the FootNet emulator

Training of the FootNet model is a supervised learning process, which requires ground truth to guide the optimization of the model parameters. Here, we use a full-physics model to generate the ground truth. We simulate footprints using the Stochastic Time-Inverted Lagrangian Transport (STILT) model (Lin et al., 2003; Fasoli et al., 2018), a Lagrangian particle dispersion model. STILT simulations are conducted for two regions: the Barnett Shale region in Texas and the SF Bay Area in California (see Fig. 2). These two regions are chosen because one has a simple topography (the Barnett Shale) and the other is topographically complex (the SF Bay Area). As such, these regions represent limiting cases for the construction and evaluation of the emulator. Further, the combination of the two regions will help prevent overfitting of the model to a single location. For the SF Bay Area, STILT simulations are run from 2018 to 2020 with receptors located at realistic sites deployed in BEACO₂N (see http://beacon.berkeley.edu, last access: 6 March 2025 and Shusterman et al., 2016). Footprints for the Barnett Shale region are generated from a 1-week WRF-STILT simulation in 2013 (Turner et al., 2018). All STILT runs are conducted within 400 × 400 km² domains at 1 × 1 km² spatial resolution (see Fig. 1). The footprints are integrated 72 h backwards from the measurement time because of the 400 km × 400 km domain used by the FootNet model. The time integration period could change depending on the spatial scales and timescales of the inversion systems.

Figure 1 The top row shows a schematic diagram of the FootNet model. A detailed structure of FootNet is shown at the bottom. The orange boxes indicate 3 × 3 convolutional layers. The red boxes represent 2 × 2 max-pooling layers. The light-blue boxes are 2 × 2 transposed convolutional layers. The dark-blue boxes represent the latent vectors concatenated from previous layers (shown as parallel arrows on top).

The output of FootNet is a source–receptor relationship (i.e., footprint H), which represents the sensitivity of atmospheric concentrations at a receptor site to emissions upwind of the receptor. This relationship between the measured concentrations and the emissions in the upwind area can be formulated as

$$\begin{array}{}\text{(1)}& \mathit{y}=\mathbf{H}\mathit{x}+\mathit{b},\end{array}$$

where y is the measured concentration, x is the emission flux in a domain around the measurement location, and b is the background concentration upwind of the domain. The units of y and x can be expressed as dry air mixing ratio (ppb) and flux rates (nmol m⁻² s⁻¹), respectively (Lin et al., 2003). The source–receptor relationships, $\mathbf{H}=\partial \mathit{y}/\partial \mathit{x}$, therefore have a combined unit (ppb (nmol m⁻² s⁻¹)⁻¹).

The calculation of measurement footprints is independent of the observed gas concentrations and could be performed using meteorological variables only. As shown in Table 1, we use four physical parameters from the NOAA High-Resolution Rapid Refresh (HRRR; Benjamin et al., 2016) model as the input variables, including the 10 m zonal wind speed (U10M), 10 m meridional wind speed (V10M), planetary boundary layer height (PBLH), and surface pressure (PRSS). The FootNet model receives input variables at the measurement time (t₀) and 6 h before the measurement time (t₀ − 6 h) to predict footprints at t₀. The choice of 6 h backwards was determined by a series of sensitivity tests of the amount of history information in the input data (see Sect. S1 in the Supplement). We found that including history information from more than 6 h could not improve the performance of FootNet further in the emulation (see Figs. S1–S3 in the Supplement). However, we note that the results from the sensitivity tests could depend on the spatial and temporal scales and the resolutions of the specific inversion problems. Evaluation of the necessary history information in other spatiotemporal regions is warranted.

Table 1 Information about the input variables of FootNet.

We scale the input variables to a similar magnitude for the stabilization of the training process (see Table 1). The output of FootNet is measurement footprints and is transformed by the natural logarithm function to reduce the skewness of the distribution of footprint values. The transformed footprints are filtered to remove values smaller than −20 and are then shifted by +20, corresponding to a scaling of the raw footprints by e²⁰. We find that including Gaussian plumes (see Fig. 1) in the input variables could significantly improve the performance of FootNet. The Gaussian plumes are calculated using the Gaussian plume model (e.g., Stern, 1976; Dobbins, 1979; Zannetti, 1990) with reversed wind fields starting from the measurement site that are used as the initial guess of the upwind areas and the measurement footprints. The Gaussian plumes can be calculated efficiently as a Hadamard product from the inputs listed above and, as such, add minimal computational expense. The Gaussian plumes also provide a localization for FootNet in that they contain the information on the measurement location and provide an initial guess for the spatial structure of the footprint. The FootNet model is trained to learn the nonlinear transformation from the idealized Gaussian plumes and measurement footprints using the meteorological fields. The input variables are interpolated to the 400 × 400 km² domain and the 1 × 1 km² spatial resolution of footprints.

The model structure underlying the footprint emulator is the U-Net model (Ronneberger et al., 2015), which is now broadly applied in the field of Earth science (Ghorbanzadeh et al., 2021; He et al., 2022a, b; Zemskova et al., 2022; Tucker et al., 2023; He et al., 2024). A schematic diagram of the model architecture is shown in Fig. 1. The model consists of four convolutional blocks and four up-convolutional blocks. Each convolutional block is a sequence of two convolutional layers with 3×3 kernels and one 2×2 max-pooling layer. In each convolutional layer, the input images will use the convolution calculation with 3×3 kernels that will scan all the images to generate output images. In max-pooling layers, the input images will be downsampled by taking maximum values in each 2×2 region in the images. Similarly, each up-convolutional layer has one 2×2 up-convolutional layer, followed by two 3×3 convolutional layers. Up-convolutional layers perform transposed convolution operations with 2×2 kernels scanning input images. The outputs from convolutional layers are all transformed by the rectified linear unit (ReLU) function to increase nonlinearity in predictions. In the training process, the entries of 3×3 convolutional kernels and 2×2 up-convolutional kernels will be optimized along the partial gradients of a loss function that measures the difference between the truth and FootNet predictions. More details about deep-learning architectures can be found in Goodfellow et al. (2016).

Figure 2 Locations of receptors for simulations of measurement footprints using the STILT model. Receptors in the SF Bay Area are located at sites in BEACO₂N (Shusterman et al., 2016). Receptors in the Barnett Shale region are at locations used in Turner et al. (2018). The map tiles are from © Stamen Design under a Creative Commons Attribution (CC BY 3.0) license.

We train and evaluate the emulator using a dataset with 10 000 natural log-transformed footprints (logH) from the Barnett Shale and 10 000 footprints from the SF Bay Area as the truth. We apply natural logarithm transformation to the measurement footprints because their values are often highly skewed, which could be challenging for the FootNet model to learn during the training process. The combined dataset is randomly split into 85 % as the training dataset and 15 % as the test dataset. The test dataset is independent of the training process; 15 % of the training dataset is used as a validation dataset during the training process to prevent overfitting. We use mean squared error as the loss function and the Adam optimization algorithm.

We use intersection over union (IoU) to measure the accuracy of the area of footprints predicted by FootNet, which is defined as follows:

$$\begin{array}{}\text{(2)}& \text{IoU}={\displaystyle \frac{|Y\cap \widehat{Y}|}{|Y\cup \widehat{Y}|}}.\end{array}$$

Here, Y and $\widehat{Y}$ stand for the ground truth (footprints simulated by STILT) and the FootNet predictions, respectively. The absolute value bars ($|\cdot |$) here refer to the area of a region. Specifically, the intersection, $|Y\cap \widehat{Y}|$, calculates the area of the region where both the truth and FootNet predictions show nonzero footprints. Similarly, the union, $|Y\cup \widehat{Y}|$, represents the area of the region where either the truth or FootNet predictions show nonzero footprints. IoU is widely used to evaluate the ability of deep-learning models to make accurately localized predictions. We also compute Pearson correlation coefficients (r) for footprints in the intersection areas between the truth, as simulated by STILT, and the corresponding FootNet predictions to help assess the performance.

Ultimately, we are interested in better understanding what drives the predictions of the FootNet model. As such, we use the PaP method to calculate the importance of input variables for footprint emulation, which provides some interpretability of the FootNet model (Fisher et al., 2019). The PaP method estimates variable importance by permuting each input variable with different data samples, and the subsequent performance change represents FootNet's sensitivity to the permuted variable. We estimate variable importance by calculating performance changes in the correlation, the IoU, and the root mean square error (RMSE) of the predicted footprints.

3 Evaluating the performance of the FootNet emulator

Figure 3 demonstrates the evolution of FootNet predictions during the training process and the overall performance of FootNet after the training converges. Figure 3d shows a footprint simulated by the STILT model from the test dataset, where the footprint is highly nonlinear with a change in direction near the receptor. The corresponding FootNet predictions are shown in Fig. 3a–c. After iteration A (shortly after the training starts), FootNet predicts measurement footprints around the receptor with a large negative bias and a low correlation coefficient of 0.49. Iteration B is about halfway in the training process, after which the FootNet prediction captures the general shape of the footprint better and the correlation is improved to 0.61. The training stops after iteration C. The final FootNet prediction has enriched details and attains a correlation coefficient of 0.75. Compared to the truth in Fig. 3d, the IoU of FootNet predictions improves from 0.28 after iteration A to 0.51 after iteration B and attains a final IoU of 0.76 (see Fig. 3e). Figure 3f shows the comparison between the truth and FootNet predictions for all footprints in the test dataset. FootNet predictions show a slight negative bias compared to footprints simulated using the full-physics STILT model. The overall Pearson correlation coefficient (r) between FootNet predictions and STILT simulations is 0.58. We conclude that FootNet is able to emulate the source–receptor relationship under both simple (Barnett Shale, TX) and complex (SF Bay Area, CA) meteorological conditions with high fidelity. However, it is worth mentioning here that we find some performance degradation using an alternative splitting of the data based on different time periods. Because the training dataset used to construct version 1 of FootNet has a relatively small size, similarities between samples are hard to avoid fully by randomly selecting training data samples, which could lead to generalizability issues when using FootNet version 1 over regions and time periods too different from the training dataset. This generalizability issue could be largely mitigated by increasing the volume of the training dataset in the future (Dadheech et al., 2024).

Figure 3 Convergence of the training process and evaluation of the model performance on the independent test dataset. Footprints are transformed using the natural logarithm (ppb (nmol m⁻² s⁻¹)⁻¹ before the transformation). (a–c) FootNet predictions from three stages in the training process, corresponding to the truth in panel (d). The blue arrows represent the wind vectors, and the green stars show the locations of the receptors. (e) Comparison between footprints simulated by STILT and FootNet predictions in panels (a)–(c). (f) Two-dimensional histogram of all natural log-transformed footprint (logH) values simulated by STILT and the corresponding predictions made by FootNet from the test dataset.

We then evaluate the performance of FootNet in predicting individual footprints for the two regions. Figure 4 shows footprints from STILT and FootNet for the two regions: the Barnett Shale and the SF Bay Area. Figure 4e shows results from the simple case (Barnett Shale, TX), where the footprint is similar to an idealized Gaussian plume with time-reversed winds. FootNet captures both the magnitudes and spatial patterns of the footprint well, with an IoU of 0.73 and a correlation coefficient of 0.54. Figure 4b and f demonstrate a more complicated meteorological scenario in the Barnett Shale region. The IoU metric and correlation coefficient between the STILT footprint and the FootNet prediction are 0.71 and 0.61, respectively, for this more complex scenario.

Figure 4 Evaluation of individual FootNet predictions from the test dataset. Footprints are transformed using the natural logarithm (ppb (nmol m⁻² s⁻¹)⁻¹ before the transformation). (a–d) Footprints simulated by the full-physics STILT model for the Barnett Shale region and the SF Bay Area. (e–h) Footprint predictions made by FootNet corresponding to panels (a)–(d). The blue arrows represent the wind vectors, and the green stars show the locations of the receptors. (i–l) Comparison and correlation between the truth and predictions for the four examples.

Atmospheric transport in the SF Bay Area is decisively more complex because the region includes steep topography, air–sea interactions, and numerous valleys and deltas. Figure 4c and d show results from the full-physics model for the SF Bay Area. Emulation of footprints in the Bay Area is more challenging and has an overall degraded fidelity as compared to the Barnett Shale region. Figure 4c and g show a receptor with the bulk of the footprint in the northwestern quadrant of the domain as a result of the typical summertime meteorology of the SF Bay Area, with westerly flow bringing air masses into the SF Bay Area past the Golden Gate Bridge. The shape and magnitude of the footprint are predicted by FootNet with an IoU of 0.53 and a correlation coefficient of 0.83. Figure 4d and h show a more complex meteorological scenario where the FootNet prediction has an IoU of 0.56 and a correlation coefficient of 0.78 as compared to STILT.

There have been other methods developed to improve the efficiency of footprint calculations. For example, Roten et al. (2021) used nonlinear weighted averaging to interpolate footprints from locations near the receptors. Fillola et al. (2023) developed a similar footprint emulator based on gradient-boosted regression trees (GBRTs) at a coarse spatial resolution (20–30 km at the middle latitudes) and 10 grid cells around the measurement location. Compared to previous work, the FootNet model reproduces the full-physics model with high fidelity at high resolution. This is remarkable given that the complex topography and meteorology of the regions studied here could complicate transport at the kilometer scale and the emulation of footprints. Moreover, FootNet only takes meteorological fields and the idealized Gaussian plume as its inputs. No additional LPDM simulations are needed to generate footprint predictions after the training process.

Emulation of footprints using the FootNet model brings co-benefits for computational efficiency and storage cost and better facilitates the application of LPDM-based flux inversion systems to dense observing systems. To conduct kilometer-scale emission inversions using 1 d of observations made at the 40 BEACO₂N sites in the SF Bay Area (approximately 650 observations per day), it takes the full-physics STILT model about 640 core hours to generate the required footprints. The generation of each footprint prediction takes ∼ 1 s on a 32-core compute node, which can be reduced further to 0.08 s on an NVIDIA A2 graphics processing unit (GPU). Only 6 min are required for FootNet on an A2 GPU node to generate the footprints for 1 d of BEACO₂N measurements. The storage requirement also makes it impractical to use full-physics models in high-resolution flux inversions with dense observations. Hourly footprints for 1 week of BEACO₂N measurements would require 4 TB storage space for future reuse. With FootNet, footprints could be generated in near real time, and there is no need to store the computed footprints.

Figure 5 shows the ranking of variable importance for FootNet calculated using the PaP method on 1000 randomly selected data samples. Overall, the most important meteorological variables are the 10 m wind speeds, which lead to a 0.2–0.3 decrease in correlation, and the IoU drops by 0.1–0.2 after being permuted. Permuting Gaussian plumes degrades the correlation and IoU of FootNet predictions by 0.1 and 0.03, respectively. We find less sensitivity of FootNet predictions to surface pressure and planetary boundary layer height than other input variables. This is because we only have training data from two locations in the current version of the model, and these two meteorological fields show much lower variability than the wind fields in the training dataset. We still include surface pressure and PBLH as input variables because they are essential information for the generation of measurement footprints. We expect to see greater importance for surface pressure and PBLH in the future for a general version of FootNet trained using footprints from more locations. Figure 5 also shows that input variables from t₀ − 6 h have consistently greater importance than t₀.

Figure 5 Rankings of variable importance estimated using the permute-and-predict (PaP) method on 1000 data samples. (a-c) Variable importance shown as a drop in correlation, a drop in the IoU, and an increase in the RMSE after permuting the 10 input variables. Orange lines show the medians. Boxes indicate ranges from the first quartiles to the third quartiles. Whiskers are the 1.5 interquartile ranges (IQRs) from the boxes. Circles are outliers.

The PaP method only provides a rough estimate of variable importance, and the intercorrelation between input variables can lead to an inflation of the feature importance (Hooker et al., 2021). Nevertheless, the estimated variable importance for FootNet is in alignment with our understanding of the calculation of footprints in a full-physics model, which relies on the advection of particles driven by precomputed wind fields. The Gaussian plume is also identified as highly important, because it is the only input field providing information about the locations of receptors.

4 Conclusions

We described the development of a machine-learning-based emulator of surface measurement footprints, FootNet. The footprint emulator can be used to improve computational efficiency when estimating high-resolution GHG fluxes using measurements made by dense observing systems. The FootNet model was trained and evaluated using footprints simulated by the STILT full-physics model for the SF Bay Area and the Barnett Shale region. We showed the convergence of FootNet predictions with the STILT truth as the training iterates. The overall correlation between FootNet predictions and the STILT truth in the test dataset reaches 0.58 after full convergence. The emulator predicts both the extents and magnitudes of footprints well, with high fidelity. We estimated the importance of input variables using the PaP method to improve the interpretability of the FootNet model. We found that 10 m wind speeds and Gaussian plumes have the greatest importance for the emulation of footprints. Emulation of footprints using FootNet brings co-benefits for computational efficiency and reduced storage cost, which will make it feasible to deliver high-resolution estimates of GHG fluxes in near real time using proliferated dense observing systems in the future.

Due to the computational cost required by the generation of high-resolution footprints, we only included footprints generated from previous studies for the two locations in training version 1.0 of FootNet. We are actively generating new footprints at 1 km from a broader region in order to further improve the emulator's performance, especially in regions with meteorological conditions different from the two locations used in this study (Dadheech et al., 2024). Generalizing this source–receptor emulator to other regions is being tackled in the next version of FootNet.