.


42 http://journals.cihanuniversity.edu.iq/index.php/cuesj CUESJ 2019, 3 (1): 42-49

ReseaRch aRticle

Prediction the Groundwater Depth using Kriging Method and 
Bayesian Kalman Filter Approach in Erbil Governorate
Kurdistan Ibrahim Mawlood, Paree Khan Aabdulla Omer

Department of Statistics, College of Administration and Economics, Salahaddin University-Erbil, Kurdistan Region - F.R. Iraq

ABSTRACT

The aim of this research is using the kriging method as one of geostatistics interpolation methods on the measured value of the specific part 
and Bayesian Kalman filter to identifying the depth of Groundwater in Erbil. Geostatistics is a tool which is developed for statistical analysis 
of any continuous data that can be measured at any location in the space. The Kalman filter is the Bayesian optimum solution to the problem 
of estimating the unknown state of a dynamic system from noisy data and is more efficient than computing the estimate directly from the 
entire past observed data. The main goal of this work is to predict anew value at the unmeasured location by kriging method and Bayesian 
Kalman filter and compare these two methods. The dataset is the observed values of the (295) wells that had been taken from a known 
specific place which called Shaqlawa – in Erbil Governorate. The comparison was done by calculating mean absolute error (MAE) and root 
mean square error (RMSE) for the value of the depth of groundwater in the eara of the study. The values of (MAE and RMSE) of each models 
are compared and the smaller values of them are the better interpolation as it shown in analyzing to evaluate the precision of the prediction.

Keywords: Bayesian estimation, covariance function, Gaussian random field, groundwater-surface interpolation, Kalman filter, kriging 
interpolation method

INTRODUCTION

The Kriging is a spatial interpolation Geostatistical method used for the first time in meteorology, geology, environmental sciences, agriculture, and others fields. This method is used 
to find the best estimator under the assumption of the second 
order stationarity. Geostatistics is a set of tools and models that are 
developed for statistical analysis of any continuous data that can 
be measured at any location in the space. Verify three data feature 
in statistical continuous data analysis: Dependency, stationery, and 
distribution. With these features, you can proceed to the modeling 
of the geostatistical data analysis like simple kriging. In addition, 
the goal of this work is to predict a new value at the unmeasured 
location by Gaussian Semivariogram function (model) and 
compare the results of this model based on the simple kriging 
method and understanding their spatial variability with another 
approach called Bayesian Kalman Filter.

Kalman filter is an optimal linear estimator which provides 
the estimation of signals in noise. Kalman used the state transition 
models for the dynamic system in the estimation process.

GROUNDWATER-SURFACE 
INTERPOLATION

In general, spatial statistics and geographic information system 
(GIS) rely on each other in many ways. Arc GIS is software 
which can be used to create covariates for inclusion in all 
statistical models and to bring out the results from statistical 

models. The work in this study might be very important to 
evaluate the results from different models in simple kriging 
interpolation approaches. This kind of comparison presents a 
relevant meaning for the variability of a physical model which 
used as a reference to validate the interpolation results.

Groundwater depth of the wells of any sample point data 
for generates surface need to be evaluated and preprocessed 
before interpolation. The locations and values of sample point 
data will impact the interpolation result. First, all the collected 
data should come from the same type of wells in the same 
aquifer. The well information should be carefully evaluated 
to make sure the data reflect the dynamics of groundwater 
in the target aquifer, not other aquifers. Second, the spatial 
distribution of sample point data should be carefully 
considered. The clustered data and sparse data in one area 

Corresponding Author: 
Kurdistan Ibrahim Mawlood, Department of Statistics, College 
of Administration and Economics, Salahaddin University-Erbil, 
Kurdistan Region - F.R. Iraq. 
E-mail: kurdistan.mawlood@su.edu.krd

DOI: 10.24086/cuesj.v3n1y2019.pp42-49

Copyright © 2019 Kurdistan Ibrahim Mawlood, Paree Khan Aabdulla Omer. This is 
an open-access article distributed under the Creative Commons Attribution License.

Cihan University-Erbil Scientific Journal (CUESJ)

Received: Apr 27, 2019 
Accepted: May 29, 2019 
Published: Jun 30, 2019

https://creativecommons.org/licenses/by-nc-nd/4.0/


Mawlood and Omer: Prediction the groundwater depth

43 http://journals.cihanuniversity.edu.iq/index.php/cuesj CUESJ 2019, 3 (1): 42-49

will cause different interpolation results. After collecting 
suitable data for interpolation, raster calculator in GIS is used 
to calculate the elevation of groundwater table in each well by 
dividing the surface elevation (digital elevation model) with 
the groundwater depth.[1]

GAUSSIAN RANDOM FIELD

Used Gaussian Random Field Z(s) to identify the spatial 
correlation structure, any field of spatial is a set of random 
variables which parameterized by some set (D⊂Rd). The 
simplest stochastic process form is as follows:

Z s s D Rd( ) :{ ∈ ⊂ }  (1)
Where:
Z(s): Field of random spatial.
s: Coordinates spatial random variable.
D: Domain of spatial random variable.
Rd: d-dimensional Euclidean space.

Any finite collection{Z(s
1
),Z(s

2
).,Z(s

k
).,} is multivariate 

normal:

1 1( ) ( )

. .
~ , ( ( ), ( ))

. .

( ) ( )

i j

k k

Z s s

N Cov Z s Z s

Z s s

µ

µ

    
    
           
            

If the following assumptions hold, then a spatial random 
field is called second-order stationary: [2]

( ( ))E Z s s Dµ= ∀ ∈  (2)

 (3)

Where:
E: Expected value. Cov: Covariance function of two locations.
h=s

2
−s

1
: Vector distance between Z(s

1
) and Z(s

2
).

However, we know that the covariance function can be 
expressed as follows:

( ) ( )
( )( )
[ ]

1 2 1 1

1 1 1 1

2
1 1

( ), ( ) ( ), ( ) ( )

( ) ( )

( ) ( )

s s h

Cov Z s Z s Cov Z s Z s h C h

E Z s Z s h

E Z s Z s h

 


+

 = + = 
 = − + − 

= + −
 (4)

When the expected value is equal to zero.

( ( )) 0,E Z s s D = = ∀ ∈  (5)

Then,

Cov Z s Z s Cov Z s Z s h

C h E Z s Z s h

( ), ( ) ( ), ( )

( ) ( ) ( )

1 2 1 1

1 1

( ) = +( ) 
= = +[[ ]  (6)

When the mean of a second-order stationary of spatial 
random field is equal to zero over the D and the covariance 
function of the locations does not depend on s

1
 and s

2
 but the 

vector h.

And for Cov h Cov Z s Z s h V Z s( ) ( ), ( ) ( )= = +( )  = [ ]0

If variance or covariance function does not exist, the 
intrinsic hypothesis and the spatial random field are called 
(intrinsic stationery) if the following assumption holds as:

E Z s or E Z s Z s s D

V Z s Z s s s s

( ( )) ( ) ( )

( ) ( ) ( )

= −[ ] = ∀ ∈

−[ ] = − ∀

µ

γ

1 2

1 2 1 2 1

0

2 ,, s D2 ∈
 (7)

Where,
V: Variance, γ: Semivariogram, 2γ: Variogram

2 1 2 1 1

1 1

2

 = −[ ] = − +[ ]
= + −[ ]( ) −
V Z s Z s V Z s Z s h

E Z s h Z s E Z

( ) ( ) ( ) ( )

( ) ( ) (ss h Z s

E Z s h Z s

1 1

2

1 1

2

+ −[ ]( )
= + −[ ]( )

) ( )

( ) ( )

 (8)

In addition, if the covariance function C(s
1
−s

2
)=C(h) 

or semivariogram γ(s
1
−s

2
)=γ(h) depends only on separation 

distance between s
1
 and s

2
, h=||s

1−
s

2
||, then the spatial 

random field is called isotropic.[3]

The Kriging

Kriging is a spatial interpolation geostatistical method used 
for the 1st time in meteorology, geology, environmental 
sciences, agriculture, and other fields. This method is 
used to find the best estimator under the assumption of 
the second-order stationarity. A geological process may 
not be stationary in reality. In case, where the process is 
non-stationary, we could use non-linear functions.[4] Here, 
the classifying Geostatistical techniques are as follows 
[Table 1].

Assumptions

In some of the interpolation methods especially 
Geostatistical methods, they have their own assumptions 
as follow:[5]

1. Stationarity
2. Intrinsic hypothesis
3. Isotropy and anisotropy
4. Unbiased.

Semivariogram Function and Covariance 
Function

The semivariogram function is a structure of intrinsically 
stationary spatial random field which describes a broader 
class of Erath phenomena. In the case of the second-order 
stationary spatial random processes, there is an equivalence 
between covariance function and semivariogram function as 
follows:

Table 1: Classification of geostatistical techniques

Model Stationary Nonstationary

Linear Ordinary/simple kriging Universal kriging; 
kriging using IRF-K

Nonlinear Disjunctive kriging 
simulation

Simulation of IRF-K

IRF-k: Intrinsic random functions of order k

( ) ( ) ( )1 2 1 1 2 2
1 2 1 2

( ), ( ) ( ) ( )

( ) ,
s sCov Z s Z s E Z s Z s

C s s s s D

µ µ = − − = 
− ∈



Mawlood and Omer: Prediction the groundwater depth

44 http://journals.cihanuniversity.edu.iq/index.php/cuesj CUESJ 2019, 3 (1): 42-49

V Z s h Z s V Z s h V Z s

Cov Z s h Z s

V Z s

( ) ( ) ( ) ( )

( ), ( )

( )

+ −[ ] = +[ ] + [ ]
− +[ ]
= [
2

2 ]] − +[ ]
= −[ ] = ÷

= −

2

2 0 2 2

0

Cov Z s h Z s

C C h h

h C C h

( ), ( )

( ) ( ) ( )

( ) ( ) ( )




 (10)

The semivariogram function is a measure of dissimilarity 
between pairs of locations or observed value Z(s+h) and Z(s). 
There are three parameters of semivariogram function for the 
spatial second-order stationary processes: Nugget effect (c

0
), 

range (a), and partial sill (c) which are shown in Figure 1. 
The sum (c

0
+c) is called the sill (Marcin and Marek, 2010), 

(Sluiter, 2009).

The Gaussian Semivariogram Model

The Gaussian model is commonly used to represent events with 
a small scale spatial structure,[6] the equation for this model is 
similar to the normal cumulative distribution function, and it 
is given by [Figure 2]:

 z h C C
h
a

h( ) exp ,= + −
−

















>0
2

0
2

1 0  (11)

The Simple Kriging Interpolation Method

Simple kriging can deliver the value at any unmeasured 
location (s

0
) using a linear estimator λ

i
 for the measured 

values at locations (s
1
,s

2
.,s

n
):

Z s Z si i
i

n

i
i

n

( ) ( ),0
1 1

1= =
= =
∑ ∑   (12)

As the sum of all linear estimators equals one, the 
unbiased prediction is ensured. Using the intrinsic hypothesis, 
the estimation variance is calculated by the formula:

( )2 20 0 0 0
0

1 1 1

ˆ( ) ( ( ) ( ( ) ( ))

2 ( ) ( )

e

n n n

i i i j i j
i i j

s Var Z s E Z s Z s

s s s s



    
= = =

= = −

= − − −∑ ∑∑
 (13)

The semivariance (s
j
−s

0
) and(s

i
−s

0
) are taken from the 

variogram model. The goal is to minimize σ2(s) under the 
unbiased conditions and to find the corresponding weight.

The weights are chosen to fulfill: 0 0ˆ( ( ) ( )) 0E Z s Z s− =
according to the assumption that E[Z(s

1
)−Z(s

2
)]=0 which 

called intrinsic stationery and to solve this the Lagrange 

multiplier (μ) is introduced, not easy solving this model by any 
other ways:

Z s s s s s so i i j
i

n

j( ) ( ) ( ) ( )= − + = −
=
∑λ γ µ γ

1
0 0  (14)

The solutions to this linear equation system are the values 
of the linear estimators. This system is called kriging system.[7] 
The kriging variance is determined by:

σ λ γ µk i i
i

n

s s s2 0
1

0= − +
=
∑ ( ) ( )  (15)

KALMAN FILTERING

The Kalman filter gives unbiased, linear, and minimum 
variance recursive estimate to the state of a dynamic system 
from noisy data taken in separate real-time.[10] KF used in the 
formulation of a dynamic model for linear dynamical systems 
provides a recursive solution to the problem of the optimal 
linear filter. The recursive solution in every state update 
estimate is calculated from previous estimates and data entry; 
only the previous estimate requires storage. KF has been widely 
used in the fields of signal processing and modern control and 
airborne surveillance systems, radar signal processing and 
control, and adaptive controls.

Here, we will provide only the equations needed to 
develop discrete recursion KF, discrete state equations are 
given as linear dynamic for signal θ

t
 and observation Y

t
.[11] We 

have two equations:

Figure 1: The relationship between the γ(h) and C(h) for the second-order stationary of a spatial random variable

Figure 2: The Gaussian semivariogram model with parameters



Mawlood and Omer: Prediction the groundwater depth

45 http://journals.cihanuniversity.edu.iq/index.php/cuesj CUESJ 2019, 3 (1): 42-49

1. Observation equation: Y Ft t t t= +θ υ

2. System equation: θ θ ωt t t tG= +−1
The initial state vector has a mean θ

0
 and a covariance 

matrix P
0
 that is:

0 0 0 0
ˆ[ ] = and [ ] = E Var P  

Y
t
 is the vector of noisy observation at time t, of dimension 

m×1.
F

t
 is the known matrix at time t, of dimension m×p relates the 

state vector to the observation vector Y
t
.

υt is the vector of observation error (stochastic error term) 
at time t, having normal distribution with zero mean 
and covariance V

t
, of dimension m×θt is the vector of 

unknown state parameters (signal) at time t, of dimension 
p×1, system states.

G
t
 is the matrix of time-varying state transition matrix 

dimension pxp.
ωt is the system error (stochastic error term) at time t, having 

normal distribution with zero mean and covariance W
t
, of 

dimension p×1.
Where, V

t
 and W

t
 are variance-covariance measurement and 

process noise matrices, respectively. The noise vectors ω
t
 

and υ
t
 process and measurement noise, respectively, are 

uncorrected white noises (have zero means) they are 
independent of each other and have normal probability 
distributions.
The disturbances υ

t
 and ω

t
 uncorrected with the initial 

state that as:

E( ) = 0υ θt
T
0  and E( ) = 0ω θt

T
0  for t=1, 2, 3,…,T

Bayesian Estimation of Dynamic Linear 
Models

Kalman filter is used for estimating a state of a dynamic 
system which gives the best-unbiased estimator using previous 
measurement knowledge. Previous algorithms stand from the 
use of all the previous information to estimate the state of the 
system at the time of the next step.

This dynamic linear model is the simple model 
represented as:

Observation equation Y
t
=μ

t
+ε

t

System equation μ
t
=μ

t−1+
σ

μ2

Where the errors ε
t
 and ω

t
 are mutually independent and 

independent of initial information probability (μ
0
|D

0
)

Initial information probability (μ
0
|D

0
) represents 

forecaster’s probabilistic information of about the μ
0
 at time 

t=0. The mean m
0
 and the variance c

0
 are a point estimate 

and the associated uncertainty of μ
0
. D

t
 consists all the 

information available up until time t, including D
0
, the values 

variances values { , : }σ σε µ
2 2 0t > , and the the observations 

Yt,Y
t−1

,…,Y
1
 or, the only new information becomes available 

at any time t is the observed value Y
t
, where D

t
={Y

t
,D

t−1
}.

Here, we start expressing initial information concerning 
the parameter μ

0
 was described in the form of a normal 

probability distribution with mean m
0
 and variance c

0
 is:

0 0 0~ ( , )N m c

Using a mathematical process and Bayes’ theorem at the 
time (t−1) of the parameter μ

t−1
 is:

1 1 1 1( | ) ~ ( , )t t t t tD N m c − − − −

Where, D
t−1

 is D Yt
t

−
−=1

1 2 2{ , , }σ σε µ .

We aim now to find a final distribution for the parameter 
μ

t
 at last time is the time appointed by our knowledge 

of all available data, or the time that we want to filter. We 
recompense an extremely presence of the information is:

D Y Dt t t= −{ , }1 .

Then Bayes’ theorem is

Posterior ∝ Observed likelihood×Prior

P D P Y P Dt t t t t t( | ) ( | )* ( | )µ α µ µ −1  (16)

Where

P(Y
t
|μ

t
) represents the likelihood distribution function at 

time t

P(μt|D
t−1

) represents the prior distribution at time t

Both of them are distributed normally

Now to find all of the weighting function P(μ
t
|D

t−1
) and 

the prior probability before observation Y
t
, using dynamic 

linear model and probability distribution as:

E Y E

E E
t t t t t

t t t t

t

( | ) ( | )

( | ) ( | )

µ µ ε µ
µ µ ε µ

µ

= +
= +
=

And

V Y V

V V
t t t t t

t t t t

( | ) ( | )

( | ) ( | )

µ µ ε µ
µ µ ε µ

σε

= +
= +

= 2

Or

( | ) ~ ( , )Y Nt t tµ µ σε
2

And from the system equation:

E D E

E D E D

m

t t t

t t t

t

t

t

( | ) ( )

( | ) ( | )

µ µ σ

µ σ
µ

µ

− −

− −

−

= +

= +

=

1 1

1 1

1

And

V D V

V D V D

c R

t t t

t t t

t t

t

t

t

( | ) ( )

( | ) ( | )

µ µ σ

µ σ

σ

µ

µ

µ

− −

− −

−

= +

= +

= + =

1 1

1 1

1
2

Or

1 1( | ) ~ ( , )t t t tP D N m R − −
Then by substituting each of (Yt|μt), P(μt|D

t−1
) in 

Equation (2.10.1) we get



Mawlood and Omer: Prediction the groundwater depth

46 http://journals.cihanuniversity.edu.iq/index.php/cuesj CUESJ 2019, 3 (1): 42-49

P
t
D
t

e Y
t t R

t
t

m
t

( | ) [ ( ) ( )]µ α
σ
ε

µ µ− − + −
−

1
2

1
2

1
1

 (17)

We can find that:

P
t
D
t

e
c
t

t
m
t

( | ) [ ( )]µ α µ− −
1
2

1

 (18)

Now we get mean m
t
 and variance ct of final distribution, 

respectively, are:

m
Y

R
m

R

c
R
Rt

t
t

t

t
t

t

=
+

+
=

+
+

−
1 1

1 1

2 1

2

2

2

σ

σ

σ
σ

ε

ε

ε

ε

, and 

Then, the final distribution for parameter μt at time t is

( | ) ~ ( , )t t t tP D N m c

We can write the mean and variance of posterior 
probability distribution after simplifying them as:

m m K Y mt t t t t= + −− −1 1( )  (19)

c Kt t= *σε
2  (20)

Where

K R Rt t t= +
−( )σε

2 1  (21)

Where the Equation (19) called Kalman filter, sometimes 
called for two Equations (19), (20) together Kalman filter, and 
Equation (21) is Kalman gain.

Forecasting Dynamic Linear Models

We mean by forecasting finding the h-step ahead for each 
distribution given the data dynamic linear model y1,y2.,yt. 
For a dynamic linear model, the h-step-ahead forecasting 
distributions, for states and observations, are obtained as a 
product of the Kalman filter.[12,13] From the observations y

1:t
 to 

Y
t
+h we noticed that the data Y

1:t
 provide information about 

θt, which, in turn, gives information about the future state 
evolution up to θ

t+k
 and consequently on Y

t+h
.[13,14]

Then, the general dynamic linear model for h step is:

Observation equation:

Y F N Vt h t h t h t h t t+ + + += +θ υ υ ~ [ , ]0  (22)

System equation:

θ θ ω ωt h t h t h t h t tG N W+ + + − += +1 0~ [ , ]  (23)

Let a0=mt, R0=Ct, then for h≥1, for forecasting at time t, 
the forecaster requires the h-step ahead marginal distributions, 
(Y

t+h
|Dt) and (θ

t
+h|Dt), predictions start from the posterior 

estimate of the state, obtained from the Kalman filter, on the 
time (t) on which the forecast is made.

The mean and variance for the observation equation are:

1 1 1

ˆ ( | )

( | )
t h t h t

t t t t

k
t

Y E Y D

E F D

FG m

 
+ +

+ + +

=
= +

=

 (24)

Q t Var Y D

FG C G F FG W G F V

t t h t

h
t

k T T i i T T

i

h

( ) ( | )

( ) ( )

=

= + +

+

=

−

∑
0

1  (25)

APPLICATION, RESULTS AND DISCUSSION

The Study Area and Data Collection

Spatial variability of any regional variable is a result of complex 
processes which is working at the same time and over long 
periods of time. Variation of a regional variable has never been 
an easy task or work. Many regional random variables vary not 
only horizontally but also with depth such as wells of water, 
oil, and gas. The spatial attribute includes approximately all 
wells at fields that exist which spatially distributed across the 
study area to represent the fluctuations of levels from place to 
another place in whole area. 

The source of dataset is the observed values of the 295 
wells that had been taken from the known specific place which 
called Shaqlawa in Erbil Governorate. These observations 
are collected by GPS by the ministry of agriculture, Figure 3 
expresses the location of each well, and each point of the 
dataset has its name and properties. Furthermore, each point 
has its goal for drilling the well, some of them for drinking or 
irrigation and agriculture.

Figure 3 shows the geographical location of the data 
points of Shaqlawa or all wells. They are shown as a surface 
which can be represented by the most probably prediction 
map and also by estimated prediction when we are applying 
a model of the simple kriging interpolation. The figure shows 
cycles that are many points closer together tend to be more 
alike than things that are farther apart (quantified here as 
spatial autocorrelation).

Results of Groundwater-Surface 
Interpolation

Figure 4 can be described as continuous data and represented 
the random field. It shows a histogram of the observed 
values and how they distributed in the region. It represents 
the curve of the observations values, and the results after 
taken the suitable transformation that has a small standard 
deviation equal to 0.37924, the skewness of dataset is equal 
to −0.61614, near to zero and the kurtosis of the dataset is 
equal to 4.5715 not near to three. These indicated that the 
distributions of depth in sample data points were approximate 
to normal distribution.

Results of Gaussian Semivariogram 
Simple Kriging Surface Interpolation

Spatial dependency can be detected in this dataset using 
several tools available in geostatistical analysis exploratory 
spatial data analysis and geostatistical wizard in GIS software. 

Figure 3: Geographical location of the dataset



Mawlood and Omer: Prediction the groundwater depth

47 http://journals.cihanuniversity.edu.iq/index.php/cuesj CUESJ 2019, 3 (1): 42-49

In geostatistics semivariogram is called spatial modeling 
(called structural analysis or variography).

Table 2 shows the model called the Gaussian 
semivariogram. This gives a quantitative description of 
variability about the same regional variation. The important 
part of the variogram is the rang which describes the distance. 
It shows the isotropy of the Gaussian semivariogram model 
for (295) dataset of groundwater of wells in Shaqlawa, 
which all points have equal directional of variability as being 
north-south and east-west. Furthermore, the anisotropy of 
Gaussian semivariogram model expresses the same data 
(554) that change with the direction which described an 
ellipsoid, and this ellipsoid specified by the length of two 
orthogonal axis (Major and Minor) with its orientation 
Angle θ.

Figure 5 shows Gaussian semivariogram model for 
isotropy and anisotropy surfaces, these used to predict of 
random field (spatial field) of Groundwater well Z(s

0
) unknown 

value in the same area, which depends on only maximum (5) 
neighbors value of the measured values. Table 2 contains all 
cross-validation of isotropy and anisotropy of simple kriging 
Gaussian semivariogram model and the prediction value of the 
depth of unknown new value of well.

Table 3 shows all information about simple kriging by 
Guassian semivariogram model to predict a new location of 
groundwater in the same area for both of isotropy and anisotropy. 
The prediction value of unknown is measured also depending 
on (5) maximum and (2) minimum value of neighbors from 
measured value for the depth of new well with longitude and 
latitude. In anisotropy Gaussian semivariogram model is equal 
to 173.2582 and its greater than the depth of the isotropy 
exponential semivariogram model which equal to 132.3405.

Estimation the Depth of Groundwater 
using Kalman Filter

The estimation process of parameters is made using the initial 
values of the depth of groundwater mean and variances are 
being of an estimate of it is values. Equations (19), (20), and 
(21) are used to estimate parameters, Error Covariance, and 
Kalman Gain for univariate simple dynamic linear model for 
each observation of the groundwater depth.

Table 4 gives the actual value, estimated groundwater depth 
parameter, Error Covariance, and Kalman Gain. Examining 
the table reveals the following: The Error Covariance is 
convergence at a time point (t=18) with the value (802.1407). 
Kalman Gain is convergence at a time point (t=16) with the 
value (0.383983). There is two convergences time and the last 
time convergence is the best estimate for the same time.

A compression between groundwater depth resulting 
from the estimate and forecasting groundwater depth values 
for all observations is shown in Figure 6. The results show how 

Figure 6: Difference between estimate values and forecast values for 
groundwater depth

Figure 4: The Histogram of the dataset

Figure 5: (a and b) Surfaces of the Gaussian semivariogram model

a

b



Mawlood and Omer: Prediction the groundwater depth

48 http://journals.cihanuniversity.edu.iq/index.php/cuesj CUESJ 2019, 3 (1): 42-49

closely the estimated model matches the forecasting electricity 
groundwater depth. The blue line is estimate values and the 
red line is forecasting values.

Evaluating Accuracy (or Error) Measures 
for Estimation

For selections, the best approach (method) for prediction the 
groundwater depth two accuracy criteria were used, which 
are (mean absolute error [MAE] and root mean square error 
[RMSE]). They are based on the error estimate, which is the 
difference between the estimated and forecast groundwater 
depth values. The results are shown in Table 5.

Table 3 represents that Kalman Filter is the better 
approach based on the two criteria, where each of them has a 
small value than simple kriging using Gaussian semivariogram 
model for both of isotropy and anisotropy.

CONCLUSIONS

The major results of performing the analysis of two approaches 
(Simple Kriging and Kalman Filter), the following main 
conclusions have been achieved:

1. The data follow approximately to normal distribution 
after taken log transformation for variability and second-
order stationary to remove the trend.

2. The variability of the depth of groundwater wells elevation 
in the region changes from north to south or from west-
north to south-east, this is due to the spatial dependency 
in the area which is one of the reasons of the depth of 
well.

3. For unmeasured value used simple kriging with Gaussian 
semivariogram model, in this model the predicted value 
by isotropy semivariogram model is better than the 
anisotropy semivariogram model depending on the value 
of the depth of groundwater [Table 2].

4. When applied the Kalman filter for estimating and 
forecasting the depth of groundwater, we found that 
the KF is suitable for estimation and yield good results, 

because in the presence of white Gaussian noise because 
when the normality condition holds, the KF will give 
optimum and sufficient results.

5. The result shows that the Bayesian Kalman filter approach 
has the best estimation and forecasting results compared 
to simple kriging using the accuracy criteria (MAE and 
RMSE).

Table 2: Results of gaussian semivariogram model

Simple kriging exponential 
semivariogram

Isotropy Anisotropy

Nugget (C
0
) 0.0714 0

Rang (α) 2 2

Major rang 2554.326 4998.32

Minor rang 2554.326 1949.16

Partial sill (C) 0.081173 0.093287

Lag size (h) 416.5267 416.5267

Number of lag 12 12

Table 3: Accuracy criteria for each approach

Methods MAE RMSE

Kalman filter 13.2975 18.078

Simple kriging/gaussian model 29.7364 43.719

MAE: Mean absolute error, RMSE: Root mean square error

Table 4: Results of simple kriging gaussian semivariogram model

Simple kriging exponential 
semivariogram

Isotropy Anisotropy

Maximum neighbors 5 5

Minimum neighbors 2 2

Predicted value

X 400,739 400,739

Y 421,531 421,531

Depth (m) 132.3405 173.2582

Table 5: Estimated parameters when ( ̂ 0 =126.02, 
P̂0  =2089.088, and W=500)

Convergence 
time

Actual 
values

Estimate Error 
covariance

Kalman 
gain

1 189 126.02 2089.088 0.73123

2 201 167.5177 1156.157 0.55345

3 91 133.6806 923.7829 0.442213

4 142 137.0526 846.7026 0.405315

5 100 122.529 818.8316 0.391973

6 130 125.4203 808.4435 0.387

7 160 138.7378 804.528 0.385126

8 140 139.223 803.046 0.384416

9 50 104.9482 802.4841 0.384147

10 151 122.6342 802.271 0.384045

11 183 145.8151 802.1901 0.384007

12 100 128.2224 802.1594 0.383992

13 100 117.3854 802.1478 0.383986

14 220 156.7878 802.1434 0.383984

15 193 170.6927 802.1417 0.383984

16 180 174.2665 802.1411 0.383983

17 200 184.1478 802.1408 0.383983

18 147 169.8836 802.1407 0.383983

19 202 182.2158 802.1407 0.383983

20 151 170.2295 802.1407 0.383983

. . . . .

. . . . .

291 239 164.2836 802.1407 0.383983

292 211 182.2219 802.1407 0.383983

293 73 140.2825 802.1407 0.383983

294 91 121.3589 802.1407 0.383983

295 82 106.2469 802.1407 0.383983



Mawlood and Omer: Prediction the groundwater depth

49 http://journals.cihanuniversity.edu.iq/index.php/cuesj CUESJ 2019, 3 (1): 42-49

REFERENCES
1. S. I. Gou. “Identifying Groundwater Dependent Ecosystems in the 

Edwards Aquifer Area, Zachry Department of Civil Engineering”. 
Texas: Texas A and M University, 2010.

2. R. Webster and M. Olivier. “Geostatistics Environmental Scientists”. 
2nd ed. New York, USA: John Wiley and Sons, 2007.

3. L. Marcin and K. Marek. “Simple Spatial Prediction Least Squares 
Prediction, Simple Kriging and Conditional Expectation of Normal 
Vector”. Poland: Department of Geomatics, AGH University of 
Science and Technology in Krakow, 2010.

4. D, Sarma. “Geostatistics with Application in Earth Science”. 2nd ed. 
Dordrecht: Springer, Captial Publishing Company, 2009.

5. R. Sluiter. “Interpolation Methods for Climate Data; Literature 
Review”. KNMI, R and D Information and Observation 
Technology, 2009.

6. I. Clark and W. V. Harper. “Practical Geostatistics”. Columbus, 
Ohio: Ecosse North America, 2000.

7. M. Disse and K. Amin. “Review of Various Methods for 
Interpolation of Rainfall and their Applications in Hydrology, 

Thesis of Chair of Hydrology and River Basin Management”. 
Munich: Technical University Munich, 2017.

8. R. E. Kalman. “A new approach to linear filtering and prediction 
problems”. Journal of Basic Engineering, vol. 82, pp. 35-45, 1960.

9. R. E. Kalman and R. S. Bucy. “New linear filtering and prediction 
problems”. Journal of Basic Engineering, vol. 83, pp. 95-108, 1961.

10. R. Eubank. “A Kalman Filter Primer”. United States of America: 
CRC Press Publishing, 2006.

11. B. D. Anderson and B. G. Moore. “Optimal Filtering”. New 
South Wales, Australia: Prentice-Hall Publishing, University of 
Newcastle, 1979.

12. M. Meng, D, Niu D and W. Sun. “Forecasting monthly electric 
energy consumption using feature extraction”. Journal Energies, 
vol. 4, pp. 1495-1507, 2011.

13. M. West and J. Harrison. “Bayesian Forecasting and Dynamic 
Models”. 2nd ed. New York: Springer-Verlag, 1997.

14. R. Gentleman, K. Hornik and G. Parmigiani. “Dynamic Linear 
Models with R”. New York: Springer Science and Business Media 
Publishing, 2009.