Earth Surface Processes and Landforms Airborne survey32, signals Earth Surf.radiometric Process. Landforms 1503–1515 (2007) Published online 30 January 2007 in Wiley InterScience (www.interscience.wiley.com) DOI: 10.1002/esp.1468

1503

Understanding airborne radiometric survey signals across part of eastern England B. G. Rawlins,1* R. M. Lark2 and R. Webster2 1 2

British Geological Survey, Keyworth, Nottingham, UK Rothamsted Research, Harpenden, Hertfordshire, UK

*Correspondence to: B. G. Rawlins, British Geological Survey, Keyworth, Nottingham NG12 5GG, UK. E-mail: [email protected]

Received 3 March 2006; Revised 23 October 2006; Accepted 2 November 2006

Abstract A low-level airborne radiometric survey provides data on the concentrations of gammaemitting elements including potassium (K), thorium (Th) and uranium (U) in the upper half metre of the soil. Where weathering has not penetrated much beyond this depth, as in the young soils that cover much of England and Wales, the signal is likely to be related to the soil’s clay content and its parent material. In these situations radiometric survey could be valuable for mapping soil digitally. We wished to understand how the radiometric signal relates to parent material and soil geochemistry, and to identify the spatial structure, if not the sources, of any unexplained variation. We analysed the joint spatial variation of the airborne gamma signal and highresolution soil geochemical survey data across part of eastern England by modelling their coregionalization. We also used REML to assess the joint effects of soil geochemistry and parent material on the radiometric signals of K and Th. The overall correlations of radiometric estimates with soil survey data for K and Th were large, as were the structural correlations for components of variation spatially dependent up to 49 and 16 km for K and Th respectively. This suggests that the radiometric signals for these two elements provide effective estimates of the amounts in the soil and their patterns of distribution. Although class of parent material accounted for a third of of the variance in the radiometric K signal, much of the variation within the classes is explained by geochemistry, suggesting that subtler changes can be detected. A larger proportion of the Th signal was accounted for by parent material. This supports our expectation that radiometric signals for K and Th provide information on parent material in the young landscapes of England and Wales. We are therefore confident that airborne radiometric surveys would be useful for making thematic maps of soil, particularly the soil’s texture and closely related properties across England and Wales. Copyright © 2007 Natural Environment Research Council. Published in 2007 by John Wiley & Sons, Ltd. Keywords: geochemistry; soil mapping; geostatistics; coregionalization;

REML

Introduction There is a demand for soil information by policy makers, regulators and environmental managers who must make decisions about land use that either depend on the properties of the soil or have the potential to modify those properties for better or worse. Traditionally, this requirement for soil information has been met by soil surveys. In many countries, however, there has been little or no soil survey, particularly at large scales. Even in England and Wales, for example, there is national coverage only at 1:250 000, for which the legend is of complex map units (soil associations). These map units group together more or less disparate types of soil. For many purposes we require more precise information, which requires, for conventional soil survey, a larger scale. Only about 25% of the country has been mapped at larger scales (1:63 360 and 1:25 000), and further field survey is unlikely to be funded in the near future. This situation is common, and in most countries of the world coverage at such scales is even rarer. Soil scientists are trying to meet the demand for soil information by using various kinds of remote sensing and related sources such as digital terrain models to map soil properties or profile classes from sparse observations in the Published in 2007 by John Wiley & Sons, Ltd.

Earth Surf. Process. Landforms 32, 1503–1515 (2007) DOI: 10.1002/esp

1504

B. G. Rawlins, R. M. Lark and R. Webster

field. This approach has come to be known as digital soil mapping (Lagacherie et al., 2006). The predictor variables from these sources are chosen not blindly, but because they reflect the key factors of soil formation. Jenny (1941) famously identified these factors as climate (Cl), organisms (including vegetation) (O), relief (R), parent material (P) and time (T). Data from remote sensors that provide information on net primary production are clearly useful (factor O), as are digital elevation models (R) and geological maps (P). Australian soil scientists have successfully used airborne radiometric (gamma-ray) survey data for digital soil mapping (Cook et al., 1996; McKenzie and Ryan, 1999; Wilford, 2007). Such surveys are used to estimate the concentrations of potassium (K), thorium (Th) and uranium (U) in the the upper 35 cm of the soil. These properties are not of direct interest in themselves. However, in the old landscapes of Australia the total K content of the topsoil is a good surrogate for the age of the parent material (factor T), since the longer it has been subject to weathering the smaller will be its K content. Thus McKenzie and Gallant (2006) used the radiometric K signal, along with other environmental covariates, to map contrasting soils in a catchment. Wilford and Minty (2006) note that while the total K signal declines with increasing time, the U and Th signals can increase relative to K as weathering progresses, and so combinations of these signals might be good surrogates for the age of the weathered material. Variations in the radiometric signal might also correspond to inherited differences in the mineralogical composition of the parent material. For example, the K content of acid igneous rocks is typically large, and Th is associated with granite, pegmatite and gneiss (Wilford and Minty, 2006). For this reason Wilford and Minty suggested that no soil is likely to have a unique radiometric signature everywhere, but nevertheless valuable information on the soil might be extracted from the radiometric variations within geomorphic units and regions with distinct lithology and geochemistry. We are investigating the sources of variation in the gamma radiometric signal from soil in England. In this paper we report results from our first studies over the east of the country, a much younger landscape than those where the Australians have succesfully used radiometry to predict soil properties. We want to know what factors underly variation of the gamma signal under these different conditions. Mayr et al. (2001) assessed the potential of gamma radiometry, in combination with other variables, to predict soil properties in small areas in England. They found that the radiometric data, and in particular the estimated contents of K and Th, predicted the topsoil texture well. Weathering in England is confined largely to the uppermost 1 m (in contrast to sub-tropical and tropical terrains), and the local variation in properties of the parent materials is therefore inherited in the soil. Furthermore, the predominant clay minerals (mica and mica–smectite; Loveland, 1984) in the soil have K integral to their structure. The contents of clay and K in the soil are therefore closely correlated, and so the gamma signal is a good surrogate for parent material, Jenny’s factor P. To our knowledge there are few regions where both high-resolution radiometric surveys and detailed geochemical soil data are available for comparison. One such region is part of eastern England, where geochemical data have been obtained by intense soil sampling (Rawlins et al., 2003). The data include measurements of the total K, Th and U in the soil, plus other elements. The measurement error has also been assessed. Peart et al. (2003) report the gamma radiometric survey of the region. We have analysed the joint spatial variation of the gamma signal and the soil geochemical data using multivariate geostatistical methods. This was to see whether the geochemistry is correlated with the remotely sensed signal, and whether this correlation depends on the spatial scale (i.e. whether the variations in geochemistry and radiometry at some spatial scales are more strongly correlated than at others). We then analysed the joint effect of soil geochemistry and parent material on the radiometric signal using the linear mixed model, in which these factors are fixed effects, and unexplained variation is modelled as a combination of independent identically distributed, two-dimensional random variation (part of which will be measurement error) and a spatially correlated random effect. This allows us to compare the effects of these factors on the radiometric signal, and to identify the spatial structure, if not the sources, of the unexplained variation.

Study Region and Surveys A radiometric survey was flown between May and September 1998 covering approximately 14 000 km2 of central England, bounded by British National Grid eastings 310 to 510 km and northings 315 to 385 km. In this study we focus on the north-easterly section of this region, bounded by the coordinates 450 km easting and 345 km northing (south-west corner) and 510 km easting and 385 km northing (north-east corner), an area of approximately 2400 km2 (Figure 1). We chose this area because

• •

comprehensive geochemical survey data (including K) were also available at a resolution of 1 soil sample per 2 km2 (Rawlins et al., 2003), we wished to constrain the number of types of parent material and

Published in 2007 by John Wiley & Sons, Ltd.

Earth Surf. Process. Landforms 32, 1503–1515 (2007) DOI: 10.1002/esp

1505

Published in 2007 by John Wiley & Sons, Ltd.

Earth Surf. Process. Landforms 32, 1503–1515 (2007) DOI: 10.1002/esp

Figure 1. Study region showing the major soil parent material types. The scale is approximately 1:300 000.

Airborne radiometric survey signals

1506

•

B. G. Rawlins, R. M. Lark and R. Webster

the land use and topography are more uniform here than they are further to the west. For example, across the region we studied (Figure 1) elevation varies from 3 m below mean sea level to 212 m above it, and a very large proportion (81%) of the land is under tillage; the rest is grassland (6%), suburbs (4%), deciduous woodland (3%) and small areas of other land use types (6%) (percentages taken from Fuller et al., 1994).

The soil’s parent materials comprise sedimentary formations ranging from Permian to Jurassic in age, dipping to the east, and a series of overlying Quaternary deposits including glacial till, riverine alluvium and river terrace deposits (Figure 1). The outcrops of the geological formations are aligned predominantly from north to south as a result of the inherited basement structure and active subsidence in the North Sea. The soil has developed largely since the Devensian glaciation (10 000 years BP), and therefore strongly reflects the composition of the parent material from which it formed. The soil types, in the classification of Avery (1990), as percentages of the geochemical sampling sites (n = 917, see below), include Brown Soils (41%), Surface-Water Gleys (20%), Ground-Water Gleys (17%), Pelosols (12%) and Peats (1%). Average annual rainfall across the region, assessed from data from 1941 to 1970, varies from a maximum of 800 mm over the western and eastern parts of the region to a minimum of around 560 mm over central and southern parts (Department of the Environment, 1992; MAFF, 2000). The geochemical survey of the soil was done in the summers of 1994, 1995 and 1996 in rural and peri-urban areas throughout the region. Sampling sites were chosen from alternate kilometre squares of the British National Grid by simple random selection within each square, subject to the avoidance of roads, tracks, railways, urban areas and other seriously disturbed ground. There were 917 in total. At each site soil was taken from between 35 and 50 cm from five holes augered at the corners and centre of a square with a side of lengty 20 m by a hand auger and combined to form a bulked sample. All samples of soil were dried and disaggregated. They were sieved to pass 150 µm, coned and quartered. From each a 50 g sub-sample was ground in an agate planetary ball mill. The total concentrations of major and trace elements were determined in each sample by wavelength dispersive XRFS (X-ray fluorescence spectrometry). These included K as K2O (%), Th and U (mg kg−1). These are the elements that concern us here. Sub-samples were taken from a subset of the original samples and analysed separately; the analytical errors expressed as percentages of the means were 0·2% for K, 6·4% for Th and 21% for U. The large error for U is due to a combination of abundances generally close to the limit of detection (0·6 mg kg−1) and a pronounced ‘nugget effect’; the latter term is used by analytical chemists (with a meaning that is related to the same term used by geostatisticians) to describe the large errors in concentrations of elements such as U following repeated analyses of particulate sub-samples in which the analyte occurs in a few discrete mineral phases (Stanley, 1998). For the radiometric survey, the separation between flight lines was approximately 400 m, with perpendicular tielines spaced at 1200 m intervals, except over an infill area of special interest, where flight-line and tie-line separations were reduced to 200 m and 600 m respectively. In all there were 109 919 radiometric data points in the region. The gamma-ray spectrometer was a 256-channel Picodas PGAM 1000 (model 6·11) with a cycle rate of 1 Hz (equating to about 70 m along the traverses). The gamma signals were converted into measurements of K (%), Th (mg kg−1) and U (mg kg−1) and total counts (not presented here), with 98% of the signal originating from the upper 35 cm of the soil (Minty, 1997). Based on a nominal survey altitude of 90 m, 80% of the signal originates from a circle with an approximated radius of 200 m beneath the detector (Duval et al., 1971). We assigned a code to each and every one of 25 main types of parent material. Then at the location of each radiometric data point we determined from the 1:50 000 maps of bedrock geology and overlying Quaternary deposits which code applied there. Figure 1 shows a simplified version of the result. We extracted the radiometric data nearest to each of the 917 locations of the soil samples. The mean distance between the sample sites and the corresponding radiometric measurement was 98 m, with an interquartile range of 92 m. These paired observations on soil properties and gamma emissions were then used in the analyses described in the next section.

Statistical Analyses and Interpretation Comparison of soil and radiometric survey data We explored the data by calculating summary statistics for each of the three elements (Table I). The mean concentrations of the elements in the soil geochemical survey are similar to those reported by Reimann and de Caritat (1998) for soils world-wide; results from the soil survey and the published values are, respectively, K, 2·04 and 1·14%; Th, 8·15 and 9·4 mg kg−1; U, 2·45 and 2·7 mg kg−1. All the distributions were positively skewed, but for none of the elements were the skewness coefficients consistently large enough to warrant transformations. We therefore did all subsequent analyses on the data as they were. Published in 2007 by John Wiley & Sons, Ltd.

Earth Surf. Process. Landforms 32, 1503–1515 (2007) DOI: 10.1002/esp

Airborne radiometric survey signals

1507

Table I. Summary statistics for the soil and radiometric survey data. Units are % for K and mg kg−1 for Th and U. By SS we denote the soil survey data in the G-BASE database, and by Rad the radiometric data K

Mean Median Variance Standard deviation Skewness

Th

U

SS

Rad

SS

Rad

SS

Rad

2·04 1·91 0·45 0·67 0·54

0·94 0·86 0·27 0·52 0·86

8·15 8·10 4·26 2·06 0·27

4·62 4·56 3·72 1·93 0·12

2·45 2·40 0·54 0·74 1·81

1·24 1·22 0·29 0·53 0·42

Table II. Correlation matrix for K, Th and U for the soil survey (SS) data and the corresponding radiometric data (Rad)

SS–K SS–Th SS–U Rad–K Rad–Th Rad–U

SS–K

SS–Th

SS –U

Rad–K

Rad –Th

Rad–U

1 0·09 0·10 0·82 0·36 0·10

1 0·44 0·14 0·46 0·26

1 0·07 0·14 0·14

1 0·62 0·27

1 0·54

1

The mean values from the soil survey were significantly greater than the radiometric data for each of the three elements. This difference arises because the soil survey samples were sieved, thereby removing the coarse (larger than 150 µm) fraction, which is dominated by quartz, and resulting in larger (on average) reported concentrations than the concentrations estimated for the whole soil by the radiometric data. Table II shows the correlations between the radiometry data and the geochemical measurements for the 917 soil survey data locations. There is a strong correlation between the K concentrations by the two methods (r = 0·82). By contrast there is a moderate correlation for Th (r = 0·46) and only weak correlation for U (r = 0·14), probably because of the greater analytical errors for these two elements for both XRFS (see above) and estimation based on gamma emissions (Minty, 1997). We chose to pursue the spatial relationships of K and Th determined by the two methods by investigating their coregionalization (see below). In the case of Th and U, their moderate positive correlations for analyses between methods (r = 0·44 and 0·54) is likely to result from a geochemical association. They are both tetravalent and hence substitute for Zr4+ (zirconium) in the mineral zircon, a common detrital constituent of weathered deposits, which accounts for a significant fraction of U and Th in whole rock analyses (Hoskin and Schaltegger, 2003). Another notable feature was the much stronger correlation between K and Th estimated by gamma emission (r = 0·62) than that in the soil survey data (r = 0·09). We chose to examine this relationship for the paired observations in more detail. First we subtracted the difference between the mean values for each data set (soil survey and radiometric) from each soil survey measurement to account for their larger concentrations caused by the sieving described above. We then plotted the modified values for the paired observations for K against Th on the same axes. The radiometric data with both large K (>1·5%) and Th (>4 mg kg−1) concentrations occur predominantly over the Triassic (Mercia) mudstone, resulting in a greater positive correlation between these elements than for the soil survey data. Figure 2 shows that this discrepancy is explained mainly by the larger Th concentrations reported in the radiometric data than those in the soil survey data. Given that we would expect the analytical errors associated with soil survey to be significantly less than those from estimation by gamma-ray emission, it is likely that the latter is overestimating Th concentrations at these larger concentrations. We should therefore be cautious in interpreting the strong correlation reported between K and Th in the radiometric data. Published in 2007 by John Wiley & Sons, Ltd.

Earth Surf. Process. Landforms 32, 1503–1515 (2007) DOI: 10.1002/esp

1508

B. G. Rawlins, R. M. Lark and R. Webster

Figure 2. Concentrations of K and Th in paired observations (n = 917) after subtractive correction – subtracting the difference between the means of K and Th (radiometric and soil survey) from each respective datum from the soil survey.

Coregionalization of K and Th from radiometry and soil analysis We denote by zu(x) a measurement of the concentration of an element in a soil sample at location x and by zv(x) a measurement of the concentration of the same element derived from the nearest radiometric measurement. The joint spatial variation of these two variables is described by their cross-variogram, which we estimate by

8 uv (h) =

1 2m(h)

m( h )

∑ {zu (x j ) − zv (x j j =1

}{

}

+ h) zv (x j ) − zv (x j + h)

(1)

where xj + h is the location of an observation separated from location xj by the lag vector h, and m(h) is the number of paired comparisons in the summation separated by h. The spatial variation of just one of these variables is described by its sample autovariogram, given by Equation (1) when u = v. We estimated the auto-variograms for the soil measurements of K and Th and the corresponding radiometric measurements, and the cross-variogram for the soil and radiometric measurement of each element. Although the soil sampling and radiometric measurements were not collocated, in the vast majority of cases the samples of soil were taken within the support of the latter; the mean separation of the locations is 92 m and 80% of the radiometric signal originates within a circle of radius 200 m. The effect of the separations between sampling locations on the crossvariograms (and subsequent coregionalization) will therefore be trivially small. Exploratory analysis suggested that the variation was isotropic (i.e., the variograms depend on the lag distance, but not on the direction). In these circumstances the lag becomes a scalar in distance only, so that h = |h|. To each pair of auto-variograms and the corresponding cross-variogram we fitted the linear model of coregionalization: K

γ uv (h) =

∑ buvk gk (h)

(2)

k =1

in which the gk, k = 1, 2, . . . , K, are basic variogram functions with unit a priori variance and specific distance parameters. The coefficents b kuv, in which the k are simply indices, not exponents, are (co)variances, i.e. sill variances of second-order stationary processes. In the event, we needed only two basic functions (K = 2) for a satisfactory fit. Published in 2007 by John Wiley & Sons, Ltd.

Earth Surf. Process. Landforms 32, 1503–1515 (2007) DOI: 10.1002/esp

Airborne radiometric survey signals

1509

The first, g1(h), is a nugget variance representing measurement error plus variation that is spatially correlated only over distances much shorter than the shortest distance for which we can compute 8 uv(h). The second term was an exponential function: ⎧ ⎛ h⎞⎫ γ (h) = c⎨1 − exp − ⎬ ⎝ a⎠⎭ ⎩

(3)

with sill c, and a distance parameter, a, that defines the spatial extent of the model. The function approaches the sill asymptotically and so does not have a finite range. For practical purposes it is assigned an effective range where γ is 95% of the sill variance, approximately 3a. These two components of the model therefore correspond to different components of variation of the two variables at different spatial scales. We also computed the hull of perfect correlation, the bounds on the linear model: K

k k hull [γ uv (h)] = ±∑ buu bvv gk (h).

(4)

k =1

The stronger the cross correlation the closer the linear model lies to one or other bound, and so plotting the fitted model inside the hull is a valuable diagnostic. We computed the autovariograms and cross-variogram for the radiometically estimated K and K measured in the soil to lag 30 km at intervals of 1 km with Equation (1). We then fitted a linear model of coregionalization, Equation (2), to them using the simulated annealing method described by Lark and Papritz (2003). As above, we incorporated two components in our model: a nugget g1(h) with h = 0 and an exponential function g2(h). k that minimize the sum of the squares of the residuals from the model, The fitting finds estimates of a and of the b uv weighted by the number of pair-comparisons at each lag, subject to the constraint that k k k buvk = bvu = ≤ buu bvv

for u, v and k = 1, 2.

We repeated the procedure for Th but to a maximum lag of only 20 km. We also calculated the structural correlation coefficients to aid our understanding of the changes in correlation with changes in spatial scale. These coefficients, ρ k, are analogous to the familiar Pearson correlation coefficient, but each is specific to the component of variation represented by the basic variogram function, gk(h), and so is specific to a particular spatial scale of variation (Goovaerts and Webster, 1994). They are estimated from the sill variances by

( kuv =

buvk k k buu bvv

.

(5)

Results and interpretation k Table III lists the coefficients buv for both K and Th with effective ranges 3a = 49·4 km and 3a = 16·4 km respectively. Figure 3 shows the models fitted to the experimental values. For both the soil survey and radiometric data, K shows strong spatial correlation to 49 km, with sill variances of 0·91 and 0·4 respectively. Their cross-correlation is strong, with the variogram model lying close to the hull (Figure 3(c)). The proportion of nugget variance in the cross-variogram (<1%) is much smaller than in the auto-variograms (8 and 11%) (Figure 3(a)–(c)), which suggests that the nugget variances in the auto-variograms derive largely from measurement errors. If such errors are uncorrelated then they disappear from the cross-variogram. The spatial correlation between Th by the two methods is also strong, but with larger nugget variances of 53 and 13% for the soil survey and radiometric data respectively. The cross-correlation for Th is much weaker than for K, with a range of approximately 16 km and with the fitted model lying far from the hull in the cross-variogram (Figure 3(f)). Again, the relatively small nugget variance in the cross-variogram of Th (2%, Figure 3(f)) suggests that this is largely measurement error. The significant correlation coefficients in the exponential component of the models of both K (r = 0·95) and Th (r = 0·69) demonstrate that their spatial relations, as determined by each of the two methods, are strong at this scale, and small in the nugget component (r = −0·009 and 0·04 respectively). A previous study (Rawlins et al., 2003) suggested that long-range spatial structure in topsoil geochemistry derives from the extensive outcrops of the main parent materials.

Published in 2007 by John Wiley & Sons, Ltd.

Earth Surf. Process. Landforms 32, 1503–1515 (2007) DOI: 10.1002/esp

1510

B. G. Rawlins, R. M. Lark and R. Webster Table III. Variances b kuv and structural correlation coefficients, (k, for the coregionalization of K and Th for soil survey and radiometric data

Variances K (SS) K (rad) K (SS) × K (rad) Th (SS) Th (rad) Th (SS) × Th (rad) Structural correlations (( k) K (SS) × K (rad) Th (SS) × Th (rad) a

Nugget, k=1

Exponentiala k=2

0·07 0·042 −0·005 2·083 0·045 0·045

0·841 0·358 0·52 1·877 3·696 1·817

−0·09 0·04

0·95 0·69

Effective ranges: 3a = 49·4 km for K; 3a = 16·4 km for Th.

Now it seems that the long effective range for K (49 km) has the same underlying cause, namely the variation in the K-bearing minerals between the main parent material types. In contrast, the Th-bearing accessory minerals are distributed in a more haphazard way, with a shorter correlation range (16 km). Ideally, an analysis of coregionalization would be based on data for the same depths in the soil. However, in the present case this was not possible; the soil was sampled for geochemistry between 35 and 50 cm depth according to the national protocol, which we could not alter, whereas the radiometric signal derives from the upper 35 cm. Nevertheless, the strong cross-correlation between the two methods suggests that the effect of the difference in depth on the comparison of K is small. For Th, the larger nugget variances suggest that measurement error is a greater source of error than the different depth ranges.

Modelling gamma emissions with respect to soil and parent material The linear mixed model. We explored the possible sources of variation in the passive gamma signal from the land by examining the relations between the K and Th signals, and their concentrations in the subsoil and the type of parent material. For this purpose we used the linear mixed model (LMM), which we can write as z = Xτ + Yη + ε.

(6)

Here the vector z contains our n observations of the gamma signal. Matrix X is an n × p design matrix that associates each of the n observations with a value of each of the p fixed effects (in this case a set of dummy variables that identify the parent material class at each site or the K or Th content of the subsoil or both). The vector τ contains the p fixed-effect coefficients. The vector η contains q random effects, realizations of a variable η, that are associated with the n observations by the n × q design matrix Y, which here is the identity matrix with n = q. We assume that η is a spatially correlated random variable and ε is a vector of independent random errors. These terms are independent of each other, and so we may write ⎡η⎤ ~ N ⎛ ⎡0⎤, ⎜ ⎢0⎥ ⎢⎣ε ⎥⎦ ⎝⎣ ⎦

⎡σ 2ξ R 0 ⎤⎞ ⎢⎣ 0 σ 2 I⎥⎦⎟⎠

(7)

where σ 2 is the variance of the independent error, ξ is the ratio of the variance of η to σ 2 and R is the correlation matrix of η. Note that we make an explicit assumption that the random terms are jointly Gaussian. The elements of ε, the nugget variances, represent both independent measurement errors and variation that arises from processes that are spatially dependent over shorter distances than those that separate the closest pairs of sampling points. Under the assumption that η is drawn from a second-order stationary random process, the correlation matrix R will depend only on the relative locations of our observations given some specified correlation function C(·) with one or more parameters that characterize the spatial dependence; so the elements of R are Published in 2007 by John Wiley & Sons, Ltd.

Earth Surf. Process. Landforms 32, 1503–1515 (2007) DOI: 10.1002/esp

Airborne radiometric survey signals

1511

Figure 3. Auto- and cross-variograms for (a) K (SS), (b) K (rad), (c) K (SS) × K (rad), (d) Th (SS), (e) Th (rad) and (f ) Th (SS) × Th (rad). The symbols are estimates of the experimental and cross-variances, and the solid lines show the models fitted to them. The dashed lines in the autovariograms show the variance in the data. Solid lines show the hulls of perfect correlation in the cross-variograms.

ri, j = Corr[η(si), η(sj)] = C(si − sj).

(8)

The correlation function may be one of several authorized functions, such as the spherical, or exponential function, the corresponding variogram function of which is given in Equation (3), Published in 2007 by John Wiley & Sons, Ltd.

Earth Surf. Process. Landforms 32, 1503–1515 (2007) DOI: 10.1002/esp

1512

B. G. Rawlins, R. M. Lark and R. Webster

⎧⎪ − si − s j ⎫⎪ C(si − s j ) = exp ⎨ ⎬ a ⎪⎭ ⎪⎩

(9)

is which a is a single distance parameter, which must be estimated. Note that both this function and the spherical function, which also has only the one distance parameter, describe isotropic variation where the variogram depends only on the distance between si and sj. The correlation function could be more complex with parameters that describe spatial anisotropy, but our exploratory analysis suggested that such elaboration was unnecessary in this instance. We estimated the parameters of the exponential function, which we represent by a, along with σ 2 and ξ, by residual maximum likelihood (REML). The residual log-likelihood function, Equation (10) below, has the unknown terms σ 2, ξ and a as its arguments, conditional on the data z. The residual log-likelihood is


R (σ 2 , ξ, a|z) = −

1 ⎧ 1 T T 2 −1 T ⎫ ⎨ln | H | + ln | X HX | + (n − p)σ + 2 z (I − WC W )z⎬ σ 2 ⎩ ⎭

(10)

where W = [X, Y] and H = ξ YRyT + I. The estimates of σ 2, ξ and a that maximize
R(σ 2, ξ, a|z) are found numerically. This REML solution removes dependence of the estimates of the parameters σ 2, ξ and a on the fixed effects τ, which are nuisance parameters in this problem and which would increase the bias of estimates based on maximum likelihood or method of moments (Smyth and Verbyla, 1996). Once the variance parameters are estimated we can obtain estimates of the fixed effects (τ) by generalized least squares. Lark and Cullis (2004) and Lark et al. (2006) give the details elsewhere. There are many places in the G-BASE inventory where there are both measurements of the total K or Th content of the subsoil and radiometric values. With so many data an efficient algorithm is essential for finding estimates of the variance parameters by maximization of Equation (10). We use the ASReml program (Gilmour et al., 2002), which embodies the average information (AI) algorithm of Gilmour et al. (1995). The algorithm is efficient, but it is not suitable for estimating variance parameters with a spherical correlation function. This is because the corresponding likelihood function is not smooth (Mardia and Watkins, 1989). For this reason we considered only the exponential correlation function. The following fixed effects models were considered to predict the radiometric measurement. 1. 2. 3. 4.

The The The The

mean mean mean mean

value only. and parent material. and the measured soil content (K or Th). value, parent material and soil concentration.

We compare the models by reference to their variance models for the random effect u and the independent error η. Results and interpretation. Figure 4 shows the raw data and the predicted results for each location under the fixed effect models (i.e. Xτ ) with more predictors than the overall mean. The variograms for the unexplained variation of each model were computed as ⎛ h⎞ γ (h) = σ 2 + ξσ 2 exp − ⎝ a⎠

(11)

where h is the lag distance and σ 2, ξ and a are the REML estimates of the variance parameters a. The variograms are shown in Figure 5. Note that all these variograms have very similar nugget terms (the variance of the independent error), but that the sills, the asymptotically bounding variances, differ. None of the terms considered account for any of the independent variation, which is consistent with our view that this is due almost entirely to noise in the sensor. The sill variances are largest for the model with the mean as the only fixed effect (equivalent to the variogram of the raw data). In the case of K, including parent material in the fixed effect model reduces the unexplained variation, but not as much as including the soil’s K content (Figure 5(a)). In contrast, parent material accounts for more of the variance than soil survey Th content (Figure 5(b)). These changes in the variogram show the results of including terms in the fixed effect model. Parent material and soil concentrations, the two fixed effects in the model, are likely to be correlated because the soil on each parent material will have its particular K or Th content. However, the reduction in the sill variance arising Published in 2007 by John Wiley & Sons, Ltd.

Earth Surf. Process. Landforms 32, 1503–1515 (2007) DOI: 10.1002/esp

Airborne radiometric survey signals

1513

Figure 4. Raw data for soil K (%) (a) and the predicted results for each location based on fixed effect plus parent material (b), fixed effect plus soil K (c) and fixed effect plus soil K and parent material (d).

from the inclusion of both fixed effects suggests that there are independent effects on the signal too, and the model predictions confirm this. For example, when we examine the predictions based on soil K only (Figure 4(c)) it is clear that there are variations in the region with the largest gamma signal, over the Triassic (Mercia) mudstone outcrop, that are not well predicted unless parent material is included in the model (Figure 4(d)). Similarly, the predictions based on parent material alone (Figure 4(b)) are dominated by the contrast between the large concentrations of K in the Triassic (Mercia) mudstone and the typically smaller concentrations in soil over the other types of parent material. When soil K is included (Figure 4(c), (d)) the model accounts for more of the variation observed over the Magnesian limestone (the western margin of the region) and the Jurassic limestones in the east of the region.

Conclusions The overall correlations of radiometric estimates with soil survey data for K and Th were large, as were the structural correlations for components of variation spatially dependent up to 49 and 16 km for K and Th respectively. This suggests that the remotely sensed signals for these two elements provide effective estimates of soil geochemical composition. The result in itself is unlikely to be of interest to soil surveyors (unless they want specific geochemical data), nor are the total soil K and Th contents of direct agronomic significance. However, the geochemical signal could be a useful surrogate for factors of soil formation, which our further results suggest to be the case. Large proportions of the spatial variation in the airborne signals for K and Th were accounted for by a simple classification of the parent material. For K, a further proportion was explained when the measured K in the soil between 35 and 50 cm depth was included as a predictor; so, although much of the radiometric signal does appear to be explained by parent material (Jenny’s factor P) as described in our classification, there is also a substantial proportion of the variation within these classes that can be explained by geochemistry. This suggests that subtler changes should be detectable. Much of the Th signal can be accounted for by parent material differences, a result that confirms our expectation that radiometric signals for K and Th provide information on parent material variations in the young landscapes of England and Wales. As above, we should not be surprised, because the dominant clay minerals in the soil contain K (Loveland, 1984). We are confident therefore that airborne radiometric surveys can aid the production of large scale (1:50 000) thematic soil maps, particularly those related to soil texture across England and Wales. However, the relation between airborne radiometric signals and important soil properties such as texture needs to be examined further across a larger range of parent materials and ages of landscape to see how general it is. Published in 2007 by John Wiley & Sons, Ltd.

Earth Surf. Process. Landforms 32, 1503–1515 (2007) DOI: 10.1002/esp

1514

B. G. Rawlins, R. M. Lark and R. Webster

Figure 5. Variogram models fitted to the fixed effect models in which the terms are the mean value (µ), PM (parent material) and soil concentrations for (a) Ks and (b) Ths.

An airborne radiometric survey is also likely to be effective for digital soil mapping in similar, youthful landscapes where K-bearing minerals dominate the clay size-fraction. Although we chose our study region to limit the impact of other variables, by collating data on rainfall, relief, slope length and land use in a more complex landscape model we could investigate the relative importance of Jenny’s other factors of soil formation (O, R, Cl) where time (T) is constant for much of the UK, using the airborne gamma signal.

Acknowledgements We thank all the staff of the British Geological Survey and volunteers who collected the soil samples and analysed them in the GBASE project, and staff who recorded and processed the signals from the radiometric survey. This paper is published with the permission of the Director of the British Geological Survey (Natural Environment Research Council). R. M. Lark’s contribution was supported by the Biotechnology and Biological Sciences Research Council of the United Kingdom through its core strategic grant to Rothamsted Research. Published in 2007 by John Wiley & Sons, Ltd.

Earth Surf. Process. Landforms 32, 1503–1515 (2007) DOI: 10.1002/esp

Airborne radiometric survey signals

1515

References Avery BW. 1990. Soils of the British Isles. CAB International: Wallingford. Cook SE, Corner RJ, Groves, PR, Grealish GJ. 1996. Use of airborne gamma radiometric data for soil mapping. Australian Journal of Soil Research 34: 183 –194. Department of the Environment. 1992. The UK Environment. Her Majesty’s Stationery Office: London. Duval JS, Cook B, Adams JAS. 1971. Circle of investigation of an airborne gamma-ray spectrometer. Journal of Geophysical Research 76: 8466–8470. Fuller RM, Groom GB, Jones AR. 1994. The land cover map of Great Britain: an automated classification of Landsat Thematic Mapper data. Photogrammetric Engineering and Remote Sensing 60: 553–562. Gilmour AR, Gogel BJ, Cullis BR, Welham SJ, Thompson R. 2002. ASReml User Guide, Release 1.0. VSN International: Hemel Hempstead. Gilmour AR, Thompson R, Cullis BR. 1995. Average information REML: an efficient algorithm for variance parameter estimation in linear mixed models. Biometrics 51: 1440–1450. Goovaerts P, Webster R. 1994. Scale-dependent correlation between topsoil copper and cobalt concentrations in Scotland. European Journal of Soil Science 45: 79–95. Hoskin PWO, Schaltegger U. 2003. The composition of zircon and igneous and metamorphic petrogenesis. In Reviews in Mineralogy and Geochemistry: Zircon, Hanchar JM, Hoskin PWO (eds). Mineralogical Society of America: Washington, DC; 27–62. Jenny H. 1941. Factors of Soil Formation: a System of Quantitative Pedology. McGraw-Hill: New York. Lagacherie P, McBratney AB, Voltz M (eds). 2006. Digital Soil Mapping: an Introductory Perspective, Developments in Soil Science 31. Elsevier: Amsterdam; 350. Lark RM, Cullis BR. 2004. Model-based analysis using REML for inference from systematically sampled data on soil. European Journal of Soil Science 55: 799–813. Lark RM, Cullis BR, Welham SJ. 2006. On spatial prediction of soil properties in the presence of a spatial trend: the empirical best linear unbiased predictor (E-BLUP) with REML. European Journal of Soil Science 57: 787–799. Lark RM, Papritz A. 2003. Fitting a linear model of coregionalization for soil properties using simulated annealing. Geoderma 115: 245– 260. Loveland PJ. 1984. The soil and clays of Great Britain: 1. England and Wales. Clay Minerals 19: 681–707. Mardia KV, Watkins AJ. 1989. On multimodality of the likelihood in the spatial linear model. Biometrika 76: 289–295. Mayr TR, Palmer R, Lawley R, Fletcher P. 2001. New Methods of Soil Mapping, Final Report for DEFRA Contract SR0120. Cranfield University: Bedford. McKenzie NJ, Gallant J. 2006. Digital soil mapping with improved environmental predictors and models of pedogenesis. In Digital Soil Mapping: an Introductory Perspective, Developments in Soil Science 31, Lagacherie P, McBratney AB, Voltz M (eds). Elsevier: Amsterdam; 327–349. McKenzie NJ, Ryan PJ. 1999. Spatial prediction of soil properties using environmental correlation. Geoderma 89: 67–94. Ministry of Agriculture, Fisheries and Food (MAFF). 2000. England Rural Development Programme 2000–2006: Appendix A4 – East Midlands Region. MAFF: London. Minty BRS. 1997. Fundamentals of airborne gamma-ray spectrometry. AGSO Journal of Australian Geology and Geophysics 17: 39–50. Peart RJ, Cuss RJ, Beamish D, Jones DG. 2003. The High Resolution Airborne Resource and Environmental Survey Phase 1 (HiRES-1): Background, Data Processing and Dissemination and Future Prospects, British Geological Survey Internal Report IR/03/112. British Geological Survey: Nottingham. Rawlins BG, Webster R, Lister TR. 2003. The influence of parent material on top soil geochemistry in eastern England. Earth Surface Processes and Landforms 28: 1389–1409. Reimann C, de Caritat PD. 1998. Chemical Elements in the Environment: Factsheets for the Geochemist and Environmental Scientist. Springer: Berlin. Smyth GK, Verbyla AP. 1996. A conditional likelihood approach to residual maximum likelihood estimation in generalized linear models. Journal of the Royal Statistical Society, Series B 58: 565–572. Stanley CR. 1998. NUGGET: PC software to calculate parameters for samples and elements affected by the nugget effect. Exploration and Mining Geology 7: 139–147. Wilford J. In press, 2007. Remote sensing with gamma-ray spectrometry. In Australian Soil and Land Resource Handbook: Guidelines for Surveying Soil and Land Resources (2nd ed.). McKenzie NJ, Webster R, Grundy MJ, Ringrose-Voase AJ (eds). CSIRO Publishing: Melbourne. Wilford J, Minty B. 2006. The use of airborne gamma-ray imagery for mapping soils and understanding landscape processes. In Digital Soil Mapping: an Introductory Perspective, Developments in Soil Science 31, Lagacherie P, McBratney AB, Voltz M (eds). Elsevier: Amsterdam; 207–218.

Published in 2007 by John Wiley & Sons, Ltd.

Earth Surf. Process. Landforms 32, 1503–1515 (2007) DOI: 10.1002/esp