  We propose a method for estimating the posterior distribution of a standard
geostatistical model. After choosing the model formulation and specifying a
prior, we use normal mixture densities to approximate the posterior
distribution. The approximation is improved iteratively. Some difficulties in
estimating the normal mixture densities, including determining tuning
parameters concerning bandwidth and localization, are addressed. The method is
applicable to other model formulations as long as all the parameters, or
transforms thereof, are defined on the whole real line, $(-\infty, \infty)$. Ad
hoc treatments in the posterior inference such as imposing bounds on an
unbounded parameter or discretizing a continuous parameter are avoided. The
method is illustrated by two examples, one using digital elevation data and the
other using historical soil moisture data. The examples in particular examine
convergence of the approximate posterior distributions in the iterations.
