are these methods appropriate for four dimensional modeling?

I am working on a thesis in soil science with four dimensional modeling 
as a major component.  I am using R for the interpolation of data and 
VisIT to render the data as volumes and as volumes changing in time.  I 
am using Applied Spatial Data Analysis with R (Bivand, Pebesma, 
Gomez-Rubio) as a reference.  I've also been using Chatfield's Analysis 
of Time Series, less for time series specifically than for helping to 
understand the statistical concepts.

My data are soil property data. A purely spatial component comprises 
nine soil physical properties taken across a 12 ha agriculture field. 
There are sixty locations with five measurement depths at each location 
(3D).  The spatial-temporal component is the measure of volumetric soil 
moisture status taken at the same 300 positions.  These data are 
discrete readings taken at unequal intervals.  The readings exist for 
growing seasons over three years.

The observation grid is half regular and half semi-random.  Each of 
thirty regularly placed locations has a satellite location set in random 
proximity but not too far and not too close.  The observation depths are 
at 15 cm intervals starting at -15 cm and extending down to -75 cm.

My first pass through the data was mainly concerned with the mechanics 
of the process.  I used IDW for interpolating each of the soil 
properties through the volume (myriad packages including sp, gstat, and 
others mentioned in Bivand et al), treating each property as independent 
of the others.

I then went after the temporal soil moisture data.  At each of the 300 
positions, I took the set of time-sequenced measures (measurement 
intervals varied from days to weeks during each growing season) and 
interpolated values for days not measured.  I used the Stineman 
algorithm provided by the na.stinterp function from the stinepack R 
package.  Once I had all days for all positions, I used IDW again to 
interpolate volumes for each day.

Each volume was exported as an unstructured point grid for rendering in 
VisIT.  My data analysis was limited to adjusting the IDW weighting 
power such that the density distributions of the interpolated values 
were similar to the density distributions of the observed values 
(eyeballing overlaid plots).

My interpolation grid cell size is defined at 10x10x1, modeling 10 m by 
10 m by 1 meter.  This is a little misleading.  I treat the depth 
dimension as having the same units as the areal dimension, so -15 cm 
becomes equivalent to -15 m.  Two reasons (the second not necessarily 
defensible): I need the depth to be scaled for reasonable visualization, 
and I want to decrease the influence of the depth measures on the layers 
above and below each depth.

In my next iteration through the data, I'd like to use the 
geostatistical techniques presented in Bivand et al.  I have (just) 
started with variograms with the hope of exploring the variogram in the 
context of the volume.  I have run into the problem of dimensionality. 
The plot functions are 2D, so projecting the variogram onto the 
observation spatial volume is not working.  That got me thinking.  If I 
use kriging, will I be using the influence of depth in the interpolation 
model?  Is the variogram created with the influence of depth when the 
spatial structure has three dimensions?  Is a better use to treat each 
depth as a separate layer to be interpolated in two dimensions, 
especially given that soils are often physically layered?  If three 
dimensional kriging is okay, how about four dimensional to move through 
a time series? Or would I be better to stick with the Stineman algorithm 
for the time-series interpolation?

Any feedback on this (including any warnings that I might be abusing 
things with my time series and IDW approaches) would be very much 
appreciated.  I am definitely in learning mode.

Mark