Abstract

Precipitation data are useful for the management of water resources as well as flood and drought events. However, precipitation monitoring is sparse and often unreliable in regions with complicated geomorphology. Subsequently, the spatial variability of the precipitation distribution is frequently represented incorrectly. Satellite precipitation data provide an attractive supplement to ground observations. However, satellite data involve errors due to the complexity of the retrieval algorithms and/or the presence of obstacles that affect the infrared observation capability. This work presents a methodology that combines satellite and ground observations leading to improved spatiotemporal mapping and analysis of precipitation. The applied methodology is based on space–time regression kriging. The case study refers to the island of Crete, Greece, for the time period of 2010–2018. Precipitation data from 53 stations are used in combination with satellite images for the reference period. This work introduces an improved spatiotemporal approach for precipitation mapping.

HIGHLIGHTS

  • Space-time precipitation trend model incorporating satellite measurements and elevation.

  • Application of 3D distance metric for spatiotemporal prediction.

  • Application of the non-separable Spartan space–time variogram in precipitation data.

  • Space–time residual kriging in geostatistical analysis of non-stationary data.

INTRODUCTION

Climate change is expected to have severe consequences for humans and ecosystems. The impact of climate change will be exacerbated by population growth and economic development. The Mediterranean region is among the ‘hot spots’ to be affected by climate change, which is expected to cause increased frequency and intensity of droughts and hot weather conditions. In the eastern Mediterranean, daytime temperatures are expected to increase and annual precipitation is expected to decrease (Lelieveld et al. 2012). The spatiotemporal modeling and analysis of precipitation and the accurate and reliable estimation of its spatial and temporal variability at ungauged locations can provide useful information. The latter can help to identify spatial patterns and areas with significant shortage of water availability and to assist in water resources balance analysis and recharge under different climate change scenarios. Therefore, new mapping methodologies are necessary that can exploit observations in order to model spatiotemporal variations and correlations. Such methodologies should also be capable of incorporating auxiliary information to produce accurate and reliable maps.

In the Mediterranean region, precipitation regimes exhibit significant heterogeneity and strong seasonal variability in many areas. The intensity, frequency, duration, and total amount of precipitation vary substantially over annual and decadal time scales. The patterns of precipitation and evapotranspiration are expected to be affected by climate change, and these changes will perturb the dynamic balance of the hydrological system. These effects will cause variations in the frequency and intensity of hydrological extremes as well as changes in the spatial–temporal hydrological patterns, with severe impact on water availability, economy, ecosystems, and the environment (Mathbout et al. 2018, 2020).

Precipitation in the Mediterranean region has decreased during the past decades, and the decline has been combined with significant changes in spatiotemporal variability. Increased evaporation and intense rain have been detected in most regions of the Mediterranean, whereas regional mean precipitation is either steady or declining. It has also been observed that rainfall duration is shorter and the average rainfall depth is lower nowadays in the South-Eastern Mediterranean region (Hoegh-Guldberg et al. 2018). The main rainfall events in this region tend to occur between October and May, while heavy precipitation peaks between December and May. Spring and autumn also contribute to a significant amount of precipitation (Tarolli et al. 2012).

One of the major patterns of atmospheric variability and surface climate in the North Hemisphere is the North Atlantic Oscillation (NAO) index. NAO is believed to be one of the most important modes of atmospheric circulation for the Mediterranean region (Dayan et al. 2015). Especially, Eastern Mediterranean rainfall variability is caused by anomalies associated with the North Atlantic climate variability (Eshel & Farrell 2000). Phases of strong positive NAO are associated with below normal precipitation over southern and central Europe. Opposite patterns are observed during strong NAO negative phases (Drake 2014; Hertig & Jacobeit 2014).

This work focuses on the island of Crete, which is located in the southern part of Greece, in the southeastern part of the Mediterranean basin (Figure 1). Crete is the largest island of Greece and the fifth largest island of the Mediterranean. The island has a length of 260 km, width ranging from 12 to 57 km, a coastline of 1,306 km in length, and it covers an area of 8,336 km2; the mean elevation is 460 m and the average slope is 22.8%. Crete has a dry sub-humid Mediterranean climate with long hot and dry summers and relatively humid and cold winters. There are variations among different regions of the island due to distance from the sea and altitude. The lower altitudes have a climate characterized by xerothermic conditions and considerable water deficiency. Crete is characterized by very short spring and exceptionally long summer seasons. Spring starts at the end of March and lasts until May. Summer begins in June and lasts until the end of September with the warmer months being July and August. Autumn starts in October and lasts until the end of December. Winter begins at the end of December and lasts until the end of March. Crete is administratively divided into four prefectures, namely from West to East: Chania, Rethymnon, Heraklion, and Lassithi (Special Water Secretariat of Greece 2017).

Figure 1

Global location of the island of Crete (top). Digital elevation map of Crete including the locations of the meteorological stations (bottom). Please refer to the online version of this paper to see this figure in color: http://dx.doi.org/10.2166/nh.2021.160.

Figure 1

Global location of the island of Crete (top). Digital elevation map of Crete including the locations of the meteorological stations (bottom). Please refer to the online version of this paper to see this figure in color: http://dx.doi.org/10.2166/nh.2021.160.

According to the recent Water Resources Management Plan, the historic average annual precipitation on the island of Crete is approximately 927 mm (Special Water Secretariat of Greece 2017). On the other hand, a recent study for the period 1981–2014 found an average annual rainfall of approximately 798.3 mm (Varouchakis et al. 2018). Two recent studies have assessed the temporal rainfall variability in Crete for the last 40 years (Varouchakis et al. 2018; Agou et al. 2019). Both studies have concluded that for the entire island, annual rainfall rates have, on average, remained almost stable with no significant decreasing trend. Recent studies have also identified a decreasing precipitation gradient from West to the East of the island (Naoum & Tsanis 2003; Vrochidou & Tsanis 2012; Varouchakis et al. 2018; Agou et al. 2019). The rainfall variability in Crete was found to be related to the NAO index variability: a statistically significant negative correlation of the island's rainfall variability with the NAO has been identified for the years 1981–2015 (Varouchakis et al. 2018). The number of rainy days on Crete has decreased during recent years by almost 15% compared with the early 80s, while the intensity of rainfall events has increased (Varouchakis et al. 2018). This observation agrees with a similar finding for the entire country (Nastos et al. 2013). An extensive precipitation monitoring network has been established on the island (Figure 1). This network, installed and operated from the National Observatory of Greece (Lagouvardos et al. 2017), provides the opportunity, by applying an appropriate methodology, to analyze the space–time dynamics of precipitation and to estimate its spatial variability at ungauged locations.

Space–time kriging methods have been successfully used in water resources and hydrology to estimate the spatial and/or temporal variations of physical variables. This work introduces methodology based on space–time kriging principles that combines satellite and ground-based precipitation observations. The incorporation of satellite data helps to improve the spatiotemporal mapping and analysis of precipitation. The proposed methodology employs standard geostatistical features that can enhance the spatiotemporal analysis of precipitation. Such features include the application of an anisotropic 3D distance metric to measure the data interdependence, the introduction of a novel space–time trend model that involves satellite precipitation measurements and high-resolution elevation, and the application of two well-established space–time variogram models for the joint space–time interdependence of precipitation data. The proposed methodology provides an improved model for the space–time dynamics of precipitation variations which can lead to more accurate predictions at unobserved locations and times than existing geostatistical approaches. In addition, it considers explicitly and transparently preliminary geostatistical testing and analysis of the samples, calculation of data interdependence, validation of method's estimation capability, and estimation of uncertainty, in order to deliver accurate and reliable predictions. A common concern in precipitation modeling is how to handle the non-stationarity of precipitation data (Wadoux et al. 2017; Lebrenz & Bárdossy 2019). In this study, we assume that precipitation over the area of interest and for a selected time period is a spatially stationary process (Bárdossy & Pegram 2013). This is a reasonable assumption, given the spatial extent of the island and taking into account that the selected time period is not long enough to observe climate change non-stationary effects.

Space–time kriging has been successfully applied to different datasets such as groundwater level (Hoogland et al. 2010; Júnez-Ferreira & Herrera 2013; Manzione et al. 2019; Ruybal et al. 2019; Varouchakis & Hristopulos 2019; Varouchakis et al. 2019) and precipitation. Especially for the latter, kriging-based approaches have been widely used to investigate the spatial variation of precipitation (e.g. Goovaerts 2000; Durão et al. 2010; Verdin et al. 2015; Varouchakis et al. 2018; Agou et al. 2019; Lebrenz & Bárdossy 2019). Furthermore, space–time kriging has been used to design rainfall networks and analyze precipitation or other meteorological variations in both space and time (Biondi 2013; Hu et al. 2017; Martínez et al. 2017; Raja et al. 2017; Zhang et al. 2018; Hu & Shu 2019; Yang et al. 2019; Cassiraga et al. 2020; Naranjo-Fernández et al. 2020; Takafuji et al. 2020). Rescaled space–time ordinary kriging was compared to copulas with respect to rainfall interpolation (Bárdossy & Pegram 2016). The space–time sum variogram model was applied to analyze rainfall variability (Subyani 2019). These works endorse the applicability of space–time kriging to precipitation data. Research which combines satellite precipitation data with kriging methods (Li & Shao 2010; Tang et al. 2010; Shi et al. 2015; Verdin et al. 2015; Verdin et al. 2016; Wadoux et al. 2017; Yan & Bárdossy 2019) is limited to either the spatial or temporal domains. A space–time application uses as covariates satellite data of a different variable than precipitation such as NDVI (Hu et al. 2017). This work combines ground-based and satellite-derived precipitation data in a space–time framework.

The accurate and detailed estimation of spatiotemporal variability of precipitation on the island of Crete can provide reliable and accurate maps that will be useful for hydrological and climate change studies in the region.

METHODOLOGY

Data availability and management

The available observation dataset consists of daily measurements at 53 meteorological stations (Lagouvardos et al. 2017) randomly distributed over the four prefectures of Crete during the period covering the hydrological years from 2009–2010 until 2017–2018 (see Figure 1 for spatial locations). The precipitation data include both observed cumulative and satellite precipitation measurements at the annual scale (over a hydrological year). Satellite precipitation products for the aforementioned period and elevation (50 × 50 m digital elevation map of Crete is shown in Figure 1) were used as covariates. Space–time modeling of precipitation is performed for the reference observation period 2009–2010 to 2016–2017 to estimate the space–time interdependence of the available data. The method's predictive (forecasting) capability is tested by comparing predictions with the available data for the year 2017–2018.

CMORPH CDR processing system and PERSIANN Cloud Classification System (PERSIANN-CCS) (Nguyen et al. 2019) satellite precipitation products were compared with ground precipitation measurements from the 53 monitoring stations for the aforementioned study period using correlation analysis (Figure 1). The analysis showed that CMORPH data are strongly correlated (Pearson correlation coefficient around 81%) to the ground measurements, while the correlation of the latter with PERSIANN is weaker (Pearson correlation coefficient around 71%). Thus, the CMORPH satellite data were used in this work. The CMORPH CDR processing system generates bias-corrected, integrated satellite precipitation estimates over the global domain. The spatial resolution is set by an 8 km-by-8 km grid and the temporal resolution is 30-min (Xie et al. 2019). Data were downloaded and processed to extract those for the area of Crete under a global project that studies the space–time variation of drought events (Diaz et al. 2020a, 2020b). The derived data were then processed to calculate cumulative annual precipitation for a period of 9 years, 2009/2010–2017/2018, at the aforementioned resolution. Here, the assumption that the calculated precipitation is homogeneously distributed inside each grid cell applies. The derived satellite precipitation dataset is available upon request due to potential copyright of the source info.

Geostatistical method

Space–time geostatistical modeling is based on the joint spatial and temporal dependence between observations in a probabilistic framework. Non-stationarity of observation data in kriging applications is an important concern. However, the non-stationarity present in precipitation data can be approximated using external covariates. The use of covariates, such as the Digital Elevation Map (DEM) and the average precipitation, improves the precipitation estimates compared with those derived by a stationary model; however, Wadoux et al. (2017) showed in their work that increasing the number of covariates does not clearly further improve rainfall estimation. Kriging methods, which incorporate external covariates, such as external drift kriging, universal kriging, and regression kriging, have been used in several applications that involve non-stationary geostatistical interpolation (Lebrenz & Bárdossy 2019). Detrending of precipitation produces a residual component. If the residual can be considered stationary, standard geostatistical methods can be applied. A stationarity test of precipitation time series using the Augmented Dickey Fuller (ADF) test was applied in this work (Dickey & Fuller 1979). After detrending, the geostatistical analysis follows the standard approach: (i) variogram estimation based on the residuals, (ii) theoretical variogram model fitting, (iii) jackknifing in order to validate the variogram estimates, and (iv) kriging estimation of the residuals followed by the addition of the trend model to obtain the interpolated values.

In this work, a spatiotemporal geostatistical analysis of precipitation is performed using space–time regression (or residual) kriging (STRK). The interpolation performance of the method is compared with space–time ordinary kriging. The geostatistical analysis using the proposed methods and tools was performed in Matlab® environment using original code. The spatiotemporal data observed at N space–time coordinates are represented by , while data values in space and time are based on , …, . The spatial model can be decomposed into a mean component representing the trend and a residual corresponding to fluctuations in space–time according to Equation (1). The trend function can be modeled using correlated auxiliary variables (Kitanidis 1997).
formula
(1)
In this work, two geostatistical features are used for the first time in spatiotemporal analysis. The first is a modified anisotropic distance metric instead of the 3D Euclidean distance to account for the fact that monitoring sites are located at significant variable elevation (5–1,250 m). Thus, the proposed distance metric between two points in a 3D coordinate system is applied.
formula
(2)
where is the distance (m), and is a dimensionless elevation scaling factor used to account for the difference between the horizontal distance in the xy plane and elevation . The second feature is defined from the function: and describes the application of a space–time trend (residual) model. The term denotes the regression coefficients, while and are covariates related to observed precipitation. The first covariate denotes the satellite precipitation estimates (over the reference a period of 9 years), and the second covariate is the elevation at the observation locations s. The correlation between the cumulative annual precipitation and topography is equal to 0.73. The methodological steps are presented in the flowchart (Figure 2) and in the following equations.
Figure 2

Flowchart of the proposed methodological process. STRK: space–time residual kriging; STOK: space–time ordinary kriging.

Figure 2

Flowchart of the proposed methodological process. STRK: space–time residual kriging; STOK: space–time ordinary kriging.

The experimental spatiotemporal variogram model of the residuals is given by:
formula
(3)
where = (Equation (2)), = |ti – tj|, and is the number of space–time pairs in the corresponding lags.
The space–time kriging estimator is a weighted linear combination of data values inside a specified space–time neighborhood. Using residual data notation, the estimator is given below:
formula
(4)
where is the set of sampling points in the search neighborhood of , is the estimate at an unsampled location time, is the sampled location-time neighbors, and is the corresponding space–time kriging weights. The kriging weights follow from the minimization of the mean square error of the estimate and are given by the following linear system of equations:
formula
(5)
formula
(6)
where is the number of points in , is the variogram between two sampled points and at times and , is the variogram between , and the estimation point , , and μ is the Lagrange multiplier enforcing the zero-bias constraint. Equation (6) enforces the zero-bias condition.
STRK estimates are expressed as follows:
formula
(7)
where is the estimated trend function, and is the interpolated residual obtained by means of STOK.
The variance of the spatiotemporal residual kriging estimator is given by:
formula
(8)
where denotes the space–time variogram function, while the variance component due to trend uncertainty is given by:
formula
(9)

In Equation (9), C is the sample residual covariance matrix at the observation locations, X is the sample matrix of covariates at the observation locations, x0 is the vector of covariates at the prediction locations, and c0 is the residual covariance vector between the observation and prediction locations. The variance is used to quantify the uncertainty associated with the kriging predictions.

Variography

Two broad families of variogram models commonly used in space–time geostatistics involve separable and non-separable functions. In the first, the space and time correlations are modeled separately using standard variogram functions. The space and time components are then combined under a joint sum, product, or product-sum type function (De Iaco 2010; Heuvelink et al. 2017). For the second, functions that include entangled combinations of space and time are used. These are mainly inspired from dynamic processes (e.g. diffusion equation and Gibbs law) or general mathematical laws (Kolovos et al. 2004; Hristopulos & Tsantili 2017; Varouchakis & Hristopulos 2019). Following standard practice (e.g. De Iaco 2010; Raja et al. 2017), we assess well-known space–time variograms from both categories. More specifically, the sum-metric separable space–time variogram (Equation (10)) and the non-separable Spartan space–time variogram function (Equation (11)) are studied.

The sum-metric model is given by the following function (Hoogland et al. 2010; Heuvelink et al. 2017):
formula
(10)
where , , and are, respectively, the spatial, temporal, and joint variograms, while is an anisotropic ratio parameter that approximates the space–time-scale variance. In the present work, different variogram models were used for the three terms of the sum-metric variogram to determine the optimal space–time model through cross-validation.
The Spartan space–time covariance family is provided by the following group of equations (Hristopulos & Elogne 2007; Varouchakis & Hristopulos 2019),
formula
(11)
where is the scale factor, is the rigidity coefficient, is a dimensionless wavenumber, and , , stands for the correlation length, is the normalized lag vector, is the separation distance norm, and is the variance. The spatiotemporal Spartan function is a non-separable model, which is derived by replacing the spatial lag h with the composite space–time lag , where t is the time lag and the relative weight of the time versus the spatial lag, as in Equation (10). The corresponding variogram model is obtained from the respective covariance function by means of .

Cross-validation

A cross-validation procedure is applied to assess the prediction (forecast) performance of the method. The estimates at the interpolation points are compared with the observations of the selected hydrological year (2017/2018) using the performance metrics listed below.

Bias (the optimum value is zero; positive or negative sign of the bias denotes, respectively, overestimation or underestimation):
formula
(12)
Mean absolute error (MAE) (optimum value is 0):
formula
(13)
Mean absolute relative error (MARE) (optimum value is zero, provided that the observations do not include zero values):
formula
(14)
where and are, respectively, the estimated and true values at the point .

RESULTS AND DISCUSSION

In this study, the spatial and temporal stationarity of the precipitation data was investigated. Spatially stationary data do not exhibit significant spatial trends. However, herein due to the geomorphology of Crete, there is a significant correlation (on average equal to 73%) between precipitation and elevation denoting a significant spatial trend. This means that the precipitation model is spatially non-stationary. The stationarity of the precipitation time series was tested using the ADF test. The ADF test uses the null hypothesis of time stationarity (Dickey & Fuller 1979). The estimated t-values at all monitoring stations are greater than the critical value of the Dickey–Fuller test at α = 0.05, indicating that the series have a unit root. Therefore, the precipitation time series are also non-stationary at every location (Jarvis et al. 2013; Modarres & Ouarda 2014). To account for the space–time non-stationarity, a spatiotemporal trend function, which was based on the satellite precipitation data and elevation, was removed from the ground-based precipitation data.

The assumption of stationarity for the residuals is supported by the behavior of their spatiotemporal variogram (Figure 3) which approaches a constant sill (Gregorich & Carter 2007; Varouchakis 2017; Ruybal et al. 2019; Varouchakis & Hristopulos 2019). The variogram shows a general increasing tendency with the time lag, indicating that precipitation residuals become less similar with increasing temporal distance. The study period of 2009/2010–2016/2017 includes hydrological years that exhibit highly variable rainfall which explains the dissimilarity. The precipitation residuals also become less similar with increasing spatial distance. This is due to the observed East-West and South-North precipitation gradients on the island of Crete (Vrochidou & Tsanis 2012). As the spatial lag increases, the geomorphology of the island and the gradients’ effect leads to less similar precipitation values.

Figure 3

Fit of the theoretical spatiotemporal variogram models (gray) to the experimental variogram estimates (white). The sum-metric variogram fit is presented on the left and the Spartan model fit on the right.

Figure 3

Fit of the theoretical spatiotemporal variogram models (gray) to the experimental variogram estimates (white). The sum-metric variogram fit is presented on the left and the Spartan model fit on the right.

A spatial nugget (discontinuity at zero lag) is present in the experimental variogram. The nugget term usually appears if the target variable fluctuates on scales shorter than the experimental resolution or due to measurement errors (Carrasco 2010). The resolution effect is probably present in our data due to the sparseness of the monitoring network and the complexity of the geomorphology which impacts precipitation. On the other hand, the rainfall stations are monitored by the same organization with similar monitoring systems, thus reducing the potential for measurement errors. However, the residuals are obtained after subtracting a trend function which is based on the satellite precipitation product; it is thus possible that errors inherent in the satellite data are propagated to the residuals. Hence, a nugget term is added to both theoretical space–time models and successfully captures the discontinuity of the experimental variogram at the origin.

The fit of the sum-metric model to the experimental variogram is shown in Figure 3. The optimal sum-metric space–time model is composed of a spatial Matérn model (15), a temporal spherical model (16), and a joint Exponential model (17).

The Matérn variogram function (Matérn 1960; Stein 1999; Pardo-Iguzquiza & Chica-Olmo 2008) is defined as:
formula
(15)
where >0 is the variance (sill), is the nugget variance term, is the spatial scale (or range) parameter, > 0 is the smoothness parameter, is the gamma function, is the modified Bessel function of the second kind of order , and || is the spatial separation distance.
The spherical variogram function is formed in terms of the temporal separation distance ,
formula
(16)
where is the temporal scale (or range) parameter. In addition, the joint variogram function is expressed in terms of the exponential model (17) using the joint distance metric, Equation (10), and a joint scale (or range) parameter .
formula
(17)

The sum-metric model fits well the experimental variogram, although it underestimates the sill. On the other hand, the Spartan variogram model fits (also shown in Figure 3) accurately the overall structure of the experimental variogram.

The Spartan variogram model parameters are as follows: , , , , nugget variance , and elevation scaling factor = 0.0012. The negative value of the rigidity coefficient (Hristopulos & Elogne 2007) allows the model to trace the observed dip of the experimental variogram at 0.28 units of normalized distance. The temporal correlation length is similar to a secondary return period of annual precipitation in Crete (equal to 4 years) which was identified via spectral analysis (Varouchakis et al. 2018).

The parameters of the optimal sum-metric variogram are presented in Table 1. In addition, the parameters of the space–time trend model, which was used to obtain the residuals, are presented in Table 2.

Table 1

Parameters of the sum-metric variogram models

ModelSill (mm2)Range Nugget να
Spatial Matérn 121 0.33 or 86 km 8.02 mm2 0.81 1.12 0.0012 
Temporal spherical 102 2 yr NA 
Joint exponential 116 0.81 NA 
ModelSill (mm2)Range Nugget να
Spatial Matérn 121 0.33 or 86 km 8.02 mm2 0.81 1.12 0.0012 
Temporal spherical 102 2 yr NA 
Joint exponential 116 0.81 NA 

The parameter ν is the shape coefficient of the Matern model, while the range parameter corresponds to the correlation length of the models.

Table 2

Coefficients of the multiple linear regression model for precipitation using satellite precipitation data and elevation as covariates

Trend model coefficientsβ0β1β2
2009/2010 186 0.78 0.44 
2010/2011 208 0.64 0.49 
2011/2012 145 1.12 0.64 
2012/2013 158 0.84 0.57 
2013/2014 161 0.71 0.48 
2014/2015 173 0.81 0.51 
2015/2016 201 0.76 0.54 
2016/2017 167 0.56 0.46 
2017/2018 195 0.62 0.39 
Trend model coefficientsβ0β1β2
2009/2010 186 0.78 0.44 
2010/2011 208 0.64 0.49 
2011/2012 145 1.12 0.64 
2012/2013 158 0.84 0.57 
2013/2014 161 0.71 0.48 
2014/2015 173 0.81 0.51 
2015/2016 201 0.76 0.54 
2016/2017 167 0.56 0.46 
2017/2018 195 0.62 0.39 

The performance of different spatial models with respect to estimating the annual precipitation values of the hydrological year 2017/2018 at the 53 observation locations is presented in Table 3 and in Figure 4. The incorporation of the satellite rainfall data as covariates significantly improves the performance, as STRK is shown to be more accurate than STOK. In addition, as evidenced in the results shown in Table 3, the use of the 3D anisotropic distance metric of Equation (2) significantly improves performance compared with the use of the planar Euclidean distance. It should be noted that the use of the anisotropic 3D distance metric is formally equivalent to the assumption of transverse isotropy which employs different correlation lengths in the vertical and the horizontal directions (Agou et al. 2019; Hristopulos 2020). Regarding variogram functions, the estimation metrics for the Spartan and the sum-metric variogram favor the former. In addition, the compact form of the Spartan variogram model provides faster (by about 36% running Matlab in Windows 10 environment on a computer with CPU: i7-8750, 2.2 GHz, RAM: 16G) space–time interpolation compared with the sum-metric model. Thus, the use of the Spartan variogram function is more suitable for fast mapping requirements.

Table 3

Validation metrics of different spatial models

MethodMAE (mm)MARE (%)Bias (mm)
Space–time residual kriging/sum-metric variogram/3D distance (STRK SM 3D) 46.30 0.09 − 18.20 
Space–time residual kriging/Spartan variogram/3D distance (STRK SP 3D) 39.54 0.07 − 15.78 
Space–time residual kriging/sum-metric variogram/Euclidean distance (STRK SM E) 64.71 0.16 21.70 
Space–time residual kriging/Spartan variogram/Euclidean distance (STRK SP E) 62.23 0.15 18.80 
Space–time ordinary kriging/sum-metric/3D distance (STOK SM 3D) 71.24 0.17 − 26.51 
Space–time ordinary kriging/Spartan variogram/3D distance (STOK SP 3D) 68.15 0.17 − 24.34 
MethodMAE (mm)MARE (%)Bias (mm)
Space–time residual kriging/sum-metric variogram/3D distance (STRK SM 3D) 46.30 0.09 − 18.20 
Space–time residual kriging/Spartan variogram/3D distance (STRK SP 3D) 39.54 0.07 − 15.78 
Space–time residual kriging/sum-metric variogram/Euclidean distance (STRK SM E) 64.71 0.16 21.70 
Space–time residual kriging/Spartan variogram/Euclidean distance (STRK SP E) 62.23 0.15 18.80 
Space–time ordinary kriging/sum-metric/3D distance (STOK SM 3D) 71.24 0.17 − 26.51 
Space–time ordinary kriging/Spartan variogram/3D distance (STOK SP 3D) 68.15 0.17 − 24.34 

The metrics are obtained by comparing the models’ forecasts with the observed annual precipitation values for the hydrological year 2017/2018. The averages are calculated over all locations of the monitoring network.

Figure 4

Validation metrics for the spatial models used in this study. The metrics are obtained by comparing the models’ forecasts with the observed annual precipitation values for the hydrological year 2017/2018. The averages are calculated over all locations of the monitoring network.

Figure 4

Validation metrics for the spatial models used in this study. The metrics are obtained by comparing the models’ forecasts with the observed annual precipitation values for the hydrological year 2017/2018. The averages are calculated over all locations of the monitoring network.

A further analysis of variogram suitability involved the application of a separability test (Cappello et al. 2018; Manzione et al. 2019). Following the steps described in Cappello et al. (2018), a non-separable model was determined suitable for this dataset (the null hypothesis of separability was rejected at confidence level a = 0.05 with the test statistic value: TS = 64.45, and a p-value of 0.0006). This result validates the results of this work and explains the better performance of the Spartan model in spatiotemporal interpolation. Non-stationarity tests have rarely been applied in similar works. The best-performing variogram model is typically selected based on cross-validation analysis. However, the use of non-stationarity tests can provide a priori guidance for selecting separable versus non-separable space–time variograms.

The scatterplot shown in Figure 5 compares the observations with the predicted values and illustrates the accuracy performance of STRK equipped with the Spartan variogram function and the 3D anisotropic distance metric. The plot includes a diagonal line with a 45° slope as a guide to the eye. Perfect agreement would be indicated by all the markers aligned with the line. As evidenced in Figure 5, the estimates are very close to the observed values, with a Pearson correlation coefficient equal to 0.91. In addition, Figure 6 shows the bias plotted versus the Easting coordinate of the monitoring stations, thus illustrating the geographical distribution of the estimation errors. A significant observation is that the locations with the highest bias are mainly located in the central West part of the island, where the monitoring network is sparse. Moreover, the highest errors (in absolute value) are obtained at the locations which exhibit lower correlation values (less than 70%) between ground observations and satellite data. Overall, the bias of cross-validation estimations for the entire dataset is equal to 21.50 mm.

Figure 5

Scatter plot illustrating the cross-validation performance of the STRK method using the Spartan variogram function and the 3D distance metric of Equation (2) for the hydrological year 2017/2018.

Figure 5

Scatter plot illustrating the cross-validation performance of the STRK method using the Spartan variogram function and the 3D distance metric of Equation (2) for the hydrological year 2017/2018.

Figure 6

Geographical distribution of STRK prediction bias for the hydrological year 2017/2018. STRK is applied using the Spartan variogram function and the 3D distance metric of Equation (2).

Figure 6

Geographical distribution of STRK prediction bias for the hydrological year 2017/2018. STRK is applied using the Spartan variogram function and the 3D distance metric of Equation (2).

The spatial distribution of the estimated precipitation over the entire island of Crete based on the optimal STRK model and the associated uncertainty of the estimates for the year 2017/2018 is presented in Figures 7 and 8. The estimated precipitation map, shown in Figure 7, agrees with established patterns, i.e., the relationship of the elevation with the precipitation amount, and the precipitation gradients from East to West and from South to North. The uncertainty map (Figure 8) is clearly affected by the nugget effect which sets a lower bound on the estimation uncertainty. Another notable feature is the higher uncertainty in the Eastern part of the island, especially along the South coast, which reflects the higher variability of the precipitation on the Eastern side. This result is also in general agreement with previous research (Agou et al. 2019).

Figure 7

Spatial distribution of STRK precipitation estimates for the year 2017–2018. STRK is applied using the Spartan variogram function and the 3D distance metric of Equation (2). Please refer to the online version of this paper to see this figure in color: http://dx.doi.org/10.2166/nh.2021.160.

Figure 7

Spatial distribution of STRK precipitation estimates for the year 2017–2018. STRK is applied using the Spartan variogram function and the 3D distance metric of Equation (2). Please refer to the online version of this paper to see this figure in color: http://dx.doi.org/10.2166/nh.2021.160.

Figure 8

Uncertainty of STRK precipitation estimates for the year 2017–2018. STRK is applied using the Spartan variogram function and the 3D distance metric of Equation (2). Please refer to the online version of this paper to see this figure in color: http://dx.doi.org/10.2166/nh.2021.160.

Figure 8

Uncertainty of STRK precipitation estimates for the year 2017–2018. STRK is applied using the Spartan variogram function and the 3D distance metric of Equation (2). Please refer to the online version of this paper to see this figure in color: http://dx.doi.org/10.2166/nh.2021.160.

Another metric which can be used to evaluate the method's performance is the reproduction of spatiotemporal continuity which is reflected in the variogram reproduction. Cross-validation assessments of spatial models often focus on univariate measures of performance, such as those presented above. Stochastic interpolation methods, however, allow comparing the experimental variogram with the variogram obtained from the leave-one-out cross-validation estimates as a test of spatiotemporal continuity reproduction. The spatiotemporal variogram of the leave-one-out cross-validation estimates obtained for the entire dataset (including the estimates for the validation year) using the optimal STRK approach is shown in Figure 9. This variogram exhibits a space–time structure very similar to that of the variogram obtained from the original dataset. This agreement further validates the STRK method's reliable representation of the precipitation space–time structure on the island of Crete.

Figure 9

Comparison of the spatiotemporal experimental variogram obtained from the original precipitation residuals (white) with the experimental variogram of the STRK leave-one-out cross-validation estimates (gray).

Figure 9

Comparison of the spatiotemporal experimental variogram obtained from the original precipitation residuals (white) with the experimental variogram of the STRK leave-one-out cross-validation estimates (gray).

In comparison to similar works focusing on Crete (Naoum & Tsanis 2003; Voudouris et al. 2006; Tzoraki et al. 2015; Varouchakis et al. 2018; Agou et al. 2019) which applied spatial methods of analysis, the spatiotemporal STRK approach presented here provides significantly improved results. Specifically, it leads to a reduction of MAE and MARE by at least 50%. In addition, the fusion of the satellite precipitation and elevation data with the observed precipitation by means of the proposed trend model, the employment of the 3D anisotropic distance metric, and the application of the Spartan variogram on the residuals provided estimation metrics that compete with those achieved by recent similar applications which employ spatiotemporal kriging (e.g. Raja et al. 2017; Yang et al. 2019; Naranjo-Fernández et al. 2020).

The proposed space–time geostatistical methodology can be investigated with other satellite products for different applications in hydrology and the environment to model variables such as temperature, evapotranspiration, NDVI, and soil moisture. In addition, STRK can be used to combine auxiliary information from different sources, e.g., geophysics, for the estimation of groundwater levels in aquifers. Furthermore, it could be employed in hydrological data fusion applications (e.g. Benoit 2021) and as a general framework for incorporating data from different sources (e.g. radar and sensors) and different geostatistical tools (e.g. distance metrics, variograms, data transformation methods, and simulation techniques) to improve the efficiency of spatiotemporal predictions.

CONCLUSIONS

This paper introduces a novel STRK method for the interpolation and short-term forecasting of precipitation data. The key findings of this study are as follows. (i) It is established that the correlation between precipitation and topography is significant in terms of annual average data. (ii) The Spartan spatiotemporal model fits better the precipitation residuals (after a trend function is removed) compared with the sum-metric model, leading to improved interpolation performance (bias and uncertainty estimation). The Spartan variogram models can capture (depending on the value of the rigidity coefficient) wave-like behavior of the experimental variogram. (iii) The anisotropic 3D distance metric used here enhances the predictions’ accuracy. (iv) The proposed trend model improves the variogram fit of the residuals and leads, by means of space–time residual kriging, to improved accuracy compared with ordinary space–time kriging. The aforementioned results apply to the present case study but not necessarily in areas with different geomorphological and climate characteristics.

The basic reason for the improved accuracy of the precipitation estimates achieved in this work is the incorporation of the joint space–time interdependence in the estimation procedure. The latter helps to capture the dynamic behavior of precipitation on the island of Crete in terms of the auxiliary information (satellite precipitation product) and the joint space–time variogram. In addition, the development of a trend model that comprises satellite precipitation data and elevation data as well as the use of the 3D anisotropic distance metric enhanced the predictive performance of the spatiotemporal model.

DATA AVAILABILITY STATEMENT

The precipitation data for the island of Crete can be downloaded from https://www.meteo.gr/crete/. The processed satellite CMORPH precipitation data cannot be made publicly available; readers should contact the corresponding author for details.

ACKNOWLEDGMENTS

The authors would like to specially thank the National Observatory of Athens, Greece for providing precipitation data for the island of Crete: https://www.meteo.gr/crete/.

REFERENCES

Agou
V. D.
Varouchakis
E. A.
Hristopulos
D. T.
2019
Geostatistical analysis of precipitation in the island of Crete (Greece) based on a sparse monitoring network
.
Environmental Monitoring and Assessment
191
(
6
),
353
.
Bárdossy
A.
Pegram
G.
2013
Interpolation of precipitation under topographic influence at different time scales
.
Water Resources Research
49
(
8
),
4545
4565
.
Bárdossy
A.
Pegram
G. G. S.
2016
Space-time conditional disaggregation of precipitation at high resolution via simulation
.
Water Resources Research
52
(
2
),
920
937
.
Benoit
L.
2021
Radar and rain gauge data fusion based on disaggregation of radar imagery
.
Water Resources Research
57
(
2
),
e2020WR027899
.
Biondi
F.
2013
Space-time kriging extension of precipitation variability at 12 km spacing from tree-ring chronologies and its implications for drought analysis
.
Hydrology and Earth System Sciences Discussions
10
,
4301
4335
.
Cappello
C.
De Iaco
S.
Posa
D.
2018
Testing the type of non-separability and some classes of space-time covariance function models
.
Stochastic Environmental Research and Risk Assessment
32
(
1
),
17
35
.
Carrasco
P. C.
2010
Nugget effect, artificial or natural?
Journal of the Southern African Institute of Mining and Metallurgy
110
,
299
305
.
Cassiraga
E.
Gómez-Hernández
J. J.
Berenguer
M.
Sempere-Torres
D.
Rodrigo-Ilarri
J.
2020
Spatiotemporal precipitation estimation from rain gauges and meteorological radar using geostatistics
.
Mathematical Geosciences
53
,
499
516
.
Dayan
U.
Nissen
K.
Ulbrich
U.
2015
Review article: atmospheric conditions inducing extreme precipitation over the eastern and western Mediterranean
.
Natural Hazards and Earth System Sciences
15
(
11
),
2525
2544
.
De Iaco
S.
2010
Space-time correlation analysis: a comparative study
.
Journal of Applied Statistics
37
(
6
),
1027
1041
.
Diaz
V.
Corzo Perez
G. A.
Van Lanen
H. A. J.
Solomatine
D.
Varouchakis
E. A.
2020a
Characterisation of the dynamics of past droughts
.
Science of the Total Environment
718
,
134588
.
Diaz
V.
Perez
G. A. C.
Van Lanen
H. A. J.
Solomatine
D.
Varouchakis
E. A.
2020b
An approach to characterise spatio-temporal drought dynamics
.
Advances in Water Resources
137
,
103512
.
Dickey
D. A.
Fuller
W. A.
1979
Distribution of the estimators for autoregressive time-series with a unit root
.
Journal of the American Statistical Association
74
(
366
),
427
431
.
Drake
J. B.
2014
Climate Modeling for Scientists and Engineers
.
SIAM
,
Philadelphia
,
USA
.
Durão
R. M.
Pereira
M. J.
Costa
A. C.
Delgado
J.
del Barrio
G.
Soares
A.
2010
Spatial–temporal dynamics of precipitation extremes in southern Portugal: a geostatistical assessment study
.
International Journal of Climatology
30
(
10
),
1526
1537
.
Eshel
G.
Farrell
B. F.
2000
Mechanisms of eastern Mediterranean rainfall variability
.
Journal of the Atmospheric Sciences
57
(
19
),
3219
3232
.
Gregorich
E. G.
Carter
M. R.
2007
Soil Sampling and Methods of Analysis
.
CRC Press
,
Boca Raton
.
Hertig
E.
Jacobeit
J.
2014
Variability of weather regimes in the North Atlantic-European area: past and future
.
Atmospheric Science Letters
15
(
4
),
314
320
.
Heuvelink
G. B. M.
Pebesma
E.
Gräler
B.
2017
Space-time geostatistics
. In:
Encyclopedia of GIS
(
Shekhar
S.
Xiong
H.
Zhou
X.
, eds).
Springer International Publishing
,
Cham
, pp.
1919
1926
.
Hoegh-Guldberg
O.
Jacob
D.
Taylor
M.
Bindi
M.
Brown
S.
Camilloni
I.
Diedhiou
A.
Djalante
R.
Ebi
K.
Engelbrecht
F.
2018
Impacts of 1.5 (C global warming on natural and human systems
. In:
Global Warming of 1.5° C. An IPCC Special Report on the Impacts of Global Warming of 1.5° C Above pre-Industrial Levels and Related Global Greenhouse gas Emission Pathways, in the Context of Strengthening the Global Response to the Threat of Climate Change, Sustainable Development, and Efforts to Eradicate Poverty
(
Masson-Delmotte
V.
Zhai
P.
Pörtner
H.-O.
Roberts
D.
Skea
J.
Shukla
P. R.
Pirani
A.
Moufouma-Okia
W.
Péan
C.
Pidcock
R.
Connors
S.
Matthews
J. B. R.
Chen
Y.
Zhou
X.
Gomis
M. I.
Lonnoy
E.
Maycock
T.
Tignor
M.
Waterfield
T.
, eds).
IPCC Secretariat
,
Geneva, Switzerland
, pp.
175
311
.
Hristopulos
D. T.
2020
Random Fields for Spatial Data Modeling
.
Springer/Nature
,
Dordrecht
,
The Netherlands
.
Hristopulos
D. T.
Elogne
S. N.
2007
Analytic properties and covariance functions for a new class of generalized Gibbs random fields
.
IEEE Transactions on Information Theory
53
(
12
),
4667
4679
.
Hristopulos
D. T.
Tsantili
I. C.
2017
Space–time covariance functions based on linear response theory and the turning bands method
.
Spatial Statistics
22
,
321
337
.
Hu
D.-g.
Shu
H.
2019
Spatiotemporal interpolation of precipitation across Xinjiang, China using space-time CoKriging
.
Journal of Central South University
26
(
3
),
684
694
.
Hu
D.
Shu
H.
Hu
H.
Xu
J.
2017
Spatiotemporal regression Kriging to predict precipitation using time-series MODIS data
.
Cluster Computing
20
(
1
),
347
357
.
Kitanidis
P. K.
1997
Introduction to Geostatistics
.
Cambridge University Press
,
Cambridge
.
Kolovos
A.
Christakos
G.
Hristopulos
D. T.
Serre
M. L.
2004
Methods for generating non-separable spatiotemporal covariance models with potential environmental applications
.
Advances in Water Resources
27
(
8
),
815
830
.
Lagouvardos
K.
Kotroni
V.
Bezes
A.
Koletsis
I.
Kopania
T.
Lykoudis
S.
Mazarakis
N.
Papagiannaki
K.
Vougioukas
S.
2017
The automatic weather stations NOANN network of the National Observatory of Athens: operation and database
.
Geoscience Data Journal
4
(
1
),
4
16
.
Lebrenz
H.
Bárdossy
A.
2019
Geostatistical interpolation by quantile kriging
.
Hydrology and Earth System Sciences
23
(
3
),
1633
1648
.
Lelieveld
J.
Hadjinicolaou
P.
Kostopoulou
E.
Chenoweth
J.
El Maayar
M.
Giannakopoulos
C.
Hannides
C.
Lange
M. A.
Tanarhte
M.
Tyrlis
E.
Xoplaki
E.
2012
Climate change and impacts in the Eastern Mediterranean and the Middle East
.
Climatic Change
114
(
3
),
667
687
.
Manzione
R.
Takafuji
E.
De Iaco
S.
Cappello
C.
Da Rocha
M.
2019
Spatio-temporal kriging to predict water table depths from monitoring data in a conservation area at São Paulo State, Brazil
.
Geoinformatics and Geostatistics: An Overview
10
,
2
.
Martínez
W. A.
Melo
C. E.
Melo
O. O.
2017
Median Polish Kriging for space–time analysis of precipitation
.
Spatial Statistics
19
,
1
20
.
Matérn
B.
1960
Spatial variation
.
Meddelanden från Statens Skogsförsöksanstalt Institute
49
(
5
),
1
144
.
Mathbout
S.
Lopez-Bustins
J. A.
Royé
D.
Martin-Vide
J.
Bech
J.
Rodrigo
F. S.
2018
Observed changes in daily precipitation extremes at annual timescale over the eastern Mediterranean during 1961–2012
.
Pure and Applied Geophysics
175
(
11
),
3875
3890
.
Mathbout
S.
Lopez-Bustins
J. A.
Royé
D.
Martin-Vide
J.
Benhamrouche
A.
2020
Spatiotemporal variability of daily precipitation concentration and its relationship to teleconnection patterns over the Mediterranean during 1975–2015
.
International Journal of Climatology
40
,
1435
1455
.
Modarres
R.
Ouarda
T. B. M. J.
2014
Modeling the relationship between climate oscillations and drought by a multivariate GARCH model
.
Water Resources Research
50
(
1
),
601
618
.
Naoum
S.
Tsanis
I. K.
2003
Temporal and spatial variation of annual rainfall on the island of Crete, Greece
.
Hydrological Processes
17
(
10
),
1899
1922
.
Naranjo-Fernández
N.
Guardiola-Albert
C.
Aguilera
H.
Serrano-Hidalgo
C.
Rodríguez-Rodríguez
M.
Fernández-Ayuso
A.
Ruiz-Bermudo
F.
Montero-González
E.
2020
Relevance of spatio-temporal rainfall variability regarding groundwater management challenges under global change: case study in Doñana (SW Spain)
.
Stochastic Environmental Research and Risk Assessment
34
,
1289
1311
.
Nastos
P. T.
Politi
N.
Kapsomenakis
J.
2013
Spatial and temporal variability of the Aridity Index in Greece
.
Atmospheric Research
119
,
140
152
.
Nguyen
P.
Shearer
E. J.
Tran
H.
Ombadi
M.
Hayatbini
N.
Palacios
T.
Huynh
P.
Braithwaite
D.
Updegraff
G.
Hsu
K.
Kuligowski
B.
Logan
W. S.
Sorooshian
S.
2019
The CHRS Data Portal, an easily accessible public repository for PERSIANN global satellite precipitation data
.
Scientific Data
6
(
1
),
180296
.
Raja
N. B.
Aydin
O.
Türkoğlu
N.
Çiçek
I.
2017
Space-time kriging of precipitation variability in Turkey for the period 1976–2010
.
Theoretical and Applied Climatology
129
(
1
),
293
304
.
Ruybal
C. J.
Hogue
T. S.
McCray
J. E.
2019
Evaluation of groundwater levels in the arapahoe aquifer using spatiotemporal regression kriging
.
Water Resources Research
55
,
2820
2837
.
Shi
T.
Yang
X.
Christakos
G.
Wang
J.
Liu
L.
2015
Spatiotemporal interpolation of rainfall by combining BME theory and satellite rainfall estimates
.
Atmosphere
6
(
9
),
1307
1326
.
Special Water Secretariat of Greece
2017
Integrated Management Plans of the Greek Watersheds
.
Ministry of Environment & Energy
,
Athens
.
Stein
M. L.
1999
Interpolation of Spatial Data. Some Theory for Kriging
.
Springer
,
New York
.
Subyani
A. M.
2019
Climate variability in space-time variogram models of annual rainfall in arid regions
.
Arabian Journal of Geosciences
12
(
21
),
650
.
Takafuji
E. H. d. M.
da Rocha
M. M.
Manzione
R. L.
2020
Spatiotemporal forecast with local temporal drift applied to weather patterns in Patagonia
.
SN Applied Sciences
2
(
6
),
1001
.
Tang
L.
Hossain
F.
Huffman
G. J.
2010
Transfer of satellite rainfall uncertainty from gauged to ungauged regions at regional and seasonal time scales
.
Journal of Hydrometeorology
11
(
6
),
1263
1274
.
Tarolli
P.
Borga
M.
Morin
E.
Delrieu
G.
2012
Analysis of flash flood regimes in the north-western and south-eastern Mediterranean regions
.
Natural Hazards and Earth System Sciences
12
(
5
),
1255
1265
.
Tzoraki
O.
Kritsotakis
M.
Baltas
E.
2015
Spatial water use efficiency index towards resource sustainability: application in the island of Crete, Greece
.
International Journal of Water Resources Development
31
(
4
),
669
681
.
Varouchakis
E. A.
Corzo
G. A.
Karatzas
G. P.
Kotsopoulou
A.
2018
Spatio-temporal analysis of annual rainfall in Crete, Greece
.
Acta Geophysica
66
(
3
),
319
328
.
Varouchakis
E. A.
Theodoridou
P. G.
Karatzas
G. P.
2019
Spatiotemporal geostatistical modeling of groundwater levels under a Bayesian framework using means of physical background
.
Journal of Hydrology
575
,
487
498
.
Verdin
A.
Rajagopalan
B.
Kleiber
W.
Funk
C.
2015
A Bayesian kriging approach for blending satellite and ground precipitation observations
.
Water Resources Research
51
(
2
),
908
921
.
Verdin
A.
Funk
C.
Rajagopalan
B.
Kleiber
W.
2016
Kriging and local polynomial methods for blending satellite-derived and gauge precipitation estimates to support hydrologic early warning systems
.
IEEE Transactions on Geoscience and Remote Sensing
54
(
5
),
2552
2562
.
Voudouris
K.
Mavrommatis
T.
Daskalaki
P.
Soulios
G.
2006
Rainfall variations in Crete island (Greece) and their impacts on water resources
.
Publicaciones del Instituto Geologico y Minero de Espana Serie: Hidrogeologia y aguas subterráneas
18
,
453
463
.
Vrochidou
A. E. K.
Tsanis
I. K.
2012
Assessing precipitation distribution impacts on droughts on the island of Crete
.
Natural Hazards and Earth System Sciences
12
(
4
),
1159
1171
.
Wadoux
A. M. J. C.
Brus
D. J.
Rico-Ramirez
M. A.
Heuvelink
G. B. M.
2017
Sampling design optimisation for rainfall prediction using a non-stationary geostatistical model
.
Advances in Water Resources
107
,
126
138
.
Xie
P.
Joyce
R.
Wu
S.
Yoo
S.-H.
Yarosh
Y.
Sun
F.
Lin
R.
& Noaa cdr Program
2019
Noaa Climate Data Record (cdr) of Cpc Morphing Technique (cmorph) High Resolution Global Precipitation Estimates
, 1st edn.
NOAA National Centers for Environmental Information
,
Boulder
.
Yang
X.
Yang
Y.
Li
K.
Wu
R.
2019
Estimation and characterization of annual precipitation based on spatiotemporal kriging in the Huanghuaihai basin of China during 1956–2016
.
Stochastic Environmental Research and Risk Assessment
34
,
1407
1420
.
Zhang
Y.
Zheng
X.
Wang
Z.
Ai
G.
Huang
Q.
2018
Implementation of a parallel GPU-based space-time kriging framework
.
ISPRS International Journal of Geo-Information
7
(
5
),
193
.
This is an Open Access article distributed under the terms of the Creative Commons Attribution Licence (CC BY 4.0), which permits copying, adaptation and redistribution, provided the original work is properly cited (http://creativecommons.org/licenses/by/4.0/).