1 Introduction

Climate change is happening and affects all aspects of agricultural systems and consequently food security in the future (Schmidhuber and Tubiello 2007). This is especially true in West Africa (WA), a region which is highly vulnerable to climate change due to low adaptive capacity (FAO 2007). On the basis of the latest scenarios of the Intergovernmental Panel on Climate Change (IPCC), the average temperature in this region for the year 2100 is expected to increase by 2 to 6 °C (Sylla et al. 2016) under the Representative Concentration Pathway (RCP) 4.5 and 8.5 scenarios, respectively, i.e., a temperature increase higher than the global average (Pachauri et al. 2014; Turco et al. 2015). Furthermore, Sylla et al. (2018) showed that in the future, most human population in the Sahel and the Western Sahara Desert would be at risk of experiencing heat stress conditions.

There is growing evidence that, among all agricultural sectors, climate changes will hit the livestock sector in developing countries (Thornton et al. 2009), like WA, where the sector accounts for about more than one third of agricultural Gross Domestic Productivity (GDP) (Zougmoré et al. 2016). Previous studies also demonstrated that temperature and humidity changes above the comfort zone will affect livestock in several ways. It influences their production and performance and even increases mortality rate depending on species, genetic potential, life stage, and nutritional status (Hahn 1999; Kadzere et al. 2002; West 2003). For example, heat stress decreases the forage intake by 18% and decreases milk yield by 32% at 32 °C compared with 16 °C (Chase 2006).

The most commonly used method for assessing heat stress is the temperature-humidity index (THI), which incorporates the effects of both relative humidity and temperature to describe the level of heat stress for livestock. This method was originally proposed by Thom (1959) to evaluate the degree of heat stress for humans. Later, different equations were developed for different livestock species that could be used in heat stress impact studies. During the past decade, several studies tried to use different THI indices to assess heat stress conditions on dairy cattle under future conditions across global and regional scales (e.g., Hernández et al. 2011; McCabe et al. 2016; Gunn et al. 2019). For example, Lallo et al. (2018) investigated heat stress conditions for different livestock in the Caribbean and concluded that ruminants will experience higher heat stress under a future climate as compared with present-day conditions, and that this will threaten their productivity. For South Africa, Nesamvuni et al. (2012) investigated the future heat stress for dairy cattle using three GCMs (ECHAM5/MPI-OM, MRI-CGCM2.3.2, and GFDL-ESM2M-CM2.0/2) under the IPCC’s A2 emission scenario and showed that dairy cattle are expected to experience severe stress conditions in most part of the country by 2081–2100.

Despite the economic importance of livestock production for West African countries and in view that climate change is predicted to be severe in this region, there is no information available on how and to what extent climate change might result in changes of heat stress condition in this region. Therefore, since heat stress is considered to be one of the most important variables affecting milk production, the current study aims to fill this gap and assess current conditions of heat stress and its trend during the historical period (1981–2010). We also aimed at projecting future (2021–2050 and 2071–2100) heat stress conditions for dairy cattle in WA by using three general circulation models (GFDL-ESM2M, HadGEM2-ES, and MPI-ESM-MR) under the RCP4.5 scenario.

2 Materials and method

2.1 Study area

The study area of this research is WA, which lies between 15° E–16° W and 4° N–25° N and consists of 14 countries (Benin, Burkina Faso, Côte d’Ivoire, The Gambia, Ghana, Guinea, Guinea-Bissau, Liberia, Mali, Niger, Nigeria, Senegal, Sierra Leone, and Togo) (Fig. 1). The climate of WA varies from the humid tropical monsoon conditions in the south to the arid desert conditions in southern Sahara and northern Sahel. The most important factors characterizing these climates are the altitude, the location of the ITCZ (Inter-Tropical Convergence Zone), the proximity of the tropical Atlantic Ocean, and the location of cyclonic (low atmospheric pressure) and anti-cyclonic (high pressure) circulation centers (Hayward and Oguntoyinbo 1987). The temperature conditions in WA are mostly influenced by latitude, where the maximum temperatures and temperature ranges increase from south to north. In the near-equatorial latitudes, the average temperatures do not vary much, and the diurnal cycles are also relatively low, with the mean annual range around 3 °C. This increases northward, with a mean annual range in the Sahel region in the order of 15 °C. The RH in humid climate zones is above 80% (often near 100%) in the early morning and is around 60–70% in the afternoon throughout the year. In semi-arid climate zones (e.g., Sahel), the RH is around 40 and 90% during the dry and wet season, respectively. In the afternoon, the RH could be even lower; in such a way that during the dry season, it is around 20% (or less), and during the wet season, it is around 60–70% (Nicholson 2018).

Fig. 1
figure 1

Map of West Africa displaying climate types and location of meteorological stations used in the study

The location of the study area, as well as its climatic zones, is shown in Fig. 1.

2.2 Data sources

2.2.1 Observational data

In the present study, maximum temperature (Tmax) and relative humidity (RH) observations from 50 synoptic weather stations across WA were used for testing the accuracy of simulated climate by global climate models and also for testing uncertainty of different THI indices. Weather station data were gathered from the archive of the WASCAL dataset (West African Science Service Center on Climate Change and Adapted Land Use), NOAA-NNDC (http://www7.ncdc.noaa.gov/CDO/dataproduct), and TuTiempo.net database (http://en.tutiempo.net/).

As is shown in Fig. 1, the spatial distribution of the synoptic stations used in this study is relatively uniform and can represent different climate zones in WA.

2.2.2 Simulation data

The present study makes use of the WASCAL’s high-resolution regional climate simulations (at 12 km resolution) performed by using the Weather Research and Forecasting Model (WRF, version 3.5.1), driven by the lateral boundary conditions from the ERA-Interim re-analysis data (Dee et al. 2011), the GFDL-ESM2M (Dunne et al. 2012), the HadGEM2-ES (Jones et al. 2011), and the MPI-ESM-MR (Stevens et al. 2013). The climate simulations are based on the RCP4.5 scenario, which is an intermediate scenario between optimistic and pessimistic scenarios (Van Vuuren et al. 2011; Thomson et al. 2011). More details about the overall concept of the WASCAL regional climate simulations, as well as detailed information on the WRF setup and its accuracy for simulating different climate characteristics, are given in Heinzeller et al. (2018). The WRF downscaling datasets have been evaluated against observations for the time period 1981–2010, and we discuss the expected future climate change impact for the two future 30-year periods, from 2021 to 2050 (near future) and from 2071 to 2100 (far future).

2.2.3 Livestock data

The Gridded Livestock of the World (GLW) database, as described in detail by Robinson et al. (2014), provides modeled livestock densities of the world, adjusted to match official (FAOSTAT) sub-national estimates, at a spatial resolution of at 1 km resolution for cattle, pigs, and chickens, and a partial distribution map for ducks. The cattle density layer was downloaded via the Livestock Geo-Wiki (http://www.livestock.geo-wiki.org).

About the dairy production systems in WA, it should be noted that they almost follow the same pattern as the agro-ecological zones and vary from the humid coastal zone in the South to the arid zone in the North. However, generally, traditional pastoral production systems (like migratory, transhumant, or sedentary) are dominating dairy production in this area (up to 90% of dairy cattle population (Mbogoh 1984; Olaloku and Debrah 1992)) and still supply the majority of the milk. These systems are especially important in Burkina Faso, Gambia, Mali, Niger, and Senegal, where they raise their cattle under extensive pastoral systems, the productivity of which are highly depending on climatic conditions. On the other hand, the commercial stall feeding production systems, which can still be influenced by heat stress conditions, have gradually increased in the southern parts of WA (like Northern Guinea Savanna in Nigeria) and become more dominant during the past decades (Thorpe et al. 2012). Furthermore, about the breeds commonly used for dairy production in WA, there are several pieces of research that have focused on cattle population for different countries and designed an inventory to identify the most predominant breeds within that area (e.g., Belemsaga et al. 2005; Santoze and Gicheha, 2019). According to these surveys, the “humpless” breeds (Bos Taurus) and “humped” (Bos indicus) breeds of cattle with their crosses are dominant cattle breeds in this region.

2.3 Testing reliability of climate data

In the present study, we used the observational daily climate data during the 1981–2010 period to test the ability of ERA-Interim to simulate Tmax and RH because we aimed to use this dataset as our historical dataset. Furthermore, the accuracy of GFDL-ESM2M, HadGEM2-ES, and MPI-ESM-MR global circulation models (GCM) was tested against the observational data during the historical condition in order to test their reliability for future climate projections.

To do this, we used Taylor’s diagram (Taylor 2001), which gives an overview of the statistical relations between two datasets, an estimate and a reference (REF) (by using root-mean-square deviation, spatial correlation, and normalized standard deviation (the standard deviation of simulations divided by the standard deviation of the observation)). In this diagram, each point represents the performance of each simulation, and their ability to simulate different parameters can be determined by comparing their distance to the REF point.

2.4 Thermal heat stress assessment

In this paper, we used THI as it is the most commonly used heat stress index globally. Many different formulae, however, exist for THI (see Table 1). For decades, many researchers tried to test the accuracy of the THI formula for different regions (e.g., Bohmanova et al. 2007; Berman et al. 2016; Kumar et al. 2018). For example, Dikmen and Hansen (2009) using rectal temperature of lactating cows in some sub-tropical region showed that the correlation coefficient between THI1 and rectal temperature was around 0.51. Furthermore, Jeelani et al. (2019) investigated the suitability of THI1 for measuring heat stress in crossbred dairy cattle in moderate to hot conditions and concluded that there were strong correlations between the THI1 and physiological, biochemical parameters, especially when the THI1 reached and crossed 80.

Table 1 List of THI equations used in this study and the agreement between different approaches for the detection of heat stress categories (%). (in this table, 100 shows the result of each THI that is always in perfect agreement with itself)

Since there was no THI equation developed for different cattle breeds and climates which exist in WA, we first explored the results of different THI equations for different climate zones over WA. To this end, widely used THI formulae for dairy cattle were gathered and compared during the historical period. Table 1 shows the list of the six THI equations explored for further use in this study. Following the results of previous studies, which have shown that the THI1 gives satisfactory performance for different regions, in the present investigation, the daily THI values of each grid cell in WA were calculated using this method.

Following the use by Vitali et al. (2009), Segnalini et al. (2013), and Lallo et al. (2018), in this study, Tmax was used for the dry bulb measurements and ambient temperature (Ta). Also, wet bulb temperature was calculated based on Wagner and Pruß (2002), and dew point temperature was calculated based on Stull (2011). Another reason for calculating these variables by using temperature and RH was the absence of measurements for these variables in both historical and future datasets.

After calculating the THI value for each day, the daily THI values were divided into four THI subgroups (Dash et al. 2016), using the following THI thresholds proposed by Wiersma (1990): thresholds none (THI < 72), mild (72 ≤ THI < 79), moderate (79 ≤ THI < 89), and severe/danger (THI ≥ 89). It is worth mentioning that thresholds used for identification of dairy cattle’s physiological responses to heat stress condition were mainly developed based on experiments for Holstein-Friesian breeds, and that the transferability of the methodology can be discussed. However, information for other breeds is currently not available. This means that the results presented in this study might be less severe for highly heat stress adapted breeds or cross-breads. Table 2 presents the thresholds used for subdividing the daily THI values into four heat stress categories and the expected response in dairy cattle.

Table 2 The thresholds used for categorizing the THI values based on their response in dairy cattle (based on Wiersma 1990)

As the last step, the probability of each THI subgroup during a year has been calculated for both historical (1981–2010) and future condition (2021–2050 and 2071–2100) under the RCP4.5 scenario. Furthermore, the THI trend during the historical period as well as the longest severe and danger subgroups has been investigated.

3 Results

3.1 The accuracy of downscaled data for simulating Tmax and RH during the historical period (1981–2010)

Figure 2 (right) depicts the Taylor diagram for daily Tmax indicating the spatial centred root-mean-square differences, normalized standard deviations, and correlation coefficients between the dynamically downscaled and observed daily Tmax in WA. As shown in this figure, all GCM models used in this study could reliably reproduce the observed daily Tmax during the historical period (1981–2010). Based on the obtained results, the spatial correlation coefficients between the observed and dynamically downscaled Tmax data range from 0.75 (GFDL-ESM2M) to 0.88 (ERA-Interim), and all simulations were statistically significant at the 95% confidence level. Centered root-mean-square differences range from 0.50 (ERA-Interim) to 0.71 (GFDL-ESM2M). The normalized standard deviations range from 0.87 (GFDL-ESM2M) to 1.09 (ERA-Interim), which shows that, except for ERA-Interim, the standard deviation of all GCMs (GFDL-ESM2M, HadGEM2-ES, and MPI-ESM-MR) are less than observed values.

Fig. 2
figure 2

Taylor diagrams of daily Tmax and RH simulations for the ERA-Interim, GFDL-ESM2M, HadGEM2-ES, and MPI-ESM-MR models as compared with the observations over West Africa

For RH, as shown in Fig. 2 (left), the spatial correlation coefficients, the centred root-mean-square differences range, and the normalized standard deviations of ERA-Interim, GFDL-ESM2M, HadGEM2-ES, and MPI-ESM-MR were nearly the same, and all four dynamically downscaled datasets were statistically significant at the 95% confidence level. The only point about these simulations is that the standard deviation of all simulations was a little higher than observed values.

From the results of this analysis, it could be concluded that the ERA-Interim, GFDL-ESM2M, HadGEM2-ES, and MPI-ESM-MR models showed a satisfactory ability to simulate the Tmax and RH for WA. Therefore, in this study, the results were presented and discussed for the historical period (1981–2010), based on downscaled ERA-Interim data, and expected changes under future climate conditions (2021–2050 and 2071–2100) are derived from the GFDL-ESM2M, HadGEM2-ES, and MPI-ESM-MR GCMs under the RCP4.5 emission scenario.

3.2 Comparison of different THI equations in different climate zones in WA

The THI equations used to quantify thermal heat stress conditions differ in the way they estimate humidity as some use wet bulb temperature (Tw), dew point temperature (Tdp), and or relative humidity (RH) for calculating THI. Moreover, equations also use different coefficients for temperature and humidity (see Table 1). For this reason, before using a specific THI equation, we compared different THIs with regard to the occurrence of heat stress category for the historical period (Table 1). The Table shows that results of the different calculation methods are in agreement with each other (for THI1 and THI2, THI2 and THI3, THI1 and THI3 more than 80%), with a few exceptions (e.g., THI4 and THI6 around 25% agreement to the other calculation methods). Based on these findings, and taking into account that THI1 which has been proposed by NRC (1971) is the most widely used method to analyze thermal heat stress condition for dairy cattle in tropical environment (Nascimento et al. 2019), we based all other calculations on this index.

3.3 Historical heat stress conditions in WA

The main results of studying the geographical distribution of the average heat stress frequencies for none, mild, moderate, and severe/danger conditions in WA during the historical period (1981–2010) are shown in Fig. 3. Based on the results, it can be concluded that during the 1981–2010 period, the moderate heat stress category was dominant in most of WA. The frequency of moderate heat stress category ranged between 8% (~ 29 days/year) in eastern parts of Nigeria to 69% (~ 251 days/year) in the southern half of Senegal. Across the entire WA, the average frequency of moderate heat stress conditions during the historical period was 47% (~ 171 days/year).

Fig. 3
figure 3

Frequency of different THI categories during the historical period (1981–2010)

The frequency of the none heat stress category, i.e., conditions at which dairy cattle have their optimum productive and reproductive performance, is around 2.5% (~ 9 days/year) in the near-coastal regions (around latitudes 4–6° N), but increases northward and peaks in the Sahel region where it is on the order of 45% (~ 164 days/year). In contrast, mild heat stress conditions are most prominently found in the southern half of WA (about 65% of the days during a year) (Fig. 3).

During the historical period of 1981–2010, the frequency of severe/danger heat stress categories ranges between 17% (~ 62 days/year) in the central parts of WA (Mopti, Gao, and the southern half of Tombouctou regions in Mali; Tillabéry and the northern half of Tahoua regions in Niger) to less than 1% (less than ~ 3 days/year) in the rest of the study area.

To this end, we have conducted a trend analysis to test if the frequency of severe/danger heat stress has increased or decreased over the 1981–2010 period. Figure 4 shows that over the historical period of 1981–2010, the frequency of severe/danger heat stress shows that about 21.7% of the study area shows a significant (95% confidence level), an increasing trend. Significant increasing trends in heat stress are most obvious in the southern countries, e.g., Nigeria, Benin, Togo, and Ghana.

Fig. 4
figure 4

Regional trends in the frequency of severe/danger heat stress events during the historical period of 1981–2010 (stippling indicates trend values significant at 95% confidence level)

3.4 Future heat stress condition (2021–2050 and 2071–2100)

The main results of studying the probable changes in heat stress conditions using three GCMs (GFDL-ESM2M, HadGEM2-ES, and MPI-ESM-MR) for two future periods (2021–2050 and 2071–2100) and under the RCP4.5 emission scenario are presented in Fig. 5, Fig. S1, and Fig. S2.

Fig. 5
figure 5

Frequency of different THI categories by 2021–2050 and 2071–2100 periods projected (average of three GCMs: GFDL-ESM2M, HadGEM2-ES, and MPI-ESM-MR) under the RCP4.5 scenario

For the 2021–2050 climate period, the results obtained by three different GCMs indicated that the general pattern of changes is similar for all heat stress categories. Though, heat stress periods for dairy cattle are slightly worse for predictions based on the MPI-ESM-MR model as compared with those based on climate predictions by GFDL-ESM2M and HadGEM2-ES models.

In all scenarios, the moderate heat stress category is still dominant, with frequencies increasing most in the near-coastal regions (around latitudes 4–6° N). This increase in the frequency of moderate heat stress periods is pronounced for HadGEM2-ES projections. Again, for the near-coastal regions, the frequency of the severe/danger heat stress category is expected to increase to 6.4% (~ 23 days/year), 8.2% (~ 29 days/year), and 16.5% (~ 60 days/year) of the year based on GFDL-ESM2M, HadGEM2-ES, and MPI-ESM-MR projections, respectively. In contrast, the average frequency of the severe/danger category during the historical period was approx. 5.5% in this area.

However, for the northern half of WA, heat stress conditions for cattle might decrease under future climate conditions. More precisely, in the northern parts of Mali (Tombouctou and Kidal regions) and Niger (Agadez Region), the average frequency of none heat stress category is expected to increase from 4.7% (historical period) to 4.9–5.7% (~ 18–20 days/year), though still this is a rather low number of days without heat stress.

For WA as a whole, by 2021–2050 period, the frequency of severe/danger heat stress category is expected to increase by 1.1% (~ 4 days/year), 2.7% (~ 10 days/year), and 11.1% (~ 40 days/year) as compared with current conditions based on GFDL-ESM2M, HadGEM2-ES, and MPI-ESM-MR projections, respectively.

For the 2071–2100 climate period, the overall patterns of changes of heat stress periods are in agreement independent of the GCM used, though calculations of THI based on MPI-ESM-MR show most pronounced changes. Figure 5 shows that by the end of twenty-first century, the frequency of dangerous heat stress conditions will increase, especially in the southern half of WA (pure tropical, transitional tropical, and transitional equatorial climate zones). Except for the northern parts of Mali (Tombouctou and Kidal regions) and Niger (Agadez Region), in other regions, the frequency of none heat stress category is expected to be < 5% (~ 18 days/year).

In the southern half of WA, the frequency of moderate heat stress category is expected to increase by 28% (~ 102 days/year, GFDL-ESM2M), 31% (~ 113 days/year, HadGEM2-ES), and 18% (~ 65 days/year, MPI-ESM-MR) by 2071–2100 period as compared with the historical period, respectively.

For WA as a whole, by 2071–2100, the average frequency of severe/danger heat stress category is expected to be 7.4%, 13.7%, and 22.2% based on GFDL-ESM2M, HadGEM2-ES, and MPI-ESM-MR projections, respectively. Furthermore, in the southern half of WA, there would be some area where the frequency of dangerous heat stress events (i.e., severe/danger classes) would be around 63.4% (231 days/year).

Since the livestock production and management systems in WA vary across different climate zones, results are presented and discussed for different climate zones (i.e., desert, semi-arid desert, semi-arid tropical, pure tropical, transitional tropical, transitional equatorial climate zones). In Fig. 6 and Fig. S3, the results of projected changes in the average number of days with severe/danger heat stress events in these climate zones in both historical (1981–2010) and future periods (2021–2050 and 2071–2100) are shown. Figure 6 shows also the range of uncertainties for the projected future climate due to different GCM formulations. It suggests that although all three GCMs predict the same changes in the average number of days with severe/danger heat stress events, there are still uncertainties associated with model outpus. However, according to obtained results, in desert and semi-arid desert climate zones, the average number of days with severe/danger heat stress events is expected to decrease from ~ 12 days in the historical period to ~ 4 and ~ 7 days/year by 2021–2050 and 2071–2100, respectively. Results of future projection of days with severe/danger heat stress for the semi-arid tropical climate zone differ in dependence of the GCM being used. While using MPI-ESM-MR, we see an increase (from ~ 29 to ~ 58 days) in heat stress in this climate zone, the opposite result, i.e., a decrease can be seen while using GFDL-ESM2M or HadGEM2-ES (from ~ 29 to ~ 8 and ~ 18 days, respectively).

Fig. 6
figure 6

Projected changes in average number of days and the length of consecutive days with severe/danger heat stress categories for future climatic conditions

For pure tropical, transitional tropical, and transitional equatorial climate zones, days with severe/danger heat stress increase in both future periods. This increase is expected to be greater for the 2071–2100 period. As an example, in the transitional equatorial climate zone, the average number of days with severe/danger heat stress condition might increase from ~ 10 days/year for historical conditions (1981–2010) to ~ 79, ~ 117, and ~ 159 days/year by 2071–2100, respectively.

For WA as a whole, by 2021–2050, the average number of days with severe/danger heat stress events is expected to increase from ~ 18 days during the historical period to ~ 19, ~ 27, ~ 62 days based on GFDL-ESM2M, HadGEM2-ES, and MPI-ESM-MR simulations, respectively. Furthermore, by 2071–2100, it is expected to have an average ~ 27, ~ 51, ~ 72 days with severe/danger heat stress based on GFDL-ESM2M, HadGEM2-ES, and MPI-ESM-MR simulations, respectively.

The average length of continuous periods during which the THI exceeds the severe/danger threshold is also considered to be an important factor for dairy cattle because longer exposure to dangerous heat stress conditions will result in higher animal production losses (Gantner et al. 2011; Pragna et al. 2017). Based on the results presented in Fig. 6 and Fig. S4, the average length of continuous severe/danger heat stress is expected to increase by the future period, especially in transitional tropical and transitional equatorial climate zones. For example, based on calculations using the climate predictions by the MPI-ESM-MR model, the average length of continuous dangerous heat stress periods is expected to increase from 12 days in the historical period to 50 or 65 days/year by 2021–2050 and 2071–2100, respectively. For WA as a whole, by 2021–2050, the average length of continuous dangerous heat stress periods is expected to change from ~ 3 days in the historical period to ~ 4–7 days and even to up to 10 days by 2071–2100.

3.5 West African dairy cattle population at risk

As discussed earlier, although in some regions, like the northern parts of Mali (Tombouctou and Kidal regions) and Niger (Agadez Region), the frequency of severe/danger heat stress categories is projected to decrease (by up to 10%); for other regions throughout WA, these frequencies are projected to increase. The largest increase is projected to happen in southern regions (e.g., southern parts of Côte d’Ivoire), where the results showed that the frequency of severe/danger categories might increase by more than 50%. According to the livestock population data of GLW 3 dataset (Gilbert et al. 2018) and FAOSTAT database (FAO 2013), the number of dairy cattle in WA has increased from 4 million to close to 10 million between 1981 and 2010. This population is expected to continue to increase by future periods to meet increasing demands for livestock products (Thornton 2010).

In order to estimate how many dairy cattle in WA are expected to be at risk with the expected increase in heat stress conditions, we overlaid the dairy cattle population data with the results of expected changes in the frequency of dangerous heat stress events in the future periods. The projected changes in heat stress frequency and proportion of affected dairy cattle population in WA are presented in Table 3.

Table 3 Dairy Cattle population at risk of change in frequency of severe/danger (dangerous) heat stress under future climate conditions (2021–2050 and 2071–2100) under RCP4.5 scenario

As presented in Table 3, depending on the GCM used, around 50, 61, and 73% of the current dairy cattle population are in places where the frequency of severe/danger heat stress is expected to increase by 2071–2100. More specifically, for example, based on MPI-ESM-MR projections under the RCP4.5 scenario, around 5.3% (by 2021–2050) and 9.8% (by 2071–2100) of cattle are living in places where the frequency of severe/danger heat stress is expected to increase by more than 40%.

3.6 Heat stress related impacts of climate change on Milk production in WA

Heat stress negatively affects the total milk production at different lactation stages (Aggarwal and Upadhyay 2012). In literature, it is often assumed that the milk production level stays almost constant until a certain threshold (normally until the threshold of the moderate heat stress condition) and then it starts to decrease linearly with increasing degree of heat stress (Ravagnolo and Misztal 2000; Bohmanova et al. 2007). Since different dairy cattle breeds have different susceptibility to heat stress depending on the level of milk production (Novak et al. 2009; Herbut and Angrecka 2012; Brügemann et al. 2012), we based our assumption on studies for the most predominant breeds within the West African region (i.e., Bos indicus and Bos taurus). It is worth noting that the milk production of Zebu type is naturally low (around 1000–1300 kg per lactation for pure Zebu cattle (Yilma et al. 2006) and 2700 kg per lactation for Zebu × Friesian crossbreds (Ahmed et al. 2007)) and depends on regions and production conditions.

In this study, we estmated the potential impact of heat stress on milk production based on the results of the study done by Santana et al. (2015) who found that milk yield above the moderate heat stress threshold declines by 0.099 kg milk/day per THI unit for Zebu cattle under a pastoral management system. Figure 7 shows that depending on the GCM used, milk production is expected to decrease by 200–400 kg/year by 2071–2100 in around 1, 7, or 11% of the study area. The most severe changes in milk production in future periods are projected to occur in the southern half of WA (latitudes 4–14° N).

Fig. 7
figure 7

Expected reduction in milk production due to changes in the frequency, intensity, and duration of moderate and severe/danger heat stress under future climate conditions (RCP4.5 scenarios for 2021–2050 and 2071–2100)

4 Discussion

The aim of this study was to examine frequency and trends of heat stress during the historical (1981–2010) and two future periods (2021–2050 and 2071–2100) in WA using a temperature-humidity index (THI) for categorizing heat stress periods on basis of measured/predicted daily temperature and humidity values. For this purpose, THI values have been calculated for both the historical and future conditions (under RCP4.5 scenario), using daily Tmax and RH data from 50 meteorological stations and four dynamically downscaled datasets (i.e., ERA-Interim, GFDL-ESM2M, HadGEM2-ES, and MPI-ESM-MR models downscaled by WRF). The following conclusions can be drawn from this study:

  • Investigating the accuracy of ERA-Interim, GFDL-ESM2M, HadGEM2-ES, and MPI-ESM-MR models for simulating Tmax and RH during the historical period (1981–2010) indicated that all dynamically downscaled datasets were able to capture these parameters with reasonably good accuracy. Heinzeller et al. (2018) also investigated the accuracy of these datasets over WA and showed that ERA-Interim downscaled by WRF agrees with the temperature observations. Also, with regard to the results of three GCM-driven runs (i.e., GFDL-ESM2M, HadGEM2-ES, and MPI-ESM-MR), MPI-ESM-MR downscaled by WRF fits the observed temperatures better than other simulations.

  • Comparing the results of six different THI equations in different climate zones in WA indicated that in some cases (like THI1 and THI2), the results of different calculation methods were in agreement with each other, and their differences were not statistically significant. Bohmanova et al. (2007) by comparing seven different THI equations in semi-arid and humid climates in the USA indicated that although they differ in scale, they follow a similar pattern and the largest differences between equations were observed in summer months.

  • Investigating the historical heat stress condition in WA indicated that, during the historical period of 1981–2010, most regions already have experienced an increasing trend of dangerous heat stress for dairy cattle, especially in the eastern and southeastern parts of WA.

  • Calculating the average frequency of heat stress categories in WA based on GFDL-ESM2M, HadGEM2-ES, and MPI-ESM-MR projections indicated that the frequency of dangerous heat stress events is expected to significantly increase under future climate conditions, especially in the southern regions. This issue is considered to be critical for pastoral and agro-pastoral producers in the north (who contribute in around 70% of milk production), as they historically move their cattle from the north to the more humid farming areas of the south in search of grazing area and water (Kamuanga et al. 2008). However, increasing the frequency of heat stress in the southern half of WA due to climate change would be considered as a new challenge for milk production.

  • The future climate scenario analyses indicate that the length of consecutive days of heat stress for dairy cattle with THI higher than the severe/danger threshold is likely to increase in future periods in WA. These changes are also important because the severity of heat stress impact on milk production would increase by increasing the length of heat stress period (Rojas-Downing et al. 2017; Krishnan et al. 2017).

  • Based on GFDL-ESM2M, HadGEM2-ES, and MPI-ESM-MR projections under the RCP4.5 scenario, by 2071–2100, around 2.0, 5.8, and 9.8% of dairy cattle population in WA would experience dangerous heat stress events 40% (~ 145 days) more frequent than the historical period, respectively.

  • According to our results, milk production in WA will be severely challenged by increasing the frequency, intensity, and duration of moderate and severe/danger heat stress in the future. Assuming a decline of 0.099 kg milk/day when above the moderate heat stress category, we estimated that, by 2071–2100 period, the milk production per cattle is expected to decrease up to 100, 200, and 400 kg in around 21, 27, and 11% of the study area based on MPI-ESM-MR projections. These expected changes highlight the urgent need for an adaptation of the dairy production system in this region.

  • Although our focus in this study has been on dairy cattle production system (since the thresholds used for classifying THI values were mainly proposed for this system), there exist many studies in the literature using the same THI formula and thresholds for other ruminants (e.g., for goat: Salama et al. 2014; for sheep: Srikandakumar et al. 2003; for beef cattle: Brscic et al. 2007; for buffalo: Dash et al. 2016). Therefore, the findings of the current investigation could also be extended to other livestock species in WA.

The probable changes in heat stress events in WA which have been discussed here could have significant consequences for dairy cattle production systems since heat stress exerts an important control on production quantity (mediated through reductions in growth, milk yield (Bouraoui et al. 2002; Gantner et al. 2011)), reproductive performance (Bilby et al. 2008), immune system (Das et al. 2016)), and milk quality (Pragna et al. 2017). Furthermore, previous studies indicated that changes in heat stress conditions would also affect feed intake by livestock (National Research Council. 1981; Könyves et al. 2017; Baumgard et al. 2006). As a next step, it is thus necessary to quantify the actual impact of heat stress changes on the livestock sector in terms of its productivity and associated economic and livelihood consequences. In addition, it will be necessary to assess the adaptative capacity and tolerance of existing breeds to heat stress and develop livestock management strategies for these countries to cope with heat stress (e.g., nutritional interventions, provision of shade, housing, fans, and sprinklers).

The results presented in this study were all based on RCP4.5, which assumes that our current feasible climate change mitigation strategies would lead to stabilization of radiative forcing at 4.5 W m−2 in the year 2100 (Thomson et al. 2011). It should, however, be noted that if the current trend in GHG emissions continues the RCP8.5 scenarios is more likely (Hayhoe et al. 2018). Thus, our study is probably a conservative estimate of the likely impacts of future heat stress for livestock production. However, the projected changes in heat stress conditions for dairy cattle across WA as outlined here can provide a good indication of trends in current and future conditions to manage the future risk of climate change and design future adaptation strategies.