Regional disparities in climate risk to rice labour and food security

a) British Antarctic Survey, High Cross, Cambridge, CB3 0ET, United Kingdom b) The Alan Turing Institute, 96 Euston Rd, Somers Town, London NW1 2DB, United Kingdom c) University of Bristol, University Road, Clifton, Bristol BS8 1SS, United Kingdom d) Met Office Hadley Centre, FitzRoy Road, Exeter, EX1 3PB e) University of Exeter Global Systems Institute, Laver Building, North Park Road, Exeter, EX4 4QE, United Kingdom f) Department of Computer Science and Technology, University of Cambridge, JJ Thomson Ave, Cambridge CB3 0FD, United Kingdom


Introduction
Asia produces 90% of rice, 631 million tonnes (Mt) 1,2 . A large proportion of rice cropland is tropical or subtropical (57% between 23.45 deg N and S, 95% between 35 deg N and S), see Figure 1). Agriculture employs more than 640 million rural people in the Asia and Pacific region 3 , in which rice production comprises 30% of total crop gross production value 4 .
Agricultural workers are especially vulnerable to hot and humid weather, which impacts health and productivity, and which will increase due to global warming [5][6][7] . Field studies have demonstrated the presence of heat strain and related health issues in agricultural workers in general 8,9 and rice harvest workers specifically 10 . These workers are subject to high seasonal temperatures, which can be harmful despite not necessarily being extreme statistically.
Global mean surface temperature has increased by 0.9-1.2 C relative to 1850-1900 as of 2017 11 , and if anthropogenic warming continues to follow recent trends, is projected to reach 1.5 C between approximately 2030 and 2052; the estimated rate of warming is 0.2 C per decade 12 . In accordance with the Paris Agreement under the United Nations Framework Convention on Climate Change 13 , most countries have committed to limit global mean surface air temperature (GSAT) warming to 2 C, and to make efforts to limit warming to 1.5 C. However, given pledges of emissions reductions as of 2019, warming is on course to exceed 3.2 C by 2100 14 .
Many studies focussed on the risk of occupational heat stress use wet-bulb globe temperature (WBGT), which is a heat-stress index defined by ISO 7243. WBGT is intended to combine all the factors that affect the human experience of heat, namely air temperature, radiant temperature, humidity, and air velocity 15 . As performing work generates heat, in a high WBGT environment labour must be reduced in order to maintain a safe body temperature. The ISO standard and various national regulations set limits on continuous working time at given thresholds of WBGT. Workers self-pace in order to cope with heat strain, so changes in labour capacity are a proxy for welfare, as well as being an indication of an economic impact.
Occupational heat strain depends on not only meteorological aspects, but also those factors relating to the workers themselves; including individual differences in acclimatization, work intensity, clothing, and hydration 16 . This means that estimates of heat strain are specific to an occupational context, and suggests that field data collected in specific occupational context may not be generally applicable.
Previous studies 6,[17][18][19][20] have weighted exposure to heat hazard based on human population maps, and examine either the hottest month of the year, or the average of the whole year. Rice harvests occur in certain times of year, depending on location. These are not always the hottest months of the year; however, these are months in which we know that heavy physical labour is performed. We weight the result using rice harvest maps, assuming that harvested weight is a proxy for the number of workers engaged in harvesting, in order to quantify the labour capacity effect of global warming on rice harvest workers specifically.

Methods
In this study we analyse atmospheric data from the Coupled Model Intercomparison Project, Phase 6 (CMIP6) 21 , which is 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, humidity, and surface level pressure. All CMIP6 models for which appropriate data were present in the Centre for  25 , were used as the historical observational datasets.
WBGT is designed to be measured directly using specialised equipment; however, empirical formulae for estimating it from standard meteorological variables do exist. Daily minimum and maximum wet bulb temperature (WBT) was calculated from daily mean and max near surface air temperatures, relative humidity, and pressure using the open source software 'psychrolib' 26 which implements formulae from the American Society of Heating, Refrigerating and Air-Conditioning Engineers handbook. In our WBGT calculation, we neglect the effects of solar irradiance, by assuming that the black globe temperature is approximated by the air temperature. This means that we assume work is occurring in the shade, and therefore underestimate WBGT by several degrees in sunny conditions 27 . Properly accounting for solar irradiance would require either working with sub-daily data, or making further assumptions about sub-daily variation. However, in this study we focus long term trends, which are driven by air temperature and humidity; as there is no clear long-term trend in surface downwelling shortwave flux, it does not affect our results.
Field measured WBT decreases with wind speed at low speed (<2 m/s, a light breeze), but higher wind speeds have a lesser effect 27 . The WBT calculation we used assumes that the wick is well ventilated, so variable wind speed is safely neglected.
Studies of occupational heat stress and strain under climate change often either assume a threshold above which a worker is at risk 18 , or assume that there is a simple relationship between WBGT and worker productivity. These relationships are assumed to be representative across sectors, and are based on either a regulatory or advisory standard 17 , limited field study or survey data, or an ad hoc fusion of the two 6,1928 .
Standards attempt to minimise harm, so can appear conservative when compared to the actual conditions in which people work. Field data for the effect of heat on worker productivity are relatively sparse, and cover only a few activities: studies of small numbers of workers are used to estimate productivity effects on the entire human population.
Sahu et al 10 observed a 5% per °C WBGT decrease in the labour capacity of labourers harvesting rice between 23-33 °C WBGT. Rate of collection was measured in 124 workers in groups of 10-18, and WBGT was measured in-situ, at a single location in India. We adopt this for our labour impact metric, and assume that this is representative of manual rice harvest labour. The impact is assumed to be linear in WBGT, although this assumption must break down as WBGT approaches human skin temperature. The systematic uncertainty due to these assumptions cannot be assessed without larger scale field observations. The labour loss function is -5.14*WBGT + 218, in units of %, clipped at 0 and 100. This means that 0% loss occurs at 23 C and 100% loss occurs at 42.5 C.
Daytime temperature variation is estimated by assuming the temperature is close to the daily mean temperature, daily max temperature, and the mid-point of the two for 4 hours each, following previous studies 6 . Comparing the result of this calculation using hourly ERA5 and daily ERA5 data, we found that this is a reasonable approximation (see Supplementary Material).
The RiceAtlas dataset 1,2 provides detailed data for 31 countries in Asia, broken down into location entities with an average area of 5000 km. 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. The harvest season is typically around 30 days long in this database. Figure 1 shows the distribution of rice cropped land in Asia, from the RiceAtlas dataset, with some key locations labelled.
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 enclosed multiple GCM grid-cells, the mean is taken of the enclosed cells; otherwise, the land grid-point closest to the centroid of the location was used.  Linear regression between the rice harvest labour impact and global mean surface air temperature (both as 20 year means), performed for individual harvest location-seasons, shows that a linear relationship is present at the local level too. Figure 4 shows the production-weighted histogram of the impact gradient (multi-model mean). There is a high level of agreement between models, with an average correlation coefficient between models of 0.95. Harvest location-seasons fall into two groups: those strongly affected by climate change (impact gradient ~ 4 %/C) and those that are not (impact gradient ~ 0 %/C, with large multimodel spread). The global mean (Fig 2) obscures this fact, as the average is over heterogeneous data. Fitting a normal distribution to the peak, we see that the strongly affected group is centred at 3.9 %/C, with a standard deviation of 0.3 %/C: we therefore identify the strongly affected group as all harvest locationseasons having a hazard gradient >3.0 %/C. The harvest location-seasons were split into two groups, split at 3 %/C. The high hazard gradient group comprises 50% of production in Asia. For locations having multiple harvest seasons, the proportion of production in the location that is in the more exposed group was calculated, and plotted in Figure 5. The proportion of production in each country that is exposed is listed in Table 1.

Results
Exposure is generally very high in in Southeast Asia. All production in Indonesia, Sri Lanka, and Malaysia is exposed. Vietnam's production is 84% exposed, with 100% exposure in the Mekong River Delta and less in the north.
Generally, harvests that peak in September, October, and November (SON) are the least exposed due to lower seasonal temperatures. Locations at higher latitudes also on average have lower seasonal temperatures. This largely explains the spatial variation of exposure in China. Jiangsu's harvest peaking in October is not exposed, whereas Chongqing's harvest peaks in August so is 100% exposed; multiple cropping in leads to partial exposure in Hubei, Jiangxi and Hunan. Across China, 27% of production is exposed.
All locations in Bangladesh have an exposed harvest, but overall exposure is 61% due to multi-cropping: harvests in March-May are exposed, harvests in November-December less so. Although in India, a large amount of production is exposed in West Bengal and Andhra Pradesh, production further north in Punjab and Uttar Pradesh is not identified as exposed, meaning that the overall exposure of India is just 36%. The rice harvest in Punjab and Uttar Pradesh is in October-December, a cooler time of the year; West Bengal has rice harvests in between April and July.
Due to the resolution of climate models, results are not reliable for locations closest to the Himalayas. Temperatures are averaged across a grid-cell that includes both mountains and plains, but rice agriculture mainly occurs in plains and valleys.
In some cases, especially wheat in Punjab, rice is multi-cropped with another crop not covered in this study, which may have its harvest exposed. Studies have suggested that changes to planting and harvesting dates 29 , and crop choice 30 , could be used to adapt to climate change and mitigate negative effects on yields.
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 taken into account when planning adaptive measures.
These results are comparable to those of similar previous studies, but with a different focus. Previous studies either average over the whole year or focus on the hottest part of the year. However, different parts of the year are not equally important in terms of agricultural labour. Planting and harvesting typically take up the most labour, and they cannot be delayed or displaced.
Orlov et al 31 project a decrease in work productivity of approximately 9.5% and 12% for South Asia and Southeast Asia respectively, in the crops sector by 2050 under the CMIP5 high-emission scenario (their figure   9). Assuming this corresponds to 2 C of warming relative to the present gives an impact gradient of 4.7 and 6.0 %/C respectively. Wind and surface downwelling shortwave radiation was included in the calculation of the WBGT, but the calculation used daily average irradiance withthe Sun assumed to always be at its zenith.
The general picture and conclusions of our work are unchanged when the labour impact function is changed to the one used by Orlov et al for the crops sector. A function based on fitting a curve to an ad-hoc combination of data from the Sahu et al study and a study of miners, and has a similar gradient from 29-32 C WBGT.
In the study by Dunne et al 17 , much more dramatic reductions in labour capacity were projected, as industrial safety guidelines for WBGT were used as the labour impact function. If the same function is applied in our analysis, a strong correlation is still seen between the calculated impact and the GSAT, but with a higher gradient. Clearly, the choice of labour impact function is very important and has a strong impact on the final result. The labour impact function used by Dunne for heavy labour reaches 100% at 33 C WBGT, at which WBGT the function from Sahu et al is 52%.
Our study is limited by the accuracy, coverage, and granularity of the rice crop data. However, by focussing on long-term trends with high levels of multi-model agreement, we largely 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.

Conclusion
The growing impact of heat on rice harvest labour will be very unevenly distributed, mainly falling on poor rural workers who will be least able to adapt, and therefore contributing to widening economic inequality. Given agriculture employs hundreds of millions of people in the region, a few-percent shift in their labour capacity is equivalent to the labour of millions of people.
Vulnerability may be higher in those locations with low levels of farm mechanization, but we note that high quality international data on the level of farm mechanization are lacking. Farm mechanisation in Asia over the past 50 years has been largely driven by small internal combustion engines, which have already alleviated much of the human drudgery of farming, reduced the use of draft animals, and increased yields 32,33 . However, more fatal injuries occur when using powered machinery 33 , and the use of internal combustion engines contributes to greenhouse gas emissions. Labour productivity decreases due to climate change will compound other environmental impacts on rice production, including declining yields as direct result of increasing temperature, water stress, and sea level rise 12 ; contributing to the unequal loss and damage created by climate change.
We see a strong relationship between heat hazard and global mean surface temperature change in the rice harvest seasons and locations comprising 50% of the Asian harvest, with a high level of agreement between climate models. Historical observational data shows already shows a statistically significant increase in our labour impact metric. Understanding disparities in the exposure of workers and industries around the world will be important for climate adaptation strategy, as well as estimation of loss and damage. Overall, the exposure of such a large proportion of rice agriculture, and the workers engaged therein, provides an argument for strong mitigation of climate change. Table 1: Rice production exposed to high hazard gradient by country, and Asia total. Countries with less than 1 million tonnes of annual rice production are not included. Rice production from RiceAtlas 1 .

Data and code availability statement
The datasets generated during the study will be available from the CEDA archive, and accession codes provided before publication. Computer code will be publicly available under the MIT license before publication.

Author contribution statement
CS performed the analysis and mostly wrote the paper. JSH, DM, RAB, and ES helped develop the ideas for the analysis, and advised on the wording and emphasis of the paper.