The UK’s transportation network is supported by critical geotechnical assets (cuttings/embankments/dams) that require sustainable, cost-effective management, while maintaining an appropriate service level to meet social, economic, and environmental needs. Recent effects of extreme weather on these geotechnical assets have highlighted their vulnerability to climate variations. We have assessed the potential of surface wave data to portray the climate-related variations in mechanical properties of a clay-filled railway embankment. Seismic data were acquired bimonthly from July 2013 to November 2014 along the crest of a heritage railway embankment in southwest England. For each acquisition, the collected data were first processed to obtain a set of Rayleigh-wave dispersion and attenuation curves, referenced to the same spatial locations. These data were then analyzed to identify a coherent trend in their spatial and temporal variability. The relevance of the observed temporal variations was also verified with respect to the experimental data uncertainties. Finally, the surface wave dispersion data sets were inverted to reconstruct a time-lapse model of S-wave velocity for the embankment structure, using a least-squares laterally constrained inversion scheme. A key point of the inversion process was constituted by the estimation of a suitable initial model and the selection of adequate levels of spatial regularization. The initial model and the strength of spatial smoothing were then kept constant throughout the processing of all available data sets to ensure homogeneity of the procedure and comparability among the obtained VS sections. A continuous and coherent temporal pattern of surface wave data, and consequently of the reconstructed VS models, was identified. This pattern is related to the seasonal distribution of precipitation and soil water content measured on site.

Many parts of the UK’s rail network were constructed in the mid to late 19th century, long before the advent of modern construction standards. Historic levels of low investment, poor maintenance strategies, and the deleterious effects of climate change have resulted in critical elements of the rail network being at risk of failure. A recent study conducted by Network Rail, who owns the largest network of earth structures in the UK, has revealed that 50% (or 5000 km) of their network of earthworks are considered to be in a poor or marginal condition (Network Rail, 2010). According to the Department of Transport Rail Accident Investigation, an average of 65 earthworks failures per year were recorded during the period from 2004/2005 to 2008/2009 (RAIB, 2008). The majority of these failures occurred as cutting and embankment slips, and a large portion were triggered by localized extreme weather events (Network Rail, 2010). One of the most recent examples is the Harbury landslide, in early 2015, in which 350,000 tons of material slipped along a 160 m long stretch leading to the closure of one of the main UK rail lines (Leamington Observer, 2015).

Over the last half century in the UK, there has been little change in the total winter precipitation; however, during this time there has been an increase in the proportion of winter precipitation resulting from extremely heavy rainfall events (Jenkins et al., 2009). Seasonal increases in pore water pressure (reductions in suction pressure) resulting from excessive rainfall may significantly reduce the strength of earthworks, leading to failure if triggered, for example, by an extreme weather event (Wilks, 2010). These failures are not only caused by elevated pore pressures due to prolonged and intense rainfall, but dry conditions can also lead to desiccation and cause cracking of slope surfaces, increasing their permeability and providing pathways for rainfall infiltration. When the effects of wetting and drying are combined, seasonal cyclical patterns may cause shrinking, swelling, and strain softening of the fill materials that may result in serviceability failures and eventual earthwork collapse (Kovacevic et al., 2001).

Advanced assessment and remediation of earthworks is significantly less costly than dealing with failures. It is therefore crucial to develop efficient approaches for the assessment of earthwork stability, so that mitigation measures can be targeted and, consequently, failures can be avoided. Geophysical methods, when used as complementary tools together with traditional geotechnical investigations, are ideally suited for earthwork investigations (Donohue et al., 2011, 2012; Chambers et al., 2014; Gunn et al., 2015a). In addition to being noninvasive, several techniques are cost effective and rapid, allowing detailed 2D and 3D volumetric information to be obtained, with repeated observations permitting temporal variations to be measured. This level of spatial and temporal detail cannot be captured using discrete boreholes or other forms of geotechnical investigation. Geophysical methods are particularly useful if the measured properties are calibrated with relevant geotechnical parameters.

A recognized body of research is emerging, focusing on the relationship among electrical resistivity, moisture content, and pore pressure (Chambers et al., 2014; Gunn et al., 2015a); however, to date, little work has been carried out using seismic geophysical methods to temporally and spatially monitor variations in the mechanical properties of slopes. S-wave (VS) and P-wave (VP) seismic velocities are related to the small-strain shear modulus G0 and constrained modulus M0 as follows:
where ρ is the soil bulk density and K0 is the small-strain bulk modulus. Therefore, in soft, fully saturated soils (saturation, Sr=100%), VP (unlike VS) is controlled by the bulk stiffness of the pore fluid; however, for a small reduction in saturation, VP decreases significantly; VS for any degree of saturation and VP in unsaturated or partially saturated soils (Sr<99%) are controlled by the state of stress, the degree of cementation, and processes that alter interparticle contacts, such as capillary pressures and electrical forces (Cho and Santamarina, 2001; Santamarina et al., 2005). The compacted soils used in earthworks are in an unsaturated state during construction and are likely to be in a partially saturated state during operating conditions. The geotechnical properties of these soils are therefore affected by the simultaneous presence of water and air in the pore spaces, which makes the pore fluid mixture compressible and influences the stress state through the water and air pressures (Mancuso et al., 2002). Temporal variations in the seismic velocities of partially saturated earthworks are therefore likely related to changes in capillary/suction pressures, which in turn may result from seasonally or climatically related variations in soil moisture. Several authors have explored the relationship among S-wave velocity, degree of saturation, and soil suction on laboratory samples for a range of materials (Qian et al., 1991; Marinho et al., 1995; Cho and Santamarina, 2001; Mancuso et al., 2002; Donohue and Long, 2010). The experimental results suggest that measurements of VS are sensitive to changes in saturation and soil suction. Although it is relatively easy to observe these variations in a controlled laboratory environment, little or no work has been performed to observe the consequence of climate effects on seismic velocities in the field. Despite the recognized importance of small strain behavior in predicting the performance of earthworks, little effort has thus far been devoted to the subject (Mancuso et al., 2002).

The analysis of surface waves is an established method for the reconstruction of S-wave velocity models for the near surface (Socco et al., 2010). In particular, the multichannel analysis of surface waves (MASWs) approach (Park et al., 1999; Xia et al., 1999; Zhang et al., 2004) is frequently used for a broad range of applications, from the estimation of soil stiffness (for instance, Foti, 2003), the evaluation of seismic site response (Yilmaz et al., 2006), reconstruction of complex subsurface structures (Ivanov et al., 2006), and the characterization of subsoil formations (as in the case of quick clays, e.g., Pasquet et al. [2014], or uncompacted granular materials, e.g., Gofer and Bachrach, 2012). Although less commonly used, surface wave (SW) attenuation data have been exploited for evaluating local seismic site response and soil amplification. The measurement of SW attenuation enables a near-surface material damping ratio profile to be determined (Rix et al., 2000). The developed techniques for the estimation of SW attenuation characteristics usually exploit the same testing configuration and raw seismic data as those methods used for SW phase velocity analysis. Both approaches can therefore be performed simultaneously (Rix et al., 2001; Lai et al., 2002; Foti, 2003).

Despite their widespread use, seismic geophysical techniques have rarely been used for reconstructing the time-lapse mechanical behavior of the subsurface, particularly on shore. The studies that have reported time lapse seismic measurements are limited to large-scale applications, such as hydrocarbon reservoir exploration geophysics (Landro, 2001; Lumley et al., 2003) or the monitoring of geologic storage reservoirs (Draganov et al., 2012; Picotti et al., 2012). In the field of near-surface geophysics, the potential application of seismic methods for tracking time-variant processes has been assessed (Jefferson et al., 1998; West and Menke, 2000), although at more reduced spatial and temporal scales than those presented in this paper. Recently, the use of seismic techniques has been proposed for the detection of varying levels of water table of shallow aquifers (Pasquet, 2014; Pasquet et al., 2015).

This study examines the use of SW and seismic attenuation data for tracking climate-induced variations in the geomechanical properties within the structure of an unstable railway embankment. To our knowledge, seismic methods have not previously been tested for assessing climate related, fluid-induced geomechanical property variations in slopes or earthworks.

In the following sections, we describe the surveyed site, the time-lapse acquisition of seismic data, and their processing in terms of extraction of SW dispersion and attenuation curves. The spatial and temporal variability of retrieved data sets, and the statistical significance of the latter with respect to experimental uncertainty, are also investigated. Finally, the SW inversion method is discussed and the resulting temporal and spatial variations in the S-wave velocity structure of the investigated embankment are presented and interpreted.

The railway embankment investigated in this study is part of the heritage Gloucestershire Warwickshire Railway and is located near the village of Laverton in the Cotswolds region of southern England. The embankment is approximately 5 m high, and it was constructed at the beginning of the 20th century via end tipping of local Charmouth Mudstone (Gunn et al., 2015b). The embankment comprises a 0.9 m thick upper layer of ballast fouled with fines, ash, and soil (rich in humus), underlain by approximately 4 m of a high-plasticity clay fill (Charmouth Mudstone). Several significant earthwork failures on the Gloucester Warwickshire Railway have occurred over recent years, and there is also evidence of previous movement at the Laverton embankment. Figure 1 shows a site layout map, and a stratigraphic log determined from intrusive investigation involving a combination of boring, sampling, trial pitting, and cone penetration testing (CPT). A CPT profile from the site is also included in this figure, and it indicates that the embankment may have been constructed in two lifts, with the boundary at approximately 3 m depth. In general, the upper interval is characterized by very low penetration resistance (qc, between 0.5 and 1 MPa), whereas the embankment fill below 3 m exhibits slightly greater qc (between 1 and 1.3 MPa). Samples of the clay fill extracted near CPT 3 (Figure 1) indicated high levels of moisture and weathering immediately beneath the ballast layer. Below the embankment base (approximately 5 m), qc increased considerably.

Data acquisition

The SW seismic data were acquired bimonthly along the crest of the railway embankment, from July 2013 to November 2014, for a total of nine acquisitions. For all acquisitions, we adopted the same source type, receiver array, recording equipment, and acquisition geometry/parameters, to make the acquired seismic data fully comparable. To ensure that sources and receivers were positioned consistently throughout all acquisitions, we placed permanent marks for spatial reference along the surveyed area.

The MASW (Park et al., 1999) data were acquired along the embankment crest using a landstreamer, consisting of 24, 4.5 Hz vertical geophones. The geophones were spaced at 1 m intervals and were coupled to the ground by means of cleated metal plates. Data were recorded every 5 m in a roll-along fashion, and, in total, a 100 m stretch of the embankment was investigated (21 acquisitions) (Figure 1). A sampling rate of 62.5  μs was used. For each setup, three recordings were acquired, which were subsequently stacked in the time domain to improve the signal-to-noise ratio (S/N) (Socco and Strobbia, 2004). The seismic source was a 4.5 kg sledgehammer hitting a metal plate located 2 m from the first receiver. A total acquisition time of less than 2 h was generally achieved. In addition to the SW data described in this paper, time-lapse P-wave refraction surveys were also conducted at the site, which will be described in detail in a later paper by the same authors.

For the purpose of stacking multiple single-shot data recorded at the same spatial location (Foti et al., 2015), an efficient and precise automatic triggering system was used, which consisted of a starter sensor attached to the hammer head. Synchronization of the recorded seismograms before stacking was further optimized using a control algorithm based on the crosscorrelation of the seismic traces recorded by the receiver closest to the source position (offset=2  m).

In Figure 2 (left column), three typical stacked seismograms acquired in January, March, and September 2014, which had the same source-receiver configuration, are shown (source located at x=41.5  m along the seismic line). The P-wave first-break arrivals and the SW train are marked with P and SW in the three seismic sections. Despite the similarity between the three seismograms, a difference in the slope of the SW trains can be visually identified, suggesting a temporal change in SW propagation.

Although the same equipment was used for the source-pulse generation and seismic data recording, the condition of the top surface of the embankment (dry or moist) inevitably influenced the source- and geophone-ground coupling. As a consequence, the frequency content of the source signals and recorded seismic traces may vary over time (Foti et al., 2015). Variations in the spectral amplitudes of processed seismograms, however, theoretically do not affect the temporal comparability of Rayleigh-wave phase velocity data. Phase velocity dispersion curves are related to the phase difference measured at neighboring traces of SW frequency components, rather than to their spectral amplitudes (Park et al., 1999; Xia et al., 1999; Socco and Strobbia, 2004). Geophone-ground coupling may, however, affect the reliability of SW attenuation data (see the following section for adopted countermeasures). Here again, it is not necessary to ensure a perfect reproducibility of the source signals, because attenuation curves are related to the spatial decay of SW spectral amplitudes and not to their absolute values (Rix et al., 2000, 2001; Foti et al., 2015). Attenuation data obtained using different source pulses are hence comparable. From a practical point of view, source pulses and hence seismic data acquired in time-lapse fashion should have similar frequency content distributions, to ensure comparable data quality within the frequency range of interest. We do not have direct measurements of source signals; however, tests were conducted using seismic traces recorded at the shortest offset (2 m) from the September and November 2014 acquisitions. September and November 2014 data were selected for comparison because they were representative of particularly dry and wet soil surface conditions, respectively (see the “Discussion” section). These tests (Figure 3) suggest a good reproducibility of the source pulses for single-shot acquisitions; moreover, the frequency content of the various seismic traces do not show dramatic changes over time or space with the majority of signal energy comprised within the frequency band of interest (10–60 Hz, see the following subsection).

For comparison with the geophysical data, climate sensors (manufactured by Decagon Devices Inc.) were installed on site and continuously recorded meteorological data (rainfall, humidity, solar radiation, and wind speed) over the period of interest. Point sensors (also manufactured by Decagon Devices Inc.), capable of measuring soil moisture were also installed at various depths in back-filled boreholes within the embankment flanks (Figure 1).


Rayleigh-wave dispersion curves

Stacked seismograms (Figure 2, left panels) were remapped into the frequency-wavenumber (f-k) domain by applying an f-k transform (Figure 2, right panels). Energy maxima were then picked in the f-k domain to obtain a fundamental mode SW dispersion curve for each positioning of the landstreamer (Socco and Strobbia, 2004). The analysis focused on the fundamental mode because it consistently represented the most energetic event in the f-k images, spanning approximately the same frequency band (10–60 Hz) throughout the investigated period of time (Figure 2, right panels). Therefore, the SW fundamental mode could be assumed as reference data, portraying in a reliable way the temporal and spatial variations of Rayleigh-wave propagation (in terms of phase velocity and attenuation).

For each bimonthly acquisition, a set of 21 dispersion curves spatially referenced to the midpoint of the corresponding recording array (which were regularly spaced by 5 m) were obtained. Example sets of dispersion curves acquired in the months of January, March, and September 2014 are shown in Figure 4 (left panels). The dispersion curves extend mostly in the 10–60 Hz frequency band and exhibit an approximate phase velocity range of 225 m/s (at approximately 10 Hz) to 110  m/s (at approximately 60 Hz). These values are in agreement with other SW phase velocity measurements on wet clays (Foti, 2003; Long and Donohue, 2007). Phase velocity significantly increases at approximately 17 Hz, marking the threshold between frequency components propagating only in the softer embankment (>17  Hz, roughly corresponding to wavelengths <7.25  m, exhibiting lower velocities) and those propagating also through the stiffer subsoil beneath the embankment (<17  Hz, higher velocities). Another feature in all data sets presented in Figure 4 is the gradual increase of phase velocities for frequencies >17  Hz with increasing distance from the southwest. Below 17 Hz, the opposite trend may be observed. Despite these common features, the different sets of dispersion curves also exhibit visible differences, which are analyzed in detail later in the paper (see the “Discussion” section).

Rayleigh-wave attenuation curves

The recorded seismograms were also analyzed in terms of Rayleigh-wave intrinsic attenuation. The intrinsic amplitude attenuation caused by material damping is representative of the dissipative properties of the medium (Foti et al., 2015). The Rayleigh-wave spatial amplitude decay due to intrinsic attenuation is frequency dependent and can be expressed as (Xia et al., 2002)
where Aj,i and Aj,i+1 represent the spectral amplitudes of jth frequency component at receivers i and i+1, respectively; ri and ri+1 are the offsets of the receivers; and αRj is the Rayleigh-wave attenuation coefficient for the jth frequency. From the acquired seismograms, the attenuation coefficients αRj, representing the dissipative properties of the soil, were estimated by applying a univariate regression of spectral amplitudes’ variation with offset (Foti et al., 2015). A set of attenuation curves (i.e., αR versus frequency) were hence obtained for every acquisition. Similar to the phase velocity data set, the attenuation curves are spatially referenced to the center of the corresponding recording array, i.e., at 21 locations spaced by 5 m along the seismic line. To portray the dissipative behavior of the Rayleigh-wave fundamental mode and to avoid the effects of combining modes, the attenuation coefficients were estimated in the frequency band in which the fundamental mode was dominant. Spectral amplitudes from short-offset traces were also discarded from the regression to minimize near-field effects.

The use of a landstreamer allows for rapid data acquisition, but the ground-receiver coupling is not ideal for attenuation analysis. A robust data-fitting algorithm in the regression analysis, yielding αRj coefficients, was therefore implemented to exclude anomalous spectral amplitudes from poorly coupled geophones. In Figure 4 (right panels), the sets of attenuation curves for January, March, and September 2014 are shown. The attenuation coefficient generally increases with frequency as expected (Foti et al., 2015): derived values are compatible with other Rayleigh-wave attenuation analyses conducted over wet clays (Foti, 2003). Although great care was paid in developing a robust algorithm for the reliable estimation of attenuation coefficients, the attenuation curves appear to be noisier and less smooth when compared with the dispersion curves (Figure 4). This may be attributed to the fact that the quality of attenuation data, in contrast with phase information, is considerably more dependent on the ground-receiver coupling. Due to irregularity of the attenuation curves, a definite spatial trend cannot be visually identified.

In the adopted inversion procedure, a single initial model was estimated and subsequently used for the deterministic inversion of all dispersion curve data sets (July 2013 to November 2014), following the approach introduced by Socco et al. (2009) and Boiero and Socco (2010).

Initial model

The common initial model for the parallel inversion of all dispersion curve data sets was obtained from Monte Carlo (MC) inversion of the sample dispersion curves from the March 2014 data set (Socco et al., 2009; Boiero and Socco, 2010; Konstantaki et al., 2013). March 2014 was selected as the reference data set for its central position with respect to the surveyed period of time. Three dispersion curves, each considered representative of parts of the seismic line (x<35  m, 35<x<75  m, and x>75  m), were selected and inverted. The MC inversion algorithm developed by Maraschini and Foti (2010) was used. This approach enables a computationally efficient search of the parameter space to be carried out. For each MC inversion, 106 models, with a redundant number of layers (10), were tested. Once the best fitting set of VS models were determined (Figure 5a and 5b), the initial model profile was obtained by manually reducing the number of layers over the stack of feasible models, according to a criterion of minimum parameterization (Figure 5b; Konstantaki et al., 2013). This operation was repeated for each of the three selected dispersion curves, adopting a common model parameterization. Figure 5a and 5b, for example, refer to the MC inversion of the dispersion curve located at x=55  m (from March 2014), which was chosen to represent the central (35<x<75  m) portion of the seismic line in the initial model. The initial VS model is shown in Figure 5c, and comprises three separate VS profiles (x<35  m, 35<x<75  m, and x>75  m), each one is the outcome of the MC inversion of a selected dispersion curve.

Least-squares inversion

The initial VS model (Figure 5c) was then used as input for the subsequent least-squares inversion of the nine complete Rayleigh-wave dispersion curves sets to provide a time-variant reconstruction of S-wave velocity within the surveyed embankment structure. Considering the regular and dense spatial distribution of dispersion curve data, these were inverted using a laterally constrained inversion (LCI) algorithm. The LCI, originally developed for the inversion of resistivity data (Auken and Christiansen, 2004) and later successfully applied to SW data interpretation (Wisén and Christiansen, 2006; Socco et al., 2009), is a deterministic inversion scheme involving the simultaneous inversion of spatially referenced geophysical data for a set of 1D models. These models are regularized through spatial constraints tying corresponding model parameters. This intrinsic assimilation of spatial regularization within the inversion process is particularly beneficial for SW dispersion data because the nonuniqueness of the SW method (Foti et al., 2009) is mitigated by introducing a priori information in the form of the expected spatial variability of the obtained S-wave velocity models.

In the LCI inversion of the SW dispersion curves, the selection of optimal levels of spatial regularization was of considerable importance for reconciling data compliance and spatial variability and ensuring more physically consistent seismic models (Boiero, 2009; Boiero and Socco, 2010). The selection of the level of spatial regularization was based on a set of trial inversions performed on the March 2014 reference data set, adopting an increasing strength for spatial constraints, and collating the obtained levels of data fitting (measured as normalized residuals for each experimental dispersion curve, Figure 6). Hence, we selected a form of medium constraints that provided spatial regularity within the obtained VS section, while also ensuring data compliance similar to that obtained with a condition of null constraints (Figure 6). This level of spatial smoothing was then kept constant throughout the processing of all available data sets, to ensure homogeneity of procedure and hence comparability among the obtained VS sections.

In this section, we describe in detail the results of the processing and inversion stages. We present the acquired time-lapse seismic data set comprising nine sets of Rayleigh-wave dispersion and attenuation curves. The dispersion curve data were inverted to obtain time-distributed VS images of the surveyed embankment. We illustrate the seismic model relevant to the March 2014 reference data, chosen for its central position along the time scale, and we introduce the S-wave velocity sections for the entire period of study.

Rayleigh-wave dispersion curves

Rayleigh-wave dispersion curves, exhibiting good quality and spatial consistency (see Figure 4), may be arranged into a global representation in which the spatial and temporal trend of these data can be visually assessed. Figure 7 illustrates the Rayleigh-wave dispersion data from the Laverton site, which are displayed as multiple phase velocity 2D pseudosections. These pseudosections were obtained by translating the dispersion curves from phase velocity versus frequency (Figure 4, left panels) to pseudodepth (approximated with wavelength/2; Abbiss, 1981) versus phase velocity. Phase velocity data were then interpolated over a regular grid extending along the embankment crest (considering the dispersion curves’ reference points) and pseudodepth axis to obtain a continuous image of the spatial variability of Rayleigh-wave phase velocity. A consistent color scale has been used for all panels in Figure 7 to enable the spatial and temporal variability of phase velocities to be observed. Because the embankment is the focus of this study, only the data points referring to the propagation of seismic waves through the earthwork were considered; consequently, in Figure 7, only phase velocities with a pseudodepth <5  m (height of the embankment) are displayed.

In Figure 7, spatial and temporal changes throughout the acquired data sets can be identified. All panels consistently present a common lateral pattern with increasing phase velocity values observed from southwest to northeast (see also Figure 4). In particular, a low-velocity area may be observed in all data sets approximately between 20<X<40  m. Temporal changes can also be distinguished with phase velocities generally decreasing in the period of July to November 2013. The November 2013 and January 2014 pseudosections exhibit particularly low phase velocities. In contrast, relatively high phase velocities are observed in the period from March to September 2014, with September 2014 exhibiting the highest velocity values measured over the monitoring period. A further decrease in phase velocity is observed in November 2014. Time-lapse variations of phase velocity data are analyzed in detail in the “Discussion” section.

Rayleigh-wave attenuation curves

It was not possible to effectively arrange the Rayleigh-wave attenuation data in a global representation to portray their spatial and temporal variability (as done for phase velocity in Figure 7) because the attenuation curves were of lower quality when compared with the phase velocity data (Figure 4). This was largely due to poor soil-receiver coupling, which resulted in exclusion of a significant portion of the attenuation coefficient αR data.

Hence, to focus on the temporal variability of the attenuation coefficients within the embankment, the overall distributions of αR were divided into three frequency bands (Figure 8, left panels), 37.5–42.5, 17.5–22.5, and 13–16 Hz, corresponding (according to the approximate half-wavelength criterion; Abbiss, 1981) to 1.5, 3, and 4.5 m depth. The panels in the left column of Figure 8 present the attenuation coefficient distributions for all acquired data sets, belonging to the three adopted frequency bands. Similar to the analysis of Rayleigh-wave phase velocities (Figure 7), it is possible to track a comprehensive and reasonably consistent temporal pattern of variation for attenuation data. In chronological order, values of αR generally increase from July to November 2013, attenuation coefficients between March and September 2014 reach lower values again, and finally they increase from September to November 2014. This seasonal trend is more evident and clear in the 17.5–22.5 and 37.5–42.5 Hz bands (the left panels of Figure 8a and 8b). As for the 13–16 Hz range (the left panel of Figure 8c), low values of αR in January and November 2014 do not agree with this trend.

Figure 8 also compares the temporal trend of αR with the corresponding phase velocities from Rayleigh-wave dispersion curves (right panels of Figure 8), relevant to the same three frequency bands. Phase velocities exhibit an opposite trend to that of the attenuation coefficients; i.e., an increase in αR corresponds in general to a decrease in velocity, and vice versa. The pattern of temporal variability for phase velocities is of course the same as that observed in Figure 7.

Inversion of Rayleigh-wave dispersion curves

The 2D VS pseudosection obtained from the LCI of the March 2014 dispersion curves (Figure 4b, left panel) is illustrated in Figure 9a. The pseudosection in Figure 9a is obtained from lateral interpolation of the 21 reconstructed vertical VS profiles, each of which refers to a particular dispersion curve. Figure 9c reports the data fitting, in terms of normalized residuals, between each experimental dispersion curve and the corresponding synthetic curve from the inverted VS profile. As shown, normalized residuals are relatively low and fall within the range 0.1–0.4. Some examples of experimental/synthetic curves are reported in Figure 10, which illustrate good experimental data compliance.

As far as the March 2014 VS model is concerned (Figure 9a), a discontinuity in the S-wave velocity distribution at an approximate 5 m depth is identified (i.e., the base of the embankment, which is consistent with the CPT data, Figure 1). Within the embankment VS increases with depth from 110 to 160  m/s. These values are consistent with seismic CPT (SCPT) data collected at various locations along the seismic line (Figure 9b), although at different stages of the year (January 2013 and July 2014). The increase of S-wave velocity with depth is not gradual and smooth. The upper 2 m are characterized by low VS (generally <125  m/s). Below this soft layer, S-wave velocities increase sharply toward 140–160 m/s. In particular, a relatively high velocity stiff layer at approximately 3 m depth may be identified. At some locations in the embankment, this layer exhibits greater VS than the underlying formation (between 45<x<70  m in Figure 9a), again, in agreement with the SCPT profiles. This distribution of VS with depth is thought to relate to the building technique adopted for the embankment, probably involving two successive construction phases, with an intermediate compaction stage once half of the total embankment height had been reached (Gunn et al., 2015b). In terms of lateral VS changes, the embankment and the underlying material exhibit differing velocity trends: proceeding from southwest to northeast, VS gently increases within the body of the embankment, whereas it decreases in the underlying subsoil. This is coherent with the spatial trend observed in the dispersion curve data sets (Figure 4, left panels).

These features of the 2D VS model from March 2014 (Figure 9a) are common to the S-wave velocity pseudosections obtained from inversion of the dispersion curves extracted throughout the monitoring period (Figure 11). Here again, as in Figure 7, spatial and color boundaries are adapted to focus on the temporal variations of VS within the embankment structure (depth <5  m).

The temporal variability of the acquired data sets (dispersion and attenuation curves) and the inverted VS models were investigated to quantitatively define a seasonal trend that can be correlated with the precipitation pattern and corresponding temporal variability of the volumetric water content within the embankment.

Rayleigh-wave dispersion curves

Because an identical acquisition geometry was used for each data acquisition, and given the high-quality Rayleigh-wave dispersion curves, a point by point comparison method was adopted. Therefore, in Figure 12, each phase velocity is compared with its equivalent from the March 2014 data set, which is considered as a reference data set for its central position on the time scale. As with Figure 7, only data referring to seismic waves propagating through the earthwork were considered (i.e., phase velocities with wavelength/2 <5  m). In Figure 12, the temporal variations of phase velocity observed in the previous section are confirmed and quantified (i.e., higher velocity values in spring-summer months as opposed to lower values in autumn-winter). Changes in phase velocity are mainly confined within a ±10% range. In terms of spatial distribution, it may be observed that the northeast portion of the embankment (approximately x>40  m) is generally characterized by greater variability in Rayleigh-wave phase velocity. In the northeastern half of the investigated area, two zones of relatively homogeneous temporal variations located approximately below and above a pseudodepth of 3.5 m may also be distinguished. Above 3.5 m (pseudodepth), phase velocity appears to vary more significantly and changes occur more rapidly (compare, for instance, the September to November 2014 panels), whereas a smoother variation over time is observed at pseudodepths >3.5  m. Intuitively, this can be related to the greater and more immediate susceptibility of the shallower portion of the embankment toward climate-related changes of phase velocity.

In Figure 13, the temporal variability in the Rayleigh-wave phase velocities is summarized and compared with the recorded precipitation pattern and volumetric water content measured by sensors installed within the embankment flanks. Each row in Figures 13a represents a histogram of the point-by-point changes in phase velocity relative to March 2014. Figure 13a shows a continuous and consistent seasonal trend that can be correlated with the rainfall data (Figure 13b). In Figure 13a, comparing July 2013 with the following autumn-winter acquisitions (September 2013 to January 2014), phase velocity globally decreases, which corresponds with frequent and intense autumn-winter rainfall (Figure 13b). This is in contrast with the spring-summer 2014 acquisitions (May 2014 to September 2014) in which a general increase of phase velocities and a decrease of precipitations may be observed, when compared with previous data sets. The “fastest” data set is from September 2014, acquired after the driest recorded period of time, with only 3.4 mm of rainfall in three weeks. The November 2014 data set, acquired after the start of intense autumn rainfall, shows much slower velocities, globally similar to those of autumn-winter 2013. The correlation between the seasonal trend of phase velocities and the precipitation pattern is consistent with the expected variation of VS with soil suction as a result of climate-induced soil moisture/saturation changes within the embankment clay fill (Cho and Santamarina, 2001; Santamarina et al., 2005). Intense and frequent rainfall will infiltrate the embankment, leading to an increase in soil water content/saturation. This increase in saturation will cause a corresponding rise in pore water pressure (i.e., reduction in suction pressure), leading to a reduction in VS.

Rayleigh-wave phase velocities are related to VS through Poisson’s ratio (Xia et al., 1999; 2003; Stokoe et al., 2004) and hence would be expected to exhibit a similar relationship with soil saturation. Consequently, the autumn-winter period, with more severe and frequent rainfall (Figure 13b), is characterized by lower Rayleigh-wave phase velocities. Oppositely, the spring-summer months are characterized by alternating periods of very low or null precipitation and extreme rainfall events (Figure 13b). These events generally saturate the shallower portion of the embankment fill (dried by preceding periods of low rainfall) but do not penetrate deeper into the embankment. Hence, the spring-summer months are likely to be characterized by lower values of water content that result in higher Rayleigh-wave phase velocities (Figure 13a).

In addition to precipitation, other factors may alter the water content within the embankment structure. The significant increase of air temperature (Figure 13b) in the spring/summer period will cause an increase in evaporation of aqueous phase from the soil to the atmosphere. Moreover, the greater temporal variability of velocity values at x>40  m (Figures 7 and 12), can be ascribed to the presence of thick vegetation on the flanks of this particular stretch of embankment, implying the abstraction of water from the soil by the roots of plants.

The volumetric water content values recorded at the dates of seismic data acquisition (Figure 13c and 13d) cannot be closely related to the observed temporal changes of seismic data because they were measured by sensors not directly located in the area covered by seismic surveys (i.e., on the flanks of the embankment, not below its crest; see Figure 1), and moreover they span only a portion of the investigation time period. Still, a significant temporal variability in water content within the embankment structure may be observed (from values as high as 0.55 to as low as 0.05). As expected, this variability correlates with the observed precipitation pattern (in particular, note the general decrease of water content values toward September 2014, and their considerable increase in November 2014, Figure 13c and 13d). The substantial variation of water content is compatible with the corresponding temporal changes in the seismic data.

Temporal variability versus experimental uncertainty

For Rayleigh-wave phase velocities, we verified the relevance of the observed variations over time with respect to the data experimental uncertainties. As suggested by Socco et al. (2009), we obtained the uncertainty for Rayleigh-wave dispersion curves from the stacked seismograms (Figure 4) by extracting the curves from the three relevant individual shots (see the “Data acquisition and processing” section): the velocity range identified for each frequency by these dispersion curves defines the uncertainty range of phase velocity. For the purpose of this test, only the experimental uncertainties for the September and November 2014 data sets were computed. The temporal variability between September and November 2014 was then compared with the corresponding experimental uncertainties, both individually, to ascertain whether temporal variability exceeded the uncertainty range, and globally, to compare their statistical distribution. As previously, only the data points corresponding to Rayleigh-wave propagation through the embankment were used (see Figures 7 and 12). The majority (84%) of the observed point-by-point data variations between September and November 2014 were found to be greater than the uncertainty range. In addition, the distributions of the experimental uncertainties for both data sets were compared with the point-by-point temporal variations (Figure 14). The experimental uncertainties show similar distributions for the September and November data sets: uncertainties appear to be symmetrically distributed around 0% and they lie within a relatively narrow interval (90% of uncertainties being comprised within ±1% range). The similarity between the September and November 2014 uncertainty distributions suggests that the level of uncertainty may be approximately constant throughout all data sets. Conversely, the distribution of temporal changes in the data (gray line in Figure 14) appears to be radically different, significantly skewed toward negative values (90% of phase velocity differences are in the range 12% to 0%). This indicates a global reduction in Rayleigh-wave phase velocity that is significantly outside the experimental uncertainty range; hence, the temporal variability in phase velocity may be reliably assumed to be related to an actual change in the seismic properties of the embankment.

Rayleigh-wave attenuation curves

The temporal pattern of Rayleigh-wave attenuation data was also analyzed and correlated with rainfall data. As discussed in the “Results” section and illustrated in Figure 8, αR from the three selected frequency bands follows similar temporal behavior, showing a general increase of αR from July 2013 to the successive autumn-winter months (September 2013 to January 2014); αR then declines to lower values for the spring-summer period (March to September 2014) and rises once more in November 2014. Again, this temporal pattern corresponds well with the observed rainfall distribution (Figure 13b), involving frequent and intense rains in autumn-winter and alternating periods of low rainfall and extreme precipitation events in spring-summer. As with the phase velocity data, the seasonal change in attenuation values is therefore likely to be related to the variation in the degree of saturation of the embankment clay fill. An increase of water content within the railway embankment will cause a corresponding rise in pore water pressure (i.e., reduction in suction pressure), and hence lead to an increase of Rayleigh-wave intrinsic attenuation: oppositely, drier conditions lead to an increase of effective stress and therefore to a decrease of intrinsic attenuation (Santamarina et al., 2001).

The definite temporal trends of αR (Figure 8, left column) and phase velocities (Figure 8, right column) thus confirm the sensitivity of Rayleigh-wave data, in terms of attenuation and velocity, in portraying climate-related changes in the mechanical properties of the subsoil. It should be pointed out that under the condition of VP/VS>2 (which is generally verified for wet clays), the influence of P-wave on Rayleigh-wave attenuation is negligible (Xia et al., 2002; Foti, 2004). The seasonal trend in αR may therefore be identified with the S-wave attenuation behavior.

Inverted S-wave velocity models

Detailed analysis of the inversion results obtained for the March 2014 reference data set and their agreement with available geologic and invasive survey (SCPT and CPT) data confirm the validity of the chosen inversion strategy (see the “Results” section). The inverted VS models for all data sets (July 2013 to November 2014) were then used for monitoring the seasonal seismic property variations within the embankment.

As expected, the reconstructed VS models (Figure 11) show spatial and temporal trends similar to those previously identified for Rayleigh-wave phase velocities. The main trends and features that may be observed are as follows:

  1. lateral increase of velocities in the embankment body (depths <5  m), proceeding from southwest to northeast along the seismic line

  2. a low-velocity area located at approximately 20<x<40  m

  3. general decrease of VS in the embankment body (depths <5  m) from July 2013 to January 2014, succeeded by a rise in spring-summer (March to September 2014), followed by a further decrease in November 2014

  4. greater temporal variability of VS for x>40  m

As explained in previous subsections, the variability over time of seismic velocities in fine-grained soils can be linked to the seasonal change of water content (Santamarina et al., 2001), which can in turn be explained with seasonal and climate-related phenomena, such as the distribution of precipitation over the yearly cycle and evapotranspiration.

The VS models of the embankment structure resulting from the inversion process globally display a generally consistent trend over time. Figure 15 represents the relative variations of VS with respect to their temporal average in the three shallower layers (1.5 m thick each), which model the embankment structure (depths between 0.75 and 5.25 m). The upper 0.75 m was excluded from the analysis because of insufficient Rayleigh-wave sampling (the minimum detectable wavelength being equal to the intergeophone distance of 1 m; Socco and Strobbia, 2004). Quantitatively speaking, VS values appear to vary within a range of approximately ±10%. The variability of inverted VS is consistent with that observed in VS values from SCPT data (Figure 9b), considered representative of winter (January 2013 acquisition) and summer (July 2014 acquisition) conditions. In terms of temporal variations in VS, a global trend may be identified, involving lower values in autumn-winter and higher values in spring-summer; this pattern is clear for the shallow layer (Figure 15a), but it is less evident in the deepest (Figure 15c) and particularly intermediate (Figure 15b) formations, where lateral variations in the temporal trend of VS are present. This clear temporal pattern of the upper layer can be ascribed to the more direct exposure of this layer to weather-related effects (precipitation and evapotranspiration), whereas in the deeper formations additional phenomena (e.g., a lateral variability of material properties within the embankment structure) may affect the seasonal change of moisture and hence S-wave velocity. Other considerations, linked to inherent features of SW inversion rather than to geomechanical aspects, could also explain the smoother VS pattern for the upper layer (Figure 15a) as opposed to the two deeper layers (Figure 15b and 15c). Rayleigh-wave dispersion data sample the shallowest layer more densely when compared with the deeper layers (the number of frequency components sampling a given depth decreases as depth increases; Socco and Strobbia, 2004). In addition, equivalence problems (Tarantola, 2005; Foti et al., 2009) might affect the estimation of VS for the deeper layers.

Although a careful inversion strategy such as LCI was adopted, the use of SW data, on its own, does not appear to be fully adequate for reconstructing climate-related changes in VS. To mitigate the inherent ambiguities of the SW method, a possible solution could be offered by performing laboratory calibration tests, with the aim of establishing a defined relationship between the measured seismic velocities and relevant mechanical properties of soil samples. This relationship could then be used for a detailed interpretation of the field geophysical data. In addition, joint inversion (Boiero and Socco, 2014) of different seismic data sets (such as SW and P-wave refraction data) would potentially provide a single, physically and spatially consistent seismic model. The inclusion of electrical resistivity measurements in the inversion scheme (Garofalo et al., 2015) would also likely be beneficial and would offer the possibility of linking a method inherently sensitive to temporal changes of moisture content (electrical resistivity tomography) to geophysical techniques more closely related to the mechanical behavior of the subsoil (seismic methods).

SW data were collected repeatedly along the crest of a heritage railway embankment with the aim of portraying the seasonal variations of S-wave velocity within the body of the earthwork. Acquisitions were repeated bimonthly in the period July 2013 to November 2014. Great care was given to ensuring that the same acquisition geometry and parameters were adopted. Seismic equipment, capable of swiftly acquiring seismic data and hence efficiently covering relatively long profiles, was used (sledgehammer seismic source, vehicle towed landstreamer). The acquired seismograms were processed to extract Rayleigh-wave phase velocity and attenuation data. The Rayleigh-wave dispersion curve data sets were subsequently inverted with the application of an LCI scheme with the aim of obtaining a time-variant S-wave velocity model of the embankment structure.

The time-lapse SW data sets were analyzed to identify temporal variations and to reconstruct their global trend over time. The good quality and spatial consistency of the Rayleigh-wave dispersion curves allowed a detailed “data point by data point” comparison method to be adopted. We determined consistent spatial and temporal variations for these data, observing an overall increase of Rayleigh velocities in the spring-summer months, and vice versa lower velocity values during the autumn-winter period. This temporal trend correlates well with the precipitation pattern and soil water content measured at the site, and it is consistent with the expected variation of seismic velocities with soil saturation and the corresponding pore water pressure changes. It was also verified that the observed temporal variability of phase velocity generally exceeds the range of experimental uncertainties.

Rayleigh-wave attenuation curves appear to be characterized by lesser smoothness and do not display a clear spatial pattern. Analyzing their temporal variability is therefore more challenging; however, a temporal trend in the overall distribution of attenuation coefficients is apparent (lower attenuation coefficients in spring-summer and higher values in autumn-winter), which is consistent with that of Rayleigh-wave phase velocities and with the seasonal cycle of water content.

It has therefore been shown that with well-established SW methods for the characterization of the near surface (phase velocity and attenuation analysis), it is possible to obtain reliable time-lapse data sets, capable of portraying seasonal variations of seismic velocities. The results obtained in terms of acquisition and processing of SW data are particularly relevant from an engineering point of view because they demonstrate the potential of this noninvasive, time-effective, and cost-effective approach for time-lapse monitoring of the geomechanical properties of earthworks at risk of failure.

Time-lapse variations of phase velocity data coherently influenced the VS models obtained from the successive inversion stage. Although consistent temporal variations of VS within the embankment structure were reconstructed, the level of detail was not sufficient to allow a geomechanical interpretation of such variations. This was ascribed to inherent features of SW data inversion. Further developments of the present work are hence related to the achievement of a greater consistency of the inversion and calibration of the field data through a series of laboratory tests.

This work formed part of the project GEOphysical Condition Assessment of Railway Earthworks (GEOCARE—PI S. Donohue) which was funded by the Engineering and Physical Sciences Research Council (EPSRC). The contributions to this work from B. Dashwood, S. Uhlemann, R. Swift, J. E. Chambers, and D. A. Gunn are published with permission from the Executive Director of The British Geological Survey (NERC). The authors acknowledge the support and the help in providing access to the site of the Gloucestershire-Warwickshire Steam Railway (GWR) staff. The authors would like to thank S. Foti, R. Cosentini, and V. Socco from Politecnico di Torino for their valuable discussions and suggestions.

Freely available online through the SEG open-access option.