Regional disparities and seasonal differences in climate risk to rice labour

In this study we analyse data from the Coupled Model Intercomparison Project, Phase 6 (CMIP6) [14], comprised of global climate models that have been run in a shared experimental configuration. The model outputs used here are daily mean and maximum temperature, specific humidity, and surface-level air-pressure. All 14 CMIP6 models for which appropriate data were present in the Centre for Environmental Data Analysis (CEDA) archive ⁸ were included. CMIP historical runs and ScenarioMIP future pathways were processed. ScenarioMIP runs simulate the future climate given assumptions about future emissions (especially of greenhouse gases) and development pathways called shared socioeconomic pathways (SSPs) [15]. So that models with multiple ensemble members were not given additional statistical weight, a single ensemble member was used for each model, selected in numerical order from what was available in the archive. A table of model runs used, with data citations, is included in the supplementary material available online at stacks.iop.org/ERL/16/124004/mmedia (supplementary tables 2 and 3). Emissions trajectories in ScenarioMIP are derived from socioeconomic assumptions; we only made use of the atmospheric variables output from the climate models and did not make any direct use of the socioeconomic data. The Climatic Research Unit gridded Time Series (CRUTS) 4.03 [16, 17], and the European Centre for Medium-Range Weather Forecasts Reanalysis 5 (ERA5) [18], were used as historical observational datasets.

Daily frequency CMIP6 data were used. Sub-daily data for the variables of interest were available for only a small number of models. Daytime temperature variation was estimated by assuming the temperature is close to the daily mean temperature, daily max temperature, and the mid-point of the two, for 4 h each, following Kjellstrom et al [6]. Comparing the result of this calculation using 3-hourly and daily CMIP6 data, we found that this is a reasonable approximation. Daily ERA5 data were used to maintain consistency with the CMIP6 data. CRUTS 4.03 is monthly and is only included as a crosscheck.

Wet-bulb globe temperature (WBGT) is a heat-stress metric defined by ISO 7243 and widely used for assessing hazards due to hot conditions [19–21]. WBGT is intended to combine the factors that affect the human experience of heat: namely air temperature, radiant temperature, humidity, and air velocity [20]. WBGT is designed to be measured directly using specialised equipment; in practice, statistical and empirical formulae for estimating it from standard meteorological variables must be used in the climate context. We assumed that work occurs in the shade, so that air temperature approximates black-globe temperature. Furthermore, cloud cover is one of the most uncertain aspects of global climate models [22]. In the supplementary material (supplementary section 2.1), we explore the effect of excluding radiation. Although this leads to underestimation of heat stress in some conditions, we find that long-term trends in WBGT are determined mainly by changes in air temperature. Furthermore, this assumption is used in other related studies [6, 11, 12].

WBGT was calculated as ${\text{WBGT }} = { }0.7{\text{*WBT }} + { }0.3{\text{*BGT}}$, where WBT is the wet-bulb temperature. WBT was calculated from air temperature, specific humidity, and air pressure using open-source software 'psychrolib' [23] which implements formulae from the American Society of Heating, Refrigerating and Air-Conditioning Engineers handbook. Field measured WBT decreases with wind speed at low speed (<2 m s⁻¹, a light breeze), but higher wind speeds have a lesser effect [24]. The WBT calculation used is consistent with a light breeze, and variation in wind speed is neglected.

Studies of occupational heat stress under climate change often either assume a threshold in WBGT above which a worker is at risk [25], or assume that there is a simple relationship between WBGT and worker productivity. Typically, these relationships are assumed to be representative across sectors, and are based on either a regulatory or advisory standard [11], limited field study or survey data, or an ad hoc fusion of the two [6, 12, 26]. Occupational heat-stress regulations attempt to minimise harm and are often exceeded in practice. Therefore, use of regulatory standards to estimate future labour-productivity, as in Dunne et al [11], may not be accurate. Field data measuring the effect of heat on worker productivity are sparse and cover only a few activities. Studies of small numbers of workers are often used to estimate productivity effects on the entire human population.

Sahu et al [4] observed a 5% per °C WBGT decrease in the labour capacity of labourers harvesting rice between 26 °C and 32 °C WBGT. Rate of collection was measured in 124 workers in groups of 10–18, and WBGT was measured in-situ, at an individual location in India. For our study, we assume that this is representative of manual rice harvest labour, and adopt it as our labour impact metric. By labour impact metric, we mean the amount that labour productivity is reduced by heat stress. The impact is assumed to be linear in WBGT, although this assumption must break down at high WBGT (>35 °C). The systematic uncertainty due to these assumptions cannot be assessed without larger scale field observations. The labour impact metric is calculated as $5.14{\text{*WBGT}} - { }218$, in units of %, clipped at 0 and 100, meaning that 0% loss occurs at 23 °C and 100% loss occurs at 42.5 °C. This function comes from interpretation of the Sahu et al [4] study by Gosling et al [26]. By using a labour impact derived from a study specific to manual rice harvesting, we represent this particular activity more accurately than the assumptions used for the economy in general in other studies (e.g. [6, 11, 12, 25–27]).

To facilitate comparison to other studies, results were also calculated using different assumptions about the relationship between WBGT and labour productivity: the assumptions used for agricultural labour in Dunne et al [11], and in Orlov et al [27]. Our calculations are consistent with these other studies above. In the supplementary material (supplementary section 2.2), equations for each of these and the effect of varying the assumptions are shown.

The RiceAtlas dataset [8, 9] provides detailed data on rice production, broken down into location entities with an average area of 5000 km². The distribution of rice cropland in Asia is shown in figure 1. Information such as yield, harvested area, planting and harvesting dates are included, and many entities have multiple yearly harvests. The data are representative of the years 2010–2012. We used this to identify harvest dates. Figure 2 shows the distribution of rice harvests through the year across the whole of Asia according to RiceAtlas. Figure 3 shows where harvests occur in each calendar season. Harvests peak in September-November, but there is a large amount of spatial variation. For a more complete illustration of the timing of harvest seasons in the RiceAtlas dataset, we refer to the paper describing that dataset.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Location entities in RiceAtlas vary greatly in size, and many are much smaller than the grid spacing of the climate models considered. Where a location geometry enclosed multiple climate-model grid-cells, the mean was taken of the enclosed cells; otherwise, the land grid-point closest to the centroid of the location was used. We used data only for countries in Asia; see table 1.

For each location and harvest season, climate data (ERA5 and CRUTS 4.03, and CMIP6 outputs) were selected using the harvest dates and spatial geometry provided by RiceAtlas, and WBGT calculated. Only climate data during the peak month of each harvest were selected. Then, the labour impact metric was calculated from WBGT. Independently for each climate model, and for each location and season, a linear trend was fitted between the 20 year means of the change in the global mean near-surface air temperature (GSAT) and the labour productivity metric. Twenty-year periods were defined between 1850 and 2009 inclusive for the historical runs, 2020–2099 inclusive for the SSP runs. As the labour productivity function was clipped at zero, the fitted function was too. We assumed that 1 °C of GSAT warming relative to 1850–1900 was representative of the present (2020) [28], and changes in both variables are relative to the present. Data from historical model runs and different SSPs were included together in a single fit, so that a range of GSAT values (up to 4 °C relative to 1850–1900) were included; but each model, and harvest season and location was fit independently. We defined the gradient of this fit as the 'hazard gradient', the purpose of which is to summarise the relationship between changes in GSAT and local changes in labour productivity via WBGT, to explore the spatial patterns of change. The hazard gradient means the amount by which the labour impact increases on average for 1 °C of global warming. This effectively collapses variation in time, the climate sensitivity of each model, and emissions scenario into a single dimension. The two-sided inverse Student's t-distribution was used to determine whether the hazard gradient in each harvest season and location was significant.

The multi-model mean of the hazard gradient was calculated for each harvest season and location. The labour impact was also averaged over all locations, weighted by total rice production, to characterise the regional situation.

Statistical tests were applied to the results. In order to test if there was a significant change in the historical period, linear regression was used on the time-series of the labour impact metric calculated from ERA5 and CRU-TS, and the two-sided inverse Student's t-distribution used to determine if the gradient was significantly different from zero. To test if there is a significant difference between different warming scenarios, the distributions of results across climate models at 1.5 °C, 2.0 °C, and 3.0 °C of global warming were compared using the Kolmogorov–Smirnov (KS) test and Student's t test.

Figure 4 shows the Asia-wide average of the heat impact on rice harvest labour, calculated from the historical observational and reanalysis datasets, weighted by total rice production in each harvest season and location. There is a statistically significant increase (p ≪ 0.01) over the full observational period (1900–2018) of CRU-TS, and in both datasets for the common observational period (1980–2018): the heat hazard associated with the rice harvest has already increased.

Zoom In Zoom Out Reset image size

Figure 5 shows heat impact on Asia-wide rice harvest labour plotted against GSAT warming across various climate models. Changes are relative to the present, assumed to mean 1 °C of warming over 1850–1900. The rice harvest labour impact is weighted by total rice production. In figure 5(a), each point is a 20 year mean in a single model simulation. Despite the different biases and climate sensitivities of the models, they each show a linear relationship between GSAT and the rice harvest labour impact. All models have a hazard gradient 2.0–2.7% °C⁻¹, p ≪ 0.01, R² > 0.98. The multi-model mean (standard deviation) of the hazard gradient is 2.3 (0.2)% °C⁻¹. A single trend line is shown in figure 5(a) for illustration only.

Zoom In Zoom Out Reset image size

In figure 5(b), the Asia-wide labour impact is linearly interpolated to 1.5, 2.0 and 3.0 °C of warming (relative to 1850–1900) individually for each model, and results are shown as a multi-model boxplot. There is a statistically significant difference between the 1.0 °C and 1.5 °C of warming with multi-model mean (standard deviation) 1.0 (0.3)%; between 1.5 °C and 2.0 °C of 1.3 (0.2)%; as well as between 2.0 °C and 3.0 °C of 2.6 (0.4)%. KS test and t test p-values on the multi-model distributions are small (≪0.001) for each case.

As the Asia-wide analysis is weighted according to total rice production in each location and season, the Asia-wide results should be interpreted carefully: the average is more sensitive to locations with more intensive cropping. However, the breakdown by location shows that the trends are general, and not skewed by a single area with high cropping intensity. The correlation coefficient of the hazard-gradient between models and within each harvest season and location is 0.95, suggesting that there is generally high agreement between the models.

Figure 6 shows the distribution of the multi-model mean hazard-gradient across harvest seasons and locations. In the blue histogram, the entire year is given equal weight in the heat stress calculation. Conversely, in the orange histogram, only the rice harvest season is considered in the heat stress calculation. Both histograms are weighted by rice production. Considering the time of year at which the rice harvest occurs, as opposed to the whole of the year, leads to some rice harvests having higher hazard-gradient, but other harvests not having a significant hazard-gradient. If the harvest season and locations with non-significant hazard-gradients are excluded, then the mean hazard-gradient is increased. The interquartile range of the distribution in orange is more than double that in blue. The two assumptions lead to completely different distributions of hazard gradient, showing the importance of including the seasonality of labour in projections of the effect of heat-stress on agricultural labour productivity.

Zoom In Zoom Out Reset image size

The 50% of rice production with the highest hazard gradients (hazard gradient >3.05% °C⁻¹) is identified as being exposed to a high hazard-gradient; the proportion of exposed production in each location is mapped in figure 7 and listed by country in table 1. We do this to identify regional differences in the estimated labour impact; we are not suggesting this is a limit to adaptation. Areas identified by this method include most of Southeast Asia, and coastal South Asia. Harvests in northern India (e.g. in Punjab) are not strongly affected.

Zoom In Zoom Out Reset image size

Figure 8 shows the hazard gradient for each harvest season and location, plotted against latitude in 6a and against the peak month of the harvest season in 6c. Figure 8(b) shows the hazard gradient against latitude assuming the full year is equally weighted. Two interacting effects explain most of the spatial variation in the results. Firstly, locations close to the equator have higher temperatures. Secondly, for locations that are further from the equator there is greater seasonal variation in temperature, so the time of year at which the harvest occurs is important. Figure 8(b) shows that if the time of year of the harvest is not considered, then the hazard gradient is mostly determined by latitude. Close to the equator, for example in Indonesia, hazard gradient is not very seasonally dependent as there is relatively little seasonal variation in temperature and harvests are spread through the year. By comparison, in China and India the time of year of the harvest season has a much stronger effect on the hazard gradient. Harvests in China occur June-November (see figure 3), with harvests occurring June-August more exposed than those occurring September-November. Figure 8(c) shows that, regionally, all large harvests that are not exposed to high hazard-gradients peak between September and December inclusive, whereas harvests July-August are more exposed.

Zoom In Zoom Out Reset image size

Global mean surface temperature (GMST) has increased by 0.9 °C–1.2 °C relative to 1850–1900 as of 2017 [28], and if anthropogenic warming continues to follow recent trends, is projected to reach 1.5 °C between 2030 and 2052; the estimated rate of warming is 0.2 °C per decade [29]. In accordance with the Paris Agreement under the United Nations Framework Convention on Climate Change [30], most countries have committed to limit GMST warming by 2100 to 2 °C, and to make efforts to limit warming to 1.5 °C. However, at the time of writing, pledges of emissions reductions are not sufficient to meet this goal [31].

We identify in figure 4 that there is a statistically significant increase in the estimated heat impact in the historical observational dataset. Overall, the long-term averages in the estimated labour impact are highly correlated with GSAT change (figure 5). This is not surprising as the labour impact metric is linear in WBGT above a threshold, and GSAT change (i.e. global warming) dominates large-scale, long-term trends in WBGT. This relationship holds at the Asia-wide level, but also individually for most location and seasons.

Labour inputs to rice production are seasonal, with more labour occurring around planting and harvesting [32–34]. This labour cannot be delayed or displaced to other times of year, and the timeliness of these activities affects yield and quality. Averaging across the full year assumes that labour for any given activity can be exchanged with labour at a different time of year, which is clearly not the case for planting and harvesting. In this study, we have shown the importance of considering the seasonality of agricultural labour when estimating the effect of heat-stress on labour productivity. We have focussed only on harvesting, and only on rice. However, most crops are seasonal, and therefore have a seasonal distribution of labour; this should be taken into account in estimates of the effect of heat stress on labour. Other studies into the effect of climate change on labour productivity, even those focussed on agricultural labour, do not take this into account [6, 12, 13, 27].

In some cases, for example wheat in Punjab, rice is multi-cropped with another crop that is not covered in this study, but for which harvest labour may be more exposed. Studies have suggested that changes to planting and harvesting dates [35], and crop choice [36], could be used to adapt to climate change and mitigate negative effects on yields. Changes to planting and harvesting dates, to maximise yield as temperature and precipitation patterns change, is a complex topic, and outside the scope of this study. Given that the world has already experienced 1 °C of warming, such adaptations may already be occurring. Changes to planting and harvesting dates could lead to workers being exposed to heat stress in locations where they currently are not, and this should be considered when studying adaptive measures.

Workers may be able to reduce heat exposure by starting work early in the morning when it is cooler, but this could mean working in the dark. Reducing work hours could also reduce heat exposure but could exacerbate labour shortage. We note that the workers observed by Sahu et al began their work near dawn and were still exposed to hazardous heat in their first hour of work [4].

The accuracy, coverage, and granularity of the rice crop data limits our study. However, by focussing on long-term trends with high levels of multi-model agreement, we avoid the issue of climate model bias. The assumed linear effect of WBGT on labour productivity comes from a single field study and cannot be accurate as WBGT approaches known human physiological limits, but the general conclusions are unchanged by substituting other labour-impact functions from the literature (see supplementary section 2.2). Due to the resolution of climate models, results will be less reliable for locations with steep gradients in altitude, for example Nepal.

Farming practices and machinery use mean that the required amount of labour for a given harvested area or weight will vary between places [37, 38]. Differences in mechanisation between and within LIC and LMIC are not a simple function of gross domestic product, and there is substantial variation within as well as between countries. Vulnerability may be higher in parts of sub-Saharan Africa, or for other crops, where there has been less uptake of agricultural machinery [39, 40]. This is discussed further in supplementary material section 1.

El Nino is an important driver of annual temperature variation in the region, and negatively affects crops through rainfall patterns (see e.g. Selvaraju [41]). Higher global temperatures could change El Nino patterns, but this is beyond the scope of this study. We have focussed on long term trends on the scale of 20 years rather than extreme events, but WBT extremes are also closely correlated with large-scale mean temperatures within the tropics [42].