A Flood of Droughts: New Perspectives on Droughts in the Netherlands

(1)

A Flood of Droughts: New Perspectives on

Droughts in the Netherlands

Jos Peters

s2701634

(2)

Master’s Thesis EORAS & Economics Supervisor: Drs. Vullings

(3)

Abstract

(4)

1 Introduction

2018 was a disastrous year for Dutch farmers, as extreme droughts caused many crop failures (RTL Nieuws, 2019). They were not the only victims of the extreme drought that year. A broad range of sectors, such as agriculture, shipping, en-ergy, and industry faced consequences of the lack of precipitation (Ecorys, 2019). Therefore, the Dutch government gave order to an overall economic assessment of the 2018 drought. The damage is estimated to be between 450 and 2080 mil-lion euros, largely due to the direct effect on crops and grassland (Ecorys, 2019). For precise estimates of the financial consequences, the definitions of droughts are important. Since droughts depend on multiple hydrological variables, such as evaporation and temperature, there is no uniform definition. The various definitions could yield different perspectives on the severity of droughts, and thus on the economic damage.

We will introduce several drought indices, which we apply on the Dutch drought history. In addition, we use extreme value theorems to model the extremely dry events. Next, we will estimate the economic impact based on each index. The financial consequences of droughts are estimated using two approaches. First, we relate the historical damage in characteristic years to the probability of ex-ceedance. These are years that represent a certain level of dryness. Second, we directly relate the different indices to crop yield losses.

The dynamics of droughts are not as clear as those of other meteorological concepts. They do not just depend on a lack of precipitation but the joint be-haviour of hydrological, agricultural and meteorological factors all contribute to the formation of drought (Chang et al., 2016). Another issue in quantifying droughts is its spatial and temporal diversity (Kao and Govindaraju, 2010). In order to identify and assess events of drought, much research has been devoted to the construction of indices. So far, the consensus is that none of them led to a universal broadly applied index. Incorporating more variables in an index generates a more diverse view, however, it could be less effective due to its in-creasing complexity. Therefore, each indicator gives its own view on drought and emphasizes certain factors and variables. The first research question of this study is: which indices describe Dutch droughts well?

(6)

can be measured by means of various variables, such as stream flow, precipi-tation or ground water. Therefore, they constructed an indicator which uses only one of those inputs to monitor drought. The Standardized Precipitation Index (SPI) is defined as the difference of precipitation from the mean, divided by its standard deviation (McKee et al., 1993). So far, this is the most widely used index as it can easily be applied to all different inputs related to droughts. Moreover, it is spatially consistent and easy to interpret in terms of probability, which makes it suitable for risk assessments (Guttman, 1998). The World Me-teorological Organization (WMO) recommended the SPI to be used worldwide by the national weather services (Standardized Precipitation Index User Guide, 2012). The disadvantage of this index is that it models droughts as univari-ate phenomena. In reality, droughts are complex multivariunivari-ate phenomena that require the incorporation of joint behaviour. Finally, the index uses a history of a predetermined amount of months to determine the standard deviation and mean, which does not correct for seasonality.

There are several indices that allow for the multivariate nature of droughts, such as the Palmer Drought Index (PDI), introduced by Palmer (1965). A long series of hydrologic data in a specific area is used to calculate the amount of moisture that is needed for normal conditions (Guttman, 1998). The monthly deviations from this normal situation are integrated in a single drought index. The normal conditions are based on a combination of loss of soil moisture, runoff, recharge and evapotranspiration. Therefore, it meets the call for a multivariate approach, in contrast to the SPI. The PDI of the current month is largely based on the PDI of the previous month, and therefore it exhibits a long-term memory. However, the complexity and long-term memory also make the index harder to calculate, interpret and compare (Guttman, 1998).

Over the years, both indices are widely studied and further extended. The Standardized Precipitation Evapotranspiration Index (SPEI) adds evaporation to the SPI based on temperature (Vicente-Serrano et al., 2009). This allows for the joint effect of evaporation and precipitation from the PDI, while keep-ing the comparability advantages of the SPI. Another index that makes use of evaporation is the Reconnaissance Drought Index (RDI), proposed by Tsakiris et al. (2005). The Standardized Soil Moisture Index (SSI), developed by Xu et al. (2018) is a simple index to measure soil moisture. This index is combined with the SPI to establish a Multivariate Standardized Drought Index (Hao and Aghakouchak, 2013). Examples of extensions on the PDI include the Palmer Hydrological Drought Index (PHDI) and Z index (Palmer, 1965).

(7)

(8)

Examples include van Beek et al. (2008), who calculated the financial efficiency of anti-drought measurements, and Haskoning (2007) who investigated the eco-nomic consequences of various climate scenarios. Since the input data for the model is not publicly available for all years, we use the estimated damage in characteristic years found in previous studies.

Next, we examine the relationship between the indices and crop damage. The question we would like to answer is: to what extend can the indices be used to estimate crop damage? In order to answer that question, we relate crop yield losses to the different indices using OLS regressions. This has for example been done by Garc´ıa-Le´on et al. (2019), who used OLS regressions to predict yields of Spanish cereals based on drought indices. Meyer et al. (1993) developed a crop specific drought index to monitor the relationship between droughts and crop yields. Mavromatis (2006) compared various indices to estimate wheat yields in Greece, and Quiring and Papakryiakou (2003) analysed indices used in the Canadian agricultural sector.

The physical crop damage is not enough to estimate the financial impact for farmers. Reinhard et al. (2015) argue that price and market effects have to be measured as well. They make a distinction between short and long term eco-nomic effects of drought. In the short run, macro-ecoeco-nomic changes will push prices upwards as a result of lower outputs. The magnitude of the effects de-pends on the price elasticity of demand of the corresponding crops. Therefore, we include a description of the market effects corresponding to the supply shock after droughts.

Climate change could also be a factor in the estimation of expected damage. So far, the KNMI index is based on the assumption that droughts are indepen-dent and iindepen-dentically distributed. The latest KNMI report states that there is no trend in droughts over the years (Sluijter et al., 2018). They described that both precipitation and evaporation are increasing. Since these two trends cancel out, there is no trend in the amount of droughts. However, extreme droughts in recent years led to a re-evaluation of the impact of climate change.

(9)

(10)

2 Data

The KNMI has an extensive weather monitoring network consisting of auto-matic weather stations, satellites, radars, and other measuring instruments (Waarneemnetwerk: KNMI Kennis-& datacentrum, 2020). We use data from 48 automatic weather stations, and 322 manual rain meters.

2.1 Precipitation

Both the automatic weather stations and the rain meters measure precipitation. We decided to use data from the rain meters, since they meet the international conditions set by the World Meteorological Organization (Vrijwillige neerslag-meters: KNMI Kennis- & datacentrum, 2020). They make use of the Hellmann meter, which is one of the most used precipitation instruments for official sta-tions worldwide (Sluijter et al., 2018).

Time series of all meters are freely available on the website of the KNMI (KNMI Klimatologie, 2020). The precipitation data are measured per 24-hour period, and are validated every ten days. We used data from 11 locations, that provided long and homogeneous time series between 1906-2020. The selected places are displayed in Table 1, together with the station number and the local average daily precipitation. As shown in Figure 1, the locations are evenly spread out across the country. For the meter in Heerde, data is missing during the 15th week of 1952. This data is interpolated using the mean of the other 10 stations during this week. Figure 2 shows the distribution of the average monthly precipitation of the 11 rain meters.

Location Station Number Average daily precipitation (mm)

De Bilt 550 2.239 Groningen 139 2.195 Heerde 328 2.201 Hoofddorp 438 2.251 Hoorn 222 2.138 Kerkwerve 737 2.024 Oudenbosch 828 2.132 Roermond 961 1.881 Ter Apel 144 2.045 West Terschelling 011 2.118 Winterswijk 666 2.136

(11)

Figure 1: Location of the 11 weather stations.

Figure 2: Distribution of the average monthly precipitation in the Netherlands of the 11 meters in Table 1.

2.2 Evaporation

(12)

Makkink to estimate it. This formula uses global radiation and temperature, to estimate the ’reference evaporation’. This is the amount of evaporation in a shortly cut, well saturated grass field (Sluijter et al., 2018).

The temperature data is available since 1906, however, global radiation is only measured in the Netherlands since July 1957. Therefore, we need to manually construct a variable representing radiation levels for the period before 1957. Be-low we discuss the applied methodology, folBe-lowing methods similar to Beersma and Buishand (2002).

Global radiation will be estimated based on the amount of sunshine combined with the time of the year, using the method of Frantzen and Raaff (1982). They provided formulas to estimate the radiation at specific locations. For three of those locations, all required data are available. These locations are represen-tative for specific parts of the Netherlands. The location in Eelde is used to estimate radiation for all stations in the northern part of the Netherlands. De Bilt is used for the centre, and Maastricht represents the stations in the south. Hence, the estimated radiation in each of the 11 stations is based on one of the three formulas. This is an improvement compared to the method of Beersma and Buishand (2002), who used the location of De Bilt to estimate the radiation for all stations.

The formulas used to estimate radiation in the three locations are displayed below.

Eelde: G =

45.0 + 32.9 ∗ sin(10t − 84.79) for S = 0 70.7 + 52.1 ∗ sin(10t − 85.57) + 137.5 + 110.1 ∗ sin(10t − 83.65) for S > 0 De Bilt:

G =

43.0 + 33.3 ∗ sin(10t − 86.04) for S = 0 65.1 + 47.0 ∗ sin(10t − 86.47) + 142.0 + 107.6 ∗ sin(10t − 82.29) for S > 0 Maastricht:

G =

44.3 + 31.8 ∗ sin(10t − 88.19) for S = 0 69.9 + 49.3 ∗ sin(10t − 86.59) + 143.1 + 103.3 ∗ sin(10t − 80.66) for S > 0, where G is the daily global radiation. S is the daily relative sunshine, which is the percentage of the day that the sun shined. Each year is divided in 36 parts, defined as decades t.

(13)

λ ∗ ρ ∗ ET = 0.65 ∗_s+γs ∗ K,

where λ is the enthalpy of vaporization of water, ρ the specific weight of water, ET the evaporation, γ the psychrometric constant, s the slope of the vapour pressure curve, and K the global radiation.

As a robustness check, the results are compared with the use of the formula of De Bilt only. The greatest difference is visible between the stations De Bilt and Eelde. The majority of the differences between the two relate to low evaporation values. Figure 3 shows that evaporation values between 0-0.2 occur more often in Eelde, while values between 0.2-0.4 occur more often in De Bilt. These low values do not contribute much to the precipitation deficits, since they mostly occur during cloudy days. Most of the time, these days do not lead to the construction of major precipitation deficits. Moreover, the location of Eelde is only related to two weather stations (Ter Apel and Groningen). Since the depreciation deficits are based on the average of 11 stations, and the differences mostly occur at low values, there are no major differences in the results using only De Bilt.

Figure 3: A comparison between the distributions of the daily estimated makkink evaporation in Eelde and De Bilt.

(14)

Figure 4: Comparison between the estimated and actual makkink evaporation. Frantzen and Raaff (1982) name various reasons for the deviations from the actual values. Most importantly, the estimation is based on temperature and sunshine hours, so it does not measure sunshine intensity. Therefore, the radi-ation will be underestimated if the sunshine hours are in the middle of the day, when the intensity is relatively higher. As a robustness check, we compared the values of the indices based on only the period after 1957, with the values of the whole time interval. This comparison did not yield notable differences.

2.3 Crop yields

(15)

(16)

3 Methodology

In this section, we will introduce the different indices. First, we describe the current Dutch drought index used by the KNMI. We will critically evaluate it, and propose an alternative measure based on the same inputs. Next, we introduce three new indices: SPI, JDI and MSPI. First, we describe how we calculate the SPI index in the Dutch context. Subsequently, we use the JDI and MSPI to compress multiple SPI values into a single index. This way, several time intervals can be combined. The JDI uses a copula to allow for dependence between the different intervals. The MSPI uses a principal component analysis to combine multiple time windows into a univariate index.

After introducing the indices, we will conduct an extreme value analysis. The extremes are of special importance in drought analysis, since they could lead to major damage. For example, the extreme droughts of 1976 and 2018 resulted in damages of at least 450 million euros (Projectgroep Droogtestudie Nederland 2004; Ecorys, 2019). The extremes appear to follow a different distribution than the other values. This could lead to imprecise estimates using the indices introduced above. Therefore, we use two extreme value models. The censored sample Generalized Extreme Value (GEV) distribution is fit on the most severe drought values. The Peak Over Threshold (POT) model estimates the distri-bution of the droughts above a certain level of severity.

Finally, we investigate if we can relate crop yield losses to drought indices. Therefore, the indices will be related to crop yield damage using OLS regres-sions.

3.1 KNMI index

The KNMI measures drought severity using the maximum precipitation deficit. It uses the difference between the amount of precipitation and evaporation. We will now describe how the maximum precipitation deficit is calculated. Next, we will discuss the strengths and weaknesses, and propose an alternative measure based on the precipitation deficit.

The Daily Precipitation Deficit at day t is denoted by DP Dt. It is the difference

between the amount of precipitation and evaporation during 24-hours, at day t: DP Dt= EVt− P Rt, (1)

where EVtis the total amount of evaporation during day t, and P Rtis the total

(17)

To capture the accumulation of the DP Dt values during a time interval of

multiple days, the Accumulated Precipitation Deficit (AP Dw_t) is constructed: AP Dw_t = max(AP D_tw−1+ DP Dw, 0). (2)

The time window starts at time t, and ends at time w. The AP D1wvalues are

displayed in the blue column of Table 2. The AP Dwt measures the accumulation

of the past DP Dtvalues, however, only as long as the accumulated value is

non-negative. The reason for this constraint, is the impact of drainage. If there is a abundance of precipitation, the excess water will be discharged. Therefore, if the precipitation deficit is zero, any additional precipitation will not be subtracted, until a new drought develops. For example, in Table 2 the AP D4

1 value equals

0. The additional rain on day 5 will not lead to a negative AP D5

1, instead, it

stays at value zero. By construction, this implies that long periods of moderate precipitation often have more impact on the AP Dw

t than a single day of heavy

rain. However, in reality the excess water could be stored to a certain degree. The storage capacity of the soil depends on many factors, and differs at each location. Due to the complexity of these factors and a lack of data availability, it is beyond the scope of this thesis to allow for these local differences. Also, due to extreme droughts in recent years, Dutch authorities try to make more use of the excess water by saving it. Last year, 7 million euro was invested in water storage facilities (NOS, 2020). These developments could be a reason to loosen the non-negativity constraint in further studies.

Table 2: Construction of the Accumulated Precipitation Deficit during a fictive week.

For the sake of convenience, the interval boundaries are often denoted as first day of the month. Hence, AP DJ un

M ar represents the AP D value for the time

interval between March 1st and the first of June. In other words, it captures the AP D accumulation during Spring.

The development of the AP Dw

t during the time window can be summarized

(18)

maximum value that the AP D reached between April 1st and October 1st : M AP Dtw= max(AP D t+1 t , AP D t+2 t , ..., AP Dwt) (3)

In the context of Table 2, we find:

M AP D7₁= max(4, 5, 3, 0, 0, 10, 9) = 10. (4) The advantage of the M AP DOct

Apr is that it is convenient to compare among

different years. It is always based on the AP D values during the same months, with the same seasons. The ignorance of the winter is justified by the relative absence of deficits in those months. Also, it measures the peak of the drought, which often causes the greatest damage. However, the maximum value also has limitations. By fixing the time interval, it lacks flexibility with respect to the temporal aspect of droughts. Short and medium term droughts are not incor-porated in the index, unless they contribute to the yearly peak. According to the University of Wageningen, a short term drought of 2-3 weeks could already cause damage. If this period does not result in the yearly peak, and is followed by a period of large precipitation, it is not incorporated in the M AP D_AprOct. Moreover, if the drought does not occur during the fixed period between April and October, it is not displayed in the index. Also, the index does not distin-guish between the length of the droughts. As only the peak is considered, the duration of the drought is not relevant. This could lead to underestimation of droughts with a long duration.

(19)

Figure 6: KNMI precipitation deficit construction. Considering the objections to the M AP DOct

Apr, we propose an alternative

in-dex. The Average Accumulated precipitation Deficit (AAP Dw

t) represents the

average of the AP D values during the time interval: AAP Dw_t = Pw t=iAP D i t w − t (5)

In the context of the fictive week in Table 2 we find that: AAP D₁7=4 + 5 + 3 + 0 + 0 + 10 + 9

7 ≈ 4.43. (6)

The AAP Dw

t allows us to collect more information regarding the development

of the droughts. Not only the peak, but also the duration and the intensity of the different droughts during the time window are considered. Longer droughts are weighted heavier compared to the KNMI measure. Also, the time interval is not fixed, to allow for flexible time intervals. Therefore, also droughts in fall and winter can be incorporated. Recent developments show that droughts do no longer end after September. The 2018 drought took so long, that it had impact on the start of 2019.

In order to compare AAP Dw

t values among different years, the same time

(20)

Figure 7 shows how the two different methods yield different results. The first graph shows the development of the precipitation deficit in 1958. A clear peak is visible, which leads to a relatively high M AP DOct_Apr. In contrast, the AAP D_AprOct is relatively low that year. The second graph shows the year 2000. Here we see multiple peaks. The KNMI approach only considers the highest peak, while the average accumulated value also incorporates the area under the other peaks. This is reflected in the AAP DOct

Apr, as it is higher than in 1958 although the

max-imum is lower. The last graph represents 2009. Here we see that the drought persists long at a high level. Therefore, the AAP DOct

Apr is high in relation to the

maximum.

Figure 7: Comparison between the AAP DOct

Aprand M AP D Oct

Apr, using an

illustra-tion of the precipitaillustra-tion deficits in 1958, 2000 and 2009. Droughts with longer duration and/or multiple peaks lead to relatively higher AAP DOct

(21)

3.2 SPI

The SPI index is based on a single hydrologic input variable (McKee et al., 1993). This could be any relevant measure, such as stream flow, evaporation or water storage. To construct the index, the input variable is transformed into a normally distributed index: the SPI. We will use the AAPD to construct the SPI. In this section, we describe the procedure that transforms the AAPD into the normally distributed SPI.

First, we will fit the Gamma probability distribution function to the historical AAPD observations. Based on the fitted distribution, the cumulative probabil-ity corresponding to each observation is calculated. By inserting these proba-bilities into the inverse normal distribution, we transformed them to a normally distributed variable which is called the SPI. In other words, we match the quan-tiles of the original distribution to the quanquan-tiles of the normal distribution, to construct a new normally distributed variable. This variable is directly re-lated to the original observations, but has pleasant properties that belong to the normal distribution. In this section, these steps will be described in more detail. First, we will search for a probability distribution function that provides a good fit to the AAPD data. This is necessary, because each hydrologic input variable is distributed differently. For example, the AAPD values used in this study are distributed differently than the maximum precipitation deficit of the KNMI, as displayed in Figure 8.

Figure 8: Comparison between the distribution of the KNMI measure ( M AP DOct

Apr, left), and the alternative (AAP DOctApr, right).

(22)

Figure 9: Comparison between AAP D values of different time intervals. On the left, we find the AAP DAug_{J ul}. The right displays the 12-month time interval, starting and ending in August, denoted by AAP DAug_Aug.

Although McKee et al. (1993) only fitted the Gamma distribution to the original precipitation data, further studies considered many other distribution. Accord-ing to Stagge et al. (2015) the Gamma distribution is frequently preferred due to its flexible shape parameter. However, for short accumulation periods both the Weibull and Gumbel distribution could outperform the Gamma distribution. All of the distributions above were considered, as well as the Lognormal, Gener-alized Extreme Value distribution (GEV) and three-parameter additions of the Gamma and GEV distributions, to allow for location parameters. It is impor-tant to find a suitable fit that is close to the true distribution of the hydrologic input. If the fitted distribution does not represent the underlying data well, it could lead to an under- or overestimation of the probability of droughts. For example, if we use a distribution with fat tails to model a hydrologic variable with thin tails, we overestimate the probability of extremes. A good fit at this stage, should later result in a normally distributed SPI index. Therefore, the normality of the index is used to evaluate the suitability of the different distri-bution functions.

(23)

Figure 10: QQ-plots of the Gamma, Gumbel, Lognormal and Weibull distribu-tions, fit on the AAP DJ un

M ar. They represent the AAP D values in spring.

(24)

and Wah, (2011). The test statistic W measures the degree of normality of the sample. Its distribution depends on the sample size. The exact formula can be found in the Appendix. The null hypothesis of the test states that the sample is normally distributed. The test was applied to the SPI values of all monthly time windows, with ending months March, September and December. This sums up to 36 windows in total. The most important ending month is September, as most crops are harvested in the months prior to September. We find that for approximately 80% of the time intervals in September and March, the null hypothesis is not rejected. Therefore, we find evidence in favour of the normal distribution for most time intervals. However, the time intervals with ending month December are rejected at a 10% significance level. Therefore, we need to be cautious with the interpretation of the SPI values at the end of the year. We find evidence that these values are not normally distributed, which could lead to misleading probability of exceedances of droughts. However, the other distributions, as shown in Figure 10 and 11, did not perform better for December. Therefore, we decided to use the Gamma distribution for all time windows. We are aware that the interpretation of the SPI index in the autumn could be misleading in terms of probabilities. Also, we include extreme value analysis to find probabilities corresponding to extreme observations, that are more vulnerable to deviations from the normal distribution.

After finding the appropriate parameters of the Gamma distributions for all time windows, the cumulative probability can be derived. The probability cor-responding to AAP Dw

t equals:

uw_t = F (AAP D_tw), (7) where F is the cumulative distribution function, with parameters corresponding to the AAPD data with time window starting a t, until time w.

Now we have found the cumulative probabilities corresponding to each histori-cal observation, we can insert them into the inverse normal distribution. This results in a standard normally distributed variable with variance one and mean zero. This property makes it comparable to SPI indices in other countries with varying climates.

SP I_tw= φ−1(uw_t) = φ−1(F (AAP Dw_t)). (8)

3.3 Seasonal correction

(25)

modified the SPI index. Each time interval is compared to the exact same time window in previous years, thereby comparing the same seasons. Practically, this implies that for each time interval, a separate history should be consid-ered. The cumulative probability corresponding to a certain AAPD level will then be estimated based on this unique history. For example, the probability of finding an AAPD level in May of 40mm is below 5%, since historically these values do not occur frequently. In contrast, this average is above 15% in August. The procedure is computationally intensive, as for each specific time window a corresponding history should be collected. Hence, we need to calculate the SPI values based on Gamma distributions with a wide range of parameters, based on different histories. Figure 12 below shows the transformation from the AAP D_{M ay}Apr to its corresponding modified SPI index, based on the average of the 11 weather stations. Table 3 shows the frequency of occurrence of certain SPI values.

Figure 12: Transformation from the AAP D_{M ay}Apr to the SP I_{M ay}Apr index.

Table 3: SPI and drought severity.

SPI Values Drought Category Time in Category [0 : -1) mild drought 29,82% [-1 : -1,5) moderate drought 14,03% [-1,5 : -2) severe drought 2,63% [-2,00 : -∞) extreme drought 1,75%

(26)

3.4 JDI

Most of the time, a single rainstorm does not stand on its own. Low-pressure areas could result in stormy and wet weather conditions for a prolonged pe-riod. Therefore, the amount of precipitation tomorrow, is related to the current weather conditions. To allow for this dependence, an obvious solution would be the use of auto-regressive time series methods. However, the data is clearly non-stationary, for example due to seasonality. Moreover, we would like to maintain the probabilistic character of the SPI index. Therefore, we will follow the method of Kao and Govindaraju (2010), by introducing copulas to allow for dependence between various time windows.

A copula is a multivariate cumulative distribution function, consisting of uni-form marginal distributions (Ning, 2010). In other words, copulas could be used to combine multiple univariate distribution functions into a joint distribution function. In the context of droughts we use copulas for the following application. If we find that the probability of exceeding a certain drought is 10% in January and 20% in February, we can use a copula to calculate the joint probability of the two months combined. The copula incorporates the dependence between the two periods. Hence, if January is dry, it is more likely that February is also dry, which will be reflected in the joint probability.

According to Sklar (1959), there exists a copula C : [0, 1]d _{→ [0, 1] such that}

∀(x1, ..., xd) ∈ IR,

F (x1, ..., xd) = C(F1(x1), ..., Fd(xd)), (9)

where F is a joint cumulative distribution function and F1, ..., Fd are marginal

cumulative distribution functions (McNeil et al., 2005).

We already used the marginal distribution functions in the construction of the SPI index. In (7), we used the Gamma distribution to calculate uwt. If we

insert the results of (7) in the copula formula (9), we find the joint cumulative distribution function:

F (AAP D1, ..., AAP Dd) = C(u1, ..., ud). (10)

We will now describe how we apply copulas in the context of droughts. Assume we calculated the 12 values uDec

t spanning a whole year, i.e. uDecJ an, uDecF eb. . . u Dec Dec.

Using the distribution function in (10), we can calculate the joint probability: P (U_{J an}Dec≤ uDec J an, ..., U Dec Dec ≤ u Dec Dec) = q. (11)

This represents the joint probability of exceedance of the considered year. When the beginning of the year is dry, the values uDec

J an, uDecF eb, uDecM arare relatively small.

(27)

for seasonality.

Using (11), we can find the probability q, that a certain drought is exceeded. However, there are multiple kinds of droughts which cause the same impact. For example, a short and intense drought could result in a similar value of q as a long but moderate drought. The probability of the combination uJ ul

J un< 0.4 and

uJ ul

M ay < 0.3 could be the same as the probability of u J ul

J un< 0.7 and u J ul

M ay< 0.2.

Therefore, Kao and Govindaraju (2010) decided to make the index more univer-sal. They looked for the probability of finding random variables that lead to the same level of q. In other words, which combinations of AAPD levels are equally likely. They formalized this concept, by introducing the Kendall distribution:

Kc(q) = P (C(U1, ..., Ud) ≤ q). (12)

According to Kao and Govindaraju (2010), we can also interpret this as the probability of finding a set of random variables, which satisfy the condition:

{(U1, ..., Ud) ∈ [0, 1]d|C(U1, ..., Ud) ≤ q}. (13)

Similar to the SPI, the Kendall distribution is re-scaled using the inverse normal, to obtain the universal JDI:

J DI = φ−1(Kc(q)) = φ−1(P (F (x1, ..., xd) ≤ q)). (14)

3.5 Empirical copula

There are several different types of copulas available. They all use different dependence structures between the values of uw

t. As we want to get a broad

view on drought, we consider at most 12 input variables representing a whole year. Less inputs are useful for shorter time intervals, and more will result in overlapping seasons. The latter is historically not too restrictive, as droughts often do not persist for periods longer than 12 months. In order to deal with this multivariate approach, parametric copulas are hard to fit to the data. McNeil et al. (2005) describe that there are several restrictions in constructing higher-dimensional parametric copulas. As we have sufficient observations, we can make use of a non-parametric empirical copula. The empirical copula follows from: ˆ C(ut_w1, ..., ut_wd) = 1 n n X i=1 1(U_w1ti ≤ ut w1, ..., U ti wd≤ u t wd). (15)

There does not exist a direct method to estimate the Kendall distribution cor-responding to an empirical copula. Therefore, we will make use of simulations. Using R, we will draw 10000 random realisations from the copula in (15). These represent ten thousand combinations of precipitation deficits that have the same distribution as the historical observations. Using these results, we can find how many times certain combinations occur. We use this to find the probability:

(28)

To keep the dependence structure in place, the simulation will be drawn at the copula level. As a final step, the outcome of (12) will be inserted in the inverse normal, to obtain the JDI index.

3.6 MSPI

Another index that combines multiple time intervals into a single index is the MSPI, constructed by Bazrafshan et al. (2015). To combine the time windows, they use a Principal Component Analysis (PCA). PCA is a technique that can reduce the dimensions of a dataset by analysing patterns in its variation. Eigenvectors are used to evaluate the structure of the data. Bazrafshan et al. (2015) describe that the principal components are calculated as:

P Ci= EiTX = N X n=1 eniXn, (17) where ET

i is the itheigenvector corresponding to the original data. The

eigen-vectors are sorted in descending order. The nth _{element of each eigenvector i}

is defined as eni. Xn is the nth variable of the original dataset. P Ci is the

principal component corresponding to the ith_eigenvector.

The eigenvectors are based on cross correlation between the SPI values. Figure 13 below shows that the correlation between SPI values of similar time windows is high.

(29)

components corresponding to each year based on the 12 time windows. Table 4 below shows that the first principal component indeed incorporates most of the variation in the data. The amount of additional explanation quickly reduces after the second component. As the main goal of the MSPI is to reduce the multivariate character of the data to a one-dimensional index, it is solely based on the first component. This component explains 83.4% of the variation in the original data.

Table 4: Proportions of variation explained by principal components.

Principal component 1 2 3 4 5 6 St.dev 3.164 1.112 0.670 0.409 0.330 0.140 Proportion of variance 0.834 0.103 0.037 0.014 0.009 0.002 Cumulative proportion 0.834 0.937 0.975 0.989 0.998 0.999 Principal component 7 8 9 10 11 12 St.dev 0.094 0.026 0.020 0.005 0.002 0.001 Proportion of variance 0.001 0.000 0.000 0.000 0.000 0.000 Cumulative proportion 1.000 1.000 1.000 1.000 1.000 1.000

We found that the first component captures most of the variation in SPI values between different time windows. The P C_{J an}Apr represents the first principal com-ponent for the time interval between January and April. This implies that all SPI values within this time interval are incorporated. The SP IF eb

J an, SP I M ar J an ,

SP I_{J an}Apr are compromised into the single P C_{J an}Apr. In order to compare the prin-cipal components among different time windows and locations, they need to be standardized. Bazrafshan et al. (2015) describes the following standardization:

M SP I_tw= P C w t −P C¯tw SDw t , (18)

where P Ctw is the principal component in the time window starting at t and

ending at w, and SD is the standard deviation in the same time window. P C¯w t

indicates the average principal component value during a specific time interval between t and w.

(30)

August. In contrast, the M SP I incorporate most of the variation in the interme-diate time windows as well. More precisely, the 1-month, 2-month, ..., 12-month time windows ending in August, are all represented within the M SP I_AugAug.

Figure 14: Comparison between the SP I_AugAug and the M SP I_AugAug.

3.7 Extreme value models

Extreme droughts could lead to major impact. Therefore, it is important to obtain precise estimates of the probability that these events occur. We found that the SPI index does not always satisfy the normality assumption. One of the underlying reasons is that the observations in the tail of the distribution appear to follow a different shape than the remainder of the distribution. When we fit a distribution to the entire sample, we need to make a concession with regards to the parameter values, in order to find a good fit for its entire range. The parameter estimates that yield the best fit for the tail will result in a sub-optimal fit in other parts of the distribution. Therefore, we use a separate model to explore the properties of the tail distribution. To model the tail, we can make use of two statistical theories regarding extreme values, introduced by McNeil et al. (2005). First, we model the Generalized Extreme Value (GEV) distribution for censored samples. This model is also used in the KNMI, in their 2014 drought report (Beersma et al. 2004). Thereafter, we introduce the Peaks Over Threshold (POT) model.

(31)

given by: F (x) =        exp − (1 + ξ(x − µ σ ) −1 ξ if ξ 6= 0 exp − e−(x−µ)σ if ξ = 0. (19)

The location parameters is defined by µ ∈ IR. The scale parameters is given by σ > 0, and the shape parameter is known as ξ. The latter determines the type of the distribution. McNeil et al. (2005) explain that the function reduces to the Gumbel distribution if ξ = 0. For ξ > 0, we find the Fr´echet distribution, and for ξ < 0 the type of distribution is known as Weibull.

Figure 15 shows the probability distribution function and the cumulative distri-bution function for the three different types of GEV distridistri-bution. For all three distributions the location and scale are given by µ = 0 and σ = 1.

Figure 15: PDF and CDF of the GEV distribution, for the shapes ξ < 0, ξ = 0 and ξ > 0, with location µ=0, and scale σ = 1.

We find that the Weibull distribution, with ξ < 0 has the longest tail. This means that it allows for the most extreme observations. In contrast, the Fr´echet type, with ξ > 0 has a shorter tail. In fact, the distribution is limited by an upper bound.

(32)

The fitdistcens function within the fitdistrplus package in R failed to fit the GEV distribution on a censored sample using the method of maximum likeli-hood. This problem is recognized by Wang (1990) who argued that the method of maximum likelihood often fails to fit censored sample GEV distributions, even if the observations are drawn from a true GEV distribution. Alternatively, the method of moments could be used. This method relates the population moments to the parameters that need to be estimated. The resulting expression can be solved to obtain estimates for the parameters. If the inverse form of a particular distribution is known, the Probability Weighted Moment (PWM) can be used. This is a variant of the method of moments that exploits the form of the inverse to find more convenient expressions than regular moments (Greenwood et al., 1979). As an explicit expression for the GEV distribution is available, we can make use of this method. Wang (1990) extended the PWM method to allow for censored samples. In this section, the method of Wang (1990) will be described based on his application on the GEV distribution. Probability weighted moments are defined by Greenwood et al. (1979) as:

Ml,j,k= E(XlFj(1 − F )k) =

Z 1

0

(x(F )lFj(1 − F )kdF, (20) where F = F (x) is the cumulative distribution function defined in (19). The corresponding inverse is denoted by x(F ). l, j and k are real numbers. If l is a nonnegative integer, we can recognize the conventional moment. However, sometimes there does not exist an explicit relationship between the conven-tional moment and parameters. Therefore, it could be useful to make use of probability weighted moments. For some distributions the parameters can be written as a function of the partial probability weighted moment, but not for the conventional moment. Greenwood et al. (1979) provides a table with the relationship between various distributions and an analytically convenient proba-bility weighted moments. Often, M1,j,0or M1,0,kare chosen, since in these cases

only the first power of x is taken. This results in more convenient expressions than the conventional moment.

To extend the PWM to censored samples, Wang (1990) introduced a lower bound in the integral in (20), which results in:

Ml,j,k0 = E(XlFj(1 − F )k) =

Z 1

F0

(x(F )lFj(1 − F )kdF, (21) where F0 = F (x0) is derived from the threshold x0. Wang (1990) called the

concept in (21) the partial probability weighted moment (PPWM). For compu-tational convenience, l = 1 and k = 0 were chosen, which yields:

β0_j= M_1,j,00 = Z 1

(33)

By inserting the cumulative distribution function F (x) and the inverse distri-bution function x(F ) of the GEV into equation (22), Wang (1990) found the expression for: β_j0 = (ξ +a k) 1 r + 1(1 − F r+1 0 ) − α k Γ(1 + k) (r + 1)1+kP (1 + k, −(r + 1) log(F0))), (23) where, P (1 + k, −(r + 1) log(F0)) =

Z

−(r+1) log(F0) 0 θk_e−θ_dθ Γ(1 + k)

The first three moments, β0, β1and β2, are derived by Wang (1990) using

equa-tion (23). They can be found in the Appendix. These three equaequa-tions contain the parameters ξ, α and k that need to be estimated. In order to solve for these parameters, the unknown values for β and F0 need to be found first. Wang

(1990) provided estimators for both, which can be found in the Appendix. Af-ter inserting these estimators into the three moment equations, we can solve for the three unknown parameters of the GEV distribution. In Figure 16 below, we find the QQ-plot of the censored sample GEV distribution with thresholds 150mm and 225mm.

Figure 16: QQ-plot of the censored sample GEV distribution with thresholds 150mm (left) and 225mm (right).

We find that the threshold of 150mm results in a good fit for the observations in the range 150-250mm. The observations above 250mm show a better fit using the higher threshold of 225mm.

3.7.2 Peaks Over Threshold

(34)

estimate the distribution above a threshold. The Peaks Over Threshold model can be used to model exceedances above a certain threshold. It is based on the assumptions that exceedances follow a homogeneous Poisson process in time. Moreover, the observations above the threshold are independent and identically distributed. Finally, according to McNeil et al. (2005), the distribution above the threshold follows a generalized Pareto distribution (GPD).

McNeil et al. (2005) define the GPD distribution as:

Gξ,β(x) =      1 − (1 +ξx β −1 ξ ) if ξ 6= 0 1 − e−xβ _{if ξ = 0,} (24)

where the scale parameter is defined as β > 0. The shape parameter is ξ ≥ 0 if x ≥ 0, and ξ < 0 if 0 ≤ x ≤ −β_ξ.

We use the mean excess function to test if our data indeed follows a GPD distribution above a certain threshold. The mean excess is the expected excess over threshold u, conditional on exceeding the threshold:

e(u) = E(X − u|X > u). (25)

Empirically, we can estimate this function by using the following formula de-scribed by McNeil et al. (2005):

e(u) = Pn

i=1(Xi− u)IXi>u

Pn

i=1IXi>u

. (26)

According to Kleiber and Kotz (2003), a cumulative distribution function is uniquely determined by its mean excess function. A well known example is the exponential distribution. The so called ’lack-of-memory’ property states that the distribution of a random variable is independent of its history. Consequently, the distribution does not depend on the level of the threshold, which results in a constant mean excess function. The GPD is also related to a specific mean excess form. According to McNeil (2005), the specific mean excess function corresponding to the GPD is given by:

e(u) = β + ξu

1 − ξ . (27)

(35)

Figure 17: Mean excess plot of the M AP D_AprOct.

The graph appears to be piecewise linear. The linear part between 50-150mm is less interesting, since none of the characteristic years below 150mm had any economic damage, according to Projectgroep Droogtestudie Nederland (2004). The next part between 150-225mm does not appear to be linear. Therefore, we focus on the linearity after 225mm. We fit the GPD distribution using the method of maximum likelihood, with the threshold of 225mm. The QQ-plot is shown below in Figure 18. As comparison, the QQ-plot of the censored sample GEV distribution with a similar threshold of 225mm is displayed on the right.

(36)

the model. Only the two most extreme observations appear to show a better fit using the GEV distribution model. Due to the evidence in favour of the POT model, it could be worthwhile to include this model when examining extreme droughts.

3.8 Climate change

The complexity of climate change models reaches beyond the scope of this re-search, moreover the data availability is limited. However, it is useful to examine what kind of impact climate change could have on the different indices. In par-ticular, we investigate if temperature and precipitation increases have the same effect on each index.

To reduce uncertainty and complexity, we focus on the 2030 climate scenario of the KNMI (KNMI’14-klimaatscenario’s, 2020). In contrast to longer term predictions, such as those for 2050 and 2085, the 2030-scenario consist of only one scenario. Data regarding the average change in precipitation, sunshine ra-diation, evaporation and temperature are available.

The distributions based on historical data are adjusted for the changing aver-ages, found in the KNMI scenario. An important shortcoming is the lack of data regarding the change in variation. Climate change could have impact on the variability of weather circumstances. This could lead to more extreme re-sults, which are of special importance in the study of droughts.

(37)

3.9 Expected damage

We use two different methods to examine the expected drought damage. The first method uses data from characteristic years. The second method relates the indices to crop yield damages.

3.9.1 Characteristic years

First, we make use of data from the Agricom model. The required input data is not available to run the model for all studied years. However, the damage in the characteristic years 1949, 1959, 1969 and 1976 are calculated based on the Agricom model in Projectgroep Droogtestudie Nederland (2004). Therefore, we focus on the results in these characteristic years.

Based on the indices, the probability of exceedance of the specific droughts in the characteristic years can be evaluated. We use the results to estimate a curve representing the probabilities corresponding to each damage. Consequently, the expected damage can be calculated as the area below that curve.

3.9.2 Crop yields

Second, we analyse to what extend the indices can be used to estimate crop damages. Garc´ıa-Le´on et al. (2019) used a model with crop yields as dependent variable, and various drought indices as explanatory variable. Following this method, the crop yield anomalies need to be standardized.

Zt=

Yt− ¯Y

σY

, (28)

where Yt represents the yield of the crop in kilo per hectare at time t. We

find that Zt is indeed standardized for the considered crops, with mean of 0

and standard deviation of 1. Consequently, we will estimate the following OLS regression, for each drought index:

Zt= β1∗ Itw+ β2∗ (Itw) 2_{+ u}

t, (29)

where the drought index is denoted as Iw

t. We apply the model to the SPI,

MSPI and JDI index. β1and β2 are the coefficients corresponding to the index

and u represents the residual. Note that the constant is omitted, as both the indices and Z are standardized. A quadratic term is included, since we expect that the yield damage is increasing with more severe droughts. Also, we con-sidered the possibility that crops can endure a certain level of dryness, before they physically suffer from the sub-optimal growing conditions. Therefore, we included a dummy variable that indicates if a certain level of severity is passed. However, this indicator did not result in a better models, since both the AIC measure increased and the adjusted R2 _{value decreased. Therefore, we did not}

(38)

the indices, we can find out which interval has the most explanatory power. The coefficients are estimated using OLS.

Some crops show a trend in the yield anomalies, denoted by Z. In Figure 19 below, we see the difference between potatoes without trend, and sugar beets with a clear trend.

Figure 19: Standardized yield anomalies for potatoes and sugar beets. Therefore, we include an alternative regression that allows for a trend.

Zt= β1∗ Itw+ β2∗ (Itw)2+ β3T + ut, (30)

where T represents the corresponding year. Again, the coefficients are estimated using OLS. For both models, we tested the normality of the residuals using the Shapiro-Wilk test. The results shown in Table 5 show that the null hypothesis is not rejected. Therefore we assume that they are normally distributed.

Table 5: Normality test Shapiro-Wilk test W p-value Potatoes 0.977 0.894 Onions 0.953 0.501 Sugar beets 0.990 0.999

3.10 Market effects

(39)

approach in which participation in weather insurances is encouraged through subsidies. Moreover, many farmers want to insure themselves against growing uncertainty related to climate change. Therefore, the demand for weather in-surances has increased during the last decade. The so called ‘broad weather insurance’ is supported by the government, and compensates farmers in case of damage due to extreme weather. This includes for example droughts, storm and frost. Between 2010 and 2015, the amount of farmers that closed this in-surance doubled to over 1000 (Berkhout et al., 2016). After 2015, it increased further to almost 2000 in 2019 (van de Laak, 2020). As of this year, the govern-ment exempted the insurance from insurance premium tax, to encourage more farmers to protect themselves against uncertainty associated with weather cir-cumstances (Rijksoverheid, 2020). The additional costs of this measurement is yearly 7 million euros.

As the participation of the government and farmers in the broad weather in-surance increases, it is useful to obtain more precise estimates of the damages eligible for compensation. The largest supplier of the broad weather insurance, Vereinigte Hagel (VH), uses the current precipitation deficit (AAP D_Aprw ) as criterion for drought damage (van der Boom, 2020). Insured are eligible for compensation if the threshold of 250mm is reached. This method leads to an ef-ficient and fast payment system, however, it does not necessarily reflect the true damage. According to the historical damages in characteristic years, damages occur well before the 250mm threshold. Moreover, the damage is not strictly in-creasing with the precipitation deficit. Hence, the consequences of the imprecise measurements could be twofold. If the damage is underestimated, the insurance might not be attractive enough for farmers. The goal of the government to insure as many farmers as possible might be undermined. On the other hand, if the damage is overestimated the premiums are too low, and the investments of the government could be more efficient. Therefore, this section will try to give a more detailed approach about the construction of the financial damage for farmers.

(40)

Figure 20: Supply shock in a basic supply and demand model.

Droughts could lead to a shock in the supply curve, which will ceteris paribus lead to a higher price. Reinhard et al. (2015) describe the effects of the price change on surpluses. Part A in Figure 20 explains the decrease in surplus for the farmers, when the price effect is ignored. This corresponds to the damage calculated in the Agricom model. Part B is the increase in price that customers pay as a result of the lower supply. This surplus shifts from consumers to the farmers. Part C is the consumer surplus that is lost, since it is not shifted towards the farmers. This approach does not take the effect of export into account. When only the welfare effects in the Netherlands are considered, a separation should be made between the part of the supply that is consumed at home, and abroad.

(41)

lays the price effects of droughts. These factors make it difficult to quantify the impact of droughts on the price of crops.

(42)

Figure 21: Supply shock in a basic supply and demand model for the Dutch potato market between 2017 and 2018, with data obtained from: Agrimatie, Baltussen et al., 2014, and CBS (StatLine, Akkerbouwgewassen, retrieved from: https://opendata.cbs.nl/ on 16 April 2020).

Note that this is an indication for the effects of droughts, and not a precise estimation. Although the sewed area barely changed during these years, other factors than drought might have had impact on the supply. Also, the demand curve could adjust due to factors such as the popularity of the product. Part A is the largest, with an area of 2.49. Part B has an area of 1.20, and part C is as big as 0.17. This means that the area that farmers lose based on the Agricom model is clearly the largest. However, almost half of this damage is compensated by a price increase, corresponding to area B.

(43)

(44)

4 Results

First, we will describe and compare the general results of the different drought indices. Thereafter, the expected damage will be examined based on the various indices. Finally, we will show the results regarding the relationship between the indices and crop yield losses.

In order to get a representative view on the indices, we apply them to different characteristic years. In particular, we focus on the years 1976 and 2018. These years are both considered to be extremely dry, however, the construction of the droughts differs. Therefore, the indices could yield different results. The index values are not always comparable due to different scales. Consequently, we used both the original values, and their corresponding probability of exceedance. When we interpret the values in terms of probabilities, we need to keep two limitations in mind. First, the normality assumption is not always satisfied, for example for SPI values with ending month December. Second, the return peri-ods rely on the assumption that the hydrologic variables are independent and identically distributed. However, this assumption is violated if we consider the effects of climate change. Therefore, we included a section that incorporates the effects of an increase in precipitation and evaporation following the 2030 KNMI climate scenario.

4.1 KNMI

(45)

Figure 22: Construction of the accumulated depreciation deficit, with as peak the M AP DOct

Apr.

Since the years are ranked based on the peak, we find that 1976 is drier than 2018 according to the KNMI measure. However, it appears that the drought has not been solved at the end of the graph. Also, it is hard to judge the severity of the short term droughts within this time window. To meet these objections, we will also compare the years using the other indices. We start with the SPI.

4.2 SPI

(46)

Figure 23: SPI 1976 (left) and SPI 2018 (right).

We find that the structure of the droughts differs after the summer. The be-ginning of 1976 is drier, but the drought persists longer in 2018. In 1976, the months after July are wet. The SPI value for the time window between Oc-tober and December is above 1.5, which corresponds to severe wet conditions. In contrast, in 2018 the drought persist until after the summer. The window between September and December of 2018 is extremely dry. Since the SPI index is normally distributed, it is easy to interpret the data in terms of probabili-ties. During the first three windows of 1976 the SPI index is under −3. Three standard deviations from the mean correspond to a probability of exceedance below 0.15%. Hence, the probability of finding more severe conditions is less than 0.15%. This does not mean that the probability of three consecutive win-dows with these conditions correspond to a probability of exceedance of 0.153_%.

For these calculations, we need to incorporate the dependence between the time windows. This can be done using the JDI index.

4.3 JDI

The probability corresponding to the JDI index equals 0.035 for 1976 and 0.036 for 2018. Hence, according to this index the droughts are similar in both years. The probability of finding a combination of time windows with more severe conditions than 1976 is 3.5%. The corresponding return period equals 28.57 years. This means that we expect a similar drought approximately once every 29 years. The JDI values are similar because of the following reason. The start of 1976 is exceptional, but is followed by a regular autumn. The beginning of 2018 is exceptional, but not as rare as in 1976. However, the wet fall of 2018 makes the overall combination as exceptional as 1976. Hence, if we consider the dependence between the time windows the years are similar in terms of dryness.

4.4 MSPI

(47)

If we compress the data of all time windows in a single index that incorporates most of the variation, we find that 2018 is slightly drier than 1976.

4.5 Robustness

Figure 24 shows an overview of all indices in 1976 and 2018. For comparability, they are all converted to their corresponding probabilities. The probabilistic value corresponding to the KNMI index in 2018 is not available, since the most recent update of the KNMI dates from 2004. We find that the differences between the indices are large. A possible explanation is the reliance on the nor-mality assumption. The observations in 1976 and 2018 are in the far tails of the distribution. If the underlying data is not exactly normally distributed, devia-tions in the tail could be large. Therefore, we also include a index comparison between the years 1967 and 1996 in Figure 25. These years are characterized by the KNMI as ’average dry’ and ’moderate dry’ respectively. We find indeed that the differences are relatively smaller. We also included 10 non-characteristic years at the start of each decade. Figure 26 shows that the driest year again shows the most deviations. To get a better view on the extreme droughts of 1976 and 2018, we apply the extreme value models.

(48)

Figure 25: Index comparisson between 1967 (left) and 1996 (right).

Figure 26: Index comparison between ten non-characteristic years at the start of each decade of the 19th century. It does not show the KNMI index, since this is only available for characteristic years.

4.6 Extreme value models

(49)

both models, given that the deficit is already higher than 225mm. In words, this probability measures the chance that a very dry year is also extremely dry. For 1976 the probability is below 0,1% for both models. Since these are conditional probabilities, they cannot be compared with the other indices. However, we can conclude that 1976 is drier than 2018 based on both extreme value models.

Figure 27: SPI 1976 (left) and SPI 2018 (right).

4.7 Time windows

The flexibility of the time windows will now be used to focus on specific seasons. In August most crops are harvested, therefore it is the most relevant ending month from an agricultural perspective. Therefore, we will show the SPI values from the years 1976 and 2018 with August as ending month in Figure 28 below.

(50)

procreate. In drier circumstances, these attempts are more likely to fail (Natu-urmonumenten et al., 2019). Various species such as the grayling, northern crested newts and the European badger are feared to become extinct due to the droughts (Nieuwsuur, 2019). Moreover, when less plants bloom due to water deficits, there is less nectar for insects, and therefore less food for birds. All in all, large parts of the ecosystem will be put to hold (Dunning, 2019). May is a more relevant time for the ecosystem than August, and therefore, we will look at the SPI values ending in May. The results are shown in Figure 29.

Figure 29: Short and medium term windows for the SPI index in 1976 and 2018 ending in May.

We find that 1976 was increasingly dry. The time windows starting at the begin-ning of the year are mild, but the drought increased near May. In contrast, the SPI values in 2018 do not indicate exceptional conditions. Parties with special interest in the time window between April and May, such as certain animals, are particular exposed to drought in 1976.

Finally, droughts in autumn have multiple direct and indirect effects. Peat needs a constant amount of water during the year. If it falls dry, it breaks down (Natuurmonumenten, 2020). As a result, large amounts of CO2 are released.

(51)

Figure 30: SPI values for different time intervals in fall during 1976 and 2018. We find a clear contrast in the ending of 1976 compared to 2018. In 1976, the summer drought is directly followed by wet months. In contrast, the 2018 drought persisted during fall. This has the implication that the deficits in 2018 were not refilled in the winter, but carried on till 2019.

All in all, we find that the SPI index is useful in comparing droughts during dif-ferent time intervals. They can be interpreted in terms of probabilities, however, we need to keep in mind that the normality assumption is not always satisfied.

4.8 Economic damage

The characteristic years and their corresponding damages are shown in Table 6 and Figure 31. All drought indices are used to calculate the probabilities corresponding to the characteristic years. 1985 is used as benchmark for the KNMI index. It represents a level of dryness that is expected to be exceeded every year.

Table 6: Economic damage in characteristic years, and their corresponding probabilities according to different indices.

(52)

Figure 31: Expected damage in characteristic years.

In Figure 31 we find the probability distribution of the characteristic years based on the different indices, combined with the expected damage in those years. The results would be more precise with data from all years, but only the expected damage for characteristic years are publicly available. We see that the graph is not strictly decreasing for all indices. This could mean that there are droughts with high damage, that are more likely than droughts with less damage. We also find that in most years, the KNMI index leads to the lowest expected dam-age. The lines of the other indices are often above the KNMI line, indicating that the probability of the damage is estimated to be higher according to these indices. This could be due to the construction of the KNMI index. As it only incorporates the peaks of the droughts, it might underestimate the probability of droughts with longer duration or multiple peaks.

(53)

Table 7: Economic damage based on characteristic years. Index Expected damage (million euros)

KNMI 217,60

JDI 243,52

SPI 254,88

MSPI 242,68

4.9 Climate change

The effect of climate change differs for each index. The dotted lines in Figure 32 show the index values after incorporating an increase in precipitation and evaporation based on the 2030 climate scenarios. The SPI index is most prone to the changes, as it leads to a higher expected damage in almost all years. The changes for the other indices are neglectable. Note that the climate change effects are only measured using a general increase in precipitation and evapo-ration. These results indicate that the reaction to changing weather conditions differs for each index. Further studies could be used to examine the relationship more exactly.

(54)

4.10 2020

Finally, we apply the indices on the current year. According to NOS (2020), demand for water pumps has been extremely high among farmers this year, because it has been dry for weeks. To assess this particular drought, we compare the indices in Figure 33 below. We applied the last update on 28 June 2020.

Figure 33: The SPI, JDI and MSPI indices applied to the time interval between March and July 2020.

The SPI values show an increasing drought. The overall period between April and July has an SPI value just below -1.5. This corresponds to a severe drought. For the time windows starting from March, the value drops below -2. This implies that the drought is extreme, with a probability of exceedance below 1.75%. The MSPI and JDI summarize the four time intervals into a single index. The corresponding values are between -1.5 and -2. Although the drought is extreme, we cannot yet apply the POT and GEV models. These models are applied to the yearly droughts.

4.11 Crop yields

We use an OLS regression to relate crop yield damages to the different indices. To find which index results in the best fitting model, we use two different mea-sures. First, we find the explanatory power of the models using the adjusted R2_{measure. Thereafter, we provide the AIC measure. The values are shown in}

(55)

Table 8: Adjusted R2_{values of the models specified in (29) (without trend) and}

(30) (with trend) for each index with various time intervals.

Potatoes Onions (with trend) Sugar Beets (with trend) Adjusted R2 _{Adjusted R}2 _{Adjusted R}2

SP I_{J ul}Aug 0.7889 0.7372 0.844 SP I_{J un}Aug 0.2884 0.6157 0.8084 SP I_{M ay}Aug 0.3247 0.4167 0.8353 SP I_AprAug -0.04832 0.1411 0.7573 M SP I_{J un}Aug 0.6543 0.7248 0.8243 M SP I_{M ay}Aug 0.5903 0.7085 0.8445 M SP I_AprAug 0.4463 0.578 0.8464 J DI_{J un}Aug 0.6013 0.7316 0.8199 J DI_{M ay}Aug 0.5805 0.6394 0.8189 J DI_AprAug 0.3329 0.4495 0.7997

We find in Table 8 that the adjusted R2 _{of sugar beets is stable around 0.8, for}

(56)

Table 9: AIC values of the models specified in (29) (without trend) and 30 (with trend) for each index with various time intervals.

Potatoes Onions (with trend) Sugar Beets (with trend)

AIC AIC AIC

SP I_{J ul}Aug 27.22855 30.9629 24.12266 SP I_{J un}Aug 50.31589 37.42522 28.03698 SP I_{M ay}Aug 49.31901 44.51928 25.15363 SP I_AprAug 57.67574 51.09674 32.52367 M SP I_{J un}Aug 36.59637 31.74819 26.38432 M SP I_{M ay}Aug 39.8229 32.72475 24.0679 M SP I_AprAug 45.54841 39.01802 23.8329 J DI_{J un}Aug 39.30687 31.32169 26.85555 J DI_{M ay}Aug 40.27476 36.34471 26.95968 J DI_AprAug 49.08842 43.53596 28.88202

(57)

Table 10: Crop yield estimation, SPI index.

Potatoes Onions Sugar Beets (Without trend) (With trend) (With trend) SP I_{J ul}Aug 0.72433∗∗∗ 0.38039∗∗ 0.18517∗ (0.10290) (0.15114) (0.10444) (SP I_{J ul}Aug)2 _-0.08193 _-0.25585∗∗ _-0.07857 (0.05913) (0.10638) (0.07335) T ime -0.06216∗∗ 0.17388∗∗∗ (0.02648) (0.01769) Constant 1.02982∗∗ -1.64364∗∗∗ (0.34303) (0.22331) Adjusted R2 0.7889 0.7372 0.844

Standard errors in parentheses

∗_{p < 0.1,}∗∗_{p < 0.05,}∗∗∗_{p < 0.01}

We find in Table 10, that the coefficients for the SP I_{J ul}Aug index are statistically significant at a 10% level for all crops. The coefficient is positive, which implies that drier years are related to higher crop yield losses. The quadratic term is negative for all crops. Since droughts are indicated with a negative SPI value, this implies that severe droughts have relatively more impact on crop yields than moderate droughts. We found in Table 3 that extreme droughts have an SPI value below −2. The estimated yield anomalies corresponding to these conditions, are −2 ∗ 0.72433 + (−2)2_{∗ −0.08193 = −1.77638 according to this}

(58)

Table 11: Crop yield estimation, MSPI index.

Potatoes Onions Sugar Beets (Without trend) (With trend) (With trend) M SP I_{J un}Aug 0.57918∗∗∗ 0.30015∗∗ 0.08407 (0.11699) (0.11870) (0.09228) (M SP I_{J un}Aug)2 _-0.11661∗ _-0.27748∗∗∗ _-0.11356 (0.06436) (0.08465) (0.06520) T -0.08299∗∗ 0.16652∗∗∗ (0.02694) (0.01818) Constant 1.32360∗∗∗ -1.50339∗∗∗ (0.35861) (0.22561) Adjusted R2 0.6543 0.7248 0.8243

∗_{p < 0.1,}∗∗_{p < 0.05,}∗∗∗_{p < 0.01}

(59)

Table 12: Crop yield estimation, JDI index.

Potatoes Onions Sugar Beets (Without trend) (With trend) (With trend) J DI_{J un}Aug 0.62503∗∗∗ 0.37826∗∗∗ 0.12707 (0.12595) (0.11293) (0.09098) (J DI_{J un}Aug)2 _-0.16470∗∗ _-0.35954∗∗∗ _-0.12563∗ (0.07559) (0.09041 ) (0.07103) T -0.09085∗∗∗ 0.16462∗∗∗ (0.02693) (0.01848) Constant 1.48154∗∗∗ -1.49507∗∗∗ (0.36330) (0.23203 ) Adjusted R2 0.6013 0.7316 0.8199

∗_{p < 0.1,}∗∗_{p < 0.05,}∗∗∗_{p < 0.01}

A Flood of Droughts: New Perspectives on Droughts in the Netherlands