Spatiotemporal geostatistical analysis of precipitation combining ground and satellite observations

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 km²; 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).

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.

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.

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 Equation (9), C is the sample residual covariance matrix at the observation locations, X is the sample matrix of covariates at the observation locations, x₀ is the vector of covariates at the prediction locations, and c₀ is the residual covariance vector between the observation and prediction locations. The variance is used to quantify the uncertainty associated with the kriging predictions.

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.

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.

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.

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 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.

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.

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.

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).

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.

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.

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.

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.

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/.