Application of statistical and geostatistical approaches in temporal and spatial estimations of rainfall

Predicting precipitation under the temporal and spatial distributions has great importance in water engineering, hydraulic structures, and water resources management (Krajewski & Smith 2002; Gräler et al. 2013; Volpi & Fiori 2014; Ghazvinian et al. 2020d, 2021; Ugbaje & Bishop 2020; Dadrasajirlou et al. 2023). The importance of this probabilistic forecast has increased in recent years due to the phenomenon of climate change and land-use change (Dore 2005; Fujita 2008; Flannigan et al. 2016; Konapala et al. 2020; Hu et al. 2021).

The spatio-temporal distribution of precipitation has shown a straightforward relationship for providing accessible water resources for humans and ecosystems (Milly et al. 2005; Oki & Kanae 2006; Mishra et al. 2011). Therefore, changes in precipitation due to human-caused climate change in the 21st century may result in changes in water accessibility (WA) that have implications for both humans and the biosphere (Haddeland et al. 2014; Schewe et al. 2014). Studies have examined several simultaneous changes in the P_year and P_max24h in these variables (Holdridge 1947; Thornthwaite 1948; Kottek et al. 2006). Furthermore, the characteristics of annual rainfall and 24-hour rainfall, even in arid regions, are the main basis for WA studies in existing climate classifications (Feng et al. 2013; Rajah et al. 2014). Also, changes in P_year cause drought, increase water stress in human society, and have adverse effects on agriculture and food production (Vaziri et al. 2018). Thus, 24-h rainfall impacts flood events.

Geostatistical methods are based on spatial variation functions (Eftekhari et al. 2021). Spatial variation functions are generally fitted using parametric methods, such as Kriging and inverse distance weighted (IDW). IDW and Kriging consider the influence of covariates on the predicted variable and use the spatial autocorrelation of the predicted variable and its cross-correlation information with related covariates, but they have different modeling methods and computational processes. St-Hilaire et al. (2003) evaluated the impact of the density of the network of water and meteorological stations on the estimation of annual and daily precipitation of rainfall and runoff events in the northern United States of America. The results showed that increasing the number of stations has an effect on the quality of estimating the spatial distribution of precipitation in the region.

Cheng et al. (2008) evaluated the rain gauge network using geostatistical methods in order to estimate the regional average of precipitation and point estimation in areas without stations. Variogram analysis showed that hourly rainfall has higher spatial changes than annual rainfall.

Therefore, forecasting P_year and P_max24h can help in planning and managing water resources. On the other hand, the phenomenon of climate change has increased the occurrence of marginal precipitation. In addition, due to the increase in the urban area, the runoff coefficient in catchments has increased. These have caused the catchments to be more sensitive to rainfall and, as a result, produce more floods than in previous years. This issue leads to an increasingly irreparable loss of life and property. In rainfall forecasting, a lot of research has been conducted with temporal and spatial distribution approaches and intelligent models (Nasseri et al. 2008; Mekanik et al. 2013; Solgi et al. 2016; Azad et al. 2019; Zhang et al. 2020). Also, studies such as Bayat et al. (2019) optimized rainfall estimation based on spatial variations of rainfall. In recent years, extensive scientific studies have been conducted to predict precipitation by considering the temporal and spatial distribution in different parts of the world.

P_max24h is used to calculate the maximum flow in the design of high-risk structures such as dams and nuclear power plants. Also, by using P_max24h, the maximum possible flood can be calculated. The calculation of P_max24h is based on the extreme and exceptional conditions of factors such as the continuity of rainfall, time distribution of rainfall, and infiltration, which requires a lot of data and information.

Geo-statistics is a branch of statistics that investigates natural phenomena that are both random and structural (Çadraku 2021). Geo-statistics originated from mineral resource assessment, but it has been widely used in many fields such as climate, hydrology, environment, and ecology (Çadraku 2022).

Spatial interpolation is the traditional way to transform point-wise rainfall into areal rainfall. It is a process of tapping and utilizing the spatial autocorrelation of rainfall and spatial intercorrelation between rainfall and related explanatory variables under a specific mathematical framework. The spatial distribution of precipitation and temporal distribution plays a critical role in the implementation of development projects. Using the spatial distribution of rainfall can help to determine the risk of floods or landslides in an area.

Fitting theoretical distributions using real data can be used to obtain the probability of occurrence or return period of natural phenomena. By choosing the appropriate statistical distribution, a better estimation of the return period of extreme phenomena with a low occurrence probability can be made, thus increasing the efficiency and reliability of water projects. Precipitation prediction in basins and rivers with different return periods that have statistical information can be conducted using theoretical statistical distributions (Rao & Hamed 2000). The advantage of matching the frequency distribution of existing data with a statistical distribution is developing the existing limited statistics to understand the course of events in the future. There is no consensus among hydrologists on using a particular distribution function. But choosing the appropriate statistical distributions whose results are close to reality (have the lowest error) is possible.

Therefore, recently a considerable number of research works in spatial distribution have been published. Husak et al. (2007) have studied the monthly rainfall distribution for drought monitoring. They used the Gamma distribution approach, and the results showed it has a good performance in 98% of stations. Xia et al. 2012 used rainfall data from 27 meteorological stations in China's Huai River Basin to study the trend of marginal rainfall in the region. Based on their study, three statistical distributions, namely GEV, GP, and Gamma were used. The GEV distribution was more accurate according to the Kolmogorov–Smirnov (K–S) test.

Bhavyashree & Bhattacharyya (2018) used 20 different probability distributions to fit the maximum daily rainfall in Karnataka, India. Their paper results revealed that probabilistic distributions such as Gamma, Pearson 5, and Weibull performed better.

Alam et al. (2018) examined the best probability distribution for predicting maximum monthly rainfall in Bangladesh. The results of their research indicated that in 36% of the stations, the GEV distribution and in 26% of the stations, the Pearson type 3 and Log-Pearson Type 3 distributions had better fit. The probabilistic distribution analysis of marginal precipitation in Mumbai, India, was implemented by Parchure & Gedam (2019) They considered 26 meteorological stations in their research, and the GEV distribution showed better results in 29% of the stations. Frechet and GP probability distributions also performed better in 27 and 22% of stations, respectively.

Li et al. (2020) increased the temporal coverage of precipitation data for reanalysis using Bayesian probabilistic methods. The results showed that the proposed method has a high performance. Lemus-Canovas et al. (2019) used a spatial regression model to investigate the spatial distribution of precipitation in the eastern Pyrenees. Consequently, the results revealed that there is a correlation between the estimated and observed precipitation with the coefficient of explanation (R² = 0.8).

Khan et al. (2019) studied the spatial distribution of rainfall in the peninsular region of Malaysia from 1951 to 2007. The results showed that the amount of annual precipitation did not change significantly. The trend of regional precipitation had significant changes, which could be due to the increase in design periods with an average of 4.8 and 4.9 days per decade. Iqbal et al. (2019) studied the spatial distribution of precipitation and rainfall trends from 1951 to 2007 in the Himalayas in Pakistan. The results indicated changes in precipitation in summer between 0.25 and 1.25 mm per year. Also, the number of marginal rainfalls and dry days increased. Their research pointed out that increases in these two parameters have led to increased flood and drought.

Vélez et al. (2019) studied the spatial distribution of daily rainfall concentrations at 20 Puerto Rican stations from 1971 to 2010. The results showed that the index value of annual and daily rainfall concentration varies from southeast to northwest.

In general, the Semnan province has a great variety of climates. Also, according to the studies, forecasting the temporal and spatial distribution of precipitation is of great importance in the implementation of development projects and the development of economic infrastructure. On the other hand, the study of the temporal and spatial distribution of precipitation requires a large amount of observational data, making such research works difficult in the region of the Semnan province. In such cases, the use of probabilistic distributions can help study temporal precipitation changes. Therefore, in the present study, the temporal and spatial distribution of P_year and P_max24h in the Semnan province have been analyzed. For this purpose, probabilistic distributions and different spatial interpolation methods have been used. Then, the amount of P_year and P_max24h and its spatial distribution with 10-, 25-, 50-, 100-, and 200-year return periods were predicted. Also, the relationship between P_year and P_max24h was estimated in different return periods, which is helpful when the P_max24h data were missed.

In the present research, precipitation information from Semnan, Shahroud, Damghan, Garmsar, Shahmirzad, and Meyami stations has been used to study the temporal and spatial variations of rainfall in the Semnan province. The rainfall data period of Semnan, Shahroud, Damghan, and Garmsar is from 1986 to 2018. This period for Shahmirzad and Meyami stations is from 2001 to 2018 (Table 1).

A time series is predictable if it has long-term memory. This approach is based on whether all the expected phenomena have occurred in the available time series or not (Karamouz et al. 2012).

Hurst (1951), by studying the water level of the Nile River, presented a test for extreme events, now referred to as the Hurst index. The steps of this method are as follows:

At the beginning, the data scale is normalized using Equation (1). By considering a runoff time series of we have:

HEC-SSP software was introduced in 2008 by the U.S. Army Corps of Engineers for statistical analysis of hydrological data (Brunner & Fleming 2010). This software has different sections, such as data definition and analysis (Harris et al. 2008; Root & Papakos 2010). Primarily, in the present study, the P_year and P_max24h were entered in the data definition section. Then, in the analysis and distribution fitting analysis section, different statistical distributions were fitted to the data. In the end, the best statistical distribution was selected, and the amount of precipitation was measured with different return periods for each station. It is necessary to mention that P_year is calculated based on the summing of daily rainfall in a year. Also, the P_max24h was the maximum daily rainfall in a year.

Interpolation is for finding unknown values of some points based on the known values of other points. The most important application of these methods is estimating the spatial distribution based on the discrete points. Kriging, IDW, and RBF methods are some of the most potent interpolation methods used in many studies (Agung Setianto & Tamia Triandini 2013; Gong et al. 2014; Arifin et al. 2015; Ikechukwu et al. 2017; Liu et al. 2020). Therefore in the present study, the mentioned methods were used for spatial interpolation of P_year and P_max24h in the Semnan province. In this study, the cross-validation (CV) method was used. The CV method is best known for widespread use in interpolation applications (Hancock & Hutchinson 2002). The main benefit of the CV approach is that it is a clearly defined and user-independent process. The CV approach is not acceptable for surfaces that have an insufficient number of observed input points (Jeffrey et al. 2001).

In the above equation, Z(s_i) and Z(s_i + h) are values that are at a distance h from each other, and is a set of values of observational distances. Theoretically, the observed values and the values calculated by the Kriging method should be the same. But in reality this does not happen due to the presence of a nugget effect (indicating measurement error or changes in the microscale). Therefore, the values calculated by the Kriging method will largely depend on the appropriate choice of variogram model (Gong et al. 2014). The most important main feature of the Kriging method is that the spatial smoothing is estimated through variograms. Also, the uncertainty estimates for the predicted values are presented as variance. A semi-variogram model is required to determine the expectations at which Kriging variance is calculated.The necessary calculations in this method can be performed in a defined local neighborhood (Goovaerts 1997), but we used global Kriging in this work, using all observations for spatial prediction.

SURFER (Keckler 1994) is very suitable and flexible software for drawing two-dimensional and three-dimensional models. This software has several versions. The latest version is a combination of CAD and SURFER 8 software. This software is very useful in drawing two-dimensional and three-dimensional isometric maps. Three important applications of this software in natural geography are as follows: (1) drawing topographic maps in different scales, (2) drawing isometric weather maps (isothermal, isobars, isohyetal), and (3) designing a three-dimensional model of topographic and geological maps in different scales.

Genetic algorithm (GA) is an evolutionary algorithm for optimization based on extensive and effective search in large spaces based on genes and chromosomes (Ehteram et al. 2018). This search consists of three steps. In the first stage, the initial population, consisting of a set of chromosomes, is formed. In the second stage, the value of each member is measured using the definition of the objective function. In the final step, genetic operators produce new members, which include producing offspring from selected parents and mutating the members, and finally, gradual evolution is performed (Goldberg 1989). The selection stage is based on the fitting degree of members and some of the most suitable chromosomes are selected for reproduction. Finally, genetic operators are implemented on members, and their genetic codes are modified and synthesized (Zhang et al. 2012). In this study, for the GA, the number of family members was 300, the number of offspring was 240, the number of mutant members was 90, the mutation rate was 0.04, and the number of iterative steps was 200. To select parents from family members, the Roulette Wheel Selection method was used.

Table 3 shows the results of the study of the adequacy of P_year and P_max24h based on the Hurst method. According to this table, all the data used in the present study in all stations are adequate.

The results of the goodness fit test for statistical distributions for the P_year of the studied stations in the Semnan province are presented in Table 4. According to this table, GP distribution for Semnan, Damghan, Garmsar, and Meyami stations; Gumbel distribution for the Shahroud station; and exponential distribution for the Shahmirzad station had the best fitting results. 66.66% of the studied stations followed the GP distribution, and about 16.67% followed the Gumbel distribution and Exponential distribution. This result may be due to skewness in the distribution of P_year at the studied stations.

Table 5 shows the goodness fit test results of statistical distributions for P_max24h in the studied stations. According to Table 5, in 50% of the stations, the best distribution is related to the GP distribution. About 16.67% of stations follow Log normal, Exponential distribution, and Value Generalized Extreme distributions. Table 6 shows the goodness fit test outcomes of statistical distributions for the ratio of P_max24h to P_year in the studied stations. According to Table 6, in 50% of the stations, the best distribution is related to the GEV distribution. About 16.67% of stations follow normal, triangular, and GP distribution. In Tables 4–6, the best Distribution has been highlighted, and other statistical distributions at lower ranks have been shown.

To evaluate the spatial distribution methods, primarily, estimated data were calculated for each of the observed data using the cross-evaluation method.

In other words, in all return periods of 10, 25, 50, 100 and 200, the amount of annual rainfall and 2-h rainfall decreases from the north and northeast of Semnan province to the southeast. The amount of error of each method is shown using evaluation criteria based on Table 7.

The Kavir plain is located in the south of the Semnan province and there is no city, and therefore, no synoptic station and almost all cities of this province are located in the northern half of the province. For instance, Semnan city, which is the seat of the Semnan province, and Garmsar city is closer to the southern desert of the province and have an arid and warm climate. However, Shahmirzad city has higher rainfall than other regions of the province due to locating on the Alborz mountain range heights. Shahmirzad is only 25 km away from Semnan city, but they are different climate conditions. In this study, stations were studied whose data were complete. Also, in this study, we wanted to examine only changes in latitude and longitude.

Despite the many advantages, the existing problems in fitting the models are among the disadvantages of Kriging methods (Teymourzade et al. 2019). Although the IDW method is quick and easy, it does not provide information on dependency or estimation of spatial uncertainty (Dingman 2002).

The climatic conditions of Iran have a great variety that has been created due to changes in latitude and altitude. Semnan province has a variety of climatic conditions in different stations, and these conditions have prevailed on the stations due to the locations of the stations. Statistical analyses of P_year and P_max24h and the ratio of P_max24h to P_year were performed in six stations in Semnan, Shahroud, Damghan, Garmsar, Shahmirzad, and Meyami. Based on the P_year analysis, the Shahmirzad station has the highest annual average and the lowest average P_year in the Damghan station. Then, spatial interpolation using Kriging and IDW methods and 10-, 25-, 50-, 100-, and 200-year return periods were studied. In general, both Kriging and IDW methods had good performance in spatial distributions. For all return periods, the P_year and P_max24h decreases to the southeast of Semnan province. At the same time, the closer we get to the north of the region, especially to the Shahmirzad station, the P_year increases to a P_max24h.

There is a good relationship between P_year and P_max24h. The values of R² coefficient in the 10-, 25-, 50-, 100-, and 200-year return periods were 0.988, 0.992, 0.988, 0.983 and 0.98, respectively. In the temporal analysis, the GP distribution with about 58% in the studied stations was the selected distribution in the analysis of P_year and P_max24h. Half of the regions followed the GEV distribution for the ratio of a P_max24h to P_year. Based on the results obtained in the present study, it is suggested that researchers can use statistical methods as an optimal approach for locating other meteorological stations. One of the limitations of this study is the small number of stations that have complete data. New methods can also be used to predict time with different return periods.

This article is not funded by any scientific institution.

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.