Quantifying the contribution of major carbon producers to increases in vapor pressure deficit and burned area in western US and southwestern Canadian forests

Observed monthly maximum temperature, minimum temperature, and vapor pressure for North America at a ¼ degree horizontal resolution during 1901–2021 was acquired from Williams et al (2022). We calculated monthly mean VPD at the native pixel scale by first calculating saturation vapor pressure from monthly average maximum and minimum temperature and having VPD set as the difference between actual and saturation vapor pressure following Seager et al 2015.

Monthly maximum temperature, minimum temperature, and specific humidity data from 28 coupled model intercomparison project, phase 6 (CMIP6) climate models were acquired for both the historical (1850–2014) and SSP2-4.5 (2015–2100) experiments (table S1) and interpolated to a common 1.0-degree horizontal resolution grid before calculating VPD as above. We then calculated the change in VPD over time by multiplying the linear least squares trend by the number of years in the time period of interest. Unless otherwise specified, throughout this study we used the same linear least squares approach to calculate the change in all modeled and observed climate variables over time.

We constrained our analysis to forested lands based on level III U.S. Environmental Protection Agency (EPA) ecoregions of the western United States and southwestern Canada in the area west of 103°W and between 30° and 53°N. We focused on March–September mean VPD over forested areas of this domain given the influence of antecedent and concurrent VPD in curing fuels and making them receptive to igniting and carrying fire (e.g. Abatzoglou et al 2013, Jolly et al 2015) as evident in strong relationships between seasonal aggregates of VPD and macroscale BA in the region (e.g. Abatzoglou et al 2021).

To evaluate relationships between the change in GMT and the change in March–September VPD averaged across the study area in both observations and climate models, we acquired observed GMT from NASA-GISS from 1880–2021 (Hansen et al 2010) and modeled GMT for each unique CMIP6 climate model and ensemble member. We calculated the rate of change in regional March–September VPD per degree change in GMT as the difference in VPD and GMT between 1991–2020 and 1901–1950 (equation (1)),

$$\begin{equation}\frac{{\Delta {\text{VPD}}}}{{\Delta {\text{GMT}}}} = \frac{{{\text{VP}}{{\text{D}}_{1991 - 2020}} - {\text{VP}}{{\text{D}}_{1901 - 1950}}}}{{{\text{GM}}{{\text{T}}_{1991 - 2020}} - {\text{GM}}{{\text{T}}_{1901 - 1950}}}}.\end{equation} \tag{ 1 }$$

2.2.1. Assessing the influence of internal variability

2.2.2. Assessing linear scaling

After establishing ΔVPD/ΔGMT, we simulated changes in GMT using an energy balance carbon cycle climate model so that the ΔVPD/ΔGMT relationship could be applied to the output. The model is based on (Millar et al 2017) with parameters defined by Ekwurzel et al (2017).

We first ran a 'full forcing' (i.e. including volcanic eruptions, solar activity, total fossil fuel emissions, etc) GMT simulation with energy balance and carbon cycle parameters that provided the best fit to observed changes in atmospheric CO₂ and GMT since 1750 (see Licker et al 2019, supplementary tables 1 and 2). We then ran the model again with the scope 1 and scope 3 carbon dioxide and methane emissions of the world's 88 largest carbon producers (hereafter 'the big 88') removed between 1900 and 2015 using annual, producer-specific emissions data (Licker et al 2019 table 7, Climate Accountability Institute 2020). This enabled us to calculate the annual GMT change resulting from emissions traced to the big 88, $\Delta {\text{GM}}{{\text{T}}_{{\text{big}}88}}\left( t \right).$ We refer to this as the 'X88' simulation.

Annual emissions data from the big 88 is not uniformly available for the post-2015 period—a period of significant wildfire activity in the study region. Rather than simulating $\Delta {\text{GM}}{{\text{T}}_{{\text{big}}88}}$ for 2016–2021, we hold $\Delta {\text{GM}}{{\text{T}}_{{\text{big}}88}}\left( t \right)$. Two lines of evidence suggest this approach will yield conservative results. First, annual emissions data for a subset of the big 88 between 2015 and 2018 show that each entity's emissions increased every year through 2018 (figures S2 and S3). Second, average annual global CO₂ emissions from fossil fuel and industry were higher over the 2016–2021 period than in 2015 (Friedlingstein et al 2021, Liu et al 2022). Results for the full 1900–2021 period are reported here. Results for the 1901–2015 period are reported in table S2.

Consistent with previous studies, we also ran the energy balance carbon cycle climate model with and without total fossil fuel aerosols for both the full forcing and the X88 cases (Licker et al 2019). Since the big 88 contribute roughly two-thirds of the total industrial aerosol emissions, however, excluding the aerosols underestimates the net cooling influence of that forcing. We therefore focus on the best estimate with full forcing, including total fossil fuel aerosols.

The order of removal of each of the 88 carbon producers within the model has a potential, albeit slight, influence on its respective contribution to GMT. To assess this source of uncertainty, for each carbon producer, we first ran the model with only that producer's emissions, i.e. without the emissions of the remaining 87. We then ran the model again assuming that the producer's emissions were added last in the set of 88.

We then determined the simulated difference in annual GMT between the full forcing and X88 simulation results. These annual $\Delta {\text{GM}}{{\text{T}}_{{\text{big}}88}}$ findings were then incorporated into VPD and BA time series based on the 28 CMIP6 climate models and observational data.

Forested BA data were acquired from three datasets: (a) Monitoring Trends in Burn Severity in the western US from 1984–2020 (Eidenshink et al 2007); (b) Canadian National Burned Area Composite in Canada from 1986–2020 (Hall et al 2020); and (c) MODIS BA (MCD64A1) from 2001–2021 (Giglio et al 2018). The latter was used to estimate BA for 2021 based on a linear model between BA from MODIS and each of the respective datasets during periods of overlap following Williams et al 2014. We limited our analysis to BA in the western US and southwestern Canada within forested level III ecoregions to match the geographic extent of the VPD data.

To calculate the portion of the VPD change from 1901–2021 resulting from emissions traced to the big 88, we constructed counterfactual VPD time series for the observational data set and for each CMIP6 model. The counterfactuals represent VPD change without the influence of the big 88 and were developed by multiplying the annual $\Delta {\text{GM}}{{\text{T}}_{{\text{big}}88}}$ for each year by the ΔVPD/ΔGMT identified previously for each dataset. The VPD change contributed by the big 88 was calculated as the difference between the total modeled or observed VPD change and the counterfactual (equation (2)),

$$\begin{align}{\text{VP}}{{\text{D}}_{\text{c}}}\left( t \right) & = {\text{VP}}{{\text{D}}_m}\left( t \right) - {(\Delta {\text{VPD}}/\Delta {\text{GMT}})_m} \nonumber\\ & \quad \cdot \Delta {\text{GM}}{{\text{T}}_{{\text{big}}88}}\left( t \right) \end{align} \tag{ 2 }$$

where t represents year, and VPD_c represents the counterfactual VPD for each model, m.

Following Abatzoglou and Williams (2016), we considered counterfactual BA scenarios that exclude the modeled influence of the big 88 from the observed VPD record. Because reliable BA observations are only available starting in the mid-1980s, our counterfactual BA analysis covers the 1986–2021 period. Model results indicated that the influence of the big 88's emissions grew substantially between 1901 and 1986. To avoid devaluing the importance of the big 88's emissions in shaping climate before 1986, we report VPD contributions since both 1901 and 1986, and BA contributions since 1986.

To calculate the counterfactual BA, we used a simple macroscale empirical model of annual forest BA of the form:

$$\begin{equation}\log \left( {{\text{BA}}\left( t \right)} \right) = \alpha + \beta \cdot {\text{VPD}}\left( t \right) + \varepsilon \end{equation} \tag{ 3 }$$

where t is the year, α and β are regression coefficients, and represents an error term (Abatzoglou et al 2021). This modeling exercise is based on the strong exponential VPD-BA relationship in the observational record (Juang et al 2022). The error term is drawn from residuals in the aforementioned equation when omitting and is meant to capture stochasticity in relationships not captured in VPD (e.g. fire weather conditions, human ignitions).

As above, counterfactual VPD time series excluding the influence of the big 88 and the time varying influence of the big 88 on GMT were constructed based on the observed VPD time series following equation (4):

$$\begin{align}{\text{VP}}{{\text{D}}_{{\text{c,obs}}}}\left( t \right) & = {\text{VP}}{{\text{D}}_{{\text{obs}}}} - {(\Delta {\text{VPD}}/\Delta {\text{GMT}})_m} \nonumber\\ & \quad \cdot \Delta {\text{GM}}{{\text{T}}_{{\text{big}}88}}\left( t \right).\end{align} \tag{ 4 }$$

This approach removes the time varying increase in GMT resulting from emissions traced to the big 88, leaving much of the interannual variability from the observed record unchanged. We then applied the counterfactual VPD time series (${\text{VPD}}$ _c,obs) for each model to equation (3) to estimate BA during 1986–2021 in the absence of carbon emissions from the big 88. Following previous work, we use counterfactuals from observations to best adhere with realized BA (Abatzoglou and Williams 2016). We ran each counterfactual model 100 times for each climate model to account for randomized time series of in modeled BA and present results for the median of each of the 28 CMIP6 models. The cumulative forest BA resulting from emissions traced to the big 88 was calculated by taking the difference between cumulative observed forest BA and cumulative forest BA simulated using the counterfactual VPD time series. We assessed the validity of this modeling approach by comparing modeled cumulative BA using equation (3) and the observed VPD time series data that include the influence of the big 88 (running the model 100 times to account for random terms) to the observed forest BA.

Observed regional March–September VPD as calculated by linear least squares regression increased by 0.82 hPa from 1901–2021, consistent with findings from previous studies (e.g. Williams et al 2022). This represents an 11.2% increase in VPD over the study period. We also note that 21 of the 22 years since 2000 have had March–September VPD above the 1901–1950 average. All CMIP6 models captured the historical increase in VPD, with a 28-model mean increase of 0.80 hPa (interquartile range (IQR) 0.53–1.03 hPa) over the same time period (figure 1(a)).

Zoom In Zoom Out Reset image size

The observational data showed a positive relationship between GMT and regional March–September VPD of 0.99 hPa °C⁻¹ or 13.4% °C⁻¹. All CMIP6 models also showed a positive relationship between GMT and regional March–September VPD, with a 28-model mean of 0.81 hPa °C⁻¹ (IQR 0.53 hPa °C⁻¹–1.05 hPa °C⁻¹) or 12.9% °C⁻¹ (figures 1(b) and (c)). The spread in ΔVPD/ΔGMT within the model suite likely reflects the combination of internal variability within each model and inter-model differences in physics. Estimates of the influence of internal variability on ΔVPD/ΔGMT for each model (figure S1), however, suggest that uncertainty in ΔVPD/ΔGMT from inter-model differences exceeds that from internal variability, similar to results from studies that have examined sources of uncertainty for regional temperature change (Hawkins and Sutton 2009).

The relationship between 21 year moving average changes in regional VPD and 21 year moving averages in GMT are quasi-linear from the mid-20th century onwards (figure 1(b)), suggestive of a scalable, albeit model dependent, GMT-regional VPD relationship. Linear scaling of ΔVPD/ΔGMT is further confirmed by noting that 22 of the 28 models had R² > 0.75 (all models had R² > 0.5) for a linear relationship between 21-year moving average regional VPD and GMT during the 1901–2020 period. These results suggest that linear scaling adequately captures the ΔVPD/ΔGMT response.

The period from 1986–2021, for which reliable BA data is available for the region, is characterized by steep, steady rise in VPD. As calculated by linear least squares regression over 1986–2021, March–September VPD increased by 1.09 hPa in the observational data and by 1.13 hPa (IQR 0.71–1.43 hPa) in the 28-model mean over that time period.

Consistent with the results of Ekwurzel et al (2017), we find that the influence of the big 88 on GMT grows over time. In the early years of the 20th century, $\Delta {\text{GM}}{{\text{T}}_{{\text{big}}88}}\left( t \right)$ is nearly negligible. However, $\Delta {\text{GM}}{{\text{T}}_{{\text{big}}88}}\left( t \right)$ increases exponentially to 0.49 °C by 2021. The best-estimate climate model accounts for 49.1% ± 1.3% of the observed warming between 1900 and 2015 in our X88 simulation.

Applying the established ΔVPD/ΔGMT relationship, the total contribution of emissions from the big 88 to VPD change for the 1901–2021 period is 0.47 hPa per the observational data, or 57% of total observed change in VPD. The total contribution of the big 88's emissions to VPD change per the 28-model mean is slightly lower at 0.39 hPa (IQR 0.25–0.50), or 48% (IQR 32%–63%) of the total modeled VPD change over the study period (figure 2). The models and observations thus imply that roughly half of the increase in VPD from 1901–2021 can be traced to the emissions of these major carbon producers.

Zoom In Zoom Out Reset image size

Over the 1986–2021 time period, the total contribution of emissions from the big 88 to VPD is 0.33 hPa in the observational data and 0.27 hPa (IQR 0.18–0.35 hPa) for the 28-model mean. Given this shorter time period, these results may be more prone to the influence of internal variability within the models than results from the full 1901–2021 time period.

As in previous studies, forest BA grows exponentially with increased VPD (figure 3). Using VPD time series that do not exclude the influence of the big 88, we approximate modeled cumulative BA of 217 100 km², which is comparable to the observed cumulative BA of 214 690 km². The differences between observed and modeled BA were not statistically significant.

Zoom In Zoom Out Reset image size

We modeled BA that excluded the influence of the big 88 for each of the 28 CMIP6 models using a counterfactual VPD time series that removes the model dependent ΔVPD/ΔGMT from the observed VPD data. Without emissions traced back to the big 88, we find that a central estimate of cumulative BA in the region would have been 134 735 km² (IQR 112 000–156 000 km²). We thus conclude that emissions from the big 88 contributed 37% (IQR 26%–47%) of the cumulative burned forest area from 1986–2021 (figure 4). This represents roughly 80 000 km² (nearly 20 million acres) of burned forest area, or an area larger than Ireland.

Zoom In Zoom Out Reset image size