Received: 13 January 2016
Accepted: 1 March 2017
DOI 10.1002/hyp.11163
RESEARCH ARTICLE
Cokriging for enhanced spatial interpolation of rainfall in two
Australian catchments
Sajal Kumar Adhikary1,2
|
Nitin Muttil1,2
|
Abdullah Gokhan Yilmaz3
1
College of Engineering and Science, Victoria
University, P.O. Box 14428, Melbourne,
Victoria 8001, Australia
2
Institute for Sustainability and Innovation,
Victoria University, P.O. Box 14428,
Melbourne, Victoria 8001, Australia
3
College of Engineering, University of Sharjah,
P.O. Box 27272, Sharjah, United Arab Emirates
Correspondence
Sajal Kumar Adhikary, College of Engineering
and Science, Victoria University, Footscray
Park Campus, P.O. Box 14428, Melbourne,
Victoria 8001, Australia.
Email: sajal.adhikary@live.vu.edu.au
Abstract
Rainfall data in continuous space provide an essential input for most hydrological and water
resources planning studies. Spatial distribution of rainfall is usually estimated using ground‐based
point rainfall data from sparsely positioned rain‐gauge stations in a rain‐gauge network. Kriging
has become a widely used interpolation method to estimate the spatial distribution of climate variables including rainfall. The objective of this study is to evaluate three geostatistical (ordinary
kriging [OK], ordinary cokriging [OCK], kriging with an external drift [KED]), and two deterministic
(inverse distance weighting, radial basis function) interpolation methods for enhanced spatial
interpolation of monthly rainfall in the Middle Yarra River catchment and the Ovens River catchment in Victoria, Australia. Historical rainfall records from existing rain‐gauge stations of the
catchments during 1980–2012 period are used for the analysis. A digital elevation model of each
catchment is used as the supplementary information in addition to rainfall for the OCK and
kriging with an external drift methods. The prediction performance of the adopted interpolation
methods is assessed through cross‐validation. Results indicate that the geostatistical methods
outperform the deterministic methods for spatial interpolation of rainfall. Results also indicate
that among the geostatistical methods, the OCK method is found to be the best interpolator
for estimating spatial rainfall distribution in both the catchments with the lowest prediction error
between the observed and estimated monthly rainfall. Thus, this study demonstrates that the use
of elevation as an auxiliary variable in addition to rainfall data in the geostatistical framework can
significantly enhance the estimation of rainfall over a catchment.
KEY W ORDS
cross‐variogram model, digital elevation model, kriging with an external drift, ordinary cokriging,
positive‐definite condition, variogram model
1
|
I N T RO D U CT I O N
in the field because the number of stations in a network is often restricted
by economic, logistics, and geological factors (Goovaerts, 2000). As a result,
Rainfall data provide an essential input for many hydrological investiga-
point rainfall data are generally accessible from a limited number of
tions and modelling tasks. Accuracy of various hydrological analyses such
stations. These limitations increase the need for using suitable spatial esti-
as water budget analysis, flood modelling, climate change studies, drought
mation methods to obtain the spatial distribution of rainfall and generate
management, irrigation scheduling, and water management greatly
rainfall map from the point rainfall values. Moreover, the network is often
depends on the correct estimation of the spatial distribution of rainfall
not deployed on a regular grid and rainfall data may not be available in the
(Delbari, Afrasiab, & Jahani, 2013; Moral, 2010). This usually requires a
target location where it is most required (Adhikary et al., 2016a). In such
dense rain‐gauge network with a large number of stations (Adhikary,
cases, spatial interpolation plays a vital role to simulate rainfall in areas with
Muttil, & Yilmaz, 2016a). However, the rain‐gauge network is often sparse
no stations based on the observed rainfall values in the surrounding areas.
This is an open access article under the terms of the Creative Commons Attribution‐NonCommercial‐NoDerivs License, which permits use and distribution in any
medium, provided the original work is properly cited, the use is non‐commercial and no modifications or adaptations are made.
© 2017 The Authors. Hydrological Processes Published by John Wiley & Sons Ltd.
Hydrological Processes. 2017;1–19.
wileyonlinelibrary.com/journal/hyp
1
2
ADHIKARY ET
Several interpolation methods have been frequently employed to
AL.
(e.g., Di Piazza et al., 2011; Feki et al., 2012; Hevesi et al., 1992; Lloyd,
interpolate rainfall data from rain‐gauge stations and produce the spatial
2005; Martínez‐cob, 1996; Moral, 2010; Phillips, Dolph, & Marks,
distribution of rainfall over a catchment. Examples of such methods range
1992; Subyani & Al‐Dakheel, 2009).
from simple conventional (e.g., Thiessen polygons [Thiessen, 1911],
A wide variety of spatial interpolation methods have been devel-
isohyet mapping [ASCE, 1996], simple trend surface interpolation
oped for the interpolation of spatially distributed point rainfall values.
[Gittins, 1968]), and deterministic methods (i.e., inverse distance
However, it is often challenging to distinguish the best interpolation
weighting (IDW; ASCE, 1996), radial basis function (RBF; Di Piazza, Conti,
method to estimate the spatial distribution of rainfall for a particular
Noto, Viola, & Loggia, 2011) to complex stochastic or geostatistical
catchment or study area. The reason is that the performance of an
methods (i.e., ordinary kriging [OK], ordinary cokriging [OCK] and kriging
interpolation method depends on a number of factors such as catch-
with an external drift [KED; Goovaerts, 2000]). Although the conven-
ment size and characteristics, sampling density, sample spatial distribu-
tional and deterministic methods have been improved over time, their
tion, sample clustering, surface type, data variance, grid size or
limitations continue to exist. These limitations have been described elab-
resolution, quality of auxiliary information to be used as well as the
orately in Goovaerts (2000) and Teegavarapu and Chandramouli (2005).
interactions among these factors. Moreover, it is unclear how the afore-
Geostatistical methods have been shown superior to the conven-
mentioned factors affect the performance of the spatial interpolation
tional and deterministic methods for spatial interpolation of rainfall
methods (Dirks, Hay, Stow, & Harris, 1998; Li & Heap, 2011). Hence,
(Goovaerts, 1997; Isaaks & Srivastava, 1989). Several studies have
the best interpolation method for a particular study area is usually
reported that rainfall is generally characterised by a significant spatial
established through the comparative assessment of different interpola-
variation (e.g., Delbari et al., 2013; Lloyd, 2005). This suggests that inter-
tion methods (e.g., Delbari et al., 2013; Dirks et al., 1998; Goovaerts,
polation methods, which are explicitly able to incorporate the spatial
2000; Hsieh, Cheng, Liou, Chou, & Siao, 2006; Mair & Fares, 2011;
variability of rainfall into the estimation process should be employed.
Moral, 2010). The comparison among different interpolation methods
In view of that, kriging has become the most widely used geostatistical
is made through a validation procedure. The interpolation method that
method for spatial interpolation of rainfall. The ability of kriging to
provides better results with lower bias and higher accuracy in rainfall
produce spatial predictions of rainfall has been distinguished in many
estimation is identified as the best interpolation method.
studies (e.g., Adhikary et al., 2016a; Goovaerts, 2000; Jeffrey, Carter,
In the past, many studies have been devoted to the comparison
Moodie, & Beswick, 2001; Lloyd, 2005; Moral, 2010; Yang, Xie, Liu, Ji,
and evaluation of different deterministic and geostatistical interpola-
& Wang, 2015). The major advantage of kriging is that it takes into
tion methods in a range of different regions and climates around the
account the spatial correlation between data points and provides
world. Dirks et al. (1998) compared four spatial interpolation methods
unbiased estimates with a minimum variance. The spatial variability in
using rainfall data from a network of 13 rain‐gauges in Norfolk Island
kriging is quantified by the variogram model that defines the degree
concluding that kriging provided no substantial improvement over
of spatial correlation between the data points (Webster & Oliver, 2007).
any of the simpler Thiessen polygon (TP), IDW, or areal‐mean
Another key advantage of kriging over the conventional and deter-
methods. Goovaerts (2000) employed three multivariate geostatistical
ministic methods is that while providing a measure of prediction stan-
methods (OCK, KED, simple kriging with varying local means [SKVM]),
dard error (also called kriging variance), it is capable of complementing
which incorporate a DEM as secondary variable and three univariate
the sparsely sampled primary variable, such as rainfall by the correlated
methods (OK, TP, and IDW) that do not take into account the elevation
densely sampled secondary variable, such as elevation to improve the
for spatial prediction of monthly and annual rainfall data available at 36
estimation accuracy of primary variable (Goovaerts, 2000; Hevesi,
rain‐gauge stations. The comparison among these methods indicated
Istok, & Flint, 1992). This multivariate extension of kriging is referred
that the three multivariate geostatistical methods gave the lowest
to as the cokriging method. The standard form of cokriging is the
errors in rainfall estimation. Martínez‐cob (1996) compared OK, OCK,
OCK method, which usually reduces the prediction error variance and
and modified residual kriging to interpolate annual rainfall and evapo-
specifically outperforms kriging method if the secondary variable (i.e.,
transpiration in Aragón, Spain. The results indicated that OCK was
elevation) is highly correlated (correlation coefficient higher than .75)
superior for rainfall estimation, reducing estimation error by 18.7%
with the primary variable (i.e., rainfall) and is known at many more
and 24.3% compared with OK and modified residual kriging, respec-
points (Goovaerts, 2000). The KED is another commonly applied
tively. Hsieh et al. (2006) evaluated OK and IDW methods using daily
cokriging method, which can incorporate the exhaustive secondary var-
rainfall records from 20 rain‐gauges to estimate the spatial distribution
iable (i.e., elevation) to give an enhanced estimation of rainfall when
of rainfall in the Shih‐Men Watershed in Taiwan. The results demon-
dealing with a low‐density rain‐gauge network. Thus, cokriging includ-
strated that IDW produced more reasonable representations than
ing the OCK and KED methods has been the increasing preferred
OK. Moral (2010) compared three univariate kriging (simple kriging
geostatistical methods all over the word. As highlighted by Goovaerts
[SK], universal kriging, and OK) with three multivariate kriging methods
(2000) and Feki, Slimani, and Cudennec (2012), rainfall and elevation
(OCK, SKVM, and regression kriging) to interpolate monthly and
tend to be related because of the orographic influence of mountainous
annual rainfall data from 136 rain‐gauges in Extremadura region of
topography. Therefore, topographic information such as the digital ele-
Spain. The results showed that multivariate kriging outperformed uni-
vation model (DEM) can be used as a convenient and valuable source of
variate kriging and among multivariate kriging, SKVM and regression
secondary data for the OCK and KED methods. The efficacy of incorpo-
kriging performed better than OCK.
rating elevation into the interpolation procedure for enhanced estima-
Ly, Charles, and Degré (2011) used IDW, TP, and several kriging
tion of rainfall has been highlighted in many studies across the world
methods to interpolate daily rainfall at a catchment scale in Belgium.
ADHIKARY ET
3
AL.
The results indicated that integrating elevation into KED and OCK did
2003; Schreider, Jakeman, Pittock, & Whetton, 1996; Yu, Cartwright,
not provide improvement in the interpolation accuracy for daily rainfall
Braden, & de Bree, 2013). The Yarra River catchment is a major source
estimation. OK and IDW were considered to be the suitable methods
of water for more than one third of Victoria's population in Australia
as they gave the smallest error for almost all cases. Mair and Fares
(Barua, Muttil, Ng, & Perera, 2012). Although the catchment is not
(2011) compared TP, IDW, OK, linear regression, SKVM to estimate
large with respect to other Australian catchments, it produces the
seasonal rainfall in a mountainous watershed concluding that OK pro-
fourth highest water yield per hectare of the catchment in Victoria
vided the lowest error for nearly all cases. They also found that incor-
(Adhikary et al., 2015). There are seven storage reservoirs in the catch-
porating elevation did not improve the prediction accuracy over OK for
ment that supports about 70% of drinking water supply in Melbourne
the correlation between rainfall and elevation lower than 0.82. Delbari
city (Barua et al., 2012). The Ovens River catchment is another impor-
et al. (2013) used two univariate methods (IDW and OK), and four mul-
tant source of water in northeast Victoria, which forms a part of the
tivariate methods (OCK, KED, SKVM, and linear regression) for map-
Murray‐Darling basin (Yu et al., 2013). The Ovens River is consid-
ping monthly and annual rainfall over the Golestan Province in Iran.
ered one of the most important tributaries of the Murray‐Darling
They reported that KED and OK outperformed all other methods in
Basin due to the availability of sufficient volume of water with accept-
terms of root mean square error (RMSE). Jeffrey et al. (2001) derived a
able quality and its good ecological condition. The average annual flow
comprehensive archive of Australian rainfall and climate data using a
of the river constitutes approximately 7.3% of the total flow of the state
thin plate smoothing spline to interpolate daily climate variables and
of Victoria (Schreider et al., 1996). The catchment contributes approxi-
OK to interpolate daily and monthly rainfall. The aforementioned
mately 14% of Murray‐Darling basin flows in spite of its relatively small
studies on spatial interpolation of rainfall indicate that each method
catchment area of less than 1% of the total Murray‐Darling basin area
has its advantages and disadvantages and thus performs in a dissimilar
(EPA Victoria, 2003; Yu et al., 2013). Thus, both the catchments have sig-
way for different regions. There is no single interpolation method that
nificant contribution towards the sustainable development of Victoria's
can work well everywhere (Daly, 2006). Therefore, the best interpolator
economy. However, high rainfall variation and diverse water use
for a particular study area or catchment should essentially be achieved
activities in these catchments has complicated the water management
through the comparative assessment of different interpolation methods.
tasks, which has further created strong burden on the water managers
To date, many studies have been conducted on spatial interpola-
and policy makers for effective water resources management. Therefore,
tion of rainfall at a regional and national scale in Australia (Gyasi‐Agyei,
enhanced estimation of rainfall and its spatial distribution is important,
2016; Hancock & Hutchinson, 2006; Hutchinson, 1995; Jeffrey et al.,
which could be beneficial for effective water supply and demand man-
2001; Johnson et al., 2016; Jones, Wang, & Fawcett, 2009; Li & Shao,
agement, and sustainable agricultural planning in both the catchments.
2010; Woldemeskel, Sivakumar, & Sharma, 2013; Yang et al., 2015).
The rest of the paper is arranged as follows. First, a brief descrip-
However, none of these studies was conducted at a local or catchment
tion of the study area and data used are presented, which is followed
scale. Likewise, elevation and rainfall relations locally have been rela-
by the detailed description of the methodology adopted in this study.
tively little studied in Australia and no such studies have been under-
The results are summarized next, and finally, the conclusions drawn
taken specifically within the Yarra River catchment and the Ovens
from this study are presented.
River catchment in Victoria, Australia. Sharma and Shakya (2006)
highlighted that any analysis of hydroclimatic variables should be car-
2
STUDY AREA AND DATA USED
|
ried out at the local scale rather than at a large or global scale because
of the variations of hydroclimatic situations from one region to
another. Therefore, the main aim of this study is to assess if relatively
2.1
|
The study area
more complex geostatistical interpolation methods that take into
This study considers two catchments in Victoria as the case study area,
account the elevation and rainfall relation provide any benefits over
which includes the Middle Yarra River catchment and Ovens River
simpler methods for enhanced estimation of rainfall within the Yarra
catchment in south‐eastern Australia. Figure 1 shows the approximate
River catchment and the Ovens River catchment in Victoria, Australia.
location of the case study area. The Yarra River catchment is located in
The specific focus is to evaluate the effectiveness of the cokriging
northeast of Melbourne covering an area of 4,044 km2. The water
methods including OCK and KED that make use of elevation as a sec-
resources management is an important and multifaceted issue in the
ondary variable over those methods including OK, IDW, and RBF that
Yarra River catchment because of its wide range of water uses as well
do not make use of such information to estimate the spatial distribu-
as its downstream user requirements and environmental flow provi-
tion of rainfall within the catchment. This study is expected to provide
sions (Barua et al., 2012). The catchment significantly contributes to
an important contribution towards the enhanced estimation of rainfall
the water supply in Melbourne and has been playing an important role
in the aforementioned two Australian catchments using the cokriging
in the way Melbourne has been developed and grown (Adhikary et al.,
methods by incorporating elevation as an auxiliary variable in addition
2015). The Yarra River catchment consists of three distinctive sub‐
to rainfall data. One specific contribution of this paper is in explaining
catchments including Upper, Middle, and Lower Yarra segments based
how rainfall varies with elevation from catchment to catchment.
on different land use activities (Barua et al., 2012). This study concen-
The Yarra River catchment and the Ovens River catchment in
trates on the Middle segment of the Yarra River catchment, which
Victoria are selected for this study because they are two important
forms part of the case study area in Figure 1. The Upper Yarra segment
water resources catchments in Australia in terms of water supply and
includes mainly the forested and mountainous areas with minimum
agricultural production (Adhikary, Yilmaz, & Muttil, 2015; EPA Victoria,
human settlement, which is mainly used as a closed water supply
4
ADHIKARY ET
FIGURE 1
AL.
Location and topography of (a) Middle Yarra River catchment and (b) Ovens River catchment in Victoria with existing rain‐gauge stations
catchment for Melbourne and has been reserved for more than
The Ovens River catchment in northeast Victoria is also consid-
100 years for water supply purposes. The MiddleYarra segment, starting
ered as a part of the case study area in this study, which is shown in
from the Warburton Gorge to Warrandyte Gorge, is notable as the only
Figure 1. The catchment covers an area of 7,813 km2 (Yu et al.,
part of the catchment with an extensive flood plain, which is mainly used
2013), which extends from the Great Dividing Range in the south
for agricultural activities. The Lower Yarra segment of the catchment,
to the Murray River in the north, with the Yarrawonga Weir forming
which is located downstream of Warrandyte, is mainly characterized
the downstream boundary. It is considered to be one of the least
by the urbanized floodplain areas of Melbourne city (Adhikary, Muttil,
modified catchments within the Murray‐Darling basin. The catchment
& Yilmaz, 2016b). Most of the land along rivers and creeks in the middle
contributes approximately 14% to the average flows of the Murray
and lower segments has been cleared for the agricultural or urban
River in spite of its relatively small size (0.75 percent of the total
development (Barua et al., 2012; Melbourne Water, 2015).
Rainfall varies significantly through different segments of the Yarra
Murray‐Darling Basin area; EPA Victoria, 2003; Yu et al., 2013). The
Ovens River is the main river in the catchment, which originates on
River catchment. The mean annual rainfall varies across the catchment
the northern edges of the Victorian Alps and flows in a north‐west-
from 600 mm in the Lower Yarra segment to 1,100 mm in the Upper
erly direction until its junction with the Murray River near Lake
Yarra segment (Daly et al., 2013). The Middle Yarra segment (part of
Mulwala. The riverine plains and alluvial flats are primarily cleared
the case study area in Figure 1) covers an area of 1,511 km2 and consists
for agricultural use, while the hills and mountains are covered by
of three storage reservoirs. Decreasing rainfall patterns in the catchment
forests with native plant species (Yu et al., 2013). Total average water
will reduce the streamflows, which in turn will lead to the reduction in
use in the catchment is about 30,000 million litres per year, 64% of
reservoir inflows and hence impact the overall water availability in the
which is diverted from the Ovens River and its tributaries. A major
catchment. Moreover, the reduced streamflows may cause increased risk
part of this water use is irrigation, which constitutes more than
of bushfires. Conversely, increasing rainfall patterns and the occurrence
16,000 million litres annually (Schreider et al., 1996). The river itself
of extreme rainfall events (as reported in Yilmaz and Perera [2014] and
provides natural conditions suitable for many significant native fish
Yilmaz, Hossain, and Perera [2014]) will result in excess amount of
species, particularly the endangered Murray Cod (EPA Victoria,
streamflows that may cause flash floods in the urbanized lower segment
2003). Thus, the catchment is considered to be important, not only
and makes it vulnerable and risk‐prone. The urbanized lower part of the
at a regional scale, but also at the national scale in terms of its water
catchment is also dependent on the water supply from the storage
supply volume for domestic and agricultural production, and high
reservoirs mainly located in the middle and upper segments of the
environmental value.
catchment (Adhikary et al., 2015). Therefore, accurate spatial distribution
The climate of the Ovens River catchment varies considerably
of rainfall in the middle and upper segments of the catchment could be
with topography and elevation (Yu et al., 2013). The average annual
useful for accurate estimation of future streamflows for optimal reservoir
rainfall varies from 1,127 mm in the Alpine region at Bright to
operation and effective flood control in the urbanized lower part.
636 mm on the alluvial plains in Wangaratta with most rainfall
ADHIKARY ET
5
AL.
occurring in winter months (Yu et al., 2013). Approximately 45% of the
part in the study area. September is the wettest month (rainfall amount
annual rainfall occurs during the winter (June to September) whereas
equals to 112.5 mm) with the highest variation in rainfall. The driest
the summer is warm and dry. Winter snowfalls are common at alti-
month is February (rainfall amount equals to 56.4 mm) with the second
tudes above 1,000 m (EPA Victoria, 2003). Therefore, enhanced esti-
highest variability. On the other hand, the annual average rainfall for
mation of rainfall and its spatial distribution could be useful for the
the same period in the Ovens River catchment varies from 231 to
effective management of water supply and agricultural activities in
2,473 mm with a mean value of 913 mm. The wettest month is July
the catchment.
(rainfall amount equals to 112.9 mm) with the highest variability and
February (rainfall amount equals to 51.3 mm) appears to be the driest
2.2
|
month with the third highest variation in rainfall.
Data used
For the OCK analysis, a DEM of both the catchments with 10 m
In this study, historical rainfall data from existing rain‐gauge stations in
resolution (shown in Figure 1) is collected from the Geoscience
the Middle Yarra River catchment and the Ovens River catchment
Australia. The elevation of the Middle Yarra River catchment varies
(Figure 1) for 1980–2012 period are considered. There are 19 rain‐
from 25 m (lowest‐mainly in central, north‐western, and western part)
gauge stations in the Middle Yarra River catchment, whereas the
to 1,243 m (highest‐mainly in northern, north‐eastern, and eastern
Ovens River catchment includes 42 rain‐gauge stations operated by
part) with a mean elevation of 621 m above mean sea level. Whereas,
the Australian Bureau of Meteorology. Daily rainfall data are collected
the elevation of the Ovens River catchment varies from 124 m (lowest‐
from the Scientific Information for Land Owners climate database
mainly in the upper north‐western part) to 1903 m (highest‐mainly in
(http://www.longpaddock.qld.gov.au/silo/) and compiled to generate
the lower southern, south‐eastern and eastern part) with a mean eleva-
monthly and annual rainfall, which are then used for the analysis. Sum-
tion of 874 m above mean sea level. Monthly and annual rainfall gen-
mary statistics of monthly rainfall data are given in Table 1. The annual
erally tend to increase with the higher elevations caused by the
average rainfall in the Middle Yarra River catchment for the aforemen-
orographic effect of mountainous terrain (Goovaerts, 2000). Several
tioned period varies from 710 mm to 1422 mm with a mean value of
studies have revealed that rainfall usually shows good correlation with
1082 mm. The southern and south‐eastern part experiences the
elevation. For example, Goovaerts (2000) showed that a good to signif-
highest rainfall, whereas the lowest rainfall occurs in the north‐western
icant correlation exists between the monthly rainfall and elevation,
TABLE 1
Summary statistics for monthly rainfall data of Middle Yarra River catchment and Ovens River catchment
Month
Mean (mm)
Minimum (mm)
Maximum (mm)
Standard deviation (mm)
Correlation coefficient
(versus elevation)a
Middle Yarra River catchment
January
67.3
1.8
201.4
31.23
0.79
February
56.4
0.0
269.5
53.40
0.69
March
66.3
9.2
217.8
36.89
0.77
April
84.7
15.2
246.0
45.81
0.74
May
88.3
10.2
239.7
42.47
0.67
June
106.4
13.8
300.2
46.18
0.70
July
102.1
17.9
303.9
48.95
0.73
August
108.6
21.1
289.7
50.67
0.72
September
112.5
25.3
350.3
56.98
0.67
October
101.6
4.0
333.9
50.86
0.69
November
96.1
15.3
258.7
47.60
0.77
December
91.3
8.2
301.2
52.88
0.74
Ovens River catchment
a
January
55.8
0.0
364.7
51.41
0.49
February
51.3
0.0
421.9
65.66
0.26
March
53.7
0.3
418.1
51.01
0.49
April
54.5
1.0
234.6
38.63
0.72
May
74.6
2.0
360.7
56.71
0.61
June
95.4
1.4
457.0
62.99
0.65
July
112.9
6.0
622.4
71.92
0.60
August
105.2
5.8
468.0
68.94
0.61
September
86.3
5.2
484.1
55.44
0.65
October
73.4
0.1
410.8
61.61
0.67
November
73.6
0.6
290.8
47.70
0.68
December
64.1
0.2
374.6
54.79
0.68
Linear correlation coefficient between rainfall and elevation data.
6
ADHIKARY ET
AL.
which varies from .33 to .83. Subyani and Al‐Dakheel (2009) found
that good correlation ranging from .34 to .77 exists between seasonal
rainfall and elevation in the Southwest Saudi Arabia. Moral (2010)
identified a good correlation ranging from .33 to .67 between monthly
and annual rainfall and elevation in the southwest region of Spain.
As can be seen in Table 1, the correlation coefficient (CC) between
the monthly average rainfall and elevation for the Middle Yarra River
catchment varies from .67 to .79, where 8 months have the CC values
greater than .70. This indicates that a strong correlation exists between
the monthly rainfall and elevation in the catchment, suggesting that
elevation may enhance the monthly rainfall estimates when used as a
secondary variable in the OCK analysis. On the contrary, the CC
between the monthly average rainfall and elevation for the Ovens
River catchment varies from .26 to .72, where 6 months exhibit the
CC values higher than .65. Apart from the driest month of February,
the correlation ranges from .49 to .72 in the Ovens River catchment.
Thus, it seems beneficial taking into account this exhaustive secondary
variable (elevation in this study) into the enhanced estimation and
mapping of rainfall in both the catchments. Goovaerts (1997) men-
2 Methodological framework adopted in this study.
IDW = inverse distance weighting; KED = kriging with an external
drift; OCK = ordinary cokriging; OK = ordinary kriging; RBF = radial
basis function
FIGURE
tioned that use of multiple elevation data other than the colocated
positions of rain‐gauges can lead to unstable cokriging systems
because the correlation between near elevation data is much greater
than the correlation between distant rainfall data. Therefore, the
colocated elevation data are used for the OCK analysis in this study,
which are extracted at the same positions of rain‐gauge stations from
n
bOK ðs0 Þ ¼ ∑ ωOK Zðsi Þ
Z
i
the DEM of the catchments.
n
with
i¼1
3
METHODOLOGY
|
The methodological framework adopted in this study for spatial inter-
¼1
∑ ωOK
i
(1)
i¼1
bOK ðs0 Þ is the estimated value of variable Z (i.e., rainfall) at target
where Z
indi(at which estimation is to be made) unsampled location s0; ωOK
i
cates the OK weights linked with the sampled location si with respect
polation of rainfall includes three kriging‐based geostatistical (OK,
to s0; and n is the number of sampling points used in estimation. While
OCK, and KED) and two deterministic (IDW and RBF) interpolation
giving the estimation at target location, OK provides a variance mea-
methods, which is shown in Figure 2. A brief description of these
sure to signify the reliability of the estimation.
methods is presented in this section. The variogram and its estimation
OK is known as the best linear unbiased estimator (Isaaks &
technique are also summarised with each of the kriging methods
Srivastava, 1989). It is linear in the sense that it gives the estimation
because it is a key component of kriging. For a more detailed descrip-
based on the weighted linear combinations of observed values. It is
tion of the methods used in the current study, readers are referred to
best in the sense that the estimate variance is minimized while interpo-
several recent geostatistical textbooks including Journel and
lating the unknown value at desired location. And it is unbiased
Huijbregts (1978); Isaaks and Srivastava (1989), Goovaerts (1997),
because it tries to have the expected value of the residual to be zero
Chilès and Delfiner (1999), Wackernagel (2003), and Webster and
(Adhikary et al., 2016b). The weight constraint in Equation 1 ensures
Oliver (2007).
the unbiased estimation in OK. For OK, the kriging weights are determined to minimize the estimation variance and ensure the unbiasedness of the estimator.
3.1
|
Ordinary kriging
Kriging refers to a family of generalized least‐squares regression
methods in geostatistics that estimate values at unsampled locations
using the sampled observations in a specified search neighborhood
(Goovaerts, 1997; Isaaks & Srivastava, 1989). OK is a geostatistical
interpolation method based on spatially dependent variance, which
The OK weights ωOK
i can be obtained by solving a system of (n + 1)
simultaneous linear equations as follows:
n
∑ γ si −sj ωOK
þ μOK
i
1 ¼ γ sj −s0
i¼1
n
for j ¼ 1; … … … …; n
(2)
∑ ωOK
¼1
i
i¼1
gives unbiased estimates of variable values at target location in space
using the known sampling values at surrounding locations. The unbi-
where γ(si − sj) is the variogram values between sampling locations si
asedness in the OK estimates is ensured by forcing the kriging weights
and sj, γ(sj − s0) is the variogram values between sampling location sj
to sum to 1. Thus, the OK estimator may be stated as a linear combina-
and the target location, s0, and μOK
1 is the Lagrange multiplier parame-
tion of variable values, which is given by
ter. Equation 2 indicates that OK highly depends on a mathematical
ADHIKARY ET
7
AL.
function called variogram model that indicates the degree of spatial
functional forms these variogram models are given in Table 2. The
autocorrelation in datasets.
three models are fitted to the experimental variogram using regression
For OK interpolation of variables, first an experimental variogram
b
γðdÞ is derived by
by noting the residual sum of squares (RSS) between the experimental
b
γðdk Þ and modelled γ(dk) variogram values (Mair & Fares, 2011) with a
trial‐and‐error approach for different lag sizes and lag intervals
1 NðdÞ
b
γðdÞ ¼
∑ ½Zðsi þ dÞ−Zðsi Þ2
2NðdÞ i¼1
(3)
(Goovaerts, 1997) such that the RSS is minimum. RSS is given by
K
2
RSS ¼ ∑ ½b
γðdk Þ−γðdk Þ
pling locations si and (si + d), respectively, located at d distance apart
and N(d) is the number of data pairs. A variogram cloud is initially generated using Equation 3 for observations at any two data points, in
which all semivariance values are plotted against their separation distance. The experimental variogram is computed from the variogram
cloud by subdividing it into a number of lags and taking an average
of each lag interval (Johnston, VerHoef, Krivoruchko, & Lucas, 2001;
Robertson, 2008). A variogram model γ(d) is then fitted to the experimental variogram. A typical variogram cloud based on Equation 3 and
a typical experimental variogram with a typical fitted model is shown in
Figure 3.
Exponential, Gaussian, and spherical are the most commonly used
variogram models for kriging applications in hydrology (Adhikary et al.,
2015), which are also used to model the experimental variogram. The
(4)
k¼1
where Z(si) and Z(si + d) are the variable values at corresponding sam-
RSS in Equation 4 provides an exact measure of how well the
variogram model fits the experimental variogram (Robertson, 2008).
Lag sizes and number of lags are varied based on a general rule of
thumb, in which the lag size times the number of lags should be about
half of the largest distance among all data pairs in the variogram cloud
(Johnston et al., 2001, p. 66). The variogram parameters (nugget, sill,
and range) are also iteratively changed to obtain the best fitted model.
The model with its corresponding parameters that minimizes RSS is
selected as the best variogram model and finally used in OK analysis.
Variogram model fitting is performed using GS+ geostatistical software
(Robertson, 2008) and OK is implemented through ArcGISv9.3.1 software (ESRI, 2009) and its geostatistical analyst extension (Johnston
et al., 2001).
3.2
|
Ordinary cokriging
OCK method is a modification of the OK method. The key advantage of
OCK is that it can make use of more than one variable rather than using
only a single variable in the estimation process. The OCK method is
used to enhance the estimation of primary variable by using secondary
variable assuming that the variables are correlated to each other (Isaaks
& Srivastava, 1989). In this study, rainfall and elevation are considered,
respectively, as the primary and secondary variables in the OCK
method. Like OK method, the aim of the OCK method is to estimate
the primary variable. The OCK estimator (Goovaerts, 1997) considering
one secondary variable (i.e., elevation), which is cross‐correlated with
the primary variable (i.e., rainfall) may be written as
m
n
bOCK ðs0 Þ ¼ ∑ ωOCK Zðsi Þ þ ∑ ωOCK V si
Z
i2
i1
1
2
with
(5)
i2 ¼1
i1 ¼1
n
m
i1 ¼1
i2 ¼1
∑ ωOCK
¼0
¼ 1; ∑ ωOCK
i2
i1
bOCK ðs0 Þ is the estimated value of primary variable at target
where Z
unsampled location s0, ωOCK
and ωOCK
are the kriging weights associated
i1
i2
with the sampling locations of the primary and secondary variables Z
TABLE 2
Commonly used positive‐definite variogram models
Model name
Exponential
Gaussian
Spherical
(a) a typical variogram cloud for a finite set of discrete lags,
and (b) a typical experimental variogram based on the variogram cloud
fitted by a typical variogram model with its parameters
FIGURE 3
Model equation
γðdÞ ¼ C0 þ C1 1− exp − 3d
a
h
i
2
γðdÞ ¼ C0 þ C1 1− exp − 3d
a2
h 3 i
γðdÞ ¼ C0 þ C1 32 da − 12 da3 ;
=C0 + C1 ,
d<a
d≥a
C0 = nugget coefficient, C0 + C1 = Sill, a = range of variogram model.
d = distance of separation between two locations.
8
ADHIKARY ET
AL.
and V, respectively, n and m are the number of sampling points for the
(elevation) variables in the OCK method are obtained by fitting with
primary and secondary variables.
an experimental cross‐variogram that is given by
The OCK weights are obtained by solving a system of (n + 2) simultaneous linear equations (Goovaerts, 1997) that can be given by
m
n
þ μOCK
¼ γzz sj1 −s0
þ ∑ γzv si2 −sj1 ωOCK
∑ γzz si1 −sj1 ωOCK
i2
1
i1
i2 ¼1
i1 ¼1
n
m
∑ γvz si1 −sj2 ωOCK
þ ∑ γvv si2 −sj2 ωOCK
þ μOCK
i1
i2
2
i1 ¼1
n
∑ ωOCK
¼1
i1
i2 ¼1
¼ γvz sj2 −s0
for j1 ¼ 1; ……; n
for j2 ¼ 1; ……; m
1 NðdÞ
∑ ½Zðsi þ dÞ−Zðsi Þ×½V ðsi þ dÞ−V ðsi Þ
2NðdÞ i¼1
(7)
It is important to note that the variogram models must satisfy the
positive‐definite condition (PDC) in kriging. For a single variable (rainfall) in the OK method, the condition is satisfied by using the posi-
i1 ¼1
m
tive‐definite variogram model functions given in Table 2. However,
¼0
∑ ωOCK
i2
the OCK method considering two variables (rainfall and elevation) con-
i2 ¼1
(6)
where γzv si2 − sj1 and γvz si1 − sj2 are the cross‐variogram values
between sampled Z and V values, and
b
γzv ¼ b
γvz ¼
μOCK
1
and
μOCK
2
are the Lagrange
multiplier parameters accounting for the two unbiased conditions.
The elementary step in the OCK method is to establish an appro-
sists of one cross‐variogram and two direct variograms, and additional
requirement for satisfying the PDC (Goovaerts, 1999). In order to
make sure that the cross‐variogram model is positive‐definite (all
eigenvalues are positive), an indicator called the Cauchy‐Schwarz
inequality (Journel & Huijbregts, 1978; Phillips et al., 1992) must be
satisfied for all distance values d, which is given by
priate model for cross continuity and dependency between the primary (rainfall) and secondary (elevation) variable. This positive
1
γzv ðdÞ≤½γzz ðdÞγvv ðdÞ2
(8)
correlation between variables is referred to as the cross‐regionalization
or coregionalization (Goovaerts, 1997; Wackernagel, 2003), which can
where γzv(d) is the cross‐variogram model between primary and sec-
be quantified by cross‐variogram or cross‐covariance. These models
ondary variables, and γzz(d) and γvv(d) are the direct variogram models
are used to define the cross continuity and dependency between
for primary and secondary variables, respectively. Based on the indica-
two variables in the OCK method (Subyani & Al‐Dakheel, 2009). The
tor shown in Equation 8, Hevesi et al. (1992) suggested a graphical test
cross‐variogram models between the primary (rainfall) and secondary
of PDC for the fitted models as follows:
Experimental variograms and fitted variogram models for the monthly rainfall data of the Middle Yarra River catchment used in the
ordinary kriging interpolation method
FIGURE 4
ADHIKARY ET
9
AL.
1
n
(9)
PDC ¼ ½γzz ðdÞγvv ðdÞ2
bKED ðs0 Þ ¼ ∑ ωKED Zðsi Þ
Z
i
i¼1
The model is considered positive‐definite if the absolute value of
the cross‐variogram model γzv(d) in Equation 8 is less than the slope
and corresponding absolute value of PDC curve in Equation 9 for all
distance values d. In the OCK method, direct and cross‐variogram
models are fitted as linear combination of the same set of basic models
given in Table 2 such that the RSS value by Equation 4 is minimum
under the requirement of PDC. Variogram model fitting is performed
using GS+ geostatistical software (Robertson, 2008) and OCK is implemented through ArcGISv9.3.1 software (ESRI, 2009) and its
n
¼1
with ∑ ωKED
i
The KED weights ωKED
can be obtained by solving a system of
i
(n + 2) simultaneous linear equations as follows:
n
∑ γR si −sj ωKED
þ μKED
þ μKED
i
0
1 V sj ¼ γR sj −s0
|
for j ¼ 1; … …; n
i¼1
n
¼1
∑ ωKED
i
;
i¼1
n
V ðsi Þ ¼ V ðs0 Þ
∑ ωKED
i
i¼1
(11)
geostatistical analyst extension (Johnston et al., 2001).
3.3
(10)
i¼1
where γR(si − sj) is the residual variogram values between sampling
Kriging with an external drift
locations si and sj, γR(sj − s0) is the residual variogram values between
KED is a particular type of universal kriging that gives the estimation of
sampling location sj and the target location, s0, and μKED
are
and μKED
1
0
a primary variable Z, known only at a small number of locations in the
the Lagrange multiplier parameters.
study catchment, through a secondary variable V, exhaustively known
OCK and KED differ in the way the secondary variable V is used.
in the same area (Feki et al., 2012). The trend or local mean of the pri-
The secondary variable (e.g., elevation in this study) gives only the
mary variable is first derived using the secondary variable (Goovaerts,
trend information in KED, whereas estimation with OCK is directly
1997; Wackernagel, 2003) and then simple kriging is carried out on
influenced by it (Delbari et al., 2013). In case of KED, the primary
residuals from the local mean. The KED estimator (Wackernagel,
and secondary variables should exhibit a linear relationship. In addition,
2003) is generally given as
estimation with KED requires the secondary variable values at all the
TABLE 3
Results of fitted variogram models for monthly rainfall data for using in the OK interpolation method
Variogram parameters
Month
Model name
2
Nugget, C0 (mm )
2
Cross‐validation statistics
Sill, C0 + C1 (mm )
Range, a (km)
SM
SRMS
Middle Yarra River catchment
January
Spherical
0.15
163.65
25.20
0.068
0.998
February
Spherical
0.20
58.70
24.75
0.059
1.000
March
Spherical
0.00
205.38
25.25
0.059
0.863
April
Spherical
1.25
270.75
24.85
0.054
0.897
May
Spherical
1.10
327.20
25.25
0.039
0.916
June
Spherical
0.10
727.60
26.85
0.041
0.847
July
Spherical
0.10
849.60
27.00
0.031
0.795
August
Spherical
4.00
839.00
26.11
0.035
0.851
September
Spherical
1.00
816.00
27.07
0.045
0.883
October
Spherical
0.00
390.00
24.18
0.049
0.903
November
Spherical
0.00
292.50
21.53
0.058
0.861
December
Spherical
0.00
320.00
22.00
0.060
0.863
Ovens River catchment
January
Spherical
92.00
455.00
109.50
‐0.048
0.992
February
Spherical
83.00
360.00
70.00
0.007
1.024
March
Spherical
85.00
501.00
107.10
‐0.021
0.991
April
Gaussian
93.00
793.00
162.29
‐0.060
1.000
May
Spherical
109.00
1016.00
108.60
‐0.047
1.016
June
Spherical
50.00
1240.00
90.00
‐0.039
1.001
July
Gaussian
180.00
2030.00
80.50
‐0.083
0.992
August
Spherical
58.00
1420.00
65.00
‐0.040
0.990
September
Spherical
25.00
937.00
67.00
‐0.028
0.997
October
Spherical
102.00
519.00
75.40
‐0.027
1.001
November
Spherical
19.00
430.00
68.00
0.001
0.988
December
Spherical
35.00
296.00
89.10
‐0.028
0.988
Note. OK = ordinary kriging; SM = standardized mean error; SRMS = standardized root mean square error.
10
ADHIKARY ET
AL.
estimation grid nodes as well as all the sampling locations si. The resid-
analysis of spatial correlation unlike kriging. When interpolating a grid
ual variogram models are fitted based on the basic models in Table 2
node, the basis kernel functions define the optimal set of weights to be
such that the RSS value between the experimental and modelled
used with the sampling points. There are several radial basis functions
variogram values by Equation 4 is minimum. Variogram model fitting
available (Johnston et al., 2001) However, thin plate spline is the most
and estimation with KED are performed using GS+ geostatistical soft-
commonly used radial basis function for interpolation (e.g., Boer, de
ware (Robertson, 2008).
Beurs, & Hartkamp, 2001; Di Piazza et al., 2011; Hutchinson, 1995).
In this study, thin plate spline is also used as the radial basis function,
3.4
|
which is given by
Inverse distance weighting
IDW interpolation method (ASCE, 1996) gives a linear weighted aver-
∅ðr Þ ¼ ðcrÞ2 lnðcr Þ
age of several neighbouring observations to estimate the variable value
(14)
at target location. This method assumes that each observation point
has local influence that diminishes with distance. IDW assigns greater
where c is the smoothing parameter, which is obtained through cross‐
weights to observation points near to the target location, and the
validation process. The optimal value of the smoothing parameter is
weights diminish as a function of distance (Johnston et al., 2001).
selected based on the lowest RMSE value between the observed and
The estimation by IDW can be written as
estimated values. RBF interpolation method is performed by
ArcGISv9.3.1 software (ESRI, 2009) and its geostatistical analyst
n
bðs0 Þ ¼ ∑ ωi Zðsi Þ
Z
where ωi ¼
i¼1
d−k
i0 =
n
∑
i¼1
d−k
i0
(12)
bðs0 Þ is the estimated rainfall value at desired location s0, Z(si) is
where Z
extension (Johnston et al., 2001).
3.6
|
Assessment of interpolation methods
the Z value at location si, ωi is the weight assigned to observation
The performance of different interpolation methods (OK, OCK, KED,
points, di0 is the distance between the sampling point at locations si
IDW, and RBF) used in this study are evaluated and compared through
and s0, n is the number of sampling points, and k is a power, which is
cross‐validation process. The cross‐validation is a simple leave‐one‐out
referred to as a control parameter.
As “k” approaches zero and the weights becomes more similar,
IDW estimates approach the simple average of the surrounding observations. However, the effect of the farthest observations on the
estimated value is diminished with the increase of k. The value of k
ranges from 1 to 6 (Teegavarapu & Chandramouli, 2005). Several studies have investigated with variations in a power to examine its effects
on the spatial distribution of information from rainfall observations
(Chen & Liu, 2012). Therefore, k value is varied in the range of 1 to 6
with an increment of 0.1 in the current study. The optimal k value is
selected based on the lowest RMSE value between the observed and
estimated values. All rain‐gauge stations are considered in the search
neighbourhood in the estimation process. IDW interpolation is performed by ArcGISv9.3.1 software (ESRI, 2009) and its geostatistical
analyst extension (Johnston et al., 2001).
3.5
|
Radial basis function
RBF (Chilès & Delfiner, 1999) is an exact interpolation method, which
stands for a diverse group of interpolation methods. The RBF estimator can be viewed as a weighted linear function of distance from grid
point to data point plus a bias factor, which is given by
n
bðs0 Þ ¼ ∑ ωi ∅ðksi − s0 kÞ þ μ
Z
(13)
i¼1
where ∅(r) is the radial basis function (r = ‖si − s0‖), r is the radial distance from target point s0 to sampling points si, ωi are the weights
and μ is the Lagrangian multiplier. The weights are obtained by solving
of a system of (n + 1) simultaneous linear equations.
The basis kernel functions in the RBF method are analogous to
variograms in kriging, which makes it similar to geostatistical interpolation methods. However, RBF does not have the advantage of a prior
Experimental variograms and fitted variogram models
based on the collocated elevation data for (a) Middle Yarra River
catchment, and (b) Ovens River catchment
FIGURE 5
ADHIKARY ET
11
AL.
validation procedure (Haddad, Rahman, Zaman, & Shrestha, 2013) in
with the lowest MBE and RMSE values and the highest R2 value is cho-
which observations are removed one at a time from the dataset and
sen as the best interpolation method.
then re‐estimated from the remaining observations using the adopted
As has been mentioned earlier, kriging gives the prediction
model. Cross‐validation provides important evidence of the perfor-
standard error while giving the estimation of unsampled variables,
mance measures for the interpolation methods. In this study, the per-
the adequacy of the variogram model for kriging and cokriging estima-
formance of all the interpolation methods for rainfall estimation is
tion should also be tested to produce correct interpolation results
compared based on mean bias error (MBE), RMSE, and coefficient of
(Johnston et al., 2001; Phillips et al., 1992). Therefore, two additional
determination (R2) values between the observed and estimated rainfall
standardized cross‐validation statistics are used in this study, which
values, which are given by Equations 15–17.
are standardized mean error (SM) and standardized root mean square
error (SRMS) as given by Equations 18–19
i
1 n h
bðsi Þ
MBE ¼ ∑ Zðsi Þ−Z
n i¼1
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
i2
1 n h
bðsi Þ
RMSE ¼
∑ Zðsi Þ−Z
n i¼1
h
n
oi2
bðsi Þ−Z
bm
∑ni¼1 fZðsi Þ−Zm g Z
2
R ¼
n
o2
bm
bðsi Þ−Z
∑n fZðsi Þ−Zm g2 ∑n Z
i¼1
(15)
(16)
(17)
i¼1
bðsi Þ are the observed and predicted
In Equations 15–17, Z(si) and Z
bm are the mean of the observed and predicted values,
values, Zm and Z
and n is the number of sampled data points. The interpolation method
h
i
b
1 n Zðsi Þ−Zðsi Þ
SM ¼ ∑
bðsi Þ
n i¼1
σ
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
"
#2ffi
u
u1 n Z ðs Þ−Z
bðsi Þ
i
t
SRMS ¼
∑
bðsi Þ
n i¼1
σ
(18)
(19)
bðsi Þ is the prediction standard error for location si. SM should
where σ
be close to 0 if the estimates using the adopted variogram model are
unbiased. SRMS should be close to 1 if the estimation variances are
consistent and the variability of the prediction is correctly assessed
(Adhikary et al., 2015; Johnston et al., 2001).
FIGURE 6 Experimental cross‐variograms with the fitted cross‐variogram models and positive‐definite condition curve based on the monthly
rainfall and collocated elevation data for the Middle Yarra River catchment used in the ordinary cokriging interpolation method
12
ADHIKARY ET
4
RESULTS A ND DIS CUS SION
|
AL.
For convenience in this study, results obtained for the Middle
Yarra River catchment are presented and discussed elaborately and
4.1
|
compared with that obtained for the Ovens River catchment.
Variogram models for OK analysis
Figure 4 shows the experimental variograms and fitted variogram
The OK analysis requires the estimation of the direct variogram models
models with optimal variogram parameters (i.e., nugget, sill, and range)
for rainfall data. In this study, an isotropic experimental variogram is
for monthly rainfall data of the Middle Yarra River catchment, which
estimated from the rainfall dataset for each month assuming an identi-
are used in OK analysis. The best fitted variogram models are selected
cal spatial correlation in all directions and neglecting the influence of
based on the minimum RSS values using a trial‐and‐error process with
anisotropy on the variogram parameters. Isotropy is assumed for the
different lag intervals. The variogram parameters are iteratively
methodological simplicity. Isotropy is a property of a natural process
changed to get the best fitted model, which yields the minimum RSS.
or data where directional influence is considered insignificant and spa-
As can be seen from Figure 4, the spherical model is the best fitted
tial dependence (autocorrelation) changes only with the distance
variogram model for all months having spatial structure of
between two locations (Johnston et al., 2001). Under the isotropic con-
.66 < R < .97 (results of R not shown in Figure). It can be also noted that
dition, the semivariance is assumed the same for a given distance,
10 months (except November and December) have R values greater
regardless of direction. Initially, the directional experimental variograms
than .75. The optimal variogram parameters for each monthly rainfall
are estimated from each monthly rainfall dataset. However, the direc-
dataset for both the catchments are provided in Table 3. As can be
tional variograms are found noisy because of the less number of rain‐
seen from Table 3, the ratio of nugget coefficient to sill of the
gauge stations in the study area. Therefore, the directional influence
variogram model is small for all months for both the catchments. This
is ignored in the experimental variogram calculation. The experimental
evidently indicates that a strong spatial correlation exists between
variogram is then fitted with three predefined variogram model func-
the monthly mean rainfall and the spatial distribution of the rain‐gauge
tions (exponential, Gaussian, and spherical as given in Table 2) to obtain
stations over the study area. This supports the use of geostatistical
the variogram models for each monthly rainfall dataset.
interpolation methods such as OK, OCK, and KED, which consider
TABLE 4
Results of fitted cross‐variogram models between monthly rainfall and elevation for using in the OCK interpolation method
Variogram parameters
Month
Model name
2
Nugget, C0 (mm )
2
Cross‐validation statistics
Sill, C0 + C1 (mm )
Range, a (km)
SM
SRMS
1.001
Middle Yarra River catchment
January
Gaussian
1.0
1387.0
12.65
0.053
February
Gaussian
40.0
801.0
16.19
0.067
1.000
March
Gaussian
1.0
1545.0
14.36
0.055
1.006
April
Gaussian
150.0
1653.0
18.35
0.057
1.004
May
Gaussian
80.0
1651.0
19.40
0.042
1.019
June
Gaussian
300.0
1874.5
19.75
0.055
1.001
July
Gaussian
55.0
3081.2
13.50
0.008
1.002
August
Gaussian
57.5
2872.9
12.10
0.007
1.006
September
Gaussian
179.4
2746.0
14.25
0.042
1.002
October
Gaussian
114.5
1860.9
15.85
0.053
1.004
November
Gaussian
3.4
1804.4
12.00
0.034
1.000
December
Gaussian
1.7
1774.7
15.66
0.096
1.002
Ovens River catchment
January
Gaussian
485.2
2197.9
60.98
‐0.020
0.989
February
Gaussian
525.8
2540.3
50.56
0.001
1.119
March
Gaussian
635.1
2349.5
52.45
‐0.012
1.093
April
Gaussian
373.7
2209.3
36.01
‐0.019
0.904
May
Gaussian
595.0
3170.7
32.99
‐0.015
1.090
June
Gaussian
342.6
3961.7
52.48
‐0.022
1.117
July
Gaussian
652.2
4175.6
71.62
‐0.030
1.028
August
Gaussian
919.2
4223.1
96.48
‐0.024
1.095
September
Gaussian
733.6
3717.8
66.49
‐0.018
1.000
October
Gaussian
256.0
2048.1
29.73
‐0.010
1.010
November
Gaussian
413.5
1925.5
64.09
‐0.013
1.050
December
Gaussian
123.5
1935.7
32.60
‐0.021
1.007
Note. OCK = ordinary cokriging; SM = standardized mean error; SRMS = standardized root mean square error.
ADHIKARY ET
13
AL.
the spatial correlation in the estimation process. The range of influence
4.2
|
Cross‐variogram models for OCK analysis
and the sill of the variogram model vary from one month to another,
but the variogram exhibit a same spherical structure in all months. This
The OCK analysis requires the simultaneous estimation of the direct
may be caused due to the control of the relief on the spatial distribu-
and cross‐variogram models for the rainfall and elevation variables.
tion of rainfall (Delbari et al., 2013). The range of influence is lowest
The three variogram models are fitted as a linear combination of the
for November (21.53 km) and highest for September (27.07 km). Fur-
same set of standard models given in Table 2 so that the RSS value is
thermore, the cross‐validation statistics in Table 3 for both the catch-
minimum under the constraints of PDC (Goovaerts, 1999). Figure 6
ments confirm that the fitted variogram models for all monthly
shows the experimental and fitted cross‐variogram models for the
rainfall data satisfy the unbiased condition and thus can be used for
Middle Yarra River catchment. The number of data pairs in each lag size
the OK analysis.
For elevation data, an isotropic experimental variogram is com-
is the same for all the three direct and cross‐variogram models.
puted ignoring the directional influence. The experimental variogram
to examine the positive‐definiteness criteria of the cross‐variogram
Figure 6 also shows the PDC curve computed based on Equation 9
is then fitted with the aforementioned three variogram model func-
models obtained for the catchment. Additionally, the cross‐validation
tions. The best fitted variogram model is selected using the same pro-
statistics are used for identifying the adequacy and final selection of
cedure described above. The Gaussian variogram model gives the best
the adopted cross‐variogram model for the OCK analysis. The cross‐
fitted model for both the catchments with the lowest RSS value, which
validation results obtained using all the adopted cross‐variogram
is shown in Figure 5. The optimal variogram parameters for both the
models for both the catchments are presented inTable 4. The results in
catchments are also shown in the figure. In order to avoid possible
Table 4 indicate that the cross‐variogram models of all monthly
inconsistencies in the subsequent modelling of direct and cross‐
datasets are
variograms in OCK analysis (Goovaerts, 1997, 2000), the colocated
neighbourhoods for both the catchments.
suitable for
the OCK analysis considering all
elevation data (see Figure 1) are used for estimating the variogram of
As can be seen from Figure 6, the Gaussian variogram model fits
elevation, not the entire DEM of the catchment. Therefore, the
well for all monthly datasets of the Middle Yarra River catchment,
cokriging method adopted in this study is referred to as the colocated
which also satisfy the PDC criteria defined by Equation 9. Further-
OCK method (Wackernagel, 2003).
more, the correlation between monthly rainfall and elevation for all
Experimental residual variograms and fitted residual variogram models for the Middle Yarra River catchment used in the kriging with an
external drift interpolation method
FIGURE 7
14
ADHIKARY ET
months in Table 1 indicates that elevation will contribute to enhance
4.3
|
AL.
Residual variogram models for KED analysis
the monthly rainfall estimation in the catchment. The figure also shows
that the values of the sample cross‐variogram increase for distances
In order to implement the KED analysis, experimental residual
from 0 to 25 km (more than half of the maximum interstation distance)
variograms are estimated based on the residuals obtained from linear
for almost all months. This indicates that a positive spatial cross‐
regression between rainfall and elevation data neglecting the influence
correlation exists between rainfall and elevation in the catchment. This
of anisotropy on the variogram parameters. The experimental residual
wide ranges may be due to the high correlation (.67 < R < .79) between
variogram is then fitted using the three standard models given in
the monthly rainfall and elevation (Table 1). Such high correlation con-
Table 2. Figure 7 shows the experimental and fitted residual variogram
firms that the monthly rainfall in the Middle Yarra River catchment is
models for all monthly datasets of the Middle Yarra River catchment. It
mainly caused by the orographic effects.
Figure 6 also shows the PDC curves, which are computed to
can be seen from the figure that the spherical model gives the best
fitted model for all monthly datasets. The optimal variogram parame-
examine the positive‐definite conditions of the cross‐variogram
ters and the corresponding cross‐validation statistics of the selected
models of the catchment. It is worth pointing out that the PDC curve
residual variogram models for both the catchments are presented in
may give a qualitative indication for the degree of correlation. As can
Table 5. As can be also seen from Figure 7 and Table 5, the residual
be observed from Figure 6, the plotted PDC curve for most of the
variogram models exhibit relatively smaller sills than those obtained
months showed a close fit to the cross‐variogram model for smaller
from the actual rainfall datasets (see Figure 4) but they follow very sim-
distances with few exceptions (February, April, May, and June
ilar structure. This is not unexpected because the residual variograms
months). For example, the PDC curve is closer to the cross‐variogram
from the linear regression represents variation, which remains after
model in the case of January, March, July, November, and December
removing the trend (Lloyd, 2005). The cross‐validation statistics shown
months depending on the degree of correlation. This conclusion holds
in Table 5 also indicate that the residual variogram models of all
true based on the higher correlation for these months as given in
monthly datasets for both the catchments are satisfactory for the
Table 1.
KED analysis.
TABLE 5
Results of fitted residual variogram models for using in the KED interpolation method
Variogram parameters
Month
Model name
2
Nugget, C0 (mm )
2
Sill, C0 + C1 (mm )
Cross‐validation statistics
Range, a (km)
SM
SRMS
Middle Yarra River catchment
January
Spherical
0.10
43.14
8.86
0.029
1.022
February
Spherical
0.01
23.11
13.10
0.025
0.981
March
Spherical
0.10
56.69
9.35
0.060
0.999
April
Spherical
0.10
75.70
10.76
0.010
1.000
May
Spherical
0.10
115.80
11.61
‐0.034
1.031
June
Spherical
39.90
298.70
27.03
0.005
0.981
July
Spherical
5.70
257.70
11.59
‐0.024
0.990
August
Spherical
0.10
265.00
11.57
‐0.013
0.984
September
Spherical
0.10
283.50
11.77
0.027
0.982
October
Spherical
0.10
139.40
10.60
‐0.010
1.004
November
Spherical
0.10
84.94
9.69
0.054
0.992
December
Spherical
0.10
100.70
10.61
0.058
0.990
Ovens River catchment
January
Spherical
38.90
239.40
42.90
‐0.017
1.012
February
Spherical
60.10
320.90
32.72
0.004
1.029
March
Spherical
13.40
246.20
26.80
‐0.011
0.996
April
Gaussian
67.20
188.00
103.05
‐0.043
1.008
May
Spherical
217.00
597.00
108.40
‐0.026
1.071
June
Spherical
60.00
841.00
70.20
‐0.044
1.028
July
Spherical
5.00
1713.00
106.00
‐0.092
1.005
August
Gaussian
250.00
2084.00
90.54
‐0.045
1.001
September
Spherical
1.00
975.00
100.85
‐0.027
1.029
October
Spherical
42.10
383.60
87.00
‐0.018
1.041
November
Spherical
6.00
367.30
92.30
‐0.012
0.989
December
Spherical
31.00
192.50
88.75
‐0.028
0.993
Note. KED = kriging with an external drift; SM = standardized mean error; SRMS = standardized root mean square error.
ADHIKARY ET
4.4
|
15
AL.
Spatial prediction of rainfall
the catchment. As can be seen from the table, geostatistical (OK,
OCK, and KED) interpolation methods perform better than determinis-
In this study, different geostatistical and deterministic interpolation
tic (IDW and RBF) interpolation methods for monthly rainfall estima-
methods including OK, OCK, KED, IDW, and RBF are adopted to
tion in the study area. The OCK method gives the best results for
estimate the spatial distribution of monthly mean rainfall in the
rainfall estimation over the study area for all months when considering
Middle Yarra River catchment and the Ovens River catchment in
all the performance measures. The KED method gives the second best
Australia. Several performance measures including MBE, RMSE, and
results, which is very close to the performance of the OCK method but
R2 are frequently used to indicate how accurately an interpolator
performs better than the OK method for both the catchments. IDW
predicts the observed data. Smaller values of MBE and RMSE with
and RBF give similar performance with higher error in rainfall estima-
a higher R2 value of an interpolator indicate better prediction by
tion over the study area. For the Middle Yarra River catchment,
the corresponding method. In case of the scatter plot, the better pre-
Table 6 also shows that in some months, the RBF method performs
diction is that if all scattered points lay close to the 450 line with the
better than the OK method for rainfall estimation. However, no
highest R2 value between the predicted and observed values
remarkable differences are seen between them when considering all
(Adhikary et al., 2016a).
the performance measures.
Table 6 presents the different performance measures of the
For OK, OCK, KED, IDW, and RBF methods, the average RMSE
adopted interpolation methods (or interpolators) for estimating
values (Table 6) for the Middle Yarra River catchment are 11.93,
monthly rainfall over both the study catchments. The different interpo-
10.29, 10.85, 12.66, and 12.22 mm, respectively, whereas the average
lation methods are quantitatively compared based on these perfor-
RMSE values for the Ovens River catchment are 17.42, 16.15, 16.65,
mance measures in order to identify the best interpolator for each of
17.54, and 23.41 mm, respectively. For OK, OCK, KED, IDW, and
TABLE 6
Performance of different interpolation (OK, OCK, KED, IDW, and RBF) methods for monthly rainfall estimation in the study area
MBE (mm)
Month
OK
OCK
R2
RMSE (mm)
KED
IDW
RBF
OK
OCK
KED
IDW
RBF
OK
OCK
KED
IDW
RBF
Middle Yarra River catchment
January
0.98
0.63
0.73
1.76
1.56
8.92
6.96
7.99
8.96
9.99
0.40
0.66
0.56
0.42
0.31
February
0.51
0.38
0.41
0.73
0.85
5.25
3.80
4.45
5.26
5.46
0.42
0.71
0.59
0.42
0.42
March
0.97
0.65
0.70
1.67
1.78
8.71
6.79
8.41
9.58
9.50
0.54
0.73
0.58
0.45
0.46
April
1.03
0.83
0.89
1.70
2.27
9.80
8.26
8.75
10.68
10.04
0.54
0.67
0.64
0.46
0.55
May
0.78
0.72
‐0.40
1.43
2.41
11.55
10.78
11.98
11.58
11.80
0.49
0.56
0.51
0.49
0.55
June
1.19
1.26
1.17
2.25
3.29
15.15
13.77
13.44
17.26
15.15
0.61
0.67
0.69
0.49
0.63
July
0.89
0.11
‐0.36
1.73
2.77
15.50
14.63
14.84
16.85
15.01
0.65
0.69
0.68
0.59
0.68
August
1.05
0.03
‐0.20
2.08
3.04
16.75
15.94
14.67
18.19
16.12
0.59
0.63
0.69
0.52
0.64
September
1.38
0.87
0.92
2.69
3.52
16.79
15.73
15.38
18.46
17.17
0.57
0.62
0.63
0.48
0.57
October
1.07
0.79
‐0.11
1.86
2.32
12.34
10.70
11.56
12.05
12.48
0.53
0.64
0.60
0.56
0.56
November
1.23
0.55
0.58
2.16
2.28
11.00
8.25
9.10
11.43
11.72
0.50
0.73
0.66
0.47
0.46
December
1.31
1.07
1.08
2.25
2.28
11.43
7.93
9.58
11.65
12.17
0.51
0.78
0.66
0.50
0.48
Average
1.03
0.66
0.45
1.86
2.36
11.93
10.29
10.85
12.66
12.22
0.53
0.67
0.62
0.49
0.54
‐0.62
‐0.56
‐0.25
‐1.78
‐1.13
12.86
12.39
12.45
13.68
18.42
0.41
0.44
0.42
0.34
0.18
0.11
‐0.01
0.08
‐0.96
1.26
17.51
16.32
17.34
17.53
25.64
0.15
0.23
0.19
0.11
0.02
‐0.28
‐0.22
‐0.15
‐1.26
‐0.25
13.13
12.30
12.87
14.44
17.14
0.44
0.59
0.51
0.33
0.32
Ovens River catchment
January
February
March
April
‐0.62
‐0.38
‐0.41
‐0.53
‐0.58
10.31
9.19
9.58
11.00
14.91
0.56
0.71
0.66
0.49
0.31
May
‐1.07
‐0.59
‐0.60
‐1.94
‐1.56
23.18
21.16
22.89
23.27
34.19
0.18
0.29
0.25
0.15
0.01
June
‐0.92
‐0.87
‐0.94
‐1.19
‐0.81
22.13
21.98
22.04
22.21
32.50
0.52
0.66
0.61
0.57
0.30
0.44
July
‐1.94
‐1.36
‐2.01
‐3.53
‐2.09
23.16
22.89
22.92
23.01
31.35
0.62
0.68
0.65
0.64
August
‐0.94
‐0.92
‐1.05
‐1.35
‐1.16
22.24
21.92
21.97
22.27
30.01
0.66
0.70
0.68
0.68
0.49
September
‐0.62
‐0.55
‐0.62
‐0.46
‐0.87
20.39
20.01
20.08
20.44
28.12
0.50
0.57
0.54
0.56
0.32
October
‐0.53
‐0.31
‐0.35
‐0.63
‐0.41
19.37
17.65
19.01
19.64
26.24
0.23
0.39
0.32
0.19
0.04
November
0.01
‐0.33
‐0.05
0.18
0.30
10.52
9.37
9.48
12.17
11.61
0.73
0.80
0.78
0.65
0.70
December
‐0.27
‐0.21
‐0.11
‐0.41
‐0.22
9.67
8.59
9.12
10.82
10.83
0.63
0.71
0.69
0.54
0.61
Average
‐0.64
‐0.53
‐0.54
‐1.15
‐0.63
17.42
16.15
16.65
17.54
23.41
0.47
0.57
0.53
0.44
0.31
Note. IDW = inverse distance weighting; KED = kriging with an external drift; MBE = mean bias error; OCK = ordinary cokriging; OK = ordinary kriging;
R2 = coefficient of determination; RBF = radial basis function; RMSE = root mean square error.
16
ADHIKARY ET
AL.
RBF methods, the average R2 values (Table 6) for the Middle Yarra
coefficient to sill) of the direct variogram models for rainfall (Table 3)
River catchment are .53, .67, .62, .49 and .54, respectively whereas
and elevation (Figure 5), and the cross‐variogram models for rainfall‐
the average R2 values for the Ovens River catchment are .47, .57,
elevation (Table 4) for both the catchments are found very small in
.53, .44, and .31, respectively. The higher R2 value in the OCK and
all months. This results in the improvement in the rainfall estimation
KED methods indicate that using elevation as a secondary variable
by the OCK method, which is thus selected as the best interpolator
brings more information in the rainfall estimation process under the
for the study area in this study. Therefore, the OCK method (the best
kriging‐based geostatistical analysis framework.
interpolator) is used to generate a continuous rainfall dataset of the
As explained by Delbari et al. (2013), using elevation as a second-
monthly average rainfall for each of the catchments, which are shown
ary variable may not always improve the prediction accuracy through
in Figure 8 and Figure 9. The created datasets are expected to be very
the OCK analysis if the spatial continuity of elevation is weaker than
useful in various hydrological and water resources planning studies
that of rainfall despite a high correlation exists between rainfall and
within the Yarra River catchment and the Ovens River catchment
elevation. In this study, the relative nugget effect (i.e., ratio of nugget
in Australia.
Spatial distribution of monthly rainfall in the Middle Yarra River catchment using the ordinary cokriging (the best interpolator in this
study) interpolation method
FIGURE 8
ADHIKARY ET
AL.
17
Spatial distribution of monthly rainfall in the Ovens River catchment using the ordinary cokriging (the best interpolator in this study)
interpolation method
FIGURE 9
5
|
C O N CL U S I O N S
show that the geostatistical methods outperform the deterministic
methods for spatial interpolation of rainfall over the study area. Specif-
In this study, three kriging‐based geostatistical (OK, OCK, and KED)
ically, the performance of the cokriging methods (OCK and KED) is
and two deterministic (IDW and RBF with thin plate spline) interpola-
better than that of other geostatistical methods. The performance of
tion methods are used for estimating spatial distribution of monthly
the RBF with thin plate spline is found practically as good as the ordi-
mean rainfall in the Middle Yarra River catchment and the Ovens River
nary kriging method for rainfall estimation, whereas the IDW method
catchment in Victoria, Australia. The objective is to compare the per-
is shown to have the worst results for the study area. OCK performs
formances of these interpolation methods to select the best interpola-
the best among all the interpolators and gives the improved estimates
tion method for generating a high quality continuous rainfall dataset in
of rainfall in all months for both the catchments. It provides the lowest
the form of a rainfall map for the study area. The elevation data
estimation errors and the highest correlations between the estimated
obtained from a DEM of the study area is used as a supplementary var-
and observed monthly average rainfall. Thus, ordinary cokriging is
iable in addition to rainfall records for the cokriging analysis using the
identified as the best interpolator in this study for estimating spatial
ordinary cokriging and kriging with an external drift methods. Results
distribution of rainfall in both the catchments. The obtained results
18
indicate that making use of elevation as an auxiliary variable in addition
to rainfall data can enhance the estimation of rainfall in a catchment
with the mountainous and/or complex terrain. This study thus recommends the use of cokriging for the generation of continuous rainfall
map especially in catchments with high spatial variation of rainfall as
ADHIKARY ET
AL.
Journal of Hydrology, 208(3–4), 187–193. https://doi.org/10.1016/
S0022‐1694(98)00155‐3
EPA Victoria. (2003). Environmental condition of rivers and streams in the
Ovens catchment. Environment Report of EPA Victoria, Publication No.
909. Available at: <http://www.epa.vic.gov.au/~/media/Publications/
909.pdf> (Accessed on: 20 November, 2016).
well as elevation.
ESRI (2009). ArcGISv9.3.1 (computer software product of ESRI). Redlands,
CA, USA: Environmental Systems Research Institute (ESRI).
ACKNOWLEDGEMEN TS
Feki, H., Slimani, M., & Cudennec, C. (2012). Incorporating elevation in rainfall
interpolation in Tunisia using geostatistical methods. Hydrological
Sciences Journal, 57(7), 1294–1314. https://doi.org/10.1080/
02626667.2012.710334
The authors acknowledge the financial support for this study through
an International Postgraduate Research Scholarship (IPRS) provided
by the Government of Australia and Victoria University, Melbourne.
The authors are thankful to the Australian Bureau of Meteorology
(BoM) and Geoscience Australia for providing the necessary data used
in this study. The authors also wish to thank the Editor and three anonymous reviewers for their valuable comments and suggestions, which
have improved the quality of this paper.
RE FE R ENC E S
Adhikary, S. K., Yilmaz, A. G., & Muttil, N. (2015). Optimal design of rain‐
gauge network in the Middle Yarra River catchment, Australia. Hydrological Processes, 29(11), 2582–2599. https://doi.org/10.1002/hyp.10389
Adhikary, S. K., Muttil, N., & Yilmaz, A. G. (2016a). Ordinary kriging and
genetic programming for spatial estimation of rainfall in the MiddleYarra
River catchment, Australia. Hydrology Research, 47(6), 1182–1197.
https://doi.org/10.2166/nh.2016.196
Adhikary, S. K., Muttil, N., & Yilmaz, A. G. (2016b). Genetic programming‐
based ordinary kriging for spatial interpolation of rainfall. Journal of
Hydrologic Engineering, 21(2) 04015062. https://doi.org/10.1061/
(ASCE)HE.1943‐5584.0001300
ASCE (1996). Hydrology handbook (2nd ed.). American Society of Civil Engineers (ASCE): Reston, VA.
Barua, S., Muttil, N., Ng, A. W. M., & Perera, B. J. C. (2012). Rainfall trend
and its implications for water resource management within the Yarra
River catchment, Australia. Hydrological Processes, 27(12), 1727–1738.
https://doi.org/10.1002/hyp.9311
Boer, E. P. J., de Beurs, K. M., & Hartkamp, A. D. (2001). Kriging and thin
plate splines for mapping climate variables. International Journal of
Applied Earth Observation and Geoinformation, 3(2), 146–154. https://
doi.org/10.1016/S0303‐2434(01)85006‐6
Chen, F. W., & Liu, C. W. (2012). Estimation of the spatial rainfall distribution using inverse distance weighting (IDW) in the middle of Taiwan.
Paddy and Water Environment, 10(3), 209–222. https://doi.org/
10.1007/s10333‐012‐0319‐1
Chilès, J. P., & Delfiner, P. (1999). Geostatistics: modeling spatial uncertainty.
New York: John Wiley & Sons.
Daly, C. (2006). Guidelines for assessing the suitability of spatial climate
data sets. International Journal of Climatology, 26(6), 707–721. https://
doi.org/10.1002/joc.1322
Daly, E., Kolotelo, P., Schang, C., Osborne, C. A., Coleman, R., Deletic, A., &
McCarthy, D. T. (2013). Escherichia coli Concentrations and loads in an
urbanised catchment: The Yarra River, Australia. Journal of Hydrology,
497, 51–61. https://doi.org/10.1016/j.jhydrol.2013.05.024
Delbari, M., Afrasiab, P., & Jahani, S. (2013). Spatial interpolation of monthly
and annual rainfall in northeast of Iran. Meteorology and Atmospheric
Physics, 122(1), 103–113. https://doi.org/10.1007/s00703‐013‐0273‐5
Di Piazza, A., Conti, F. L., Noto, L. V., Viola, F., & Loggia, G. L. (2011). Comparative analysis of different techniques for spatial interpolation of rainfall data
to create a serially complete monthly time series of precipitation for Sicily,
Italy. International Journal of Applied Earth Observation and Geoinformation,
13(3), 396–408. https://doi.org/10.1016/j.jag.2011.01.005
Dirks, K. N., Hay, J. E., Stow, C. D., & Harris, D. (1998). High resolution studies of rainfall on Norfolk Island part II: Interpolation of rainfall data.
Gittins, R. (1968). Trend surface analysis of ecological data. Journal of
Ecology, 56(3), 845–869. https://doi.org/10.2307/2258110
Goovaerts, P. (1997). Geostatistics for natural resources evaluation. New
York: Oxford University Press.
Goovaerts, P. (1999). Geostatistics in soil science: State‐of‐the‐art and perspectives. Geoderma, 89, 1–46. https://doi.org/10.1016/S0016‐7061
(98)00078‐0
Goovaerts, P. (2000). Geostatistical approaches for incorporating elevation
into the spatial interpolation of rainfall. Journal of Hydrology, 228(1–2),
113–129. https://doi.org/10.1016/S0022‐1694(00)00144‐X
Gyasi‐Agyei, Y. (2016). Assessment of radar‐based locally varying anisotropy on daily rainfall interpolation. Hydrological Sciences Journal, 61
(10), 1890–1902. https://doi.org/10.1080/02626667.2015.1083652
Haddad, K., Rahman, A., Zaman, M., & Shrestha, S. (2013). Applicability of
Monte Carlo cross validation technique for model development and
validation using generalized least squares regression. Journal of Hydrology, 482, 119–128. https://doi.org/10.1016/j.jhydrol.2012.12.041
Hancock, P. A., & Hutchinson, M. F. (2006). Spatial interpolation of large climate data sets using bivariate thin plate smoothing splines.
Environmental Modelling & Software, 21(12), 1684–1694. https://doi.
org/10.1016/j.envsoft.2005.08.005
Hevesi, J. A., Istok, J. D., & Flint, A. L. (1992). Precipitation estimation in
mountainous terrain using multivariate geostatistics (part 1 and 2). Journal of Applied Meteorology, 31(7), 661–688. https://doi.org/10.1175/
1520‐0450(1992)031<0661:PEIMTU>2.0.CO;
Hsieh, H. H., Cheng, S. J., Liou, J. Y., Chou, S. C., & Siao, B. R. (2006). Characterization of spatially distributed summer daily rainfall. Journal of
Chinese Agricultural Engineering, 52, 47–55.
Hutchinson, M. F. (1995). Interpolating mean rainfall using thin plate
smoothing splines. International Journal of Geographical Information Systems, 9(4), 385–403. https://doi.org/10.1080/02693799508902045
Isaaks, H. E., & Srivastava, R. M. (1989). Applied geostatistics. New York:
Oxford University Press.
Jeffrey, S. J., Carter, J. O., Moodie, K. B., & Beswick, A. R. (2001). Using spatial interpolation to construct a comprehensive archive of Australian
climate data. Environmental Modelling & Software, 16(4), 309–330.
https://doi.org/10.1016/S1364‐8152(01)00008‐1
Johnson, F., Hutchinson, M. F., The C, Beesley, C., & Green, J. (2016).
Topographic relationships for design rainfalls over Australia. Journal of
Hydrology, 533, 439–451. https://doi.org/10.1016/j.jhydrol.2015.12.035
Johnston, K., VerHoef, J. M., Krivoruchko, K., & Lucas, N. (2001). Using
ArcGIS geostatistical analyst. Redlands, CA, USA: ArcGIS Manual by
ESRI.
Jones, D. A., Wang, W., & Fawcett, R. (2009). High‐quality spatial climate
data‐sets for Australia. Australian Meteorological and Oceanographic
Journal, 58(4), 233–248.
Journel, A., & Huijbregts, C. (1978). Mining geostatistics. New York:
Academic Press.
Li, J., & Heap, A. D. (2011). A review of comparative studies of spatial
interpolation methods in environmental sciences: Performance and
impact factors. Ecological Informatics, 6(3–4), 228–241. https://doi.
org/10.1016/j.ecoinf.2010.12.003
ADHIKARY ET
19
AL.
Li, M., & Shao, Q. (2010). An improved statistical approach to merge satellite
rainfall estimates and rain‐gauge data. Journal of Hydrology, 385(1–2),
51–64. https://doi.org/10.1016/j.jhydrol.2010.01.023
Lloyd, C. D. (2005). Assessing the effect of integrating elevation data into the
estimation of monthly precipitation in Great Britain. Journal of Hydrology,
308(1–4), 128–150. https://doi.org/10.1016/j.jhydrol.2004.10.026
Ly, S., Charles, C., & Degré, A. (2011). Geostatistical interpolation of daily rainfall at catchment scale: The use of several variogram models in the Ourthe
and Ambleve catchments, Belgium. Hydrology and Earth System Sciences,
15(7), 2259–2274. https://doi.org/10.5194/hess‐15‐2259‐2011
Mair, A., & Fares, A. (2011). Comparison of rainfall interpolation methods in a
mountainous region of a tropical island. Journal of Hydrologic Engineering,
16(4), 371–383. https://doi.org/10.1061/(ASCE)HE.1943‐5584.0000330
Martínez‐cob, A. (1996). Multivariate geostatistical analysis of evapotranspiration and precipitation in mountainous terrain. Journal of Hydrology,
174(1–2), 19–35. https://doi.org/10.1016/0022‐1694(95)02755‐6
Melbourne Water. (2015). Port Phillip and Westernport Regional River health
strategy: Yarra catchment. Available at: <http://melbournewater.com.au/
aboutus/reportsandpublications/key‐strategies/Documents/Port%
20Phillip%20and%20Westernport%20Regional%20River%20Health%
20Strategy%20‐%20Yarra%20catchment.pdf> (Accessed on: 20
October, 2015).
Moral, F. J. (2010). Comparison of different geostatistical approaches to
map climate variables: Application to precipitation. International Journal
of Climatology, 30(4), 620–631. https://doi.org/10.1002/joc.1913
Phillips, D. L., Dolph, J., & Marks, D. (1992). A comparison of geostatistical
procedures for spatial analysis of precipitation in mountainous terrain.
Agricultural and Forest Meteorology, 58(1–2), 119–141. https://doi.
org/10.1016/0168‐1923(92)90114‐J
Robertson, G. P. (2008). GS+: geostatistics for the environmental sciences.
Gamma Design Software: Plainwell, Michigan, USA.
Schreider, S. Y., Jakeman, A. J., Pittock, A. B., & Whetton, P. H. (1996). Estimation of possible climate change impacts on water availability,
extreme flow events and soil moisture in the Goulburn and ovens
basins, Victoria. Climate Change, 34(3–4), 513–546. https://doi.org/
10.1007/BF00139304
Sharma, R. H., & Shakya, N. M. (2006). Hydrological changes and its impact
on water resources of Bagmati watershed, Nepal. Journal of Hydrology,
327(3–4), 315–322. https://doi.org/10.1016/j.jhydrol.2005.11.051
Subyani, A. M., & Al‐Dakheel, A. M. (2009). Multivariate geostatistical
methods of mean annual and seasonal rainfall in Southwest Saudi
Arabia. Arabian Journal of Geosciences, 2(1), 19–27. https://doi.org/
10.1007/s12517‐008‐0015‐z
Teegavarapu, R. S. V., & Chandramouli, V. (2005). Improved weighting
methods, deterministic and stochastic data‐driven models for estimation of missing precipitation records. Journal of Hydrology, 312(1–4),
191–206. https://doi.org/10.1016/j.jhydrol.2005.02.015
Thiessen, A. H. (1911). Precipitation averages for large areas. Monthly
Weather Review, 39(7), 1082–1084. https://doi.org/10.1175/1520‐
0493(1911)39<1082b:PAFLA>2.0.CO;2
Wackernagel, H. (2003). Multivariate geostatistics: An introduction with
applications (3rd ed.). Berlin: Springer‐Verlag.
Webster, R., & Oliver, M. A. (2007). Geostatistics for environmental scientists
(2nd ed.). Chichester, United Kingdom: John Wiley & Sons.
Woldemeskel, F. M., Sivakumar, B., & Sharma, A. (2013). Merging gauge
and satellite rainfall with specification of associated uncertainty across
Australia. Journal of Hydrology, 499, 167–176. https://doi.org/
10.1016/j.jhydrol.2013.06.039
Yang, X., Xie, X., Liu, D. L., Ji, F., & Wang, L. (2015). Spatial interpolation of
daily rainfall data for local climate impact assessment over greater
Sydney region. Advances in Meteorology, 2015, 563629. https://doi.
org/10.1155/2015/563629
Yilmaz, A. G., & Perera, B. J. C. (2014). Extreme rainfall nonstationarity
investigation and intensity–frequency–duration relationship. Journal of
Hydrologic Engineering, 19(6), 1160–1172. https://doi.org/10.1061/
(ASCE)HE.1943‐5584.0000878
Yilmaz, A. G., Hossain, I., & Perera, B. J. C. (2014). Effect of climate change
and variability on extreme rainfall intensity–frequency–duration relationships: A case study of Melbourne. Hydrology and Earth System
Sciences, 18, 4065–4076. https://doi.org/10.5194/hess‐18‐4065‐
2014
Yu, M. C. L., Cartwright, I., Braden, J. L., & de Bree, S. T. (2013). Examining
the spatial and temporal variation of groundwater inflows to a valley‐
to‐floodplain river using 222Rn, geochemistry and river discharge: The
Ovens River, Southeast Australia. Hydrology and Earth System Sciences,
17(12), 4907–4924. https://doi.org/10.5194/hess‐17‐4907‐2013
How to cite this article: Adhikary SK, Muttil N, Yilmaz AG.
Cokriging for enhanced spatial interpolation of rainfall in two
Australian catchments. Hydrological Processes. 2017. https://
doi.org/10.1002/hyp.11163