Engineering, Technology & Applied Science Research Vol. 8, No. 4, 2018, 3213-3217 3213 www.etasr.com Laghari et al.: Analysis of Mapping Techniques for Mountain Precipitation: A Case Study … Analysis of Mapping Techniques for Mountain Precipitation: A Case Study of Alpine Region, Austria Abdul Nasir Laghari Department of Energy and Environment Engineering, Quaid-e-Awam University of Engineering, Science and Technology, Nawabshah, Pakistan a.n.laghari@quest.edu.pk Gordhan Das Walasai Department of Mechanical Engineering, Quaid-e-Awam University of Engineering, Science and Technology, Nawabshah, Pakistan valasai@quest.edu.pk Daddan Khan Bangwar Department of Civil Engineering, Quaid- e-Awam University of Engineering, Science and Technology, Nawabshah, Pakistan skb_khan2000@yahoo.com Aftab Hameed Memon Department of Civil Engineering, Quaid-e-Awam University of Engineering, Science and Technology, Nawabshah, Pakistan aftabm78@hotmail.com Abdul Hannan Shaikh Department of Mathematics, Quaid e Awam University of Engineering, Science, and Technology, Nawabshah, Pakistan hanangul12@yahoo.co.uk Abstract—Truly representative precipitation map generation of mountain regions is a difficult task. Due to poor gauge representativity, complex topography and uneven density factors make the generation of representative precipitation maps a very difficult task. To generate representative precipitation maps, this study focused on analyzing four different mapping techniques: ordinary kriging, spline technique (SP), inverse distance weighting (IDW) and regression kriging (RK). The generated maps are assessed through cross-validation statistics, spatial cross-consistency test and by water balance approach. The largest prediction error is produced by techniques missing information on co-variables. The ME and RMSE values show that IDW and SP are the most biased techniques. The RK technique produced the best model results with 1.38mm and 72.36mm ME and RMSE values respectively. The comparative analysis proves that RK model can produce reasonably accurate values at poorly gauged areas, where geographical information compensated the poor availability of local data. Keywords-mountain regions; poor gauge representativity; spatial interpolation techniques I. INTRODUCTION Scientists generally agree that the earth is undergoing critical climate changes [1-5]. Climatic changes based on the hydrology, development and management of water resources have been under major attention over the years. The distributed hydrological models are gaining huge standing in analyzing and investigating the overall impacts on mountain regions and their environment [6-8]. Distributed hydrological models require input variables like estimates of climatic variables at regular and continuous intervals as pre-requisites for their proper functioning [5, 6, 9, 10]. The amount of rainfall is the most vital parameter for any distributed hydrological model. Nevertheless, the amount of rainfall is a matter of various uncertainties like measurement errors, systematic errors during applying interpolation and stochastic errors resulting from the random nature of rainfall. The performance of the models depends heavily upon the accurate estimation of precipitation over specific area and time. The results can be highly compromised [11]. This challenge of reliable and accurate rainfall estimates increases in mountain regions where the geography is complex and measuring stations are scattered over vast areas and concentrated in the valleys [12-14]. The measurement of accurate data for the mountainous range is a very difficult task, resulting in poor representation in the model that analyzes the various rainfall patterns. Resultantly, in these types of situations, when no single method is optimal nor the accuracy of a specific interpolation technique is proven, the performance relies on the variable under study, spatial configuration, and the assumptions used in the estimation [15- 17]. The accuracy of measured data under a certain technique can be verified by comparing and analyzing by applying different techniques to the same data. In order to achieve this objective, the current study analyzes the range of stochastic and deterministic mapping techniques to estimate the values at un- gauged locations. II. DATA AND MODELS A. Study Area The study was carried out at the Alpine area of Kitzbühl Ache region situated in the Austrian Eastern Province of Tyrol with an area of about 2000km2. Complete details about the study region are given in [18]. The available 30 year time series (1960-1990) of mean annual precipitation of 14 gauge stations are taken into consideration for spatial interpolation. The catchment topography is highly rugged with elevation that ranges from 400m to 2400m above m.s.l as shown in Figure 1. The catchment shows strong seasonal precipitation behaviors. B. sto dis pre kri spa is [17 [20 wh C. cli of est dis wh ob loc po sig ass inf far D. dim tec of cur the cal cur ten Engineerin www.etasr Spatial Inter Mapping tec ochastic appro stance weighin ecipitation, w iging [19-21]. atial correlatio able to remov 7, 26]. This a 0, 26-28]. Th hereas the deta Fig. 1. S IDW Techni It is form o matic data by each unknown timated have d stance from kn ˆ( ) i T X   1 ( , ) i pd x x   i xd S ,( 1 here )(ˆ XT = in served values cation, d(x, xi) int, and P= gnificantly inf signs more w formation, wh raway samples SP Techniqu It is a de mensional cu chnique assign neighboring rvature. This p e sample point lculated samp rvature. The te nsion. The firs ng, Technology r.com rpolation Tech chniques are c oaches. The fo ng and spline whereas, the la The regressio on of the differ ve the drifts w approach has hese techniqu ailed discussio Study area with lo ique of a determin y averaging sa n point [31-33 decreasing inf nown locations * ( )i iT X .p S p ix ) 1 , follows nterpolated val at the ith locat i)= distance f weighting po fluence the in weightage to c hile lower po s with little inf ue eterministic te urves on 3-d ns the mathem points that produces an ev ts. This is alik ple points whi echniques com st one generat y & Applied Sci hniques categorized int ormer techniq and is used ater technique on is carried rent variables with validated been used fo ues are brief n is referred to ocations of 14 rain nistic techniqu ample values i 3]. It’s assume fluence or wei s. The formula 1i i  . lue at location tion, T(Xi)= ob from known p ower. The cho nterpolation re loser points a ower assigns m formation. echnique, wh dimensional s matical function t decrease t ven surface w ke to fit a rubb ilst decreasing mmonly used tes the uniform ience Research Laghari et al.: to deterministi que includes in for interpolati e includes ord out on the ba [22-25]. The m intrinsic hypo or climatic ma fly presented o [20, 29, 30]. n gauge stations ue that interp in the neighbo ed that the var ight with incre a is given as: (1) (2) n x, λi= the wei bserved value point to interp oice of powe esults. High p and results in more weighta hich represen surfaces [34]. ns to certain nu the whole su while passing a ber piece acro g the entire su are regularize m changing su h V : Analysis of M ic and nverse ion of dinary asis of model othesis apping here, polates orhood riables easing ight of s at ith preted er can power more age to nts 2- . The umber urface around oss the urface d, and urface with whi valu insp regu foll f whe unk solu kno fun coe E. and neig bas of k for by D whe obs valu the loca arou esti kno wei by: for sem sele (T̂ equ Vol. 8, No. 4, 20 Mapping Techniq h values that ile the other o ues confined pection and r ularized splin lowing formul ˆ( ) ( ,iT X f x { 21 ( ) 2 4 2 0 r R r r K              f(xi, yi)=a1+a2 ere )(ˆ XT = in known smooth ution of a syst own point to nction, C= the efficients found OK Techniqu This techniqu d works the s ghboring mea ed on distance known points. kriging, after D.L. Krige: 0 1 ˆ( ) n i T X    ere )(ˆ 0XT = th served value a ues at the ith l number of ation and spat und the pre imator varianc own as “best ights sum to u 1 ( , n i i j i X X    all j, where φ mivariogram a ected as a bes )( 0X , and t uation: bAc .1 , 018, 3213-3217 ques for Mount may remain one influences by the sam results through e method was a is used for su 1 , ) ( N i i i y R r   } 2 ln 2 ln 2 r c r c                   2x+a3y nterpolated spl h function, λ tem of linear interpreted p e constant eq d by a solution ue ue uses a semi ame way as I asured values. es but also dep . Author in [3 a first use of t . ( )i iT X he interpolated at ith location, location. The measured poi tial relationsh ediction locat ce and ensure t linear unbi unity [36]. Est ) ( ,j jX   φ=Lagrangian analysιs, the st-fit model fr the weights λi c 7 tain Precipitati outside the s s the rigidity o mple data ra h cross-valida s selected from urface interpo )ir } 1   line value at lo i= the coeffic equations, ri= point, K0= the qual to 0.577 n of a system o ivariogram for IDW works, However, we pend on total s 35] developed the theory in t d value at loc and λi= the w λi depends on ints, distance hips among th tions. Krigin es the unbiase iased estimato timator varianc 0, )X n multiplicator Gaussian va from all other i are calculate        i  3214 ion: A Case Stu sample data r of the surface ange. After v ation statistics m both types lation: (3) (4) (5) ocation x, f(xi, cients found = the distance e modified B 7215, and αi= of linear equat r spatial predi which weight eights are not spatial arrange d a general for the mining ind (6) cation x, T(Xi) weight of obse n: Fitting mod to the predi he measured v ng minimizes ed estimation or (BLUE)” ce can be ach (7) r. After perfor ariogram mod models to pr ed from the m (8) udy … range, with visual s, the . The , yi)= by a from Bessel = the tions. iction ts the only ement rmula dustry = the erved del to iction values s the often with ieved rming del is redict matrix Engineering, Technology & Applied Science Research Vol. 8, No. 4, 2018, 3213-3217 3215 www.etasr.com Laghari et al.: Analysis of Mapping Techniques for Mountain Precipitation: A Case Study … where matrix A contains the semivariances of all data point pairs, b is a vector containing semivariances between the location of interest and observed point. λi is the weight to be calculated [20, 37-41]. F. RK Technique Authors in [17] conclude that ordinary kriging does not produce representative precipitation values in mountain regions all the time. Alternative techniques might fully utilize the relationship between predicted variable and co-variables for variability analysis. Two of these techniques are cokriging and regression kriging. Due to poor cross-covariance between precipitation and any topographical variable, the former application did not produce good results. However, in such situations regression kriging seemed the natural choice, which is commonly used in hydro-sciences [42, 43]. The technique defines the relationship between target and co-variables in order to predict the values at grid nodes through linear regression. The auxiliary variables are easy to measure, they provide an alternative to target variable at the under sampled locations to model and quantify the existing patterns. To quantify the existing trend of the variable and its variability in a regression model, we preferred to use multiple linear regressions - the further extension of straightforward linear regression with the variety of descriptive variables. As altitude alone flopped to signify the variation of precipitation sufficiently, we tested other variables like slope, topographic index, aspect, hill shade, curvature etc. These extra explanatory variables were obtained from the elevation model, in order to make better predictions of the target variables at grid nodes of the DEM. We adopted stepwise procedures to select the most crucial variables and the subsequent regression equation to predict the target variable at un-sampled locations. The regression led to a three parameter equation, significant at 5% level, explaining 62% of variability of precipitation.      1 ˆ : 1075.94 42.59 * 0.26 * 3.98 * R Hill shade DEM Slope T      (9) The resulting regression residuals (ε) are further kriged at grid nodes by fitting variogram models, and then finally both values are summed up to predict the target variable values. ˆ ˆ ˆ( ) ( ) ( )RK R OKT X T X X  (10) where RKT̂ = Combined predicted values of target variable, RT̂ = Regression predicted values of target variable, OK̂ = Kriged values of regression residuals. III. VALIDATION METHODS The mapping techniques performance was assessed by following three steps: A. Cross-Validation Statistics The cross validation method depends on eliminating one sample location (measurement station) from the data set at a time and calculating the value of the removed sample with the remaining data points. This routine was followed for each measurement station. The comparative indices were then used as a measure of prediction quality by the ME and the RMSE, which are defined as follows: 1 1 ˆ( ) n i i i ME T T n    (11) 0.5 2 1 1 ˆ( ) n i i i RMSE T T n          (12) where n is the number of validation points, and ˆ &T Ti i are the predicted and observed values at location i. The ME criterion is used to check the conditional bias property, while the RMSE criterion assesses the precision quality. A smaller value of RMSE indicates higher accuracy and vice versa. Cross validation statistics can be used to find the optimal mapping technique, however, the presence of short range correlations in data may raise questions regarding the reliability of its statistical results [35]. B. Cross-Consistency A second step to analyze reliability and consistency of predictions, spatial cross-consistency approach was adopted [44]. All statistical parameters of different calculated precipitation mapping estimates were compared with a referenced precipitation map (RPM). This RPM was carefully produced during a 4 year project (www.waterpool.org) in which different experts from different institutes were involved, and results were consistent with water balance estimates. C. Water Balance Approach Finally all calculated precipitation maps including the RPM were evaluated by means of a general water balance approach.        – – Q discharge P precipitation ET evapotranspiration S Storage changes   (13) Gridded discharges were calculated for each mapping technique as a result of subtraction of actual evapotranspiration grid estimates from interpolated mapping precipitation grid estimates. The actual evapotranspiration values are obtained from the hydrological Atlas of Austria. Storage changes can be ignored, as for long-range mean annual water balances; it was assumed that there is no sensible change in the water contents of different reservoirs, e.g. groundwater, snow cover [44, 45]. The difference between the calculated discharges with observed discharges gives a measure for the reliability and consistency of the precipitation. IV. RESULTS AND DISCUSSION The ME values clearly indicate the superiority of RK technique over the other techniques, showing almost 42% less bias than OK. Whereas IDW and SP techniques produce higher bias, almost 2 to 5 times higher than OK technique. The RMSE results also reveal the primacy of RK technique over the other ones. However, all techniques yield high uncertainty in calculated values. Thus, considering co-variables into account certainly improves the performance by decreasing RMSE values from 111 to 72. The spatial cross-consistency tests are also conducted by computing precipitation maps with ref ado Fig stud dif pro inc per ref the ele po sta of cov zon Fig com reg map Fig com Neg valu ann Engineerin www.etasr ferenced preci opted in [46]. g. 2. Cross va dy region are the The statisti fferences. Ho omising resul corporation o rcentile varia ferenced map. eir performanc evation zones sition of gaug ations are locat basin area. Zo vering 39% a ne-4 possesses g. 3. Summary mputed precipitat gion. Negative val p value, while po g. 4. Compara mputed precipitati gative values ind ues indicate the o The basin le nual values in ng, Technology r.com ipitation raster The comparis alidation statistics ME and RMSE ical paramete owever, the R lts except in of topographi ation in RK To further an ces, first at 4 c of the whol ge stations are ted below 100 ones 2 and 3 c and 26% of a s only 3% area y statistics of p tion estimates by lues indicate a lo ositive values indi ative statistics of ion estimates by dicate the less val opposite. evel values fo n the region. y & Applied Sci r. The same s on results are s of four mappin ers clearly RK technique n maximum v ical informati map is clos nalyze all tech catchment leve e basin (Figu e highly biase 00m (Zone-1), ontain only on area respectiv a without any percent differenc IDW, Spline, OK ower value with r cate the opposite. f percent differen IDW, Spline, OK lue from the refer ollow the sam RK method ience Research Laghari et al.: strategy was e given in Figur ng techniques app show the v e yields the values due t ion. Similarly sely matched hniques we ana els and then th ures 3 and 4) ed, more than covering only ne station each vely. The top gauge station. ces between RP K and RK for th respect to the refe . nces between RP K and RK at four renced map and p me pattern of is found to p h V : Analysis of M earlier re 2. plied in visible most to the y the d with alyzed hrough ). The n 80% y 32% h, with most PM and e entire ferenced PM and basins. positive mean present bett othe tech just con altit OK altit How not and high perf dep zon com The imp map was ord esti atla map gau with The tech Fig. com elev map pro valu at h uno pro leve [1] [2] Vol. 8, No. 4, 20 Mapping Techniq ter result in th er techniques hniques like tifies the co nclusions are s tude. Without K techniques p tude which co wever, at the be effective, d 33% at zon her sparse dat formance of pendence resul nes 2 and 3 mparatively, al e better resu portance of g pping at highe s adopted for der to validate imates of the w as of Austria. ps and evapotr uged runoff an h RPM, 2.7% ese outcomes hnique in the h 5. Comparati mputed precipitati vation zones. Neg p and positive valu Comprehensiv duces reason ues. The prim high altitude sp obtainable. Ho duce good re el height with B. Klove, P. A Kværner, T. M Velasco, M. groundwater an 518B, pp. 250-2 H.-M. Fussel, A A. Bigano, S. C 018, 3213-3217 ques for Mount hree basins, w have better es IDW, spline onclusions of strengthened b t elevation in performed com ontains 86% o higher altitud where the con e-4 resulting ta zones. Duri OK at high lted to be hig 3. For analy ll stations wer ults from RK geographical er sparse data z the four inte e the RPM. M whole region w Total runoff w ranspiration es nd computed ru % with RK and s validate th higher altitude ive statistics of ion estimates by gative values indi ues indicate vthe V. CO ve analysis c nably accurate macy of this te parse data zon owever, the sp esults over ge a wide networ REFER Ala-Aho, G. Bertr Muotka, H. Mykr Pulido-Velazqu nd dependent eco 266, 2014 A. Jol, A. Marx, M Castellari, M. Erh 7 tain Precipitati whereas in Kit stimates. The s over geosta f many othe y further analy nformation the mparatively w of the total sta des these three nsistency rang in relegation ing semivariog zone was po gher at zone 1 yzing the ma re included in K at higher information zones. A wate erpolation pre Mean annual were taken fro was computed stimates. The unoff found o d 13-20% with he overall su e region with sp percent differenc y IDW, Spline, icate the less val opposite. ONCLUSION learly indicat e average an echnique is pa nes where map line and IDW eostatistical te rk of gauge sta RENCES rand, J. J. Gurdak ra, E. Preda, P. uez, “Climate osystems”, Journ M. Hilden, A. Apa ard, B. Georgi, C 3216 ion: A Case Stu tzbüleher Ach superiority of atistical techn er studies. T ysis on the ba e IDW, spline well below 10 ations under s e techniques c ges at 8% at zo of their abili gram modelin oor, as the sp 1 and weaken apping techn the final map zones prove in estimating er balance appr ecipitation ma evapotranspir om the hydrolo d from precipit difference bet out to be only h other techni uperiority of parse data. ces between RPM OK and RK a lue from the refe es that RK m nnual precipit articularly obse pping variable W model results echniques at a ations. k, H. Kupfersber Rossi, C. B. U change impact nal of Hydrology aricio, A. Bastrup Climate change, im udy … he the usual niques These sis of e and 000m study. could one-2 ity at ng the patial ned in niques pping. e the g the roach aps in ration ogical tation tween 0.4% iques. f RK M and at four renced model tation erved es are s also a low rger, J. Uvo, E. ts on y, Vol. p-Birk, mpacts Engineering, Technology & Applied Science Research Vol. 8, No. 4, 2018, 3213-3217 3217 www.etasr.com Laghari et al.: Analysis of Mapping Techniques for Mountain Precipitation: A Case Study … and vulnerability in Europe 2016-An indicator-based report, EU Publications, 2017 [3] M. Burke, J. Dykema, D. B. Lobell, E. Miguel, S. Satyanath, “Incorporating climate uncertainty into estimates of climate change impacts”, Review of Economics and Statistics, Vol. 97, No. 2, pp. 461- 471, 2015 [4] C. B. Field, V. R. Barros, K. Mach, M. Mastrandrea, Climate change 2014: impacts, adaptation, and vulnerability. Vol. 1, Cambridge University Press, 2014 [5] A. Laghari, D. Vanham, W. Rauch, “To what extent does climate change result in a shift in Alpine hydrology? A case study in the Austrian Alps”, Hydrological Sciences Journal, Vol. 57, No. 1, pp. 103-117, 2012 [6] D. Viviroli, M. Zappa, J. Gurtz, R. Weingartner, “An introduction to the hydrological modelling system PREVAH and its pre-and post- processing-tools”, Environmental Modelling & Software, Vol. 24, No. 10, pp. 1209-1222, 2009 [7] R. Cibin, P. Athira, K. Sudheer, I. Chaubey, “Application of distributed hydrological models for predictions in ungauged basins: a method to quantify predictive uncertainty”, Hydrological Processes, Vol. 28, No. 4, pp. 2033-2045, 2014 [8] Y. Chen, J. Li, H. Xu, “Improving flood forecasting capability of physically based distributed hydrological models by parameter optimization”, Hydrology and Earth System Sciences, Vol. 20, No. 1, pp. 375-392, 2016 [9] K. X. Yu, L. Gottschalk, L. Xiong, Z. Li, Li, “Estimation of the annual runoff distribution from moments of climatic variables”, Journal of Hydrology, Vol. 531, pp. 1081-1094, 2015 [10] Z. K. Tesemma, Y. Wei, M. C. Peel, A. W. Western, “Including the dynamic relationship between climatic variables and leaf area index in a hydrological model to improve streamflow prediction under a changing climate”, Hydrology and Earth System Sciences, Vol. 19, No. 6, pp. 2821-2836, 2015 [11] K. Beven, “How far can we go in distributed hydrological modelling?”, Hydrology and Earth System Sciences, Vol. 5, No. 1, pp. 1-12, 2001 [12] E. Lepuschitz, “Geographic information systems in mountain risk and disaster management”, Applied Geography, Vol. 63, pp. 212-219, 2015 [13] S. Takaoka, “Origin and geographical characteristics of ponds in a high mountain region of central Japan”, Limnology, Vol. 16, No. 2, pp. 103- 112, 2015 [14] N. Boers, B. Bookhagen, N. Marwan, J. Kurths, “Spatiotemporal characteristics and synchronization of extreme rainfall in South America with focus on the Andes Mountain range”, Climate Dynamics, Vol. 46, No. 1-2, pp. 601-617, 2016 [15] J. Creutin, C. Obled, “Objective analyses and mapping techniques for rainfall fields: an objective comparison”, Water Resources Research, Vol. 18, No. 2, pp. 413-431, 1982 [16] D. D. Weber, E. J. Englund, “Evaluation and comparison of spatial interpolators II”, Mathematical Geology, Vol. 26, No. 5, pp. 589-603, 1994 [17] C. Prudhomme, D. W. Reed, “Mapping extreme rainfall in a mountainous region using geostatistical techniques: a case study in Scotland”, International Journal of Climatology, Vol. 19, No. 12, pp. 1337-1356, 1999 [18] D. Vanham, E. Fleischhacker, W. Rauch, “Technical Note: Seasonality in alpine water resources management? a regional assessment”, Hydrology and Earth System Sciences, Vol. 12, No. 1, pp. 91-100, 2008 [19] M. Knotters, D. Brus, J. O. Voshaar, “A comparison of kriging, co- kriging and kriging combined with regression for spatial interpolation of horizon depth with censored observations”, Geoderma, Vol. 67, No. 3-4, pp. 227-246, 1995 [20] I. O. A. Odeha, A. B. McBratney, D. Chittleborough, “Further results on prediction of soil properties from terrain attributes: heterotopic cokriging and regression-kriging”, Geoderma, Vol. 67, No. 3-4, pp. 215-226, 1995 [21] I. O. A. Odeha, A. B. McBratney, D. J. Chittleborough, “Spatial prediction of soil properties from landform attributes derived from a digital elevation model”, Geoderma, Vol. 63, No. 3-4, pp. 197-214, 1994 [22] V. Chaplot, C. Walter, P. Curmi, “Improving soil hydromorphy prediction according to DEM resolution and available pedological data”, Geoderma, Vol. 97, No. 3-4, pp. 405-422, 2000 [23] I. D. Moore, P. E. Gessler, G. A. E. Nielsen, G. A. Peterson, “Soil attribute prediction using terrain analysis”, Soil Science Society of America Journal, Vol. 57, No. 2, pp. 443-452, 1993 [24] J. L. Richardson, W. J. Edmonds, “Linear regression estimations of Jenny's relative effectiveness of state factors equation”, Soil Science, Vol. 144, No. 3, pp. 203-208, 1987 [25] J. A. Thompson, J. C. Bell, C. A. Butler, “Digital elevation model resolution: effects on terrain attribute calculation and quantitative soil- landscape modeling”, Geoderma, Vol. 100, No. 1-2, pp. 67-89, 2001 [26] M. R. Holdaway, “Spatial modeling and interpolation of monthly temperature using kriging”, Climate Research, Vol. 6, pp. 215-225, 1996 [27] A. Martinez-Cob, “Multivariate geostatistical analysis of evapotranspiration and precipitation in mountainous terrain”, Journal of Hydrology, Vol. 174, No. 1-2, pp. 19-35, 1996 [28] D. L. Phillips, J. Dolph, D. Marks, “A comparison of geostatistical procedures for spatial analysis of precipitation in mountainous terrain”, Agricultural and forest meteorology, Vol. 58, No. 1-2, pp. 119-141, 1992 [29] P. Goovaerts, Geostatistics for natural resources evaluation. Oxford University Press, 1997 [30] T. Hengl, G. B. Heuvelink, A. Stein, Comparison of kriging with external drift and regression kriging, ITC, 2003 [31] P. A. Burrough, R. A. McDonnell, C. D. Lloyd, Principles of geographical information systems, Oxford University Press, 2015 [32] D. R. Legates, C. J. Willmott, “Mean seasonal and spatial variability in global surface air temperature”, Theoretical and Applied Climatology, Vol. 41, No. 1-2, pp. 11-21, 1990 [33] C. Stallings, R. Huffman, S. Khorram, Z. Guo, Linking gleams and GIS, Paper-American Society of Agricultural Engineers (USA), 1992 [34] M. Hutchinson, P. Gessler, “Splines—more than just a smooth interpolator”, Geoderma, Vol. 62, No. 1-3, pp. 45-67, 1994 [35] G. Matheron, Le krigeage universel: cahiers du Centre de Morphologie Mathematique, École nationale supérieure des mines de Paris, 1969 (in French) [36] R. Webster, M. A. Oliver, Geostatistics for environmental scientists, John Wiley & Sons, 2007 [37] J. C. Davis, Statistics and Data Analysis in Geology, John Wiley & Sons, 1986 [38] C. V. Deutsch, A. G. Journel, GSLIB Geostatistical Software Library and User’s Guide, Oxford University Press, 1992 [39] G. W. Heine, “A controlled study of some two-dimensional interpolation methods”, COGS Computer Contributions, Vol. 3, No. 2, pp. 60-72, 1986 [40] N. S. N. Lam, “Spatial interpolation methods: a review”, The American Cartographer, Vol. 10, No. 2, pp. 129-150, 1983 [41] A. G. Royle, F. L. Clausen, P. Frederiksen, “Practical universal kriging and automatic contouring”, Geoprocessing, Vol. 1, pp. 377-394, 1981 [42] S. Ahmed, G. De Marsily, “Comparison of geostatistical methods for estimating transmissivity using data on transmissivity and specific capacity”, Water Resources Research, Vol. 23, No. 9, pp. 1717-1737, 1987 [43] J. Delhomme, “Kriging in the hydrosciences”, Advances in Water Resources, Vol. 1, pp. 251-266, 1978 [44] J. Hofierka, J. Parajka, H. Mitasova, L. Mitas, “Multivariate interpolation of precipitation using regularized spline with tension”, Transactions in GIS, Vol. 6, No. 2, pp. 135-150, 2002 [45] M. Kuhn, H. Escher-Vetter, “Die Reaktion der österreichischen Gletscher und ihres Abflusses auf Änderungen von Temperatur und Niederschlag”, Osterreichische Wasser-und Abfallwirtschaft, Vol. 56, No. 1-2, pp. 10-16, 2004 [46] S. G. Custer, P. Fames, J. P. Wilson, R. D. Snyder, “A Comparison of Hand and Spline-Drawn Precipitation Maps for mantainous Montana”, Journal of the American Water Resources Association, Vol. 32, No. 2, pp. 393-405, 1996