Spatiotemporal geostatistical analysis of precipitation combining ground and satellite observations

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.


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. ).
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. , ).
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. ). 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. ).
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. ). Especially, Eastern Mediterranean rainfall variability is caused by anomalies associated with the North Atlantic climate variability (Eshel & Farrell ). 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 ; Hertig & Jacobeit ).
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 2 ; 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.  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 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. ; Lebrenz & Bárdossy ). 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 ). 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. 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 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.
where d s i,j is the distance (m), and ϑ is a dimensionless elevation scaling factor used to account for the difference between the horizontal distance in the x-y plane and elevation φ. The second feature is defined from the The experimental spatiotemporal variogram model of the residuals is given by: where r s ¼ kd s i,j k (Equation (2)), r t ¼ |t it j |, and Nðr s , r t Þ 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: where S 0 is the set of sampling points in the search neighborhood of (s 0 , t 0 ),Ẑ 0 (s 0 , t 0 ) is the estimate at an unsampled location time, Z 0 (s i , t i ) is the sampled location- where N 0 is the number of points in S 0 , γ Z 0 (s i , s j ; t i , t j ) is the variogram between two sampled points s i and s j at times t i and t j , γ Z 0 (s j , s 0 ; t j , t 0 ) is the variogram between s j ,t j and the estimation point s 0 , t 0 , and μ is the Lagrange multiplier enforcing the zero-bias constraint. Equation (6) enforces the zero-bias condition.
STRK estimates are expressed as follows: where m Z (s 0 , t 0 ) is the estimated trend function, andẐ 0 (s 0 , t 0 ) is the interpolated residual obtained by means of STOK.
The variance of the spatiotemporal residual kriging estimator is given by: where γ ST denotes the space-time variogram function, while the variance component due to trend uncertainty is given by: In Equation (9) (10)) and the non-separable Spartan space-time variogram function (Equation (11)) are studied.
The sum-metric model is given by the following function where γ ST (r s , 0), γ ST (0, r t ), and γ ST (h) are, respectively, the spatial, temporal, and joint variograms, while α is an anisotropic ratio parameter that approximates the space-timescale 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 ; Varouchakis & Hristopulos ), where η 0 is the scale factor, η 1 is the rigidity coefficient, and ω 1,2 ¼ (jη 1 ∓ Δj=2) 1=2 , Δ ¼ jη 2 1 À 4j 1=2 , ξ stands for the correlation length, h ¼ r=ξ is the normalized lag vector, h ¼ jhj is the separation distance norm, and σ 2 Z is the variance. The spatiotemporal Spartan function is a nonseparable model, which is derived by replacing the spatial 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 γ ST (h) ¼ C Z (0) À C Z (h).

Cross-validation
A cross-validation procedure is applied to assess the predic- Bias (the optimum value is zero; positive or negative sign of the bias denotes, respectively, overestimation or underestimation): Mean absolute error (MAE) (optimum value is 0): Mean absolute relative error (MARE) (optimum value is zero, provided that the observations do not include zero values): whereẐ(s i ) and Z(s i ) are, respectively, the estimated and true values at the point s i .

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 hypoth- where σ 2 Z >0 is the variance (sill), σ 2 nug is the nugget variance term, ξ r > 0 is the spatial scale (or range) parameter, ν > 0 is the smoothness parameter, Γ(Á) is the gamma function, K ν (Á) is the modified Bessel function of the second kind of order ν, and |r s | is the spatial separation distance.
The spherical variogram function is formed in terms of the temporal separation distance r t , where ξ t > 0 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) 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  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)   The parameter ν is the shape coefficient of the Matern model, while the range parameter corresponds to the correlation length of the models.     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   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 ) 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

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.