1. Introduction

In the first part of this series of two articles [1], the process of frequency analysis of extreme daily maximum rainfall data for Madagascar using the Block Maxima method was studied for the 1970-2023 and 2000-2023 time series, and the results for the six provincial capitals were presented as an example. In particular, these results showed that stationary GEV models produced unrealistic results and that the return periods for the same extreme event were different because they depended on the size of the time series considered. In addition, despite the fact that the two time series were nested, it was also found that the trends, reflecting climate change, were different and even in the opposite direction in some cases.

For better statistical modelling, the longest series (1970-2023) had then been studied and, to take account of climate change, non-stationary GEV distribution models were applied to this time series with variation as a function of time of the position parameter and the scale parameter but keeping the shape parameter constant. This resulted in eight possible non-stationary models depending on the values of the intercepts and and the values of the slopes and . These models also included the special case of the Gumbel distribution when the slopes were zero, i.e. and . This led to much more realistic results, but in which the one-to-one relationship between the return period and the return level was called into question by the presence of a third parameter interpreted as the risk level.

Even if this process of frequency analysis were absolutely necessary to determine the return period and/or the return level at a given spatial point, it is extremely cumbersome and tedious and may not be within everyone's reach, particularly for engineers in charge of establishing design flows for hydraulic and road infrastructures, etc. In this case, a map of frequency isohyets is much more practical because, for a return period and an estimated level of risk, all you have to do is identify the spatial point concerned on the map and read the value of the corresponding return period. This would considerably reduce the time spent on this activity in projects and studies in various fields such as flood risk assessment and mitigation, environmental impact studies, land-use planning, river restoration and the design of hydraulic structures, etc. [2].

Such maps are produced using spatial interpolation methods because of the presence of ungauged areas, which is the case for virtually all of Madagascar, and the entire frequency analysis process described above was carried out beforehand on the 132 measurement points shown in Figure 1. These points were chosen at random but so that there was at least 1 measurement point in all 119 districts of Madagascar.

Figure 1. Locations of rainfall data sites on the island of Madagascar

As already indicated in the first article, the overall aim of the study is to update this type of isohyet map, whose most recent values for the whole of Madagascar, to our knowledge, date from 1993 [3]. There are practically no publications on Madagascar either on isohyets or on spatial interpolation and the only recent publication we have found on this subject is an article dating from 2018 which deals with the use of historical normals to improve the modern monthly climate normal surfaces for Madagascar [4].

The literature describes several possible spatial interpolation methods, including purely geometric methods such as Thiessen polygons, deterministic methods such as Inverse Distance Weighting and Thin Plate Spline surfaces, and especially geostatistical methods such as the various forms of kriging. The vast majority of these studies have generally concluded that kriging methods are superior, in various fields such as groundwater studies, extreme or non-extreme rainfall, etc., and with various time steps, monthly or annual [5-16]. In line with this conclusion, this article considers only kriging interpolation methods.

As already indicated above, the use of non-stationary GEV distribution models brings into play, for each return period, a third parameter interpreted as being the notion of risk [17] and which corresponds to the uncertainty limits of the return levels derived from the a posteriori probabilities [18]. As such, these uncertainty limits vary from 0% to 100% and it is not possible to draw up isohyet maps for a return period for the full extent of these uncertainty limits. In this article, therefore, we will confine ourselves to calculating and then plotting the frequency isohyets for two risk levels: the median level (probability = 50%) and the very high level (probability = 90%), even though some authors suggest that the probability to be assigned to a given return period for a very high risk level is the third quartile (probability = 75%) for the final return level [17].

The isohyet maps envisaged for this work concern return periods of 10, 20, 50, 100 and 200 years, but with two risk levels (M: median and VH: very high) per return period, this leads to 10 isohyet maps. However, for lack of space, only the procedures and results for the 10-year return period VH will be presented in this article, but the procedures are strictly identical for the other return periods, with obviously different results. The 10-year return period was chosen because it is generally the contractual return period for engineering works in rural Madagascar.

2. Materials

2.1. Rainfall Data

The rainfall data are those used in the first article, i.e. daily rainfall data from 1970 to 2023, combined with GLDAS (1970 to 2000), TRMM (2001 to 2019) and GPM (2020 to 2023) data. They were obtained from gridded data, but reduced to the measurement points (Figure 1) and downloaded from the https://giovanni.gsfc.nasa.gov/website.

2.2. Digital Elevation Model (DEM)

The DEM used for this study was a set of DEMs produced by SRTM (Shuttle Radar Topography Mission) and available at https://earthexplorer.usgs.gov/. After being assembled, the various tiles were projected onto the CSR EPSG: 29702, which is the CSR covering the entire island of Madagascar.

2.3. Spatial Interpolation Grid

The spatial interpolation grid used for this study was a grid with square cells of dimensions 5000 m × 5000 m, which provides sufficient accuracy given the large size of the country (1700 km × 500 km) while maintaining a certain fluidity in the calculations.

2.4. Calculation Tools

All calculations relating to this study were carried out using the free software R 4.3.2 (https://www.r-project.org/ ) with the corresponding packages, mainly the gstat package [19] and its dependencies. Other packages were also used for the various operations required.

3. Methodology

3.1. Preliminary Analysis

In order to determine the return levels corresponding to the different return periods = 10, 20, 50, 100 and 200 years, the entire procedure described in the first article was applied to each of the 132 measurement points shown in Figure 1 with non-stationary GEV distribution models.

3.2. Exploratory Analysis and Spatial Structure of the Data

In this analysis, the aim was to examine both the statistical distribution of the various values of (24-hour rainfall with a return period of years) and the spatial distribution of these values in order to detect a spatial trend. Indeed, as we propose to carry out an interpolation by kriging, this is only possible if the process that generated the values at the observation points is a random process following a normal distribution and if, spatially, they have a stationary structure (Gaussian fields) [10]. To verify the normality of the data, summary statistics were calculated to account for the asymmetry and dispersion of the data. The presence of a spatial trend was also highlighted using a map of observed points weighted by their values, as well as the polynomial regression of values against X and Y coordinates.

3.3. Normalisation of the Data

For geostatistical modelling to be possible, the data must be (approximately) normally distributed and stationary, i.e. the mean and variance must not vary significantly in space. To make the of approximate to a normal distribution, the well-known Box-Cox transformation was used [20,21]:

(1)

Following this transformation, the y data was replaced by , which follows a normal distribution. The parameter was determined using the maximum likelihood (ML) method. These transformed data were then used for spatial interpolation. To return to the initial data the inverse transformations are immediate by inverting equations (1):

(2)

where is the estimated value of the parameter .

3.4. Spatial Interpolation

Generally speaking and whatever the method used, the objective of spatial interpolation is to determine the value at a non-gauged location, knowing the values at locations where is a location and and are the coordinates in geographical space. Thus, at an ungauged point this value is given by

(3)

where is the weight for gauged station

However, as already indicated above, the spatial interpolation considered in this study is spatial interpolation by kriging, despite the fact that there are other forms of spatial interpolation, such as the very popular IDW [12,22], which we did not use because of the subjectivity involved in choosing the power parameter to be applied to the inverse distance. The kriging interpolation method differs in that it provides an unbiased estimate by considering the autocorrelation of all available data on the basis of various semivariograms [23].

3.4.1. Empirical Semivariogram

Under the assumption of intrinsic stationarity, it is assumed that the process which generated the samples at the measured points is a random function composed of a mean and a residual [24,25]:

(4)

as well as a variogram:

(5)

Under this assumption, the variance of is constant and the spatial correlation of does not depend on the location but only on the separation distance It is then possible to form several pairs whose separation vectors are (almost) identical and estimate the correlation from these. The empirical semivariogram obtained from the transformed data in equation (1) is then defined according to equation (6):

(6)

where is the measured semivariance for distance is the number of data pairs separated by distance and is the value of the variable at location

In addition to what has been done in the exploratory analysis, the empirical semivariogram is a valuable exploratory tool for assessing the presence of spatial correlation in the data. It should also be noted that, under the assumption of intrinsic second-order stationarity, the empirical semivariogram expressed by equation (6) is an unbiased estimator of the theoretical semivariogram expressed by equation (7) below because the unknown mean of the random function is not involved in the expression of the empirical semivariogram in equation (6).

3.4.2. Theoretical Semivariogram Models

As the empirical semivariogram only covers the observed points, it has to be adjusted by a theoretical model to cover the entire space of points to be interpolated. Under the assumption of second-order stationarity and the intrinsic assumption, the semivariogram shows how the dissimilarity between and changes with distance [26]:

(7)

During this study, after a visual examination of the empirical semivariogram, three of the most common theoretical models were tested:

• the exponential model:

(8)

where is the value of the semivariogram when it reaches the sill, and is the range. For this model, the practical range is defined as the distance for which

(9)

i.e. the distance above which the covariance represents only 5% of the value recorded at the origin. This means that

• the spherical model:

(10)

with the same definitions for the sill and the range .

• the Gaussian model:

(11)

with the same definitions for the sill and the range .

3.4.3. Anisotropy

When the theoretical models given by equations (5), (7) and (8) are isotropic, the overall omnidirectional semivariogram depends only on the length of vector , i.e. and being the components of vector projected onto X and Y. However, anisotropy has an important effect on kriging weights because, if it exists, the weight of points located in the main direction of the anisotropy ellipse will increase: it is therefore essential to analyse the presence of anisotropy in the data. When geometric anisotropy exists in the stationary case, it implies that the thresholds of the directional semivariograms, although equal, are reached with two different distances (ranges). The anisotropy ratio is defined as follows [12]:

(12)

where is the largest range in the direction of the major axis and is the smallest span in the direction of the minor axis. The spatial variation between any pair and is then calculated by the distance in the maximum continuity axis variogram model according to the rotation angle and the anisotropy ratio where is given by equation (11) below [12,27]:

(13)

Instead of the omnidirectional empirical semivariogram will then be used to be fitted by the best of one of the theoretical models (equations (7), (9) and (10)) for kriging interpolation.

3.4.4. Ordinary Kriging

We have Ordinary Kriging (OK) when, in equation (4) above, the mean m is constant (or zero) and the kriging is then done on the residuals . For Ordinary Kriging (OK), the weights in equation (3) are obtained by solving the following system of linear equations [12,28]:

(14)

where is the value of the theoretical semivariogram function for the distance between and is the value of the theoretical semivariogram of the distance between and is the Lagrange parameter. The second equation in (14) ensures that the interpolation is unbiased.

3.4.5. Universal Kriging (UK)

When, in equation (4), the mean m of the measured data is no longer stationary but depends on the X and Y coordinates, this confirms the presence of a trend that needs to be removed, and interpolation then needs to be performed on the remaining residuals. Interpolation by kriging can only be performed on random fields that are substantially Gaussian [10]. The predictions at the ungauged locations are then obtained from:

(15)

The Universal Kriging equations are then identical to those of Ordinary Kriging (equations (14)) but with the weights replaced by from equation (15). As for the mean which includes the spatial tendency to be removed, it is generally expressed in the form of a first or second degree polynomial [25,29]:

(16)

or

(17)

where are regression coefficients to be determined and are regression residuals. In the case of equation (16), the trend surface is an affine surface while in the case of equation (17), the trend surface is a quadratic surface.

3.5. Assessment of Model Quality

Given the large amount of modelling to be carried out in this work, it was necessary to have criteria and/or procedures for evaluating the quality of these models in order to choose the best model. The performance of the final theoretical semivariogram model, and hence of the kriging predictions, is assessed using leave-one-out cross-validation [30], a procedure in which the observed values of the data are eliminated one by one and then predictions are made for each of these data using the remaining data [31,32]. The predicted values are then compared with the observed values.

The quality of the interpolation fit was also assessed by comparing the following diagnostic metrics [22,33,34]: the mean error ME, which measures bias, the mean absolute error MAE and the root mean square error RMSE:

(18)

(19)

(20)

where and are the observed and predicted values, respectively, at the same point is the number of observation points. The closer the values of ME, MAE and RMSE are to zero, the better the fit.

4. Results and Discussion

As already indicated, the results shown in the figures are for the 10-year return period T10VH (Very High), but the tables show all the results.

4.1. Results of the Preliminary Analysis

In a non-stationary context and after applying the procedures described in the first article for the 1970-2023 series of annual maxima, the distributions found for the 132 data points are shown in Figure 2. For the GEV distributions, the shape coefficient was such that while for the Gumbel distributions, we had for this shape coefficient.

Traditionally, the Gumbel distribution had always been the distribution used in Madagascar. However, the results show that this is not entirely correct, since 28.7% of the measurement points here follow a GEV distribution, even though the vast majority (71.3%) follow a Gumbel distribution.

Figure 2 also shows that there is no preferred location for the nature of the distribution (GEV or Gumbel) either in relation to coastal areas or in relation to elevation.

Figure 2. Distribution models on a terrain relief background (elevation in metres). Black circles: Gumbel model (87), red squares: GEV model (35)

4.2. Results of the Exploratory Analysis and the Spatial Structure of the Data

The skewness coefficients for the data leading to the 10 maps considered are shown in Table 1.

Table 1. Asymmetry coefficients

The spatial structure of the data for the return period T10VH is shown in Figure 3 below. Figure 3, left-hand side, shows that there is an overall increase in T10VH values in the west-east direction and also in the south-north direction, although there are a few points that do not conform to this observation. This is proof of the presence of a spatial trend, which is confirmed by the polynomial regressions on the right-hand side of Figures 3 and 4. These regression curves show that along X, the polynomial regression curves of order 2 and order 3 are practically identical, whereas along Y, the regressions of order 1 and order 2 are identical. The results of the evaluations of this spatial trend are given below.

Figure 3. Left: Map of proportional values of T10VH. Right: Polynomial regressions of T10VH as a function of X and Y coordinates

4.3. Results of Trend Surface Calculations

As the presence of a trend was highlighted above in Section 4.2, it is important to confirm its existence in the data. Note that if the raw data contains a trend, the same will also be true for the Box-Cox transformed data (BC data). The parameters of the trend surface equations (16) and (17) had been determined by ordinary least squares and are shown in Tables 2 and 3.

Table 2. Parameters of the 1st order trend surfaces of the raw data

Table 3. Parameters of 2nd order trend surfaces of raw data

Graphically, and for example for T10VH data, these trend surfaces are shown in Figures 4 and 5.

Figure 4. Map of the trend surface predicted by the 1st order model and map of the variances of its residuals (T10VH)

Figure 5. Map of the trend surface predicted by the 2nd order model and map of the variances of its residuals (T10VH)

The maps for the other data also show similar patterns (but with different values).

To determine which of the two models to adopt, the residuals from each of these two models were tested for independence (Durbin-Watson test), normality (Shapiro-Wilk test) and homoscedasticity (Breusch-Pagan test) of variances. For all the data, both models passed the independence and homoscedasticity tests, but only the 1st order models passed the normality tests, while the 2nd order models did not. Consequently, the 1st order models were adopted. This is confirmed by Table 3, which shows that the coefficients, to 5, of the second-order terms are all less than 10^-9 for the second-order surface. This decision is consistent with the principle of parsimony in statistical modelling.

Visually, if we refer to Figures 4 and 5, we can see that the high variances are found outside the island, which is normal since these are non-sampled areas.

4.4. Results of Data Normalisation

Given the results in Table 1, a Box-Cox transformation was used to normalise the data according to equations (1) and the λ values found are summarised in Table 4.

Table 4. Values of parameter λ of the Box-Cox transformation

It is noting that many of these values of λ are fairly close to zero, which means that a logarithmic transformation could have been performed [35]. However, in the present study, these values of λ were kept unchanged and it was the first form of equations (1) and (2) that was used during the transformation and back-transformation.

4.5. Anisotropy Study Results

As it is on the BC transformed data that the kriging interpolations will be applied, the examination of the anisotropy was carried out on this BC transformed data. The directional empirical semivariograms for T10VH and for the four main directions in space, 0, 45, 90 and 135 degrees, are shown in Figure 6.

Figure 6. Directional empirical semivariograms of BC data for T10VH. The dotted horizontal line corresponds to the sill (variance) of the BC data

Examination of the ranges of the directional empirical semivariograms in Figure 6 shows that the major axis of the anisotropy ellipse is along the 0° direction (south-north direction) while the minor axis is along the 45° direction (west-east direction). All the other data in the study show patterns similar to Figure 6, but with different values for sills and ranges. The anisotropy ratios (equation (12)) for all the data are summarised in Table 5.

Table 5. Anisotropy ratios

Thus, in accordance with equation (13) and with , the distance to be considered between any pair and in semivariogram calculations should be:

(21)

Nevertheless, the application of this distance did not give any significant differences with the omnidirectional semivariograms taken without consideration of the anisotropy. This can be explained by the fact that the anisotropy may be a consequence of the non-stationarity, even intrinsic, of the random function studied and the differences between the different directional semivariograms may be due to a lack of data in certain directions [26]. This lack of data is clearly visible in Figure 1, where 119 of the 132 points observed were centred on the centroids of the district territories, which are much closer together in the central part than in the northern and western parts of the island. As these stationarity problems are supposed to disappear with a Box-Cox transformation and the removal of the trend, anisotropy was not taken into account in the spatial interpolation.

4.6. Results of Theoretical Semivariogram Models

After removing the 1st order trend, by comparing the sums of squared errors, the fits of the theoretical models to the omnidirectional semivariograms of the BC data all gave as the best model a spherical model with a zero nugget among the theoretical models tested (exponential, Gaussian and spherical). The characteristics of these models are shown in Table 6 for Ordinary Kriging and Table 7 for Universal Kriging.

Table 6. Characteristics of spherical theoretical models of semivariograms of BC data (nugget = 0) for Ordinary Kriging

Table 7. Characteristics of spherical theoretical models of semivariograms of BC data (nugget = 0) for Universal Kriging

Figure 7 shows, for example, the omnidirectional empirical semivariograms and the fitted theoretical models for BC T10VH data.

Figure 7. Empirical semivariograms and fitted theoretical model for BC T10VH data. Left: Ordinary Kriging. Right: Universal Kriging

These results are not entirely surprising, at least for the Gaussian model that has been excluded, since the Gaussian model exhibits parabolic behaviour near the origin, which means that it can be associated with infinitely differentiable second-order stationary functions. However, such regularity is practically non-existent in practical applications, and the Gaussian model is considered unrealistic under normal circumstances [26].

The exponential model, although also having a linear behaviour near the origin like the spherical model, was not retained because its sum of squared errors was higher than that of the spherical model. In addition, the zero value of the nugget reflects a better quality of the data, probably after it had been transformed according to Box-Cox.

4.7. OK and UK Results

The theoretical semivariogram models found in Section 4.6 had been applied for Ordinary Kriging (OK) and Universal Kriging (UK) on the interpolation grid described in Section 2.3 (5000 m × 5000 m). For both Ordinary Kriging and Universal Kriging, this interpolation was carried out on the residuals of the BC data and, for Universal Kriging, the trend surface (affine) was directly taken into account in the model in the form ∼ X + Y and the interpolation was carried out on the residuals of this model.

For both Ordinary Kriging and Universal Kriging, it can be seen that the variances of the prediction residuals are minimal at the observed points, which clearly confirms the unbiased nature of spatial interpolation by kriging. Furthermore, the maximum values of these variances are found outside the island of Madagascar, which is normal since these are unsampled areas.

With the exception of the values, the mapping for the other data in the study shows similar patterns to Figures 8 and 9 and, to summarise, the maximum variances of these predictions are given in Table 8.

Table 8. Maximum prediction variances (mm2) for OK and UK

Figure 8. Ordinary Kriging for T10VH. Left: predicted values. Right: variances of predicted residuals

Figure 9. Universal Kriging for T10VH. Left: predicted values. Right: variances of predicted residuals

Table 8 shows three things: firstly, the maximum variances are lower for interpolation by OK than by UK. Secondly, the values are extremely accurate because, even for the highest value (UK.var = 1.782 mm²), the standard deviation is only 1.335 mm. Thirdly and finally, these maximum values are outside the area of interest, i.e. in the unsampled part of Madagascar. Figures 8 and 9 show, for example, that within Madagascar, the variances for OK and UK lead to standard deviations of less than 1 mm.

4.8. Quality Assessment of the Two Kriging Models

Using the leave-one-out cross-validation (LOOCV) method, the two kriging models were compared: this method compares the predicted values with the observed values and the conclusions on the variances are similar to those obtained for Table 8 above, i.e. Ordinary Kriging presented the lowest variances at the observed points. In relation to the observed points, the two kriging models show exactly the same prediction behaviour, as shown for example in Figure 10 for the T10VH data.

Figure 10. Prediction behaviour of OK and UK models compared with observed values. Red line: model

Figure 10 shows that the model (red line) overestimates the predicted values above the limit, which is 225 mm for the T10VH data example. For all the data, we observed the same overestimation behaviour above a certain limit, but with different values for this limit. These limit values are given in Table 9 below.

Table 9. Limit values beyond which the kriging model overestimates the predicted values

Nevertheless, these overestimation differences can be neglected (of the order of a millimetre), as we have seen with the variances of the residuals of each of these two models, which are outside the area of interest. The quality of the two models was also assessed using diagnostic metrics (equations (18), (19) and (20)), the results of which are shown in Table 10.

Table 10. Comparison of diagnostic metrics for OK and UK models

Table 10 shows that the ME values that measure the bias are all very close to zero, with the values for Ordinary Kriging being the lowest, and the maximum value being 0.0040 mm (T25VH for UK) which confirms once again that kriging is an unbiased interpolator. The MAE and RMSE values are also very close to zero. Therefore, in addition to the comments on prediction variances already made above, the best model for spatial interpolation is the Ordinary Kriging model and this is the one that was used for the final plotting of the frequency isohyets.

4.9. Frequency Isohyet Plot

Based on the rasters produced by Ordinary Kriging, the frequency isohyet maps, i.e. the contour lines, were drawn using the free software QGIS. A simplified version of these maps for the T10VH data is shown in Figure 11.

Figure 11. Overview of the T10VH frequency isohyet map

Figure 11 is shown just to give an idea of the final version of the isohyet maps produced in this work. Full maps with additional graticules and city markers to aid user location are available in PDF format and at larger scales for all data processed. Other useful information is also included on these definitive maps. These maps can be obtained by sending an e-mail to the author.

5. Conclusions

This work on extreme rainfall in Madagascar consists of two articles, the final objective of which was to update the frequency isohyet maps, the most recent of which dates back to 1993. The first article dealt mainly with frequency analysis, which led to a completely different view of the notion of return period, with non-stationary models to take account of the effects of climate change.

In this second article, the focus was on spatial interpolation, where a preliminary analysis of the data had shown that the traditional use of the Gumbel distribution to model extreme annual daily rainfall no longer made sense when the effects of climate change were taken into account. Analysis of the spatial structure had also shown the presence of spatial trends in the data, trends modelled by an affine surface or a quadratic surface depending on the coordinates. However, the various tests showed that this trend was best expressed by an affine surface, which was also in line with the principle of parsimony in statistical modelling.

Spatial interpolation was then carried out on the Box-Cox transformed data, as the raw data showed non-negligible skewness coefficients, and two competing spatial interpolation models were considered: Ordinary Kriging and Universal Kriging. To carry out this geostatistical interpolation, after removing the trend, the directional and omnidirectional empirical semivariograms were analysed and it was shown that, despite the presence of geometric anisotropy, this had no significant influence on the interpolation results. Finally, the best theoretical models fitted to the omnidirectional semivariograms were all spherical models without nuggets, for both Ordinary Kriging and Universal Kriging.

After comparing the variances of the interpolation residuals, leave-one-out cross-validation (LOOCV) and diagnostic metrics, it was found that Ordinary Kriging was better than Universal Kriging and that the final frequency isohyet maps were drawn using the rasters produced by Ordinary Kriging. The cross-validation also showed that the interpolation models found overestimated the predicted values beyond a threshold value, although the deviations were extremely small and outside the area of interest, i.e. outside Madagascar.

In view of recent disasters in various fields (floods, bridge and road failures, etc.), we hope that this work will help to mitigate or even eliminate these disasters by taking into account the correct return levels and/or return periods to be applied in the design of infrastructures or in the management of areas at risk.