Palaeochannels in Northern New South Wales: inversion of electromagnetic induction data to infer hydrologically relevant stratigraphy.
Abstract: Palaeochannels, or prior streams, are strings of sandier sediments that occur frequently in the irrigated alluvial plains of Northern New South Wales, Australia. These landscape features have been recognised as locations of substantial deep drainage losses, and are therefore target areas for water use efficiency. Electromagnetic induction (EM) measurements were used to identify the width and the depth of the palaeochannel sediments in a 2-dimensional transect. Three different inversion techniques, Tikhonov regularisation, the McNeill layered earth model, and an optimal linear combination of EM measurements, were applied to a combination of EM-38 and EM-34 data. Using various kriging techniques, the resulting conductivity profiles were interpolated to soil property data and transformed to saturated hydraulic conductivities using pedotransfer functions. There were distinct differences in the resulting stratigraphies depending on the inversion and interpolation method employed. Trend kriging of the sampled soil property data using the Cook and Walker and Tikhonov inversion data as a trend surface gave the most consistent hydraulic conductivity values compared to the sampled soil property data. However, differences between inversion and interpolation methods were negated by uncertainties in the pedotransfer functions.

Additional keywords: geophysics, vadose zone hydrology, kriging.
Subject: Geophysics (Research)
Authors: Vervoort, R.W.
Annen, Y.L.
Pub Date: 05/01/2006
Publication: Name: Australian Journal of Soil Research Publisher: CSIRO Publishing Audience: Academic Format: Magazine/Journal Subject: Agricultural industry; Earth sciences Copyright: COPYRIGHT 2006 CSIRO Publishing ISSN: 0004-9573
Issue: Date: May, 2006 Source Volume: 44 Source Issue: 3
Topic: Event Code: 310 Science & research
Product: Product Code: 8523400 Geophysics NAICS Code: 54171 Research and Development in the Physical, Engineering, and Life Sciences
Geographic: Geographic Scope: Australia Geographic Code: 8AUST Australia
Accession Number: 144605694
Full Text: Introduction

Palaeochannels, prior streams or relict river channels, have long been recognised throughout the alluvial plains of eastern Australia (Butler 1950; Stannard and Kelly 1968). Pal~eochannels are aggregated relic landscape features, stemming from ancient river activity, which generally no longer convey water. On the surface, they appear as bands of lighter textured soil, and can often be identified as such from aerial photographs. Sometimes a distinction is being made between features with surface expression (deemed prior streams) and features without surface expression (deemed palaeochannels), but since the meaning of the Greek word 'palaeos' only translates to 'old', this paper will define all these features as palaeeochannels.

The spatial extent and depths of the larger palaeochannel systems have been broadly mapped in Australia (e.g. Page et al. 1996; Young et al. 2002), partly due to their significance in terms of groundwater storage (Rogers et al. 2002) and mineral exploration. But, due to the complex and multi-branched nature of the palaeochannel networks in the alluvial plains, smaller palaeochannels are a common feature in irrigated paddocks (Rogers et al. 2002; Young et al. 2002). Considering that studies have identified palaeochaunels and associated soils as areas with a high risk of deep drainage (e.g. Stannard and Kelly 1968; Rogers et al. 2002; Triantafilis et al. 2003), understanding of the hydrology of palaeochannels is essential for efficient irrigation and effective recharge management.

Generally speaking, 2 main features would govern the hydrology and recharge of palaeochannels: the size (spatial extent and depth) of the feature and the stratigraphy (layering of soil properties) within the profile. Traditionally, the size of palaeochannels has been identified from aerial photographs and the profile stratigraphy has been investigated by deep drilling into the sediments (e.g. Stannard and Kelly 1968). However, drilling is costly, and more rapid non-invasive techniques would be more useful if extensive mapping for hydrological simulation is needed.

Electromagnetic (EM) induction has become increasingly popular as an indirect method for identifying changes in soil properties. Electromagnetic methods measure a depth-weighted average of the bulk soil electrical conductivity (EC), otherwise known as the apparent electrical conductivity (E[C.sub.a]). The bulk soil EC is generally a function of the soil water and salt content, but in low salinity (<4 dS/m) environments, EC is also a function of the texture (clay content), mineralogy, and structure (Friedman 2005).

The horizontal spatial extent of palaeochannels can be successfully identified using various EM instruments (Williams and Hoey 1987; Triantafilis et al. 2003), because there tend to be clear differences in texture and possibly wetness and salt content between the palaeochannel and surrounding sediments. Repeated EM measurements at different coil spacings or different heights above the ground can also be used to identify vertical changes in conductivity through the soil profile (McNeill 1980; Corwin and Rhoades 1984; Cook and Walker 1992; Hendrickx et al. 2002) and this could help identify layering in a profile. While other geophysical methods such as electrical imaging (Acworth 1999) can give very detailed stratigraphy information, the cost, time investment and limited spatial extent of these methods quickly decreases the value at larger scales compared with EM measurements. In addition, inversion methods related to these techniques are often complex.

Most studies employing EM have concentrated on finding linear or optimal combinations of measurements to predict soil property observations (such as EC profiles) using regression techniques (e.g. Corwin and Rhoades 1984; Triantafilis and Lesch 2005). Cook and Walker (1992) extended this approach; they developed an optimal linear combination of depth weighing functions to calculate the ECa for a specific depth using least squares minimisation.

Because E[C.sub.a] is a convolution of the depth response function of the instrument and the soil EC profile, inversion of the data is actually needed to be able to resolve the true EC in the profile (e.g. Hendrickx et al. 2002). Inversion of EM data involves solving the following matrix equation for the unknown soil EC profile ([sigma]) (Borchers et al. 1997):

(1) K[sigma] = d

where d is the vector of apparent conductivities or instrument readings, and K is a matrix of depth response functions of the instrument at different heights (Borchers et al. 1997). Non-linear inversion of geophysical data, such as EM and resistivity, is well established and several commercial software packages are available for developing solutions (e.g. Auken and Christiansen 2004; Monteiro Santos 2004). However, the strength of the EM method is that it is approximately linear (McNeill 1980; Hendrickx et al. 2002), allowing non-numerical inversion schemes, being equally effective as the non-linear approximations and using less computer resources (Hendrickx et al. 2002). McNeill (1980), for example, suggests a linear inversion method based on the layered earth model using the depth response function of the instrument, but this is limited to even-determined problems (Menke 1989). Under-determined linear problems can be resolved using a regularisation method, such as Tikhonov regularisation (Botchers et al. 1997; Hendrickx et al. 2002).

Most of the EM inversion studies have concentrated on finding 1-dimensional EC-depth functions (McNeill 1980; Corwin and Rhoades 1984; Cook and Walker 1992; Hendrickx et al. 2002). Only one study attempted to reconstruct soil layers (horizons) from EM measurements (McBratney et al. 2000). This study indicated that the Tikhonov regularisation method was useful, but that much depended on the regularisation order without a-priori information. Borchers et al. (1997) found, compared with the method developed by Cook and Walker (1992) and regression methods, that the second-order Tikhonov regularisation was best for estimating EC profiles as it predicts smooth changes in EC down the profile.

Apart from inversion of the EM profile, further data manipulation steps need to be undertaken to predict hydrologically important parameters, such as hydraulic conductivity ([K.sub.S]), from the EM signal. While the EM signal has been correlated to soil properties (Triantafilis et al. 2003), this needs to be taken one step further to predict [K.sub.S]. For example, pedotransfer functions, which have been developed as a means of converting soil property data in to hydraulic properties, could be used (Minasny and McBratney 2002).

From this review of literature it appears that different inversion techniques could lead to different stratigraphies due to differences in the mathematical approach to the inversion problem. If the purpose of the EM inversion were to describe the hydrology of a palaeochannel for management or simulation purposes, this would lead to different hydrological behaviours, due to different hydraulic properties of the profile inferred from soil properties contained in the different stratigraphies. The main aims of this study therefore are:

(i) identifying whether a combination of soil coring and EM measurements, subsequent inversion interpolation will assist in developing a more complete picture of the stratigraphy of palaeochannels; and

(ii) investigating whether a combination of the method of inversion and interpolation influences the final stratigraphy and hydraulic properties of the system.

Materials and methods

The study site was located on a property in the lower Gwydir River valley (north-western New South Wales, Australia). The field is generally used for furrow-irrigated cotton but was planted in a rotational wheat crop during the time of the study. Waterlogging caused by the formation of a perched water table associated with a palaeochannel has been a persistent problem in this field (Triantafilis et al. 2003).

A 370-m transect was set up from the edge of the field (Peg 1) directly across the palaeochannel (Fig. 1). The transect was not aligned with the furrows, and was made long enough to extend beyond the area affected by waterlogging. Pegs were placed at 10-m spacings along the transect, with Peg 1 at 0 m (-29.272[degrees]N, 149.743[degrees]E) and Peg 38 at 370 m (-29.274[degrees]N. 149.746[degrees]E). The precise extent to which the palaeochannel sediments occurred, both laterally across the field and vertically down the profile, was initially not known. Drilling to obtain soil samples was performed at 6 locations (Pegs 1, 6 13, 17, 25, 38) along the transect to a depth of 13 m, except at Peg 25 where a depth of 15 m was reached. The selection of these sampling sites was mainly based on a visual interpretation of aerial photography and maps generated in earlier work (Triantafilis et al. 2003). Soil was sampled from each metre and a small subsample was sealed in a tin in order to determine the soil wetness. Bulk soil samples were air-dried for approximately 3 weeks. A fraction of each soil sample was ground and sieved to 2 mm. Particle size analysis was performed for each sample using the pipette method (Gee and Bauden 1986) to determine the clay (<2 [micro]m), silt (2-20 [micro]m), fine sand (20-200 [micro]m), and coarse sand (200-2000 [micro]m) fractions.

The gravel fraction of any samples containing gravel was determined by weighing approximately 100 g of the original uncrushed soil sample and then wet sieving the soil through 2 sieves to determine the fine (2-4.75 mm) and coarse (>4.75 mm) gravel fractions. The sieved gravel was oven-dried at 105[degrees]C and weighed. Given the destructive nature of the drilling process, bulk density samples could not be obtained. The soil EC of each sample was determined using a 1:5 soil/water suspension (Rayment and Higginson 1992). Gravimetric water contents were determined by oven drying the sample at 105[degrees]C and calculating the loss in weight.

EM survey and raw apparent electrical conductivity (E[C.sub.a]) data analysis

The day before the drilling commenced, an electromagnetic induction survey was performed along the transect. The GeonicsTM EM-38 and EM-34 electromagnetic induction devices were used in order to obtain 6 different depth-weighted E[C.sub.a] readings. Measurements were taken using both vertical and horizontal dipole configurations for both instruments. The EM-38, which has a fixed intercoil spacing of 1 m and operates at a frequency of 14600 Hz, was used to provide measurements to approximately 0.75 m in the horizontal dipole configuration and 1.5 m depth in the vertical dipole configuration. Measurements were taken at 10-m spacings along the transect (i.e. at every peg), resulting in a total of 75 measurements.

The EM-34 is designed for coil separations of 10m, 20m, and 40 m with corresponding frequency settings of 6400, 1600, and 400 Hz, respectively. The suggested depth of exploration is 0.75 x coil separation for the horizontal dipole configuration and 1.5 x coil separation for the vertical dipole configuration. In this study we only used the 10- and 20-m coil separation, since the suggested depth of exploration for the 40-m spacing would be too high for this study.

EM-34 readings at the 10- and 20-m separation were taken at every peg (i.e. at 10-m intervals) in both horizontal and vertical mode, with the 2 coils centred around the peg. This resulted in a total of 152 readings taken using the EM-34, with 4 different measurements taken at every peg.

Characterising the stratigraphy of the palceochannel profile

We assumed that we needed hydrologically important information at depths of 1 m down to a depth of approximately 10 m. Three different approaches were used to infer the EC values for the profile. The first inversion method is based on McNeill (1980) using an even-determined inversion of Eqn 1 (Menke 1989). The K matrix consists of n linear combinations of response functions, where n is the number of observations:


R[(z).sub.i] is either the vertical or horizontal cumulative response function of the instrument, depending on the coil configuration to obtain [d.sub.i]. Note that in this case the number of depth layers needs to equal the number of observations (n). The different EM measurements provided an E[C.sub.a] reading for the bulk soil approximately to the different exploration depths. For the layered earth model, we assumed that each layer is considered to have a thickness that equals the difference in exploration depth for 2 consecutive EM readings (McNeill 1980). The top 2 layers were determined from the EM-38 readings and each has a thickness of 0.75 m. The other layers varied in thickness from 6 to 15 m. The final bottom layer is assumed to be infinitely thick.

The inversion problem was solved in R using the optimisation procedure 'optim' (R Development Core Team 2004). We applied an upper bound of 5 dS/m and a lower bound of 0.10 dS/m to avoid negative solutions. To project the derived EC values to the l-m layers, a further cubic smoothing spline was fitted to the EC depth curves, using the 'smooth.spline' function with the smoothing parameter (spar) set at 0.5.

The second inversion method follows the approach taken by Cook and Walker (1992) to predict ECa for the 1-m layers. The minimisation problem was again solved in R using 'optim' by minimising the response above and below the depth of each layer and the responses applied to the measured E[C.sub.a] values. This method does not invert the EM signal and the results are therefore still apparent conductivities rather than bulk soil conductivities. However for ease of notation and for better understanding we have dropped the 'a' subscript from the resulting data.

The third inversion method follows the Tikhonov regularisation method described by Borchers et al. (1997), which involves finding the solution to the minimisation problem:

(3) min [[parallel]K[sigma] - d[parallel].sup.2] + [[lambda].sup.2] [[parallel]L[sigma][parallel].sup.2]

using non-negativity constraints. Here K, d, and [sigma] are the same as in Eqn 1, L can be chosen to be the discrete version of the second-derivative operator and the factor X controls the trade-off between minimising the norm of the residual and minimising [parallel]L[sigma][parallel] (Borchers et al. 1997), which is in fact a smoothing operation on the data. The MATLAB code developed by Borchers et al. (1997) was transcribed to R for this study. However, we chose to automatically select the optimal value of [lambda] from the derived L-curve using the general linear cross-validation method (Hansen 1992).

The analysis was limited to 10 m depth, since increased gravel and a possible aquifer exists below that depth, causing anomalous EM readings and very low correlations with soil properties.

We split the measured soil property data into a validation and a calibration set using Latin Hypercube sampling (McKay et al. 1979). This method of sampling for the dataset was chosen because the distribution of soil properties is very skewed, and sampling a small subset would give different results depending on whether the set includes the outer values of the distribution. By using Latin Hypercube sampling, we can ensure that the subset has the same distribution as the original dataset. The calibration set was used for all the subsequent analysis and this set consisted of approximately 65% of the original soil data.

The clay, silt, sand, and fine sand data from the calibration dataset were interpolated onto a grid 10m by 1 m using ordinary kriging in the package 'gstat' in R (Pebesma 2004). This spatial interpolation represents the soil property and stratigraphy information gained from sampling soil without performing additional EM measurements. This will be referred to as the interpolated soil data.

One way to improve the spatial interpolation of soil data would be to use the EM data as a trend surface and to interpolate the soil property data using universal kriging with trend. All 3 derived EC datasets were used to interpolate the calibration soil dataset using universal kriging with trend in 'gstat' (Pebesma 2004). The z-data were transformed to account for the strong anisotropy in the dataset. The dataset is much denser in the z-direction (sampled every 1 m) than in the x-direction (sampled every 10m). These datasets will be referred to as the 'trend kriged' datasets.

A second interpolation method involved regression kriging (Odeh et al. 1995). The derived EC data were correlated with the particle size fractions at the sampled locations to develop predictive equations for these soil properties using linear regression. The derived EC data were log-transformed before regression, as it should be expected that the prediction would be bounded by:

(4) EC [not equal to] 0 for Clay = 0 EC < [infinity] for Clay = Maximum

Other transformations could be the reciprocal and the square root of the EC values, since these would also conform to the boundaries in Eqn 4.

The EC data were converted to soil properties based on the linear regressions, after which the residuals were interpolated over the grid using ordinary kriging, and added to the predicted soil property values (Odeh et al. 1995). The z-data were again transformed to account for the strong anisotropy in the dataset. Only the regressions with the highest correlation coefficient ([r.sup.2]) and the lowest root mean square error (RMSE) were used to calculate particle size fractions from the EC data and remaining particle size fractions were calculated by difference. This dataset will be referred to as the 'regression kriged' datasets.

Because the hydrological properties of the palaeochannel system were the focus of this study, we calculated hydraulic properties from the soil properties using the 'Australia wide' pedotransfer functions in the program NEUROTHETA (Minasny and McBratney 2002). While NEUROTHETA is capable of predicting both the soil water characteristic curve and the K(h) relationship, we will limit our discussion to the saturated hydraulic conductivity ([K.sub.S]).

Results and discussion

Sampled soil properties

For most sampled sites, clay content increased below 2 m and then slowly decreased between 8 and 12m depth (Fig. 2). An exception was at 120m along the transect, where the palaeochannel is located (Peg 13 in Fig. 1). Silt fractions were fairly uniform along the transect. At all sampled sites, the sand distribution (Fig. 2) mirrored the clay distribution. The increase in sand at a depth of 2m at Peg 13 (120m along the transect) marks the top of the buried palaeochannel, while the decrease in the sand fraction and increase in the clay fraction at 7 m mark the bottom of the palaeochannel.


The palaeochannel in this field therefore appears to be a discrete entity consisting of 4 m of sand from a depth of 2 to 6m from the surface. The sand content within the palaeochannel bed is around 80%, with fine sand dominating from 2 to 4 m and coarse sand dominating in the bottom 2 m. This is most likely to be a depositional feature, where the coarser material settled at the bottom of the active channel and progressively finer sand was deposited as stream flow waned.

A lighter textured clay layer (containing 40-50% clay and 30-50% sand) occurred in the top 2m along most of the transect (Fig. 2). This is most likely to be a result of mixing of the coarser palaeochannel sediments with the clay during field levelling and cultivation, but could also be a result of different depositional events. However, the mixing hypothesis may also explain the extent of the lighter coloured soil seen in the aerial photograph (Fig. la) well outside the palaeochannel bed. Peg 38 (at 370 m along the transect) is located beyond the palaeochannels sediments, in a grey Vertosol (Isbell 1996) (see Fig. 1b).


A gravel bed (Fig. 2) containing various amounts of gravel (9-62%) and very coarse sand (6-57%) underlies the entire transect below 10m depth. We think that this gravel may define the top of an upper alluvial aquifer, which is reported to occur between 15 and 25 m in this region (Water Conservation and Irrigation Commission 1966). This gravel bed is also visible at varying depths in borelogs in the surrounding area.

During the drilling it was noted that the palaeochannel was saturated between 3 and 6 m below the surface due to a perched water table that had formed in the sand above the underlying clay. The field had been irrigated prior to sowing a wheat crop, around 2 months prior to the sampling. A perched water table is regularly observed in the field at the end of each cotton season (T. Richards, Agronomist, pers. comm.). The underlying gravel, in contrast, was very dry compared with the rest of the profile.

Electrical conductivities were low (<3 dS/m) around the palaeochannel (at 120 m) and in the lower part of the profile. This is expected as salts are more easily leached out of the more permeable sandy and gravelly soil than out of the surrounding clays.

Inversion of the EM signal

The different EM signal inversion methods gave quite different results (Fig. 3). The Tikhonov regularisation method produced the highest estimated conductivities (ranging from 0.26 to 3.52 dS/m) followed by the McNeill inversion method (ranging from 0.06 to 3.15dS/m) and lowest estimated conductivities were predicted by the Cook and Walker method (ranging from 0.46 to 1.26 dS/m). The lowest results from the Cook and Walker method are somewhat expected since this method does not invert the EM signal. However, it can be seen from the 2-D plots of the data that all methods reproduce an area of lower conductivity at the location of the palaeochannel (120-150m, Fig. 3), as well as possibly a bit earlier (around 50-70 m in the McNeill method), and at the end of the transect (around 320 m). Outside the palaeochannel sediments, the Cook and Walker method is quite constant with depth with mainly lateral variations. The check pattern in the Cook and Walker method is a feature of the method. It creates similar horizontal weighting bands for each depth, which are subsequently scaled by the predominant signal at each horizontal location.

The McNeill inversion produces an increase in conductivity with depth, with a rather sharp boundary at around 9m depth. Above this depth the conductivity is generally quite low, and below this depth the predicted conductivity is much higher. In contrast, the Tikhonov inversion indicates a band of higher conductivity between about 5 and 10 m depth, with lower conductivies above and below these depths. The soil data also indicated an initial gradual increase in clay content to 4 m depth and a drop in clay content below 8-10 m depth, coinciding with the Tikhonov inversion results. There are also differences in the way the 3 methods represent the sandier sediments, particularly in the width and lateral variation. Such differences could have a major impact on the prediction of water flow through the channel.

Below 15 m all EC values dropped sharply (data not shown). We believe that this sudden change in conductance is due to the presence of the gravel layer. Since drilling only reached 15 m as a maximum and most of the gravel found was dry, it remains unclear whether this is due to (most likely) a (saline) watertable at lower depths or different electromagnetic properties of the sediments and gravel beyond that depth.

It is important to note here that the original data were spread unevenly in the z-direction. Because this study only employed the EM-34 and EM-38, a major gap in the data existed between approximately 1.5 and 7.5m, which is particularly an issue for the even-determined problem from the McNeill inversion. As the palaeochannel appears to be most distinct at 2-7 m, this lack of data should have had an effect on the further analysis. In hindsight, using an EM-31 (for which the exploration depth can range from 2 to 6 m) for this study would have been more useful. However, this instrument was not available to us at the time of this study.

Correlation of inverted EM values with soil properties

The overall EM signal is a combination of the different soil components (Friedman 2005). However, stepwise linear regressions of the different soil properties with the EM-derived EC data did not give evidence of any strong relationships ([r.sup.2] range 0.32-0.47), but all regressions gave significant relationships (Table 1). The highest correlation was, surprisingly, for the relationship of the inverted McNeill data and included a combination of clay, sand, soil wetness, EC (1:5), and fine gravel contents as predictors. Even though the correlations were not great, these relationships give at least some confidence in the prediction of soil conductivity using the inversion process. However, since we are interested in predicting single soil properties from the derived EC data, such multiple regressions were not useful in this study.

Correlations from the single regressions between soil properties and the EC data were quite low, ranging from almost zero for the prediction of some properties for the Cook and Walker method to 0.44 for the prediction for clay content for the McNeill method (Table 2). Such low correlations lead to quite significant root mean square errors (RMSE), which make the relationships not very useful for prediction of these soil properties at locations where only an EM signal is available.

Although the salt content in the soil solution (as reflected by EC water) is generally regarded as the strongest determinant of the EM signal (Friedman 2005), EC (1:5 in water) exhibited only low correlations with the inverted EM signal (Table 2).

All regressions are quite strongly influenced by the variation in the data related to the sampling location (at 120 m) within the palaeochannel (closed squares in Fig. 4). The samples in, below, and above the palaeochannel are characterised by low EC (Fig. 2), while the sandy-textured samples in the palaeochannel (at depth 2-6 m) were also quite wet, leading to different EM responses. Although the Tikhonov and Cook and Walker method data have a lower correlation with the soil property data than the McNeill data, the distribution of the data is in this case more even. The McNeill data regressions, in contrast, are more strongly influenced by the palaeochannel sediments.


Overall the prediction of the soil properties at the, presumed unsampled, validation sites, using any of the inversion and interpolation methods based on the EM signal, is not very strong (Table 3). This is indicated by lower correlation coefficients and higher RMSE values for the validation sites in relation to the calibration sites. But it appears to be more a function of the interpolation method (kriging) than of the inversion method, as the decrease in prediction accuracy occurs also at the interpolated (ordinary kriged) soil dataset.

Although the Tikhonov inversion method appeared to have better correlations, the McNeill inversion method tended to have lower RMSEs at the calibration locations for the regression kriging method (Table 3). This indicated that while the trend was better predicted by the Tikhonov inversion dataset, the locational prediction was better for the McNeill inversion dataset. For the trend kriging interpolation the Tikhonov inversion dataset had the highest [r.sup.2] values and the lowest RMSE values for predicting soil properties in both the validation and calibration datasets.

Overall, trend kriging appeared to be a better predictor of the soil property data than regression kriging as evidenced by the lower RMSE and higher [r.sup.2] values at the validation locations. This is probably due to the original low correlations between the derived EC data and the soil properties (Table 2), which forms the basis for the regression kriging method. The Cook and Walker dataset, which had the lowest correlation with soil properties (Table 2), also exhibits the lowest predictive power in the regression kriging method.

Trend kriging based on the Tikhonov data as a trend surface also outperforms ordinary kriging of soil property data. This is encouraging, as the soil sample dataset in this study was relatively dense. For larger areas, where the soil sampling density is lower, EC data derived from inversion of more densely measured EM data would be a valuable addition to infer stratigraphies.

Prediction of hydraulic properties

We attempted direct correlation of the hydraulic conductivity with the E[C.sub.a] data, but this indicated very low correlations. We therefore decided not to present this as an alternative. This lack of correlation was somewhat surprising as we had expected that the EM signal would be sensitive to the same soil properties that determine [K.sub.S]. This probably indicates the importance of soil structure and bulk density in prediction of [K.sub.S] and the EM measurement could be relatively insensitive to these properties. However, it could also be related to the uncertainty of the prediction of hydraulic properties using pedotransfer functions.

Visually, there were only small differences between the interpolation methods and the different EM inversion methods in terms of prediction of the stratigraphy of log saturated hydraulic conductivity for the transect (Fig. 5). From this figure, the similarities between the interpolation methods for the same inversion method are immediately apparent. For example the Tihonov trend kriged results are more similar to the Tikhonov regression kriged results than to any other method. It is difficult to judge the goodness of the final predicted hydraulic properties, due to the fact that the predicted results include uncertainties due to sampling density, inversion and interpolation techniques, and pedotransfer functions. We have chosen to use [K.sub.S] values derived from all the soil data (at both the validation and calibration sites) as a benchmark, because there are no measured [K.sub.S] values for the transect.


Overall correlations between predicted hydraulic properties using the regression kriging and trend kriging methods at all the sampled locations range from 0.69 to 0.87 (Fig. 6) as is also reflected in the spatial patterns (Fig. 5). There was high correlation (>0.99) between kriging methods based on the same inversion method. This appears to indicate that any differences in the spatial patterns are mainly due to inversion methods and less due to the interpolation techniques, and that the final resulting hydraulic conductivity field was very similar for all methods.


Comparison between predicted [K.sub.S] values using the different EM inversion and interpolation methods and the [K.sub.S] values derived from all sampled soils data indicated that the difference between the Tikhonov inversion and Cook and Walker method is only small in terms of [r.sup.2] and deviation of the 1:1 line (Fig. 6); however, the McNeill inversion method appears less accurate after all interpolation. In contrast to this study, Borchers et al. (1997) and Hendrickx et al. (2002) found large differences in predicted EC profiles between the Cook and Walker method and Tikhonov inversion.

We can focus on how these different methods predict hydraulic conductivities needed for the simulation of recharge and soil water processes at a specific location along the transect (Fig. 7). Again, it seems the prediction of hydraulic properties is relatively even across the methods. Most of the methods can predict smooth transitions in soil properties and hydraulic conductivities reasonably well (i.e. at the 200-m sample site). Particularly the interpolated McNeill inversion method has some difficulties with sharp transitions in the soil properties or hydraulic conductivities (i.e. at 4-6m depth at the 50-m and 120-m sample sites). Other researchers predicting soil conductivities have found similar results (Cook and Walker 1992; Borchers et al. 1997), and this reinforces the strength of the other methods.


The prediction of hydraulic properties in NEUROTHETA includes a considerable uncertainty, which is indicated by the 95% error bars (Fig. 7), because the prediction is based on a database of soils and the program predicts the average value. This means that while the average hydraulic conductivities at each depth appear different for the different methods, the values are not significantly different. While this might also translate into no significant difference in water flow, this is not certain. Using stochastic modelling approaches, which include input and model uncertainty (Krzysztofowicz 1999), these differences could be assessed. From this study it seems that the uncertainty created by the prediction of hydraulic properties far outweighs the uncertainty in the inversion and interpolation methods.

In addition to the uncertainty in the prediction of hydraulic properties in NEUROTHETA, there is also uncertainty in the prediction of the soil properties from EC data and subsequent interpolation and the initial inversion of the EM data. The regression kriging method incorporates some of the uncertainty in the prediction of soil properties from EM data through inclusion of the regression errors in the interpolation method. The trend kriging method does not include this, as it only uses the EC data as a co-predictor. Uncertainties could be calculated, however, by using small perturbations in the EC field. Uncertainty due to the inversion method is unknown at this stage, but could be calculated in the Tikhonov inversion process by varying the regularisation parameter [lambda], in a narrow range around its optimal value and calculating the different realisations of the inversion. This would again allow calculation of 95% confidence intervals.

In this study we have limited our analysis to 3 inversion techniques and 2 interpolation methods. Other versions of the L operator in the Tikhonov inversion can be chosen (i.e. McBratney et al. 2000) and higher dimensional non-linear inversion techniques (i.e. Auken and Christiansen 2004; Monteiro Santos 2004) could be applied. The simplicity of the linear inversion, however, makes the method easy to handle and available to a wide range of practitioners. The R-scripts used in this study are available on request from the corresponding author.

The success of the trend kriging method in predicting hydraulic conductivity also implies another possible interpolation method. The EC data could be interpreted as a scaling factor field for the [K.sub.S] data and variations in [K.sub.S] could be calculated based on a reference [K.sub.S] function for the profile or for identified distinct sediments. This would circumvent the problem of predicting hydraulic properties from soil property data using pedotransfer functions.

An extension on this study is the implementation of the derived [K.sub.S] fields into a hydrological model and comparison of the final hydrological response of the different inversion and interpolation methods. We are currently attempting this for another dataset.


This research has shown that, although the overall patterns of hydraulic conductivity (and thus the stratigraphy) can be inferred from the combination of EM inversion, soil sampling, and pedotransfer functions, the absolute magnitude of the hydraulic conductivity cannot easily be predicted. This is partly due to uncertainties in the inversion techniques but mainly due to uncertainties in the conversion of soil property data into hydraulic properties. In this research it was shown that the combination of EM data with sparse sampling would be more beneficial at larger scales to characterise hydrologically important stratigraphies. The different inversion methods, McNeill (McNeill 1980) and Tikhonov regularisation (Borchers et al. 1997) and the method developed by Cook and Walker (1992), caused only small differences in stratigraphies. After interpolation and conversion, the Cook and Walker and Tikhonov methods appeared to represent the existing [K.sub.S] field more accurately than the McNeill inversion method. Trend kriging as an interpolation method for soil property data, using the derived EC data as a trend surface, seemed to be superior to regression kriging based on correlations between derived EC data and sampled soil property data. Differences between interpolation and inversion methods were mostly negated in the final prediction of hydraulic properties due to uncertainties in the prediction using pedotransfer functions.


This research was partly funded by the Australian Cotton CRC through an honours research scholarship for YLA. The authors would like to thank Auscott Ltd for access to the field site and ongoing research support, Dr B Borchers for making his MATLAB code available on the worldwide web, and Dr B Minasny for assisting in converting some of the MATLAB code to R.


Acworth RI (1999) Investigation of dryland salinity using the electrical image method. Australian Journal of Soil Research 37, 623-636.

Auken E, Christiansen C (2004) Layered and laterally constrained 2D inversion of resistivity data. Geophysics 69, 752-761. doi: 10.1190/1.1759461

Borchers B, Uram T, Hendrickx JHM (1997) Tikhonov regularisation of electrical conductivity depth profiles in field soils. Soil Science Society of America Journal 61, 1004-1009.

Butler BE (1950) A theory of prior streams as a causal factor of soil occurrence in the riverine plain of south-eastern Australia. Australian Journal of Experimental Agriculture 1, 231-252. doi: 10.1071/AR9500231

Cook PG, Walker GR (1992) Depth profiles of electrical conductivity from linear combinations of electromagnetic induction measurements. Soil Science Society of America Journal 56, 1015-1022.

Corwin DL, Rhoades JD (1984) Measurement of inverted electrical conductivity profiles using electromagnetic induction. Soil Science Society of America Journal 48, 288-291.

Friedman SP (2005) Soil properties influencing apparent electrical conductivity: a review. Computers and Electronics in Agriculture 46, 45-70. doi: 10.1016/j.compag.2004.11.001

Gee GW, Bauden JW (1986) Particle size analysis. In 'Methods of soil analysis. Part 1'. (Ed. A Klute) (American Society of Agronomy, Soil Science Society of America: Madison, WI)

Hansen PC (1992) Analysis of discrete ill-posed problems by means of the L-curve. SIAMReview 34, 561-580. doi: 10.1137/1034115

Hendrickx JMH, Borchers B, Corwin DL, Lesch SM, Hilgendorf AC, Schlue J (2002) Inversion of soil conductivity profiles from electromagnetic induction measurements: Theory and experimental verification. Soil Science Society of America Journal 66, 673-685.

Isbell RF (1996) 'The Australian Soil Classification.' (CSIRO Publishing: Collingwood, Vic.)

Krzysztofowicz R (1999) Bayesian theory of probabilistic forecasting via a deterministic hydrologic model. Water Resources Research 35, 2739-2750. doi: 10.1029/1999WR900099

McBratney AB, Bishop TFA, Teliatnikov IS (2000) Two profile reconstruction techniques. Geoderma 97, 209-221. doi: 10.1016/ S0016-7061(00)00039-2

McKay MD, Beckman RJ, Conover WJ (1979) A comparison of three methods for selecting values of input variables in the analysis of output from a computer code. Technometrics 21,239-245.

McNeill JD (1980) Electromagnetic terrain conductivity measurement at low induction numbers. GEONICS INC. Technical Note TN-6.

Menke W (1989) 'Geophysical data analysis: Discrete inverse theory.' (Academic Press: San Diego, CA)

Minasny B, McBratney AB (2002) The neuro-m method for fitting neural network parametric pedotransfer functions. Soil Science Society of America Journal 66, 352-361.

Monteiro Santos FA (2004) 1-D laterally constrained inversion of EM34 profiling data. Journal of Applied Geophysics 56, 123-134. doi: 10.1016/j.jappgeo.2004.04.005

Odeh I, McBratney AB, Chittleborough DJ (1995) Further results on prediction of soil properties from terrain attributes: heterotropic cokriging and regression-kriging. Geoderma 67, 215-226. doi: 10.1016/0016-7061(95)00007-B

Page K, Nanson G, Price D (1996) Chronology of Murrumbidgee River Palaeochannels on the Riverine Plain, southeastern Australia. Journal of Quaternary Science 11, 311-326. doi: 10.1002/(SICI) 1099-1417(199607/08) 11:4 < 311 ::AID-JQS256 > 3.3.CO;2-T

Pebesma EJ (2004) Multivariable geostatistics in S: the gstat package. Computers and Geosciences 30, 683-691. doi: 10.1016/ j.cageo.2004.03.012

R Development Core Team (2004) 'R: A language and environment for statistical computing.' (R Foundation for Statistical Computing: Vienna) (URL

Rayment GE, Higginson FR (1992) 'Australian laboratory handbook of water chemical methods.' (Inkata Press: Melbourne, Vic.)

Rogers MP, Christen EW, Khan S (2002) Aquifer identification and characterisation for salinity control by shallow groundwater pumping: A case study of the Coleambally Irrigation Area. CSIRO Land and Water, Griffith, Technical Report 16/02.

Stannard ME, Kelly ID (1968) 'The irrigation potential of the Iower Gwydir Valley.' (Water Resources Commission of New South Wales: Sydney)

Triantafilis J, Huckel AI, Odeh IOA (2003) Field-scale assessment of deep drainage risk. Irrigation Science 21, 183-192.

Triantafilis J, Lesch SM (2005) Mapping clay content variation using electromagnetic induction techniques. Computers and Electronics in Agriculture 46, 203-237. doi: 10.1016/j.compag.2004.11.006

Young RW, Young ARM, Price DM, Wray RAL (2002) Geomorphology of the Namoi alluvial plain, northwestern New South Wales. Australian Journal of Earth Sciences 49, 509-523. doi: 10.1046/j. 1440-0952.2002.00934.x

Water Conservation and Irrigation Commission (1966) Water resources of the Gwydir Valley, Survey of Thirty NSW River Valleys. Water Conservation and Irrigation Commission Report No. 5.

Williams BG, Hoey D (1987) The use of electromagnetic induction to detect the spatial variability of the salt and clay contents of soils. Australian Journal of Soil Research 25, 21-27. doi: 10.1071/SR9870021

Manuscript received 21 March 2005, accepted 11 October 2005

R. W. Vervoort (A,B) and Y. L. Annen (A)

(A) Faculty of Agriculture, Food and Natural Resources, The University of Sydney, NSW 2006, Australia, and Australian Cotton Cooperative Research Centre, Narrabri, NSW 2390, Australia.

(B) Corresponding author. Email:
Table 1. Results of the stepwise multiple linear regressions between
measured soil properties and inverted soil conductivity profiles

The results indicate overall low correlations but significant
relationships at P = 0.05

Inversion method    [r.sup.2]     RMSE              Predictors

Tikhonov              0.34        0.55     Soil wetness and fine sand
McNeill               0.47        0.36     Sand, gravel and fine gravel
                                             contents, soil wetness and
                                             EC 1 : 5 in water
Cook and Walker       0.37        0.12     Clay, sand silt, and fine
                                             gravel contents

Table 2. Correlations between soil properties and log-transformed
derived EC values using the different inversion methods, indicating
overall low correlation between the derived EC values and soil
properties, with the Cook and Walker (1992) method being the least
related to sampled soil properties as indicated by lowest [r.sup.2] and
highest root mean square error (RMSE)

                              McNeill                  Tikhonov

Soil property          [r.sup.2]    RMSE (%)    [r.sup.2]    RMSE (%)

Clay content            0.44 *       11.17       0.28 *       12.72
Silt content            8.0E-3        5.03       0.03 *        4.98
Sand content            0.43 *       11.87       0.20 *       14.04
Fine sand content       0.32 *        9.77       0.37 *        9.42
Coarse sand content     0.14 *        8.65       0.08 *        9.34
EC (1 : 5) in water     0.25 *      293.7        0.22 *      300.1

                          Cook and Walker

Soil property          [r.sup.2]    RMSE (%)

Clay content            0.13 *       13.92
Silt content            1.6E-5        5.05
Sand content            0.13 *       14.66
Fine sand content       0.11 *       11.19
Coarse sand content     0.03 *        9.20
EC (1 : 5) in water     4.3E-5      339.8

* Relationship significant (P < 0.05).

Table 3. Comparison between observed and predicted soil properties for
the calibration and validation datasets

The [r.sup.2] values reflect the level of explained variation of the
relationship, while the root mean square error (RMSE) gives and
indication of how well the points are predicted on average. As can be
seen from the results, specific points in the domain are not very well
predicted for any of the inversion and extrapolation methods, except
the kriging of soil properties. Abbreviations used are C, calibration
dataset; V, validation dataset; RK, regression kriging; TK, trend
kriging; Tikh, Tikhonov inversion; McNeill, McNeill inversion; C&W,
Cook and Walker (1992) method; Soil, interpolated soil dataset

                  RK         Tikh         RK        McNeill
               [r.sup.2]     RMSE      [r.sup.2]     RMSE
                              (%)                     (%)

C Clay          0.64         9.02       0.59          7.90
C Sand          0.59         9.88       0.53          8.79
C Fine sand     0.70         6.07       0.70          5.48
V Clay          0.08 *      12.86       0.18 *       11.78
V Sand          0.05 *      14.19       0.12 *       13.34
V Fine sand     0.17 *       7.57       0.28          8.18

                  RK         C&W        TK        Tikh
               [r.sup.2]    RMSE     [r.sup.2]    RMSE
                             (%)                   (%)

C Clay          0.53        10.69      0.89       4.22
C Sand          0.48        12.02      0.90       4.10
C Fine sand     0.63         6.86      0.85       3.83
V Clay          0.15 *      16.45      0.44       8.32
V Sand          0.10 *      17.66      0.50       7.72
V Fine sand     0.31         9.15      0.67       5.64

                  TK        McNeill       TK         C&W
               [r.sup.2]     RMSE      [r.sup.2]     RMSE
                              (%)                    (%)

C Clay           0.82        4.77        0.86        5.11
C Sand           0.83        4.69        0.86        5.41
C Fine sand      0.85        3.64        0.87        3.69
V Clay           0.38        9.58        0.33       11.54
V Sand           0.37        9.88        0.37       11.29
V Fine sand      0.58        7.13        0.64        6.44

               [r.sup.2]    Soil

C Clay           0.92       3.43
C Sand           0.92       3.65
C Fine sand      0.89       3.14
V Clay           0.31       9.11
V Sand           0.37       8.90
V Fine sand      0.65       5.53

* Non-significant relationship (P > 0.05).
Gale Copyright: Copyright 2006 Gale, Cengage Learning. All rights reserved.