Inequitable distribution of risks associated with occupational heat exposure driven by trade

Analytical framework

This study uses a multi-regional input-output model with a labour satellite account to estimate labour hours embodied in the global production networks, which is further combined with the gridded population to obtain grid-level labour in production networks. By combining this with a climate model that provides grid-level heat exposure, as well as with other data such as average working hours, labour’s exposure to heat stress can be accounted for. How these modules are integrated is illustrated in Supplementary Information (Supplementary Figs. 1, 2).

Labour embodied in trade

The input-output (IO) model has been widely adopted to investigate the interconnections between different sectors. By applying multi-regional input-output (MRIO) models, researchers can uncover how countries/regions are connected bilaterally through trade. These MRIO models can be extended with satellite accounts to trace emissions, resources, pollution, water, land, as well as many other aspects transferred through the international production system^42,64,65. Here, we use a global multi-regional input-output model to trace labour hours embodied in trade. Assuming that there are \(m\) regions and \(n\) sectors involved in international trade, the equilibrium between product supply and demand can be expressed in Eq. (1).

where, \({{\bf{A}}}\) is a direct input coefficient matrix, \({{\bf{X}}}\) is the total output matrix, and \({{\bf{Y}}}\) is the final demand matrix. The Leontief model links final demand with gross outputs, which is shown in Eq. (2).

where \({{\bf{B}}}\) is the Leontief inverse matrix. Denote the labour intensity matrix by \({{\bf{E}}}\), which is a diagonalized matrix representing labour employment per unit of gross output for each region-sector, we can get Eq. (3).

\({{{\bf{L}}}}_{{mn}\times {mn}}\,\) is the global labour employed, in which its element \({l}_{{isjr}}\) is the labour employment in region \(i\) sector \(s\) created by region \(j\) sector \(r\). Thus, by summing the matrix \({{\bf{L}}}\) by rows, we obtain production side labour employment of each region-sector; summing the matrix \({{\bf{L}}}\) by columns, we obtain consumption side labour induced by each region-sector.

Matching labour employment data to MRIO

There are 26 sectors in the Eora multiregional input-output tables, and 14 economic activities in the International Labour Organisation (ILO) labour employment database. To match the labour employment categories with Eora economic sectors we follow a widely adopted approach, namely matching these two data sources by using a concordance table (Supplementary Table 2).

Measuring extreme heat exposure hours

To measure labour force exposure hours to extreme heat, we combine climate models with socio-economic information. The climate model is based on the ERA5⁶⁶, from which researchers have developed the universal thermal climate index (UTCI)³⁶. UTCI is a human biometeorology parameter—an equivalent temperature (°C)—that measures human physiological response to the thermal environment. By considering how the human body experiences atmospheric conditions (air temperature, humidity, wind, and radiation), the UTCI describes synergistic heat exchanges between the thermal environment and the human body. Four variables from the ERA5 are necessary to build the UTCI index, including 2 m air temperature, 2 m dew point temperature (or relative humidity), wind speed at 10 m above the ground level and mean radiant temperature (MRT). There are 10 UTCI thermal stress categories that correspond to specific human physiological responses to the thermal environment. The categories related to UTCI degree Celsius (°C) are as follows: above +46: extreme heat stress; + 38 to + 46: very strong heat stress; + 32 to + 38: strong heat stress; + 26 to + 32: moderate heat stress; + 9 to + 26: no thermal stress; + 9 to 0: slight cold stress; 0 to − 13: moderate cold stress; − 13 to − 27: strong cold stress; − 27 to − 40: very strong cold stress; below − 40: extreme cold stress (Supplementary Table 3). Thus, working hours each day can be classified into different heat stress levels (Supplementary Figs. 35–222). Defining the exposure to heat exposure as the number of hours exposed to extreme/very strong/strong heat stress, the UTCI index can be converted to a binary heat exposure indicator, which is shown in Eq. (4).

where \({{{\rm{UTCI}}}}_{{gydh}}^{c}\) is the converted UTCI index for the gth grid on hour \(h\) in day \(d\) of year \(y\), by subtracting 273.15 from the original UTCI index stored in Kelvin (K). Then we make two assumptions. First, following the idea of 12 h exposure time in previous studies, we set the maximum daily working hour to be 12 h, which is equivalent to 2400 annual working hours for 200 work days/year and 3000 annual working hours for 250 work days/year. Second, as different grids may differ in annual work hours and average daily work hours, we assume that the annual work hours are evenly distributed throughout the 12 h available work time per day. Then, the daily and annual heat exposure cap can be summed by using Eq. (5).

Thus, the average exposed work hours to heat can be obtained by using Eq. (6).

where \({\omega }_{{gy}}\) is the average work hours of the labour force in gth grid for year \(y\), which is obtained from region-specific database and estimation.

Linking grid cell labour embodied with grid cell exposure

We use a population density module to link the interregional embodied labour flow with grid cell annual average per capita work hours exposed to extreme heat. Thus, the labour force within each grid cell can be linked from a production network perspective, which is shown in Eq. (7).

where \({p}_{{gy}}\) is the population density of the gth grid in year \(y\), \({a}_{g}\) is the area of the gth grid, and \({G}_{i}\) is all grid cells in country/region \(i\). Then the labour force exposure to extreme temperatures (LET) linked with the production network on the grid cell level and country/region level can be expressed in Eqs. (8, 9).

Here, \({{{\rm{LET}}}}_{{isjr}}\) represents the embodied heat exposure of the labour force, \({l}_{{isjr}}\) represents the embodied labour employment, and \(\frac{{\sum }_{g\in {G}_{i}}{w}_{{gy}}\times {p}_{{gy}}{\times a}_{g}}{{\sum }_{g\in {G}_{i}}{p}_{{gy}}\times {a}_{g}}\) reflects the average per capita exposure. Summing \({{\rm{LET}}}\) by \(i\) and \(s\), we obtain the total heat exposure induced by sector \(r\) in region \(j\), which is a consumption side exposure; while summing \({{\rm{LET}}}\) by \(j\) and \(r\), we obtain the total heat exposure induced by sector \(s\) in region \(i\), which is a production side exposure.

Structural decomposition analysis

We use an SDA to attribute the role of each driving factor in the total changes in labour force heat exposure (Supplementary Method 1). Based on Eqs. (3) and (9), we quantify \({{\bf{LET}}}\) based on the input-output framework in a matrix form,

where \({{\bf{H}}}\) is the diagonal matrix for average heat exposure index per labour employed in each country, \({{\bf{E}}}\) is a diagonal matrix for labour intensity per gross output for each country-sector, \({{\bf{B}}}\) is the Leontief inverse matrix representing the production structure, \({{{\bf{Y}}}}_{s}\) is a \({mn}\times m\) matrix with the sum of all elements equal to 1, representing the final demand structure, and \({{{\bf{Y}}}}_{t}\) is the diagonal matrix with all diagonal elements as the same value—the deflated global total demand. Thus, the changes \({{\bf{LET}}}\) can be attributable to its driving factors, including the climate change effect \(\Delta {{\bf{H}}}\), the labour intensity effect \(\Delta {{\bf{E}}}\), the production structure effect \(\Delta {{\bf{B}}}={\Delta {{\bf{B}}}}^{d}+{\Delta {{\bf{B}}}}^{e}\), the final demand structure effect \(\Delta {{{\bf{Y}}}}_{s}=\Delta {{{{\bf{Y}}}}_{s}}^{d}+{{{{\bf{Y}}}}_{s}}^{e}\), and the total final demand growth effect \(\Delta {Y}_{t}\),

Then, we apply the SDA method on constant price input-output tables obtained by taking the Producer Price Indices (PPIs) as deflators and the Logarithmic Mean Divisia Index (LMDI) approach to estimate the role of each driving factor between time periods (Supplementary Fig. 7).

Robustness analysis

To test the robustness of the results, this study further conducts several sensitivity analyses. First, a comparison between this study and previous studies was conducted (Supplementary Table 4). Second, we use alternative data sources by replacing the input-output table from Eora with Inter-Country Input-Output tables compiled by the Organisation for Economic Co-operation and Development (OECD-ICIO), which consists of 67 economies and 45 sectors from 1995 to 2018. The results drawn from these two data sources are compared on the global and national levels both on the production and consumption sides from 1995 to 2015 (Supplementary Figs. 223–225). Third, we also use Monte Carlo simulation to test the robustness of the original heat exposure estimation, based on whether it falls within the expected range of values derived from the simulations. We conduct 1000 Monte Carlo simulations, where the heat stress drawn from the ERA5 is set to a normal distribution with a mean equal to the baseline values and a standard deviation set to 10% of each respective baseline value. The results indicate that the original estimation falls well within the range of expected values generated by these simulations, and the global labour force’s total heat exposure remain robust from 1995 to 2020 when heat stress varies within an acceptable level (Supplementary Figs. 226–228).

Data

Multiregional Input-output tables. To measure trade flows between countries/regions, we need to use information from MRIO tables. There are several available MRIO databases, including Eora⁶⁷, World Input-Output Database (WIOD)⁶⁸, OECD-ICIO⁶⁹, and Exiobase⁷⁰, as well as other newly emerging databases. In this study, the global MRIO table is collected from the Eora global supply chain database. Eora provides high-resolution multiregional IO tables with matching environmental and social satellite accounts for 190 economies and 26 sectors. Its full geographical coverage, especially with detailed information on tropical and low-income economies, allows all the global labour employment to be analysed. In addition, it has a long temporal coverage, which helps uncover the trends and compare them across different time periods.

Labour employment. Labour employment data were drawn from the International Labour Organisation (ILO) via ILOSTAT explorer. The ILO-modelled estimates provide a complete set of internationally comparable labour statistics, which is a balanced panel data set with consistent country/region coverage. The data is mainly based on nationally reported observations, and the missing data is imputed using the ILO model. ILO evaluates existing self-reported data, selects only those observations deemed sufficiently comparable across countries/regions and runs the model to obtain the ILO-modelled estimates. For the working-age population who is at least 15 years old, the labour employment data from ILO provides a breakdown of total employment by economic sector and gender. The economic activity classification is based on the International Standard Industrial Classification (ISIC) Rev.4, which includes 14 sectors.

Working Hours. The average number of annual working hours for the labour force in each country/region is originally obtained from Penn World Table version 10.01. Penn World Table is a database with information on relative levels of income, output, input and productivity, covering 183 countries/regions between 1950 and 2019⁴¹. We use variable avh in the database—average annual work hours by labour engaged—as the work hours for each country/region. For those countries that have missing data for a given year, we use the random forest to fill in the missing values (Supplementary Method 2 and Supplementary Fig. 229). For the year 2020, the relevant data from 2019 was applied in the analysis.

Gridded population density data. To capture heterogeneous levels of exposure to extreme heat geographically, we need to split those country-level employment data to the grid level. To do this, we use the population density in each grid from the fourth version of the Gridded Population of the World data (GPWv4). In GPWv4, population input data were collected at the detailed spatial resolution available from the results of the 2010 round of Population and Housing Censuses and then extrapolated to produce population estimates for years 2000, 2005, 2010, 2015, and 2020. For the year 1995, which is not covered by GPWv4, the population density of the year 2000 is used instead. The population density raster data sets are gridded with an output resolution of 30 arc-seconds (~ 1 km at the equator).

Climate data. To obtain gridded extreme heat exposure, we draw climate data from the ECMWF ERA5 reanalysis. The ERA5 provides UTCI—an equivalent temperature that measures human physiological response to the thermal environment. To obtain the daily extreme heat hours in each grid, we separate the UTCI index by categories: (1) above + 46: extreme heat stress; (2) + 38 to + 46: very strong heat stress; (3) + 32 to + 38: strong heat stress; (4) + 26 to + 32: moderate heat stress; (5) + 9 to + 26: no thermal stress; (6) + 9 to 0: slight cold stress; (7) 0 to − 13: moderate cold stress; (8) − 13 to − 27: (9) strong cold stress; (10) − 27 to − 40: very strong cold stress; below − 40: extreme cold stress.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.