Global supply chains amplify economic costs of future extreme heat risk – Nature

Our methodology, in essence, combines three modules of climate, health and economy with full validation (Extended Data Fig. 1). The integrated model links climate module (estimating future climate parameters including surface air temperature and relative humidity and so on), demographic and health module (simulating future world population dynamics and exposure–response functions to warming) and economic module (dynamic footprint of heat-induced labour loss on global economy and supply chain).

Climate module

Fourteen GCMs involved in the framework of CMIP6 (Extended Data Table 1) with ten bias-corrected models from ISIMIP3b^49,50 are used to estimate the modelled heat stress projection for the end of the twenty-first century. Five models were randomly averaged several times from the climate model ensemble as a Monte Carlo uncertainty analysis. ERA5 re-analysis data⁵¹ from 1985 to 2022 are used for bias-correction and validation. Climatic parameters such as maximum and average temperature and relative humidity on a daily scale are integrated, which are closely related to future working environment (Supplementary Fig. 9).

Many institutes, including International Standards Organization (ISO) and US National Institute for Occupational Safety and Health (NIOSH), use WBGT to quantify different amounts of heat stress and define the percentage of a typical working hour that a person can work while maintaining core body temperature. To facilitate the long-term calculation, we use¹⁸ simplified WBGT, which approximates WBGT well using temperature (T_a) and relative humidity (RH)^52,53 as parameters such as solar radiation and wind speed have higher uncertainty and weaker effects at the global scale. To take into account indoor heat exposures for industrial and service sector workers, we used the approximation that indoor WGBT_indoor = WBGT_outdoor − 4, based on a deduction of the radiation exposure factor from the formula below¹⁸:

We also calculated the spatial and temporal evolutionary trends in the occurrence of future heatwaves to calculate excess mortality. There is no consistent definition for heatwave worldwide because people may have acclimatized to their local climatic zones and different studies have applied various temperature metrics^54,55. Heatwaves are usually defined by absolute or relative temperature threshold in consecutive days⁵⁶. There are various ways to define a heatwave. For example, the IPCC defines heatwave as “a period of abnormally hot weather, often defined with reference to a relative temperature threshold, lasting from two days to months”, whereas the Chinese Meteorological Administration defined heatwave as “at least three consecutive days with maximum temperature exceeding 35 °C”. Others³¹ identified heatwave using the TX90p criterion, that is, when the 90th percentile of the distribution of regional maximum temperatures spanned by data from the period 1981–2010 was exceeded for at least three consecutive days. In our study, two or more consecutive days above the 95% threshold of the 1985–2015 ERA5 daily mean temperature^51,57 were defined as a heatwave, which is considered to be a moderate estimation and is widely used in epidemiological studies^36,58,59. Several definitions, such as four or more consecutive days above the 97.5% threshold, are used as sensitivity analysis. Considering certain amounts of climate adaptation of the local resident along the warming climate, dynamic heatwave thresholds⁶⁰ are defined as part of the uncertainty analysis in this study; that is two or more consecutive days above the 95% threshold of the daily mean temperature between 1985 and the year before the target year were defined as a heatwave (ERA5 data are used for 1985–2014; climate projection data are used after 2015). The use of a dynamic threshold based on both historical and climate projections data helps to incorporate the human adaptation of heat stress in a long-term warming scenario, as reported in recent studies^61,62,63,64.

Health costs related to heat exposure

Some studies have shown that the health impact of heatwaves could vary substantially with location^65,66. Few studies have investigated the heatwave-induced mortality risk at a global scale^41,67. A primitive health risk function associating heatwave mortality risks with four different climate zones was established by ref. ³⁶ on the basis of a comprehensive study using data from 400 communities in 18 countries/regions across several years (1972–2012). Here, we used the relative risk coefficients (Extended Data Table 2) from figure 4 of ref. ³⁶ for four different climate zones (Extended Data Fig. 3) to estimate potential heatwave-related death due to climate change on a global scale. The simplified four-climate-zone-based estimation may neglect subregional characters and should be interpreted with caution, as further factors affecting heat-induced death (such as air condition accessibility⁶⁸, age^69,70,71,72 and humidity⁷³) are not included in this study.

The number of excess deaths D_hw during a heatwave period was calculated at each grid cell level (0.5°) with the following equation:

POP is the population at the given location consistent with the SSPs⁷⁴. MR is the average daily mortality rate (2009–2019) at the country level obtained from the World Bank⁷⁵. For 37 countries with large territory and more refined data (for example, European Union (including UK), Russia, Ukraine, China, the USA, Canada, Brazil, South Africa, India and Australia), we used state/provincial statistics based on data from national statistical offices (Source, World Bank; state/province level data for European Union, Eurostat⁷⁶; Russia, The Russian Fertility and Mortality database⁷⁷; China, China Statistical Yearbook 2019⁷⁸; the USA, National Institutes of Health⁷⁹; Brazil, Fundação Amazônia de Amparo a Estudos e Pesquisas⁸⁰; Canada, Statistics Canada⁸¹; Australia, Australian Bureau of Statistics⁸²; India, Ministry of Finance Economic Survey⁸³). RR is the relative risk of mortality caused by heatwaves. HWN is the number of heatwave days for the given year and location (Extended Data Fig. 2).

The calculated excess deaths are translated to a social-economic loss on the basis of the value of statistical life (VSL). The concept of VSL is widely used throughout the world to monetize fatality risks in benefit–cost analyses. The VSL represents the individual’s local money–mortality risk tradeoff value, which is the value of small changes in risk, not the value attached to identified lives. The country-based VSL estimation used in this research is adopted from the global health risks pricing study by ref. ⁸⁴. The estimation is based on the estimated VSL in the USA (US$₂₀₁₉11 million) and coupled with an income elasticity of 1.0 to adjust the VSL to other countries using the fixed-effects specification. A similar health valuation method has been adopted in past studies^85,86 and was recommended in the report of the World Bank⁸⁷. Moreover, a sensitivity test is conducted under the assumption that all life would be valued equally across the world (Supplementary Figs. 2 and 3). For such a test, an averaged VSL is calculated by summing up each country’s income-based VSL times its population then dividing by the total population of the world.

Expose function of labour productivity

The increase in daily temperatures affects the efficiency of workers and reduces safe working time. A compromise in endurance capacity due to thermoregulatory stress was already evident at 21 °C. Different studies used similar methods to evaluate the labour loss function. The form of logistic function with ‘S’ shape has become the consensus of the academic community but the specific functional equation and parameters are various in different studies. The loss functions used in mainstream research include exponential function⁸⁸ as equation (4), cumulative normal distribution function^5,41 as equation (5) and so on. In this research, we adopt the cumulative normal distribution function (equation (5)) as our benchmark function because it was extensively applied and case proven in 3-year reports of the Lancet Countdown on health and climate change^5,41,89,90. Because the Hothaps function (equation (4)) is subject to parameter uncertainty as a result of being based on a few empirical studies, we use it to test for the sensitivity of our estimates (Supplementary Figs. 4 and 5). Our methodology identifies three ISO standard work intensity amounts: 200 W (assumed to be office workers in the service industry, engaged in light work indoors), 300 W (assumed to be industrial workers, engaged in moderate work indoors) and 400 W (assumed to be construction or agricultural workers, engaged in heavy work outside). For example, to calculate workability loss fraction in India’s food production sector (300 W, indoor), we bring the corresponding parameters (Extended Data Table 3) and WBGT_indoor into equation (5). Previous studies have tended to ignore indoor workforce loss, assuming that the indoor workforce was very low under current climate condition or protected by air conditioning⁹¹. However, a growing number of studies have proved that future indoor labour losses cannot be underestimated³¹. For example, only 7% of households in India possess an air conditioner, despite having extremely high cooling needs. Considering the severe adaptation cooling deficit in emerging economies⁹², indoor labour losses must be fully considered in global-scale studies. This study uses the climate–income–air conditioner usage function published by ref. ⁹³ to assess the rate of air conditioning protection in conjunction with the per capita income of each country under each SSP scenario. Higher per capita income in each country leads to higher air-conditioning penetration, whereas the climate base determines the rate and trend of increase in air-conditioning penetration (elasticity of penetration to income). In our study, we improved the function by replacing cooling degree days (CDDs) with indoor WBGT, as CDDs only consider temperature neglecting humidity. Only the indoor workforce under air conditioning, will be protected from heat-induced loss.

Of which the parameters for a given activity level (Prod_mean and Prod_SD, defined as the amount of internal heat generated in performing the activity) are given in Extended Data Table 3, and ERF is the error function defined as:

To calculate average daily impacts, we use an approximation for hourly data based on the 4 + 4 + 4 method implemented by ref. ¹⁴. We assume that 4 h per day is close to WBGT_max and 4 h per day is close to WBGT_mean (early morning and early evening). The remaining 4 h of a 12 h daylight day is assumed to be halfway between WBGT_mean and WBGT_max (labelled WBGT_half). The analysis above gives the summer daily potential workability lost in each grid cell at each amount of work intensity and environment (200–400 W, indoor or outdoor). By combining this with the dynamic population grid under each SSP scenario (see Supplementary Fig. 13 for comparison with static population setting), we aggregate to obtain country-scale labour productivity losses. In the disaster footprint model, we adopt the approach presented by ref. ⁵ which defines the timeframe for computing labour productivity losses as the warm season (June to 30 September in the Northern Hemisphere and December to 30 March in the Southern Hemisphere) to adjust the overestimation of the risk of moderate hot temperature, as the model is more applicable to sudden and strong shocks rather than moderate changes throughout the year.

Global disaster footprint analysis module

The global economic loss will be calculated using the following hybrid input–output and computable general equilibrium (CGE) global trade module. Our global trade module is an extension of the adaptive regional input–output (ARIO) model^20,94,95, which was widely used in the literature to simulate the propagation of negative shocks throughout the economy^96,97,98,99. Our model improves the ARIO model in two ways. The first improvement is related to the substitutability of products from the same sector sourced from different regions. Second, in our model, clients will choose their suppliers across regions on the basis of their capacity. These two improvements contribute to a more realistic representation of bottlenecks along global supply chains¹⁰⁰.

Our global trade module mainly includes four modules: production module, allocation module, demand module and simulation module. The production module is mainly designed for characterizing the firm’s production activities. The allocation module is mainly used to describe how firms allocate output to their clients, including downstream firms (intermediate demand) and households (final demand). The demand module is mainly used to describe how clients place orders to their suppliers. And the simulation module is mainly designed for executing the whole simulation procedure.

Production module

The production module is used to characterize production processes. Firms rent capital and use labour to process natural resources and intermediate inputs produced by other firms into a specific product. The production process for firm i can be expressed as follows,

where x_i denotes the output of the firm i, in monetary value; p denotes type of intermediate products; \({z}_{i}^{{\rm{p}}}\) denotes intermediate products used in production processes; va_i denotes the primary inputs to production, such as labour (L), capital (K) and natural resources (NR). The production function for firms is f(·). There is a wide range of functional forms, such as Leontief¹⁰¹, Cobb–Douglas and constant elasticity of substitution production function¹⁰². Different functional forms reflect the possibility for firms to substitute an input for another. Considering that heat stress tends to be concentrated in a specific short period of time, during which economic agents cannot easily replace inputs as suitable substitutes, might temporarily be unavailable, we use Leontief production function which does not allow substitution between inputs.

where \({a}_{i}^{{\rm{p}}}\) and \({b}_{i}\) are the input coefficients calculated as

and

where the horizontal bar indicates the value of that variable in the equilibrium state. In an equilibrium state, producers use intermediate products and primary inputs to produce goods and services to satisfy demand from their clients. After a disaster, output will decline. From a production perspective, there are mainly the following constraints.

Labour supply constraints

Labour constraints during heat stress or after a disaster may impose severe knock-on effects on the rest of the economy^21,103. This makes labour constraints a key factor to consider in disaster impact analysis. For example, in the case of heat stress, these constraints can arise from employees’ inability to work as a result of illness or extreme environmental temperatures beyond health threshold. In this model, the proportion of surviving productive capacity from the constrained labour productive capacity (\({x}_{i}^{{\rm{L}}}\)) after a shock is defined as:

Where \({\gamma }_{i}^{{\rm{L}}}(t)\) is the proportion of labour that is unavailable at each time step t during heat stress; \((1-{\gamma }_{i}^{{\rm{L}}}(t))\) contains the available proportion of employment at time t.

The proportion of the available productive capacity of labour is thus a function of the losses from the sectoral labour forces and its predisaster employment level. Following the assumption of the fixed proportion of production functions, the productive capacity of labour in each region after a disaster (\({x}_{i}^{{\rm{L}}}\)) will represent a linear proportion of the available labour capacity at each time step. Take heatwaves as an example; during extreme heatwaves that last for days on end, governments and businesses often shut down work to reduce the risk of serious illnesses such as pyrexia. This imposes an exogenous negative shock on the economic network.

Constraints on productive capital

Similar to labour constraints, the productive capacity of industrial capital in each region during the aftermath of a disaster (\({x}_{i}^{{\rm{K}}}\)) will be constrained by the surviving capacity of the industrial capital^{30,96,104,105,106}. The share of damage to each sector is directly considered as the proportion of the monetized damage to capital assets in relation to the total value of industrial capital for each sector, which is disclosed in the event account vector for each region \(({\gamma }_{i}^{{\rm{K}}})\), following ref. ¹⁰⁷. This assumption is embodied in the essence of the input–output model, which is hard-coded through the Leontief-type production function and its restricted substitution. As capital and labour are considered perfectly complementary as well as the main production factors and the full employment of those factors in the economy is also assumed, we assume that damage in capital assets is directly related with production level and, therefore, VA level. Then, the remaining productive capacity of the industrial capital at each time step is defined as:

Where, \({\bar{K}}_{i}\) is the capital stock of firm \(i\) in the predisaster situation and K_i(t) is the surviving capital stock of firm \(i\) at time \(t\) during the recovery process

Supply constraints

Firms will purchase intermediate products from their supplier in each period. Insufficient inventory of a firm’s intermediate products will create a bottleneck for production activities. The potential production level that the inventory of the pth intermediate product can support is

where \({S}_{i}^{{\rm{p}}}(t-1)\) refers to the amount of pth intermediate products held by firm i at the end of time step t − 1.

Considering all the limitation mentioned above, the maximum supply capacity of firm i can be expressed as

The actual production of firm i, \({x}_{i}^{{\rm{a}}}(t)\), depends on both its maximum supply capacity and the total orders the firm received from its clients, \({{\rm{TD}}}_{i}(t-1)\) (see section on the ‘Demand module’),

The inventory held by firm i will be consumed during the production process,

Allocation module

The allocation module mainly describes how suppliers allocate products to their clients. When some firms in the economic system suffer a negative shock, their production will be constrained by a shortage to primary inputs such as a shortage of labour supply during extreme heat stress. In this case, a firm’s output will not be able to fill all orders of its clients. A rationing scheme that reflects a mechanism on the basis of which a firm allocates an insufficient amount of products to its clients is needed¹⁰⁸. For this case study, we applied a proportional rationing scheme according to which a firm allocates its output in proportion to its orders. Under the proportional rationing scheme, the amounts of products of firm i allocated to firm j, \({{\rm{F}}{\rm{R}}{\rm{C}}}_{j}^{i}\) and household h, \({{\rm{H}}{\rm{R}}{\rm{C}}}_{h}^{i}\) are as follows,

where \({{\rm{F}}{\rm{O}}{\rm{D}}}_{i}^{j}(t-1)\) refers to the order issued by firm j to its supplier i in time step t − 1, and \({{\rm{H}}{\rm{O}}{\rm{D}}}_{i}^{h}(t-1)\) refers to the order issued by household h to its supplier j. Firm j received intermediates to restore its inventories,

Therefore, the amount of intermediate p held by firm i at the end of period t is

Demand module

The demand module represents a characterization of how firms and households issues orders to their suppliers at the end of each period. A firm orders its supplier because of the need to restore its intermediate product inventory. We assume that each firm has a specific target inventory level based on its maximum supply capacity in each time step,

Then the order issued by firm i to its supplier j is

Households issue orders to their suppliers on the basis of their demand and the supply capacity of their suppliers. In this study, the demand of household h to final products q, \({{\rm{HD}}}_{h}^{q}\left(t\right)\), is given exogenously at each time step. Then, the order issued by household (HOD) h to its supplier j is

The total order received (TOD) by firm j is

Simulation module

At each time step, the actions of firms and households are as follows in Monte Carlo simulations.

Firms plan and execute their production on the basis of three factors: (1) inventories of intermediate products they have, (2) supply of primary inputs and (3) orders from their clients. Firms will maximize their output under these constraints.

Product allocation

Firms allocate outputs to clients on the basis of their orders. In equilibrium, the output of firms just meets all orders. When production is constrained by exogenous negative shocks, outputs may not cover all orders. In this case, we use a proportional rationing scheme proposed in the literature^20,108 (see section on ‘Allocation module’) to allocate products of firms.

Firms and households issue orders to their suppliers for the next time step. Firms place orders with their suppliers on the basis of the gaps in their inventories (target inventory level minus existing inventory level). Households place orders with their suppliers on the basis of their demand. When a product comes from several suppliers, the allocation of orders is adjusted according to the production capacity of each supplier.

This discrete-time dynamic procedure can reproduce the equilibrium of the economic system and can simulate the propagation of exogenous shocks, both from firm and household side or transportation disruptions, in the economic network. From the firm side, if the supply of a firm’s primary inputs is constrained, it will have two effects. On the one hand, the decline in output in this firm means that its clients’ orders cannot be fulfilled. This will result in a decrease in inventory of these clients, which will constrain their production. This is the so-called forward or downstream effect. On the other hand, less output in this firm also means less use of intermediate products from its suppliers. This will reduce the number of orders it places on its suppliers, which will further reduce the production level of its suppliers. This is the so-called backward or upstream effect. From the household side, the fluctuation of household demand caused by exogenous shocks will also trigger the aforementioned backward effect. Take tourism as an example, when the temperature is well beyond the comfort range of the visitor, the demand for tourism from households all over the world will decline significantly. This influence will further propagate to the accommodation and catering industry through supplier–client links.

Economic footprint

We define the VA decrease of all firms in a network caused by an exogenous negative shock as the disaster footprint of the shock. For the firm directly affected by exogenous negative shocks, its loss includes two parts: (1) the VA decrease caused by exogenous constraints and (2) the VA decrease caused by propagation. The former is the direct loss, whereas the latter is the indirect loss. A negative shock’s total economic footprint (TEF_i,r), direct economic footprint (DEF_i,r) and propagated economic footprint (PEF_i,r) for firm i in region r are,

and,

and,

Global supply-chain network

We build a global supply-chain network based on v.10 of the Global Trade Analysis Project (GTAP) database¹⁰⁹ and use GTAP 9 (ref. ¹¹⁰), EMERGING database¹¹¹ for robustness analysis. GTAP 10 provides a multiregional input–output (MRIO) table for the year 2014. Also, the database for the year 2011 was used for robustness testing. This MRIO table divides the world into 141 economies, each of which contains 65 production sectors (Supplementary Tables 4 and 5). If we treat each sector as a firm (producer) and assume that each region has a representative household, we can obtain the following information in the MRIO table: (1) suppliers and clients of each firm; (2) suppliers for each household and (3) the flow of each supplier–client connection under the equilibrium condition. This provides a benchmark for our model. We also used a dynamic CGE model consistent with the SSP scenarios for a parallel assessment and as part of the robustness check of the ARIO results. Specifically, the CGE model we used is a G-RDEM¹¹² with aggregated ten regions and ten sectors^113,114,115 (Supplementary Information section 1.3).

When applying such a realistic and aggregated network to the disaster footprint model, we need to consider the substitutability of intermediate products supplied by suppliers from the same sector in different regions^115,116,117. The substitution between some intermediate products is straightforward. For example, for a firm that extracts spices from bananas it does not make much of a difference if the bananas are sourced from the Philippines or Thailand. However, for a car manufacturing firm in Japan, which uses screws from Chinese auto parts suppliers and engines from German auto parts suppliers to assemble cars, the products of the suppliers in these two regions are non-substitutable. If we assume that all goods are non-substitutable as in the traditional input–output model, then we will overestimate the loss of producers such as the case of the fragrance extraction firm. If we assume that products from suppliers in the same sector can be completely substitutable, then we will substantially underestimate the losses of producers such as the Japanese car manufacturing firm. To alleviate these shortcomings in the evaluation of losses under the two assumptions, we allow for the possibility of substitution for each sector depending on the region and sector of the supplier (Supplementary Information section 1.3).

Nonetheless, our estimates of economic damages from heat stress are subject to some important uncertainties¹¹⁸ and our methods may not capture all types of economic damages. We only include economic losses caused by heat stress on human activities without considering the impacts on infrastructure, crop growth and other factors. Considering the challenges of predicting changes to socioeconomic systems globally, we have followed the approach from the literature^23,31,91,119 to simulate supply-chain indirect losses by considering the impact of future climate risks on current socioeconomic settings. We have not considered the potential substitution of labour with capital resulting from technological advances, such as mechanization. Our analysis ignores the different levels of trade openness and globalization among SSP narratives, as well as the role of dynamic factors such as technology and price. Again, although we have conducted robustness tests for different degrees of trade substitutability, the relevant parameter is set randomly in the Monte Carlo simulation rather than derived through a general equilibrium model. The results should therefore be interpreted with caution as indicating potential future climate change risks to the existing economy rather than as quantitative predictions, given that the static representation of the economic structure in our model inevitably skews the assessment in the long run.