A comparison of small strain stiffness in till as measured by seismic refraction and barometric loading response

Soil stiffness can vary over several orders of magnitude depending on the actual range of strain imposed by testing, or as a result of operational strains in geotechnical structures. Soil stiffness changes rapidly with strain level at low strain levels (0.01–0.1%) and the variation with strain is not linear. Characterization of the in situ small strain stiffness of stiff soils is important in geotechnical design; however, analyses of the mechanical behaviour of these soils is confounded by stiffness values that vary with strain level. Stiff till cuttings are susceptible to progressive failure as a result of strain softening. As a consequence, the evolution of stiffness during progressive failure is both a key parameter in characterizing pre-failure slope deformations and a key diagnostic of softening. Changes in strength (due to softening) should be reflected in commensurate temporal and spatial changes in stiffness; consequently, real-time, in situ measurements of stiffness would better define the progression of softening. Seismic surveys, which create small compression and shear strains, have been used to estimate in situ small strain elastic moduli. These spatially extensive measurements can be correlated to temporal variations in stiffness from the monitoring of barometric loading efficiency. In this latter method, the pore pressure response of a grouted (sealed) piezometer to barometric pressure fluctuations is used to measure the compressibility (stiffness) of the formation. This article summarizes the results of field trials within a cutting in stiff till in Northern Ireland in which these two techniques were used to characterize small strain stiffness.


INTRODUCTION
The road and rail network in Northern Ireland (NI) encompasses a significant number of large cuttings in heavily overconsolidated, stiff, lodgement till. The stress-strain regime of these cut slopes are complex, with different principle stress directions at different positions along the potential failure plane. For example, loading may be primarily in extension near the toe of the slope, whilst compressive loading is predominant at the crest of a slope. Cuttings in overconsolidated clay are known to be susceptible to progressive failure, with softening of the toe and the development of a rupture surface into the slope, as is well documented for London Clays (Potts et al., 1997).
The mechanism of progressive, or delayed, failure has been attributed to many factors including the dissipation of pore water pressures generated during an excavation and the strain softening associated with climatically driven pore pressure cycling (Potts et al., 1997;Hughes et al., 2007). If the long-term stability of cuttings in lodgement tills is susceptible to this failure mechanism, methods of assessing the change in in situ stiffness over time would provide a valuable indicator of slope condition for geotechnical asset owners.
The measurement of small strain elastic properties of soils are most often carried out in the laboratory using bender elements installed in triaxial cells (Woods, 1994;Lings and Greening, 2001;Donohue and Long, 2010;Bonal et al., 2012); however, even 'undisturbed' samples of soils tested in this manner are subject to sufficient disturbance to compromise the measurements. Laboratory testing on reconstituted samples of tills has been carried out to investigate the time-dependent behaviour of the NI tills . However; it was found that high and variable gravel contents makes it impractical to obtain good quality undisturbed samples from the field for accurate laboratory testing. This, combined with the time frames required for sample preparation and testing, suggests that alternative methods are required to determine the in situ stiffness of the till.
On site, pressuremeters can be used to measure in situ strength and stiffness parameters of soils and rocks, which cannot be measured with typically used push equipment (for example, Standard Penetration Tests and Cone Penetration Tests). Pressuremeters enter the ground by pushing, by pre-boring a hole into which the probe is placed, or by self-boring, whereby the instrument makes its own hole. Increments of pressure are applied to the inside of the membrane, forcing it to press against the surrounding ground, and so loading a 'cylindrical cavity'. The test will report on a series of pressure readings and the respective horizontal displacements of the cavity wall. The method of installation causes irrecoverable disturbance; however, the self-boring method reportedly causes disturbance which is small enough to lie within the elastic range of the material. It should be noted that these instruments are complex by conventional site investigation standards, can only be operated by trained personnel, inappropriate analysis can lead to misleading parameters, and it is comparatively expensive to other site investigative techniques.
Seismic surveys induce very small strains and it is generally accepted that they provide results relevant to the linear-elastic phase of soil deformation (Whiteley, 1994;Matthews et al., 2000;Foti, 2003;Soupios et al., 2006;Long and Menkiti, 2007;Clayton & Heymann, 2001). It has also been shown that laboratory stiffness levels are similar to field seismic measurements (Michaels, 1998;Clayton, 2011).
Barometric loading efficiency can also be used to quantify the stiffness of stiff confined aquifers or aquitards. The barometric loading efficiency is the ratio of the changes in pore-pressure to the changes in barometric pressure. It is a simple function of the stiffness of the formation as well as the porosity and stiffness of the fluid. Since the latter are known, the stiffness of the formation can be calculated directly from the barometric loading efficiency (Rojstaczer, 1988;van der Kamp and Schmidt, 1997;Anochikwa et al., 2012, Smith et al., 2016Smith et al., 2013).

Research Site
The research site used for this study is near Loughbrickland in County Down, Northern Ireland ( Figure 1). Loughbrickland cutting is a 25 m high, east facing slope, with an average slope angle of approximately 26 °. It was excavated in 2004 through drumlin topography in order to improve the horizontal realignment of the A1 carriageway. Borehole logs (Clarke, 2007) were used to infer the topography of a cross-section through the cutting ( Figure 2). The bedrock is completely to moderately weathered Greywacke sandstone with completely weathered slaty mudstone interbeds, typical of the Palaeozoic, Silurian Gala Group bedrock geology of the area (Anderson, 2004). The bedrock is overlain by a lower and upper till layer, differentiated primarily by observed differences in hydraulic conductivity  with an upper weathered zone primarily associated with rooting and seasonal wetting and drying cycles ( Figure 2).
The reader is refered to Figure 2 and 3 for detailed figures of the site, including the location and identification numbers of the installed piezometers at the research site. A number of itmsoil heavy duty vibrating wire piezometers (VWP) have been installed using the fully grouted method at the site at a range of depths ( Figure 2); the itmsoil VWPs measure pressure ranges from -30 to 300 kPa. Barometric pressure changes are measured using a VWP at the surface, open to atmospheric air pressure. The site is also fully instrumented to monitor surface water balance and infiltration, water table elevation and climatic conditions. For more detail on the site see the following references Clarke, 2007;Carse, 2014;McLernon, 2014 contents were calculated for the matrix after discarding the stony material retained on the 5 mm sieve. The water content with the gravel or large particles sizes included are also presented for comparison to highlight that the matrix (stone-free) water contents are generally 30-50 % higher than the intact sample water contents (Table 1). This correction was carried out as the author believes that the Till behaviour is dominated by the clay matrix and, given the variability of the stone content in the Tills, the matrix water content was more indicative of the soil behaviour. A summary of the Particle Size Distribution (PSD) ranges for the Loughbrickland site is also presented in Table 1, alongside typical PSD ranges for another cutting site along the A1 at Dromore as well as typical ranges for Dublin Boulder Clays. It can be seen that the material is variable in nature, but with clay contents between 16 and 26 %.
During excavation of the cutting, large inclusions of soils with higher clay content were also observed Kelly, 2018). These inclusions reinforce the appreciation of the highly heterogeneous nature of drumlin formations.

SURFACE WAVE SURVEY
Multichannel analysis of surface wave (MASW, Park et al., 1999) surveys were carried out at the Loughbrickland research site from May 2014 to assess the geological structure of the site (i.e. estimating the thickness of the till, identifying any inhomogeneity in its materials, mapping the surface of the underlying bedrock) and to offer an independent alternate estimation of the small strain soil stiffness. Surface wave (SW) surveys have been established as a reliable tool for the estimation of shear wave velocity (Vs) and small strain shear modulus (G) in the near surface (Nazarian and Stokoe, 1984;Gabriels et al., 1987;Foti, 2003).

Data Acquisition
Surface wave (SW) data were acquired along three transects, located at the crest of the cutting, at the toe, and on a berm crossing the slope transversely, at approximately half of the slope height ( Figure 3). Data acquisition was achieved using 24 no. 4.5 Hz vertical geophones.
Receivers were arranged in linear arrays with a regular spacing of 2 m (at the crest of the slope), 1 m (at the berm) and 0.25 m (at the toe), following the necessity to ensure greater investigation depths (proportional to the length of the seismic transect, (Socco and Strobbia, 2004) as the elevation increases. The seismic source was a 4.5 kg sledgehammer (crest and berm acquisitions) or a 450 g claw hammer (toe acquisition) hitting a metal plate placed on the ground surface. The seismic source was positioned at both ends of the recording arrays (Socco and Strobbia, 2004).
A key issue encountered during the acquisition of the seismic data was that the site is located beside a busy dual-carriageway; therefore, a high level of undesired noise due to traffic was present in the collected seismic traces. This problem was mitigated by operating the seismic source multiple times at each shooting position and then stacking (i.e. summing) the sets of data with the same source-receiver configuration, so as to maximize the ratio between the signal generated by the sledgehammer and the incoherent noise (Foti et al., 2015).

Data Processing
Recorded seismograms ( Figure 4a) were translated to the frequency-wavenumber (f-k) domain, where SW frequency components can be easily identified due to their greater energy content ( Figure 4b). Energy maxima were picked (white asterisks in Figure 4b), thus obtaining an experimental dispersion curve (blue asterisks in Figure 4c), i.e. a curve correlating each SW frequency components to its corresponding phase velocity (v), which is obtained from: The frequency-dependence of SW phase velocity shown in Figure 4c is caused by the dispersive properties of surface waves. According to this property, SW frequency components propagate in a portion of the subsurface whose thickness is approximately equal to one wavelength (Socco and Strobbia, 2004). Hence, in a vertically heterogeneous medium different frequency components having different wavelengths exhibit different velocities of propagation. Moreover, the propagation of SW is a multimodal phenomenon, so that each frequency component can exhibit more than one velocity, each of which corresponds to a different mode (as shown in Figure 4c).
The processing operations of the surface wave method were applied separately to seismic data from the crest, berm and toe seismic arrays, thus obtaining three experimental dispersion curves ( Figure 4c shows the dispersion curve for the transect at the crest of the slope). The dispersion curve constitutes the set of experimental observations from which the parameters of the postulated subsurface models are estimated, in an operation called inversion (Lines and Treitel, 1984). During the inversion process, the experimental dispersion curve is collated with a set of synthetic curves corresponding to layered subsurface models whose parameters (See Table 2: thickness, ρ, Poisson's ratio and VS of various layers) are assumed (Socco and Strobbia, 2004). The operation allowing a synthetic dispersion curve to be computed from a given subsurface model (known as forward modelling) is available in semi-analytic form (Thomson, 1950;Haskell, 1953). The models, whose corresponding synthetic curves best fit the experimental dispersion curve represent the inversion result, and hence the estimate of the properties of the investigated subsurface.
The inversion strategy adopted for the present work is a Monte Carlo inversion procedure as developed by Maraschini and Foti (2010). Each experimental dispersion curve was collated with two 106 synthetic curves, corresponding to as many layered Vs subsurface models whose parameters (thickness and S-wave velocity (Vs) of various layers) were randomly selected.
Since SW data are poorly sensitive to the values of Poisson's ratio and ρ (Xia et al., , 2003, these were assigned according to available geotechnical information and data from literature (Table 2: Clarke, 2007;Carse, 2014;McLernon, 2014), as is common practice in SW inversion (Socco and Strobbia, 2004). The inversion algorithm employed by Maraschini and Foti (2010) allows a rapid computation of the misfit between experimental and synthetic curves; moreover, it efficiently includes all available SW dispersion modes in the inversion considerably improving the robustness of the results.

Discussion of Data
One of the main sources of error in the characterisation of near surface via the SW methods is in the incorrect identification of the different modes on the experimental dispersion curves Zhang and Chan, 2003;Maraschini and Foti, 2010). Additionally, as suggested by Karray and Lefebvre (2012) a correct identification of SW modes can be hampered by mode jumps, caused by sharp velocity contracts or velocity inversion (Maraschini and Foti, 2010). The Monte Carlo code adopted for the inversion (Maraschini and Foti, 2010) tackles these issues by not requiring any prior attribution of dispersion curve points to the SW modes, searching all subsurface models that can explain the experimental curve, considering the possibility of mode jumps or apparent dispersion modes. In the case of the Loughbrickland site, the lack of sharp velocity contrast between the lower till formation and the underlying bedrock can be ascribed to the fact that the intact bedrock is overlain by a weathered bedrock layer. Figures 4c displays the comparison between the experimental and the best fitting synthetic curves; Figure 4d shows the corresponding subsurface models, selected as statistically equivalent to the lowest misfit Vs profile by applying a Fisher test with a 95 % confidence level (Socco and Boiero, 2008;Maraschini and Foti, 2010). These profiles constitute the inversion result (i.e. the estimation of VS structure) for the seismic array placed at the crest of the slope. Despite the SW method (Foti et al., 2009), the selected VS profiles, extracted from a population of two 10 6 different models, show a good reciprocal consistency, thus confirming the reliability of the obtained solution.

Results
The best fitting profiles form Figure Carr et al., 1998;Long and Menkiti, 2007). At around 25 m, Vs sharply increases again towards values of 1250 m/s, witnessing the location of the upper boundary of the weathered bedrock as identified in borehole logs (Clarke, 2007).
From the VS profiles, it is possible to derive an estimate of the shear modulus at small strains (G) by applying: where ρb is the soil bulk density. From Vs models in Figure 4d, a single representative profile was derived from the S-wave velocity models shown in Figure 4d and used to calculate a G0 modulus profile. Since the various profiles of Figure 4d belong to the same statistical class (according to the applied Fisher test, Socco and Boiero, 2008), the single VS model was obtained by manually drawing a profile following the global trend of accepted models, as is common practice in Monte Carlo results interpretation (e.g. Konstantaki et al., 2013). The resultant Vs profile (blue dashed line in Figure 4d) Carse, 2014;McLernon, 2014). The resultant G0 model was used to deduce an elastic stiffness (E) modulus profile using the relationship: In order to deduce an accurate estimate of Poisson's ratio, using the relationship between Vs and Vp profiles correlates to a Poisson's ratio (Richart et al., 1970). The Vs/Vp ratio was obtained by collating the Vs profiles from the surface wave survey conducted in this study (Figure 4), and the Vp model from a P-wave refraction survey previously performed at the same site (Carse, 2014 (Santamarina, 2001), or the presence of layers originating from different depositional processes (as seems to be the case for the Loughbrickland site).

BAROMETRIC LOADING EFFICIENCY
Changes in barometric pressure at the ground surface will act as a change in a uniform, aerially extensive, change in total stress. This total stress change will be transmitted across the entire soil depth and will have to be picked up by the soil skeleton (e.g. effective stress change) and the pore water pressure (Anochikwa et al., 2012;van der Kamp and Maathuis, 1991; Kamp and Schmidt, 1997; Bardsley and Campbell, 2007). This pore pressure response is similar to that created by the construction of an aerially extensive surface fill that induces a pore water pressure response within an underlying saturated formation (Bishop, 1954;Skempton, 1954). This phenomenon was recognised early on within both the geotechnical (Skempton, 1954) and the hydrogeology disciplines (Jacob, 1940). For example, Jacob (1940) observed similar responses within a confined aquifer as a result of tidal loading. More recently, attempts have been made to also use this observation to track changes in surface loading as a result of changes in soil moisture loading as a result of changes in stored water within the soil profile (Anochikwa et al., 2012;Marin et al., 2010, van der Kamp andMaathuis, 1991;Barr et al., 2000).

Theoretical Background
The theoretical explanation for the pore pressure response to loading is based in the generalised constitutive relationships for coupled stress-strain and pore water pressure response of saturated, isotropic, linearly elastic, porous media as presented by Biot (1941) and Nur and Byerlee (1971). This formulation was described more simply for the case of barometric loading by van der Kamp and Gale (1983) for a laterally constrained domain, subject to a spatially continuous surface loading, such as barometric pressures (Anochikwa et al., 2012). These authors showed that the pore water pressure response to an instantaneous change in surface loading, referred to as the loading efficiency, can be described as follows: where ̅ is the loading efficiency, better known as Skempton's ̅ coefficient (Skempton, 1954), u is the pore water pressure, and σB is the barometric pressure. Theoretically, the loading efficiency for a linear, elastic, isotropic, saturated soil is related to the stiffness of the soil through the following equation: where Ec is the drained constrained modulus of elasticity, n is the porosity of the soil formation and Ew is the bulk modulus of water (kPa). In hydrogeology, these are better known as the reciprocals of the constrained formation compressibility (α, kPa -1 ) and the bulk compressibility of water (β = 4.8 x 10 -7 kPa -1 at 25 °C). Jacob (1940) and Skempton (1954) both developed similar expressions.
If it is assumed that the soil response to barometric loading is undrained then rearranging Equation (5) gives: The value of Young's Modulus (E) can be obtained using the relationship between E, Ec and ν (Poisson's Ratio) using the relationship developed by Poulos and Davis (1974) and later used by van der Kamp and Gale (1983): At the Loughbrickland site, Poisson's ratio (ν=0.45) was estimated from the relationship between wave velocities; see Section 2.4.

Estimating Loading Efficiency
There are a range of methods available in literature for the estimation of loading efficiency from pore water pressure records and barometric loading records. For example, when monitoring field data, there is a lag between the water level fluctuation and the barometric observations. Rasmussen and Crawford (1997) therefore used regression deconvolution to estimate a barometric response function (BRF) as the technique is useful in estimating how a parameter in a system responds to a stimulus when the response is not instantaneous, and the magnitude of the response changes with time (Toll and Rasmussen, 2006). Rasmussen and Crawford (1997) and Spane (2002) document the regression deconvolution method and its use as a diagnostic tool. This paper however presents the manual adjustment method, where the pore pressure data from the piezometers are corrected for barometric response at monthly time intervals. The barometric pressure was subtracted from the raw pore water pressures, and multiplied by a correction factor which was adjusted by trial and error to eliminate barometric effects as determined by visual inspection of the data (e.g. Smith et al., 2013;Anochikwa et al., 2012). The correction factor represents the loading efficiency 'B-bar' ( ̅ ); a worked example is presented in the subsequent section. It must be noted, that due to this method utilising trial and error, it therefore is subject to a level of uncertainty. Therefore, this was repeated on a monthly basis, and the data presented in Tables 3 and 4 presents an average Loading efficiency.

Loading Efficiency Analysis Example
Barometric pressure was measured at the research site using a sheltered piezometer at the logger box (18/01/2010 -02/08/2011) and was measured later using a Davis Vantage Pro weather station (30/04/2015 -22/11/2015). The raw pore water pressure data from the piezometers were corrected for barometric response by subtracting the changes of barometric pressure, multiplied by a correction factor that was adjusted by trial and error to eliminate barometric effects as determined by visual inspection of the corrected data, as described earlier in the article.   (Clarke, 2007). The pressure in each packer was controlled at one bar of pressure by an air reservoir which was inspected at regular intervals. On occasion the reservoir was found to have deflated and on rare occasions leaks were found in the air tubing; these sections of data were excluded from the analysis presented in this paper.

Discussion
The assessment of the in situ small strain stiffness of till deposits carried out using seismic refraction surveys and loading efficiency analysis in this research gives results which are comparable across three research sites in Northern Ireland (Carse, 2014), and are also consistent with literature ( Figure 8). For example, work by Clayton and Heymann (2001) places the small-strain undrained stiffness (E) of London Clay at 240 MPa. Long and Menkiti (2007) analysed areas of Dublin Boulder Clay using a number of laboratory and geophysical techniques and found its small strain shear modulus (Gmax) between 250 kPa and 1500 MPa (approx. E of 750-4200 MPa, with an assumed υ of 0.4). Figure 9 shows the typical strain ranges under which the majority of geotechnical structures operate (10 -5 to 10 -2 ), which fits well within the range of strain levels induced by geophysical tests, thus providing reliable estimates of small strain stiffness.
An interesting observation from the data presented in Figure 5 is the depositional layering at the Loughbrickland research site, showing the location of the upper and lower till. It is also worth noting that the stiffness at the toe of the slope is relatively low, which is to be expected from the construction of the cutting, and has likely been reduced due to the development of artesian conditions during construction in 2004 . Harley et al. (2016) also demonstrated how stiff till cuttings are susceptible to progressive failure as a result of strain softening. As a consequence, the evolution of stiffness during progressive failure is both a key parameter in characterising pre-failure slope deformations and a key diagnostic of softening. Changes in strength due to softening, is reflected in commensurate temporal and spatial changes in stiffness; consequently, real-time, in situ measurements of stiffness is a novel means of monitoring a slope's condition.
This research has shown how simultaneously carrying out barometric loading efficiency analysis, and seismic surveys are complimentary; for example, seismic surveys allow an accurate estimation of Poisson's ratio with depth, which is a key assumption in calculating the compressibility of a formation. Figure 6 demonstrated the reliability of the two methods adopted for this research, as both methods yielded similar estimates of elastic modulus (E) with depth at each of the three transect locations.

CONCLUSION
Both Multichannel Analysis of Surface Wave (MASW) and Barometric Loading Efficiency analysis gives similar and reasonable estimates of in situ small strain stiffness. The estimation of stiffness using the MASW analysis is dependent on a reasonable estimation of the bulk density (γ), and the estimation of stiffness using the barometric loading method is dependent on a good estimate of Poisson's Ratio (υ). This study showed that both methods were complimentary, in that Poisson's Ratio can be determined in situ with depth from the ratio of surface and compression wave velocities (Vs and Vp), and that manual dips in boreholes on the day of data acquisition allowed an accurate report of the water In situ soil stiffness profiles are required if reliable predictions of displacement patterns around new and existing infrastructure are to be made, particularly in urban areas. In industry, the measurement of soil stiffness profiles for a range of applied civil engineering research applications would be beneficial, for example deep foundations and retaining wall construction; as aforementioned, the small strain elastic modulus of stiff soils can vary over many orders of magnitude for small operational strains, which makes predicting ground displacement patterns due to tunnelling and deep excavations a difficult task. Derivation of the small strain constrained modulus, in conjunction with stiffness degradation data, will also be suitable for tunnelling and slope stability assessment. However, further research is required, which will require close links with industrial project partners from engineering consultants and contractors, to allow researchers to directly feed outputs of the research into applied research in practice, through existing projects. Suitable existing piezometer data, or re-activated monitoring sites under industrial partner ownership would allow data exploitation to gain valuable stiffness information.
It should be noted, that methods of measuring in situ elastic properties of soils and soft rocks are required in the accompanying fields of geography, geology, hydrogeology and resources extraction, and would therefore be beneficial across numerous sectors internationally. For example, Argillaceous sediments (clay-rich aquitards, commonly referred to as shales) with low hydraulic conductivities (<10 -8 ms -1 ) make up 2/3 of all sedimentary rocks on Earth (Smith et al., 2016). It is these Shales that typically control recharge and chemical transport to adjacent aquifers, as well as act as isolating units to protect shallow groundwater from contamination by fluid or gas migration from deeper formations. Management and protection of these groundwater resources is often dependent on accurate determinations of the geotechnical properties of formations, therefore Smith et al. (2016) utilised a barometric loading efficiency method to understand the stress behaviour of Argillaceous aquitards.
It should be noted that these are evolving techniques used in cutting-edge practice in Canada, for applications including accurate determination of hydrogeological properties of aquitard formations, and compressibility of aquitard formations which can influence yield to adjoining aquifers, subsidence from fluid extraction, propagation of stress changes (Smith et al., 2013) and in calculating specific storage (Cook et al., 2017). Hendry et al. (2017) utilise the barometric loading method to assess a site with dynamic groundwater conditions due to fluctuating river levels, which illustrates the robustness of the method in a dynamic environment which historically was difficult to characterise. Both techniques are currently underutilised, and could be more widely used in the UK and Ireland.    . Surface wave data from transect at the top of the slope: acquisition, processing and inversion. a) Seismic traces acquired by the array placed at the crest of the slope; b) Normalised f-k spectrum obtained from the seismic recordings at the crest of the slope (coloured background) and experimental dispersion curves (white asterisks, located on the energy maxima of the spectrum); c) Experimental dispersion curve (blue asterisks) in f-v domain compared with the best fitting synthetic dispersion curves (coloured lines, the colour depends on the corresponding misfit; d) Corresponding VS profiles: here again the colour refers to the misfit. The blue dashed profile represents the interpreted final model from which the G modulus, and therefore the E modulus profile is derived (Figures 5 and 6).