Biswarup Mukherjee
email: biswarup.mukherjee@mines-paristech.fr

Marcelo De Castro Fernandes

Luigi Vanfretti

A PMU-Based Control Scheme for Islanded Operation and Re-synchronization of DER

Keywords: Power systems modeling, PMU, Controlled islanding, Island operation, Re-synchronization, Power grid, DER, Modelica, OpenIPSL

This article proposes a novel synchrophasor-based control scheme enabling controlled islanding, islanded operation and automatic re-synchronization of a distributed energy resource (DER) in a distribution network. The performance of the proposed control system is studied using a test power system model. The DER controller uses a centralized architecture, receiving phasor measurement unit (PMU) measurements from both the transmission and distribution grids. In simulation, the frequency control function inside the proposed controller uses frequency estimates calculated using a new formula that exploits the bus voltage in rectangular form (real and imaginary values). The performance of the proposed frequency computation method is studied and compared with the conventional washout filter (WF) approach used by most power system software tools. The study also discusses why unwrapped bus angles are necessary to perform the automatic re-synchronization process. The performance of the proposed controller is evaluated using both deterministic and stochastic load models, allowing the assessment of variability in distribution grids. The implementation of the proposed control scheme and the simulation of the test system is carried out leveraging rich features of Modelica language and the Open-Instance Power System Library (OpenIPSL).

Introduction

Motivation

The concept of operating a smart grid includes the application of smart or intelligent equipment that makes power transmission and distribution secure, reliable and flexible. An important and desired function of smart grids is that of facilitating the autonomous operation of micro-grids [START_REF] Guoping | An overview of microgrid planning and design method[END_REF][START_REF] Zhou | An overview on microgrid technology[END_REF][START_REF] Morstyn | Control strategies for microgrids with distributed energy storage systems: An overview[END_REF], where distributed generators (DGs) are connected in the vicinity of, or co-located with, loads. Providing these two functions in today's power system is an important task considering the increasing penetration of DERs.

Nowadays, islanding and re-synchronization operations need to be executed occasionally by transmission system operators (TSOs) and/or distribution system operators (DSOs) to maintain proper balance between power generation and consumption. This is often carried out either with proprietary re-synchronizers that lack interoperability, or manually (i.e. without automation) where the operators coordinate changes at the DER site with plant engineers.

This article investigates how to automate the controlled islanding, islanded operation and re-synchronization processes for a synchronous generator-based DER located at the distribution network (DN). To this end, the article presents a centralized control architecture, models the necessary control functions and trip signals that need to be sent to protection relays. The control functions are envisioned to make use of measurements from PMUs [START_REF] Aminifar | Synchrophasor measurement technology in power systems: Panorama and state-of-the-art[END_REF][START_REF] Mahmood | Real-time reduced steady-state model synthesis of active distribution networks using pmu measurements[END_REF][START_REF] Perić | Optimal pmu placement for power system ambient data-based mode estimation applications[END_REF], located at transmission and distribution networks. This allows the monitoring of synchronization variables across the traditional TSO/DSO operation boundaries and with interoperability assuming they adopt standardized protocols such as IEEE C37.118.2 [START_REF]IEEE standard for synchrophasor measurements for power systems[END_REF].

Literature review and previous work

In current power grids, manual and autosynchronizers are widely used to perform re-synchronization in any islanded power network having multiple generators and loads connected to the grid [START_REF] Thompson | Fundamentals and advancements in generator synchronizing systems[END_REF]. Although manual systems are still used, they are mostly kept as a back-up and automatic systems are preferred due to its inherent greater precision to assess the criteria for synchronization [START_REF] Thompson | Fundamentals and advancements in generator synchronizing systems[END_REF][START_REF] Belyaev | Generator to grid adaptive synchronization technique based on reference model[END_REF][START_REF] Manson | Case study: Smart automatic synchronization in islanded power systems[END_REF]. Indeed, in order for synchronization to occur without unwanted consequences, it is necessary to closely monitor synchronization variables: voltage quantities (magnitude, phase-angle) and frequency, and to assure that they are within acceptable limits for the connection of two isolated systems to occur [START_REF] Belyaev | Generator to grid adaptive synchronization technique based on reference model[END_REF][START_REF] Manson | Case study: Smart automatic synchronization in islanded power systems[END_REF].

The problem of re-synchronization is not new and many different creative solutions have been proposed in the literature recently, since the topic has gained attention due to the bolstering deployment of DERs. In fact, a virtual synchronous machine (VSM)-based control strategy is presented in [START_REF] D'arco | A synchronization controller for grid reconnection of islanded virtual synchronous machines[END_REF], allowing microgrids to operate in both grid-connected and islanded without control loop reconfiguration. The proposed strategy acts on reducing frequency error from a VSM reference and phase-angle difference for a seamless re-connection of islanded grids. In [START_REF] Piardi | A control structure for smooth transfer from grid-connected to islanded operation of distributed synchronous generators[END_REF], a novel strategy for the smooth transition of DERs from grid-connected to islanded operation modes is investigated. In that paper, the proposed control action is independent on the islanding detection time. Moreover, distributed control of islanded microgrids are presented and investigated in [START_REF] Xin | Control of island ac microgrids using a fully distributed approach[END_REF] and in [START_REF] Dehkordi | Distributed resilient adaptive control of islanded microgrids under sensor/actuator faults[END_REF] for increased resiliency. The distributed approach in [START_REF] Xin | Control of island ac microgrids using a fully distributed approach[END_REF] allows microgrids to perform primary, secondary and tertiary control actions locally, reducing its dependency on external measurement. On the other hand, [START_REF] Dehkordi | Distributed resilient adaptive control of islanded microgrids under sensor/actuator faults[END_REF] uses the distributed architecture but proposes adaptive protocols that can compensate any adverse effect that might occur to the sensor.

It is important to observe that the previous solutions only contemplated the usage of local measurements. This is mainly because most of the previous work use inverter-based energy resources (IBRs) which are able to derive local frequency measurements from their phased-locked loops (PLLs), so that controlled islanding and post reconnection procedures can be performed quickly and precisely [START_REF]System protection schemes in power networks[END_REF][START_REF] Yang | Slow-coherency-based controlled islanding-A demonstration of the approach on the august 14, 2003 blackout scenario[END_REF]. However, there are different studies in the literature that present how remote sensing can be used for islanding and reconnection procedures, and this becomes attractive for more conventional generators that would also need to operate under such conditions. These remote sensing quantities may come from PMUs and, studies on the application for synchrophasors for these purposes has been proposed in [START_REF] Phadke | Synchronized phasor measurements in power systems[END_REF][START_REF] La Ree | Synchronized phasor measurement applications in power systems[END_REF][START_REF] Kyriakides | Synchronized measurements in power system operation: International practices and research issues[END_REF]. In [START_REF] Franco | Using synchrophasors for controlled islanding -a prospective application for the uruguayan power system[END_REF], a study analyzes how synchrophasors could be used to create a controlled islanding in the Uruguayan power system as a protective measure. As an example of a microgrid-related work, a strategy leveraging the remote sensing of both voltage and frequency quantities is proposed in [START_REF] Assis | Automatic reconnection from intentional islanding based on remote sensing of voltage and frequency signals[END_REF]. The control actions for automatic re-connection after an intentional islanding of a microgrid are based on conventional synchronism check relays and they are demonstrated to result in a smooth procedure.

The usage of both local signals and synchrophasor measurements would, then, be essential for a centralized architecture to function properly. In fact, it is absolutely necessary to adequately model the measurements, especially the frequency. A common approach in power systems tools is to use a WF in order to model the estimate of the bus frequency [START_REF] Ortega | Comparison of bus frequency estimators for power system transient stability analysis[END_REF]. However, this well-known solution is prone to numerical issues [START_REF] Ortega | Comparison of bus frequency estimators for power system transient stability analysis[END_REF] and, therefore, it is necessary to find an alternative approach. In [START_REF] Milano | Frequency divider[END_REF] the frequency divider approach is presented and it is compared to the WF approach in [START_REF] Ortega | Comparison of bus frequency estimators for power system transient stability analysis[END_REF], while this method is attractive, it is not well suited for generic modeling tools and has only been implemented in the tool reported in [START_REF] Ortega | Comparison of bus frequency estimators for power system transient stability analysis[END_REF]. In this paper, a more practical and simpler approach is proposed based on rectangular components of the voltage phasor. This is attractive because it becomes straightforward to implement in general purpose tools, but also, it follows the same principle that actual PMUs use to estimate the frequency internally, i.e. they estimate frequency based on the voltage angle [START_REF]IEEE standard for synchrophasor measurements for power systems[END_REF]. Furthermore, in previous work [START_REF] Mukherjee | Modeling of PMU-based islanded operation controls for power distribution networks using Modelica and openIPSL[END_REF], the authors proposed a control system architecture for islanded-operation, based on PMU measurements; however, that architecture did not provide the means to perform automatic resynchronization functions and to provide trip signals to relays to perform the required re-closing of the circuit breaker. In addition, the authors modeled an automatic re-synchronization function for DER generators [START_REF] Vanfretti | Automatic resynchronization controller analysis within a multi-domain gas turbine and power system model[END_REF], but it did not include an islanded-operation control function that would allow the islanded system to self-regulate its frequency. Therefore, this work proposes an integrated control architecture for both controlled islanding, islanded operation and automatic re-synchronization. The proposed controller, if adopted in the field, can operate taking measurements from synchrophasor data, hence providing an interoperable solution as compared to today's common practices.

Contributions

The main contributions of this article are:

(i) For modeling and simulation purposes, a new frequency estimate computation method based on the bus voltage phasor data in rectangular form is proposed. This aims to emulate frequency measurements, such as that obtained from PMUs, when used in a closed loop controller. For comparison, the Modelica [START_REF] Fritzson | Principles of object-oriented modeling and simulation with Modelica 2[END_REF] implementation of the conventional WF model [START_REF] Ortega | Comparison of bus frequency estimators for power system transient stability analysis[END_REF][START_REF] Milano | Frequency divider[END_REF] for frequency computation used in most power system tools, is presented.

The computed frequency values obtained using both the techniques are analyzed and compared.

(ii) An integrated control system architecture is proposed, which is interfaced in cascade with a turbine's governor and excitation controls (i.e. voltage regulator), without replacing the existing generator control structure. Within this control architecture, a control scheme for controlled islanding, islanded operation and autonomous re-synchronization is proposed.

(iii) Taking advantage of PMU measurements, an angle difference control function is proposed to influence the phase angle difference during the grid re-connection process. It is modelled inside the proposed control architecture, and compliments traditional voltage and frequency control functions.

(iv) The use of advanced object-oriented equation-based modeling language features from Modelica are exploited and illustrated, giving an example of how such advanced language features can facilitate design and analysis of new control functionalities, i.e. controlled islanding, islanded operation and re-synchronization.

(v) The performance of the control system is evaluated considering both deterministic, and stochastic load models. This analysis allows to understand the feasibility of the proposed control system's operation under different levels of time-domain variability.

(vi) Finally, to extend the analysis for effectiveness and robustness of the proposed control scheme, aspects related to the use of the PMUs for feedback are considered. The performance of the proposed scheme is evaluated for different values of the PMU reporting rate and PMU data transmission delay.

Paper structure

The reminder of this paper is organized as follows: In Section II, for simulation purposes, a new frequency computation approach is proposed and described along with the conventional method used by power system simulators. Section III presents the proposed control architecture and its Modelica implementation.

The power system model and simulation set-up is presented in Section IV. Four case studies are presented in Section V, while the conclusions summarizing the paper results are presented in Section VI.

Frequency Computation

When the distribution grid is disconnected from the bulk transmission grid to become an islanded system, the phase angle of the bus voltage in the distribution grid starts deviating from the main transmission system reference. As a result, the distribution network frequency also deviates from the transmission grid frequency. The frequency at the distribution grid can be controlled by a DER, if equipped with a controller for this purpose. If the frequency at the distribution grid is monitored with a PMU, this measurement can be used to control the distribution network's frequency. However, PMUs do not directly measure frequency. Instead, frequency is a derived quantity, obtained from estimates that use the bus voltage's rate of change in their phase angles [START_REF] Ieee | IEEE/IEC international standard -measuring relays and protection equipment -part 118-1: Synchrophasor for power systems -measurements, IEC/IEEE[END_REF]. The frequency f (t) estimation can be expressed as:

f (t) = 1 2π × dθ(t) dt = f 0 + 1 2π × dφ(t) dt , (1)
where the phase function θ(t) = 2πf 0 t + φ(t) includes both frequency and phase angle.

Theory and implementation

In this work, it is assumed that PMUs estimate frequency from the bus voltage angles. The PMUs are assumed to be located at both the transmission and distribution network will provide frequency signals to the proposed control system. However, this has to be carefully considered in modeling and simulation [START_REF] Ortega | Comparison of bus frequency estimators for power system transient stability analysis[END_REF][START_REF] Milano | Frequency divider[END_REF], because if the bus frequency is calculated from the bus angle directly, angle wrapping will produce sharp discontinuities in the derived frequency values. To overcome this issue and to provide useful and adequate frequency signals for control purposes, a new method is proposed herein. Here we have considered only two PMUs, one at the transmission grid and the other one at the distribution grid where the DER is located. Note that for the islanded operation, only one PMU is necessary in theory, but the re-synchronization requires both the PMUs. Consequently, we assume no other PMUs are installed throughout the system. In practice, at a minimum, at least two PMUs at the transmission grid and two PMUs at the distribution grid, all at different locations, should be utilized.

Let V i and V r be the imaginary and real parts of complex bus voltage V , then we have that

V = V 2 i + V 2 r , (2)
and the bus angle θ can be expressed as in [START_REF] Morstyn | Control strategies for microgrids with distributed energy storage systems: An overview[END_REF].

θ = arctan V i V r (3)
Therefore the angular frequency ω = 2πf can be estimated taking the first order derivative of the bus angle θ, as presented in [START_REF] Aminifar | Synchrophasor measurement technology in power systems: Panorama and state-of-the-art[END_REF].

θ = ω = V r Vi -V i Vr V 2 i + V 2 r ( 4
)
Note that for the quantity ω to be in per unit, it is necessary to divide it by the base frequency (i.e. 2πf 0 ). Equation ( 4) was implemented using Modelica language within a frequency computation model, as follows: Note that the derivative function blocks, denoted by der1 and der2, use a gain value of k=1 and a time constant of T= 0.01 sec. Since the proposed computation method relies on the real and imaginary values of complex bus voltage, the method is independent of rate of change of voltage phase angle. This makes the proposed method advantageous over the WF model, presented in the next section.

model f r e q u e

Washout Filter implementation in Modelica

In traditional power systems tools [START_REF] Siemens | pss/e-high-performance transmission planning and analysis software[END_REF], the WF model is used to estimate the grid frequency from the bus voltage phase angle. This model computes bus frequency by passing the voltage phase angle through a derivative block and a first order filter. This paper presents the Modelica implementation of the WF model shown in Fig. 1, which is deployed in the power system model. The parameter used for the WF model are as follows: k= 1, T f= 1 sec and T w = 2 sec, where T f and T w are the time constants of the derivative and first order filter blocks, respectively. The gain k remains same for both the derivative block and the first order filter models.

In this work we have performed the simulations using phasor time-domain modeling framework, it is to exploit both the components (real and imaginary) of the complex bus voltage phasor to compute the frequency, similar to the way this is performed by PMUs which derive the frequency from the bus voltage angle. Note that this approach is much different than those used in power electronic systems which derive the frequency from a PLL. Hence, because of the modeling framework and to be consistent with how PMUs derive frequencies, we propose the frequency computation using bus voltage phasor components.

Hence, in Section 2.1 we propose an alternative and more practical approach over the WF model.

Proposed Controller Architecture and Functions

Control architecture

The proposed control system architecture and internal functions for control islanding, islanded operation and automatic re-synchronization is presented in Fig. 2. The overall structure including sensing (PMUs), data communications and the controller's analog and digital outputs is shown in Fig. 2a, while the internal functions are presented in Fig. 2b. The controller has two main modes of operation: (1) controlled islanding and islanded operation (highlighted by green dotted line), ( 2) automatic re-synchronization (highlighted by orange dotted line).

The mode for islanded operation requires a PI controller, that can be used to restore the frequency deviation from the nominal value after the system is separated from the main grid. This function can be activated as soon as the island is detected.

The control mode for automatic re-synchronization consists of three major functions or units: computation unit, activation unit and control unit. The synchronizing variables are determined from two PMUs each at different remote sub-stations, one located at the transmission network and one in the distribution network. The synchronization variables are: voltage difference (∆V), frequency difference (∆f) and angle difference (∆θ). These three synchronization variables are determined by the computation unit within the centralized controller at the DER. The activation unit plays a key role in activating three different controllers (voltage controller, frequency controller and an angle-difference controller) following a sequential decision making process, that checks thresholds for three different synchronization variables (∆V, ∆f, ∆θ). When all three synchronization variables are below their thresholds after a given amount of time, the controller will send a trip signal to a relay to close the corresponding circuit breaker and to complete the re-synchronization process. quency control functions can be modeled using PI controllers, as reported in [START_REF] Assis | Automatic reconnection from intentional islanding based on remote sensing of voltage and frequency signals[END_REF]. To control ∆θ, which is the third re-synchronization variable, an angle 245 difference controller that requires a PID block is proposed herein. In the sequel, each individual control function is summarized, followed by a discussion on how they are used in the sequential control logic. Implementation aspects exploiting the Modelica language are summarized also.

Modeling of the automated re-synchronization control mode

Modeling 250

The voltage controller output is defined by Eq. ( 5), where x represents the Boolean input signal, (start voltage) that activates the controller when the re-synchronization process initiates.

y(x) =    K P V ∆V + K I V t 0 ∆Vdt , if x = true 0, otherwise (5)
where K P V and K I V are the proportional and derivative gains of the voltage control function.

y(x) =    K P f ∆f + K I f t 0 ∆fdt , if x = true 0, otherwise (6)
where K P f and K I f represent the proportional and integral gains of the frequency control functions, respectively. The Boolean input signal x activates the frequency controller after the voltage bounds are checked successfully. When all three (∆V, ∆f and ∆θ) bounds are checked successfully, the frequency controller is deactivated.

To monitor and control the phase angle of the bus voltage an angle-difference controller is used. The angle controller uses a PID function whose input is the unwrapped angle difference estimated from the transmission and distribution side bus voltage angles. The output of this controller corresponds to the following equation:

y(x) =    K P θ ∆θ + K I θ t 0 ∆θdt + K D θ d∆θ dt , if x = true 0, otherwise (7)
where K P θ , K I θ and K D θ represent the proportional, integral and derivative gain of the angle-difference control function respectively. When the re-synchronization process begins, a signal activates the control logic inside the activation unit. The outputs of the re-synchronization controller can be expressed using the following equations:

V rs = K P V ∆V + K I V t 0 ∆V dt (8a)
P rs = f P + f θ (8b) f P = K P f ∆f + K I f t 0 ∆f dt (8c)
f θ = K P θ ∆θ + K I θ t 0 ∆θdt + K D θ d∆θ dt (8d)
where V rs is applied to the supplementary input signal of the AVR, P rs is applied to the turbine's power set-point, while f P and f θ are applied to the governor's speed set-point. The controllers will change these outputs with the goal of minimizing the synchronization variables.

When the synchronization variables (∆V, ∆f and ∆θ) are below their thresholds for a pre-defined amount of time, the Boolean output signal is sent to a protection relay that would close the circuit breakers of the tie-line, which finishes the re-synchronization of the islanded distribution network with the transmission network.

The sequential control logic in Fig. 3 was implemented within activation unit, where limit checking blocks are modeled to monitor the synchronizing variables, with the following thresholds:

∆V min ≤ ∆V ≤ ∆V max (9a) ∆f min ≤ ∆f ≤ ∆f max (9b) ∆θ min ≤ ∆θ ≤ ∆θ max (9c)
where the ∆x min and ∆x max values for x : V, f & θ need to be determined via simulation. The minimum value ∆θ min must be larger than the angle of the equivalent impedance between the tie-line linking the transmission and distribution network. This implies that ideally, the re-connection should be performed when there is minimal power being fed through that line, resulting in a smooth coupling between the systems. Similarly, ∆V max would be the largest possible voltage drop through that line, while ∆f min, max can be targeted to have a low error, for example ±0.01 Hz.

Modelica implementation

To avoid the numerical issues during simulation, advanced Modelica language features are used to handle phase angle wrapping before it was applied to the ∆θ controller. The controller uses the updated difference value theta diff new estimated from the angular difference between the transmission and distribution network bus voltages. The calculation of theta diff new uses different Modelica operators such as homotopy, smooth, noEvent [START_REF] Tiller | Introduction to physical modeling with Modelica[END_REF][START_REF]Modelica operators & expressions[END_REF], as follows: For the homotopy operator the argument simplified = theta diff implies that the initialization process uses this simplified expression, thereby making the initialization process more robust. Meanwhile, the smooth operator is provided an argument 0 to indicate that the variable theta diff new is discontinuous, i.e. non-differentiable. This allows to model the actual behavior for this control input variable as seen in the field [START_REF] Vanfretti | A phasor-data-based state estimator incorporating phase bias correction[END_REF], where the measured phase angles will be discontinous when the system undergoes a line trip and similar discrete changes. Finally, the noEvent operator helps with event suppression. This minimizes the number of states events that can arise during simulation. State events require reinitialization of the differential-and-algebraic equation model, and thus, degrade simulation performance [START_REF] Henningsson | DAE solvers for large-scale hybrid models[END_REF]. Therefore, noEvent helps avoiding such events, improving simulation performance.

Modeling of the controlled islanding and island operation mode

Once the distribution network is islanded, both power balance and frequency would need to be controlled by at least one of the DERs in the network. In this work, only the case when a single DER is located in the islanded portion of the network is considered. It is assumed that the DER consists of a conventional synchronous generator with both speed and voltage controllers. To achieve autonomous operation as an island, the DER is equipped with the proposed controller architecture described above. To achieve autonomous island operation, the proposed controller at the DER would be equipped with an ICL (Intelligent Control Logic) for islanded operation, which that can use real-time frequency measurements from PMUs that can be placed in any location within the distribution grid where the DER is present.

The ICL is modeled within the distribution network generator (G22) (see Figs. 5 and6) that uses a PI controller in order to maintain zero frequency deviation when the islanding occurs. The output of the PI controller within ICL can be expressed by Equation [START_REF] Manson | Case study: Smart automatic synchronization in islanded power systems[END_REF], where K P I and K I I are the proportional and integral gain respectively of this PI function and ∆Ω L represents the error in the distribution network frequency from its nominal value. The implemented algorithm within ICL disables the PI controller modeled inside when the automatic re-synchronization process begins and, at the same time, the ResynchUnit initiates the voltage control action. Hence, the re-synchronization maneuvering process for an isolated distribution grid is carried out autonomously.

Ω rs =    K P I ∆Ω L + K I I t 0 ∆Ω L dt , if x = true 0, otherwise (10)

Modeling of the control mode transitions

The transitions between different control modes used by the proposed controller are shown in Fig. 4. Note that it shows labels indicating how the scheme executes the transitions from one control mode to another.

For illustration, assume that a DER with the proposed control system is islanded. The control system will either detect an island or the island will be created and activate the islanded operation controller (see label (A) in Fig. 4). After the system has been operated in islanded mode under acceptable conditions, the re-synchronization controller can be activated as shown by (B) in Fig. 4, which concurrently disables the islanded operation controller (see (C) of Fig. 4). As soon as the re-synchronization controller starts operating, the sequential control mode logic (see Fig. 3) is activated. Naturally, once the re-synchronization process is completed, the re-synchronization controller is deactivated, leaving only the conventional controls (AVR and Governor) to drive the DER.

Power System Modeling and Simulation Set-up Implementation

The power system model, shown in Fig. 5, was implemented using Modelica language and the OpenIPSL [START_REF] Vanfretti | itesla power systems library (iPSL): A Modelica library for phasor time-domain simulations[END_REF][START_REF] Baudette | OpenIPSL: Open-instance power system library -update 1.5 to iTesla power systems library (iPSL): A modelica library for phasor time-domain simulations[END_REF] in Dymola [START_REF] Henningsson | DAE solvers for large-scale hybrid models[END_REF][START_REF]The Dynamic Modeling Laboratory software (Dymola)[END_REF] in order to study the performance of the proposed automatic re-synchronization controller. The controller architecture was placed inside the DER generator model (G22) at the distribution network, as shown in Fig. 6. The transmission network generator model (G1) uses a GENSAL block that corresponds to the synchronous generator model with salient poles and a simple excitation control system represented by the SEXS model. To start the islanding process, the entire system frequency was perturbed by introducing a speed change in the hydro-turbine and governor system model (HYGOV) inside G1. The HYGOV is a standard hydro turbine governor model, and it is the same as the one in PSS/E.

For the DER at the distribution network the different components of the plant are shown in Fig. 6. An IEEESGO block was used to model the gas turbine and governor system within G22. The IEESGO model is an IEEE standard model for turbine-governor systems, which can be used to represent both reheat steam turbines or simplified gas turbines. To model the excitation control system the SEXS block was used. The other two components of G22 are the ResynchUnit and the ICL blocks, as explained in Section III.

When the islanding begins, the controlled circuit breaker (CB2) islands the distribution grid from the transmission network and at the same time a signal is sent to the ICL, activating a PI block within this model. The algorithm to automate the switching from islanded mode to re-synch mode and the operation Modeling of the stochastic variability of the aggregate load was carried out by injection of white noise to the PQ load model, as highlighted in green in Fig. 5. Note that white noise generates a signal having normal distribution characterized by a mean and variance.

A simulation set-up block1 was used for two purposes. First, it helps to create the conditions that result in the islanding of the distribution grid eight seconds after the simulation begins. Second, to create the conditions to perform the automatic re-synchronization operation of the same distribution network with the transmission grid using the proposed control scheme. The frequency computation blocks are modeled inside the simulation set-up, that act as PMUs to compute both transmission and distribution network bus frequencies using the proposed method explained in Section II.

Case Studies

We considered the power system model as shown in Fig. 5 for the following case studies were performed using Dymola [START_REF]The Dynamic Modeling Laboratory software (Dymola)[END_REF] software. For all simulations, DASSL [START_REF]A description of DASSL: A differential/algebraic system solver[END_REF] numerical solver was chosen as an integration algorithm with a tolerance of 0.0001 having an output/reporting interval of 20 ms.

Case 1: Comparative analysis between the WF model and the proposed frequency computation technique

In Section II the Modelica implementation of the WF model was presented. This section compares the performance of the proposed frequency computation approach with the traditional frequency computation approach using the WF filter model while performing the automatic re-synchronization process.

The main differences between the WF filter model and the proposed model is that the WF uses the bus (B6) voltage phase angle while the proposed technique uses the real and imaginary parts of the same bus (B6) voltage to compute the distribution network frequency. Figure 7 shows that the WF model introduces transients in the frequency calculation (red trace) due to bus angle switching, whereas the proposed technique provides a smooth frequency output (blue trace).

For both the models the outputs of the frequency controller are plotted in Fig. 8. The plot shows, due to frequency calculation provided by WF model, the frequency controller fails to operate successfully and introduces spikes or transients (red trace) in its output. On the other hand, the blue trace shows that the controller operates successfully by producing the expected steady state error signal, when it has the input frequency computation using the proposed approach. This study shows the advantage of using the proposed frequency computation technique over the conventional WF method for frequency computation.

Case 2: Controlled scheme performance analysis

This case study was performed to evaluate the performance of the proposed control scheme for both PQ load and an infinite bus connected to the transmission network side of the developed power system model, discussed in the previous section. Simulation results are plotted in Fig. 9-16 for 5 MW power dispatch from the distribution network generator. Figure 9 plots the distri- bution network frequency deviation for both a controlled and an uncontrolled islanding event executed before the automatic re-synchronization process by ResynchUnit. Observe from Fig. 9 that the distribution network frequency deviation returns to zero (the red trace) for the controlled islanding operation, while in case of uncontrolled islanding the control architecture fails to do so (shown by the blue dotted trace). This demonstrates the functionality of the islanded operation controller as shown in Fig. 4. Figure 10 depicts the frequency deviation in both transmission and distribution network when a PQ load is connected to transmission grid at bus B5. This is performed to illustrate the difference of modeling the reminder of the network with an infinite bus and how that influences the system behavior. Figure 10 clearly exhibits how frequency controller acts to minimize the frequency difference when the re-synchronization process initiates. Note that when the re-synch starts, the control functions are applied to first minimize the voltage difference (∆V) (see Fig. 11), followed by the ∆f control. Note that because the transmission side (TS) does not have automatic generation control (AGC), the frequency deviates according to G1's droop, and the ∆f controller drives the G22 to match the distribution network frequency with the TS frequency instead of returning it to a fixed reference frequency. This is important as the TS frequency is never exactly at its reference value (i.e. 50 or 60 Hz), which means that the proposed controller will be able to track the TS instead of a fixed frequency reference.

Figures 12 and13 show the frequency deviation and voltage profile of the power system due to the re-synchronization operation when the transmission grid remains connected to an infinite bus. Observe from Fig. 12 that as the transmission network frequency swings, but returns to 50 Hz after both the islanding and the re-synchronization. Hence, this is an ideal scenario; because the transmission grid includes an infinite bus that establishes the reference frequency in the TS, and, as a result, the distribution grid frequency returns at 50/60 Hz before and after the re-synchronization process.

From Fig. 13 it can be observed that as soon as the islanding begins, the bus voltage levels at B4 and B6 are affected, as a result there exist the voltage difference until the voltage control action begins. After voltage control action ends, the ∆V≈ 0 as it can be observed from the green trace. This ensures that the voltage controller reduces the control error to the desired threshold. It can also be observed from the simulation result that the re-synchronization process does not affect the bus voltages or introduce switching transients, this is due to the stiffness of the infinite bus at the transmission network. trol scheme for both IEESGO and HYGOV models. The results validate the efficacy of the proposed scheme when two different turbine and governor systems are modeled within the distribution generator. Figure 14 shows the frequency deviation in the distribution grid for different gains (i.e. K P V and K I V ) of the voltage controller within the re-synchronization function. It can be observed for IEESGO the transient performance can be improved (as seen from the red trace, the overshoot is reduced ≈ 16%) by tuning the control parameters. A similar performance can be noticed in Fig. 15 when the transmission network is connected to an infinite bus. Here we observe that the change in gain parameters of the voltage controller results in different re-synchronization times. Frequency deviation (Hz)

Hygov (k P V = 0.05, k I V = 2) IEESGO (k P V = 0.05, k I V = 2) IEESGO (k P V = 0.1, k I V = 2)
re-synchronization begins (see B in Fig. 4) frequency control begins (see Fig. 3) re-synchronization executed For different permanent droop of HYGOV turbine and governor system Fig. 16 plots both transmission and distribution network frequency deviations. It is interesting to observe that when the droop is 1, the control loop becomes unstable, and as a result, the controller fails to re-synchronize two grids; however, with a droop value 0.5 the controller successfully re-synchronizes two islanded grids. These results show that it is necessary to tune all control parameters, not only the ones of the new controller, but other control systems in cascade and involved in the re-synchronization process. Hence, for each kind of power plant, the parameters need to be carefully calibrated considering the turbine and governor time constants and dynamics.

Case 3: Effect of stochastic load model on control action

This case study is performed to analyze the performance of the proposed controller for a stochastic load model connected to the distribution network using the Modelica noise modeling [START_REF]Modelica Noise[END_REF]. The stochastic load model simulates the behaviour of consumer demand, showing variability when load forecasting is performed [START_REF] Zhao | Probabilistic voltage stability assessment considering stochastic load growth direction and renewable energy generation[END_REF]. This also allows considering the time-domain variability of behind-the-meter and other non-dispatchable energy sources at the distribution network level. Therefore it is necessary to investigate whether the islanded and automatic re-synchronization control actions can be executed together under aggregate load variability conditions. The stochastic load model has been simulated with standard deviation (σ) values 0.00004 and 0.000025; and a period of 20 ms.

For simulation execution we consider that the noise generation begins at 5 s. As a result both the transmission and distribution network frequency becomes stochastic as long as they are not isolated from each other. As shown in Fig. 17, after the distribution network is intentionally islanded the control action begins and successfully restores the distribution network frequency deviation to zero.

It is also interesting to observe that the transmission network frequency is deterministic in nature as long as it remains isolated from the islanded distribution grid which includes the stochastic load model. Re-synchronization occurs at 225 s, and, again, it is observed from Fig. 17 that both transmission and distribution network frequency becomes stochastic, which is as expected. Therefore this case study verifies that the proposed control architecture will successfully operate in the case of load uncertainties or non-deterministic variability of nondispatchable energy sources, which are becoming more common in distribution networks.

Case 4: Effect of delay control action

In this section, we analyze the impact of both PMU reporting rate and data transmission delay on the performance of the proposed control scheme. Here we consider an HYGOV turbine and governor system within the DER model and, the transmission grid remains connected to the infinite bus.

To model different PMU reporting rates a ZOH from MSL (Modelica Standard Library) [START_REF]Modelica Standard Library[END_REF] was used. Observe from Fig. 18 that delays from 50 ms to 150 ms have no major impact on the controller's performance; this is because the frequency dynamics being controlled are much larger than typical PMU reporting rates. This is a positive result, as typical PMU reporting rates are ≤ 16 s, i.e. 10 to 60 samples per second.

On the other hand to mimic the aggregate time-delay from a PMU device a "fixedDelay" block from MSL was used in the distribution network before this frequency was applied to the controller. It is interesting to observe in Fig. 19 that as the data transmission delay increases (the blue trace) the control action gets delayed which eventually increases the re-synchronization time. These results are positive, as typical synchrophasor systems only incur delays in the order of 100s of milliseconds, up to a few seconds.

Conclusion 530

An integrated and centralized control system architecture for DER operation was proposed in this work and tested in a power system model. For modeling and simulation purposes, a frequency computation approach was presented and used to emulate frequency measurements from PMUs that are utilized as input to different functions of the proposed scheme.

535

The proposed centralized control system architecture is capable of both operating the DER in an island thereby avoiding the need of disconnecting it, and also capable of performing re-synchronization by executing the maneuver-ing process autonomously. Modeling the proposed control system architecture within one DG we aimed to show that, using such control architecture, it is possible to use remote measurements from synchrophasor data from different sub-stations to re-synchronize the DER autonomously. Moreover, it will only require the new control system to be installed in cascade with existing technology. In contrast, using existing technologies, it is only possible to re-synchronize the system with a separate system called a re-synchronizer [START_REF] Thompson | Fundamentals and advancements in generator synchronizing systems[END_REF]. Similar to the author's proposal, other methods for controlled islanding [START_REF] Franco | Using synchrophasors for controlled islanding -a prospective application for the uruguayan power system[END_REF] and synchronization [START_REF] Almas | A hybrid synchrophasor and goose-based power system synchronization scheme[END_REF][START_REF] Assis | Automatic reconnection from intentional islanding based on remote sensing of voltage and frequency signals[END_REF] have only been explored in simulation or lab experiments and have not been used in the field, and moreover, they also only address the control of a single generator/ DER.

Here we have considered the PMUs located in both the transmission and the distribution network. The best locations for these PMUs in terms of their placement and delay are out of the scope of this work and will be subject to future research. By using measurements from two parts of the grid, we show how an integrated control system architecture that is interfaced in cascade with a turbine's governor and excitation controls, can effectively perform two functions for which separate systems are proposed in the literature, all while at the same time without requiring the replacement the existing generator controls, but instead, complimenting them.

In this work we have only considered the grid constraints as the control/synchronizing variables which are useful for controller modeling; however, additional grid constraints (e.g. line ampacities) and cyber-physical aspects have not been considered in the modeling context.

Finally, the method can be expanded to consider multiple DERs in the distribution network, as well as distributing the control architecture to support such case. However, expanding beyond the control of a single DER is more complex. It would require an additional layer of coordination that allocates the control effort to each DER. Such work while of interest, it is outside of the scope and intent of this paper and is subject for future research.

Further work may involve tuning controllers inside the control scheme. For this purpose, the controllers inside the control scheme should be re-designed for robustness to improve the overall performance.

  n c y C a l c u l a t i o n B l o c k M o d e l i c a . B l o c k s . I n t e r f a c e s . R e a l I n p u t u ; M o d e l i c a . B l o c k s . I n t e r f a c e s . R e a l I n p u t u1 ; M o d e l i c a . B l o c k s . I n t e r f a c e s . R e a l O u t p u t y ; M o d e l i c a . B l o c k s . C o n t i n u o u s . D e r i v a t i v e d e r 1 ( k=1,T=0.01) ; M o d e l i c a . B l o c k s . C o n t i n u o u s . D e r i v a t i v e d e r 2 ( k=1,T=0.01) ; equation y = ( u * ( d e r 2 . y )-u1 * ( d e r 1 . y ) ) / ( ( u ˆ2 ) +(u1 ˆ2 ) ) ; connect ( u 1 , d e r 2 . u ) ; connect ( u , d e r 1 . u ) ; end f r e q u e n c y C a l c u l a t i o n B l o c k ;

Figure 1 :

 1 Figure 1: Washout Filter (WF) implementation in Modelica.

Figure 3

 3 summarizes the algorithm for the sequential control logic for resynchronization. During the re-synchronization process, the voltage and fre-Centralized controller architecture at the DER using remote PMU measurements. islanding and operation Automatic re-synchronization (b) Internal control functions.

Figure 2 :

 2 Figure 2: Controller Architecture and Functions.

Figure 3 :

 3 Figure 3: Sequential control logic for re-synchronization.

  equation t h e t a d i f f = (-B 6 . a n g l e ) + B 4 . a n g l e ; t h e t a d i f f n e w = homotopy ( a c t u a l = smooth ( 0

Figure 4 :

 4 Figure 4: Control scheme transitions between modes of operation. Grey lines and fill colors indicate existing conventional generator control system functions, while new functions and signals are shown in black and other colors.

Figure 5 :

 5 Figure 5: The power system model.

Figure 6 :

 6 Figure 6: Centralized control structure implementation within the generator model (G22).

Figure 7 :

 7 Figure 7: Frequency computation comparison between the proposed and the WF approaches.

Figure 8 :

 8 Figure 8: Output of frequency controller with both proposed and WF method.

  in the frequency deviation

Figure 9 :

 9 Figure 9: Controlled islanding vs. uncontrolled islanding.

Figure 10 :

 10 Figure 10: Frequency deviation in both transmission and distribution network for a PQ load model connected to the transmission network.

Figure 11 :

 11 Figure 11: Transmission and distribution network bus voltages for a PQ load model connected to the transmission network.

Figure 12 :

 12 Figure 12: Frequency deviation in both transmission and distribution network for an infinite bus connected to transmission grid.

Figure 13 :

 13 Figure 13: Transmission and distribution network bus voltages for an infinite bus connected to transmission grid.

Figures 14 -

 14  are plotted to compare the performance of the proposed con-

Figure 14 :

 14 Figure 14: Distribution network frequency deviation for both IEESGO and HYGOV, for a PQ load model connected to the transmission network.

Figure 15 :

 15 Figure 15: Frequency deviation in both transmission and distribution network for an infinite bus connected to transmission grid, with an IEESGO model in DER.

Figure 16 :

 16 Figure 16: Frequency deviation in both transmission and distribution network for an infinite bus connected to transmission grid, with an HYGOV model in DER.

Figure 17 :

 17 Figure 17: Transmission and distribution network frequency deviation for a stochastic load model in the distribution network.

Figure 18 :Figure 19 :

 1819 Figure 18: Distribution network frequency deviation for different PMU reporting rates.

  The output of the voltage controller is applied to the auto-fed to the PI block. If this Boolean signal is false, the output of the controller becomes zero.

	The output of the frequency controller modulates the turbine's power error
	signal. Similar to voltage controller, the frequency controller output can be
	defined by Equation (6):

255 matic voltage regulator (AVR) error signal, which controls the generator field voltage. If the Boolean input signal, start voltage, is true, the computed volt-age difference is

For brevity the detailed modeling of this block is not presented here.

Acknowledgements B. Mukherjee is supported by ARMINES, at Centre PERSEE of MINES ParisTech through the EVA project. EVA is funded in the framework of the joint programming initiative ERA-Net Smart Energy systems' focus initiative Integrated, Regional Energy Systems, with support from the European Union's Horizon 2020 research and innovation program under grant agreement No 775970. M. Fernandes and L. Vanfretti are supported by the New York State Energy Research and Development Authority (NYSERDA) and the Engineering Research Center Program of the National Science Foundation and the Department of Energy under Award EEC-1041877, in part by the CURENT Industry Partnership Program.

Authorship contribution statement

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.


Keywords: Mixture model, Generalized pareto distribution, Drought hazard, Extreme value theory JEL Codes C15, Q54

A critical stage in drought hazard assessment is the definition of a drought event, and the measure of its intensity. Actually, the classical approach imposes to all climatic region the same set of thresholds for drought severity classification, hence resulting in a loss of information on rare events in the distribution tails, which are precisely the most important to catch in risk analysis. In order to better assess extreme events, we resort to an extreme value mixture model with a normal distribution for the bulk and a Generalized Pareto distribution for the upper and lower tails, to estimate the intensity of extreme droughts and their occurrence probability. Compare to the standard approach to drought hazard, which relies on a standardized precipitation index and a classification of drought intensity established from the cumulative standard normal distribution function, our approach allows the drought threshold and the occurrence probability of drought to depend on the specific characteristics of each precipitation distribution. An application to the West Africa region shows that the accuracy of our mixture model is higher than that of the standard model. The mixture performs better at modelling the lowest percentiles and specifically the return level of the centennial drought, which is generally overestimated in the standard approach.

Introduction

This study takes place in the strand of literature dedicated to drought risk assessment (e.g. Carrao et al 2016). In sub-Saharan Africa, drought is a major and recurrent threat for people, economic activity and environment. Rainfed agriculture is the dominant production system so that the agricultural sector, which contributes for a large part to the gross domestic product, is directly exposed to dry spells. In these low-income countries, droughts have long-lasting effects that can be felt long after weather conditions have returned to normal. Human and physical capital losses caused by drought often result in protracted food crisis, and contribute to poverty trap. Drought events have also spillover effects on regions or countries that are not directly affected by the hazard but are connected by trade. Food price spikes resulting from production shortfalls tend to propagate to trading partners that import or export the food products. Therefore, in sub-Saharan Africa, severe droughts are associated to food insecurity, even famine and human loss, conflicts and economic slowdown (O'Grada 2007). This was the case in the beginning of the 70s and 80s, and more recently in 2010-2011, three periods of extreme drought that have marked the continent's history.

Despite the pervasive impact of droughts in Africa, drought risk management strategies have long been limited to ex-post emergency responses [START_REF] Wilhite | Managing Drought Risk in a Changing Climate: The Role of National Drought Policy[END_REF]. This situation led the United Nations Convention to Combat Desertification to plaid for urgent actions at the global, regional and national level, to establish a strategic framework for drought risk management in Africa (UNCCD, WMO and FAO 2018). To define a consistent framework for drought risk management, and enhance resilience to drought, one prerequisite is to identify and understand the drought risk components [START_REF] Hayes | Reducing Drought Risk: Bridging Theory and Practice[END_REF]).

The risk associated with a natural hazard is defined as the potential loss from the hazard. It is commonly view as the combination of three components: hazard, exposure and vulnerability [START_REF] Dilley | Natural Disaster Hotspots: A Global Risk Analysis[END_REF][START_REF] Ipcc | Managing the Risks of Extreme Events and Disasters to Advance Climate Change Adaptation. A Special Report of Working Groups I and II of the Intergovernmental Panel on Climate Change[END_REF]). Therefore, depending on the vulnerability of exposed socioeconomic and natural systems, a natural hazard may become a disaster.

Hazard, the first component of risk, can be defined as the probability of occurrence of a potentially damaging event [START_REF] Rajsekhar | Integrated Drought Causality, Hazard, and Vulnerability Assessment for Future Socioeconomic Scenarios: An Information Theory Perspective[END_REF]. The difficulty in measuring drought hazard comes from the multiplicity of drought definitions. There is no unique definition of drought, and drought assessment is necessarily region and impact specific [START_REF] Wilhite | Understanding: The Drought Phenomenon: The Role of Definitions[END_REF]. As other natural hazards, droughts can be characterized in terms of their severity, location, duration and timing [START_REF] Wmo | Integrated Drought Management Programme (IDMP), Integrated Drought Management Tools and Guidelines Series[END_REF]. But droughts differ from other natural hazards in many respects. Droughts are slow-onset, large-scale and long-duration phenomena. They result from large-scale anomalies in atmospheric circulation, require a minimum of two to three months to become established, and persist for months, seasons or years. Therefore, the onset and end of droughts, their intensity and spatial characteristics are difficult to assess.

Drought is also a socioeconomic construct, which definition relates to its impact. The climatological community has defined four types of drought: meteorological drought, agricultural drought, hydrological drought, and socioeconomic drought [START_REF] Wilhite | Understanding: The Drought Phenomenon: The Role of Definitions[END_REF][START_REF] Wmo | Integrated Drought Management Programme (IDMP), Integrated Drought Management Tools and Guidelines Series[END_REF]. Meteorological drought is defined as a prolonged dry period in the natural climate cycle. Thus, meteorological drought is a relative, rather than absolute phenomenon, specific to a region, that can occur in any part of the world whether arid or wet [START_REF] Wilhite | Drought as a natural hazard: concepts and definitions[END_REF]. It is measured by departure from normal precipitations over some specified period of time. Agricultural droughts originate from soil water depletion during a specific growing stage of the crop leading to reduced yields. A commonly used index to assess agricultural droughts is the Standardized Precipitation Evapotranspiration Index (SPEI), which combines precipitation and evapotranspiration indicators to catch soil moisture [START_REF] Beguería | A Multiscalar Drought Index Sensitive to Global Warming: The Standardized Precipitation Evapotranspiration Index[END_REF]. Hydrological drought refers to persistent low surface or subsurface water supply, in streams, rivers, reservoirs, or groundwater. The Palmer Drought Severity Index [START_REF] Palmer | Meteorological Drought. Research Paper No. 45[END_REF]) is used to assess this type of long-term drought, which usually results from a prolonged meteorological drought. The index combines temperature and water balance indicators. Socioeconomic drought occurs when the demand for water exceeds the supply for domestic use or economic activity (for instance hydroelectric energy production or irrigation).

All types of drought originate from a deficiency of precipitation, i.e. a meteorological drought that propagates to other types of drought (agricultural, hydrological and socioeconomic) [START_REF] Wilhite | Drought as a natural hazard: concepts and definitions[END_REF]. Whatever the type of drought, the intensity of drought is not evaluated upon its social or economic impact but in relation to the climatologic norm for the region or to the demand for water in the region or sector of activity. Droughts intensity being generally the main, or even the sole, criterion used to assess drought hazard, the measurement of drought intensity plays a crucial role in estimating drought risk The reference classification system for drought severity is that of [START_REF] Mckee | The relationship of drought frequency and duration to time scales[END_REF]. It has been established for a standardized index, namely the Standardized Precipitation Index (SPI). The SPI measures cumulated precipitation deviations from their long-term mean. It can be computed over several intervals of time e.g. 3, 6, 12 months or more. Based on precipitations, the SPI specifically assesses the intensity of meteorological droughts but depending on the considered time scale, it can also be informative on other types of drought. The raw precipitation data are fitted to an incomplete Gamma distribution and then transformed to Gaussian equivalents with zero mean and variance equal to one (see for more computational details [START_REF] Guttman | Accepting the Standardized Precipitation Index: A Calculation Algorithm 1[END_REF][START_REF] Edwards | Characteristics 20th century drought in the United States at multiple time scales[END_REF]. Therefore, SPI values are dimensionless, comparable across regions with different climatic regimes, and they offer a readily interpretable measure of drought intensity. [START_REF] Mckee | The relationship of drought frequency and duration to time scales[END_REF] distinguished four types of droughts from mild to severe, according to the SPI values and a fixed set of threshold values.

Working with a standardized index of drought, such as the SPI, and using the classification of [START_REF] Mckee | The relationship of drought frequency and duration to time scales[END_REF], the standard procedure to compute a drought hazard index (DHI) consists in giving weights to each drought category. Then the DHI is the weighted sum of the frequency of occurrence of each type of drought. Alternatively, each drought category can be further split into subcategories corresponding to different occurrence probabilities, from low to high, with increasing ratings. Thus, the DHI is calculated by combining the weights and ratings of the various subcategories of droughts. The sub-classification can be arbitrarily chosen [START_REF] Shahid | Drought Risk Assessment in the Western Part of Bangladesh[END_REF] or identified using a natural break classification method (Carrao et al 2014;[START_REF] Rajsekhar | Integrated Drought Causality, Hazard, and Vulnerability Assessment for Future Socioeconomic Scenarios: An Information Theory Perspective[END_REF]He et al 2013).

The main advantage of this approach to drought classification and hazard estimation, based on a standardized index, is to offer a comprehensive and consistent framework allowing comparing precipitation depth across regions. However, the standardization process imposes a common distribution to the transformed data and the same set of thresholds for drought classification. Therefore, each class of drought has the same probability of occurrence and the same intensity whatever the geographical location (Carrao et al 2014). More importantly, standardization results in a loss of information on the characteristics of precipitation deficit, especially on rare events in the distribution tails, which are precisely the most important to catch in risk analysis.

In this paper, we focus on extreme events, and propose an alternative method to drought hazard estimation based on extreme value theory. Compared to previous ones, our approach does not impose a unique set of thresholds for drought classification nor standardization, or any other transformation, of the precipitation data. The drought threshold and the probability of occurrence of a drought are estimated, location by location, using an extreme value mixture model with a normal distribution for the bulk, and a Generalized Pareto distribution (GPD) for the upper and lower tails of the rainfall distribution. The great advantage of this approach is to provide estimation of droughts' intensity and of their probability of occurrence that are specific to each geographical location. As a consequence, the drought hazard index can be measured more precisely, using the estimated intensity of each drought event instead of a somehow arbitrary system of weights and rating. Another advantage of our method is to provide an estimation of the drought threshold in the same unit than the raw data, allowing setting up an operational drought monitoring system.

To test the relevance of our approach to drought hazard estimation we work on a gridded set of rainfall data covering West Africa, from January 1901 to December 2016. Although the study area is limited, it includes a wide range of climatic zones. To facilitate comparisons with the standard approach of drought classification and DHI evaluation, we focus on meteorological droughts measured from cumulated annual rainfall. For ease of interpretation most of the results are provided in the form of maps.

The remainder of the paper is organized as follows. Section 2 presents the extreme value mixture model. Section 3 describes the studied area and the database. Section 4 provides the estimation results. Section 5 compares our estimation of drought hazard with that obtained using the standard approach. Section 6 concludes.

The extreme value mixture model

We focus on extreme precipitation events in the lower tail of the rainfall distribution. Drawing on the extreme value theory, we define the drought threshold as the precipitation value below which precipitation distribution can be approximated by an extreme value model. Symmetrically, we consider that heavy precipitation, above an upper threshold, can also be approximated by an extreme value model.

To estimate these thresholds separating the bulk of the distribution from the tails, mixture models are of great interest. Thanks to their ability to adapt to a wide variety of random phenomena, mixture models play an important role in statistical data analysis. They have been successfully used in a myriad of applications and statistical data engineering, for instance in the fields of econometrics as well as biology and epidemiology, meteorology, and more recently Machine learning. The general idea of the model is to express the underlying distribution as a weighted sum of a given family of reference distributions. In our context, we are interested only in finite mixture models i.e. finite weighted sums. Let us denote 𝑋 1 , … , 𝑋 𝑇 a random sample from a finite mixture of 𝑚 arbitrary distributions possibly 𝑑-vector valued; the density function 𝑓 of 𝑋 𝑘 can be defined, for all 𝑥 𝑖 ∈ ℝ 𝑑 , as follows:

𝑓 Θ (𝑥 𝑖 ) = ∑ 𝜋 𝑗 𝜑 𝑗 (𝑥 𝑖 ) 𝑚 𝑗=1 , 𝑤ℎ𝑒𝑟𝑒 ∑ 𝜋 𝑗 𝑚 𝑗=1 = 1.
Denote Θ = (Π, 𝝓) = (𝜋 1 , … , 𝜋 𝑚 , 𝜑 1 , … , 𝜑 𝑚 ) an object describing the mixture density of the sample. We assume that the densities 𝜑 𝑗 are drawn from some family ℱ of multivariate density functions. Mixtures are said to be parametric or semi-parametric when the elements of ℱ are completely described by a finite or a discrete number of scalars. Gaussian mixtures are, by far, the most popular and most used parametric mixture models in literature either in Gaussian regressions or in clustering problems. In this case:

ℱ ℊ = {𝜑(⋅ |𝜇, 𝜎 2 ) = density of 𝒩(𝜇, 𝜎 2 ), (𝜇, 𝜎 2 ) ∈ ℝ × ℝ + * }
where 𝒩(𝜇, 𝜎 2 ) is the normal distribution of mean 𝜇 ∈ ℝ and variance 𝜎 2 > 0, with probability density function defined for all 𝑥 ∈ ℝ by

𝜑(𝑥|𝜇, 𝜎 2 ) = 1 𝜎 𝑘 √2𝜋 exp (- 1 2 ( 𝑥 -𝜇 𝑘 𝜎 𝑘 ) 2 )
The general representation of the model reduces so that the parameters become

Θ = (Π, 𝝓) = (𝜋 1 , … , 𝜋 𝑚 , 𝜇 1 , 𝜎 1 , … , 𝜇 𝑚 , 𝜎 𝑚 ).
Despite their success and popularity, Gaussian mixtures are not convenient to model rare and extreme events that are intrinsically linked to precipitation data. In the aim to model these events we draw on the "Peak Over Threshold" (POT) approach of Pickands (1975). Under the conditions of the Pickands-Balkema-de Haan theorem, the asymptotic distribution of observations X exceeding a sufficiently high threshold is a Generalized Pareto distribution (GPD). Therefore, we adopt the generalized Pareto family to model the precipitation tails:

ℱ 𝓅 = {𝜑(⋅ |𝑢, 𝜎, 𝜉) = density of 𝐺𝑃𝐷(𝑢, 𝜎, 𝜉) for (𝑢, 𝜎, 𝜉) ∈ ℝ × ℝ + * × ℝ}
where 𝐺𝑃𝐷(𝑢, 𝜎, 𝜉) is the Generalized Pareto distribution with location parameter 𝑢 ∈ ℝ, scale parameter 𝜎 > 0 and shape parameter 𝜉 ∈ ℝ. For all the rest of the paper we will take parametric mixtures: ℱ = ℱ ℊ ∪ ℱ 𝓅 We consider a particular mixture model, denoted GNG, with m=3: a generalized Pareto distribution for the upper and lower tails, and a normal for the bulk of the distribution. The cumulative distribution function for the GNG mixture is given by 1 : [START_REF] Hu | evmix: Extreme Value Mixture Modelling, Threshold Estimation and Boundary Corrected Kernel Density Estimation[END_REF].

𝐹(𝑥) = {  𝑙 𝐺 𝑙 (𝑥|𝑢 𝑙 , 𝜎 𝑙 ,  𝑙 ) 𝑥 < 𝑢 𝑙  𝑙 + (1 - 𝑙 - 𝑢 ) [ 𝐻(𝑥) -𝐻(𝑢 𝑙 ) 𝐻(𝑢 𝑢 ) -𝐻(𝑢 𝑙 ) ] 𝑢 𝑙  𝑥  𝑢 𝑢 (1 - 𝑢 )+ 𝑢 𝐺 𝑢 (𝑥|𝑢 𝑢 , 𝜎 𝑢 ,  𝑢 ) 𝑥 > 𝑢 𝑢 1 See
H(x) is the normal cumulative distribution function, G(x) the conditional GPD cumulative distribution function. The parameters to be estimated are: the location parameters 𝑢 𝑙 and 𝑢 𝑢 , defining respectively the lower and upper tail thresholds;  l the lower tail GPD scale parameter; 𝜉 𝑙 the lower tail GPD shape parameter;  u the upper tail GPD scale parameter; 𝜉 𝑢 the upper tail GPD shape parameter. The tail fractions  𝑢 and  𝑙 represent, respectively, the probability of being above the upper threshold and the probability of being below the lower threshold. In our estimations they are given by:  𝑙 = 𝐻(𝑢 𝑙 ) and  𝑢 = 1 -𝐻(𝑢 𝑢 ) so that F(x) reduces to:

F(𝑥) = { 𝐻(𝑢 𝑙 )𝐺 𝑙 (𝑥) 𝑥 < 𝑢 𝑙 𝐻(𝑥) 𝑢 𝑙  𝑥  𝑢 𝑢 𝐻(𝑢 𝑢 )+(1 -𝐻(𝑢 𝑢 )) 𝐺 𝑢 (𝑥) 𝑥 > 𝑢 𝑢
The conditional GDP cumulative distribution function up to the lower threshold is given by:

𝐺 𝑙 (𝑥|𝑢 𝑙 , 𝜎 𝑙 ,  𝑙 ) = Pr(𝑋 < 𝑥|𝑋 < 𝑢 𝑙 ) = { ⌊1 + 𝜉 𝑙 ( 𝑢 𝑙 -𝑥 𝜎 𝑙 )⌋ -1 𝜉 𝑙 ⁄ if 𝜉 𝑙  0, exp [-( 𝑢 𝑙 -𝑥 𝜎 𝑙 )] if 𝜉 𝑙 = 0 for all 𝑥 ∈ ]𝑢 𝑙 + 𝜎 𝑙 𝜉 𝑙 ; 𝑢 𝑙 [ if 𝜉 𝑙 < 0 ; for all 𝑥 < 𝑢 𝑙 if 𝜉 𝑙  0
The conditional GDP cumulative distribution function above the upper threshold is given by:

𝐺 𝑢 (𝑥|𝑢, 𝜎 𝑢 ,  𝑢 ) = Pr(𝑋 < 𝑥|𝑋 > 𝑢 𝑢 ) = { 1 -⌊1 + 𝜉 𝑢 ( 𝑥 -𝑢 𝑢 𝜎 𝑢 )⌋ -1 𝜉 𝑢 ⁄ , if 𝜉 𝑢  0 1 -exp [-( 𝑥 -𝑢 𝑢 𝜎 𝑢 )] , if 𝜉 𝑢 = 0 for all 𝑥 ∈ ]𝑢 𝑢 ; 𝑢 𝑢 - 𝜎 𝜉 𝑢 [ if 𝜉 𝑢 < 0 ; for all 𝑥 > 𝑢 𝑢 if 𝜉 𝑢  0.
The shape parameter () determines the tail behavior: if  = 0, the distribution has an exponential, light tailed, distribution and belongs to the Gumbel family; if  < 0 the distribution has short bounded tail and belongs to the Weibull family; if  > 0 the distribution has a heavier tail and belongs to the Fréchet family. As an illustration, Fig. 1 depicts the density of a Gaussian law and of a mixture with two Weibull distributions for the tails, and a normal for the bulk.

Figure 2 shows the sensitivity of the left threshold to the shape parameter of the GPD: the threshold is close to the mean when 𝜉 𝑢 is positive, and moves away from it when 𝜉 𝑢 is negative.

.0000

.

0004

.0008

.0012

.0016

.0020

Gaussian GNG Density

Left threshold

Right threshold

Fig. 1 Results of simulations: probability density of a GPD -Normal -GPD mixture (red) and of Gaussian distribution (blue). The two distributions have the same mean.

.0000

.

0002

.0004

.0006

.0008

.0010

.0012

.0014

.0016

1,000

Region and Data

The region under study is located between latitudes 0 and 18°N, and longitude 20°W and 20°E (Fig. 3). This area is centered on West Africa but includes also a part of Central Africa. Rainfall patterns in this region are determined by the monsoon system and topography. The region encompasses three major climatic areas and several climate transition zones, which cross the region from east to west, forming virtually parallel strips from the Tropic of Cancer to the Equator.

In the north of the studied area, above latitude 18°N, the climate is semi-desert, of northern Sahelian type, or desert. The rainy season is concentrated over 1 or 2 months, average annual rainfall varies from 300 mm to less than 100 mm. Further south, the tropical climate is characterized by very distinct dry and wet seasons. The length of the wet season and the amount of annual rainfall allows distinguishing 3 types of tropical climates. The Sahelian climate, or semi-arid tropical, lies between latitudes 12°N and 18°N; the dry season lasts 7 to 10 months, and the average annual rainfall varies between 300 and 700 mm. The Sudanian climate, or pure tropical, corresponds to a dry season of 5 to 6 months, with an average annual rainfall between 700 and 1000 mm (1200 mm on the Senegalese coast). The Guinean climate, or transitional tropical, corresponds to a shorter dry season, from 4 to 5 months, and an average annual rainfall above 1000 mm. The transitional equatorial and equatorial climate lie below 10°N over the countries of the Gulf of Guinea. Annual precipitations are abundant, generally above 1000 mm with a bimodal distribution (two dry seasons). However, the rainfall distribution is quite heterogeneous in this area. The climate is very humid, with annual rainfall above 2500 mm, in the coastal regions of Guinea, Cameroon and Gabon, which are characterized by mountainous relief. By contrast, the Dahomey Gap, which includes parts of Ghana, Togo and Benin, with annual rainfall total of 1000 to 1200 mm, is abnormally dry for the region (see Fig. 3). The rainfall data comes from the Climatic Research Unit at the University of East Anglia. We use the gridded time-series of monthly rainfall total at 0.5° × 0.5° latitude/longitude resolution, for the period January 1901 -December 2016. The data was extracted from the revised version of January 2016 of the dataset referred to as CRU TS3.10. The construction of the CRU dataset is detailed in [START_REF] Harris | Updated high-resolution grids of monthly climatic observations -the CRU TS3.10 Dataset[END_REF]. Before proceeding to the estimations, the database had to be cleaned. Prior to 1930, numerous outliers have been detected, especially for cells located in the northern part of the region corresponding to the Sahara border. Therefore, the area above the latitude 18°N has been excluded from the sample. We also dropped from the dataset all total annual rainfall that repeat during two or more consecutive years.

We computed annual rainfall from monthly data. Annual precipitation total is commonly used to assess meteorological drought. A decline in annual total is also a rough indicator of agricultural drought especially in the Sahelian and tropical part of the region where the rainfall season lasts 2 to 6 months and the annual rainfall total is equivalent to the seasonal cumulated rainfall.

Estimation results

The mixture model is estimated cell by cell, for annual rainfall total over the period 1901-2016. As a rule, we present only the estimation results for the left tail of the distribution, corresponding to drought events.

Goodness of fit

According to the Kolmogorov-Smirnov (KS) goodness of fit test, the estimated mixture model provides a satisfactory representation of the data for the vast majority of cells (Fig. 4). The KS test rejects the equality of the empirical and the mixture cumulative distribution function (CDF) in 3.5% of cases only. This good result is confirmed by the chi-square test, which rejects the null hypothesis of no significant difference between the observed and expected value in 1.7% of cases only. Results are not presented, but are available upon request.

We also compared the accuracy of the GNG mixture and the gamma distribution, which is used to compute the SPI, on the basis of the Root Mean Square Error (RMSE). The difference in the RMSE of the GNG and the gamma shows that the GNG mixture model better fits the data in 62.9 % of cases. To formally test whether the prediction errors from the GNG and the Gamma are significantly different, we run the test of Diebold and Mariano (1992) corrected for small sample-bias [START_REF] Harvey | Tests for Forecast Encompassing[END_REF]. Using the absolute error loss criterion, the test results provided in Fig. 5, show that GNG mixture has better predictive accuracy than the Gamma in 39.4 % of cases. Conversely, the Gamma distribution has better predictive in 23.2 % of cases only. We can conclude that the accuracy of the GNG model is as high as or higher than that of Gamma in 83% of locations. When using the square error loss criterion, the results of the Diebold and Mariano test are not as good but still favorable the GNG model, which accuracy is as high as or higher than that of Gamma in 77% of locations. Brown cells: lower accuracy of the GNG is rejected and higher accuracy of the GNG is not rejected. Orange cells: higher accuracy of the GNG is rejected and lower accuracy of the GNG is not rejected. White cells: the forecasts accuracy of the GNG and of the Gamma are not significantly different.

Spatial pattern of the GPD parameters

The location parameter of the lower GPD, defining the extreme drought threshold, displays the same latitudinal trend than annual precipitations (Fig. 6 and7). To allow spatial comparisons, the location parameters are expressed as z-score:

(𝑥-𝜇) 𝜎 with  and  the mean and standard deviation of annual rainfall over the 1981 -2010 period. Therefore, the standardized values indicate the intensity of the precipitation deficit below which a dry spell can be classified as an extreme drought: the lowest the standardized threshold value, the higher is the intensity of extreme droughts (Fig. 8).

A break appears in the latitudinal organization of the standardized thresholds: they tend to be higher above the 10° latitude and lower below (Fig. 9). Three main clusters of particularly low values appear whose epicenter is located, respectively, in the western part of Nigeria, the southwest and southeast part of the Central African Republic (Fig 8).

The latitudinal trend is less obvious for the scale parameter ( 𝑙 ) but there is some evidence of a decreasing trend in sigma above the latitude 12 (Fig. 12 and13). Below, the spatial distribution of sigma is quite heterogenous. The shape parameter ( 𝑙 ) displays the most erratic spatial pattern (Fig. 14 and15). In more than 77% of cases, the estimated shape parameter is negative, meaning that the precipitation distribution has a short bounded left tail. The mixture estimation also gives an implicit parameter, the probability of being below the location parameter or tail fraction ( 𝑢 ) (Fig. 10). Results show that the probability of nonexceedance the extreme drought threshold is concentrated around its mean (11.4%), with no clear spatial pattern (Fig. 11).

The Drought hazard and hazard exposure maps

Focusing on extreme droughts we consider that rainfall totals below the extreme value threshold, given by the location parameter of the left GDP, constitute a homogenous category and do not distinguish different classes of droughts below this threshold. Therefore, the drought hazard index (𝐷𝐻𝐼 𝑔𝑛𝑔 ) is computed as the product of the minimum intensity of extreme droughts, given by the standardized left threshold value (𝑍 𝑢𝑙 ), and the probability of being below the threshold ( 𝑢𝑙 ):

𝐷𝐻𝐼 𝑔𝑛𝑔 = |𝑍 𝑢𝑙 |.  𝑢𝑙 (1)
Figure 16 displays the DHI, which has been rescaled using the min-max normalization method.

Figure 17 displays the Gi* statistic of [START_REF] Getis | The Analysis of Spatial Association by Use of Distance Statistics[END_REF]. The G* statistic evidences four clusters of large values of the DHI (hot spots), higher than average. The main one is located in a region covering southwestern Nigeria, southern Benin and Togo. Two other hot spots of drought hazard are located, respectively, in the border area of Cameroon and the Republic of Central Africa (RCA), and in the southeast part of RCA. Last, another spatial area of high values of DHI emerges in the northern edge of the study area, straddling Niger and Chad. By contrast, numerous cold spots, of DHI values lower than average appear. Most of them are located above the 10 th parallel.

Figure 18 highlights a U-shape relationship between DHI and latitude: DHI is the lowest between parallel 10 and 15, and higher above and below this band. Clusters of high and low values of DHI alternate from east to west, giving a sinusoidal shape to the relationship between the DHI and the longitude (Fig. 19).

Comparison with the standard approach of drought hazard

To assess the relevancy of our approach to drought hazard estimation, we compare our estimation of the DHI with that obtain following the standard approach, and compare the ability of the GNG mixture and the SPI to model extreme droughts.

The drought hazard index

Drought hazard is commonly defined as a combination of the intensity of droughts and their frequency, and estimated from the SPI and the drought classification of McKee et al (1993) (Table 1). Droughts are classified according to a fixed set of thresholds, and the probability of occurrence of a drought of a given intensity is the same whatever the geographical location and the statistical characteristics of the raw data. Therefore, the intensity of droughts is caught through a fixed weighting, and the DHI is given by:

𝐷𝐻𝐼 𝑠𝑝𝑖 = ∑ 𝑊 𝑖 × 𝑓 𝑖 4 𝑖=1 (2)
𝑊 𝑖 is the weight assigned to droughts of class 𝑖 (Table 1, col.2); 𝑓 𝑖 is the frequency of droughts of type 𝑖, namely the number of drought events of type 𝑖 recorded during the period over the number of observations in the period (table 1, col.5). The 𝐷𝐻𝐼 𝑠𝑝𝑖 map is given by Fig. 20. and7). We measure the intensity of a class of drought by the standardized value of the upper threshold of each class. Thus, the DHI is given by:

Drought
𝐷𝐻𝐼 𝑔𝑛𝑔 = ∑ |𝑧 𝑖 | × 𝑓 𝑖 3 𝑖=1
(3)

With 𝑧 1 : the standardized value of the 16 th percentile; 𝑧 2 : the standardized value of the 7 th percentile; 𝑧 3 : the standardized value of the 2 nd percentile of the GNG distribution. The corresponding map is given by Fig. 21. The two maps, Fig. 20 and Fig. 21 give a very different picture of the spatial distribution of the drought hazard. According to the standard approach (Fig. 20), the draught hazard is higher in the Sahelian band and on the western side of the studied area: Senegal, Guinea Bissau, Guinea, Burkina Faso, Northern Nigeria, and Southern Chad. By contrast, according to our approach, the drought hazard is higher in the southern quadrant of the region (Fig. 21).

These differences are mainly due to the standardization process. In the standard approach, the thresholds, and consequently the intensity of the various drought categories (Wi), are the same for all series. In our approach, the drought classification is based on percentiles of the theoretical distribution of the raw data. The thresholds and consequently the intensity of the various drought categories, z i , vary according to the characteristics of each rainfall distribution.

Another source of divergence, of lesser importance, comes from the unequal quality of adjustment of the data. While the empirical frequency of the different types of droughts (fi) should the same in the two approaches, and closed to their theoretical probability of occurrence, some discrepancies appear (Table 1, col 5 and 8). Unsurprisingly, the GNG model performs better to predict extreme drought. This point is explored below.

Distribution of precipitation in the lower tail

The drought classification in We can see that the GNG mixture performs better than the SPI to model the left tail of the rainfall distribution, with 29.2% of observations below the 30 th percentile compared to 28.1% for the SPI. When looking at the distribution of observations below the 30 th percentile, the GNG is better at modeling the highest and lowest percentiles corresponding, on the one hand, to the "abnormally dry" and "moderate drought" categories and, on the other hand, to "exceptional drought". Both the SPI and the GNG perform poorly in detecting "severe droughts" between the 10 th and 5 th percentiles.

The centennial drought

The centennial drought has a return period of 100 years and a return level equals to the first percentile of the distribution. Therefore, the 100 years return level is the amount of annual rainfall total that can be expected on average in 1 year out of 100 years. Considering these centennial droughts, the GNG mixture model appears more efficient than the standard approach based on the SPI to model these events.

The number of dry spells with an annual rainfall total less than or equal to that of the 100-year drought, differs greatly depending on whether it is calculated using the GNG or the SPI distribution (Table 3). According to the SPI approach, almost 12% of locations did not experience a rainfall episode equal to or below the 100-year drought level while 19% of locations experienced more than 2 drought years below the 100-years level, in the period 1901 -2016. Results are more consistent when using the GNG distribution: less than 1% of locations did not experienced a drought spell below the 100-years drought level, and only 2.7% of locations experienced more than 2 drought episodes below the 100-years level. Figure 22 shows the intensity of the centennial drought calculated from the GNG distribution, and measured by the standardized return value. The comparison with the intensity of the centennial drought calculated from the gamma distribution, shows that in a majority of cases (63.4%), the Gamma minimizes the severity of the centennial drought (Fig. 23). In other words, if the GNG is a good approximation of the rainfall distribution, using the Gamma, will lead to underestimate the severity of the centennial drought.

Conclusion

To design efficient drought risk management and risk coping strategies, at the national or regional level, we need above all information on the temporal and spatial distribution of precipitation. This information is necessary to assess drought hazard and hazard prone locations. A critical stage in drought hazard assessment is the definition of a drought event, and the classification of droughts. While the standard method is based on an ad hoc fixed set of thresholds, the method we propose aims at providing an objective definition and a measure of drought, based on the statistical properties of the rainfall distribution.

We consider the class of droughts that are potentially damaging, which distribution can be approximated by an extreme value model. The extreme drought threshold is not set exogenously but estimated using the mixture model with a generalized Pareto distribution for the upper and lower tails, and a normal for the center of the distribution. The extreme drought threshold is given by the location parameter of the lower tail GPD and the probability of precipitation below the threshold is given by the lower tail fraction. Therefore, the two components of drought hazard, the drought threshold and the drought probability, are idiosyncratic, dependent on the specific characteristics of precipitation data.

Another advantage of the proposed method is its operational aspect. The drought threshold can be easily estimated and the estimation accuracy tested. The drought threshold is measured in the same unit than the rainfall data, in mm per year in our case, so that a rainfall deficit can be readily classified and appropriate responses quickly triggered.

The application conducted on the West Africa region and an indicator of meteorological drought, shows the great sensitivity of results to the method used to estimate the drought hazard. The use of the SPI to identify extreme drought leads to underestimate the intensity of extreme droughts, and to an abnormally high frequency of extreme events with a return period of 100 years or more. This result potentially challenges many research studies that use the SPI to identify rainfall shocks and estimate their economic and social impact.

Last the map of drought hazard index offer a different picture of the spatial distribution of the hazard prone locations. Our drought hazard measure shows more heterogeneity across locations than the standard one. More important, our results show that the climatic regions more favorable to agricultural production are also the more exposed to extreme droughts. Considering that the annual rainfall is a rough indicator for agricultural drought, the results should lead to pay a greater attention to the wetlands, which are the most agriculture intensive and where the intensity of extreme droughts is the highest.

This analysis should be considered as exploratory. Its scope could be extended to a larger geographic area, and to other indicators of drought. In drought hazard assessment, the assumptions made on the rainfall distribution are crucial and the relevancy of the drought hazard estimates depends on the goodness of fit of the assumed distribution. The tests and the distribution of drought events show that the mixture model better fit the data than the standardization process used to calculate the SPI. However, a still better fitting model can be looked for in further research.

Fig. 2

 2 Fig. 2 Drought threshold for different values of the shape parameters of the left tail GPD ( 𝑙 )

Fig. 3

 3 Fig.3Average annual rainfall over the period1901-2016 (mm)

Fig. 4

 4 Fig.4 Kolmogorov-Smirnov test. Red cells: bootstrap p-value < 0.05, the GNG does not fit the data.

Fig. 5

 5 Fig.5 Diebolt and Mariano Test. Loss function: Absolute Error

Fig. 6 Fig. 8

 68 Fig. 6 Extreme drought threshold value in mm (𝒖 𝒍 ) Fig. 7 Extreme drought threshold by latitude

Fig. 16

 16 Fig. 16 Rescaled DHI calculated from equation (2).

Fig. 17

 17 Fig. 17DHI 1 . Getis Ord statistics. Red cells: hot spot (high values) of DHI; yellow cells: cold spot (low values)

Fig. 20

 20 Fig. 20 Rescaled DHI calculated using the SPI and equation (2)

Fig. 22

 22 Fig. 22 GNG 100-year drought level: standardized values

Table 1 .

 1 Drought type designation and thresholds. DHI. W: weight For comparison purposes, we recalculate our DHI measure, based on the GNG mixture, distinguishing four categories of drought intensity with occurrence probabilities closed to that of McKee (Table 1, col. 6

	category	W	Threshold	prob. of	Empirical	Thresholds	prob. of	Empirical
			values for the	occurrence	frequency:	for the GNG	occurrence	frequency:
			SPI		SPI	in percentiles		GNG
	Mild drought	0	] -1 ; 0 ]	34 %	0.339	]16% ; 50% ]	34%	0.354
	Moderate drought 1	] -1.5; -1 ]	9.2%	0.083	]7% ; 16% ]	9%	0.088
	Severe drought	2	] -2; -1.5 ]	4.4%	0.040	] 2%; 7% ]	5%	0.043
	Extreme drought	3	 -2	2.3%	0.029	 2%	2%	0.020

Table 2 .

 2 Table 2, col.1, is that commonly used by the US Drought Monitor (USDM)(Sovboda et al 2002). Thresholds are given by percentiles of the rainfall distribution (col. 2). Column.3 gives the corresponding SPI values i.e. the percentiles of the Gaussian law. Column 4 gives the theoretical probability of occurrence of each type of drought while columns 5 and 6 give the mean frequency of each drought type using alternatively the SPI values or the GNG percentiles as thresholds. Probability of occurrence and empirical frequency of drought categories. Empirical frequency is the mean frequency over the whole sample, which includes 2981 locations and 263754 annual rainfall observations.

	Drought categories	Percentiles	SPI values	Theoretical	Empirical	Empirical
	(USDM)			prob. of	frequency:	frequency:
				occurrence	SPI	GNG
	Abnormally dry	]20% ; 30%]	] -0.84; -0.53]	10%	0.093	0.105
	Moderate drought	]10% ; 20%]	] -1.29; -0.84]	10%	0.091	0.098
	Severe drought	]5% ; 10%]	] -1.65; -1.29]	5%	0.044	0.044
	Extreme drought	]2% ; 5%]	] -2.05; -1.65]	3%	0.027	0.025
	Exceptional drought	2%	-2.05	2%	0.026	0.020
	Total			30%	0.281	0.292

Table 3 .

 3 Frequency of 100-years drought episodes in the period 1901 -2016

Acknowledgment

The authors acknowledge the support received from the Agence Nationale de la Recherche of the French government though the IDEX-ISITE initiative (16-IDEX-0001 CAP 20-25) and the LABEX IDGM+ (ANR-10-LABX-14-01), within the program ''Investissements d'Avenir''.