Escalating environmental summer heat exposure—a future threat for the European workforce

Heat exposure constitutes a major threat for European workers, with significant impacts on the workers’ health and productivity. Climate projections over the next decades show a continuous and accelerated warming over Europe together with longer, more intense and more frequent heatwaves on regional and local scales. In this work, we assess the increased risk in future occupational heat stress levels using the wet bulb globe temperature (WBGT), an index adopted by the International Standards Organization as regulatory index to measure the heat exposure of working people. Our results show that, in large parts of Europe, future heat exposure will indeed exceed critical levels for physically active humans far more often than in today’s climate, and labour productivity might be largely reduced in southern Europe. European industries should adapt to the projected changes to prevent major consequences for the workers’ health and to preserve economic productivity.


Introduction
Mean surface temperature as well as the frequency and duration of heat waves have increased in the last decades in Europe (IPCC 2013), in particular in the south-east of the continent (Pogačar et al. 2018; Barriopedro et al. 2011;Morabito et al. 2017). Also, future heat waves are very likely to be more frequent and longer-lasting (IPCC 2013), mainly as a direct consequence of the increase in mean temperatures (Schär, et al. 2004;Fischer and Schär 2010). Those changes relate to increasing environmental heat exposure throughout the twenty-first century (Willett and Sherwood 2010;Zhao et al. 2015;Knutson and Ploshay 2016;Coffel et al. 2018;Li et al. 2018;Matthews 2018) which in turn might have an effect on mortality, well-being and labour productivity (Dunne et al. 2013;Kjellstrom et al. 2018;Mora et al. 2017;Flouris et al. 2018;Levi et al. 2018;Moda et al. 2019). The combination of environmental heat exposure and internal heat production generated from metabolic processes results in heat stress (Xiang et al. 2014). Air temperature is an important aspect, but all relevant environmental factors should be considered in order to provide improved early warning systems and future projections of heat exposure (Li et al. 2018;Matthews 2018). Under high-airtemperature situations, the only means for the body to remain within healthy temperature limits is through loss of heat via sweat evaporation (Kenny and Flouris 2014). However, high air humidity and lack of wind limit sweat evaporation and, hence, heat dissipation that, enhanced by intense solar radiation, leads to a rise in core body temperature (Parsons 2014). This may increase the risk of heat-caused illnesses such as heat cramps, heat syncope, heat exhaustion and heat stroke (Koppe et al. 2004).
Excessive heat exposure greatly affects not only the more vulnerable population groups (such as children and elderly people) but also those performing physically demanding activities like physical exercise (Junge et al. 2016) or work (Kjellstrom et al. 2009a). Improper clothing (due to personal work protection or cultural reasons) and/or absence of cooling systems exacerbates the adverse heat conditions even further (Kjellstrom et al. 2009b). Environmental heat is a key factor for the workforce, and it is associated with a reduction of physical productivity (Wyndham 1969;Kjellstrom et al. 2009b;Sahu et al. 2013;Ioannou et al. 2017). Therefore, a better knowledge of environmental heat exposure is decisive for protecting population health and is a key management information for industries to plan labour activity and anticipate changes in productivity.
There is a large number of documented indices (see, e.g. Coccolo et al. 2016) that quantify heat exposure by combining several meteorological drivers (most of them based just on air temperature and humidity, and a few include wind speed and solar radiation). In the present work, heat exposure is represented by the wet bulb globe temperature (WBGT), since (1) it is the most widely used index to assess heat stress on working people, (2) it can be calculated from standard meteorological parameters including also the effect of solar radiation and (3) it can be interpreted via international labour standards (ISO 1989(ISO , 2017. Employing the purely meteorology-based WBGT, this work paves the way for further studies based on other heat indices that potentially involve physiological criteria and that are relevant for other impact analyses (e.g. mortality, Mora et al. 2017).
Previous works on climate projections of environmental heat accounted for shaded conditions only and were developed with a global perspective (Willett and Sherwood 2010;Zhao et al. 2015;Knutson and Ploshay 2016;Coffel et al. 2018;Li et al. 2018;Brouillet and Joussaume 2019), highlighting the dangerous conditions in densely populated areas (e.g. Southeast Asia, Africa), but the detailed impact in Europe is often overlooked. In order to assess future changes in environmental heat exposure on a pan-European level, the comprehensive and state-of-the-art regional climate model (RCM) ensemble of the EURO-CORDEX initiative Kotlarski et al. 2014) is exploited in the present work. Model simulations are statistically adjusted to the local, site-specific climate at more than a thousand locations in Europe. The effect of the projected changes of environmental heat on occupational settings is shown as an example of application.
This work is structured as follows. The data and study methods are described in 'Data and methods'. 'Results' shows climate change projections of heat exposure and impacts on labour productivity in Europe. The main conclusions are summarized in 'Discussion and conclusions'.

Wet bulb globe temperature and derived indices
Two implementations of the WBGT are used in order to account for shaded (WBGT shade ) and sunny conditions (WBGT sun ). WBGT sun (Liljegren et al. 2008) takes into account air temperature, dew point temperature, wind speed and solar radiation, whereas WBGT shade is a simplified version based on air temperature and dew point temperature only (Bernard and Pourmoghani 1999), assuming a wind speed of 1 m/s, which is the apparent wind created by limb and torso for actively working people (i.e. equivalent to a slow walk), and no heat from radiation. The reader is referred to Lemke and Kjellstrom (2012) for a comparison of WBGT calculations and the detailed formulations. Those two versions of WBGT have been implemented in the R package HeatStress, 1 under license GPL-3.
In this work, we target daily maximum heat exposure, obtained from the daily values of the input variables that enhance it, i.e. daily maximum air temperature and solar radiation and daily mean dew point temperature wind speed (see Annex A and Fig. S1 in the Supplementary Material). Several derived indices accounting for heat exposure risk are obtained from the daily WBGT values and presented for the summer season for the sake of brevity (Table 1). We acknowledge that the chosen metric (daily mean or maximum) and temporal aggregation (monthly, seasonal) have implications on the resulting indices, as discussed in the Supplementary Material (Annex B and Table S1).
Several organizations (ISO 2017; American Conference of Governmental Industrial Hygienists (ACGIH) (2016)) have defined reference and critical values of WBGT for different classes of metabolic rate (resting, low, moderate, high and very high), clothing, and considering whether a person is acclimatized or not. The application of the different WBGT reference values allows analysing heat stress for resting conditions, elderly people or vulnerable groups of the population. In this work, for the sake of brevity, we focus on exposure limits assuming moderate work intensity (300-350 W) for an unacclimatized (26°C) and acclimatized (28°C) worker (see NIOSH (2016) for a comparison of thresholds for different exposure limits as adopted by several international institutions). This way the uncertainty due to the acclimatization state of the workers is also considered. These reference values are a conservative choice since they imply protecting most of the workforce (depending on the individual conditions, one worker can be vulnerable to suffer from heat stress under lower or higher WBGT values); therefore, they should be understood as an example of application. For the sake of conciseness, we will refer to WBGT as overall environmental heat exposure and conditions with WBGT above 26°C as moderate heat risk and above 28°C as high heat risk throughout the manuscript.
There are a few approaches to estimate productivity loss at country and individual levels. For instance, Sahu et al. (2013) and Junge et al. (2016) quantify productivity as production and power output, respectively, whereas Ioannou et al. (2017) introduced a time-motion analysis to assess the effects of workplace heat on productivity and labour effort. Regardless of the method used, a graduate loss of productivity (UNDP 2016) and economic impact (Orlov et al. 2019) with increasing thermal exposure is evident. Exposure-response relationships are usually established between hourly heat exposure and productivity (e.g. WBGT above 31°C under moderate work implies 25% reduction of labour productivity in one hour for specific epidemiological studies, Kjellstrom et al. 2018). We approximated hourly values with daily mean and maximum WBGT values following the 4 + 4 + 4 method , which is a good approximation compared to more complex and computationally intense temperature models (Bilbao et al. 2002). This method assumes the daily maximum WGBT value during the 4 central hours of the day (12-16 h), the daily mean WBGT during 4 h in the early morning (8-10 h) and the early evening (18-20 h), and the remaining 4 h in between are approximated with the average of daily mean and maximum values (10-12 h and 16-18 h). Daily mean values of WBGT are derived from daily mean values of all the input parameters. Combining the hourly WBGT with the exposure-response curve from ISO (1989), we obtain the percentage of hours lost due to heat exposure, assuming that the workers are potentially active 8 h/day (from 9 to 17 h).

Observational data
Due to the lack of dense observational datasets for all meteorological parameters needed in the WBGT calculation, daily station data from different data sources are used. First, we consider the European Climate Assessment & Dataset (ECA&D, Klein Tank et al. (2002)), which is the basis for the E-OBS gridded product (Haylock et al. 2008) and contains series of daily observations at meteorological stations throughout Europe and North Africa. Daily series for maximum and mean temperature, relative humidity and wind speed were downloaded in June 2016 from www.ecad.eu.
Second, due to the low station density for non-standard parameters (e.g. humidity, wind), we combine the available ECA&D stations with the Global Surface Summary of the Day (GSOD, Smith et al. 2011) dataset which is based on data exchanged under the framework of the World Meteorological Organization (WMO) World Weather Watch Program. Daily series of daily mean and maximum temperature and daily mean dew point temperature and wind speed were downloaded from ftp://ftp.ncdc.noaa.gov/pub/data/gsod/ in July 2016.
For both products, only the stations with less than 20% of missing data in the period 1981-2010 for the three variables were considered (Fig. S2) after an additional filter for outliers. The two datasets were combined into a single one, and for the station locations where ECA&D and GSOD are available, priority was given to ECA&D. We additionally made use of the station set from the SwissMetNet 2 (SMN) dataset over Switzerland, which is provided and maintained by MeteoSwiss. Summarizing, the finally merged dataset considered in this work consists of 1370 European stations, coming from 351 ECA&D stations (Fig. S2, green), 976 GSOD stations (blue) and 43 SMN stations (red).
In order to cover all European stations for which the rest of the variables was available, we consider satellitederived daily downward surface solar radiation data (SARAH from the CMSAF surface radiation dataset, Pfeifroth et al. 2017) that are based on measurements in the visible range from the Meteosat First and Second Generation Satellites (Posselt et al. 2012;Posselt et al. 2014). They operated from 1983 to 2015, with best quality from 1994 on (Posselt et al. 2014). The satellite data are provided on a grid of 0.05°horizontal resolution (approximately 5 km). Daily mean and maximum radiation were calculated from hourly mean values retrieved for the analysis period  in February 2017. In order to combine these gridded data with the station data, the closest grid box from the satellite data to each of the 1370 stations was considered.

Model data
Regional climate model (RCM) simulations from the Coordinated Regional Downscaling Experiment (CORDEX, Giorgi et al. (2009) In this work, we use 84 simulations (GCM-RCM chains) in total performed by 11 RCMs driven by 10 different global climate models (GCMs), at two horizontal resolutions (0.11°a nd 0.44°, approx. 12 and 50 km) and assuming three Representative Concentration Pathways (RCP2.6, RCP4.5 and RCP8.5). Historical runs for the period 1981-2005 and future projections until 2099 were considered here. Since we use the reference period 1981-2010, the years 2006-2010 from the projection runs were merged with the historical runs. Most of the simulations were extracted from the ESGF in May 2017, and additional EURO-CORDEX-compliant simulations from ETH Zurich were also included in the analysis (Sørland et al. 2018). The reader is referred to Table S2 for a summary of the simulations used.

Bias correction
Despite their continuous improvement, RCMs are prone to systematic biases and their resolution is still too coarse for use in many sectoral climate change impact applications (Christensen et al. 2008). A common approach to account for systematic biases and to further downscale RCM results is to use distribution-based statistical transfer relations that adjust/correct systematic model biases and that might also include a downscaling component (see, e.g. Déqué (2007); Piani et al. (2010)). Among the available bias correction methods, we use empirical quantile mapping (QM) that consists of matching the simulated and observed distributions by establishing a quantile-dependent correction function, between the observed and simulated quantiles (Panofsky and Brier 1968). QM is nowadays a well-established method and was used to produce localized climate information from national climate change scenarios, e.g. in Switzerland (CH2018 2018) and Austria (Formayer et al. 2015). The implicit assumption of QM is that a climate model can sufficiently project ranked categories of the variable of interest, i.e. quantiles, but not its actual values (Déqué 2007). Here we use the implementation from Rajczak et al. (2016), where the 99 empirical percentiles of daily data are corrected and linear interpolation is used for the values between two percentiles. Constant extrapolation is applied for values outside the calibration range, i.e. the correction function of the last (first) percentile is applied to all the values above (below) it. The calibration is carried out independently for each day of the year with a moving window of 91 days, considering the full period 1981-2010, for the closest grid box from the RCMs to each of the 1370 stations. The correction functions are then applied to the historical period  and the future period  or transiently for 1981-2099 (Fig. 1). Each of the input parameters of the WBGT (daily maximum temperature and solar radiation and daily mean dew point and wind speed) is corrected independently prior to the index calculation. Wilcke et al. (2013) showed that univariate bias correction is able to retain the temporal structure and the inter-variable relationships of the uncorrected data. Different to multivariate bias correction methods (Vrac and Friederichs 2015;Cannon 2016), univariate techniques are less complex and easier to interpret and the parameter estimates of the transfer functions are consequently more robust. The correction of the individual variables prior to the calculation of a multivariate index or impact model is still the common and preferred approach to downscale climate indices (Teutschbein and Seibert 2012;Casanueva et al. 2018), even though it assumes statistical independence of the variables.
The EURO-CORDEX RCMs do not provide daily maximum solar radiation, and its approximation from the available variables in the model is not straightforward. To circumvent this problem, we bias-correct daily mean solar radiation from the RCMs against daily maximum observed solar radiation. By doing this, the distribution of the mean solar radiation of the models is mapped onto the distribution of the maximum observed counterpart. Thus, the correction function is not a pure bias correction, but also translates daily mean solar radiation into daily maximum solar radiation. The suitability of this approach was tested with the SMN data (Fig. S3). Although results show an underestimation of the daily variability in the quantile mapped maximum solar radiation with respect to the observed counterpart, the effect on the WBGT is small, probably because solar radiation does not play a major role in the heat exposure index (compared to temperature and  Fig. S4), and WBGT sun can be fairly well reproduced.

Evaluation of the bias correction
One of the potential limitations of QM is that the inter-variable dependencies may be modified (Ehret et al. 2012;Vrac and Friederichs 2015) since it is a univariate bias correction method and does not explicitly correct for RCM biases in the physical, temporal or spatial relationships among variables. Therefore, it is crucial to verify if the WBGT calculated with the bias-corrected input variables agrees with the observed counterpart. Figure S5 shows the Q-Q plots of the observed vs. simulated WBGT (only summer days) for shaded and sunny conditions, for the raw (black) and bias-corrected data (blue) for five representative stations in Europe. They are representative in the sense that they cover different climates and a north-to-south gradient along the continent. Moreover, they are important European cities, with many inhabitants, which might be suffering from an increase of environmental heat exposure in the future. It is shown that after the independent correction of the input variables, the heat indices agree very well with the observations (i.e. the blue dots lie approximately on the diagonal). Thus, the inter-variable relationships are obviously retained after the correction of the individual parameters (Wilcke et al. 2013). Note that results may be less accurate in a cross-validation framework, although our approach can be considered independent since the evaluated aspect (here: multivariate consistency) is not directly tackled by QM.

Future increase in summer maximum heat exposure
Previous works have shown that hot temperatures (on single days) increase faster with global warming than average temperatures (Schär et al. 2004;Fischer and Schär 2009;Cattiaux et al. 2015); thus, we first focus on the most extreme case of summer maximum WBGT. Our results show that, overall, environmental heat exposure is projected to increase in the course of the twenty-first century in Europe. For an illustrative set of European cities (Fig. 1), summer maximum WBGT (WBGTx3d, i.e. maximum of moving 3-day average WBGT; see Table 1) observed during the past 30 years ranged from 22°C in Oslo to 30°C in Seville for shaded conditions. Due to the effect of solar radiation, WBGT in the sun is systematically higher than WBGT in the shade, with differences of around 2-3°C or even more for the Mediterranean stations. Until 2050, heat exposure is projected to increase similarly for all RCPs considered, while after 2050, the changes according to the three RCPs show distinct behaviour, especially for the more southerly locations. At the five illustrative stations and especially for RCP8.5, heat exposure might reach far beyond the observed range with summer maximum WBGT up to 26°C in Oslo and 34°C in Seville in the shade. As for the observational period, projected WGBT values in the sun are systematically higher by 2-3°C. Model uncertainty together with interannual variability (shaded range) is larger at the northern stations (4-5°C) and smaller for the southern locations (2-3°C). The uncertainty range increases in the course of the century for all of the stations and all RCPs. At Northern and Central European stations, the selected thresholds for moderate and high heat risk (orange-shaded areas in Fig. 1) are occasionally reached in the present but those conditions will occur much more often in the future.
Single-day summer maximum WBGT (WBGTx1d) increases by 1-4°C by the end of the twenty-first century with respect to the historical reference period depending on the specific station and on the emission scenario (Fig.  S6). Changes of WBGTx1d are more uniform across the continent than air temperature changes which are higher in Southern Europe than in the north Kröner et al. 2017). This might be due to the projected decrease in relative humidity in the Mediterranean region (Ruosteenoja and Räisänen 2013), which counterbalances the strong temperature increase and hence results in a smaller WBGT change (note the effect of increasing temperature and decreasing dew point temperature on WBGT in Fig. S4), in agreement with Brouillet and Joussaume (2019). In general, future environmental heat exposure is enhanced in those regions where it is already a problem in today's climate, resulting in a clear northsouth heat exposure gradient (Fig. 2). For instance, in Southern Europe, WBGTx1d in the shade reaches 27-29°C in present-day climate and might amount to 28-29°C, 28-30°C and 30-32°C for the three RCPs, respectively. Accordingly, WBGTx1d in the sun is projected to rise from 30-32°C to 31-33°C, 32-34°C and 34-36°C, respectively. In Central Europe, where nowadays WBGTx1d ranges between 24 and 26°C (in the shade) and 27 and 29°C (in the sun), it may reach up to 30°C and 32°C under the strongest emission scenario (RCP8.5) in many locations.
The climate change signal for summer mean WBGT (WBGTmean) is very similar to that for the summer maximum (WBGTx1d, see Fig. S7), unlike the amplification found in temperature extremes compared to mean temperature changes (Fischer and Schär 2010). A possible reason for this could be the role of humidity (dew point temperature) which contributes to WBGT and, thereby, modulates the pure temperature signal. Dew point temperature shows very similar changes of summer mean and maximum values (Fig. S7, blue points).
High heat exposure to become more frequent in the future The described increase of summer mean and maximum WGBT is associated with more frequent situations with high heat risk (Fig. 3). In the shade, the number of summer days with high heat risk (WBGTg28) currently reaches up to 20 days on average in a few locations in Spain, Italy, Greece and Cyprus (Fig. 3, left), whereas in the sun, the area with such conditions spans over Southern Europe and reaches 30-50 days in some Mediterranean locations. Frequencies of high heat risk for both shaded and sunny conditions are projected to increase by 5-15, 15-30 and 30-50 days per year in the Mediterranean region by the end of the century with respect to today's climate for the three RCPs, respectively (Fig. S8). Hence, on average, there will be 5-15, 10-40 and 20-70 days with high heat risk per summer season for shaded and 20-50, 30-60 and 50-80 for sunny conditions for the three RCPs in the Mediterranean area (Fig. 3, middle and right).
In some Southern European locations, workers might suffer high heat risk at the end of this century for almost the entire summer in sunny conditions (values close to 100%) and more than half of the summer days under shade conditions ( Table 2). The comparison between the frequencies of WBGT above 26°C and 28°C (see 'Wet bulb globe temperature and derived indices') partly samples the uncertainty due the acclimatization state of the worker for the specific case under moderate work intensity (NIOSH 2016;ISO 2017). Although the frequency of heat risk is lower for acclimatized workers (Table 2, right columns: WBGTg28), future changes with respect to present-day climate might be up to 2 times larger than for unacclimatized workers (WBGTg26) in Mediterranean locations for WBGT in the sun. In these locations, the frequency of high heat risk is substantial even for the low-emission scenario. While high risk in shaded conditions is confined to the Mediterranean area, it extends northwards for WBGT in the sun, affecting the entire continent except Scandinavia and the British Isles by 2100 (cf. Fig. S8). High Observed and projected number of days with high heat risk. As Fig. 2, but for the frequency of days with high heat risk (WBGTg28). The two maps refer to shaded conditions heat risk (WBGTg28) in the sun is projected to occur during up to 10% of summer days in many locations in Central Europe for the low-and medium-emission scenarios and even on 20-30% of summer days for RCP8.5 ( Fig. 3 and Table 2). Thus, many workers in Europe need to adapt to a future with more frequent heat situations, even when performing moderate work.

Model uncertainty in the climate change signal
Uncertainties in climate projections often represent a challenge for climate change communication and its interpretation by stakeholders. Climate indices, such as heat indices, offer the possibility to describe the climate in a more approachable way. The joint assessment of uncertainties using multivariate indices can potentially increase or compensate model uncertainties depending on the case (Fischer and Knutti 2013). The uncertainty in the climate change signal of WBGT would be modulated by the model uncertainty of all input variables (air and dew point temperature, wind speed and radiation) and the relationships among them (Fischer and Knutti 2013). Here, the ratio between the multi-model mean change (signal) and standard deviation of the multi-model changes (uncertainty) is used as an indicator to compare the model agreement relative to the projected change for summer m aximum WBGT (WGBTx1d) in the shade and the input variables (Fig. 4). This ratio is large for most of the stations in the Mediterranean region and Central and Southeast Europe, especially for RCP8.5 for which larger signals are projected (Fig. 4i). That means that the uncertainty is small compared to the signal. Signal-to-noise ratios are larger for WGBT and dew point temperature than for maximum temperature over many regions, implying that the incorporation of dew point temperature increases the robustness (Fig. 4). This result is consistent with previous findings (Fischer and Knutti 2013) and highlights the potential source of predictability provided by relative humidity, which has been found in climate change assessments (Scoccimarro et al. 2017) and seasonal forecasts . The north-south pattern in the signal-tonoise ratios of WBGT in the shade (Fig. 4g-i) may explain the differences in the model uncertainty by the end of the twenty-first century in Fig. 1, since in Northern Europe there is larger uncertainty for air temperature projections than in the south.
More robust results for WBGT than air temperature also hold for RCP2.6 although with smaller signal-to-noise ratios, due to the smaller signals, and large uncertainties in Northern Europe where two simulations project a much larger temperature increase than the RCP2.6 ensemble mean.
Other sources of uncertainty beyond model uncertainty are discussed in the Supplementary Material (Annex C). Table 2 Observed and projected frequency (% of summer days) with moderate (WBGTg26) and high (WBGTg28) heat risk in the sun and in the shade, for five exemplary stations and the three RCPs in the period 2070-2099. A value of 100% would mean that the entire summer is at risk. Projections represent the multi-model median and the 80% of the multi-model range (10th to 90th percentile, in brackets)

Productivity losses due to heat exposure applying ISO recommendations
The workers' natural response to protect themselves against the risk of heat-related illnesses is to slow down work intensity and/or limit working hours, thus minimizing body heat production and reducing heat exposure, respectively (Ioannou et al. 2017). As a consequence of these strategies, labour productivity and economic output are reduced (Kjellstrom et al. 2009a;Kjellstrom, et al., 2009b). Staying below certain exposure limits (e.g. ambient temperature below 35°C (Kjellstrom et al. 2009b) and core body temperature below 38°C (Kjellstrom et al. 2009a)) reduces the risk of heatrelated illness, but does not preclude the possibility of other adverse effects such as a loss of productivity. For that purpose, specific epidemiological studies (Wyndham 1969;Sahu et al. 2013) have been used to estimate relationships between exposure (heat level) and response (productivity) . For moderate work intensity (300-350 W), when WBGT is above 31°C, the hourly work capacity is reduced b y 2 5 % f o r a n a v e r a g e w o r k e r a c c o r d i n g t o these epidemiological studies .
International standards (ISO 1989) would be more restrictive than that (70% of losses for WBGT above 31°C, Kjellstrom et al. 2018), because they are based on heat levels that protect most of the workforce. We quantify the productivity losses (reduction in available work time) that would occur if the international standard recommendations are enforced (ISO 1989, see 'Wet bulb globe temperature and derived indices'). If working people do not reduce their work as much, the losses would be less but at the expense of other effects on workers' well-being (Flouris et al. 2018). These ISO recommendations apply to acclimatized workers; thus, they might be more appropriate for Southern European cities. In present-day climate, there are only few locations in the Mediterranean area where up to 20% of summer working hours are lost under sunny conditions (Fig. 5,  left). By the end of the twenty-first century, under RCP8.5, Southern Europe will experience a widespread loss of working hours by at least 15%, reaching more than 50% in some locations in Spain, Italy, Greece and Cyprus (Fig. 5, middle). The different estimates of productivity losses vary nonlinearly with the three emission scenarios, and substantial differences are projected between RCP4.5 and RCP8.5 for specific hot spots such as Rome and Seville (Fig. 5, right).

Discussion and conclusions
Changes in environmental heat exposure due to climate change are highly relevant to stakeholders and policy makers in order to respond to possible impacts on the main European industries. At the end of the twenty-first century, many European workers will very likely be affected by heat stress, not only due to the increase of heat exposure but also because such situations will become more frequent in large areas of the continent. This will affect not only the south of the continent but-especially for workers active in the sun-also regions in Central and Northern Europe, where heat exposure has a smaller effect in present-day climate. Even if stronger global mitigation actions are implemented (RCP2.6), high heat risk is found for large parts of Southern Europe during the twenty-first century, underlining the need for reconsidering international working regulations for European industries to reduce heat exposure. A direct consequence of environmental heat is the loss in labour productivity especially for sunny conditions, which can result in a reduction of 15-60% of the working hours in the Mediterranean area under the strongest emission scenario by the end of the twenty-first century. Model uncertainty may be large, but the current study shows a strong agreement among climate models in terms of an increased heat risk. Furthermore, climate change projections of summer maximum heat exposure are more robust than those for maximum air temperature, mainly due to the higher model agreement in dew point temperature. We focus on the effects of summer heat exposure on labour productivity, but a broader picture would consider non-linear effects (Burke, et al., 2015), e.g. future warming might increase outside labour productivity during the winter months in northern Europe, thus compensating the loss in summer due to high temperatures.
In light of the results presented here, mitigation and adaptation measures should be taken to prevent productivity losses in the course of the twenty-first century. Some protective strategies to alleviate heat exposure might be heat wave monitoring and warning, a reduction of sources of heat in workplaces, a reduction of physical work intensity, personal protection through movable personal microclimate cooling and sophisticated technical developments in clothing based on cooling Fig. 5 Productivity loss due to environmental heat exposure. Observed (left) and projected (middle) percentage of summer working hours lost due to heat exposure under sunny conditions. Projections show the multimodel ensemble median (over 39 model chains) for the strongest emission scenario RCP8.5 and the period 2070-2099. Boxplots (right) summarize the projected percentage of working hours lost by the end of the century for the two WBGT implementations (shade and sun) and three selected European stations. Each boxplot represents the multi-model uncertainty range for each RCP and location. To obtain the percentage of summer working hours lost due to heat exposure, we approximate hourly values of WBGT by combining daily mean and maximum values of WBGT  and apply the exposure-response relationship from ISO (1989) (see 'Wet bulb globe temperature and derived indices') with ventilation and phase change materials that absorb or release latent heat when they change phases (Gao et al. 2018). Since reducing heat exposure by the workers (i.e. limit working hours, increase pauses) derives in productivity losses, there might be, however, a tradeoff between reducing heat exposure and preventing productivity losses that needs to be assessed case by case. In order to cope with this situation, specific solutions need to be applied for each case, including feasible, sustainable and effective solutions (Morris et al. 2020).
The analysed RCM simulations do not include a dedicated module for urban climates, and the implemented land use is static. Moreover, meteorological stations as those used here for the bias correction are usually located outside the cities, such as at airports or in rural environments. Thus, the urbanization effect is underestimated and not fully accounted for in our analysis. Further research is needed to assess the urban heat island effect and to quantify future heat exposure in urban areas on a European scale. The urban heat island effect is, for instance, associated with warmer nights in urban compared to rural areas (Oleson et al. 2011;Parlow et al. 2014). In a warmer climate, a higher frequency of high-heat-stress nights and tropical nights in urban compared to rural sites is projected (Fischer et al. 2012; Burgstall 2019 and references therein). This is highly relevant since people in cities might not recover from the daytime heat and might subsequently not be able to handle any extreme heat the following day (Perkins, 2015). So, the results presented in this study should be regarded as a lower bound of future heat risk in urban areas.
Many opportunities arise based on this work. The climate change projections of environmental heat exposure could be combined with demographic data and economic models in order to quantify economic losses due to heat (as foreseen in the HEAT-SHIELD project; see www.heat-shield.eu). At shorter time scales, a heat-warning system has been developed in order to allow stakeholders timely and precise prevention strategies and better planning of the work activities up to four weeks in advance (see http://heatshield.zonalab.it). Accordingly, the development and dissemination of heathealth planning and warning systems is now among the priorities of the World Meteorological Organization (WMO) and the World Health Organization (WHO). All these examples are certainly a good test bed for developing effective climate services in the context of extreme heat. CH2018 group (www.ch2018.ch) for technical support and scientific discussions. Financial support for this work is provided by the HEAT-SHIELD Project (European Commission HORIZON 2020, research and innovation programme under the grant agreement 668786). The authors wish to thank the Swiss National Supercomputing Centre (CSCS) for providing the technical infrastructure.
Finally, the authors are also grateful to the editor and two anonymous reviewers who helped to improve the original manuscript.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.