Time-lapse monitoring of fluid induced geophysical property variations within an unstable earthwork using P-wave refraction

A significant portion of the UK ’ s transportation system relies on a network of geotechnical earthworks (cuttings and embankments) that were constructed more than 100 years ago, whose stability is affected by the change in precipitation patterns experienced over the past few decades. The vulnerability of these structures requires a reliable, cost- and time-effective monitoring of their geomechanical condition. We have assessed the potential application of P-wave refraction for tracking the seasonal variations of seismic properties within an aged clay-filled railway embankment, located in southwest England. Seismic data were acquired repeatedly along the crest of the earthwork at regular time intervals, for a total period of 16 months. P-wave first-break times were picked from all available recorded traces, to obtain a set of hodocrones referenced to the same spatial locations, for various dates along the surveyed period of time. Traveltimes extracted from each acquisition were then compared to track the pattern of their temporal variability. The relevance of such variations over time was compared with the data experimental uncertainty. The multiple set of hodocrones was subsequently inverted using a tomographic approach, to retrieve a time-lapse model of V P for the embankment structure. To directly compare the reconstructed V P sections, identical initial models and spatial regularization were used for the inversion of all available data sets. A consistent temporal trend for P-wave traveltimes, and consequently for the reconstructed V P models, was identified. This pattern could be related to the seasonal distribution of precipitation and soil-water content measured on site.

Time-lapse monitoring of fluid-induced geophysical property variations within an unstable earthwork using P-wave refraction INTRODUCTION A significant portion of the UK's transportation system relies on a network of geotechnical earthworks (cuttings/embankments) that are more than 100 years old and were not built to modern construction standards. The present condition of these earthworks has been further compromised by poor maintenance strategies and by climate variations experienced over recent decades (Jenkins et al., 2009). Wetter winters and drier summers may cause seasonal variations of pore-water pressure, which cyclically deteriorate the stability of earthworks (Kovacevic et al., 2001). Additionally, more frequent severe weather events can trigger earthwork slope failures (Wilks, 2010). In 2012 (the second-wettest year on record in the UK), Network Rail, who operate the majority of the UK rail network, reported 144 earthwork failuressignificantly more than in previous years. The vast majority of these failures were associated with heavy rainfall (Network Rail, 2013).
The vulnerability of these earthwork structures requires reliable, time-and cost-effective monitoring of their geomechanical condition. Noninvasive geophysical methods, which are capable of estimating physical properties of significant volumes of the subsurface, are well-suited for this application (Donohue et al., 2011). Electrical resistivity tomography (ERT) (Loke and Barker, 1996), for example, has, over recent years, been used for near-surface time-lapse imaging of natural slopes and landslides (e.g., Miller et al., 2008;Uhlenbrook et al., 2008;Cassiani et al., 2009)  the topic of the present work, monitoring the internal condition and temporal dynamics of geotechnical assets (e.g., Sjodahl et al., 2008, Wilkinson et al., 2011Chambers et al., 2014;Gunn et al., 2015a). Favored by the sensitivity of resistivity data to changing hydrogeological conditions, ERT, however, requires an additional stage of calibration, often site dependent, to invert for parameters more relevant to the mechanical description of the investigated target (e.g., moisture content, as in Cassiani et al., 2009;Chambers et al., 2014). Seismic geophysical methods, which enable seismic velocity models of the subsurface to be reconstructed, are more directly related to the mechanical properties of the subsurface. In this work, we particularly focus on the P-wave refraction technique (Redpath, 1973) that is regarded as an established tool for the reconstruction of subsurface V P models (ASTM International, 2011). It is based on the identification of P-wave first-arrival times and on the subsequent tracing of corresponding rays to derive a V P model of the subsurface. The P-wave refraction method has a wide range of uses at several scales of application, from seismological studies (e.g., Zelt and Smith, 1992) to the nondestructive evaluation of cultural heritage artefacts (e.g., Sambuelli et al., 2011). At the engineering scale, P-wave refraction is frequently used for the reconstruction of seismic models for the shallow layers in hydrocarbon exploration geophysics (e.g., Macrides and Dennis, 1994), imaging of complex structures in the subsurface (Nyquist et al., 1996;Ramos Martinez et al., 1997), monitoring of environmental issues (Lanz et al., 1998;de Iaco et al., 2003), landslide monitoring (e.g., Uhlemann et al., 2015), hydrogeological studies (e.g., Ayers, 1990;Donohue et al., 2015), and detection of underground cavities (Sloan et al., 2015;Toth et al., 2015).
Despite their widespread application, the use of seismic techniques for time-lapse monitoring of geomechanical properties has, however, thus far been mainly restricted to large-scale applications, e.g., hydrocarbon reservoir exploration geophysics (among others, Landro, 2001;Lumley et al., 2003) or the monitoring of geologic storage reservoirs (e.g., Draganov et al., 2012;Picotti et al., 2012). In near-surface geophysics, the potential use of seismic techniques for tracking time-variant processes has been assessed (e.g., Jefferson et al., 1998;West and Menke, 2000), although at reduced spatial scales and spanning shorter time intervals (from hours to some days), when compared with those involved in this paper. Recently, the use of seismic techniques has been proposed for monitoring water table levels in shallow aquifers (Pasquet, 2014;Pasquet et al., 2015aPasquet et al., , 2015b. In this paper, seismic techniques are not used for detecting the condition of full saturation, but they are instead used for evaluating the fluctuating seismic properties of a railway earthwork in response to varying levels of water content, usually below saturation. The seismic behavior of compacted soils used in earthworks is influenced by changing levels of water content, affecting capillary/suction pressures, which in turn influence the state of stress of the solid matrix of the soil (Mancuso et al., 2002). Many authors have highlighted the sensitivity of seismic velocities to changes in soil suction (e.g., Qian et al., 1991;Marinho et al., 1995;Cho and Santamarina, 2001;Mancuso et al., 2002;Donohue and Long, 2010); however, these studies are limited to laboratory measurements and do not comprise field surveys.
This study proposes and tests the use of P-wave refraction to track these variable properties within the structure of a heritage railway embankment; it complements the work by Bergamo et al. (2016), focused on surface-wave data and methodology. In the following sections, we describe the surveyed site and the periodical acquisition of seismic data (for a total period of 16 months) and their processing in terms of extraction of multiple sets of P-wave first-break arrivals. The trend of seasonal variability of the obtained traveltimes is reconstructed, and its relevance with respect to data experimental uncertainty is also investigated. The multiple sets of P-wave traveltime data are then inverted, using a tomographic approach, to obtain a time-variant model of V P distribution within the embankment structure.

SITE DESCRIPTION
The geotechnical structure investigated in this study is a stretch of railway embankment of the heritage Gloucestershire-Warwickshire Railway. The site is located near the village of Laverton, Gloucestershire, in southwest England. The embankment was constructed in the early 20th century (approximately 1900-1906); the investigated portion is approximately 5 m high and approximately 9 m wide at the top (approximately 30 m at the base), and it runs from southwest to northeast (Figure 1). The structure comprises approximately 4 m of compacted high plasticity clay fill (locally derived Charmouth Mudstone), and it is topped by approximately 0.9 m of humus-rich ballast fouled with fines, ash, and soil (Gunn et al., 2015b). In previous years, several significant earthwork failures have occurred along the Gloucestershire-Warwickshire Railway. Close to the Laverton site, there is visual evidence of previous movement along the earth-  Figure 1 displays a stratigraphic log obtained from intrusive investigation involving a program of boring, sampling, trial pitting, and cone penetration testing (CPT). A sample CPT profile from the site is also presented (upper left panel in Figure 1), obtained from a location (CPT3 in Figure 1) within the area covered by the seismic surveys. The behavior of cone end resistance (q c ) with depth provides evidence of two layers within the earthwork: an upper layer (depths 1-3 m) characterized by lower q c values (0.5-1 MPa) overlies a deeper portion of fill exhibiting slightly greater q c (1-1.3 MPa). This suggests that the embankment may have been constructed in two lifts, with the boundary at approximately 3 m depth. Below the embankment base (approximately 5 m), q c increases significantly, as the cone penetrates into the stiffer formation underlying the embankment.

DATA ACQUISITION AND PROCESSING
P-wave refraction data (Redpath, 1973) were recorded bimonthly, from July 2013 to November 2014, along the crest of the railway embankment, for a total number of nine data sets. For all acquisitions, the same source type, receiver array, recording equipment, and acquisition geometry/parameters were adopted, to make the acquired seismic data fully comparable. To ensure reproducibility of the acquisition geometry, permanent marks were placed for spatial reference along the surveyed area, indicating the position of sources and receivers.
The configuration of the recording array adopted in all nine acquisition campaigns consisted of 24, 4.5 Hz, spike coupled, vertical geophones, positioned along a straight line on the crest of the embankment ( Figure 1) at 2 m intervals. The seismic source was a 4.5 kg sledgehammer striking a metal plate. The source was positioned at 16 different locations along the seismic line. For each of these source positions, three recordings were acquired and then stacked in the time domain, to obtain seismic sections with greater signal-to-noise ratio (Socco and Strobbia, 2004). A sampling rate of 62.5 μs and a 1 s record length was used. P-wave first-break arrival times were manually picked on stacked seismograms. The upper panels in Fig black crosses indicate the manually picked P-wave first breaks. Parallel to P-wave refraction data acquisition (i.e., at the same dates and spatial location; Figure 1), multichannel analysis of surface waves (MASW) surveys were also conducted; these surveys and their results are described by Bergamo et al. (2016). The MASW data were acquired using a landstreamer, with 24 plate coupled geophones spaced at 1 m. During each acquisition campaign, the geophone array was successively moved by 5 m intervals, to cover a 100 m stretch of the railway embankment, overlapping in its central section with the refraction line ( Figure 1). The seismic source, geophone type (4.5 Hz vertical receiver), sampling rate, and number of stacks were the same adopted for refraction data acquisition. To in-crease data coverage, P-wave first-break arrival times were also picked from the 10 MASW seismic sections overlapping with the P-wave refraction line. Despite the varying conditions of the soil surface (affecting source-and receiver-ground coupling), tests conducted on MASW data demonstrated a good reproducibility of the source signal and its frequency content, without dramatic changes over space (along the acquisition line) and time (among different acquisition campaigns; Bergamo et al., 2016). More importantly, to ensure a fair comparison of traveltimes extracted at different times during the surveyed period, a reliable, automatic triggering system was adopted, which was based on a starter sensor placed on the hammer head. This efficient triggering system was also essential for the operation of stacking multiple single-shot seismograms in the time domain, ensuring synchronous, and hence, superimposable acquisitions (Foti et al., 2015). The optimal prestack synchronization of recorded seismograms was further improved using a control algorithm based on the cross-correlation of the seismic traces recorded closest to the source position.
The output data set obtained from each acquisition campaign after the processing (picking) stage was a set of P-wave first-break traveltimes, referenced to the same spatial locations. An example set of picked traveltimes for the January, March, and September 2014 acquisitions are shown in Figure 2d-2f, in the time-versus-offset domain (this data configuration is henceforth referred to as hodocrone). A clear change in slope at an offset of approximately 16 m can be observed. This value coincides with the approximate boundary, identified in the subsequent inversion stage (see the "Inversion method" section), between the seismic rays propagating exclusively through the embankment (offset <16 m, steeper slopes in Figure 2d-2f indicating lower V P ), and the rays penetrating also into the soil underlying the earthwork (offset >16 m). Despite this common feature, differences in terms of arrival times (and hence To assess the reliability of the manual picking procedure and to evaluate the level of experimental uncertainty affecting the hodocrone data sets, the ranges of uncertainty for the traveltimes from two sample data sets (September 2014 and November 2014) were computed. Following a procedure similar to that adopted by Bergamo et al. (2016) for surface-wave data, the uncertainty ranges for the arrival times from stacked seismograms (Figure 2a-2c) were estimated by additionally picking P-wave first breaks on the relevant three individual shots. The time interval determined by the picks from the individual shot seismograms defines the uncertainty range for the corresponding traveltime identified on the stacked seismic section. Figure 3 displays the uncertainties obtained for the September and November 2014 P-wave hodocrone data sets, represented as function of source-receiver offset. The two data sets exhibit similar behavior: uncertainties appear to be evenly distributed around 0% and lie within a relatively narrow interval. The similarity between September and November 2014 suggests that the level of uncertainty may be approximately constant throughout all data sets.
In addition to the seismic surveys, climate data were also recorded, using permanent devices installed on site, which were capable of continuously measuring rainfall, humidity, solar radiation, and wind speed for the period of interest. Several point sensors were installed at various depths in back-filled boreholes on the embankment flanks, to provide a continuous measurement of volumetric water content ( Figure 1). All of these sensors were manufactured by Decagon Devices Inc.

INVERSION METHOD
Multiple sets of P-wave hodocrones, acquired in time-lapse fashion over the crest of the surveyed railway embankment, were inverted to retrieve a V P model of the investigated area, and to reconstruct seasonal changes of P-wave velocity within the earthwork structure.

Initial model
The initial model for P-wave traveltime inversion was derived from the 1D, intercept time interpretation (Reynolds, 2011) of the March 2014 hodocrones (Figure 2e), selected as reference data set for its central position with respect to the surveyed period of time. The obtained 1D V P model comprises a shallow layer (0 −4.5 m depth) with V P ¼ 420 m∕s, an intermediate layer (4.5-9 m depth) with V P ¼ 1190 m∕s, and overlying a half-space with V P ¼ 1800 m∕s. This 1D model was then adapted to the discretization of the surveyed section of embankment into a set of regular cells (2 × 1.5 m each). The layered structure of the 1D V P model was hence modified into a stack of thinner blocks (1.5 m thick), whose P-wave velocities increased from 420 (shallower cells) to 1800 m∕s (cells at depths ≥ 10 m; Figure 4).
The discretization of the investigated section of earthwork as an array of 29 × 9 rectangular cells ( Figure 4) is functional to the following stage of tomographic inversion. The design of such discretization reconciles several needs: (1) covering the whole subsurface section surveyed by seismic rays with a sufficient degree of spatial resolution, (2) ensuring a well-conditioned, predominantly overdetermined inversion problem (the total number of traveltimes for every data set, approximately 515, exceeding the total number of unknowns), and (3) providing regular-sized cells with a height and width of comparable size (GeoTomCG, 2004;Jones 2010). Additionally, the vertical discretization of the embankment structure is compatible with the one adopted for surface-wave data inversion (Bergamo et al., 2016), so that the inversion results of surface-wave dispersion curves and P-wave traveltimes can be adequately compared (see the "Discussion" section).

Tomographic inversion
The adopted inversion scheme involved a tomographic interpretation of P-wave traveltimes, using the commercial software for tomographic analyses GeoTomCG®. The code performs inversions with the simultaneous iterative reconstruction technique (Lytle et al., 1978;Peterson et al., 1985), modifying a conveniently discretized initial velocity model through an iterative process. Each iteration comprises three successive steps: (1) computation of model traveltimes, (2) calculation of residuals, and (3) application of velocity corrections to produce an updated velocity model. A curved raytracing method with a revised form of ray bending is adopted  for forward modeling (GeoTomCG, 2004), derived from the Um and Thurber (1987) technique.
Considering the limited extent of the surveyed section of the embankment (approximately 50 m) and the available geologic information on its structure (a relatively homogeneous construction technique implying reduced lateral variations), a form of lateral regularization was adopted in the inversion. To choose the optimal level of spatial smoothing, a series of trial inversions of the reference data set (March 2014) were carried out, with different degrees of lateral smoothing. By comparing the distributions of normalized residuals from the various inversions ( Figure 5), a medium level of spatial regularization was identified as optimal, leading to realistically homogeneous V P sections without significantly penalizing the data fitting (Boiero, 2009;Boiero and Socco, 2010). Further increases in the degree of spatial smoothing caused an appreciable increase of misfit ( Figure 5). The selected level of regularization was used for the inversion of all nine P-wave hodocrone data sets.

RESULTS
The V P section obtained from tomographic inversion of the March 2014 data set (Figure 2e) is displayed in Figure 6a. The V P increases smoothly with depth (from approximately 400 to 650 m∕s) in the shallower portion of the embankment (depths 0-3.75 m), and higher V P values (approximately 850 m∕s) are observed in the lower layer of the earthwork structure (depths 3.75-5.25 m). This range of velocities (400 − 850 m∕s) is compatible with the expected values of V P in unsaturated clays (Knight and Endres, 2005). Below the base of the embankment (depths >5.25 m), P-wave velocity sharply increases to approximately 1100 m∕s, then reaching 1700 m∕s at 10 m depth. Both characteristics, i.e., a stiffer lower portion of the embankment and a sharp discontinuity at the base of the earthwork, are consistent with the results of cone penetration testing (CPT) (Figure 1). In terms of lateral variations, an area with slightly lower V P is located at approximately 30 < x < 40 m. Figure 6b and 6c displays the corresponding normalized residuals, obtained from the comparison between simulated and experimental traveltimes. Both panels illustrate a satisfactory experimental data compliance, with 78% of residuals falling below the threshold of 5% (Figure 6c). Similar features in the distribution of V P are present in most of the sections obtained from inversion of all other traveltime data sets (Figure 7). In Figure 7, the color and spatial scale were chosen to focus on the temporal changes of Pwave velocity within the body of the embankment (depths <5 m), which is the target of this study. The reconstructed variations over time of V P are analyzed in detail in the "Discussion" section.

DISCUSSION
The temporal variability of the P-wave traveltime data sets was accurately investigated to quantitatively define a consistent seasonal trend. The variations of P-wave hodocrone data sets coherently influenced the variability of the V P sections obtained from the inversion process. The pattern of variation of P-wave traveltimes and reconstructed velocities was also collated with the results of the analysis and inversion of the timelapse surface-wave data set acquired in parallel to P-wave refraction data (Bergamo et al., 2016).

P-wave traveltime data sets
As the hodocrone data sets were of good quality and spatial consistency (Figure 2d-2f), they could be arranged into a global representation, where their overall spatial and temporal trend could be visually assessed. All of the acquired traveltime data are displayed in Figure 8 as 2D sections referred to the same spatial coordinates (source and receiver position along the seismic line). Because this study focuses on the temporal variation of seismic velocities within the embankment structure, only the data points referring to the propagation of P-waves through the earthwork are displayed (offset <16 m; see the "Data acquisition and processing" section). A common feature in the spatial distribution of traveltimes is an area of higher values of arrival times approximately located at 30 < x < 40 m (narrower blue band at short offsets), particularly evident   in September 2013, January, March, September, and November 2014. This feature indicates an area with lower P-wave velocity that is in agreement with observations from the analysis of surface-wave data (Bergamo et al., 2016). As far as temporal variations are concerned, P-wave traveltimes appear to increase (indicating a corresponding decrease in P-wave velocities) from July to November 2013, and remain relatively high until January 2014 (thinner blue area at short offsets in Figure 8). Oppositely, in March-September 2014, a decrease of traveltimes with respect to the preceding winter period (wider blue area at short offsets) is observed, suggesting an increase of V P within the embankment structure. Successively (November 2014), the arrival times increase again, representing a decrease of P-wave velocities. Hence, traveltimes appear to follow a seasonal variation pattern.
Because seismic data were collected using an identical acquisition geometry, to quantitatively evaluate the temporal data variations, a "point-by-point" comparison method was adopted; i.e., every traveltime was compared with its corresponding data point from the reference data set (March 2014). The results of this "data point by data point" comparison are illustrated in Figure 9. This confirms and quantifies the observations regarding the temporal variations of traveltime data identified in Figure 8. Relative changes of arrival times are mainly within a AE20% variability range. The relatively short length of the P-wave refraction line limits the possibility of identifying clear spatial patterns in the traveltime data sets; however, sharper changes are generally identified in the northeastern half of the refraction line (approximately x > 30 m), similar to what was observed for surface-wave data (Bergamo et al., 2016).
In Figure 10, the temporal variability of P-wave traveltimes is statistically summarized (Figure 10a) and compared with the temporal trend of the parallel time-lapse set of surface-wave data (in particular, fundamental mode dispersion curve phase velocities, Figure 10b; Bergamo et al., 2016). Finally, the data temporal variability is related to the recorded precipitation pattern and earthwork volumetric water content measurements (Figure 10c-10e). Each row in Figure 10a and 10b represents a histogram of the distributions of point-by-point changes, relative to March 2014. Figure 10a and 10b shows a symmetrical distribution for the temporal changes of P-wave traveltimes as compared with Rayleigh-wave phase velocities, suggesting a reciprocally consistent seasonal trend for surface-wave and P-wave velocities (an increase in traveltimes corresponds to a decrease of Pwave velocities and vice versa).
This seasonal trend can be correlated with the pattern of rainfall data recorded at the site (Figure 10c). Comparing July 2013 with the following autumn-winter data sets (September 2013 to January 2014), P-wave arrival times globally increase (Figure 10a), whereas phase velocities reduce (Figure 10b), simultaneously with the frequent and intense autumn-winter rainfall (Figure 10c). The opposite effect may be observed in the spring-summer period of 2014 (May 2014 to September 2014), in which a decrease in P-wave traveltimes, an increase of phase velocity, and a corresponding decrease in rainfall are generally recorded, when compared with previous data sets. Data from November 2014, acquired after the start of intense autumn precipitations, show much slower phase velocities/higher Pwave arrival times, similar to those of autumn-winter 2013. A partial exception to this pattern is observed in May 2014, when a marked global decrease of arrival times corresponds to a slight reduction of phase velocities.
The correlation between the symmetrical seasonal trends of traveltimes, phase velocities, and precipitation is consistent with the expected variation of seismic velocities with saturation in soils Santamarina et al., 2005). Intense and frequent rainfall is likely to correspond to an increase of water content within the railway embankment, which in turn causes a decrease of interparticle capillary pressures (i.e., a reduction in suction pressure). This leads to a decrease in the effective stress acting on the soil particles, and hence to a reduction of seismic velocities in the subsoil (i.e., V S and V P below full saturation; Santamarina et al., 2001). Therefore, the autumn-winter period with more severe and frequent rainfall (Figure 10c) is characterized by higher P-wave traveltimes and lower Rayleigh-wave phase velocities (Figure 10a and 10b). Vice versa, the spring-summer months are characterized by alternating periods of very low or null precipitation and extreme rainfall events (Figure 10c). These events generally saturate the shallower portion of the soil (dried by preceding periods of low rainfall), but they do not penetrate deeper into the subsoil. Moreover, greater evapotranspiration occurs when temperatures are higher. Hence, spring-summer months are likely to be characterized by lower values of water content, which lead to a rise of interparticle capillary pressure and to a corresponding increase in effective stress. As a consequence, spring-summer periods bear higher Rayleigh-wave phase velocities and P-wave velocities (implying lower P-wave traveltimes; Figure 10a and 10b). The discrepancy represented by the May 2014 data can be explained by the fact that those seismic data were acquired on a day with particularly heavy rainfall (22 mm, highest daily value during data acquisition). This event may have locally caused the shallower portion of the embankment to become fully saturated, leading to an increase of V P  and hence to a decrease in traveltimes, particularly at short offsets (Figure 9, May 2014 panel). On the contrary, Rayleigh-wave phase velocities, being more sensitive to the S-than P-wave velocity structure (Xia et al., 1999), are less likely to be affected by the condition of full saturation. In addition to precipitation, the variability of water content within the earthwork may be also influenced by the presence of vegetation. For example, the relatively larger changes in traveltimes observed in the northeastern half of the line (Figure 9) can be ascribed to the presence of thick vegetation on the flanks of this par-ticular stretch of embankment, implying the extraction of water from the soil by the roots of plants.
The volumetric water content values recorded at the dates of seismic data acquisition (Figure 10d and 10e) cannot be systematically 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; Figure 1), and, moreover, they span only a portion of the investigation time period. Nevertheless, a significant variability over time of water content within the embankment structure is observed (from values as high as 0.55 to as low as 0.05). As expected, this variability correlates with the observed precipitation pattern (see in particular the general decrease of water content values toward September 2014, and their significant increase in November 2014). This observed variation of volumetric water content is indicative of markedly changing levels of saturation within the embankment body, and it is therefore compatible with the substantial temporal changes of seismic data.

Data temporal variability versus experimental uncertainty
The relevance of the observed variations of P-wave arrivals with respect to the data experimental uncertainty was analyzed. For the purpose of this test, the temporal variations between the September and November 2014 traveltime data sets were collated with the corresponding experimental uncertainties (Figure 3), both individually, to ascertain whether temporal changes exceed the uncertainty intervals, and globally, to compare their statistical distributions. As previously, only the data points corresponding to P-wave propagation through the embankment were considered. The majority (67%) of the observed point-by-point data variations between the September and November 2014 traveltimes are larger than the corresponding uncertainty range. Also, the global distributions of experimental uncertainties and temporal changes differ significantly ( Figure 11). Experimental uncertainties for the September and November data sets show similar patterns: Uncertainties appear to be symmetrically distributed around 0% and lying within relatively narrow intervals (90% of uncertainties being comprised within a AE10% range). Conversely, the distribution of temporal changes in the data appears to  be considerably different, significantly skewed toward positive values (90% of traveltime differences lie within −5% to 35%, denoting a general decrease of P-wave velocities). The temporal variations of P-wave traveltimes are therefore statistically significant with respect to their uncertainty intervals, and they reliably represent a seasonal change in the seismic properties within the earthwork structure.

Inverted P-wave velocity models
The analysis of the inversion results obtained for the March 2014 reference data set (see the "Results" section) suggests a good reliability of the adopted inversion strategy, confirmed by the consistency between the reconstructed V P model and the outcome of invasive surveys (CPT data) and by the satisfactory data compliance. Therefore, the V P sections derived from the inversion of all available data sets were used to track and quantify the seasonal trend in seismic properties of the embankment.
As expected, the V P models being obtained from the inversion of the presented sets of hodocrones, the reconstructed P-wave velocity sections exhibit a seasonal variability directly derived from the temporal pattern observed in the traveltime data sets (previous subsection). The V P models (Figure 7) generally show higher velocities within the embankment body (depths <5 m) in the spring-summer period (July 2013 and March to September 2014); oppositely, autumn-winter months (September 2013 to January 2014 and November 2014) are characterized by lower V P values. Figure 12a-12c quantitatively represents the relative variations of V P with respect to their temporal average in the three shallower layers (each 1.5 m thick), which model the embankment structure (0.75-5.25 m depth). The shallower layer of cells (depths <0.75 m) were excluded from the analysis due to poor ray coverage. It is again possible to globally identify a seasonal trend in V P , with a good lateral consistency observed among cells located along the seismic line (color scale). The V P values appear to vary within a range of approximately AE10%, with the exception of the upper layer, where some peaks reach AE20% (Figure 12a). This wide range of variation could be related to the greater exposure of this layer to weather effects (precipitation and evaporation) and hence to sharper changes of moisture content. However, it is also possible that this is related to the depth dependence on the influence of picking error of P-wave traveltimes.
As shown in Figure 3, greater relative experimental uncertainties are observed for short offset arrival times (mainly related to the estimation of V P of shallower layers), whereas smaller relative values are obtained for longer offset data (also related to the estimation of V P of deeper layers). The offset dependence of relative experimental uncertainty (and hence of the expected experimental error) can be ascribed to the error introduced by the manual picking of P-wave traveltimes having a greater relative impact on the small values of P-wave arrival times at short offsets. Hence, the larger values of experimental error for the short offset traveltimes are likely to produce greater relative estimation errors in the reconstruction of P-wave velocities for the shallower layer. This could explain the sharper trend and the higher peaks appearing in Figure 12a.
For comparison, the lower panels of Figure 12 (Figure 12d-12f) represent the temporal behavior of the reconstructed V S values for the same layers (depths 0.75-5.25 m), obtained from inversion of the parallel time-lapse data set of surface-wave data (Bergamo et al., 2016). The seasonal cycles observed in P-and S-wave velocities show a similar overall pattern within the same spatial range (i.e., with a reference point at 7 < x < 55 m, in the brown to orange colors in Figure 12). The embankment section covered by surfacewave data (for which V S models are available) is wider than the area surveyed with P-wave refraction (Figure 1). Significant discrepancies between the temporal behavior of V P and V S are present, and they can be ascribed to the different spatial resolution and sensitivity of the surface-wave and P-wave refraction methods (e.g., Pasquet et al., 2014Pasquet et al., , 2015b or to inversion artifacts. Another possible explanation is the occurrence of conditions of local soil saturation, which, at the investigation scale of P-wave refraction surveys, would result in increasing P-wave velocities that are not quite as high as those which would normally be expected for fully saturated media . This is probably the cause of the disagreement observed for the shallower layer in the May 2014 models: V P values increase with respect to March 2014, whereas V S values decrease (left panels of Figure 12).
As previously discussed, the May 2014 acquisition took place on a day with heavy rainfall, which may have caused local saturation of areas within the shallower portion of the embankment. The reconstructed temporal V P (and V S ) models of the earthwork structure suggest a cyclical change of its seismic properties, globally in agreement with the expected variability of soil moisture and resultant pore-water pressure; however, the achieved level of analysis does not enable conclusions to be drawn regarding the geo-mechanical mechanism taking place within the embankment fill. There are ambiguities regarding the lateral and vertical distribution of seismic velocity variations, and the coherence between reconstructed V P and V S variability inevitably suffers from the dissimilar subsurface coverage of surface wave versus P-wave traveltime data and the different sensitivity of the inversion algorithms.
To achieve an adequate understanding of the geomechanical behavior of the embankment fill, a set of laboratory tests on soil samples could be performed to establish definite relationships between the geophysical and geotechnical properties of the material. The resultant geomechanical relationships could then be conveniently extended to the whole surveyed section of the embankment using time-variant seismic models derived from inversion of the geophysical data. CONCLUSIONS P-wave refraction data were collected repeatedly along the crest of a railway embankment, with the aim of acquiring a time-lapse P-wave traveltime data set able to portray the seasonal variations of V P within the body of the earthwork. Seismic data were acquired bimonthly for a total period of 16 months (July 2013 to November 2014). To ensure the acquired seismic data were fully comparable, for each acquisition the same source type, receiver array, recording equipment, and acquisition geometry/parameters were adopted. Pwave first-break arrival times were extracted from all acquired data, to obtain multiple sets of hodocrones referenced to the same spatial location but temporally distributed over a considerable period of time. These traveltime data sets were then thoroughly analyzed to identify their temporal variability. The high quality and spatial consistency of the extracted P-wave hodocrones enabled a detailed "data point by data point" comparison method to be performed. A consistent pattern of data variability over time was observed, including an overall decrease of P-wave traveltimes in the springsummer months, and oppositely, higher arrival times during the autumn-winter period. This temporal trend correlates well with the seasonal distribution of precipitation measured at the site. The correlation can be explained by considering the seasonal variation of moisture content and the resultant pore-water pressure in the embankment fill and its effect on the seismic properties of the embankment. It was also verified that the observed temporal traveltime changes generally exceed the range of experimental uncertainties. The seasonal trend of P-wave arrival times is similar to that identified in a simultaneously acquired surface-wave data set.
Inversion of the P-wave hodocrones was also carried out to obtain a time-variant V P model of the embankment structure. A tomographic inversion approach was adopted, using the same initial model and level of spatial regularization for each data set. The reconstructed temporal behavior of V P within the surveyed earthwork is coherently influenced by the seasonal cycle of P-wave traveltimes, in global agreement with the expected trend of water content and pore-water pressure. There are also similarities between the retrieved temporal pattern of V P and V S models, the latter being obtained from the inversion of surface-wave data.
This work illustrates the potential application of P-wave refraction, for the time-lapse monitoring of seasonal effects in vulnerable geotechnical assets. In particular, the acquisition and processing stages have yielded promising outcomes, proving a significant and reliable sensitivity of P-wave traveltimes toward seasonal changes of geo-mechanical properties in geotechnical structures. Further work is, however, required to refine the final, inversion phase, where the achieved level of analysis does not allow a robust interpretation of the reconstructed temporal variability of seismic velocities in terms of geomechanical behavior. In this perspective, the necessity to carry out geotechnical tests to establish definite relationships among water content/seismic/mechanical properties of the subsoil materials and the need for an integrated interpretation of different types of seismic data (e.g., joint inversion) are proposed.