Spatial attribution of aircraft mass balance experiment CO₂ estimations for policy-relevant boundaries: New York City

A high proportion of the U.S. population (82.3%) resides within Census-defined urban regions (US Census Bureau, 2022), projected to reach 89.2% by 2050 (United Nations, 2019). A recent study has shown that ≈44% of carbon dioxide (CO₂) emissions in the United States originate from the urban sector for 2015 (Gurney et al., 2020). U.S. cities are currently employing ambitious climate mitigation strategies (Sharp et al., 2011; Sethi et al., 2020). For example, the New York City (NYC) Council has passed the Climate Mobilization Act, which is a series of legislations with the aim to reduce greenhouse gas (GHG) emissions and includes periodic energy audits for large commercial buildings (New York City Council, 2019). As a result, it is necessary for owners to consider potential carbon emissions when reviewing real estate capital projects.

Assessing the progress of urban emission mitigation policies requires accurate monitoring of GHG emissions. Cities evaluate their policies by following protocols utilizing self-reported emission or “bottom-up” inventories drawing on a range of fuel consumption and activity data. This approach begins with emission estimates from statistically analyzed activity data together with emission factors along politically defined boundaries. The uncertainties of such an approach alone are difficult to quantify and often result in inconsistent reporting when compared to peer-reviewed works (Turnbull et al., 2015; Quilcaille et al., 2018; Lauvaux et al., 2020; Gurney et al., 2021). Urban emission determinations using top-down approaches (based on atmospheric mole fraction measurements) are particularly challenging due to the environmental and meteorological complexities surrounding cities, in addition to the large spatial distribution of emissions stemming from multiple sources and sinks (Pataki et al., 2006). For this reason, the Northeast Corridor (i.e., northeast metropolis area from Boston to Washington, DC) is a region of interest to act as a test bed for the assessment of urban emissions among major cities including Boston (McKain et al., 2015; Sargent et al., 2018), Baltimore/Washington, DC (Lopez-Coto et al., 2017; Mueller et al., 2018; Ren et al., 2018; Ahn et al., 2020; Karion et al., 2020; Lopez-Coto et al., 2020; Lopez-Coto et al., 2022), Philadelphia (Anderson et al., 2021), and the New York-Newark urban region (Plant et al., 2019; Floerchinger et al., 2021; Hajny et al., 2022; Pitt et al., 2022).

Assessments must be performed through measurable and verifiable means of GHG emission quantification. Significant advancements in space-based retrievals of GHG signatures have been achieved through the development of high-resolution measurement satellites providing fine spatial column averages with which to conduct long-term monitoring of emission trends for megacities (Demetillo et al., 2020; Kiel et al., 2021; Lauvaux et al., 2021; Plant et al., 2022). However, most satellite instruments only operate in the presence of sunlight and do not have the ability to acquire data in regions with heavy cloud cover and variable surface albedo. Therefore, satellite measurements require complementary techniques to help support space-based observations.

Emission fluxes can be quantified using atmospheric transport models to trace the magnitude and pattern of the emitted species to the observing system (in situ or remotely sensed column averages). This approach allows for the calculation of simulated enhancements within a region of interest given gridded emission map information. Statistical estimators of emissions (i.e., inverse modeling approaches) can then be employed to obtain an improved estimate of the model emission map by comparing against the observed signal. Ground-based measurements from tower networks provide an effective method to conduct long-term monitoring of measured GHG mole fractions (Davis et al., 2017; Verhulst et al., 2017; Huang et al., 2019; Lauvaux et al., 2020).

The aircraft-based mass balance experiment (MBE) approach has also been used to measure and quantify GHG emission from several cities (Cambaliza et al., 2015; Heimburger et al., 2017; Ren et al., 2018; Pitt et al., 2019; Ahn et al., 2020) and point sources (Lavoie et al., 2015; Conley et al., 2016; Ryoo et al., 2019; Hajny et al., 2023). This approach enables large spatial coverage, albeit at the cost of limited temporal information. Typically, the aircraft follows a raster pattern directly downwind of the source plume and perpendicular to the mean wind direction. A series of stacked horizontal transects are performed resulting in a 2D sampling plane downwind of the city. The mass flux through this 2D plane is calculated from the measured enhancements above the background mole fractions. The background mole fraction can be acquired from the upwind sampling of the city or from the downwind measurements along the plume edge (Cambaliza et al., 2014; Heimburger et al., 2017; Mueller et al., 2018; Pitt et al., 2019; Balashov et al., 2020). This method is particularly effective for point sources and isolated cities with negligible outside influence (Karion et al., 2013; Cambaliza et al., 2015; Hajny et al., 2019). However, it is more challenging to assess the background for cities with complex surroundings (Pitt et al., 2019) due to multiple emission sources (e.g., cities along the Northeast Corridor) and, thus, traditional aircraft-based MBE techniques only provide bulk emissions that are difficult to unambiguously attribute to a precisely defined spatial area.

In this study, we developed a new methodology to spatially allocate CO₂ emissions estimated from MBEs to policy-relevant boundaries and applied this approach to NYC. We conducted a series of aircraft measurements downwind of NYC during the nongrowing season and employed the conventional MBE technique to estimate the bulk downwind CO₂ emissions. To accurately represent emissions within the specific area of interest (AOI), we employed an ensemble of CO₂ emission inventories and transport models to calculate the fraction of the MBE emission rate that represents the NYC 5 boroughs (i.e., the Bronx, Brooklyn, Manhattan, Queens, and Staten Island).

Flights were performed with a modified twin-engine Beechcraft Duchess, Purdue University’s Airborne Laboratory for Atmospheric Research (ALAR) (Cambaliza et al., 2014; Heimburger et al., 2017; ALAR, 2022). Details on the design and structural configuration of the aircraft payload can be found in the study of Cambaliza et al. (2014). ALAR is equipped with a global positioning and inertial navigation system (GPS/INS), and a 9-port pressure probe (Best Air Turbulence, BAT) providing 3D winds at 50 Hz (Crawford and Dobosy, 1992; Garman et al., 2006; Garman et al., 2008). A Picarro Cavity Ringdown Spectrometer (CRDS) provided in situ measurements of CO₂, methane (CH₄), and water vapor (H₂O_v) mole fractions at 0.5 Hz; however, this work only investigates CO₂ observations. The Picarro CRDS was calibrated with an in-flight 3-point calibration system using National Oceanic and Atmospheric Administration certified standard cylinders for CO₂ (ref. scale: WMO-CO₂-X2007) (Tans et al., 2017) and CH₄ (ref. scale: WMO-CH₄-X2004A) (Dlugokencky et al., 2005) during each flight mission. Additional details of instrument calibrations are discussed in detail in our recent works (Hajny et al., 2022; Pitt et al., 2022).

Figure 1a shows the compiled flight tracks from all experiments discussed in this work and overlaid on a state map with the AOIs highlighted, as the urban region (i.e., Greater New York-New Jersey area, pink) and NYC (teal). The aircraft was flown perpendicular to the wind direction and performed a series of horizontal transects downwind of NYC at varying altitudes extending to the top of the planetary boundary layer (PBL), as illustrated in Figure 1b. All flight experiments were conducted between 11:00 and 16:00 (local time) allowing for a fully developed and consistent PBL height during the sampling period. The PBL height was identified from the vertical profiles of trace gas species, H₂O_v, potential temperature, and variance in the vertical winds. Two sets of vertical profiles were flown during the upwind and downwind curtains whereby the aircraft follows an ascending/descending spiral maneuver from the near-surface to ≈3,000 m above ground level (a.g.l) into the free troposphere. Upwind transects were flown during most flights prior to performing the downwind passes, allowing for the identification of the city plume due to elevated GHG mole fractions. In addition, the downwind horizontal transects were sampled beyond the horizontal extent of the observed plume allowing for simultaneous measurement of regional background mole fractions along the transect edges and outside the influence of the direct urban emissions. The 7 MBE flights downwind of NYC analyzed in this work were conducted between November 2018 and March 2020 (all pre-COVID impacts) with all flights occurring in the months of November, February, or March to minimize the effects of CO₂ exchange with the biosphere.

The aircraft-based MBE approach was used to estimate the total CO₂ flux through a perpendicular plane downwind of NYC. This approach assumes a static PBL with consistent emissions over the time frame of the experiments, no horizontal flow divergence/convergence, and a constant wind direction throughout the sampling period (Tadić et al., 2017; Ryoo et al., 2019). The net GHG mass flow (in kmol/s) along this crosswind plane is calculated following equation:

where [0, z] is the PBL depth, −x and +x are the effective horizontal boundaries within the plume, $U⊥,ij$ is the perpendicular component of the horizontal wind speed (in m/s) at the downwind location (x, z), P_ij is the measured pressure (in atm), T_ij is the measured temperature (in K), and R is ideal gas constant (in $atm⋅m3mol⋅K$⁠), C_ij is the instantaneous measurement of the downwind dry air CO₂ mole fraction (in $µmolCO2mol dry air$⁠, or ppm), and lastly, $CijBG$ is the background dry air CO₂ mole fraction (in ppm) at the same location (x, z) (Cambaliza et al., 2014; Heimburger et al., 2017; Ren et al., 2018). Note that i and j correspond to x and z dimension, respectively. Note that the downwind transects were extended beyond the edges of the plume such that the concentrations return to background levels. Certainly, the choice of wind data and background can impact the flux estimations (Cambaliza et al., 2014; Ryoo et al., 2019). For example, Ryoo et al. (2019) have shown that when utilizing the perpendicular wind component for the flux calculation, the sensitivity of the emission estimate to the selection of background was minimal (<1%). In addition, if the substantial background mole fractions of trace gas species like CO₂ are not properly considered, the measurement uncertainty of volume-integrated horizontal wind divergence can result in significant uncertainty in the estimation of flux. However, choosing the minimum observed value as the background on each vertical level resulted in significantly different calculated fluxes at the urban scale (Ryoo et al., 2019). For this work, we used the 10s rolling-averaged perpendicular wind speed to derive the CO₂ flux estimations consistent with our previous works (Cambaliza et al., 2015; Heimburger et al., 2017; Hajny et al., 2019). We also defined the background as a linear fit anchored along the edges of the urban plume for each transect as shown in Figure S1. The bounds of the plume for each downwind transect are defined as the CO₂ mole fractions that were above the lowest 5th percentile mole fraction on the negative and positive distance side, respectively (dashed lines in Figure S1b). The fluxes derived from the CO₂ enhancements (measurements after subtracting the background), and perpendicular winds were then interpolated using a kriging approach to fill in the spatial gaps between actual measurements to generate a 2D gridded plane (Matlab-based EasyKrig3.0) (Chu, 2004) as shown in Figure 1c and d. The linear variogram model used for EasyKrig3.0 was chosen to fit the empirical semivariogram of the data, determined through a visual examination of the experimental variogram. Three key parameters were used to fit the theoretical variogram including the nugget (the variance or variability at distances smaller than the minimum separation distance between data points), sill (the maximum level of variability or variance that is observed as the lag distance approaches infinity), and range (the distance at which the dependence diminishes to a negligible level). The quality of the interpolation was evaluated using the Q1 (deviation distribution for the mean) and Q2 (deviation distribution for the standard deviation) criteria (Chu, 2004). In addition, the kriging resolution of 10 m in the vertical direction and 100 m in the horizontal direction was used for all 7 flights. The calculated bulk emission rate (kmol/s) is then obtained by integrating the net mass flow of the trace species, which is proportional to the total enhancements observed downwind (Equation 1). Additional details on the aircraft-based mass balance kriging approach can be found in the study by Cambaliza et al. (2014). In addition, we performed a sensitivity test using the “gstat: Spatial and Spatio-Temporal Geostatistical Modelling, Prediction and Simulation” (Pebesma, 2004; Gräler et al., 2016) R package as an independent analysis tool, using 4 variogram models “Sph” (Spherical), “Exp” (Exponential), “Mat” (Matern), and “Ste” (Matern, M. Stein’s parameterization) (Stein, 1999). Note that vertical profile measurements along the downwind curtain were removed prior to kriging analysis as it conflicts with the calculation due to overlapping or closely adjacent data points as vertical profile spirals are performed.

Common limitations of the aircraft-based MBE approach are fuel, air traffic control, Federal Aviation Administration regulations (minimum of 1,000 ft above surface structures in urban environments), cloud cover, and terrain restricting the number of downwind transects within the PBL, and consequently, being unable to fully capture the full vertical extent of the downwind enhancements. Figure 1d illustrates the potential underestimation of the total flux through the downwind curtain due to these limitations. Previous works have attempted to represent and include extrapolation of surface emission below the lowest flight level using various near surface concentrations and extrapolations to the surface (Gordon et al., 2015; Ryoo et al., 2019). We employed 3 extrapolation methods (Hajny et al., 2023) to fill in the spatial gaps both near the surface level and the PBL height, as shown in Figure S2. Note that our flights took place under either neutral or unstable conditions, and the aircraft transects were carried out downwind of the primary sources. Generally, these transects spanned over 2 eddy turnovers from the sources, allowing sufficient time for vertical mixing of concentrations, at least until the lowest transect. The first extrapolation method repeats the observed CO₂ measurements along the lowest and highest passes, at the bottom and top of the boundary layer, respectively. In other words, we assumed that the lowest altitude transect is the best mathematical guess to represent ground level while the highest altitude transect represented the mole fractions at the PBL height (herein refer to as synthetic transects) (Figure S2b). The second method addresses these large gaps without measurements by taking the average of the 2 passes at each extreme and assumes the data within those gaps are constant with height (Figure S2c). Lastly, we calculate the average flux across all transects and extrapolate to the surface and PBL height (Figure S2d). Henceforth, the average of the 3 extrapolation methods is used to calculate the emission rate with an average relative standard deviation (σ) of ≈10% (1σ) among the 7 flights, that represents the uncertainty resulting from these data gaps. Note that ≈10% is referring to the mean of the standard deviation across the 7 fights.

In complex emission scenarios such as urban regions, the AOI is surrounded by multiple sources that can obscure the signal, making it difficult to quantify emissions within policy-relevant boundaries. In this case, we wish to estimate emissions from NYC. However, emission sources outside the NYC boundary influence the CO₂ mole fractions measured at both within-plume and background locations, thus impacting the derived MBE emission rate. In some cases, emission sources from within the NYC boundary can also influence the measured background. Therefore, traditional MBE analysis cannot distinguish between emissions within NYC and emission sources outside the AOI.

To better understand the situation and address this issue with the MBE approach, we utilize a transport model (using footprints and emission inventories). Figure 2a shows the fluxes under the modeled footprints of the entire downwind transect (left), as well as the fluxes under the footprints of the plume edges used as background (center and right), for the flight conducted on March 4, 2020. In a situation like this, where the transect is sensitive to sources upwind and downwind of the AOI, the mole fraction observed will reflect the superposition of these sources. On the other hand, the non-zero fluxes at the edges, larger in the southern part in this case, elevate the mole fractions at the background points differently. A linear interpolation between the edges of the transect will provide a background that will account for some of the superposition of sources (approach used in the MBE). However, how this linear interpolation relates to the “true” background (i.e., the mole fraction that would have been measured at each location in the absence of NYC emissions) depends on the relative strength of these 2 competing effects.

To quantify this impact, we define 2 tracers: (1) Total CO₂ (Enh_TOT), which are the modeled enhancements produced by all CO₂ sources in a large regional domain up to ≈800 km upwind of the city, and (2) NYC only (Enh_NYC), which are the modeled enhancements produced by sources only within the AOI, that is, the 5 boroughs of NYC in this case. The left panel in Figure 2b shows (black solid line) the total modeled enhancements (Enh_TOT) for the transect, which reflect the elevated mole fraction of the south edge with respect to the north edge. We can now mimic the background determination used in the MBE by linearly interpolating between the edges (black dashed line) and compared to the “true” background (red dashed line), that is, the enhancements produced by all sources outside of the AOI (Enh_TOT – Enh_NYC). In this particular case, we can see how the linear background is accounting for a large portion of enhancements produced elsewhere (enhancements below the black dashed line). However, it still fails to account for a portion of enhancements produced outside NYC (enhancements above the black dashed line and below the red dashed line). The right panel shows the enhancements after subtracting the linear background (Enh_MBE = Enh_TOT – linearBG, black dashed line) (hereafter referred to as the MBE-like modeled enhancements) and the enhancements due to sources only within NYC (Enh_NYC, red dashed line) (hereafter referred to as NYC enhancements). From the curves, it is clear that the integral enhancement produced by NYC is smaller than the integral enhancement calculated using the linear background approximation used in the MBE for this transect. In fact, we can quantify the fraction (Φ_j) of NYC integral enhancements to the integral MBE-like modeled enhancements for this transect leg (j) as in the following equation:

This fraction (Φ) gives us a measure of the enhancements that we are interested in measuring relative to the enhancements we would infer from MBE by applying the linear background. Note that j corresponds to the ensemble member of inventory (3), transport (8) and transect (29) within the boundary layer. This results in a total of 696 calculated modeled values of Φ derived across the 7 research flights. If the fraction is smaller than 1 (Φ < 1), as in the example shown in Figure 2, the MBE-like modeled enhancements carry contributions from sources inside and outside of the AOI because the linear background is underestimating the contribution from sources outside the AOI. On the other hand, if the fraction is larger than 1 (Φ > 1), the MBE-like modeled enhancements are missing contributions from inside the AOI because the transect edges are still sensitive to the AOI and thus the linear background partially removes contributions from sources inside the AOI (Figure S3). However, in some cases, the fraction can be less than 0 (Φ < 0) due to poorly defined MBE-like modeled enhancements along the transect (i.e., Σ(Enh_TOT – linearBG) < 0). Note that 6 of the 696 calculated modeled fractions were found to have Φ < 0 and were excluded in the ensemble results.

Since in the MBE calculation, the emission rate is proportional to the integral enhancements (emission rate ≈ Σ(ΔC)), the ratio between integral enhancements (i.e., fraction, Φ) is also the ratio between emission rates. Thus, we can relate the mass balance bulk emission rate (MBE_bulk; using the linear background) calculated using the observations with the source-apportioned mass balance emission rate for the city (MBE_C) using the model-derived enhancement fraction (Φ) as shown in the following equation:

Since this approach is dependent on the spatial distribution of sources and the atmospheric transport, we used an ensemble of 3 prior emission inventories and 8 transport models. We used as prior CO₂ emissions the annual mean of ACES v2.0 (Gately and Hutyra, 2017), Vulcan v3.0 (Gurney et al., 2020), and EDGAR v5.0 (Crippa et al., 2019; European Commission, 2019). In order to account for biogenic CO₂ emissions into the modeled enhancements, we used the vegetation photosynthesis and respiration model (VPRM) (Mahadevan et al., 2008; Gourdji et al., 2022). VPRM estimates the net ecosystem exchange, which encompasses both respiration and photosynthesis, by considering various factors such as the enhanced vegetation index, land surface water index, temperature, and carbon uptake from photosynthesis. We use the Gourdji (2021) and Gourdji et al. (2022) version that includes a reformulation of the respiration equation functional form that helps correct a respiration low bias in the original version. Also, model parameters were specifically optimized for the Northeast U.S. Enhancements simulated using these biogenic fluxes were then combined with (added to) the anthropogenic enhancements to calculate the total enhancements. Furthermore, the sensitivity of measurements to surface fluxes (gridded footprints) along the sampling trajectory of the aircraft was calculated using the Hybrid Single Particle Lagrangian Integrated Trajectory (HYSPLIT) model v5.0.0 (Stein et al., 2015; Loughner et al., 2021) which incorporates features from the Stochastic Time-Inverted Lagrangian Transport (STILT) model (Lin et al., 2003; Fasoli et al., 2018). The HYSPLIT gridded footprints were calculated using 4 different input meteorological (MET) models for 2 domains, a coarser, regional domain at 0.8° resolution (bounds 34.4°N 44.4°N, 83.7°W, and 69.7°W) and a higher resolution domain at 0.02° resolution (bounds 39.2°N 42.0°N, 75.7°W, and 72.1°W). The enhancement calculation used both nested domains. The input MET models used for this work include (1) the Global Forecast System (GFS) model, (2) the European Centre for Medium Range Weather Forecasts Fifth Reanalysis (ERA5), (3) the North American Mesoscale Forecast System (NAMS) model, and (4) the High-Resolution Rapid Refresh (HRRR) model. In addition, 2 mixing parameterizations were employed: kblt2 (Kantha and Clayson, 2000) and kblt5 (Hanna, 1982), yielding 8 different transport models to generate the calculated footprints, from which to assess the variability derived from the model uncertainty. Note that by applying our model-derived fractions to the MBE results, we are incorporating transport model information into the MBE_C calculations Additional details of model parameters are shown in Table S1 and discussed in our other works including Pitt et al. (2022) and Hajny et al. (2022).

The observed measurements are influenced by CO₂ emissions from all sources including point sources such as fossil fuel power plants and incinerators; mobile combustion sources within and around NYC; and fossil fuel combustion attributed to commercial and residential buildings. Therefore, the aircraft-based MBE approach estimates the bulk CO₂ emissions rate representing the net emission from all sources along the background-subtracted downwind footprint.

Figure 3 shows the kriged CO₂ flux distribution within the PBL and along the downwind curtain for the 7 research flights conducted on the following days: (1) November 9, 2018, (2) March 1, 2019, (3) March 27, 2019, (4) November 10, 2019, (5) November 15, 2019, (6) February 16, 2020, and (7) March 4, 2020. For visual purposes, data in Figure 3 are plotted to only include positive fluxes. However, the calculated MBE_bulk emission rate includes negative flux values, which can be considered as a component of background noise (especially at city scales) and has been a standard procedure in similar works (Cambaliza et al., 2014; Heimburger et al., 2017). Henceforth, the reported MBE_bulk emission rate is calculated from the average of all 3 kriged extrapolation approaches as discussed in the Methods section. Large spatial gaps between the measurement transects is a limiting factor when performing kriging analysis for aircraft-based MBEs resulting in additional uncertainty in the calculations. An extreme case of this can be observed for March 1, 2019 having an insufficient number of downwind transects to accommodate measurements up to the PBL height. In this case, the original extrapolation to the PBL top as described in the Method section provided an emission rate of 151 kmol/s. Since most of it was derived from the extrapolation, we used the vertical profile measurements to estimate the remaining CO₂ values between the highest transect at 400 m and PBL height at 1,050 m, instead of relying on a well-mixed assumption. This was performed for the March 1, 2019 research flight as shown in Figure S4 to account for the limited number of downwind transects. The values of CO₂ along the synthetic transects were then estimated and scaled based on a fit applied to the vertical profile measurements, as shown in Figure S4b. Prior to the addition of synthetic transects, the total emission rate for March 1, 2019 research flight was calculated to be 52.4 kmol/s compared to 112 kmol/s when estimating for the CO₂ mole fractions within the whole PBL based on the CO₂ vertical profile.

The mean MBE_bulk for CO₂ emissions for the 7 research flights was calculated to be 68 ± 31 (1σ) kmol/s. The significant variability in MBE_bulk (44%) can be attributed to the day-to-day variability of emissions, the sampling footprint, and methodological uncertainties. Cambaliza et al. (2014) examined sources of uncertainty for the aircraft-based MBE approach and found that the largest uncertainty was attributed to the variability in the background mole fraction of CO₂, and the PBL height. In addition, the authors estimated that the interpolation of flux measurements throughout the downwind curtain amounted to between 8% and 12% uncertainty. However, improvements in precision of campaign mean aircraft-based MBE can be achieved through repeated measurements and averaging (Heimburger et al., 2017).

As discussed, the bulk CO₂ emissions were calculated from the interpolation of the CO₂ flux by kriging the downwind measurements using the Matlab-based EasyKrig3.0 program package (Chu, 2004) using a linear variogram model. The sensitivity test using the gstat (Pebesma, 2004; Gräler et al., 2016) based kriging (R-Krig) using the synthetic transects case is illustrated in Figure S5 where the best-performing variogram model (out of the 4 tested) was selected separately for each flight. A detailed description of the variogram modeling scheme is discussed in the study of Pebesma et al. (2004). We then compared the total CO₂ emission rates derived from EasyKrig3.0 to the total emission rate calculated from R-Krig, as shown in Figure S6. Among the 4 variogram models tested in R-Krig, the calculations yielded an average difference of 0.7%. Similarly, the calculated emission rates from R-Krig and the emission rate derived from EasyKrig3.0, which uses a linear variogram model, were different by only 1.2% on average, with a maximum difference of 6.2% for November 10, 2019. Lastly, a linear regression between both sets of results indicated that both methods are statistically indistinguishable (Figure S6b).

Figure 4 shows the distribution of fractions (Φ) for each research flight using the ensemble of transport models and prior fluxes. Note that separate Φ values were calculated for each transect sampling the plume (i.e., the boxplots consist of individual values representing a single transect-model-prior combination). As discussed earlier, Φ of less than 1 implies that the bulk MBE estimate includes sources outside NYC in addition to those within the city. On the other hand, Φ greater than 1 indicates that the MBE_bulk estimate is smaller than the total emissions produced in NYC as the background incorporates NYC emissions. In general, for 5 out of 7 flights the average Φ was less than 1, ranging between 0.52 and 0.93. Research flights November 9, 2018 and March 27, 2019 had a fraction greater than 1, with an average Φ value of 1.2 and 2.3, respectively. Values of Φ >1 are due to the sensitivity of the plume edges to NYC itself. We can then relate this Φ value of modeled enhancement to the ratio of emission rates from MBE calculations to obtain a source-apportioned mass balance CO₂ emission rate for NYC, MBE_C, as described in Equation 3.

Figure 5 breaks down the source-apportioned MBE_C emission rate across the different ensemble prior inventories (top panel), transport models (middle panel), and research flights (bottom panel). The average NYC apportioned MBE_C emission rate for CO₂ was (57 ± 24) (1σ) kmol/s with 1σ representing 1 standard deviation from the mean MBE_C emission rate among the 7 research flights. The variability among the daily ensemble mean values (i.e., daily variability) of 42.1% is caused by both emissions variability and methodological uncertainties and is similar to the calculated 44% variability shown by the MBE_bulk estimates. The significant day-to-day variability is the consequence of the irregular spatiotemporal sampling of temporally and spatially varying emissions (aliasing of emission sampling) due to changing meteorological conditions between days, and thus flight plan orientations, as has been discussed in related works (Hajny et al., 2022; Pitt et al., 2022), as well as other methodological uncertainties. As shown in Figure 5a and b, the campaign-level variability, that is, the comparison of the campaign means as a function of the different ensemble components, induced by the priors and transport models was only 2.0% and 3.9%, respectively. However, the variability is larger for each individual flight or transect, as illustrated in Figure S7. As such, we calculated the mean variability, which refers to the prior, transport, and transect variability calculated for each day, averaged across days. This is, for each day, we calculate the variability induced by the transport (or the prior) and then we average this variability for all days in the campaign with the aim of quantifying the expected uncertainty in a single day estimation, on average. The mean variability across the 7 days were calculated to be 7.6%, 14.1%, and 19.1% for prior, transport, and transect, respectively. We note that the transect variability stems from applying the spatial allocation Φ calculated by transect to the bulk MBE estimation calculated using kriging, which uses all transects at once. The MBE_C is consistently greater (≈3.8%) when using kblt2 (Kantha and Clayson, 2000) vertical mixing parameterization rather than kblt5 (Hanna, 1982), indicating that kblt2 is systematically more dispersive than kblt5, as also shown by Pitt et al. (2022) and Hajny et al. (2022). It is worth highlighting that we rely on the ensemble of several latest generation models to quantify the uncertainty introduced by the transport model. Relying on the mean of a plausible group of transport models mitigates potential errors in any one particular transport model while also allowing us to estimate the contribution of transport uncertainty, which is quantified as the spread in the posterior emission rate across the transport models, to the overall posterior emission rate uncertainty.

Finally, biogenic emissions were accounted for using VPRM to allow for an appropriate comparison between the measured and modeled enhancements. In general, the emissions modeled with VPRM during these flights tend to be small. The nearly identical emission rates calculated using only anthropogenic simulated enhancements and total simulated enhancements (i.e., anthropogenic + biogenic) (Figure S8) indicate that the simulated biogenic fluxes exhibited a certain level of spatial homogeneity during the flights. Consequently, their impact on the background was found to be minimal for this set of flights. However, it is important to note that the emission rate calculated for the AOI remains to be the total emission rate as it includes biogenic fluxes within the AOI. Nevertheless, these fluxes make only a negligible contribution for these set of flights, as argued by Pitt et al. (2022). In addition, a recent work has shown that the biosphere influence has a <10% impact for MBE style experiments conducted outside of summer and that utilize the downwind edges to account for background influence (Lal and Kort, 2023).

To evaluate the performance of the MBE source apportionment approach, we compared the same flights with complementary approaches in quantifying GHG emissions for the NYC area. This includes a nested Bayesian inversion (IV) framework (Pitt et al., 2022) and a spatially explicit scaling factor (SF) approach (Hajny et al., 2022).

Figure 6 shows the results from the 7 research flights comparing traditional MBE approach kriging emission rate values and MBE_C derived in this work along with the posterior emission rates calculated from SF and IV approaches, for comparison. It is important to note that both Hajny et al. (2022) and Pitt et al. (2022) analyzed 2 additional research flights (e.g., March 26, 2019 and March 7, 2020) that were not investigated in this work due to a limited number of downwind transects resulting in highly variable kriging outputs. As mentioned earlier, our approach uses the same transport models and prior emissions as the SF and IV approaches allowing for a direct comparison of the calculated CO₂ emission rate between the different estimation schemes. The percent difference in emission rates with respect to the average SF and IV improved from 24.5% using the MBE_bulk kriging to 6.5% with the source apportioned MBE_C. As shown in Figure 6, a significant improvement in the campaign average CO₂ emission rate was observed when comparing the emission rate derived from only kriging (blue line) and the NYC-specific source apportionment MBE_C approach (black line) to the average posterior emission rate reported between SF and IV schemes (purple line). We also compared the standard deviation among methods, and the correlation between the MBE_C daily emission rates and those from the IV and SF, as summarized in Table 1. First, we compared the performance of traditional MBE kriging and source-apportioned MBE_C to the reported emission rate from SF and IV. The source-apportioned MBE_C consistently shows a lower root mean squared error (RMSE relative to the emission rates from SF and IV, respectively) compared to traditional MBE kriging suggesting better agreement of the daily means and temporal variability with the complementary approaches. In fact, the coefficient of determination (R²) between the daily estimates also shows a large improvement; from 0.50 to 0.58 when comparing the source-apportioned MBE_C to the SF approach and from 0.39 to 0.73, when comparing the source-apportioned MBE_C to the IV approach. Considering all these statistical evaluations between the traditional MBE kriging and the source-apportioned MBE_C, and comparing to the independent methods, we conclude that the source apportionment method can be used to apply MBE estimates to a specific AOI, and with less variability in results.

The aircraft mass balance method has been used in several studies at different scales and sampling patterns (Kalthoff et al., 2002; Mays et al., 2009; Karion et al., 2013; Cambaliza et al., 2014; Gioli et al., 2014; Karion et al., 2015; Heimburger et al., 2017; Tadić et al., 2017; Hajny et al., 2019; Pitt et al., 2019; Ryoo et al., 2019) due to its simple execution that does not require atmospheric transport modeling efforts. However, the traditional MBE method is limited and impacted by complex regional emissions and by the chosen background definition, and the spatial representativeness of the observed MBE_bulk is difficult to define. The source apportionment MBE approach discussed in this work account for the challenges in determining an appropriate background definition by scaling the MBE results based on the model-derived fractions to estimate a spatially allocated emission rate allowing us to consistently compare mass balance estimations with spatially explicit modeling results from both the Bayesian inversion (i.e., probabilistic emission estimates) and the SF approach (i.e., scaling up the inventory based on the model measure comparison). This enables the observed emission enhancements from the MBE to be isolated for a defined area and provides a solution to interpret MBE_bulk calculations for specific AOIs in complex emissions scenarios.

The proposed method here introduces a unique aspect by incorporating atmospheric transport information into the mass balance technique, which offers a conceptually different approach of estimating emissions and calculating the background. While the 3 methods are using measured enhancements as the basis for their estimation, both IV and SF methods utilize simulated enhancements to estimate the emissions based on the model-data mismatch. On the other hand, the source-apportioned MBE_C approach uses an integral approach based on the mass conservation equation in steady state to estimate the emissions. After the emissions are estimated, we then use the model to estimate the fraction of emissions originating in the AOI. In this approach, model-data mismatch is not used as the basis of the emissions estimation.

The top-down aircraft-based MBE method is an effective way to quantify urban emissions due to the efficiency at capturing the entire plume of a large/area emission source within hourly timescales. However, this approach assumes negligible influence from emission sources outside of the AOI, and therefore it is a challenge to attribute the estimated emissions to a particular area under complex emission environments.

In this work, we developed a new methodology to spatially allocate CO₂ bulk emissions estimated using the MBE method to a specific AOI in a complex emissions scenario, using the 5 boroughs of NYC as a highly relevant example. A campaign-average emission rate of (57 ± 24) (1σ) kmol/s from 7 flights between 2018 and 2020 was calculated using this source apportionment MBE approach. The results obtained using the method proposed in this work were comparable to the NYC 5 boroughs emission rates derived from ACES ([45 ± 9 ] kmol/s) and Vulcan ([52 ± 10] kmol/s) bottom-up inventories for a timeframe representative of our flights, as well as to other top-down estimation methods including the SF (Hajny et al., 2022) and IV schemes (Pitt et al., 2022), reducing the difference from 24.5% when utilizing the conventional MBE approach to 6.5%. In contrast to the traditional MBE, this source apportionment MBE approach can be used for cities surrounded by multiple and spatially complex emission sources. However, the method inherently relies on prior emission and transport and dispersion modeling, resulting here in a transport model uncertainty of 14.1% and a prior uncertainty of 7.6%, on average at a daily scale. In addition, the daily flight-to-flight variability accounted for 42.1%. Therefore, it is recommended to accompany this method with an ensemble of prior and transport models (as performed here) rather than a single model run to partially mitigate dependency on particular modeling choices and to enable assessment of the variability in the results caused by transport model and inventory uncertainties.

Aircraft measurement data used in this study are available from the Harvard Dataverse at https://doi.org/10.7910/DVN/UECQXX. In addition, the publicly available VPRM data can be found in Gourdji (2021), Biospheric CO₂ surface flux estimates from the Vegetation Photosynthesis and Respiration Model (VPRM) in the eastern USA and Canada: November 2016 to December 2020, National Institute of Standards and Technology, https://doi.org/10.18434/mds2-2382 (accessed July 11, 2023).

The supplemental PDF file for this article includes the following:

Figure S1. Linear fits for each transect in the research flight on 04 March 2020.

Figure S2. Extrapolation methods to characterize the downwind curtain between surface and PBL height for the research flight on 04 March 2020.

Figure S3. Modeled footprint along a single downwind transect and comparison of modeled Enh_TOT, Enh_NYC, and linear background subtracted Enh_TOT. An example case for fraction (Φ) greater than 1.

Figure S4. Estimating CO₂ measurements above the highest altitude observations using an exponential growth function for research flight 01 March 2019

Figure S5. Comparison of the kriged CO₂ flux distribution as a function of height above ground and horizontal distance for all relevant fights in this work for the synthetic transect extrapolation case.

Figure S6. Comparison of emission rates derived from independent kriging analysis including the Matlab-based EasyKrig3.0 (EK) with the linear variogram is compared to those derived from the R-based R-Krig (RK) with the best fit of the 4 variogram models tested for the synthetic transect extrapolation case.

Figure S7. Source attributed mass balance emission rate of NYC (i.e., MBE_C) for the different research flights (rows) separated by prior, transport, and individual transects (columns).

Figure S8. Comparing the role of biogenic fluxes in mass balance estimation of NYC (MBE_C) for flights conducted during the non-growing season.

Table S1. Temporal and spatial resolution of the inventory data and transport models used in this work.

The authors would like to thank the Purdue University Jonathan Amy Facility for Chemical Instrumentation for their help and support in maintaining equipment in ALAR and the National Institute for Standards and Technology (NIST; Award No. 70NANB21H021) for support of this work. The authors would also like to thank Anna Karion for her guidance and support. The authors acknowledge the support from air traffic control in and around NYC for their cooperation and support in completing our desired flight patterns.

The commercial equipment and instrumentation used in this work is limited to specificity purposes of the experimental methods. The National Institute of Standards and Technology (NIST) does not intend to endorse the equipment to be the best available option for the purpose of the study.

Purdue University and Stony Brook University groups acknowledge funding from NIST via Award Nos. 70NANB19H167 and 70NANB21H021.

The authors have no competing interests. Of note, Paul B. Shepson is an Associate Editor at Elementa, however, he was not involved during the peer-review process for this article.

ILC and PBS devised this project. JMT, ILC, KDH, JRP, RK, BHS, TJ, CRF, RC, and PBS contributed to flight measurements and flight design. JRP, KDH, and ILC provided modeling and inventory data. JMT, ILC, KDH, JRP, and PBS contributed to the analysis and interpretation of the data. JMT, ILC, and PBS wrote the manuscript with contribution from all coauthors.

How to cite this article: Tomlin, JM, Lopez-Coto, I, Hajny, KD, Pitt, JR, Kaeser, R, Stirm, BH, Jayarathne, T, Floerchinger, CR, Commane, R, Shepson, PB. 2023. Spatial attribution of aircraft mass balance experiment CO₂ estimations for policy-relevant boundaries: New York City. Elementa: Science of the Anthropocene 11(1). DOI: https://doi.org/10.1525/elementa.2023.00046

Domain Editor-in-Chief: Detlev Helmig, Boulder AIR LLC, CO, USA

Knowledge Domain: Atmospheric Science

Associate Editor: Lori Bruhwiler, Global Monitoring Division, National Oceanic & Atmospheric Administration Earth System Research Laboratory, Boulder, CO, USA