# Global assessment of storm disaster-prone areas

Nazzareno Diodato, Pasquale Borrelli, Panos Panagos, Gianni Bellocchi

► To cite this version:

Nazzareno Diodato, Pasquale Borrelli, Panos Panagos, Gianni Bellocchi. Global assessment of storm disaster-prone areas. PLoS ONE, 2022, 17 (8), pp.e0272161. 10.1371/journal.pone.0272161. hal-04160823

HAL Id: hal-04160823

https://hal.inrae.fr/hal-04160823

Submitted on 12 Jul 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

RESEARCH ARTICLE

# Global assessment of storm disaster-prone areas

Nazzareno Diodato¹, Pasquale Borrelli²,³, Panos Panagos¹,⁴*, Gianni Bellocchi¹,⁵

1 Met European Research Observatory—International Affiliates Program of the University Corporation for Atmospheric Research, Benevento, Italy, 2 Department of Science, Roma Tre University, Rome, Italy,

3 Department of Biological Environment, Kangwon National University, Chuncheon, Republic of Korea,

4 European Commission, Joint Research Centre (JRC), Ispra, Italy, 5 Université Clermont Auvergne, INRAE, VetAgro Sup, UREP, Clermont-Ferrand, France

* Panos.PANAGOS@ec.europa.eu

## Abstract

### Background

Advances in climate change research contribute to improved forecasts of hydrological extremes with potentially severe impacts on human societies and natural landscapes. Rainfall erosivity density (RED), i.e. rainfall erosivity (MJ mm hm⁻² h⁻¹ yr⁻¹) per rainfall unit (mm), is a measure of rainstorm aggressiveness and a proxy indicator of damaging hydrological events.

### Methods and findings

Here, using downscaled RED data from 3,625 raingauges worldwide and log-normal ordinary kriging with probability mapping, we identify damaging hydrological hazard-prone areas that exceed warning and alert thresholds (1.5 and 3.0 MJ hm⁻² h⁻¹, respectively). Applying exceedance probabilities in a geographical information system shows that, under current climate conditions, hazard-prone areas exceeding a 50% probability cover ~31% and ~19% of the world's land at warning and alert states, respectively.

### Conclusion

RED is identified as a key driver behind the spatial growth of environmental disruption worldwide (with tropical Latin America, South Africa, India and the Indian Archipelago most affected).

## Introduction

Although there is a growing need to assess ecosystem responses to climate change-induced disaster risk reduction, there is a lack of research on sensitive areas and on coastal zones, drylands and watersheds, particularly in Global South low-income countries [1, 2]. In order to advance global climate change studies, climate impact indicators are needed to help develop

<figure><img src="image_3.png" /></figure>

**OPEN ACCESS**

**Citation:** Diodato N, Borrelli P, Panagos P, Bellocchi G (2022) Global assessment of storm disaster-prone areas. PLoS ONE 17(8): e0272161.
https://doi.org/10.1371/journal.pone.0272161

**Editor:** Chun Liu, Jinan University, CHINA

**Received:** March 4, 2022

**Accepted:** July 13, 2022

**Published:** August 24, 2022

**Peer Review History:** PLOS recognizes the benefits of transparency in the peer review process; therefore, we enable the publication of all of the content of peer review and author responses alongside final, published articles. The editorial history of this article is available here: https://doi.org/10.1371/journal.pone.0272161

**Copyright:** © 2022 Diodato et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data Availability Statement:** The data can be requested using the request form through the following link: https://esdac.jrc.ec.europa.eu/content/global-rainfall-erosivity.

**Funding:** The authors received no specific funding for this work.
**Competing interests:** The authors have declared that no competing interests exist.

guidelines for landscape conservation planning [3, 4] and support decision-making based on hydrological damage data [5, 6]. Rainfall erosivity density (RED in MJ hm⁻² h⁻¹, equivalent to the most common MJ ha⁻¹ h⁻¹), i.e. rainfall erosivity per unit of rainfall, is an important climatic indicator of floods and soil erosion [7, 8]. RED effects are apparent in catastrophic weather events, due to the intensification of daily storms in the northern Hemisphere (Fig 1A) and erosive precipitation in other parts of the world (Fig 1B). The effects of more intense extreme rainfall worldwide [9–12] include intensified sub-daily precipitations [13, 14], the local occurrence of flash floods [15] and the incidence of the erosive force of rainfall [16–19].

Damaging hydrological events are extreme phenomena, the source of multiple hazardous events with potentially serious impacts on human societies [21]. Among these, floods and events causing flooding are identified as hydrological disasters, which also include meteorological disasters like thunderstorms. Like other extreme phenomena, hydrological disasters typically leave behind socio-economic damage, the severity of which depends on the resilience of the affected population and the available infrastructures [22]. Potentially severe natural events are typically not classified as natural disasters if they occur in areas without a vulnerable population, e.g. deserts [23]. Information on the spatial distribution of RED can help delineate areas prone to multiple damaging hydrological events [24, 25] and, in turn, can transfer important prescriptions for disaster response planning. However, geoinformation and hazard mapping have been and remain critical in modern visual communication science practice [26], particularly in landscape decision-making, where the uncertainty component must be included in the mapping [27, 28]. In particular, it is difficult to represent disaster-prone areas by RED, as different parts of terrestrial ecosystems respond differently to uneven, often nonlinear and universal forcing agents, with threshold-like features [29, 30]. Landscapes are, in fact, highly responsive and non-linear systems to both external dynamics, such as climatic and non-climatic factors, with a combination of gradual changes coupled with infrequent high-magnitude events has led to dramatic landscape responses throughout Earth's history [31]. Then, the effects of these historical landscape responses to climate extremes and threshold processes are key parameters affecting the geomorphological impacts of extreme hydrological events [32, 33].

Threshold mapping, and the processes of developing probability mapping, is a challenge in geographic information science and spatial downscaling [34]. Scaling and integrating the relative uncertainty of these thresholds, for which storms drive surface flows, including

<figure><img src="image_4.png" /><figcaption>**Fig 1. Emerging catastrophic weather events over the last four decades.** a) Annual mean daily maximum storm depth (mm d⁻¹) in the Northern Hemisphere (mean data of areal maxima from NCEP/NCAR Reanalysis [20]), and b) Global evolution of the erosive force of rainfall (rainfall erosivity) in MJ mm hm⁻² h⁻¹ yr⁻¹ (arranged from Bezak et al. [19]). Graph a) covers the period 1981–2018, and graph b) covers the period 1981–2020.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0272161.g001
geomorphic effectiveness and erosional soil degradation, affect the spatial pattern of damaging hydrological events, and ultimately improve our understanding of ecological responses to climate extremes and thresholds [35].

The apparent naturalness of a landscape measures the degree to which it is free from storm damage (Fig 2A). Under erosive climatic conditions, the strength, intensity and frequency of a given rainfall influence hydrological processes, even though they may still maintain the shape of the land and the equilibrium of the environmental system in a resilient landscape (Fig 2B). On the contrary, a change in the hydrological regime, especially when the thresholds of an acceptable level of hydrological disturbance are exceeded, can have harmful consequences for the landscape (Fig 2C).

The prerequisite for RED modelling is that the environmental system adapts to changes in the natural hydrological regime. Usually, the prediction of precipitation and its extremes is performed using physical models, mainly due to the high spatial variability and nonlinearity of the problem [36]. The limitations of physical models (which are mainly computational at the global level) are encountered with very extreme rainfall, when the prediction has to be performed at time-scales less than 30 minutes over many years and for a spatial resolution of a few kilometres. Recently improved convection models can predict extreme rainfall, but simulations are not currently within reach due to their computational cost and degree of uncertainty [37]. Retrieving extreme rainfall from satellite data would be an option to improve the estimation of rainfall erosivity [38], but satellite data require corrections and need further assessment, especially on a global scale [39]. With the advancement of global atmospheric reanalysis data, numerical weather prediction techniques are becoming an encouraging mean of estimating rainfall erosivity [40]. In particular, the production of global RED maps is a challenge due to two conflicting conditions, namely that the analysis requires precipitation data with both high temporal resolution and global coverage. With a large volume of data and uneven distribution of stations, geostatistical methods provide reasonable estimates of what the variable of interest would be at intermediate locations. Geostatistics offers different approaches to deal with this issue and provides attracting results when experimentally determined rainfall erosivity data are available at both regional [41, 42], continental [43, 44], or even global scales [5]. However, geostatistical estimation of RED [17, 45, 46], and its spatial patterns above given thresholds [47], have generally received limited attention [48].

<figure><img src="image_5.png" /><figcaption>**Fig 2.** Landscape under different weather regimes with changing RED thresholds for a given level of hydrological disturbance in: a) Storm-free landscape, b) Landscape under storms not exceeding thresholds, and c) Stressed landscape where thresholds are exceeded resulting in damaging hydrological events (image arranged from MeteoBlue, https://static.meteoblue.com/assets/images/crosslinks/yearcomparison.svg).</figcaption></figure>

https://doi.org/10.1371/journal.pone.0272161.g002
The value of geostatistical probability mapping for geographic information systems (GIS) lies in providing reliable interpolation methods in error assessment and scaling, which can be used in environmental modelling and decision-making [49]. The spatial uncertainty associated with RED hazard over a range of scales is, in any case, an open issue, coupled with uncertainties in downscaling methods and a lack of primary information data in many areas. The uncertainty of RED thresholds actually poses challenges for geospatial assessment, as the worst storms can fall at locations not well-covered by weather recording stations [50]. Global quantification of storms at fine time scales also remains challenging, as hourly or sub-hourly rainfall data of sufficient length are poor, especially for critical and vulnerable regions such as the tropics [51].

Here, we present for the first time the use of the Global Rainfall Erosivity Database (GLOREDa) to estimate RED worldwide [5]. For this study, this database was updated as GLOREDa-V2 and contains erosivity values estimated as RED from 3,625 stations in 63 countries with time resolutions from 1 to 60 minutes. The results obtained are based on climate data and a probabilistic approach to proceed under a soft geovisualisation in order to mitigate the uncertainty involved in downscaling and geocomputational tracking. For this purpose, we used a parametric kriging technique, which provides great flexibility for modelling environmental data [52]. In particular, GIS and log-normal ordinary kriging coupled with output probability mapping (LNOKpm) were used to continuously delineate the spatial uncertainty of RED thresholds and predict areas prone to damaging hydrological events on a global scale.

# Methods

## Rainfall erosivity data

We refer to annual rainfall erosivity data from the Global Rainfall Erosivity Database (GLOREDa, here updated to GLOREDa-V2), which covers 3,625 precipitation stations from 63 countries with temporal resolutions of 1 to 60 minutes. It is the result of an extensive collection of high temporal resolution rainfall data from as many countries as possible in order to have a representative sample across climatic and geographical gradients (Fig 3).

The number of stations varies greatly from continent to continent, with no station-data available above 70° North and below 47° South. However, the latter is not an issue. In fact, the heaviest rains have a low probability of occurring over the northern limit, where rainfall erosivity values are close to zero, while large areas below the southern limit are open sea water or ice-covered (Antarctic ice cap). Precipitation time-series range from a minimum of 5 years to maximum of 52 years (on average, 16.8 years).

GLOREDa contains the best available global datasets of annual rainfall erosivity (RE, MJ mm hm⁻² h⁻¹ yr⁻¹), in the form of (R)USLE-R factor [53], calculated on a monthly (j = 1, ..., 12) basis (RE_m, MJ mm hm⁻² h⁻¹ month⁻¹), from which we obtained rainfall erosivity density (RED, MJ hm⁻² h⁻¹) as the ratio between rainfall erosivity and precipitation amount [54]:

$$ RE_j = \frac{1}{n} \cdot \sum_{i=1}^{n} \sum_{k=1}^{m_j} \{ I_{30} \cdot \sum_{r=1}^{m} [0.29 \cdot (1 - 0.72 \cdot e^{-0.05 \cdot i_r})] \cdot v_r \} \cdot k \quad (1) $$

$$ RE_y = \sum_{j=1}^{j=12} RE_j \quad (2) $$

$$ RED = \frac{RE}{\frac{1}{n} \cdot \sum_{i=1}^{n} P_y} \quad (3) $$

where n is the number of years recorded; m_j is the number of erosive events during a given
<figure><img src="image_6.png" /><figcaption>Fig 3. Geographical extent of the global network of stations (red circles) where rainfall erosivity density (RED) data were available.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0272161.g003

month *j*; *k* is the index of a *k*th single event; *v*<sub>r</sub> is the volume of rainfall (mm) during the *r*th period of a storm, which splits into *m* parts; *I*<sub>30</sub> is the maximum 30-minute rainfall intensity (mm h<sup>-1</sup>); *i*<sub>r</sub> is the rainfall intensity during the time interval (mm h<sup>-1</sup>); *P*<sub>y</sub> is the amount of rain (mm) during a given year *y*. The numerator in Eq (3) represents the long-term mean rainfall erosivity (MJ mm hm<sup>-2</sup> h<sup>-1</sup> yr<sup>-1</sup>), while the denominator is the relative mean amount of precipitation (mm) over the same number of years.

## Warning and alert thresholds

We use concepts of a “warning state” and an “alert state” [55] for erosivity density, as regions with high RED are at risk of flooding and water scarcity owing to their infrequent but very intense and erosive rainstorm [19, 56]. As such, RED reflects not only the component of climate forcing reproduced in the aggressiveness of the storm (rainfall erosivity), but also the damaging hazard and its associated hydrological risk [18]. In particular, RED values > 3.0 MJ ha<sup>-1</sup> h<sup>-1</sup> indicate an increased risk of erosive rainstorms, soil erosion and flooding [57]. This critical threshold value was popularised in storm geomorphology by Dabney et al. [58] as the runoff increases when RED exceeds 3 MJ hm<sup>-2</sup> h<sup>-1</sup>, thus leading to an increasing erosive hazard as storm erosivity represents a large proportion of the rainfall amount in an intense event. Mostly set regionally on a monthly basis [44], this hydrological threshold (alert) could help detecting areas prone to erosion- and overland flow [59, 60]. Abstractions such as thresholds of change and strength are all essential as landscapes may be able to counter or assimilate pulses of change as a form of sensitivity or stability [61]. The use of this abstraction for the way geomorphological systems react to climate variability and change is still an important topic since its conception by Allison and Thomas [62] and Phillips [63]. Given the global scale and annual resolution of this study, a second, smaller RED threshold value (warning) was also used to better capture interactions between changing hydrologic variability and ecological responses. For instance, major tropical forest regions like Amazon and African regions are less vulnerable to hydrological hazards than Southeast Asia [64]. When RED values are > 1 MJ ha<sup>-1</sup> h<sup>-1</sup>, only a certain amount of precipitation can cause relatively high erosivity [57]. Therefore,
we adopted two threshold values as follows: RED exceeds the threshold of 1.5 MJ hm⁻² h⁻¹ in the warning state and exceeds 3.0 MJ hm⁻² h⁻¹ in the alert state. The two values correspond, respectively, to the median and the 3rd quartile of the distribution of the RED data (Fig 4A). These values also mark the range of critical RED values identified by Diodato et al. [8] with a 50-year return period in north-western Italy. The advantage of using RED threshold values compared to other rainfall aggressiveness indicators are [54]: a) greater propensity to classify geomorphological hazards; b) greater stability, obtained with a shorter and heterogeneous period of recording; and c) easier mapping over large areas due to its independence on altitude up to 3000 m a.s.l. [65].

## Log-normal ordinary kriging probability map

Kriging is a generic name covering a range of spatial least-square prediction methods [66]. Most kriging algorithms with references to hydrogeological applications were reviewed in Kitanidis [67]. For some kriging models, only one estimate per cell is required; for others, as in decision-making, it is necessary to know the local uncertainty associated with the estimates [68, 69]. The ordinary lognormal kriging algorithm has the potential to improve maps because such techniques can be combined with data transformation and detrending options, and can take different forms, such as maps of probability outputs, quantiles and standard errors of prediction, when the normality of the distribution is verified [52]. In contrast, indicator kriging cannot perform data-transform and detrending, when *soft* information is required within the probability mapping.

If the prediction in the unknown locations is normally distributed, then the mean and median of the predictions will be positioned centrally in the probability density distribution of each location.

The area under the distribution curve to the right of the threshold line predicts the probability that the value is greater than a threshold value. The distribution of predictions changes for each location as the mean and standard error of the predictions change. Thus, by keeping the threshold value constant, a probability map is produced for the whole area. Since the RED data show a skewed distribution, we used a lognormal ordinary kriging in the form of a probability output map (LNOKpm). The assumption of normality of the distribution was checked after the data were log-transformed. A straightforward approach is to classify as hazardous all locations where the probability of exceeding a critical threshold value, z_k, is greater than a critical RED value (1.5 MJ hm⁻² h⁻¹ for the warning state and 3.0 MJ hm⁻² h⁻¹ for the alert state). The ordinary kriging model assumes that the data are a realisation of an auto-correlated process plus an indipendent random error. For a complete analytical procedure of LNOKpm, refer to Krivoruchko et al. [70].

## Results and discussion

### Exploratory data analysis and transformation

Exploratory data analysis is important to inspect and explore data statistically before deciding whether and how to transform them for analysis and to illustrate what can be achieved by transforming data (e.g. into logarithms) of single variates and calculating principal component analysis of multivariate data [71]. The first step in the spatial analysis is to check the raw data for drifts and outliers [72], and finally for normal distribution. Outliers and spatially-skewed data can be detected by the frequency distribution and the third standardised moment or skewness (g), whose value is 0 for a normal distribution and any other symmetric distribution with finite third moment.
<figure><img src="image_7.png" /><figcaption>**Fig 4.** Exploratory data analysis and verification of normality of rainfall erosivity density (RED) data. a) Distributional frequency on the original RED data and b) after log-normal transformation; c) QQ-Plot of the theoretical and estimated distribution of the original RED data and d) after log-normal transformation. Warning and alert thresholds (MJ hm⁻² h⁻¹) are shown in a).</figcaption></figure>

https://doi.org/10.1371/journal.pone.0272161.g004

A lognormal transformation was used to approximate the skewed distribution of the RED data (Fig 4A), with $g = 2.244$ to a normal distribution (Fig 4B), with $g = 0.011$.

The QQ-plot shows the theoretical and estimated distribution deviating from the normal distribution as the RED values increase (Fig 4C), while a normalisation of the distribution is evident after the logarithmic transformation (Fig 4D). With the normality of the distribution restored, it is possible to extend the log-normal ordinary kriging to probability mapping [73].

The drift analysis demonstrated the existence of a non-random (deterministic) component in spatial distribution of the data: the moderate gradient of RED data occurs from north to south regions of the world. However, we considered that the stationarity hypothesis does not hold for the global, but only locally. Thus, in order to make a robust assumption of homogeneity of variances, the concept of process stationarity has been replaced by a stationarity of governing influence regarding local hydrological processes and nearby local anisotropy. In such a situation, ordinary kriging is recommended for interpolation [74].

## Spatial structural modelling

To instruct the kriging interpolation, a regionalisation model was fitted using an iterative procedure developed by Johnston et al. [73], which consists of two steps. Step 1 assumes an isotropic model, and performs an initial run of the experimental spatial structure on the standard deviation-scaled data $z(s_0) = (z(s_0) - \bar{z}) \cdot \sigma^{-1}$, where $z(s_0)$ is used to denote the $j^{th}$ measurement of a variable at the $\alpha^{th}$ spatial locations $s_0$, and $\sigma$ is the sample standard deviation. With step 2 any parameter is calibrated interactively, such as: lag number (assumed equal to 7), lag size $h$ (assumed equal to 10 km for the warning state and 20 km for the alert state), range $a$
representing the limit of spatial dependence (equal to 112 km for the warning state and 70 km for the alert state, which are comparable to the ranges of spatial dependence of extreme precipitation on a global scale, varied between 54 and 265 km [75].

Isotropic semivariograms were then modelled as a combination of two distinct spatial structures, the nugget variance and a spherical structure, as shown in Eq (4):

$$ \gamma(h) = \begin{cases} 0 & h = 0 \\ C_0 + C \left( \frac{3h}{2a} - \frac{1h^3}{2a^3} \right) & 0 < h \le a \\ C_0 + C & h > a \end{cases} \quad (4) $$

The nugget effect $C_0$, equal to 0.066 for the warning state and 0.072 for the alert state, is simply the sum of the measurement error, and the microscale variation of RED, which remains unknown due to the spatial variability associated with the distance between raingauge stations. The value that the semivariogram model reaches at the range (the value on the h-axis) is called the *sill* (partial sill + nugget), and counts 0.182 for the warning state and 0.189 for the alert state, while $h$ is the distance between the unknown point ($h = 0$) and a generic point-station.

In this way, Fig 5 shows the experimental unidirectional semivariogram computed from the 3,625 data of RED, with spherical admissible models fitted for threshold-values $z_k$ (RED > 1.5, and > 3.0 MJ hm⁻² h⁻¹).

Semivariogram values increase with separation distance, reflecting the assumption that RED data that are close tend to be more similar than data that are farther away. In particular, the spherical semivariogram model fluctuates around the *sill* value at 1.05° (~112 km) for the warning state (Fig 5A, violet curve) and at 0.70° (~70 km) for the alert state (Fig 5B, violet curve), suggesting that the phenomenon recorded at alert state operates on a smaller spatial scale than in the warning state. This is physically correct as RED events affect a smaller area as they become more extreme, and the physical mechanisms of extreme storm events are size-dependent [76].

The errors involved in the transfer of information from the point to the landscapes via LNOKpm were assessed by means of quantitative standard error estimation and cross-validation [77], re-estimating the RED data at raingauge locations after removing one RED value at a time from the dataset [78]. The difference between the estimated and the corresponding actual indicator value is the experimental error. Thus, repeating this estimation for the number of the experimental data $n = 3,625$, the cross-validation statistics were calculated, as mean error and root mean square errors (RMSE).

<figure><img src="image_8.png" /><figcaption>Fig 5. Modelling of spatial dependence to instructing kriging interpolation. Experimental semivariogram (dots) with permissible spherical model estimates (violet curve) at a) threshold $z_k > 1.5$ MJ hm⁻² h⁻¹ and b) at threshold $z_k > 3.0$ MJ hm⁻² h⁻¹ (b). Units of the semivariance $\gamma$ are multiplied by 10.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0272161.g005
Cross-validation estimates the proportion of neighbourhood values that are below or above the threshold value. The high proportions of low values not exceeding the threshold and high values exceeding the threshold are a measure of the success of the kriging estimates Fig 6.

This is an estimate of the proportion of the values in the neighbourhood that are above the threshold value. The mean error values of -0.005 (warning state) and -0.019 (alert state), and RMSE = 0.262 (warning state) and RMSE = 0.215 (alert state), are close to zero, making it clear that there are no systematic errors.

## Global spatial pattern of RED-kriged probability estimation

Fig 7A presents the map of the hydrological hazard-prone areas at the warning state, highlighting the fact that about 31% of the world's land area has a greater than 50% probability of exceeding the $z_k$ threshold (RED > 1.5 MJ hm⁻² h⁻¹). The map indicates that the phenomenon taken into account by LNOKpm is not smooth (i.e. RED values change strongly with distance). The most affected hydrological hazard-prone areas are Africa and the southern Asian continents, southern Saudi Arabia, Australia, almost all of the USA with an offshoot to western Canada. These figures are consistent with statistics indicating that about 10 million hectares of cropland are lost each year due to soil erosion worldwide [79], particularly in Asia, Africa and South America, where erosion is more severe [80, 81].

China's far southeast, where the probability of exceeding RED at warning state is high, has experienced significant upward trends in rainfall erosivity over the period 1950–2010 [82]. Similarly, in many parts of Africa, where the warning state is expected to be exceeded over a wide area, soil erosion is becoming a major problem due to the high sediment production in tropical mountain streams [83]. The Mediterranean also has a high probability of reaching a warning state, where the aggressiveness of rainfall [8], in tandem with RED, seems to show a propensity to increase in recent decades [18]. However, in addition to the frequent high erosivity of precipitation regimes and human disturbances, it can be argued that much of the Mediterranean landscape is naturally vulnerable to soil erosion processes [84, 85]. Fig 7B, on the other hand, shows the map of the hydrological hazard-prone areas at alert state, highlighting the fact that ~19% of the world's land area has a greater than 50% chance of exceeding the $z_k$ threshold (RED > 3.0 MJ hm⁻² h⁻¹).

With respect to the alert threshold (Fig 7B), regions that have become free hazard-areas include the Mediterranean lands, almost all of the USA, Japan, Pakistan, northern and central

<figure><img src="image_9.png" /><figcaption>**Fig 6.** Cross-validation for warning and alert states. Scatterplots between actual rainfall erosivity density (RED) values above the given threshold and LNOKpm probability for the thresholds a) $z_k$ (RED) > 1.5 MJ hm⁻² h⁻¹ and b) $z_k$ (RED) > 3.0 MJ hm⁻² h⁻¹. The white vertical lines in both graphs represent the respective RED thresholds (a, warning state; b, alert state). The cross-validation scatter diagrams (a and b) show that the actual RED values below and above the given thresholds at the warning and alert states are in agreement with the respective kriged probability.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0272161.g006
<figure><img src="image_10.png" /><figcaption>**Fig 7.** Global spatial patterns of kriged-probability map over the period 2002–2011. Exceedance of the rainfall erosivity density (RED) threshold-value at a) warning state: $z_k$ (RED) > 1.5 MJ hm$^{-2}$ h$^{-1}$, and b) alert state: $z_k$ (RED) > 3.0 MJ hm$^{-2}$ h$^{-1}$.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0272161.g007

Africa and south-central Australia. Canada, Greenland and Eurasia still remain below the $z_k$ threshold in warning and alert states, because these regions have a climate less exposed to the critical level of RED, and where rainfall erosivity is also lower with a decrease in rainfall intensity as latitude increases northwards. For the areas of the USA exposed to hydrological hazards at the alert stage, the map in Fig 7B roughly overlaps with the map of the country highlighting the hurricane-prone areas [86].

The southernmost part of China, and parts of the southern USA, central and southern Africa, Latin America and India have a very high probability of exceeding the warning threshold. These regions reach soil erosion hotspots of > 20 Mg hm$^{-2}$ yr$^{-1}$, and are among that most intensely eroded areas in the world [87]. In particular, the observation in India follows the small convective systems that dominate throughout the Western Ghats region [88], but large events are also more intense [89].

These findings are consistent with the results of Medeiros et al. [90], who showed that the most perturbed areas in the tropical region are associated with convective storms that can have smaller radii (~10 km), as suggested by the smallest kriged-range obtained with the
semivariogram function exceeding the $z_k$ threshold (RED > 3.0 MJ hm⁻² h⁻¹) at alert state. These are the highest RED values found in most of South America and central and southern Africa, which are characterised by a complex property in the transfer of erosive energy to land. In these regions, sub-grid scale convection and intensification of rainfall generation processes are very hazardous, as long-lived mesoscale convective systems are well organised at these latitudes and contribute disproportionally to extreme tropical precipitation, with ~40% of days with more than 250 mm of rain over land being associated with convective systems lasting more than 24 hours [91].

## Kriged-probability mapping validation at continental spatial scale

The geostatistical approach used is robust to outlier effects even when the number of experimental data is relatively small and irregularly distributed [92], as the few network stations over some regions (e.g. Russia and the African continent). Then, a main drawback that can occur is when rainfall decreases significantly after heavy rainstorms [93]. However, our approach has proven to give satisfactory results when comparing LNOK-based probability maps with the hydrological disasters recorded by the Munich Reinsurance Company (Munich Re, https://www.munichre.com, on a continental scale (data not shown). Thus, for the Asian continent, it was possible to compare the kriged-probability map, for RED > 1.5 MJ hm⁻² h⁻¹, with the damaging hydrological events (storm + erosional soil degradation) that occurred between 1981 and 2018 (Fig 8).

The orange and red colours of the kriged map covering Saudi Arabia, Yemen, Oman, Pakistan, Nepal, Bangladesh, Burma, Vietnam, Thailand, Cambodia, the Indonesian Archipelago, southern Japan and south-eastern China (Fig 8A, red areas) are roughly overlapping with the hazards in multiple locations with the aggregate hydrological impact of the Munich Re dataset (not shown).

The average number of weather, climate and water hazards per decade has increased over the period 1970–2019 in Asia (Fig 8B). They have increased from, on average, one disaster every 15 days to one every three days [94], with a higher proportion of floods and storm events (Fig 8C). On the central and northern Asian continent, the low probability of high RED values (Fig 7), generally associated with a low probability of damaging hydrological events (Fig 8A), is consistent with the rainfall minima of desert areas (e.g. in Mongolia and in northwestern China) and Siberia [95].

For the African continent, spatial dependence may be more difficult to detect because RED data are scarce. Although broad climatic patterns can be identified across the African continent, there are many local variations from place to place, with the most important differentiating climatic factor being rainfall [96]. The continent is most affected by both a continental tropical air mass to the north and by maritime tropical and equatorial air masses to the south (meeting in the Inter-Tropical Convergence Zone) [97]. Essentially, the equatorial maritime air mass is unstable and brings rain, while the tropical maritime air mass, when fully developed, is stable and generally does not bring rain unless forced over a high mountain [98].

However, with the exception of some northern coastal countries and an overestimation for some southern countries, it is still possible to roughly identify the areas most susceptible to damaging hydrological events for African lands below 15° N (Fig 8D, red areas), compared to the flood pattern of the Munich Re dataset (not shown).

This is also in line with increased flooding and storm-driven erosion in Africa (Fig 8E), and nutrient loss due to increased extreme diurnal rainfall observed in tropical eastern and southeastern Africa in the late 20th century [99], calling into question the sustainability of food security for the ~300 million people currently living in Africa south of the equator [100]. Floods account for 60% of recorded disasters, while storms account for 17% (Fig 8F).
<figure><img src="image_11.png" /><figcaption>**Fig 8.** Comparison of kriged probability maps with damaging hydrological events (Asia) and flood events (Africa) over the last four-five decades. a) Rainfall erosivity density (RED) > 1.5 MJ hm⁻² h⁻¹ with b) with trend of disasters and c) related percentage; d) Rainfall erosivity density (RED) > 3.0 MJ hm⁻² h⁻¹ with e) trend of disasters and f) related percentage. The graphs a) and d) refer to the period 2002–2011, graphs b), c), e) and f) to 1971–2019.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0272161.g008

## Conclusion

Detailed rainfall datasets and appropriate modelling approaches are necessary to establish the mechanisms underlying the complex set of physical processes that govern the response of hydrological cycles to a changing climate. Indeed, hydrological hazard responses to climate change remain difficult to quantify with existing modelling frameworks. Although in some areas of the world (such as the African continent) the network of weather stations is still
insufficient, with the widespread availability of high-temporal resolution rainfall records for large areas and modern advances in climate modelling, new opportunities are opening up for the use of geostatistical methods for large-scale planning, hazard management and risk prevention. This paper presents a geostatistical modelling framework for the proper interpretation of spatially-explicit hydrological hazard, which assumes a set of quantitative data for the location of interest (precipitation extremes and rainfall erosivity) and probabilities associated with ranges of rainfall erosivity density above critical values. They include the best available precipitation data from a set of stations (3,625) worldwide (63 countries) and rainfall erosivity data, as provided by an updated version of the Global Rainfall Erosivity Database.

In this study, we have analysed for the first time the spatial pattern of hydrological hazard associated with rainfall erosivity in a global-scale visualisation. The results indicated that about 31% and 19% of the world's land area have a greater than 50% probability of exceeding the warning and alert thresholds of 1.5 and 3.0 MJ hm⁻² h⁻¹, respectively, with the most affected regions being tropical Latin America, South Africa, India and the Indian Archipelago. The geostatistical modelling, designed for a spatial resolution of ~100 km, is compatible with the vast majority of countries in the world (167 out of 234 having an area >10,000 km², https://www.worldometers.info/geography/largest-countries-in-the-world) and within the aggregation range of most environmental and biodiversity models [101]. Highlighting the potential of probabilistic geostatistical modelling, our results suggest the possibility of using geostatistical spatial modelling to determine the probability of exceeding thresholds of high erosivity density and to generate probability maps to delineate the most sensitive areas, which may lead to catastrophic regime shifts related to the occurrence of damaging hydrological events. This soft-computing modelling represents a paradigm shift on how to provide timely, accurate and actionable information on hydrological hazards [102–104]. Without giving the value at each point but returning a probability map, this approach offers the possibility to obtain information also where no data or measurements exist, as it identifies the hydrological hazard associated with the probability of exceeding an erosivity density threshold. In this way, the approach can support decision-making. We thus offer these results as a springboard to support policymakers, local authorities and civil protection in planning medium- and long-term actions to reduce hydrological disasters [105–106]. As rainfall erosivity is projected to increase by at least 35% globally by 2070 [107], the probability of hydrological hazards will also show similar trends.

A more careful validation of the global map of hydrological disaster-prone areas is certainly needed. However, this is a promising first step, and the global probability map was well suited to hydrological disasters in regions where data coverage was substantial. Quantifying the probability of exceeding threshold values of erosivity density in a way that enables meaningful comparisons with hydrological records is an important topic of study, and our article is a step forward towards this goal. In fact, geostatistical methods can be practically implemented to create spatially explicit probabilistic maps at the country level, and can be useful in the study of erosive hazards, which, however, depend on several interacting factors, such as complexography, large-scale air flow and teleconnection patterns. Then, purely geostatistical findings do not produce an explicit mechanistic modelling of rainfall erosivity density, which, however, is data demanding and can be accompanied by a large amount of uncertainties in the estimates. Future studies should assess the results in diverse physical geographic conditions and socioeconomic situations, taking into account that population density, infrastructures, plant density and other factors also influence the occurrence of damage. In addition, probability calculations have also to take into account particularly long periods of low rainfall intensity, which are not erosive but can lead to deadly flooding and landslides.
**Acknowledgments**

The authors would like to acknowledge the following services for proving their data: Bureau of Meteorology in Australia, New Zealand Institute of Water and Atmospheric Research (NIWA), Japan Meteorological Agency (JMA), Korea Meteorological Administration (KMA), National Meteorological Information Center (NMIC) of China, India Meteorological Department (IMD) of the Ministry of Earth Sciences (MoES), Iranian Meteorological Organization (IMO), Directorate General of Civil Aviation (DGCA)—State of Kuwait, Lomonosov Moscow State University (LMSU), Israel Meteorological Service (IMS), Turkish Ministry of Forestry and Water Affairs (MoFWA), U.S. Climate Reference Network (USCRN) and National Oceanic and Atmospheric Administration (NOAA), Comisión Nacional Del Agua (CONAGUA) —Servicio Meteorologico Nacional (SMN) of Mexico, Meteorological Service of Jamaica (MSJ), University of Costa Rica (UCR), Centro Nacional de Investigaciones de Café (CENICAFFÉ) of Colombia, General Directorate of the Water Resources (GDWR) of Chile, Meteorological Department Suriname (MDS), Mauritius Meteorological Services (MMS), Algerian National Agency for Hydraulic Resources (ANRH).

**Author Contributions**

**Conceptualization:** Gianni Bellocchi.

**Data curation:** Nazzareno Diodato, Gianni Bellocchi.

**Formal analysis:** Nazzareno Diodato, Gianni Bellocchi.

**Funding acquisition:** Panos Panagos.

**Investigation:** Nazzareno Diodato.

**Methodology:** Pasquale Borrelli, Panos Panagos.

**Project administration:** Panos Panagos.

**Resources:** Nazzareno Diodato, Pasquale Borrelli.

**Software:** Nazzareno Diodato.

**Supervision:** Pasquale Borrelli, Panos Panagos, Gianni Bellocchi.

**Validation:** Nazzareno Diodato, Pasquale Borrelli.

**Visualization:** Nazzareno Diodato.

**Writing – original draft:** Nazzareno Diodato, Pasquale Borrelli, Panos Panagos, Gianni Bellocchi.

**Writing – review & editing:** Nazzareno Diodato, Pasquale Borrelli, Panos Panagos, Gianni Bellocchi.

**References**

1. Rineau F, Malina R, Beenaerts N, Arnauts N, Bardgett RD, Berg MP, et al. Towards more predictive and interdisciplinary climate change ecosystem experiments. Nat Clim Chang 2019; 9:809–16.

2. Sudmeier-Rieux K, Arce-Mojica T, Boehmer HJ, Doswald N, Emerton L, Friess DA, et al. Scientific evidence for ecosystem-based disaster risk reduction. Nat Sustain 2021; 4:1–8.

3. Xu Z, Pan B, Han M, Zhu J., Tian, L. Spatial–temporal distribution of rainfall erosivity, erosivity density and correlation with El Niño–Southern Oscillation in the Huaihe River Basin, China. Ecol Inform 2019; 52:14–25.

4. Heikkinen RK, Leikola N, Aalto J, Aapala K, Kuusela S, Luoto M, et al. Fine-grained climate velocities reveal vulnerability of protected areas to climate change. Sci Rep 2020; 10:1–11.
5. Panagos P, Borrelli P, Meusburger K, Yu B, Klik A, Lim KJ, et al. Global rainfall erosivity assessment based on high-temporal resolution rainfall records. Sci Rep 2017; 7:1–12.

6. Hu X, Li Z, Nie X, Wang D, Huang J, Deng C, et al. Regionalization of soil and water conservation aimed at ecosystem services improvement. Sci Rep 2020; 10:1–10.

7. Sofia G, Nikolopoulos EI. Floods and rivers: a circular causality perspective. Sci Rep 2020; 10:1–17.

8. Diodato N, Gómara I, Baronetti A, Fratianni S, Bellocchi G. Reconstruction of erosivity density in north-west Italy since 1701. Hydrol Sci J 2021; 66:1185–96.

9. Fischer EM, Knutti R. Observed heavy precipitation increase confirms theory and early models. Nat Clim Chang 2016; 6:986–91.

10. Donat MG, Angélil O, Ukkola AM. Intensification of precipitation extremes in the world's humid and water-limited regions. Environ Res Lett 2019; 14:065003.

11. Tabari H. Climate change impact on flood and extreme precipitation increases with water availability. Sci Rep 2020; 10:1–10.

12. Brunner MI, Swain DL, Wood RR, Willkofer F, Done JM, Gilleland E, et al. An extremeness threshold determines the regional response of floods to changes in rainfall extremes. Commun Earth Environ 2021 2:1–11.

13. Förster K, Thiele LB. Variations in sub-daily precipitation at centennial scale. npj Clim Atmos Sci 2020; 3:13.

14. Lesk C, Coffel E, Horton R. Net benefits to US soy and maize yields from intensifying hourly rainfall. Nat Clim Chang 2020; 10:819–22.

15. Fowler HJ, Lenderink G, Prein AF, Westra S, Allan RP, Ban N, et al. Anthropogenic intensification of short-duration rainfall extremes. Nat Rev Earth Environ 2021; 2:107–22.

16. Zhao X, Li Z, Zhu Q. Change of precipitation characteristics in the water-wind erosion crisscross region on the Loess Plateau, China, from 1958 to 2015. Sci Rep 2017; 7:1–16.

17. Janapati J, Seela BK, Lin PL, Wang PK, Kumar U. An assessment of tropical cyclones rainfall erosivity for Taiwan. Sci Rep 2019; 9:1–14.

18. Diodato N, Ljungqvist FC, Bellocchi G. Fingerprint of climate change in precipitation aggressiveness across the central Mediterranean (Italian) area. Sci Rep 2020; 10:1–13.

19. Bezak N, Borrelli P, Panos P. Exploring the possible role of satellite-based rainfall data to estimate inter- and intra-annual global rainfall erosivity. Hydrol Earth Syst Sci 2021; preprint (https://hess.copernicus.org/preprints/hess-2021-417).

20. Kalnay E, Kanamitsu M, Kistler R, Collins W, Deaven D, Gandin L, et al. The NCEP/NCAR 40-Year Reanalysis Project. Bull. Am. Meteorol Soc 1996; 77:437–72.

21. Petrucci O, Polemio M. Damaging hydro-geological events: definitions, methods and case studies. In: Laughton RH, editor. Aquifers: formation, transport and pollution. Hauppauge: Nova Science Publishers; 2011, p. 235–57.

22. Bankoff G, Hilhorst D, Frerks G. Mapping vulnerability: disasters, development, and people. London: Earthscan; 2004.

23. Cannon T. A hazard need not a disaster make: Vulnerability and the causes of 'natural disasters'. In: Merriman PA, Browitt CWA, editors. Natural disasters: protecting vulnerable communities. London: Thomas Telford; 1993, p. 92–105.

24. Diodato N, Ljungqvist FC, Bellocchi G. A millennium-long reconstruction of damaging hydrological events across Italy. Sci Rep 2019; 9:1–10.

25. Petrucci O, Papagiannaki K, Aceto L, Boissier L, Kotroni V, Grimalt M, et al. MEFF: the database of Mediterranean flood fatalities (1980 to 2015). J Flood Risk Manag 2019; 12:e12461.

26. MacEachren AM, Robinson A, Hopper S, Gardner S, Murray R, Gahegan M, et al. Visualizing geospatial information uncertainty: What we know and what we need to know. Cartogr Geogr Inf Sci 2005; 32:139–60.

27. San BT, Ulusar UD. An approach for prediction of shoreline with spatial uncertainty mapping (SLIP-SUM). Int J Appl Earth Obs Geoinf 2018; 73:546–54.

28. Marrec R, Moniem HEA, Iravani M, Hricko B, Kariyeva J, Wagner HH. Conceptual framework and uncertainty analysis for large-scale, species-agnostic modelling of landscape connectivity across Alberta, Canada. Sci Rep 2020; 10:1–14.

29. Arnell NW. Uncertainty in the relationship between climate forcing and hydrological response in UK catchments. Hydrol. Earth Syst Sci 2011; 15:897–912.
30. Seddon AWR, Macias-Fauria M, Long PR, Benz D, Willis KJ. Sensitivity of global terrestrial ecosystems to climate variability. Nature 2016; 531:229-32. https://doi.org/10.1038/nature16986 PMID: 26886790

31. Naylor LA, Spencer T, Lane SN, Darby SE, Magilligan FJ, Macklin MG, et al. Stormy geomorphology: geomorphic contributions in an age of climate extremes. Earth Surf Process Landforms 2017; 42:166-90.

32. Phillips JD, Van Dyke C. Principles of geomorphic disturbance and recovery in response to storms. Earth Surf Process Landforms 2016; 41:971-79.

33. Yellen B, Woodruff JD, Cook TL, Newton RM. Historically unprecedented erosion from Tropical Storm Irene due to high antecedent precipitation. Earth Surf. Process. Landforms 2016; 41:677-84.

34. Srivastava SK. Threshold concepts in Geographical Information systems: A step towards conceptual understanding. J Geogr High Educ 2013; 37:367-84.

35. Kayler ZE, De Boeck HJ, Fatichi S, Grünzweig JM, Merbold L, Beier C, et al. Experiments to confront the environmental extremes of climate change. Front Ecol Environ 2015; 13:219-25.

36. Foresti L, Pozdnoukhov A, Tuia D, Kanevski M. Extreme precipitation modelling using geostatistics and machine learning algorithms. In: Atkinson PM, Lloyd CD, editors. geoENV VII-Geostatistics for Environmental Applications. Dordrecht: Springer; 2010, p. 41-52.

37. Zhang X, Zwiers FW, Li G, Wan H, Cannon AJ. Complexity in estimating past and future extreme short-duration rainfall. Nat Geosci 2017; 10, 255-59.

38. Talchabhadel R, Nakagawa H, Kawaike K, Prajapati R. Evaluating the rainfall erosivity (R-factor) from daily rainfall data: an application for assessing climate change impact on soil loss in Westrapti River basin, Nepal. Model Earth Syst Environ 2020; 6:1741-62.

39. Kim J, Han H, Kim B, Chen H, Lee JH. Use of a high-resolution-satellite-based precipitation product in mapping continental-scale rainfall erosivity: A case study of the United States. Catena 2020; 193:104602.

40. Dai Q, Zhu J, Zhang S, Zhu S, Han D, Lv G. Estimation of rainfall erosivity based on WRF-derived rain-drop size distributions. Hydrol Earth Syst Sci 2020; 24:5407-22.

41. Borrelli P, Diodato N, Panagos P. Rainfall erosivity in Italy: a national scale spatio-temporal assessment. Int J Digit Earth 2016; 9:835-50.

42. Singh J. & Singh O. Assessing rainfall erosivity and erosivity density over a western Himalayan catchment, India. J Earth Syst Sci 2020; 129:1-22.

43. Panagos P, Ballabio C, Borrelli P, Meusburger K, Klik A, Rousseva S, et al. Rainfall erosivity in Europe. Sci Total Environ 2015; 511:801-14. https://doi.org/10.1016/j.scitotenv.2015.01.008 PMID: 25622150

44. Yue T, Yin S, Xie Y, Yu B, Liu B. Rainfall erosivity mapping over mainland China based on density hourly rainfall records. Earth Syst Sci Data Discuss 2022; 14:665-82.

45. Ballabio C, Borrelli P, Spinoni J, Meusburger K, Michaelides S, Beguería S, et al. Mapping monthly rainfall erosivity in Europe. Sci Total Environ 2017; 579:1298-1315. https://doi.org/10.1016/j.scitotenv.2016.11.123 PMID: 27913025

46. Schmidt S, Alewell C, Panagos P, Meusburger K. Regionalization of monthly rainfall erosivity pattern-sin Switzerland. Hydrol Earth Syst Sci 2016; 20:4359-73.

47. Diodato N, Borrelli P, Panagos P, Bellocchi G, Bertolin C. Communicating hydrological hazard-prone areas in Italy with geospatial probability maps. Front Environ Sci 2019; 7:193.

48. Zhu D, Xiong K, Xiao H. Multi-time scale variability of rainfall erosivity and erosivity density in the karst region of southern China, 1960-2017. Catena 2021; 197:104977.

49. Burrough PA. GIS and geostatistics: Essential partners for spatial analysis. Environ Ecol Stat 2001; 8:361-77.

50. Willmott CJ, Legates DR. Rising estimates of terrestrial and global precipitation. Clim Res 1991; 1:179-86.

51. Moustakis Y, Onof CJ, Paschalis A. Atmospheric convection, dynamics and topography shape the scaling pattern of hourly rainfall extremes with temperature globally. Commun Earth Environ 2020; 1:1-9.

52. Krivoruchko K. Using linear and non-linear kriging interpolators to produce probability maps. Red-lands: Environmental Systems Research Institute; 2001.

53. Wischmeier WH, Smith DD. Predicting rainfall erosion losses: a guide to conservation planning. Agri-culture Handbook no. 537. Washington DC: United States Department of Agriculture; 1978.
54. Renard KG, Yoder DC, Lightle DT, Dabney SM. Universal Soil Loss Equation and Revised Universal Soil Loss Equation. In: Morgan RPC, Nearing MA, editors. Handbook of erosion modelling. Hoboken: Wiley-Blackwell; 2010, p. 137–67.

55. Furmanczyk K, Dudzinska-Nowak J. Effects of extreme storms on coastline changes: a southern Bal- tic example. J Coast Res 2009; 56:1637–40.

56. Panagos P, Ballabio C, Borrelli P, Meusburger K. Spatio-temporal analysis of rainfall erosivity and ero- sivity density in Greece. Catena 2016; 137: 161–72.

57. Ponjiger TM, Lukić T, Basarin B, Jokić M, Wilby RL, Pavić D, et al. Detailed analysis of spatial-tempo- ral variability of rainfall erosivity and erosivity density in the central and southern Pannonian Basin. Sustainability 2021; 13:13355.

58. Dabney SM, Yoder DC, Vieira DAN, Bingner RL. Enhancing RUSLE to include runoff-driven phenom- ena. Hydrol Process 2011; 25:1373–90.

59. Kinnell PIA. Event soil loss, runoff and the Universal Soil Loss Equation family of models: A review. J Hydrol 2010; 385:384–97.

60. Panagos P, Borrelli P, Spinoni J, Ballabio C, Meusburger K, Beguería S, et al. Monthly rainfall erosiv- ity: Conversion factors for different time resolutions and regional assessments. Water 2016; 8:1–18.

61. Huggett RJ. Terrestrial catastrophism: causes and effects. Prog Phys Geogr 1988; 12, 509–32.

62. Allison RJ, Thomas DSG. The sensitivity of landscapes. In: Thomas DSG, Allison RJ, editors. Land- scape sensitivity. Chichester: John Wiley & Sons; 1993, p. 1–5.

63. Phillips JD. Evolutionary geomorphology: Thresholds and nonlinearity in landform response to envi- ronmental change. Model Appl Chaos Theory Mod Sci 2011; 10:196–214.

64. Malhi Y, Adu-Bredu S, Asare RA, Lewis SL, Mayaux P. African rainforests: past, present and future. Phil Trans R Soc B 2013; 368:20120312. https://doi.org/10.1098/rstb.2012.0312 PMID: 23878339

65. Foster GR, Toy TE, Renard KG. Comparison of the USLE, RUSLE1. 06c, and RUSLE2 for application to highly disturbed lands. In: Renard KG, McElroy SA, Gburek W, editors. First Interagency Confer- ence on Research in Watersheds. Washington DC: United States Department of Agriculture; 2003, p. 154–60.

66. Webster R, Oliver MA. Geostatistics for environmental scientists. Chichester: John Wiley & Sons, 2007.

67. Kitanidis PK. Introduction to geostatistics: applications in hydrogeology. Cambridge: Cambridge Uni- versity Press; 1997.

68. Goovaerts P. Geostatistical modelling of uncertainty in soil science. Geoderma 2001; 103:3–26.

69. Burrough PA, McDonnell RA. Principles of geographical information systems: Spatial information sys- tems and geostatistics. Oxford: Oxford University Press; 1998.

70. Krivoruchko K, Gribov A, Ver Hoef J. Predicting exact, filtered, and new values using Kriging. In: Yarus J, Chambers R, editors. Stochastic modeling and geostatistics 2; Tulsa: AAPG Computer Applica- tions in Geology; 2000.

71. Webster R. Statistics to support soil research and their presentation. Eur J Soil Sci 2001; 52:331–40.

72. Kishné AS, Bringmark E, Bringmark L, Alriksson A. Comparison of ordinary and lognormal kriging on skewed data of total cadmium in forest soils of Sweden. Environ Monit Assess 2003; 84:243–63. https://doi.org/10.1023/a:1023326314184 PMID: 12807263

73. Johnston K, ver Hoef JM, Krivoruchko K, Lucas N. Using ArcGis Geostatistical Analyst. Redlands: Environmental Systems Research Institute; 2001.

74. Journel AG. Nonparametric estimation of spatial distributions. J Int Assoc Math Geol 1983; 15, 445– 68.

75. Tan X, Wu X, Liu B. Global changes in the spatial extents of precipitation extremes. Environ Res Lett (2021); 15:054017.

76. Wasko C, Sharma A. Global assessment of flood and storm extremes with increased temperatures. Sci Rep 2017; 7:1–8.

77. Journel AG, Huijbregts C. Mining geostatistics. London: Academic Press; 1978.

78. Isaaks EH, Srivastava RM. An introduction to applied geostatistics. Oxford: Oxford University Press; 1989.

79. Pimentel D. Soil erosion: A food and environmental threat. Environ Dev Sustain 2006; 8:119–37.

80. Mondal A, Khare D, Kundu S, Mukherjee S, Mukhopadhyay A, Mondal S. Uncertainty of soil erosion modelling using open source high resolution and aggregated DEMs. Geosci Front 2017; 8:425–36.
81. Li X, Ye X. Variability of rainfall erosivity and erosivity density in the Ganjiang River Catchment, China: Characteristics and influences of climate change. Atmosphere 2018; 9:1–18.

82. Qin W, Guo Q, Zuo C, Shan Z, Ma L, Sun G. Spatial distribution and temporal trends of rainfall erosivity in mainland China for 1951–2010. Catena 2016; 147:177–86.

83. Kroese JS, Batista PVG, Jacobs SR, Breuer L, Quinton JN, Rufino MC. Agricultural land is the main source of stream sediments after conversion of an African montane forest. Sci Rep 2020; 10:1–15.

84. Woodward JC. Pattern of erosion and suspended sediment yield in Mediterranean river basins. In: Foster I, Foster IDL, Gurnel A, Webb B, editors. Sediment and river quality in river catchments. Chichester: Jonh Wiley & Sons; 1995, p. 365–89.

85. Panagos P, Ballabio C, Poesen J, Lugato E, Scarpa S, Montanarella L, et al. A soil erosion indicator for supporting agricultural, environmental and climate policies in the European union. Remote Sens 2020; 12:1–21.

86. Zachry BC, Booth WJ, Rhome JR, Sharon TM. A national view of storm surge risk and inundation. Weather Clim Soc 2015; 7:109–17.

87. Borrelli P, Robinson DA, Fleischer LR, Lugato E, Ballabio C, Alewell C, et al. An assessment of the global impact of 21st century land use change on soil erosion. Nat Commun 2017; 8:2013. https://doi.org/10.1038/s41467-017-02142-7 PMID: 29222506

88. Bhat GS, Kumar S. Vertical structure of cumulonimbus towers and intense convective clouds over the South Asian region during the summer monsoon season. J Geophys Res 2015; 120:1710–22.

89. Nikumbh AC, Chakraborty A, Bhat GS. Recent spatial aggregation tendency of rainfall extremes over India. Sci Rep 2019; 9:1–7.

90. Medeiros B, Clement AC, Benedict JJ, Zhang B. Investigating the impact of cloud-radiative feedbacks on tropical precipitation extremes. npj Clim Atmos Sci 2021; 4:1–10.

91. Roca R, Fiolleau T. Extreme precipitation in the tropics is closely associated with long-lived convective systems. Commun Earth Environ 2020; 1:1–6.

92. Pardo-Igúzquiza E. Optimal selection of number and location of rainfall gauges for areal rainfall estimation using geostatistics and simulated annealing. J Hydrol 1998; 210:206–20.

93. Seo DJ, Smith JA. Characterization of the climatological variability of mean areal 665 rainfall through fractional coverage. Water Resour Res 1996; 32:2087–95.

94. Douris J, Kim G. The atlas of mortality and economic losses from weather, climate and water extremes (1970–2019). Geneva: World Meteorological Organization; 2021.

95. Dando WA. Asia, climates of Siberia, Central and East Asia. Encyclopedia of World Climatology. Dordrecht: Springer; 2005.

96. Mitchell T. Africa rainfall climatology. Seattle: University of Washington; 2001.

97. Smith DE, Oliver JE. Tropical and equatorial climates. Dordrecht: Springer; 1987.

98. Miller AA. Air mass climatology. Geography 1953; 38:55–67.

99. Donat MG, Lowry AL, Alexander LV, O’Gorman PA, Maher N. More extreme precipitation in the world’s dry and wet regions. Nat Clim Chang 2016; 6:508–13.

100. Pohl B., MacRon C. & Monerie P. A. Fewer rainy days and more extreme rainfall by the end of the century in Southern Africa. Sci Rep 2017; 7:6–12.

101. Amatulli G, Domisch S, Tuanmu M-N, Parmentier B, Ranipeta A, Malczyk J, et al. A suite of global, cross-scale topographic variables for environmental and biodiversity modeling. Sci Data 2018; 5:180040. https://doi.org/10.1038/sdata.2018.40 PMID: 29557978

102. Diodato N. Geostatistical uncertainty modelling for the environmental hazard assessment during single erosive rainstorm events. Environ Monitor Assess (2005); 105:25–42. https://doi.org/10.1007/s10661-005-2815-x PMID: 15952510

103. Diodato N. Spatial uncertainty modeling of climate processes for extreme hydrogeomorphological events hazard monitoring. J Environ Eng 2006; 132:1530–38.

104. Alipour ZT, Mahdian MH, Pazira E, Hakimkhani S, Saeed M. The determination of the best rainfall erosivity index for Namak Lake basin and evaluation of spatial variations. J Basic Appl Sci Res 2012; 2:484–94.

105. de Medeiros ES, de Lima RR, de Olinda RA, Costa dos Santos CA. Modeling spatiotemporal rainfall variability in Paraíba, Brazil. Water 2019; 11:1843.

106. Guardiola-Albert C, Díez-Herrero A, Cuervo-Arango MA, Bodoque JM, García JA, Naranjo-Fernández N, et al. Analysing flash flood risk perception through a geostatistical approach in the village of Nava-luenga, Central Spain. J Flood Risk Manag 2020; 13:e12590.
107. Panagos P., Borrelli P., Matthews F., Liakos L., Bezak N., Diodato N. et al., 2022. Global rainfall erosivity projections for 2050 and 2070. Journal of Hydrology, 610, Art.no.127865