1 Introduction

Quantitatively assessing the responsibility of individual nations for the impacts of climate change is necessary to develop and evaluate equity-based mitigation agreements (Okereke & Coventry, 2016) and inform legal discussions regarding climate liability (Burger et al., 2020; Stuart-Smith et al., 2021). Furthermore, precisely quantifying whose emissions have caused damage to whom, and how much, is critical to advance emerging legal discussions about how courts might apportion liability among the world’s emitters (Burger et al., 2020). Climate change attribution work has long been motivated by the possibility of being used for liability claims (Allen, 2003; Marjanac & Patton, 2018) and scientific efforts to assess national responsibility for climate change have been undertaken in earnest since Brazil’s 1997 proposal that mitigation requirements be proportional to each country’s historical responsibility for warming (M. Den Elzen & Schaeffer, 2002; Prather et al., 2009).

Despite the history of work on national attribution, individual nations’ emissions have not been quantitatively linked to the downstream economic impacts of warming, which requires an integrated analytical framework. Many previous studies have analyzed individual steps in the causal chain from emissions to impact. Existing national attribution work has quantified contributions to global emissions (Ciais et al., 2013; M. G. Den Elzen et al., 2013; Höhne et al., 2011; Matthews, 2016), GHG concentrations (Ciais et al., 2013; Wei et al., 2012), global radiative forcing (Li et al., 2016; Murphy & Ravishankara, 2018; Ward & Mahowald, 2014), or global mean temperature change (Höhne et al., 2011; Matthews et al., 2014; Skeie et al., 2017; Wei et al., 2012). At the same time, there has been a broadened appreciation of how varied regional temperature responses can be in response to global mean warming (Deser et al., 2012), which has complicated the attribution of climate change at local scales by requiring large ensembles of simulations with a single model. Separately, advances in applied econometrics have allowed scientists to tie historical warming to large effects on observed national economic productivity (Burke et al., 2015; Dell et al., 2012; Diffenbaugh & Burke, 2019). However, no study has combined these three literature bases in a coherent framework to quantify changes in economic output attributable to individual emitters. The central contribution of this study is to combine historical data with climate models of varying complexity to quantify the contribution of each nation’s emissions to the historical economic effects of warming in every other nation. This analysis represents an important step forward for at least four reasons.

Firstly, the legal requirement to prove standing in order to sue for climate-related harms means that plaintiffs must credibly link the damage they have suffered to a specific actor responsible for that damage (Burger et al., 2020). Attributing economic damage to global warming writ large is not considered legally sufficient to prove the causal role of a specific emitter (Marjanac & Patton, 2018), so existing studies that have quantified historical losses from climate change without connecting them to specific emitters do not support climate liability claims. The lack of scientific evidence linking individual emitters to the downstream impacts of warming has been a principal evidentiary gap in climate litigation (Stuart-Smith et al., 2021), so an analysis extending attribution of economic losses from warming to specific emitters is a key step forward.

Secondly, previous studies linking emissions to warming or warming to damages have often used different data, modeling approaches, or spatiotemporal resolutions, making it difficult to quantify national responsibility for economic impacts from existing literature. For example, many studies linking individual nations’ emissions to climate change have done so at the global (M. G. Den Elzen et al., 2013; M. Den Elzen & Schaeffer, 2002; Höhne et al., 2011; Matthews, 2016; Matthews et al., 2014; Skeie et al., 2017) or grid-cell (Wei et al., 2012) scale, as opposed to the district or country scale used in empirical climate-economy studies (Burke et al., 2015; Burke & Tanutama, 2019; Dell et al., 2012; Kalkuhl & Wenz, 2020). Two recent studies connected individual emitters to regional (Lewis et al., 2019, p. 201) or country-level (Beusch et al., 2022) temperature changes, but both did so only for the top few emitters and focused on the next several decades instead of the historical period, making their analyses unsuited for quantifying national contributions to historical losses from warming and thus present liability claims.

Thirdly, nonlinearities in the coupled climate-economy system make any inference of national responsibility for climate damages from existing studies irresponsible (den Elzen et al., 2005; Höhne & Blok, 2005; Trudinger & Enting, 2005). The efficiency of carbon sinks in the land and ocean decreases with increasing emissions, the economic growth penalty of rising temperatures is more severe at higher temperatures (Burke et al., 2015), and reductions in growth compound exponentially over time (Moore & Diaz, 2015). Identical GHG emissions from different countries at different times therefore may have different economic effects. As a result, studies showing the effect of overall historical warming on the global economy (Diffenbaugh & Burke, 2019), do not, on their own, position an assessment of individual nations’ responsibility for those effects (Burger et al., 2020). A key contribution of our work is to interactively simulate the carbon cycle and climate response to individual nations’ emissions as part of an analytical pipeline that connects these responses to country-level temperature changes.

Finally, uncertainties compound at each step in the causal chain from emissions to impact (Prather et al., 2009). Uncertainty can provide emitters with a veil of plausible deniability against the damages suffered by claimants. As such, responsibly assessing the culpability of an individual emitter for the climate damages suffered by another country (i.e., the signal) requires a robust estimate of the uncertainty in that culpability (i.e., the noise). Uncertainties in the carbon cycle response to emissions complicate attribution of global warming to individual emitters (Höhne et al., 2011; Karstensen et al., 2015; Wei et al., 2014). Internal climate variability increases uncertainty further, as many local (or country-level) temperature changes are consistent with the same magnitude of global temperature change (Deser et al., 2012). Uncertainty in the economic response to warming has dominated uncertainty in overall climate damage estimates (Burke et al., 2018), but it is not obvious which sources of uncertainty are most important in calculations of damages attributable to individual nations. Further, explicitly comparing uncertainties within the same framework will provide information to decision-makers about where to focus scientific efforts in support of climate liability discussions. Such work has not yet been performed, as it requires an integrated end-to-end attribution analysis, such as the one we present here.

Here we assess the economic effects of the warming attributable to each country’s emissions. We use a simple carbon-climate model to estimate the contribution of countries’ emissions to global mean surface temperature (GMST) change, propagate these global contributions to the country scale using a pattern scaling approach trained on existing global climate model (GCM) output, and apply existing temperature-growth relationships to calculate the economic effects of these country-level temperature changes. To estimate uncertainty at each step in the causal chain, we run a series of perturbed-parameter simple climate model simulations, use both single- and multi-model GCM ensembles to propagate global forcing to the country level, and resample the temperature-growth relationship (Methods). This methodology yields a distribution of two million possible values of damage suffered by one country attributable to another country in each year, a distribution which we leverage to test the statistical significance of this damage (Methods). While the burden of proof for demonstrating liability is high and nonlinearities in the climate system make any such attribution difficult, our work is an important contribution to the effort to develop evidence that can be used to make claims for legal standing.

2 Data

We use 1960–2014 temperature observations from Berkeley Earth (Rohde & Hausfather, 2020), the 20th Century Reanalysis (Compo et al., 2011), and the University of Delaware (Willmott, 2000) (UDel), regridded to a common 1°-by-1° resolution and averaged together. We similarly average precipitation from Global Precipitation Climatology Centre (Schneider et al., 2011) and UDel at a 0.5°-by-0.5° resolution. We average multiple datasets to reduce observational uncertainty associated with any single gridded data product, as ensembles of datasets have been shown to outperform individual datasets (Gottlieb & Mankin, 2021). However, we note that we find very similar regression results using each individual data product, so this choice does not substantially affect our results (Supplementary Fig. 1). While our analysis focuses on the effects of temperature, we use precipitation data as a covariate in the statistical model (Methods).

This observational data is aggregated to population-weighted country means using year-2000 population data from the Gridded Population of the World (Center for International Earth Science Information Network, 2016) and merged with 1960–2014 gross domestic product per capita (GDPpc) data from the World Bank (Bank, 2016) and United Nations World Population Prospects country population data. Population weighting of the climate data focuses the analysis on populous regions within large, geographically heterogeneous countries. All spatial interpolation is performed with bilinear interpolation using Python’s “xarray” package.

We use two sets of climate model data to sample the uncertainty in model responses to forcing and their representation of internal climate variability: the first realization of each of 12 models from the sixth phase of the Coupled Model Intercomparison Project (Eyring et al., 2016; Gillett et al., 2016) (CMIP6, Supplementary Table 1) and 20 realizations of the CESM1 Large Ensemble (Deser et al., 2020; Kay et al., 2015) (CESM1-LE). The CMIP6 models sample uncertainty across model structure, while CESM1-LE allows us to isolate the role of internal climate variability in propagating global forcing to the local level (Hawkins & Sutton, 2009).

We use the CMIP6 “historical” and “historical-nat” simulations from 1850 to 2020; the historical-nat simulations extend to 2020, whereas the historical simulations extend only to 2014, so we splice them with the first 6 years of each model’s Shared Socioeconomic Pathway 5–8.5 (SSP5-8.5) simulation to extend these simulations to 2020 (Schwalm et al., 2020). From CESM1-LE, we use the single-forcing experiments named “ALL” and “XGHG” over 1920–2020 to represent a world with and without GHG emissions. For simplicity, we refer to all these experiments as “historical” and “natural.” GMST has risen faster in the CESM1-LE ensemble mean than in the CMIP6 ensemble mean (Fig. 1a), which may be due to structural model differences and to the fact that the CMIP6 models incorporate aerosols, which have a cooling effect. However, when considering individual models and realizations, the CMIP6 ensemble range encompasses the CESM1-LE realizations (Supplementary Fig. 2), so this discrepancy does not substantially affect our results.

Fig. 1
figure 1

Attribution approach. (a) FaIR GMST change from historical emissions (red) and historical-minus-United States emissions from 1990 to 2014 (blue). Shading shows the 95% ensemble range. Black lines show the ensemble mean GMST change from two GCM ensembles: the CMIP6, which provides an estimate of structural uncertainty across models, and the CESM1-LE, which provides an estimate of uncertainty from internal climate variability. GMST change is calculated by subtracting each simulation’s corresponding natural simulation. (b) Mean national temperature change from pattern scaling using FaIR and GCM output, with FaIR-based historical in red and GCM output in black. Shading shows the 95% ensemble range. Map shows each country’s pattern scaling coefficient (Methods). (c) Relationship between annual temperature and GDPpc growth relative to the optimum temperature. Shading shows the 95% bootstrap range (Methods). Map shows the marginal effect of a 1 °C increase in annual mean temperature for each country. (d) Change in India’s GDPpc attributable to historical emissions (red), historical-minus-United States (U.S.) emissions (blue), and the U.S. alone (yellow). Boxes span the 25–75th percentiles and whiskers extend to the 2.5–97.5th percentiles. Upper text denotes the P-value of a Kolmogorov–Smirnov (K-S) test on the two distributions, with N denoting the sample size

We use monthly near-surface air temperature (“tas”) from all ensembles, resampled to annual means and interpolated to a common 2°-by-2° grid. While our calculation of economic damages extends only to 2014, we use climate model simulations through 2020 to enable us to smooth these data with a centered running mean in the pattern scaling analysis (Methods).

To simulate the effects of individual nations’ emissions, we use the Finite amplitude Impulse Response (FaIR) v1.3 model (Millar et al., 2017; Smith et al., 2018). As inputs to FaIR, we use 1850–2014 national fossil fuel emissions data from the Community Emissions Data System (Hoesly et al., 2018) (CEDS), merged with national land use and land cover change (LULCC) emissions data from Houghton and Nassikas (2017). The CEDS data are the official emissions data for use in CMIP6 simulations, which avoids discrepancies between our FaIR simulations and the GCM output. These emissions data use “territorial” accounting, where emissions are assigned to the countries in which they occur. We also compare these results to 1990–2014 “consumption-based” emission data from the Global Carbon Project (Peters et al., 2012) (GCP) to account for emissions embodied in international trade (Davis & Caldeira, 2010), though consumption-based data is only available for fossil fuel CO2 emissions and not LULCC emissions.

3 Methods

Our approach is threefold: (1) Estimate national contributions to GMST change using a simple emissions-driven climate model; (2) Propagate these GMST changes into country-level temperatures using pattern scaling methods trained on output from coupled climate models; and (3) Calculate the resulting economic damages using empirical temperature-growth relationships.

3.1 National contributions to global temperature change in a simple climate model

To simulate the contributions of individual countries to global climate change, we use FaIR in a set of “leave-one-country-out” experiments. We run FaIR from 1850 to 2014 with total historical emissions to represent the historical scenario (Fig. 1a). Then, for each country in the CEDS data (N = 174), we simulate the same historical period with all emissions except those from that country (Fig. 1a). In each of the 174 sets of simulations, we remove a single country’s contributions to fossil fuel and LULCC carbon dioxide, methane, and nitrous oxides; we leave other chemical species, along with solar and volcanic forcing, unchanged. We also run a “natural” simulation where all emissions are fixed at 1861–1880 averages (Smith et al., 2018).

When we use consumption-based emissions, we scale the original CEDS territorial fossil fuel CO2 emissions by the ratio of consumption-to-territorial emissions from GCP. This procedure generates consumption-based emissions data that are internally consistent with the CEDS data, which avoids propagating absolute differences in emissions amounts between GCP and CEDS into our results.

To match our calculations of economic damage from temperature change, which are performed from 1990 to 2014 and from 1960 to 2014, our leave-one-out simulations with FaIR are forced with total historical emissions from 1850 until either 1990 or 1960, at which point the forcing removes a nation’s emissions. A start date of 1990 examines the counterfactual in which countries had phased out fossil fuels once the scientific consensus on climate change became clear; a start date of 1960 examines the counterfactual in which postwar economic development had occurred with renewable energy rather than fossil fuels. Note that global emissions include processes with ambiguous national affiliations, such as international shipping and aviation. We include these emissions in the global total, but CEDS does not assign them to individual nations, so they are never subtracted from that global total. As such, our leave-one-out simulations may not precisely sum to global totals.

Further discussion of the appropriateness of the “leave-one-out” method, explanation of the choice to use total emissions instead of per capita emissions, and details on how we implement the perturbed-parameter FaIR simulations are available in the Supplementary Information.

3.2 Propagating national contributions to warming to the country level

To extend the results of our FaIR simulations to the country level, we use a “pattern scaling” approach that relates GMST to country-level temperatures (Mitchell, 2003; Santer, BD et al., 1990; Tebaldi & Arblaster, 2014). We first calculate annual GMST from each model by averaging surface temperature, weighting by the square root of the cosine of latitude. We then take the difference between the historical and natural simulations for each model for both GMST and each country’s annual temperature, and we smooth these differences with an 11-year centered running mean to reduce interannual variability, a method which has been previously used to reduce the influence of internal variability in pattern scaling (Mitchell, 2003). The centered running mean means that the value for 2014 represents the average from 2009 to 2019, inclusive, which is why we splice the historical CMIP6 simulations with the SSP5-8.5 projections to allow this calculation.

We then linearly regress the simulated 1960–2014 country-level temperature change values onto the simulated 1960–2014 GMST change values for each model (Beusch et al., 2020; Mitchell, 2003; Tebaldi & Arblaster, 2014), motivated by the strong linear relationship between GMST change and local temperature change (Seneviratne et al., 2016). Other pattern scaling approaches use a “delta” method that compares epoch differences rather than fitting a regression (Tebaldi & Arblaster, 2014), but the regression method has been found to outperform delta method approaches (Lynch et al., 2017; Mitchell, 2003). This approach yields a regression coefficient for each country that describes its sensitivity to changes in GMST; coefficients range from 0.6 to 0.7 °C in the tropics, indicating slower warming than the global mean, to > 1.5 °C in the poles, indicating faster warming than the global mean and consistent with the classic “Arctic amplification” pattern (Fig. 1b, inset). We perform this regression over the model-simulated 1960–2014 period to maintain consistency in time periods with the rest of the analysis.

We then predict the time evolution of country-level temperatures in each FaIR realization using the FaIR GMST values and the above coefficients, in both the historical and leave-one-out simulations (Fig. 1c). Note that the predicted time series do not contain realistic interannual variability; in the next section, we show how we reference these predicted time series to observed temperatures to impute realistic variability to the counterfactual time series.

3.3 Attribution of climate damages to each country

The final step in our analysis is to extend these reconstructed country-level temperature changes into income changes. Here we apply existing methods for empirically estimating the effect of temperature changes on economic growth rates. Our main analysis uses the global nonlinear approach developed by Burke et al. (2015), though we also test several alternative functional forms (Supplementary Information). We apply the regression specification of Burke et al. (2015) to our data, rather than relying on their published parameter estimates, so that each step of our analysis uses the same data and spatiotemporal scale and so that we can easily propagate uncertainty through our analysis. Specifically, we estimate the following model with ordinary least squares:

$${g}_{ct}={\upbeta }_{1}{T}_{ct}+{\upbeta }_{2}{T}_{ct}^{2}+{\uppi }_{1}{P}_{ct}+{\uppi }_{2}{P}_{ct}^{2}+{\upgamma }_{c}+{\upmu }_{t}+{\uptheta }_{c1}t+{\uptheta }_{c2}{t}^{2}+{\upepsilon }_{ct}$$
(1)

In Eq. 1, g denotes GDPpc growth, T denotes population-weighted country-mean annual mean temperature, and P denotes population-weighted country-mean precipitation (annual mean of monthly total rainfall). The country fixed effect γ controls for time-invariant country-specific factors such as geography and the year fixed effect μ controls for common global shocks. θ denotes country-specific linear and quadratic time trends, which account for secular trends in output due to factors such as technological change. To sample uncertainty in the β parameters, we use 250 bootstrap iterations to re-estimate the parameters after resampling with replacement by country, which accounts for autocorrelation in errors within each country (Burke et al., 2018). Our results may differ from the exact parameter values found in Burke et al. (2015) due to different data and a slightly longer period of analysis, but the quadratic relationship between temperature and growth we find is very similar to previous estimates (Fig. 1c).

To calculate the marginal effect of an additional degree in the annual mean temperature (Fig. 1c, inset), we differentiate Eq. 1 with respect to temperature.

$$\frac{\partial {g}_{ct}}{\partial {T}_{ct}}={\upbeta }_{1}+2*{\upbeta }_{2}*{T}_{ct}$$
(2)

We use this temperature-growth relationship to compare observed growth in each country with the counterfactual growth that would have occurred without each other country’s emissions. ∆Thist denotes the historical temperature change from FaIR and ∆Thist−a denotes the temperature change induced by the leave-one-out simulations for country a, with distinct values for each GCM, FaIR realization, country, and year.

For each country, GCM, and FaIR realization, we construct two counterfactual temperature time series:

$$\updelta {T}^{\mathrm{hist}}=T-\Delta {T}^{\mathrm{hist}}$$
(3)
$$\updelta {T}^{a}=T-\left(\Delta {T}^{\mathrm{hist}}-\Delta {T}^{\mathrm{hist}-a}\right)$$
(4)

Here, T refers to observed temperatures, δThist refers to temperatures that would have occurred without any anthropogenic climate change, and δTa refers to temperatures that would have occurred without country a’s emissions. Supplementary Fig. 3 provides a schematic of this calculation. These counterfactual temperatures are referenced to the observed time series to: (1) impute realistic interannual variability to these smooth counterfactual time series; and (2) bias-correct the model output by using differences rather than absolute temperatures.

To calculate the resulting economic damage for each country, year, FaIR realization, GCM, and damage function bootstrap, we first calculate the effect of temperature change on economic growth ∆g using the parameters β1 and β2, where ∆g is positive when countries would have grown faster without climate change:

$$\Delta {g}^{\mathrm{hist}}=\left[{\upbeta }_{1}\left(\updelta {T}^{\mathrm{hist}}\right)+{\upbeta }_{2}{\left(\updelta {T}^{\mathrm{hist}}\right)}^{2}\right]-\left[{\upbeta }_{1}\left(T\right)+{\upbeta }_{2}{\left(T\right)}^{2}\right]$$
(5)

Calculating the change in growth this way subtracts out the non-climate determinants of growth and isolates the change associated with temperature. This difference in growth is then applied to the beginning of the observed GDPpc time series and integrated to calculate a counterfactual GDPpc time series (cGDPhist) for each country, model, FaIR realization, and damage function bootstrap:

$$cGD{P}_{t}^{\mathrm{hist}}=cGD{P}_{t-1}^{\mathrm{hist}}+\left[cGD{P}_{t-1}^{\mathrm{hist}}*\left({g}_{t}+\Delta {g}_{t}^{\mathrm{hist}}\right)\right]$$
(6)

Economic damage D is calculated as the difference between a country’s actual GDPpc time series and this counterfactual time series. GDPpc values are multiplied by population to convert to absolute GDP where appropriate.

We first perform the above calculation using the historical counterfactual time series δThist and the observed time series T, yielding total historical damage Dhist. We then compute the leave-one-out damage value Dhist−a, which refers to the damage done by all emissions other than those of country a, by repeating the counterfactual growth and GDP calculation using δThist and δTa (Supplementary Fig. 3). The economic damage attributable to country a (Da) for each country, year, GCM, FaIR realization, and damage function bootstrap is the difference between these damage values:

$${D}^{a}={D}^{\mathrm{hist}}-{D}^{\mathrm{hist}-a}$$
(7)

Note that Da can take on both negative and positive values, representing economic damages and benefits resulting from country a’s emissions, respectively. Once we determine total dollar effects in individual nations, we aggregate them to total global income changes; where appropriate, we show cumulative attributable damages, which refer to the sum over time of attributable income changes.

One key caveat in our work relates to the emission of anthropogenic aerosols, which have a cooling effect. Based on analysis presented in detail in the Supplementary Information, the varied treatment of aerosols across models does not drive model differences in attributable damages (Supplementary Fig. 4).

3.4 Statistical tests of significant national contributions to damages

To determine whether attributable damage for country a is statistically significant, we compare the distributions of damage with (Dhist) and without (Dhist−a) the emissions of country a. For each country and year, we use a Kolmogorov–Smirnov (K-S) test to determine whether these two distributions are statistically distinguishable (Fig. 1d). P-values less than 0.05 are considered significant, indicating that the “with a” and “without a” distributions are not likely drawn from the same underlying distribution. When aggregating total attributable damages, we only aggregate from countries and years that pass this significance test.

We use a K-S test both because it is nonparametric and because it incorporates both the location and shape of the damage distributions, allowing us to leverage our uncertainty quantification analysis. However, repeating our analysis using a Student’s t-test yields similar results, so our results are not sensitive to this choice (Supplementary Fig. 5).

3.5 Uncertainty partitioning

Our analysis yields an array Da, which refers to the cumulative effect of country a’s emissions on country c, with uncertainty from the FaIR carbon cycle parameters, the global-to-local pattern scaling, and the empirical temperature-growth parameters. Each type of uncertainty is calculated as the standard deviation of Da across that dimension and the mean of all other dimensions, for example, the standard deviation of Da across all FaIR carbon cycle parameters and the mean across the pattern scaling coefficients and temperature-growth parameters. We split out CMIP6 and CESM1-LE in this calculation to represent global-to-local pattern scaling uncertainty from model structure and internal variability, respectively.

We then follow previous uncertainty partitioning work by calculating additive 95% uncertainty ranges around the mean (Hawkins & Sutton, 2009; Lehner et al., 2020). The component of the 95% range due to carbon cycle uncertainty R, for example, is calculated as:

$$\pm \frac{1.96*R}{F}$$
(8)

where F is calculated as:

$$F=\frac{R+B+MS+IV}{\sqrt{{R}^{2}+{B}^{2}+M{S}^{2}+I{V}^{2}}}$$
(9)

Here B denotes temperature-growth parameter uncertainty (i.e., the standard deviation of the damage estimates across the temperature-growth parameters with the other dimensions held at their mean values), R denotes FaIR carbon cycle uncertainty, MS denotes uncertainty in the pattern scaling from model structure, and IV denotes uncertainty in the pattern scaling from internal variability. The sum of all four sources of uncertainty is shown as the range of the bars in Fig. 2b, c, with the components of total uncertainty colored by type.

Fig. 2
figure 2

source of uncertainty (Methods). Note the different y-axes in b and c

U.S.-attributable climate damages. (a) Ensemble mean GDPpc changes in each country attributable to U.S. emissions, over 1990–2014 with territorial emissions accounting and a short-run (contemporaneous) damage function. Missing data (white countries) denotes countries without continuous GDPpc data from 1990 to 2014. b, c U.S.-attributable damages in the five countries with the greatest GDPpc percent decreases (b) and percent increases (c). The black lines show the mean, the boxes denote the 95% ensemble range, and the colored portions denote the additive fraction of each 95% range due to each

4 Results

The spatial pattern of economic damages from warming attributable to U.S. emissions is shown in Fig. 2. Income changes are positive in the high latitudes, mild in the midlatitudes, and negative in the tropics (Fig. 2a), consistent with each country’s placement relative to the optimum temperature for growth (Fig. 1c). Losses are concentrated around 1–2% of GDPpc across nations in South America, Africa, and South and Southeast Asia, where temperature shocks may damage economic factors such as labor productivity (Dunne et al., 2013) and agricultural yields (Schlenker & Roberts, 2009). On the other hand, gains exceed 3–4% of GDPpc in Canada, Russia, and Scandinavia; cold baseline temperatures mean that warming makes economic output easier in these countries.

Our integrated framework allows us to quantify not just the magnitude of losses attributable to individual countries’ emissions, but also the various uncertainties in these estimates (Methods), which has not been possible in studies that only address part of the causal chain from emissions to impact. Our uncertainty decomposition reveals significant roles for the responses of the global carbon cycle and climate to emissions, the model response to global forcing as represented by CMIP6 (“GCM structure”), internal climate variability in the response to global forcing as represented by the CESM1-LE (“internal var.”) and the empirical temperature-growth parameters. In the five most-damaged countries, most of the uncertainty comes from GCM structure and internal variability (Fig. 2b). Variation in damages is therefore highly sensitive to uncertainty in the local warming that is consistent with any one emissions trajectory. Internal variability alone can play a significant role in pattern scaling uncertainty (Giorgi, 2008); this source of uncertainty has been highlighted as crucial for evaluating policy success and failure (Mankin et al., 2020). However, despite these uncertainties, the confidence intervals on the U.S. effect on these highly damaged countries do not include zero, indicating high confidence that U.S. emissions have caused harm (Fig. 2b). Uncertainty grows at higher latitudes, with major contributions from all four sources, especially GCM structural uncertainty (Fig. 2c).

Our results contrast with uncertainty decompositions for the economic costs of overall global warming, where uncertainty in the temperature-growth regression is largest (Burke et al., 2018). These previous analyses did not explicitly consider internal climate variability, which we find to be a major contributor to uncertainty in attributable damages. Our results may also differ because we consider uncertainty in country-level damages attributable to other countries, rather than uncertainty in total global losses. Uncertainty in total losses is dominated by the temperature-growth regression in part because many high-income countries are near the temperature optimum, so small changes in the optimum can flip the sign of the effect of climate change in these countries (Burke et al., 2018). However, for tropical countries who lose out in all temperature-growth parameter samples (Fig. 2b), uncertainty in the propagation of global forcing to local climate can dominate. A country-level uncertainty decomposition is important because climate liability is often addressed through the lens of local decision-makers seeking redress for a particularized harm (Burger et al., 2020), and our results demonstrate that decision-makers in tropical countries seeking this redress should focus on reducing physical climate uncertainty rather than economic uncertainty, for example.

One source of uncertainty we do not consider is uncertainty in emissions inventories, which is poorly constrained and an area of active research (Fowlie & Reguant, 2018; Hoesly et al., 2018; Höhne et al., 2011; Karstensen et al., 2015). However, emissions uncertainty is largest further back in time and in countries whose emissions are primarily LULCC-driven (Höhne et al., 2011), meaning that damages caused by the largest emitters in recent years are likely not subject to substantial emissions uncertainty (Karstensen et al., 2015). Moreover, our use of the official CMIP emissions inventory positions an appropriate analysis using data accepted by the World Climate Research Programme in support of the UNFCCC, even if that data is provisional and subject to uncertainty.

Performing this same damage attribution exercise for every emitting country and aggregating these country-level damages to global sums yields the cumulative global damages and benefits caused by the emissions of each nation (Methods, Fig. 3). We show losses and benefits attributable to each country separately. The reason for this presentation choice is simply that adding the damages and benefits attributable to any one country to calculate net damages would assume that the benefits accruing in one place from that country’s emissions cancel the losses suffered in another from those same emissions, which is not the case. Claims for climate liability require concrete and particularized injuries suffered by specific plaintiffs (Burger et al., 2020), so assessing damages and benefits independently is essential.

Fig. 3
figure 3

Global attributable climate damages and benefits. (a) Cumulative income losses attributable to each country over 1990–2014, using territorial emissions accounting and a short-run (contemporaneous) damage function. (b) Cumulative income gains attributable to each country using the same accounting choices as in (a). (c, d) Changes in each country’s attributable losses (c) and gains (d) over 1990–2014 when the accounting shifts from territorial emissions to consumption-based emissions, as countries that import (export) carbon-intensive goods and services cause more (less) damages and benefits. Changes are calculated relative to (a) and (b); that is, panel (c) shows the change in attributable damages when using consumption-based accounting relative to the attributable damages shown in (a). Positive values mean that consumption-based accounting results in greater attributable damages than territorial accounting. (e, f) Changes in each country’s attributable losses (e) and gains (f) when using territorial emissions accounting and an earlier start date of 1960 rather than 1990 as in (a)–(d). In all maps, only countries with statistically significant damages and/or benefits (Methods) are filled. Countries with statistically insignificant damages and/or benefits are shown in gray

Global income changes attributable to the U.S. and China’s emissions over 1990–2014 each exceed US$1.8 trillion in both losses and benefits; losses and benefits induced by Russia, India, and Brazil each individually exceed US$500 billion (Fig. 3a, b). The US$6 trillion in cumulative losses attributable to those five countries alone is comparable both to some 11% of annual world GDP and to the economic losses associated with warming the planet to 2 °C rather than 1.5 °C (Burke et al., 2018). Large emitters make disproportionate contributions to climate damages; the top 10 most damaging countries are together responsible for more than 67% of losses and 70% of benefits (Supplementary Fig. 6). The U.S. contributes the most, responsible for 16.5% of losses and 18% of benefits, followed by China, responsible for 15.8% of losses and 16.8% of benefits; every other country individually contributes less than 10%.

The compounding uncertainties in each step of the analysis—from emissions to GHG concentrations, from GHG concentrations to temperature changes, and from temperature changes to economic damages—imply that country-level attributable damages can be highly uncertain. To address this, our analysis incorporates a significance test to ensure that we only consider statistically significant damages (Methods). Figure 3 reflects this significance test by only aggregating and showing statistically significant damages. The compounding uncertainties involved mean that many of the smallest emitters do not cause statistically detectable effects (gray countries in Fig. 3). Large emitters, on the other hand (the 32% of countries with data in Fig. 3), tend to have damage signals substantial enough to be significant when considered against all sources of uncertainty. Our incorporation of a significance test also means that small differences in emissions can lead to major differences in attributable damages. For example, Pakistan and Bolivia have similar CO2 emissions, with Pakistan emitting an average of 35 MtC year−1 of CO2 over 1990–2014 and Bolivia emitting an average of 32 MtC year−1. But that slight difference means that while Pakistan can be tied to $130B in statistically significant losses, Bolivia cannot be tied to any. So, while the result that heavy emitters have caused significant attributable damages may seem obvious, it is not possible to perform a quantitative attribution of precisely how much damage and to whom, nor how confident we can be in such estimates, absent an integrated framework that propagates uncertainty from each step in the analysis.

Our main analysis (Fig. 3a, b) uses territorial emissions accounting from 1990 to 2014, but other accounting choices are equally valid (Skeie et al., 2017). Consumption-based accounting incorporates emissions embodied in international trade (Davis & Caldeira, 2010) and shifts the spatial pattern of attributable damages, in particular diminishing the income changes attributable to countries in the Eastern Hemisphere (Fig. 3c, d). Under consumption-based accounting, attributable damages increase by 1.5% for the U.S., increase by 10–20% for some European countries, and decrease by 15% for Russia and 9% for China (Fig. 3c), expanding the gap in responsibility between the U.S. and all other countries (Supplementary Fig. 6).

Additionally, because of nonlinearities in the responses of the climate to emissions and of the economy to the climate, the timing of emissions changes their economic effects (Höhne & Blok, 2005). Emissions released early in the twentieth century, for example, may have more severe effects because of the longer time period over which the growth effects of such emissions accumulate; at the same time, carbon sinks are stronger at lower emissions levels and the growth effects of warming are milder when baseline temperatures are cooler, so the opposite could be true. To disentangle these effects, we repeat our analysis with a start date of 1960 rather than 1990 (Fig. 3e, f). Such a shift results in larger effects attributed to all countries due to a longer period of accumulation for both GHGs and the growth effects of warming. However, the increase is not globally uniform. Developing countries’ attributable damages and benefits only increase modestly, while many countries in Europe experience > 400% increases, given their large shares of pre-1990 emissions. The sensitivity of attributable damages to the accounting period is crucial to how these results should be interpreted: damages and thus potential liabilities depend on when emissions accounting begins and on what emissions are considered. These choices are ultimately political, but our analysis provides a quantitative grounding to inform these political choices.

Our main analysis quantifies parametric uncertainty in the temperature-growth relationship, but there is also structural uncertainty in the functional form of that relationship. We test two alternative formulations: one alternative functional form which examines the lagged effects of temperature shocks to determine whether those shocks are transient or persistent (Supplementary Figs. 710) and one alternative functional form which allows the responses of high- and low-income countries to differ (Supplementary Figs. 11 and 12). Including lagged responses increases the fraction of uncertainty attributable to the damage function, but the magnitude of damages and their statistical significance are similar to our main analysis (Supplementary Material). Differentiating the responses of high- and low-income countries increases both the damages and benefits attributable to major emitters (Supplementary Material). These analyses demonstrate that our results are generally robust to uncertainty in the functional form of the temperature-growth relationship.

The spatial pattern of attributable damages (e.g., Fig. 2a) has clear distributional implications. Countries that lose income are warmer and poorer than the global average, generally located in the tropics. Countries that gain income are cooler and wealthier than the global average, generally located in the midlatitudes. Global warming to date has amplified, and will continue to amplify, this extant pattern of global economic inequality (Diffenbaugh & Burke, 2019).

Our results provide an additional dimension to this globally unequal pattern: The cool, relatively wealthy countries that have gained from anthropogenic warming are also those that (1) have emitted the most and (2) caused the most damage to other countries from their emissions (Fig. 4). Nearly all the high-emitting nations in North America and Eurasia are in the top two GDP per capita income quintiles over 1990–2014, though China, India, and Indonesia are exceptions (Fig. 4, inset map). These top several income quintiles have caused income losses in the poorest two quintiles, while they have caused income gains for themselves that exceed those losses in magnitude (Fig. 4). On the other hand, countries in the lowest income quintile, primarily in Africa and central and south Asia, have caused nearly zero effects on other countries while suffering the greatest disadvantages from the emissions of larger economies. While uncertainty is large for the gains experienced by higher income quintiles (see also Fig. 2c), uncertainty is relatively low for damages in the bottom two quintiles. As such, the distributional picture of culpability and equity in climate damages is evident, despite the myriad complexities involved in linking individual nations’ emissions to climate impacts.

Fig. 4
figure 4

Income distribution of damages caused and experienced. Bar heights show the average GDPpc change experienced by countries in each income quintile that are attributable to the average emitter in each income quintile (colors). Error bars show + / − one standard deviation of the mean across the distribution of pattern scaling coefficients, FaIR realizations, and damage function parameters (Methods). Inset map shows each country’s income quintile, calculated using 1990–2014 average GDPpc

The results shown in Fig. 4 use territorial emissions accounting over 1990–2014 as in our baseline analysis, but many high-income, high-emitting countries have also imported additional emissions through their demand for products from the developing world (Davis & Caldeira, 2010). Shifting to consumption-based emissions accounting would therefore amplify the distributional picture illustrated in Fig. 4 (primarily in Europe; cf. Figure 4, Fig. 3c, d), as would shifting to an earlier accounting period (Fig. 3e, f).

These results are consistent with previous work that shows increases in global inequality from historical warming (Diffenbaugh & Burke, 2019). However, our findings emphasize that the culpability for warming rests primarily with a handful of major emitters, and that this warming has resulted in the emitters’ enrichment at the expense of the poorest people in the world: we quantify those costs to individual countries and who precisely is responsible for them. Our work shows that anthropogenic warming constitutes a substantial international wealth transfer from the poor to the wealthy.

5 Discussion

Our analysis has shown that GHG emissions from high-emitting countries have caused substantial economic losses in low-income, tropical parts of the world and economic gains in high-income, midlatitude regions. Critically, these economic changes are attributable to the largest emitters despite the substantial uncertainties at each step in the causal chain from emissions to impact. Our results are robust despite the wide range of carbon cycle parameters and climate sensitivities, global-to-local forcing strengths, and temperature-growth specifications we test.

These results have two key implications. Firstly, they illustrate that physical climate uncertainty may constitute the dominant source of uncertainty in losses in tropical countries that are attributable to major emitters. While uncertainty in the relationship between the climate and economy is the dominant uncertainty in global losses from warming (Burke et al., 2018), our results demonstrate that this does not hold at the country level. In the low-income tropical countries that are most vulnerable to warming, internal climate variability and differences in model structure can produce a wide range in damages attributable to major emitters like the U.S. Scientific efforts to narrow uncertainty in regional climate change may therefore pay large dividends for countries seeking legal recourse for climate damages.

Secondly, our results show that the actions of specific emitters can be tied to the downstream monetary implications of climate change. Emerging discussions about climate liability have been limited to date by a lack of scientific evidence supporting causal linkages between individual countries’ emissions and the consequent local impacts (Burger et al., 2020; Stuart-Smith et al., 2021). Our framework shows that such linkages can be quantified using state-of-the-art climate models and empirical approaches and that we can process-trace exactly who has caused economic losses from their emissions, and how much. While previous studies have illustrated the economic harms of global warming, our work shows that these harms can be assigned to individual emitters in a way that rigorously accounts for the compounding uncertainties at each step of the causal chain from emissions to local impact. Finally, it is worth noting that our approach can be generalized to other actors, such as individual firms (Ekwurzel et al., 2017; Heede, 2014; Licker et al., 2019), or to other harms, such as the economic losses suffered by farmers due to extreme heat (Diffenbaugh et al., 2021). These results therefore contribute to resolving a key barrier to climate liability efforts and advance these critical emerging discussions.