Natural Hazards and Earth System Sciences Indices of precipitation extremes in Southern Portugal – a geostatistical approach

Most of the actual studies and previews of future rainfall patterns, based on past observed records for Mediterranean climate areas, focus on the decline of the rainfall amounts over the years, and also on the increase of the frequency of heavy/intense rainfall events particularly in the winter season. These changes in heavy rainfall events may have severe implications and impacts on soil erosion resulting in increased soil degradation risks. The objective of the present work is to evaluate the spatial distribution of extreme precipitation events in Southern Portugal, using a geostatistical approach to assess the relationships between spatial and temporal extreme rainfall patterns. The used dataset comprises a set of 105 stations’ records of daily precipitation within the period 1960–1999. Two indices of extreme precipitation were selected to be computed based on the daily precipitation observation series: one representing the frequency of extremely heavy precipitation events (R30) and another one characterizing flood events (R5D). The space-time patterns of the precipitation indices were evaluated and simulated using a geostatistical approach. Despite no significant temporal trends were detected on the calculated indices series, the space-time decadal patterns are becoming more continuous in the last two decades than the previous ones.


Introduction
In Mediterranean climate regions, precipitation patterns are highly variable concerning time, space, amount and duration of the events.Most of the studies and predictions of future Correspondence to: R. Durão (rmdurao@ist.utl.pt)precipitation patterns, based on past observed records and climate change scenarios for the Mediterranean Basin, indicate a decrease of rainfall amounts over the years, and an increase of the frequency of heavy/intense rainfall events in autumn and winter seasons, particularly in the winter (Brunetti et al., 2001;Kostopoulou and Jones, 2005).In some regions, the increasing number of heavy precipitation events increases the flood risk (Hidalgo et al., 2003).Issues such as drought and erosive rainfall have been raising concern about the risks of land degradation and desertification (Lázaro et al., 2001).These changes in heavy precipitation events may also have severe implications and impacts on soil erosion resulting in increased soil degradation risks.
A common tool to understand and assess the precipitation patterns over a region is to use extreme precipitation indices based on daily precipitation series, as indicators of climate change, (Jones et al., 1999;Karl et al., 1999;Brunetti et al., 2001).A considerable number of extreme precipitation indices are described and analyzed in the literature.Generally, these indices can be split in two main categories: one involves arbitrary fixed thresholds, such as the number of days with daily precipitation exceeding a specific amount or threshold (in mm) (e.g.Klein Tank and Können, 2003;Kostopoulou and Jones, 2005) and the other category is based on statistical quantities such as percentiles, which are more appropriate for regions that contain a broad range of climates (Haylock and Nicholls, 2000;Klein Tank and Können, 2003).The use of extreme precipitation indices are recommended by the joint working group on climate change detection of the World Meteorological Organization -Commission for Climatology (WMO-CCL), the Research Programme on Climate Variability and Predictability (CLIVAR) (Peterson et al., 2001), and the CLIVAR/GCOS/WMO workshop on indices and indicators for climate extremes (Karl et al., 1999) that has released a list of recommended indices.
From this list, two precipitation indices representing wet conditions (R30 and R5D) have been selected in this work.
The R30 index is based on the count of days per year crossing a fixed threshold (30 mm) and R5D which is a durationbased index, defined as the highest consecutive 5-day precipitation total in each year.The use of annual indices greatly simplifies the analysis of extremes and provides useful measures for impact analysis as they can be related with extreme events that affect human society and the natural environment (Klein Tank and Können, 2003).
The majority of studies on extreme precipitation indices only focus on temporal trends rather than space-time patterns and trends, while for some studies a spatial analysis is not feasible due to the sparse number of monitoring stations over large study regions (Costa et al., 2008a;Rodrigo and Trigo, 2007;Kostopoulou and Jones, 2005;Moberg and Jones, 2005;Klein Tank andKönnen, 2003, Serrano et al., 1999).Rodrigo and Trigo (2007) presented an analysis on extreme precipitation changes across the Iberian Peninsula using a seasonal and annual aggregation of observed precipitation data and showed that the main changes occurred in the intensity of northern and southern stations, with a decreasing trend over time, suggesting that further work will be necessary to obtain a complete view of the spatiotemporal variability of daily precipitation in the Iberian Peninsula.Costa et al. (2008b) provide an insight of the spatial distribution of extreme precipitation indices in Southern Portugal, investigating yearly trends and decadal space-time patterns showing that there are no significant trends in the regional extreme indices studied.
In order to accommodate the spatial dimension to precipitation data analysis, some geostatistical studies have already been done, but mainly focused on the spatial interpolation of precipitation fields (e.g.Prudhomme and Reed, 1999;Goovaerts, 2000;Daly, 2006).However, there are very few studies on the modelling of space-time patterns of extreme precipitation indices.Hundecha and Bárdossy (2005) interpolated daily precipitation measurements on a 5 km×5 km grid through external drift kriging, using a digital elevation model as secondary information and, afterwards, several extreme precipitation indices were calculated on grids of 5, 10, 25 and 50 km.
A different approach is proposed for this work since the extreme precipitation indices are calculated directly from daily precipitation measurements and then the spatial distribution of extreme precipitation events is analysed using another geostatistical methodology.The semi-variograms of extreme precipitation series were modelled and the direct sequential simulation was performed to assess the relationships between spatial and temporal extreme precipitation patterns.
A description of the study region and data is presented in Sect.2, and the applied methodology is briefly described in Sect.3. The obtained results are presented and discussed in Sect.4, and the major conclusions are stated in Sect. 5.

Study region
The study region corresponds to the southern part of continental Portugal, which comprises two different administrative regions, Alentejo and Algarve, having different physiographic characteristics.In the far South, the relief is dominated by the two main Algarve's mountains: Monchique on the West, and Caldeirão on the East.In contrast, the Alentejo region is characterized mainly by vast flat to rolling country, the peneplain, where the average altitude is approximately 200 m.The São Mamede mountain ridge, the highest in the Alentejo region with an altitude of 1000 m, lies in the extreme North-East.These regions are characterised by a Mediterranean climate and present a high level of susceptibility to drought events and desertification phenomena (e.g., Pereira et al., 2006).
The precipitation regime in Southern Portugal is characterized by highly irregular behaviour in both the spatial and temporal domains, namely amount and distribution of rainfall (Daveau, 1977).Given the impacts of droughts and floods, the study of extreme precipitation space-time variability is of paramount importance in terms of water basins management in the south of the Iberian Peninsula (Ramos and Reis, 2002;Trigo et al., 2004).
The precipitation regime in Portugal can be explained by two different seasonal atmospheric mechanisms.In summer, the large-scale atmospheric circulation is steered by the Azores anticyclone, which is displaced towards its northwesterly position, producing northerly or northeasterly winds that bring warm and dry air into Portugal, which is either of continental or maritime origins (Trigo and DaCamara, 2000).However, in Southern Portugal, summer precipitation is sometimes associated with local convective activity (Trigo and DaCamara, 2000).On the other hand, during winter, the large-scale circulation is mainly driven by the position and intensity of the Icelandic low, and Portugal is affected by westerly winds that carry moist air and produce rainfall events mainly in northern Portugal (Trigo and DaCamara, 2000).

Data
The station network utilized is an original set of 105 monitoring stations with daily precipitation data (Fig. 1).The data set was compiled from the National System of Water Resources Information database (SNIRH -Sistema Nacional de Informac ¸ão de Recursos Hídricos, http://snirh.inag.pt)and three daily series were compiled from the European Climate Assessment (ECA) dataset (http://eca.knmi.nl).All series data were quality controlled by several procedures (e.g., Costa and Soares, 2008;Costa et al., 2008a).For each station, annual precipitation series were computed and studied for homogeneity through the application of six statistical tests, by means of the hybrid approach proposed by Wijngaard et al. (2003).Besides this, 62% of the long-term series were also checked through relative approaches (testing procedures that use records from reference stations), comprising the application of five homogeneity tests which are capable of locating the year where a break is likely (Costa and Soares, 2008;Costa et al., 2008a).
The final dataset comprises the time series of 105 monitoring stations, distributed irregularly over the study region, with daily precipitation observations within the period 1960-1999.Extreme precipitation indices are sensitive to the number of missing days, thus the selected stations have less than 16% of daily records missing in each year.Therefore, for each station, the indices for a specific year were set to missing if there were more than 16% of the days missing for that year (Haylock and Goodess, 2004).
The length of the precipitation series ranges from 3 to 40 years (Fig. 2).The average length is 23 years, and 28 stations have at least 30 years with precipitation records.Basic statistics of the daily rainfall data are presented in Table 1  in order to characterize the precipitation pattern in the study region within the period 1960-1999.The annual extreme precipitation indices (R30, R5D) were computed using the chosen daily precipitation series.
The index R30 measures the frequency of heavy precipitation events and is defined as the number of days per year with precipitation equal or over 30 mm.This threshold fits the extreme events regime of the region understudy, since it corresponds approximately to the 95% regional-average percentile of the 1961-90 climate normal.
The index R5D is defined as the highest consecutive 5-day precipitation total in each year and provides a measure (in mm) of the magnitude of strong precipitation events.

Methodology for assessing space-time patterns in extreme precipitation indices
Geostatistical estimators, known as kriging, provide statistically unbiased estimates of surface values from a set of observations at recorded locations, using the estimated spatial (and temporal) covariance model of the observed data.Consider the two dimensional problem of estimating a primary variable z at an unsampled location u 0 .Let {z(u α ), α=1, . .., n} be the set of primary data measured at n locations u α .Most of geostatistics is based on the assumption that the set of unknown values is a set of spatially dependent random variables, hence each measurement z(u α ) is a particular realization of the random variable Z(u α ).

R. Durão et al.: Indices of precipitation extremes in Southern Portugal
Kriging uses a linear combination of neighbouring observations to estimate the unknown value at the unsampled location u 0 .This problem can be expressed in terms of random variables as: (1) The optimal kriging weights λ α are determined by solving the kriging equations that result from minimizing the estimation variance while ensuring unbiased estimation of Z(u 0 ) by Ẑ(u 0 ).Kriging methods require a stationarity assumption, expressed in two parts.First, the mean of the process is assumed constant and invariant with spatial location (first order stationarity).Second, the variance of the difference between two values is assumed to depend only on the distance h between the two points, and not on their location u (second order stationarity).Stationarity assumptions on kriging are traditionally accounted for by using local search neighbourhoods so that the dependence on stationarity becomes local (Goovaerts, 1997).
When developing the kriging equations the model of spatial covariances, or variogram (inverse function of the spatial covariances), is assumed known.This is a key function of geostatistics and characterizes the variability of the spatial (and temporal) patterns of physical phenomena.Typically, a mathematical variogram model is selected from a small set of authorised ones (e.g.exponential or spherical) and is fitted to experimental semivariogram values calculated from data for given angular and distance classes.
The experimental semivariogram γ (h) is computed as half the average squared difference between data pairs belonging to a certain angular and distance class: where N (h) is the number of pairs of data locations a vector h apart.
The parameters of the semi-variogram model (sill, range and nugget) are used to assign optimal weights for spatial prediction using kriging.The nugget is determined when h approaches 0. The nugget effect results from high variability at short distances that can be caused by lack of samples, or sampling inaccuracy.
In bounded models (e.g., spherical and exponential), variogram functions increase with distance until they reach a maximum, named sill, at an approximate distance known as the range.The range is the distance h at which the spatial (or temporal) correlation vanishes, i.e. observations separated by a distance larger than the range are spatially (or temporally) independent observations.The range is one of the most important parameters because it is related with the spatial (or space-time) extent of continuity of the phenomenon.

Direct sequential simulation
Direct sequential simulation, which means without any transformation of the original variable, has been used to characterize the space-time variability of phenomena (Russo et al., 2008;Nunes and Soares, 2005), which is the main objective of this work.The aim of sequential simulation is to generate a set of equiprobable realizations of a random field, rather than the most probable realization (given, for example, by least squares interpolators), which reproduce the main spatial patterns as revealed by the spatial covariances and semivariograms inferred by the experimental data.Direct sequential simulation is based on the principle proposed by Journel (1994), and its simple algorithm succeeds in reproducing the variogram and histogram of a continuous variable.
Consider the continuous variable Z(x) with a global cumulative distribution function (cdf) and stationary variogram γ (h).The objective is to reproduce both F z (z) and γ (h) in the final simulated maps.Soares (2001) describes the sequence of the direct sequential simulation (DSS) algorithm of a continuous variable as follows: 1. Randomly select the spatial location of a node z(u 0 ) in a regular grid of nodes to be simulated; 2. Estimate local mean and variance identified with simple kriging estimator and kriging variance.Sample from the global histogram a value z s (u i ) centred in the estimated local mean and variance.
3. Return to step Eq. ( 1) until all nodes have been visited by the random path.

Space-time pattern's methodology
The geostatistical methods described above were used in this paper for the evaluation of space-time patterns of extreme precipitation indices over the study region.Spatial correlations between monitoring stations were generalized in a correlation function of distance between any two points, the variogram, which summarizes the main spatial continuity patterns of the chosen precipitation indices.
The methodology was developed in two stages.In the first one, the extreme precipitation indices were computed for each monitoring station and the space-time correlations (semi-variograms) were obtained for the given period.In the second stage, stochastic simulation (DSS algorithm) is used to evaluate the mean spatial pattern and the local variability of the extreme precipitation indices.
Direct sequential simulation is also used to illustrate the space-time precipitation patterns and the spatial random field simulations, which have the same probability distribution function and spatial semi-variograms (covariances) of the extreme precipitation indices.

Space-time continuity analysis
The R30 and R5D precipitation indices were computed for each monitoring station using data within the period 1960-1999.Having in mind the purpose of detecting spatialtemporal changes in the extreme precipitation indices behaviour, the data were grouped in seven ten-years periods: 1960-69, 1965-74, 1970-79, 1975-84, 1980-89, 1985-94, and 1990-99.Spatial semi-variograms (covariances) were then calculated to assess the spatial correlation between monitoring stations during each ten-year period.
Semi-variograms, which are inverse functions of covariances (inverse measure of the correlation for a given distance h), describe how the spatial continuity changes as a function of the distance and direction (where anisotropy is considered) between any pair of points in space and time.Variogram's values increase with increasing distance of separation until they reach a maximum, named sill, at a distance known as the range.
The semi-variogram models and their parameters are estimated taking into account the experimental semi-variograms that are computed with data values.The main objective is to capture the spatial pattern of physical phenomenon in the variogram model rather than getting a best fit of a second moment (Goovaerts, 1997).Therefore, we chose exponential models that capture the major spatial features of the attribute under study within each decade.
The exponential model is a function of two parameters: the range, here denoted by a, and the sill denoted by C. The exponential model reaches the sill asymptotically: In this type of models, a practical range a is defined as the distance at which the model value is at 95% of the sill: γ (a)=0.95C.
The selected model was fitted to the experimental semivariogram values, not only by giving more importance to the smaller lags and the ones computed from more data pairs, but also taking into account physical knowledge of the area and phenomenon.Considering the results of a thorough analysis on directional semi-variograms, only the omnidirectional semi-variograms of the spatial dimension were retained.Consequently, the spatial variability is assumed identical in all directions (i.e.isotropic) within each ten-year period.
Basic statistics of the daily rainfall data presented in Table 1 show that comparing the studied period data with the climate normal data, all the statistics are quite similar except for the studied period maximum (274.7 mm), which is 50% higher than the climate normal one (186 mm).
The statistics of R5D (Table 2) show that although the mean and median are considerably constant over time, the range of the distribution is increasing, the minimum is decreasing and the maximum is increasing over decades.Also, the variance and standard deviations show an increasing pattern over time.In what concerns R30, the analysis of the basic statistics (Table 3) does not indicate evidence of tendencies over time.Analysis of Table 3 shows that the maximum of the R30 index varies over decades therefore the variance varies over through decades too.
The semi-variogram model parameters (Figs. 3 and 4) of each index are summarized in Tables 4 and 5. Analysis of the R30 parameters (Fig. 3 and Table 4) shows that the spatial range is increasing gradually over decades, from 65 km maximum range in the 1960s to 160 km in 1990s.The continuity of the temporal pattern varies little, but it seems that it has an increasing tendency.R30's sill presents the same behaviour as the variance parameter.
Analysis of the R5D parameters (Fig. 4 and Table 5) shows that the spatial range is also increasing over decades, from 70 km maximum range in the 1960s to 165 km in the 1990s.On the other hand, the continuity of the temporal pattern seems to be slightly decreasing over decades, although this needs to be further investigated.It can be observed that the main change of this index occurs between the 1970s (70 km) and the 1980s (150 km), where the spatial continuity doubles in a decade.Hence, the obtained semi-variograms model parameters show that the spatial continuity of the two indices is increasing over time and have doubled during the studied period, while no relevant tendencies have been detected in the temporal scale.
The ranges of the exponential models fitted to the experimental semi-variograms, which express the extent of spatial continuity of the phenomena, are generally increasing over decades for both indices.This is a very important achievement that can be visualized in the space-time images of the main patterns of extreme precipitation provided by the simulations (see Sect. 4.2).This means that extreme events tend to be more spatially homogeneous through time in this region.

Space-time patterns of extreme precipitation indices
The space-time patterns of the indices can be visualised with maps generated by direct sequential simulation on 800 m×800 m grid cells, using the space-time semi-    The simulation algorithm generates a set of realizations of the spatial phenomenon with the same probability distribution function, spatial covariances and semi-variograms of the data of the R5D and R30 indices.For each ten-year period, 100 equiprobable simulated images per year were computed and the local variability was assessed through variance maps of the equiprobable simulated images set, computed at each grid node.
Since the precipitation regime in Portugal is characterized by a large inter-annual variability, the computation and analysis of average maps per decade permits to evaluate longterm trends.On the other hand, the computation of the variance maps of the equiprobable simulated image allows us to assess the variability associated to the computed average maps.Average equiprobable maps (images) give a mean image of R30 and R5D indices per decade, while the local variability maps enable the quantification of spatial variability/homogeneity of each index per decade.R30 and R5D average maps of the 1970-79, 1980-89 and 1990-99 decades and their associated local variability are presented in Figs. 5  and 6, respectively.
In Fig. 5I, the R30 average values present an identical spatial pattern over the three decades, showing the highest values in mountainous regions of Algarve (Monchique  and for all territory.As expected, the results for R5D (Fig. 6) are consistent and similar and both in terms of averages and variability.

Conclusions
The main objective of this paper is to analyse temporal trends in the spatial patterns of extreme precipitation in Southern Portugal.The spatial patterns of two extreme precipitation indices were characterized for a set of 105 monitoring stations with records within the period 1960-1999.Spatial covariances (semi-variograms) were computed for several tenyear periods, which allow characterizing the space-time continuity of the indices.
The results of the extreme precipitation data analysis show that: -spatial variability, as can be seen from the semivariogram ranges, has decreased over time since the 1960s up to the end of the 20th century.This means that extreme phenomena, characterized by the R5D and R30 indices, have become more spatially homogeneous since then; -the time trend of overall spatial variability can be better visualized from the local variability/homogeneity maps given by geostatistical simulations.
Due to the large variability that characterizes the precipitation regime in Portugal, it is difficult to detect temporal trends, but the spatial-temporal analysis provided by geostatistical tools, allowed us to further understand the precipitation regime dynamics.This study proves that the precipitation regime is becoming more homogenous over the Portuguese southern territory, but the understanding of the global and regional atmospheric processes involved in this change needs further investigation.
It is generally recognized that the NAO has a strong influence on the precipitation regime over Portugal, especially in winter, as well as several modes of atmospheric circulation variability (e.g., Corte-Real et. al., 1998;Trigo and DaCamara, 2000;Trigo et al., 2002;Goodess and Jones, 2002;Trigo et al., 2004;Santos et al., 2005).However, the association between the NAO circulation mode and the surface climate of the European continent has been shown to be non-stationary, since the correlation between the NAO index and local climate variables has changed over time, this indicates that other features at the synoptic and smaller scales (e.g., land use changes) can play an important role.Moreover, Trigo and DaCamara (2000) presented a circulation weather type (CWTs) classification and their influence on the precipitation regime, showing that the cyclonic and westerly type weather have a major contribution to the total and relative amount of precipitation at the Southern regions of Portugal.The same study also suggested that the cyclonic weather class is associated with a fairly homogeneous distribution of precipitation over most of the country, while the "rainy" classes with an Atlantic origin (W, SW and NW) are to be associated to the strong decrease in precipitation from North to South.Accordingly, our results may suggest that the contribution of the cyclonic weather class to the precipitation regime in the southern region is increasing against the westerly weather classes, but this relation needs to be proved.

Fig. 1 .
Fig. 1.Study domain and stations' locations with selected daily precipitation series.

Fig. 2 .
Fig. 2. Distribution of precipitation data series by length.

Table 1 .
Basic statistics of the daily precipitation data (in mm).

Table 2 .
Basic statistics of the R5D index (greatest consecutive 5day precipitation total per year, in mm).

Table 3 .
Basic statistics of the R30 index (number of days per year with daily rainfall above or equal to 30 mm).

Table 4 .
Parameters of the space-time semi-variograms for the R30 index, by decade.

Table 5 .
Parameters of the space-time semi-variograms for the R5D index, by decade.
and Caldeirão), and the lower values in the Alentejo region which is characterized mainly by vast flat to rolling country, the peneplain, where the average altitude is approximately 200 m.But the maps of uncertainty, Fig.5II, clearly present a decreasing variability of the R30 index over the three decades