Projection of irrigation water demand based on the simulation of synthetic crop coefficients and climate change

In the context of major changes (climate, demography, economy, etc.), the southern Mediterranean area faces serious challenges with intrinsically low, irregular, and continuously decreasing water resources. In some regions, the proper growth both in terms of cropping density and surface area of irrigated areas is so significant that it needs to be included in future scenarios. A method for estimating the future evolution of irrigation water requirements is proposed and tested in the Tensift watershed, Morocco. Monthly synthetic crop coefficients (K_c) of the different irrigated areas were obtained from a time series of remote sensing observations. An empirical model using the synthetic K_c and rainfall was developed and fitted to the actual data for each of the different irrigated areas within the study area. The model consists of a system of equations that takes into account the monthly trend of K_c, the impact of yearly rainfall, and the saturation of K_c due to the presence of tree crops. The impact of precipitation change is included in the K_c estimate and the water budget. The anthropogenic impact is included in the equations for K_c. The impact of temperature change is only included in the reference evapotranspiration, with no impact on the K_c cycle. The model appears to be reliable with an average r² of 0.69 for the observation period (2000–2016). However, different subsampling tests of the number of calibration years showed that the performance is degraded when the size of the training dataset is reduced. When subsampling the training dataset to one-third of the 16 available years, r² was reduced to 0.45. This score has been interpreted as the level of reliability that could be expected for two time periods after the full training years (thus near to 2050).

The model has been used to reinterpret a local water management plan and to incorporate two downscaled climate change scenarios (RCP4.5 and RCP8.5). The examination of irrigation water requirements until 2050 revealed that the difference between the two climate scenarios was very small (< 2 %), while the two agricultural scenarios were strongly contrasted both spatially and in terms of their impact on water resources. The approach is generic and can be refined by incorporating irrigation efficiencies.

1 Introduction

Water resources are scarce in semiarid areas and a major part is allocated to agriculture. In the southern Mediterranean region, irrigation allocation to agriculture represents 80 % of total water abstraction. It varies from 46 % in eastern countries up to 88 % in Morocco in 2010 (FAO, 2016). This percentage has been decreasing in most southern Mediterranean countries during the last decades in particular due to the limitation of available resources and the increase of the urban water demand. In parallel, the pressure on water resources led to what Margat and Vallée (2000) called a “post-dam era” or what Molle et al. (2019) called a “groundwater rush” where subterranean water is overused to satisfy the growing water demand. In recent years, overexploitation of groundwater has been facilitated by technological inventions, affordable cost of exploitation, and weak monitoring by authorities (MED-EUWI working group on groundwater, 2007). The overexploitation of aquifers can be observed in different countries over the Mediterranean watershed (Custodio et al., 2016; Le Goulven et al., 2009).

To satisfy the continuous increase in food demand associated with population growth, the agricultural sector has been asked to pursue its already initiated process of conversion toward agricultural intensification and above all towards a sharp increase in yields. This context goes hand in hand with the increase in food trade. The replacement of traditional crops by more financially attractive crops is already underway (Jarlan et al., 2016). In the “growth” scenario (which is more or less the actual trend), presented by Malek et al. (2018), the annual production of cultivated land increases by 40 % and the production of permanent crops increases by 260 %. In the “sustainable” scenario, annual crop production and tree production increase by 30 % and 38 %, respectively. As a consequence, irrigation water needs are expected to increase (Fader et al., 2016). The expansion and intensification of tree crops will also further rigidify the demand for agricultural water and increase the pressure on groundwater reservoirs (Jarlan et al., 2016) in order to keep tree crops alive during drought events (Le Page and Zribi, 2019; Tramblay et al., 2020). This study is carried on in the Tensift basin in Morocco, where the increase in the irrigated area and the intensification of irrigation during recent decades have caused a long-lasting drop in the groundwater table (Boukhari et al., 2015). A multimodel analysis of the area over the period 2001–2011 (Fakir et al., 2015) has shown that the groundwater falls from 1 to 3 m/year and that the cumulated groundwater deficit in 10 years (about 100 hm³/year since 2001) is equivalent to 50 % of the reserves lost during the previous 40 years. Among the main causes of this depletion is a reduction and higher irregularity of precipitation (Marchane et al., 2017) for crop growth and groundwater recharge, a reduction of snow water storage (Marchane et al., 2015), an increase and intensification of irrigated areas, and a progressive conversion to arboriculture due to national strategy. Since irrigation relies increasingly upon groundwater abstraction, questions are inevitably raised concerning the future of local agriculture and groundwater.

Therefore, in regions where the cropping density and surface area of irrigated areas are growing strongly (see, for example, the cases of China and India; Chen et al., 2019), it is necessary to make projections of the irrigation water demand with the actual trend in order to build alternative scenarios. In their review, March et al. (2012) defined scenario analysis as “internally consistent stories about ways that a specific system might evolve in the future. Narrative scenarios are plausible accounts of the future rather than forecasts”. Narrative scenarios also help identify the drivers of change and the implications of current trajectories as well as the options for action (Raskin et al., 1998). Scenarios are therefore halfway between facts and speculations in terms of complexity and uncertainty (van Dijk, 2012). They are commonly used as a management tool for strategic planning and for helping managers strengthen decision-making.

Various studies (Arshad et al., 2019; Lee and Huang, 2014; Maeda et al., 2011; Schmidt and Zinkernagel, 2017; Tanasijevic et al., 2014; Wang et al., 2016) have addressed the estimation of irrigation water requirement (IWR) of a region with an approach close to Eq. (1). A simplified balance between crop evapotranspiration (ET_c) and effective precipitation (P_e) is computed for each crop i and multiplied by the corresponding irrigated area. ET_c can be inferred from the crop coefficient method (K_c; Allen et al., 1998). Most of the effects of the various weather conditions are incorporated into ET₀, which accounts for the water demand of a reference crop, while K_c mainly accounts for crop characteristics. K_c varies according to four main characteristics: the crop height, the albedo of the crop and soil, the canopy resistance of the crop to vapor transfer (leaf area, leaf age and condition, and the degree of stomatal control), and the evaporation from soil, especially the soil exposed to solar radiation. Furthermore, as the crop develops, those different characteristics change during the various crop phenological stages. K_c values and stage lengths, for typical climate conditions, are considered to be well known for many crops and have been compiled in the FAO (Food and Agriculture Organization) look-up tables.

$$\begin{array}{}\text{(1)}& \begin{array}{rl}\mathrm{IWR}& =\sum _{i=\mathrm{1}}^n({\mathrm{ET}}_{{\mathrm{c}}_i}-P_{\mathrm{e}})\cdot {\mathrm{Area}}_i\\ & =\sum _{i=\mathrm{1}}^n\left[\left(K_{{\mathrm{c}}_i}\cdot {\mathrm{ET}}_{\mathrm{0}}\right)-P_{\mathrm{e}}\right]\cdot {\mathrm{Area}}_i\end{array}\end{array}$$

This approach is also a very convenient way to account for both climate and crop changes. On the one hand, P_e and ET₀ are obtained from meteorology or climatology. On the other hand, K_c values and stage lengths are taken from the tables. Most of the work consists of evaluating the future irrigated area of each crop and assessing the impact of more efficient irrigation techniques.

There is a significant amount of literature about the estimation of land use and land cover changes (Mallampalli et al., 2016; Noszczyk, 2018), with various techniques to estimate or predict them. Many land cover change approaches are based on the transition probability that was introduced by Bell (1974) and have been eventually connected to cellular automata to account for geographical interrelationships (Houet et al., 2016; Marshall and Randhir, 2008). A very interesting technique has been to combine the top-down (demand-driven) and bottom-up (local conversion) processes of land cover change by proceeding to a simplification of local processes (van Asselen and Verburg, 2013; Verburg and Overmars, 2009). Despite a huge bibliography both in climate change and land cover change, scenario analysis over the past 25 years has mostly focused on climate change projections, while the impact on land use and land cover has been neglected. Titeux et al. (2016) found that only 11 % of the 2313 studies analyzed have included both land cover and climate changes. Also based on a large review, March et al. (2012) have called this a “hegemony” of climate as a driver of change. Furthermore, the implementation of land cover change techniques appears to be tedious and does not account for the intensification of cropping patterns. The motivation of the present work is to take into account both land cover change and crop intensification in future scenarios of irrigation water demand, taking into account the impact of climate change.

The hypothesis is that the dynamics of land use change and intensity of the cropping patterns (hence irrigation water demand) can be reduced to a synthetic monthly K_c time series for every irrigated area of the studied region. Each synthetic time series would then account for the spatial variability of cropping patterns inside each irrigated area. Unless no sudden change occurs, a statistical model accounting for the monthly trend of K_c, the impact of rainfall, and the effect of saturating the land cover with tree crops should give an accurate fit that would allow for extrapolating into the next few decades. As the synthesis of K_c for separated irrigated areas would also decrease substantially the amount of information compared to a land cover approach, some information, like the amount of tree crops, should be retrieved back from the time series. This approach would prevent the need for tedious land cover classifications and reduce the difficulty of working with discrete values for developing future scenarios.

The objective of this study is to produce two scenarios of irrigation water demand for the different irrigated areas of the Tensift watershed in Morocco. One is the trend and the other one is an alternative scenario derived from a narrative scenario. To do so, a method for simulating and extrapolating K_c is proposed and is tested for the 2000–2016 period. The future time series of K_c is used with climate change scenarios to obtain an estimate of future irrigation demands. As the strategy of irrigation for tree crops is treated very differently from annual crops in the Tensift watershed, the yearly percentage of tree crops is also retrieved from the K_c time series.

The content of the article is separated into four parts and a conclusion. The first part describes the study area and the dataset. The methodology is then detailed in three parts: (1) K_c time series are obtained from remote sensing observations for the observation period, (2) a model is constructed to project the IWR based on simulating the evolution of K_c and integrating climate change, and (3) an alternative scenario is built from the trend to reflect a narrative scenario built by local water managers. In the Results section, the performance of the model and its ability to make projections are analyzed. The alternative scenario is compared to the simulated trend. In the Discussion section, other paths are discussed, in particular those which are related to irrigation management.

2 Study site and data

2.1 Study area: the Haouz Plain in Morocco

The study area includes the Haouz Plain, underlain by a large Plio-Quaternary alluvial aquifer. It is located within the Tensift watershed in southern Morocco around the city of Marrakech (Fig. 1). The plain is bordered to the south by the High Atlas range. The High Atlassian watersheds that peak at 4167 m a.s.l. receive about 600 mm/yr of precipitation of which a major portion falls as snow; 25 % of the streamflow is generated by snowmelt (Boudhar et al., 2009). The plain, which is located at altitudes ranging between 480 and 600 m a.s.l., has a semiarid climate (rainfall of 250 mm/yr and potential evapotranspiration of 1600 mm/yr) with mild, wet winters and very hot and dry summers. Irrigation is essential for good crop development and yield. The so-called “traditional” areas (purple), which originally received water from the Atlassian wadis, are more and more dependent on groundwater. The areas irrigated from lake reservoirs (green) receive water from the Lalla Takerkoust reservoir, the Moulay Youssef reservoir, and also from the Hassan I reservoirs through a 140 km canal called the Rocade Canal. The rest of the irrigation water comes from the groundwater.

Figure 1 The Haouz-Mejjat aquifer is located within the Tensift watershed (gray) and separated into four subareas. The most important traditional and dam irrigated areas are shown on the map. The N'Fis area is a mix of traditional and modern irrigation. All of these areas also use groundwater, but outside them, the irrigated areas (private irrigation) rely exclusively on groundwater.

2.2 Data

In order to compute the crop water demand, regional climate model (RCM) simulations over the Mediterranean region from the Med-CORDEX initiative (https://www.medcordex.eu (last access: 1 January 2017); Ruti et al., 2016) were utilized to provide future values of near-surface atmospheric variables, namely the 2 m air temperature and total precipitation. RCM results that are from two contrasting greenhouse-gas concentration RCP (Representative Concentration Pathway) trajectory scenarios were selected: the RCP4.5 scenario, which represents the optimistic scenario for CO₂ emissions (a stabilization of emissions) and the RCP8.5 (increasing emissions). Simulations from Centro Euro-Mediterraneo sui Cambiamenti Climatici (CMCC), Météo-France (CNRM), and IPSL-Laboratoire de Météorologie Dynamique (LMD) were used in the current study. The three RCMs give a good representation of the spread typically found in such simulations. The so-called “delta-change” or perturbation method (Anandhi et al., 2011) was used to downscale the RCM data using the Global Soil Wetness Project Phase 3 forcings (GSWP3: http://hydro.iis.u-tokyo.ac.jp/GSWP3, last access: 1 January 2017) as the reference historical data. Statistical scores (bias and daily correlation coefficients) were computed between the GSWP3 product and 11 SYNOP (surface synoptic observations) stations over northwest Africa over the 1981–2010 period. GSWP3 well represents the variability of temperature, relative humidity, and total precipitation. The correlation coefficients computed between both monthly time series are 0.99, 0.87, and 0.74 for temperature, relative humidity, and precipitation, respectively. Standard deviations are quite small for the three variables (0.01, 0.03, and 0.12, respectively). At the daily time step, GSWP3 well represents the day-to-day temperature variability and especially the relative humidity. However, GSWP3 has more errors in terms of the daily evolution of precipitation (r=0.08, SD = 14.6). In semiarid regions, where precipitation is infrequent and often convective, it is difficult to reproduce the precipitation events at the actual time. In terms of bias, the absolute temperature and relative humidity values are close to that observed: 0.22 ^∘C and 0.54 %, respectively; meanwhile, precipitation is underestimated by approximately 10 % on average, but this bias is not always consistent from one station to another. The 0.5^∘ data were resampled to a resolution of 1 km, applying a correction of altitude between the GSWP3 geopotential and a digital elevation model (GTOPO30 available at https://www.usgs.gov, last access: 1 January 2017) for air and dew point temperatures.

The resulting downscaled and rescaled (to 1 km resolution) dataset of meteorological variables (rain, temperature, wind speed, relative humidity and global radiation) extends from 2000 to 2050.

For generating vegetation time series, we used the MOD13Q1 Collection 6 product (Didan, 2015), which provide 16 d composite series of Normalized Difference Vegetation Index (NDVI) from the Moderate Resolution Imaging Spectroradiometer (MODIS) at a spatial resolution of 250 m. This product is computed from atmosphere-corrected (Justice et al., 2002), daily bidirectional surface reflectance observations, using a compositing technique based on product quality (Wan et al., 2015). We used the 16 d composites to reduce cloud coverage. The data have been compiled for the period 2000–2016.

3 Methodology

3.1 K_c assessment from remote sensing time series

Neale et al. (1990) proposed to estimate the K_c coefficient from vegetation indices obtained from remote sensing imagery. In this work, we have been using empirical linear equations where the slope (a) and intercept (b) have been previously calibrated for common crop type in local field experiments (Duchemin et al., 2006; Er-Raki et al., 2007, 2008):

$$\begin{array}{}\text{(2)}& K_{\mathrm{c},\mathrm{sat}}=a\cdot \mathrm{NDVI}+b.\end{array}$$

Equation (2) accounts for combined evaporation from the soil and transpiration from the crop in a very simple way, but as no water balance is performed, several assumptions must be kept in mind. (1) The extra evaporation due to wetting events (rainfall, irrigation) is not computed but is assumed to be represented in the calibrated NDVI $/$ K_c,sat relationship, which is related to the precipitation of the year used for calibration. The inability of this method to increase the evaporation part of evapotranspiration when frequent small rainfall events occur should be assessed and corrected for, either at the level of the ET_c computation or in the P_e computation. (2) The lower evaporation due to a reduction of the wetted fraction characteristic of micro-irrigation is not taken into account. (3) The crop water stress is not taken into account in this calculation, as the time step of 1 month causes this to be quite complex to represent. However, significant stress may impact the evolution of the NDVI signal. The approach of estimating K_c,sat is displayed in the light gray area in Fig. 2.

Figure 2 Processing steps flowchart. The light gray area of the flowchart was presented in Le Page et al. (2012). The calculations on K_c are performed on 29 irrigated areas. The input data are indicated in gray boxes. Variable extraction is indicated by a square box. Processing steps (and their results) are indicated by boxes with beveled corners. The indications between parentheses show the temporal, spatial, and scenario resolutions. The scenario generation described in this article is indicated in dark gray.

Figure 3 (a) Quantities used to perform the simulation. The example illustration on the left depicts the synthetic K_c of the N'Fis private area where a mix of tree crops (olive, orange, apricot) and cereal crops (mainly wheat) are cultivated. The peak occurs during March and April when cereal crops are at their maximum development; the valleys correspond to the summer months when there are no seasonal crops. The three curves are the K_c time series obtained from remote sensing processing (solid line), the linear estimation (K_c,lin, Eq. 3, dotted blue line), and the estimation corrected with rainfall (K_c,sim, Eq. 6, dashed red line). K_c,Trees has been set to 0.55 in Eq. (8), and minKc is the yearly minimum of K_c on K_c,sim. (b) The fitted correction curve of K_c according to the cumulated rainfalls (P_corr) is almost linear in this example.

3.2 Trend projection of K_c

The simulation of K_c was carried out as follows: a linear adjustment from the original K_c,sat is obtained from remote sensing; secondly, this first guess is corrected according to the rainfall of the year; finally, a correction is made according to the percentage of tree crops. Figure 3 illustrates the quantities used to perform this calculation. The four steps detailed below were reproduced for each of the 29 irrigation areas of the study.

1. The first step was to perform a linear adjustment (Eq. 3). For each month M, of the year y, a linear least-squares fit was made to the K_c,sat values for all years in the series; thus, we obtained a set of 12 monthly regression equations for K_c,sat as a function of the year. These 12 linear regressions thus formed the time series $K_{{\mathrm{c}}_{\mathrm{lin}}}$. An example is plotted in blue in Fig. 3:

   $$\begin{array}{}\text{(3)}& \begin{array}{rl}{K}_{\mathrm{c},{\mathrm{lin}}_M}& \left(y\right)=a{\mathrm{1}}_{M}y+b{\mathrm{1}}_M(\mathrm{where}a\mathrm{1}\mathrm{and}b\mathrm{1}\\ & \mathrm{are}\mathrm{fitted}\mathrm{to}{K}_{\mathrm{c},{\mathrm{sat}}_M}\left(y\right)\mathrm{for}M\in [\mathrm{1},\mathrm{12}]).\end{array}\end{array}$$

2. The second step consisted of correcting the first approximation of K_c,lin for the impact of rainfall. In this semiarid region, sowing and the development of the vegetation depend very much on the time distribution and accumulation of rainfall during the agricultural season. We chose March as the month most representative of these interannual differences (Le Page and Zribi, 2019). On the one hand, we calculated the difference in K_c (ΔK_c) between the observation (K_c,sat) and the linear regression (K_c,lin) in the month of March as in Eq. (4). On the other hand, the rainfall accumulation P_y was summed between September and March according to Eq. (5). This has been done for each year from 2000 to 2016. We found that the relationship between ΔK_c and P_y was best modeled using a second-degree polynomial according to Eq. (6). An example of this fit is given on the right side of Fig. 3.

   $$\begin{array}{}\text{(4)}& {\displaystyle }\mathrm{\Delta }{K}_{\mathrm{c}}\left(y\right)=K_{\mathrm{c},\mathrm{lin}}\left(y\right)-K_{\mathrm{c},\mathrm{sat}}\left(y\right),\text{(5)}& {\displaystyle }{P}_y=\sum _{m=\mathrm{sep}}^{\mathrm{mar}}{P}_m,\text{(6)}& {\displaystyle }\begin{array}{rl}{K}_{\mathrm{c},\mathrm{cor}}& \left(P_y\right)=a\mathrm{2}\cdot P_y^{\mathrm{2}}+b\mathrm{2}.P_y+c\mathrm{2}(\mathrm{where}a\mathrm{2},b\mathrm{2},\\ & \mathrm{and}c\mathrm{2}\text{ are fitted to}\mathrm{\Delta }{K}_{{\mathrm{c}}_M}\mathrm{for}{\mathrm{P}}_y).\end{array}\end{array}$$

3. In the third stage, the tree crops and their preservation were taken into account. In fact, in semiarid Mediterranean regions, evergreen tree orchards (olive, citrus) are the major crops responsible for green vegetation in the dry season of summer, although the cultivation of summer vegetables is also possible. In the Tensift region, orchards come first in order of preference for water allocation in agriculture, especially in dry years. The aim is to preserve the tree crops from drought. To ensure that the tree crop area would not be reduced from one year to the next, the minimum K_c,sim was forced to not run in the opposite direction to the trend slope (it will not be reduced when the trend is growing and inversely). In other words, the yearly rate of minimum K_c is growing; a new value is ensured to be greater or equal to the minimum value of the previous year. Finally, to ensure that the resulting synthetic K_c would be in line with the proportion of trees in the irrigated area, the percentage of tree crop has been obtained during the dry season when it is assumed that there is no other crop than trees:

   $$\begin{array}{}\text{(7)}& \mathit{\%}\mathrm{Trees}\left(y\right)=\mathrm{minKc}\left(y\right)/K_{\mathrm{c},\mathrm{Trees}},\end{array}$$

   where minKc is the minimum K_c,sim for the hydrological season (see the blue line in Fig. 3), and K_c,Trees is the maximum K_c for the tree crop considered (see the orange line in Fig. 3). For the study area, in which olive trees are the dominant tree crop, K_c,Trees has been set to 0.55 (Allen et al., 1998). The potential maximum K_c of each year (maxKc) was computed with the weighted sum as follows:

   $$\begin{array}{}\text{(8)}& \begin{array}{rl}\mathrm{maxKc}\left(y\right)& =K_{\mathrm{c},\max }\cdot (\mathrm{1}-\mathit{\%}\mathrm{Trees}(y\left)\right)+K_{\mathrm{c},\mathrm{Trees}}\\ & \cdot \mathit{\%}\mathrm{Trees}\left(y\right),\end{array}\end{array}$$

   where the non-tree-crop area has a maximum K_c of K_c,max, which was set according to winter wheat as in Allen et al. (1998). Winter wheat is a dominant crop of the study area, and its maximum K_c (1.15) is among the highest known.

4. The final estimate of monthly K_c (K_c,sim) integrated the linear trend, the yearly variation due to rainfall, and the limitation due to the coverage of tree crops (Eq. 9).

   $$\begin{array}{}\text{(9)}& \begin{array}{rl}{K}_{\mathrm{c},{\mathrm{sim}}_M}\left(yP\right)=& \max [min\left(K_{\mathrm{c},{\mathrm{lin}}_M}\left(y\right)+K_{\mathrm{c},\mathrm{cor}}\left(P\right),\right.\\ & \left.\mathrm{maxKc}\left(y\right)\right)K_{\mathrm{c},\min }],\end{array}\end{array}$$

   where the variability of K_c due to yearly rainfall (K_c,cor) was added to the first guess K_c,lin. It was ensured that this sum did not exceed the potential maximum K_c (maxKc, Eq. 8). It was also ensured that it did not go below the K_c of bare soil ($K_{\mathrm{c},\min }=\mathrm{0}$).

The set of Eqs. (3) to (9) together model the future trend of K_c. The water demand was then computed from Eq. (1), where ET₀ incorporated the weather conditions, and in particular the modifications of temperature predicted by the RCPs.

The performance of the model has been assessed using different sampling techniques to obtain calibration and validation datasets. First, calibration years are selected to be equally spaced along the time axis: 1 out of every 2 years (8 years), 1 out of every 3 years (6 years), and 1 out of every 4 years (4 years) of calibration. For the “every 2 years” sampling, there were two sets of calibration years; for the “every 3 years” sampling, there were three calibrations sets, and so on. The three different versions of the model were run for the different combinations of calibration and validation datasets. The versions tested were the linear fit (K_c,lin), the linear fit with rain correction (K_c,cor), and then with tree correction (K_c,sim). The performance indicators were averaged for each group ($\mathrm{1}/\mathrm{2},\mathrm{1}/\mathrm{3},\mathrm{1}/\mathrm{4}$). The years have also been separated into two groups: years with less precipitation are called the “dry years” (2001, 2002, 2005, 2007, 2008, 2011, 2012, 2014, 2016) and the rest are called “wet years”. Three statistical metrics have been computed: the determination coefficient r² (Eq. 10), the root mean squared error (RMSE, Eq. 11), and the standard deviation of error SE (Eq. 12):

$$\begin{array}{}\text{(10)}& {\displaystyle }{r}^{\mathrm{2}}=\sqrt{{\displaystyle \frac{\mathrm{Cov}(X,Y)}{\mathit{\sigma }X\mathit{\sigma }Y}}},\text{(11)}& {\displaystyle }\mathrm{RMSE}=\sqrt{{\displaystyle \frac{{\sum }_{\mathrm{1}}^n(X-Y{)}^{\mathrm{2}}}{n}}},\text{(12)}& {\displaystyle }\mathrm{SE}=\sqrt{{\displaystyle \frac{{\sum }_{\mathrm{1}}^n(X-Y)}{n}}}.\end{array}$$

3.3 Alternative evolution of K_c

The set of Eqs. (3) to (9) was used as a basis for simulating the long-term trend of K_c. The coefficients a1 and b1 of Eq. (3) were fitted with the data from the entire time series. The parameters a2, b2, and c2 of Eq. (6) were fitted using the rainfall data taken from the downscaled climate scenario.

The system of equations allows for the consideration of various possibilities for “bending the curve” (Raskin et al., 1998), which here means modifying the trend of the cropping scenario: (1) the parameters a1 and b1 in Eq. (3) can be modified to account for an overall change in the trend; (2) the coefficients a2, b2, and c2 in Eq. (6) can be modified to account for the influence of cumulative rainfall. (3) The hypothetical law that determines that the rate of change in tree coverage always has the same sign could be relaxed so that a decrease of tree coverage could be possible.

As the status of the Haouz aquifer is becoming critical, alternative trend scenarios have been proposed by the Tensift Hydraulic Basin Agency (Agence du Bassin Hydraulique du Tensift) with the support of GIZ (Gesellschaft für Internationale Zusammenarbeit) (ABHT GROUP AG – RESING, 2017). The main drivers defining these changes include (1) an expansion of irrigated areas through groundwater pumping; (2) the conversion of surface irrigation to drip irrigation; (3) an intensification of arboriculture; and (4) an abandonment of irrigated cereals for vegetable crops. However, neither the plans nor the interviews specify the intensity or locations of these changes. Here, we do not address all of the measures considered, only those concerning land cover changes. A reduction in the rate of increase of olive tree areas by 50 % is expected to begin in 2021. At the same time an intensification due to intercropping and the replacement of traditional Moroccan Picholine (100 to 200 trees per hectare) by Spanish Arbequina with a density of 1000 trees per hectare will take place. This is correlated with the expected production increase of 27 %. The fate of citrus trees is less clear. The recent trend is a strong expansion; however, the plan expected a stabilization of citrus orchards due to the lack of water especially in the western part of the basin. Finally, a progressive reduction of cereal crops (−722 ha/yr) is expected. Here, we proposed to only modify the main trend of the K_c curve in order to represent the idea of a pause in crop-cover expansion. The set of equations to simulate the bending of the curve are the following:

$$\begin{array}{}\text{(13)}& \begin{array}{rl}& a{\mathrm{1}}^{\prime }=a\mathrm{1}\cdot b_{\mathrm{c}},\\ & b{\mathrm{1}}^{\prime }=\left(a\mathrm{1}\cdot b_{\mathrm{c}}+b\mathrm{1}\right)-y_{\mathrm{c}}a{\mathrm{1}}^{\prime },\\ & K_{{\mathrm{c}}^{\prime },\mathrm{sim}}=a{\mathrm{1}}^{\prime }\cdot K_{\mathrm{c},\mathrm{sim}}+b{\mathrm{1}}^{\prime }.\end{array}\end{array}$$

A new slope (a1^′) and intercept (b1^′) were computed according to a bending coefficient (b_c) and applied for the desired years where y_c is the beginning year of the bending. The bending coefficient was expressed as a percentage.

In order to demonstrate the feasibility of applying scenarios with this model, it was determined that the introduction of two bending points to the trend scenario (Eq. 13) gave a good representation of the vision of the narrative scenario (−50 % on the trend in 2020 and another −50 % by 2040). Note that the narrative scenario did not give a detailed view of the changes, so those coefficients have been applied uniformly to the 29 irrigated areas.

Figure 4 Performance indicators (r², RMSE, SE) for the three variations of the model (K_c,lin, K_c,cor, K_c,sim) as a function of the splitting strategy for determining the calibration/validation dataset.

4 Results

4.1 Performance of the simulated K_c,sim

Figure 4 shows the three statistical metrics of the fit for different calibration groups. The calibration is computed against a subset of the 16 years, where one-half ($\mathrm{1}/\mathrm{2}$) means 8 calibration years for 8 validation years, one-third ($\mathrm{1}/\mathrm{3}$) means 5 calibration years for 11 validation years, and one-quarter ($\mathrm{1}/\mathrm{4}$) means 4 calibration years for 12 validation years. r² values were generally located around 0.5, which indicates an average fit, but they decreased when the calibration set size decreased. RMSE ranged between 0.1 for half of the calibration years to 0.17 for the corrected methods for only four calibration years. It must be noted that the K_c,sim and K_c,cor versions lower the performances regarding RMSE when only $\mathrm{1}/\mathrm{3}$ or $\mathrm{1}/\mathrm{4}$ of the years is used. This is also true for the dry calibration years, which is explained by the fact that the rainfall correction (K_c,cor) needs a good representativeness of extreme cumulative precipitations.

Figure 5 Examples of K_c simulations compared to the original time series for four selected areas.

The trajectory of each irrigated sector can be very different, according to its mix of crop, cropping intensity, and evolution. Figure 5 shows an example of this diversity with four contrasted irrigated areas. The figure also shows the calibrated and original time series of K_c when the calibration is carried out over the entire time series. The Chichaoua private area gave the worst result. As can be seen, the low r² was mainly due to the poor distribution of the data at very low K_c. It is very likely that most of the area was not irrigated continuously during the whole study period. The N'Fis right bank is located on the outskirts of the city of Marrakech, under high urban pressure. The negative trend was well predicted and the correlation score was satisfactory (r²=0.77) despite the probably erroneous pause in the original K_c data in 2006 and 2007. The Ourika traditional area is located at the outlet of the Ourika River. The scores were average (r²=0.646, slope =0.747). The general growing trend was reproduced, especially for the base of the curve. The maximum and minimum K_c values simulated with the amount of annual rainfall were generally accurate, except for the 2 years (2002 and 2003). Those are extremely dry years where rainfed cereals could not develop, and there were strong irrigation restrictions. The Bouidda district is an area irrigated by the Moulay Youssef reservoir in the eastern part of the study area. It obtained one of the best scores (r²= 0.863). The area is dominated by cereals, which explains the typical seasonal K_c curve. The increase in permanent tree crops was well simulated at the baseline. There was also an intensification of cereal crops that increased the peaks at the end of the study period. The lower K_c years, 2006 and 2007, were well simulated. The very dry year 2001 was simulated but did not reach the minima observed in the original curve.

The performance of the algorithm may be summarized as follows. The simple linear relationship provided an average approximation. The rainfall correction improved the average performance if the calibration years were representative of the validation years; elsewhere, it worsened the prediction. There was no visible effect of the tree correction during the period of time when observations were available. The bounding constraint was basically not reached. Finally, the model results were much better when the calibration years were representative of the rainfall variability (wet and dry years). Given that the remote sensing dataset was processed for the 2000–2016 period (16 years), the performance of the prediction at the 2050 horizon should be comparable to the $\mathrm{1}/\mathrm{3}$ test: r²= 0.5, RMSE = [0.1–0.14], and SE = [0.02–0.03].

Figure 6 Examples of long-term simulation of K_c in four different irrigated areas with the trend and alternative scenarios of K_c and the RCP4.5 and RCP8.5 climatic scenarios.

4.2 Evaluation of the trend and alternative evolution under climate change

Figure 6 shows examples of the extrapolation of K_c for four contrasted irrigated areas, with the trend and alternative K_c scenarios and RCP4.5 and RCP8.5 climatic scenarios. The region of private irrigation in Mejjat was all but covered with drylands until the early 90s. A boom of groundwater exploitation has provoked the development of irrigated areas. K_c spikes in the rare wet years, but most of the time, the amplitude of the signal is very small, showing that the area is mostly composed of bare soils and tree crops. The trend projection predicts that the K_c will tend toward an average of 0.2 by the year 2050, which means the impact of irrigation will be much greater than at the beginning of the 21st century. The alternative K_c greatly reduced the trend and stabilized around 0.15. Ourika is an irrigated area located at the piedmont of the Atlas mountains. This area has traditionally relied on spate irrigation from the Atlas wadis fed by snowmelt. This area is also located over the Haouz aquifer where it benefits from recharge from below-river flows. The full coverage with tree crops was attained around 2040 in the Ourika area. Therefore, the synthetic K_c of this area is controlled by the maximum K_c value of Eq. (8), which was set equal to the K_c of olive trees (0.55), which is the dominant tree crop of this region.

On the alternative trajectory, K_c is less saturated, meaning that there is still a place for annual crops. As stated earlier, the N'Fis right bank irrigated area, which is located near Marrakech, already had a trend toward the disappearance of irrigated crops. The two projections maintained the K_c level near zero except for wet years. The Bouidda irrigated area in the Tessaout (also spelled Tassaout) region follows a trend of reconversion toward tree crops without a trend of cropped area development, which is almost fulfilled by 2050. In the alternative scenario, this trend was controlled, and the area maintained a mix of annual and tree crops.

In those four examples, the climatic scenarios of rainfall only had a minimum impact on the K_c time series.

Figure 7 Trend and alternative scenarios for irrigation water demand in the four planning areas of the Tensift with RCP8.5.

Figure 7 shows the relative evolution of water demand compared to the year 2000 and the expected impact of the “bending” effect for the four different planning areas in the watershed (see Fig. 1). The expected temperature increase between RCP4.5 and RCP8.5 is about 1 ^∘C by 2050, leading to a 2 % difference in yearly ET₀. The impact on irrigation water demand is therefore minimal by the end of the simulation and has not been plotted in Fig. 7. However, we must keep in mind that if climate change only slightly impacts the water demand, it will probably reduce the water supply, because the runoff decrease will reduce the surface water available for irrigation (snowmelt, runoff, reservoirs). The area equipped for irrigation (AEI) also remains unchanged while the actual irrigated area may change inside those AEIs. Those simulations showed that each area had it's own dynamics (bearing in mind that subareas also had their own dynamics). The Mejjat area, which was mostly unexploited in the year 2000, has been experiencing significant growth since. We can also see that by 2040 the demands tend to stabilize in the trend lines owing to the saturation effect of the irrigated areas. In the alternative scenario, the reduction in the agricultural development rate stabilizes at a lower level. The N'Fis area showed the lowest increase, and the modeled trend toward a flat line showed that the area will be fully dominated by tree crops by 2050. The Haouz and Tessaout areas showed the most significant increases in irrigation water demand, as it almost tripled during the 50 years of simulation. However, the alternative scenario mitigated these increases.

5 Discussion

The proposed method allowed for projecting irrigation water demand by including both the anthropogenic and climatic vector of changes. The output is given for the next 30 years for each demand site at the monthly time step, which is a timescale that is commonly used by water planners. Various improvements to the method have been identified. (1) K_c was determined by a precalibrated linear relation to the NDVI time series. If an actual evapotranspiration product (see, for example, Xu et al., 2019) would be available for the whole period of study and at a resolution compatible with the size of the irrigated areas, K_c could be determined by computing the ratio between actual evapotranspiration products and reference evapotranspiration. (2) The gridded data were synthesized at the monthly time step for each irrigated area, by computing arithmetic averages, which resulted in a loss of information. If needed, other statistical indicators, such as the dominant land cover of the irrigated area, could be kept for further analysis. (3) The main temporal trend was fitted with monthly linear regressions (Eq. 3). The degree of the regression could be easily adapted if the model does not fit adequately with the observations.

The very low difference in impact between the two RCP scenarios is interesting. It can be explained by the short time range (2050) which resulted in a difference of only 1 ^∘C and almost no difference in precipitation over the study region. It should be noted that the predicted impact of climate change on precipitation is more variable in space and less certain than that for temperature. It also may be due to the absence of the impact of including increasing atmospheric CO₂ concentration prescribed by the RCP scenario in our approach. Different strategies are possible. For example, Fares et al. (2015) proposed a modification of the ratio bulk resistance over stomatal resistance $r_{\mathrm{s}}/r_{\mathrm{a}}$ in the ET₀ equation based upon the CO₂ emission scenario. The low impact of RCP scenarios could also be explained by the fact that the potential decrease in the annual crop growth cycle duration associated with temperature rise is not represented. According to Bouras et al. (2019), the reduction of the cycle can reach 30 % at mid-century under the RCP8.5 scenario for wheat in the study region. Similar results were obtained with the ISBA-A-gs model (Calvet et al., 1998) when studying the impact of climate change on Mediterranean crops (Garrigues et al., 2015). It is also interesting to note that if the water demand for seasonal crops is reduced with the crop cycle duration, the yield may also be affected by extreme temperatures (Hatfield and Prueger, 2015) so that the trend might be affected in order to reach the production objective.

Figure 8 IWR_G according to different efficiencies for 2 years of example (September 2011 to August 2013) in the case of an irrigated wheat. In the best case without drip curve, all coefficients are set to one except the efficient rainfall, which is set to a classical 0.75 (σ=1, φ= 1, α= 0.75, β= 1). In the best case with drip, The σ coefficient of localized irrigation is set to 0.7 (σ= 1, φ= 0.7, α= 0.75, β= 1). The best case with drip but without taking into account rainfall simulates the irrigation practice done to prevent soil leaching (σ= 1, φ= 0.7, α= 0, β= 1). The dominant surface irrigation simulates an irrigated sector where surface (gravity) irrigation dominates (α= 0.5) and where the system performance is good (σ= 0.9) (σ= 0.9, φ= 1, α= 0.9, β= 0.5). The worst case depicts an irrigated sector with low irrigation efficiencies (σ= 0.75, φ= 1, α= 0.75, β= 0.5).

Finally, the assessment of irrigation water demand should also take into account different losses at the system level (storage, transportation, and operating losses) and at the plot level (deep percolation and runoff). The future evolution of the irrigation framework toward more efficient systems such as drip irrigation could also be considered using Eq. (14). This equation of gross irrigation water requirement (IWR_G expressed in m³) is rewritten from Eq. (1):

$$\begin{array}{}\text{(14)}& {\mathrm{IWR}}_{\mathrm{G}}=\mathrm{1}/\mathit{\sigma }{\displaystyle \frac{\mathrm{AEI}}{\mathrm{10}}}{\displaystyle \frac{\mathit{\phi }{\mathrm{ET}}_{\mathrm{0}}{K}_{\mathrm{c}}-\mathit{\alpha }P}{\mathit{\beta }}},\end{array}$$

where AEI (hectares) is the area equipped for irrigation. The system efficiency σ is the ratio of water delivered upstream to the water that is distributed to the plots (Blinda, 2012). The coefficient α produces the efficient rainfall similar to that proposed by USDA-SCS (1970). The coefficient φ is used to account, in a simple manner, for a reduction of evapotranspiration due to micro-irrigation. β is the effective irrigation coefficient; it is low for surface irrigation (50 %–70 %) and close to 100 % for drip irrigation. The behavior of Eq. (14) is displayed in Fig. 8 for 2 sample years (September 2011 to August 2013) for the R3 irrigated sector where wheat is the dominant crop. According to those different efficiency configurations, IWR can vary from single to more than quadruple between the ideal best and the worst situation. However, note that those configurations are theoretical. The best case with drip irrigation is virtually unachievable in real conditions.

6 Conclusions

Owing to their inherent complexity, scenarios of agricultural evolution based on interviews are very difficult to translate into numbers (when and how much) and to represent spatially (where). The simple and flexible statistical approach proposed here is a possible solution for quantifying spatially distributed scenarios of agriculture evolution in the context of climate change and irrigated areas that are rapidly changing owing to socioeconomic influences. This is the case for many Mediterranean areas like the Tensift watershed, which was the case used for the current study. The performance of the model was acceptable in most irrigated areas, giving monthly correlation coefficients for K_c of up to 0.92 but with significant differences between the irrigated areas. It was shown that the prediction performance of the model was a function of the length of the calibration period. With 16 years of data (2000–2016), the prediction could only be done for two other periods of 16 years (until 2050). In order to demonstrate the flexibility of the model for scenario building, a local scenario of water resource management was reinterpreted to build an alternative scenario upon the trend scenario. Two climate scenarios (RCP4.5 and RCP8.5) were introduced into the trend as alternative scenarios. The local anthropogenic scenarios of irrigation water demand had rapidly changing dynamics and were spatially contrasted. On the contrary, the difference in impact between the two climate change scenarios appeared to be very small over the next 30 years, although the impacts of atmospheric carbon concentration on stomatal controls of the plant were neglected. Finally, the discussion showed that the approach could easily be combined with irrigation efficiency scenarios. The conjunction of the present approach with irrigation efficiencies would probably be a suitable combination for the modeling of water resources management. The approach can be applied to other regions where there is a significant growth of irrigated areas. The MODIS products are available over all the globe, so the retrieval of a 20-year-long time series of crop coefficients can be done anywhere. Different alternative techniques exist to retrieve K_c so that it is possible to switch between methods. A caveat of using these data, however, is their relatively low spatial resolution (250 m), which implies that some field-scale details might be missed or undersampled. The separation of the territory between different irrigated areas seems adequate in particular when using a water allocation model. The synthesizing algorithm is efficient and can be applied in different regions. The model itself is mainly linear, even if it is corrected by yearly precipitation and an eventual saturation due to the extension of tree crops. In other situations, for example, a late expansion of irrigated areas or where a reduction of the expansion is already noticeable, it will be necessary to switch to a nonlinear system. A simple ordinary least square adjustment has been used to fit the coefficients, but other more sophisticated techniques could be used. Those different questions are being taken into account in the frame of the AMETHYST project over the Merguellil watershed in Tunisia.