Numerical modelling of the effects of foundation scour on the response of a bridge pier

Foundation scour can have a detrimental effect on the performance of bridge piers, inducing a significant reduction of the lateral capacity of the footing and accumulation of permanent settlement and rotation. Although the hydraulic processes responsible for foundation scour are nowadays well known, predicting their mechanical consequences is still challenging. Indeed, its impact on the failure mechanisms developing around the foundation has not been fully investigated. In this paper, numerical simulations are performed to study the vertical and lateral response of a scoured bridge pier founded on a cylindrical caisson foundation embedded in a layer of dense sand. The sand stress–strain behaviour is reproduced by employing the Severn-Trent model. The constitutive model is firstly calibrated on a set of soil element tests, including drained and undrained monotonic triaxial tests and resonant column tests. The calibration procedure is implemented considering the stress and strain nonuniformities within the samples, by simulating the laboratory tests as boundary value problems. The numerical model is then validated against the results of centrifuge tests. The results of the simulations are in good agreement with the experimental results in terms of foundation capacity and settlement accumulation. Moreover, the model can predict the effects of local and general scour. The numerical analyses also highlight the impact of scouring on the failure mechanisms, revealing that the soil resistance depends on the hydraulic scenario.


Introduction
At the design stage of a bridge, the performance of its components has to be carefully evaluated from the hydraulic, structural, and geotechnical points of view. In addition, it is essential to ensure the stability and functionality of the structure during its service life, which can evolve due to both structural damage or ageing and variation of the supporting conditions [2].
Several studies have shown that the majority of bridge failures are a consequence of extreme flood events [30,63]. The presence of a structure in a channel constitutes an obstacle to the water flow, leading to the development of turbulences around piers and abutments. Such vortices increase the erosive power of the flow, which may locally exceed the shear strength of the riverbed gaining the energy to excavate and carry out the soil material. This phenomenon, usually termed local scour, keeps going on until local equilibrium is re-established [41,49]. In addition, aggradation and degradation phenomena, together with the reduction of the channel section, may result in global erosion of the riverbed-general scour. This erosive process typically takes place also around large monopile foundations, such as the ones usually employed for offshore wind turbines, affecting significantly their lateral behaviour [38].
The detrimental impact of scour on the mechanical response of bridge piers is nowadays well-recognized (e.g. [2]). Nevertheless, as noticed by Brandimarte et al. [9] bridge failures due to natural hazards and, in particular, foundation scour seem to be increasing with time. This is partly due to the large number of recent extreme flooding events, becoming increasingly frequent as a consequence of climate changes. In addition, budgetary constraints often prevent the timely implementation of risk mitigation measures, leading to an increase in vulnerability to extreme events [39]. As a consequence, the static and seismic responses of scoured foundations have to be re-assessed in the aftermath of major flood events.
In the state of practice, such assessment is typically performed through numerical analyses (e.g. [62]), often neglecting the hydraulic processes leading to foundation scour. However, as shown experimentally by Qi et al. [48] and Ciancimino et al. [15], the effects of general and local scour are markedly different. In addition, the assumed shape of the scour hole has been shown to have an impact on the predicted performance of the foundation [13]. Unfortunately, it is difficult to a priori define the most critical situation, especially when seismic analyses are performed. Consequently, it is essential to model the response of scoured foundations as realistically as possible.
In this study, the impact of scouring on the performance of a bridge pier supported on a cylindrical foundation embedded in a layer of dense Perth sand is investigated through numerical analyses. The simulations refer to the prototype bridge pier which was investigated through centrifuge modelling by Ciancimino et al. [15]. Nonlinear soil response is modelled employing the Severn-Trent constitutive model [20,22,23]. The model has been proven to be effective in predicting sand behaviour for a wide range of boundary value problems (e.g. [1,37,51]). The constitutive parameters are carefully calibrated on the basis of element tests of Perth sand samples. The calibration procedure is firstly performed by comparing the results of drained triaxial tests to the numerical predictions, obtained by imposing a ''true'' triaxial stress state on a single finite element. The procedure is then refined by modelling the laboratory tests as boundary value problems, thus including the effect of sample edge restraints. The numerical predictions are also compared to the results of resonant column tests, confirming the validity of the calibration procedure.
Finally, the carefully calibrated constitutive model is used to simulate a series of centrifuge model tests conducted at the ETH Zurich (ETHZ) Geotechnical Centrifuge Center (GCC), allowing a detailed comparison of numerical predictions to experimental results. The purpose of this comparison is twofold: (a) it allows the assessment of the reliability of the numerical simulations, serving as a validation of the numerical model, and (b) the combined use of numerical and experimental results provides validated insights on the mechanisms and the overall response of bridge piers, especially concerning the impact of scour on the failure mechanisms that develop in the soil around the foundation.

Characterization of the Perth sand
The soil used in the centrifuge model tests is a fine silica sand from Western Australia, close to Perth (termed hereafter, Perth sand) with a mean grain size d 50 = 0.23 mm and a soil particle specific gravity G s-= 2.65, characterized by a uniformity coefficient C u-= 1.79 and a curvature coefficient C c = 1.26. The limiting void ratios are e min = 0.50 and e max = 0.87. The reported values are based on previous experimental studies [10,45], except for e max for which a slightly higher value is adopted [15].
The stress-strain behaviour of Perth sand is discussed herein in the light of new element tests. Specifically, isotropically consolidated drained (CID) and undrained (CIU) compression triaxial tests were carried out on loose and medium-dense samples to estimate the position of the critical state locus. In addition, a series of resonant column (RC) tests was conducted to measure the small-strain shear modulus G 0 and the modulus reduction and damping (MRD) curves as a function of the effective confining pressure p 0 . The main findings of the laboratory tests are summarized in the next sections. More details on the testing procedures and the analysis of the results can be found in Ciancimino [14].

Triaxial tests
The main objective of the laboratory investigation is to identify the critical state line (CSL) of Perth sand in both the qp 0 and e -log p 0 planes. Defining its position is not a trivial task, especially in volumetric terms. A dense (i.e. dilative) sand sample sheared in both drained and undrained conditions will eventually reach the critical state at large axial strains, beyond the limits given by the assumptions adopted to interpret the results of laboratory tests. Moreover, the potential development of shear bands implies the necessity of defining the void ratio locally, even if the sample is sheared under undrained conditions with zero global volumetric strains (e.g. [16,60]). Conversely, a sample looser than critical tends to reach the critical state at low axial strains, without significant development of shear bands [16]. The void ratio of the sample is thus almost homogeneous throughout the specimen. Consequently, it was decided to mainly perform undrained tests on loose and medium-dense samples to study the position of the CSL, looking for total, or at least partial, liquefaction.

Experimental programme and testing procedures
The tests were carried out by using the triaxial testing apparatus of the Geotechnical Laboratory of Politecnico di Torino, specifically developed for testing sand samples under stress or strain control [12].
The specimens were prepared through the moist tamping method by adding a 5% weight of water to obtain loose, homogeneous samples with a diameter of 70 mm and a height of 140 mm. The sand was placed in a standard latex membrane (0.25 mm thick) applied via vacuum against the wall of a cylindrical split mould and levelled up to the desired height using a flat-bottom tamper. The procedure was repeated for 10 layers until the total specimen was created. After preparation, the samples were firstly saturated until B-values [57] higher than 0.97 were obtained and subsequently consolidated isotropically up to the initial mean effective pressure p 0 0 . Finally, they were sheared in compression under displacement control through a digital press at a constant rate of 5 mm/h. Stresses and strains were monitored with an acquisition frequency of 0.2 Hz.
To compute the deviatoric stress q, a cylindrical deformation of the specimen was assumed, according to the ASTM D7181-20 standard. Despite being valid for the preliminary test phases (saturation and consolidation), this assumption is not completely appropriate during shearing. Nevertheless, it was decided to assume an average q for consistency with the adopted volumetric e vol and axial e ax strains which are computed as average global values.
As summarized in Table 1, five triaxial tests were performed in total.

Results
The results of the triaxial tests are summarized in Fig. 1. The CIU tests are presented in terms of deviatoric stress q and excess pore water pressure Du w versus e ax , whereas for the CID test the results are presented in terms of q and e vol .
The first two undrained tests (TX-PoliTO-01 and 02) were performed on medium-dense samples with e 0 % 0:7, consolidated up to p 0 0 of 100 and 200 kPa, respectively. A dilative behaviour, typical of sands slightly denser than critical, is identifiable: an initial increase of Du w up to the phase transformation is followed by a decrease of Du w down to negative values and a consequent significant hardening ( Fig. 1). In TX-PoliTO-01, the denser state of the soil and the lower p 0 0 resulted in larger negative Du w compared to TX-PoliTO-02. The critical state condition for TX-PoliTO-01 is not clearly identifiable, even if the response seems to stabilize towards the end of the test. Conversely, in TX-PoliTO-02 the critical state seems to be reached for e ax % 0:08, after which Du w and q are almost stable.
Tests TX-PoliTO-03 and 05 were conducted on loose samples (e 0 = 0.79 and 0.84) consolidated up to p 0 0 = 300  and 200 kPa, respectively. In TX-PoliTO-03, the sample experiences temporary liquefaction: a rapid increase of q is followed by a softening and temporary stabilization of Du w (Fig. 1). This condition is referred to as quasi-steady state and is followed by further hardening [60]. The question of whether this is an actual material behaviour or rather a sample behaviour caused by the limitations associated with the testing procedures is still debated (e.g. [65,69]). Anyway, the two conditions (quasi-steady and steady state) are usually very close. The point was then included in the subsequent analyses for comparison purposes. Conversely, in TX-PoliTO-05 the very loose initial state of the sand leads to complete liquefaction of the specimen and an almost total loss of shear strength. After the peak, both q and Du w become stable and the sand reaches and maintains a clear critical state condition (Fig. 1).
The contractive behaviour of a loose (e 0 = 0.79) sample was investigated under drained conditions in test TX-PoliTO-04. The sand exhibits a monotonic increase of q until e ax % 0:1, after which q is almost constant and the critical state condition is reached (Fig. 1a). Such a trend is also visible in volumetric terms, where the sand experiences gradual contraction up to e vol % 0:01 (Fig. 1b). In this specific test, the initial state of the soil was probably really close to being critical, so small e vol developed during the shearing phase.
The critical state points (represented with circles in Fig. 1) are summarized in Table 1. For the two tests performed on medium-dense samples, the void ratios were not identified. Under undrained conditions, the tendency of medium-dense samples to dilate implies a decrease in the pore pressure until the critical state is reached (Fig. 1b). As observed by other studies (e.g. [26,40]), the latter could lead to cavitation of the pore fluid and, in turn, to partially drained conditions. No values of void ratios at critical state are thus reported for these tests in Table 1 and the points are included just for comparison purposes.

Critical state line
The CSL of Perth sand in the qp 0 plane is presented in Fig. 2a, also reporting the quasi-steady state for comparison. The slope of the line is obtained through least-square regression and it is equal to a critical state stress ratio M CS ¼ 1:23. The corresponding critical state friction angle is u 0 CS ¼ 30:6 AE 0:7 . An estimate of the CSL in the volumetric plane is presented in Fig. 2b. Both the quasi-steady and the critical state points are reported, together with the paths followed by the samples during the compression phase. The CSL is described by adopting the equation proposed by Been and Jefferies [8]: where v k and k are two calibration parameters. By expressing p 0 in kPa, the parameters are 1.89 and 0.02, respectively. The position of the line appears to be reasonable and the intercept is quite close to the maximum specific volume associated with the loosest possible state. However, given the limited number of points used for the regression, the proposed equation is considered just an estimate of the CSL.

Resonant column tests
Resonant column tests were carried out to investigate the nonlinear response of Perth sand in a wide strain range. The tests were conducted using the combined RC and Torsional Shear device of the Geotechnical Laboratory of Politecnico di Torino [46], which is a modified version of the free-fixed type apparatus designed at the University of Texas at Austin [31]. The tests were performed and processed according to the prescriptions given in ASTM D4015-15.

Experimental programme and testing procedures
The RC tests were performed on dense samples with e 0 as close as possible to the centrifuge tests (0:55 À 0:60). With a diameter of 50 mm and a height of 100 mm, the specimens were prepared through the dry pluviation, aiming to achieve high relative densities while maintaining an acceptable degree of homogeneity. After preparation, the samples were saturated (in the case of undrained tests) and isotropically consolidated. The loading was imposed under load control, applying torques with increasing amplitude through an electromagnetic driving system. For a given amplitude, a wide range of frequencies was analyzed to identify the first torsional resonance frequency of the sample f 0 . The shear wave velocity V s and the secant shear modulus G S were computed by applying the theory of torsional wave propagation under steady-state conditions in a linear visco-elastic medium [50]. The damping ratio D was obtained by applying the free-vibration decay method. The dynamic parameters of the soil are referred to a cyclic shear strain amplitude c c equal to 2/3 of the maximum shear strain [29,64]. A summary of the tests is given in Table 2. For each test, the experimental programme comprised several stages, each one consisting of initial isotropic consolidation to p 0 0;i , followed by a full (i.e. from small to large strains) RC test performed in undrained mode. During the test, the eventual pore water pressure build-up was recorded. At the end of each stage, before starting a new consolidation phase, the drainage was opened to allow dissipation of Du w . Test RC-PoliTO-01 comprised 4 stages carried out at p 0 0 ranging from 25 kPa up to 200 kPa; 5 stages were conducted in test RC-PoliTO-02, with a maximum p 0 0 of 400 kPa. The procedure adopted for the first two tests allowed investigation of the influence of p 0 on the modulus reduction and damping (MRD) curves of Perth sand, similarly to previous studies (e.g. [43]). However, such a procedure is not appropriate to study the evolution of G 0 with confining pressure. When a sand sample is excited in undrained conditions above the volumetric threshold, it exhibits a pore pressure build-up. Consequently, during the dissipation phase, the sample experiences permanent changes in the microstructure and further densification [56].
A third experiment (RC-PoliTO-03) was performed on a dry sample, following a different 16-stages procedure. The first 10 stages (RC-PoliTO-03-A in Table 2) were performed at low strain amplitudes, so as not to exceed the linearity threshold, thus preventing any change in the structure of the sample. Such stages allowed the definition of G 0 evolution for confining pressures ranging from 25 to 800 kPa. When the maximum desired p 0 0;i was reached, a full (from small to large shear strains) RC test was carried out. The confining pressure was then progressively decreased in further stages (RC-PoliTO-03-B in Table 2), and for each stage, the MRD curves of the sample were obtained.

Small-strain shear modulus
The small-strain shear modulus G 0 of sand mainly depends on the actual state of the material. Hardin and Black [28] firstly pointed out this dependency, suggesting to express G 0 as a function of p 0 and e. Figure 3 presents the experimental data in function of normalized confining pressure p 0 =p atm , with superimposed the relationship of Hardin and Black [28], calibrated for this data set. The obtained A G is 31.8 (with G 0 in MPa), which is in good agreement with values suggested for rounded sands, while for the exponent n G a typical value of 0.5 is adopted for the sake of simplicity.

Modulus reduction and damping curves
The nonlinear behaviour of soils is usually represented through the MRD curves, exploiting the dependency of the visco-elastic parameters on the cyclic strain amplitude c c [53]. The main parameters affecting the MRD curves of sands are p 0 and particle gradation [42,54]. The experimental curves are thus usually referred to p 0 0 , corresponding to the initial value measured at the end of the consolidation phase. Nevertheless, p 0 can significantly change throughout testing under undrained conditions, due to the pore pressure build-up. It is then crucial to normalize each point of the MRD curves to the actual p 0 i measured at resonance, in such a way to consistently compare results from different tests (see also [14]).  The normalized MRD curves resulting from the RC tests are depicted in Fig. 8 (presented in the next section) as a function of p 0 i . The increase of p 0 0 leads to a shift of the MRD curves to the right, widening the almost linear range. In addition, an increase of p 0 0 results in a decrease of the small-strain D. However, the influence of p 0 seems to be relevant for small pressures (\ 100 kPa), becoming less significant at higher pressures consistently with previous findings (e.g. [54]).

Soil constitutive model
Nonlinear sand response is modelled employing the Severn-Trent sand model, as in the version proposed by Gajo [20]. Formulated within a bounding-surface approach, the original model [22,23] can reproduce the stress-strain response of sands for a wide range of states with a single set of parameters. It relies on two surfaces: the strength surface and the yield surface, both defined through the Argyris et al. [4] criterion. The size of the strength surface is defined as a fraction of the critical state surface through ratio r: where k is a constitutive parameter and W is the state parameter, as proposed by Been and Jefferies [8]: being v and p 0 indicative of the current state of the sand. The size of the strength surface can then expand or shrink according to the evolution of the state of the material, while the position of its axis is fixed. The yield surface can undergo both isotropic and kinematic hardening and its size is defined as a fraction of the strength surface according to constitutive parameter R.
Although equivalent to the original one, the version employed herein [20] is based on the framework of the classical kinematic-hardening plasticity, adopting a hyperelastic formulation to describe the soil behaviour inside the yield surface. Such formulation is modified as proposed by Argani and Gajo [4] to improve the rate of convergence for negligibly small confining pressures. At yielding, the magnitude of the plastic strain increment is related to the distance between the current stress state and an image stress state located on the strength surface through a hyperbolic relationship, based on constitutive parameter B. The stress-strain response is smoothed in the transition from the elastic to the elastoplastic zone through the introduction of a further constitutive parameter a. The model assumes the deviatoric associativity to link volumetric and distortional strain increments through the flow rule. The latter is defined based on two constitutive parameters: the multiplier A and a parameter k d , defining the contribution of W to the plastic flow.

Model calibration
The constitutive model is completely defined after thirteen constitutive parameters are calibrated. The parameters are summarized in Table 3 along with the values representative of the Perth sand response.
The first three parameters, v k , k, and u 0 CS , define the position of the critical state surface and are selected on the basis of triaxial test results (see Fig. 2). In the present study, the inherent elastic anisotropy of the material (and its evolution) is neglected, so the initial fabric tensor B 0 is assumed constant and equal to I. The remaining parameters that control the elastic behaviour are derived by adopting Poisson's ratio t = 0.2 and G 0 defined through the Hardin and Black [28] equation (Fig. 3). Previous studies (e.g. [1,23]) assumed a fraction E R of G 0 instead of its full value as a strategy to better capture sand response under undrained conditions. The drawback of this strategy is that the drained behaviour, in the very small-strain range, could not be adequately reproduced. Bearing in mind that the target of this study is to reproduce the response observed in the centrifuge tests, it was decided to focus on the drained response, adopting G 0 . Consistently with this assumption, a small size of the yield surface R ¼ 0:01 is selected (i.e. 1/100th of the strength surface).
Parameter m defines the deviatoric shape of the yield and strength surfaces through the Argyris et al. [4] criterion. A typical value of 0.8 is adopted to ensure the convexity of the surfaces. The remaining parameters are calibrated through a trial-and-error process, based on the results of the element tests. In particular, the calibration of A, B, k, and k d is performed on the basis of drained triaxial

Simulation of drained triaxial tests
The calibration of the constitutive model is performed mainly focusing on the results of drained tests. In particular, reference is made to test TX-PoliTO-04 (Table 1) and to the drained triaxial tests conducted by Bühler et al. [11] on dense Perth sand samples with e 0 very close to those measured in the centrifuge tests (Table 4).
To interpret the test results, stresses and strains are computed in function of average external measurements, assuming cylindrical deformation of the specimen. Although this is a common assumption, it deserves an indepth discussion. The effects of external boundary conditions can have a significant impact on the measured response. The rough end plates are known to impose nonuniform stress and strain fields within the specimen, which tend to lead to barrel-shaped deformation. Such a typical barrel shape was also observed in the performed Perth sand triaxial tests. Figure 4 compares the shape at failure of an initially denser (Fig. 4a) to a looser (Fig. 4b) than critical Perth sand specimen. Strain localizations can readily be identified in the dense sample.
This issue has been analysed by several experimental and numerical studies (e.g. [44,55]). The presence of the rough end plates leads to an apparent increase in shear strength (e.g. [35]). Moreover, the end-restraint prevents full development of volumetric strains throughout the specimen, which in the case of dense sand tends to localize within the shear bands [16]. As discussed by Sheng et al. [55], such boundary effects should be taken into account for proper constitutive model calibration. Following the notation introduced by Muraro and Jommi [44], the drained triaxial tests are simulated considering both material behaviour under true triaxial stresses and strains, and sample behaviour, treating the experiments as boundary value problems. Tensor controlling the inherent anisotropy I  Fig. 4 Typical barrel-shaped deformation observed on: a dense [11], and b loose sand specimens A single CAX8R (i.e. 8-node axisymmetric quadrilateral, biquadratic displacement, reduced integration) element subjected to an initial isotropic triaxial state is sheared under drained conditions to simulate the material behaviour (Fig. 5a). Conversely, the sample behaviour is simulated by imposing the external boundary conditions on a quarter of the soil sample, modelled in axisymmetric conditions, and discretized as shown in Fig. 5b. The numerical model includes both the specimen and the steel top plate, which is assumed linear elastic with Young's modulus E S ¼ 210 GPa and Poisson's ratio t ¼ 0:3. A rough interface (infinitely rigid in the normal direction) with friction coefficient l ¼ tan u 0 CS is introduced between the specimen and the top plate. Preliminary sensitivity analyses showed that in terms of global stresses and strains, the results are not sensitive to the assumed l, neither on the mesh size [14]. After initial isotropic consolidation to p 0 0 , the sample is sheared by applying an axial displacement to the top plate under drained conditions (i.e. imposing Du w ¼ 0).

Material behaviour
In the simulations performed to evaluate the actual soil material behaviour, stresses and strains are uniform throughout a single element. The results can thus be directly extracted from the element integration points. The numerical results are compared to the measurements in Fig. 6a. The simulated deviatoric response of the dense samples is quite satisfactory for the entire range of e ax , although some small discrepancies can be observed. In particular, the simulations predict slightly lower q peak , after which a slower (compared to the experimental data) softening takes place. However, in terms of e vol the experimental results are captured only up to e ax = 0.1. For larger e ax , an apparent reduction of the e vol increment can be observed in the experimental measurements, which is attributed to strain localization (Fig. 4a).
Regarding the response of the initially looser than critical specimen, the comparison is satisfactory for the entire range of e ax in terms of progressive hardening of the material towards the critical state. In contrast to the dense sand experimental results, no significant strain localization is observed on initially loose samples (Fig. 4b). Consequently, in the case of loose sand, the interpretation of test results in terms of global quantities is considered reliable even at large e ax .

Sample behaviour
The sample behaviour is analysed considering conventional stresses and strains, computed for the entire specimen consistently with the procedures adopted to analyse the experimental results. Figure 6b compares the specimen simulation to the experimental results. In contrast to the single element analysis (material response), the experiments conducted on samples denser than critical are well captured for the entire range of e ax in terms of q peak and post-peak softening. In addition, the predicted volumetric responses exhibit an apparent plateau at large strains, in excellent accord with the experimental results. No significant differences are observed between material and sample response for the initially looser than critical specimen. It is well known that strain localization is typically associated with a pathological dependency on the FE mesh. For this reason, the numerical simulations were repeated with different mesh refinements (not shown here for the sake of shortness) which led to substantially unchanged results in Fig. 6b [14].
Direct visualization of the differences between material and sample behaviour is offered in Fig. 7, which presents numerically predicted contours of void ratio e for a specimen initially denser (TX-ETH-04, Fig. 7a) and looser (TX-PoliTO-04, Fig. 7b) than critical, for global axial strain e ax = 0.15. The confinement offered by the steel end plates leads to a pronounced bulging of the initially dense specimen, whereas a smoother barrel shape is observed for the looser specimen. In the dense specimen, a clear strain localization can be observed, starting from about a quarter of the height and developing towards the middle of the specimen. In the elements located inside the shear band, e reaches quite a high value (0.81), in contrast to the socalled dead zones close to the top (and the bottom) of the specimen, where e is slightly higher than the initial one (0.51). In stark contrast, no significant nonuniformities can be observed in the e field of the loose sample.
The results are consistent with previous experimental and numerical findings, where specimens initially denser than critical presented an X-shape localization of shear deformations resulting in the formation of shear bands [16,55]. Such shear bands formation implies a concentration of shear stresses, followed by faster softening of the global deviatoric response. On the other hand, it induces an apparent final plateau in the global volumetric strain response. The latter has to be interpreted as a manifestation  Numerically computed contours of void ratio e distribution for a global axial strain e ax = 0.15 for: a specimen initially denser than critical and b initially looser than critical of strain localization, rather than the attainment of critical state condition for the entire specimen [16]. The sand attains critical state only locally, within the developing shear bands.

Modulus reduction and damping curves
Cyclic direct simple shear (CDSS) tests are simulated to compare the numerical predictions to the experimentally measured MRD curves. The numerical analyses are performed on a single C3D8R (8-node trilinear displacement, reduced integration) element. The analysis starts with an initial consolidation phase up to p 0 0 , obtained by increasing the vertical effective stress r 0 z under oedometric conditions, followed by the shearing phase simulated by applying a cyclic horizontal displacement to the top nodes of the element under drained conditions.
The loading protocol consists of 10 cycles applied for 13 increasing displacement amplitudes, to adequately describe the range of shear strains investigated by the RC tests. The dynamic parameters of the soil are computed from the 10th load cycle for each displacement amplitude.
The numerical predictions are compared to the experimental data in Fig. 8. Although the present study focuses on the monotonic response of embedded foundations, the comparison is presented also in terms of damping curves for completeness. The smoothed variations of the curves with increasing c c are well captured by the constitutive model as well as their dependency on the confining pressure p 0 0 . The only noticeable difference to the experimental results is a slightly faster decay of G S =G 0 for c c ranging from 10 -5 to 10 -4 , and the corresponding overestimation of D at large strains, as typically predicted by nonlinear soil models.

Numerical simulations of centrifuge experiments
The thoroughly calibrated soil constitutive model is employed to simulate the centrifuge tests of Ciancimino et al. [15], focusing on the effects of scouring on the vertical and lateral response of a caisson foundation. Recognizing their ability to withstand severe earthquakes while suffering admissible damages in terms of permanent displacements and rotation [66], the seismic performance of caisson foundations has recently received growing attention (e.g. [24,25]). However, the reduction of embedment due to scouring may lead to a substantial deterioration of their performance [68].

Centrifuge modelling
An idealized prototype bridge pier is considered, consisting of a single degree of freedom system supported on a cylindrical caisson foundation of diameter D f ¼ 2 m, embedded by d ¼ 2 m in a layer of dense (D R % 80%) Perth sand (Fig. 9). To focus on the foundation response, the pier is modelled as a practically rigid slender structure. The physical model tests were conducted employing a decoupled hybrid approach [15]. The first step modelled in 1 g conditions the hydraulic process of local scour around the pier foundation, employing a miniaturized tidal generator (MTG) developed at ETHZ [33]. The morphology of the scour hole ( Fig. 10) was acquired through 3D-scanning and used to 3D-print moulds to realistically reproduce its geometry in a centrifuge model, tested in the ETHZ drum centrifuge [58] to investigate under Ng the mechanical part of the problem. The tests were performed under drained conditions and checked through pore pressure transducers embedded in the sand. The vertical response of the foundation was investigated through vertical push tests carried out at 100 g. The experimental setup consisted of an aluminium model of the foundation pushed with a constant velocity. The imposed vertical settlement w and the corresponding load V were measured to obtain the V À w response of the foundation before and after local scour (tests V1 A and V2 A , respectively). Displacement-controlled monotonic pushover tests were conducted at 50 g to study the lateral foundation response. The aluminium foundation was rigidly connected to the bridge pier, modelled using a practically rigid brass column, which was monolithically connected to a steel mass on the top representing the deck. The mass was selected so as to achieve a safety factor under purely vertical loading FS V % 8, corresponding to a moderately loaded system. The FS V is computed as follows: where: W 0 3 , W 2 , and W 1 are the weight of the footing accounting for buoyancy, the weight of the column, and the weight of the deck, respectively.
The lateral displacement was applied using an electric screw-jack actuator, connected to the model structure through a sliding hinge, equipped with a pinned connection, free to move in the vertical direction. The horizontal force was measured through a high-resolution load cell and used to obtain the overturning moment M, while the inplane movement of the pier (in terms of imposed rotation h and vertical settlement w) was monitored by three laser sensors.
The tests considered in the present study investigated the lateral monotonic response before foundation scour (test MH1), after local scour (tests MH2 and MH3), and after general scour (test MH5), corresponding to generalized erosion of a 1-m-thick soil layer (i.e. the average scour depth due to local scour). Under local scour conditions, the system response was studied along the transverse direction of the bridge, which is usually critical. Given the asymmetry of the scour hole along this direction (Fig. 10), it was necessary to carry out two tests, pushing the structure towards the upstream and downstream directions. An isotropic linear elastic constitutive law is adopted for the caisson, assuming the properties of aluminium, E = 70 GPa, t = 0.33, and c = 27.5 kN/m 3 . As in previous studies (e.g. [66,67]), the caisson is modelled as ''wished in place''. This simplification was done to focus on the effects of scouring but is also consistent with the preparation of the physical models. Moreover, the numerical simulations (and the centrifuge model tests) neglect the loading-unloading induced by progressive erosion of soil along the sides of the foundation. However, such a process is practically negligible, since the maximum thickness of eroded soil is 1 m, corresponding to Dr 0 v % 10 kPa. A purely frictional contact algorithm is employed to model the soil-caisson interface, thus allowing sliding and detachment. A ''hard'' pressure-overclosure relation is used Fig. 9 Idealization of the prototype bridge pier of the centrifuge tests (R.P. stands for Reference Point) Fig. 10 Morphology of the local scour hole (in prototype scale) obtained through 3D-scanning of the 1 g physical model tested in the MTG [15] in the normal direction, while the interface response along the tangential direction is modelled with a ''penalty'' algorithm, assuming a friction coefficient l ¼ 0:23. The latter is based on the results of a series of direct shear interface tests [52].

Numerical modelling aspects
All numerical simulations are performed accounting for geometric nonlinearities to capture, in the lateral tests, the additional moment DMðdÞ given by second-order (i.e. P À d) effects. It is worth noting that large displacements are anyway expected to occur only along the (practically rigid) bridge pier. The nonlinear geometry has therefore practically no impact on the soil response.

Simulation of vertical push tests
As demonstrated in several previous studies (e.g. [3,19,27]), the lateral response of a foundation system is strongly dependent on its vertical loading, usually described by its safety factor against purely vertical loading FS v . Therefore, simulating the response of the foundation under pure vertical loading is an essential first step.
The simulations of the vertical push tests consist of the following steps: (i) definition of initial equilibrium through a geostatic step; (ii) ''wished in place'' installation of the caisson; and (iii) vertical push under drained conditions. Consistently with the centrifuge tests, the analyses are performed under displacement control. As shown in Fig. 11, the d = 2 m embedded foundation extends 1 m above the ground surface (i.e. its total height is 3 m) to avoid numerical instabilities due to the ''collapse'' of soil along sides. The additional frictional resistance of this upper 1 m is practically negligible, and therefore the comparison is considered consistent, provided that this part of the foundation is modelled weightless.
A 2D axisymmetric FE model is employed to study the response of the foundation before and after generalized scour (Fig. 11a), while a 3D FE model is used to model the foundation subjected to local scour (Fig. 11b) due to the irregular shape of the scour hole. Taking advantage of the symmetry along the direction of the water flow (Fig. 10), half of the domain is modelled in 3D to reduce computational cost. The dimensions of the discretized soil field match the geometry of the cylindrical strongboxes used in the centrifuge tests. The bottom nodes are fixed in the vertical direction, while the side nodes are free to move vertically. After preliminary sensitivity analyses [14], the element size is equal to 20 cm close to the foundation, progressively increasing towards the model boundaries. The 3D mesh is defined through a revolving procedure built-in ABAQUS, so as to be equivalent to the axisymmetric model. The discretization of the scour hole geometry is obtained through a bottom-up manual procedure, based on the actual morphology obtained from the MTG tests (Fig. 10).
To allow reaching very high-stress gradients in the soil close to the foundation avoiding the typical numerical issues related to volumetric and shear locking, the soil and foundation materials are modelled through CAX8R elements for the axisymmetric models and C3D20R (i.e. 20-node brick, triquadratic displacement, reduced integration) elements for the 3D model.
The numerical predictions (solid lines) are compared to the centrifuge test results (dashed lines) in Fig. 12 in terms of load-settlement (V À w) and tangent vertical stiffnesssettlement (K v À w) response. The comparison also includes the numerical results for generalized scour, for which no experimental result is available. The simulations predict adequately well the evolution of K v with increasing w (Fig. 12b). The only noticeable difference is observed for the small-strain range (w \ 0.1 m), where the experimentally measured K v is significantly lower than the numerically predicted. With the increase of the imposed settlement w, the numerical curves converge towards the experimental measurements. Similar conclusions can be drawn by the V À w response (Fig. 12a), where the experimental V initially increases at a slower rate compared to the numerical prediction, progressively converging with the increase of w. A possible explanation of this discrepancy lies in the experimental procedure. To allow the caisson to freely settle during spinning up, the actuator arm is initially not in contact with the foundation. After attaining the target g-level the arm is pushed with constant velocity and contact is established. It is then possible that during the very first part of the vertical push, the contact is not yet perfect made, leading to larger measured w. The latter result in a lower initial value of K v , which apparently increases at small strains. Nevertheless, soon after the contact is established, the issue is overcome and K v starts to decrease. Therefore, the observed discrepancies in the small strain K v may be due to such initial imperfect contact: ''seating'' effect.
As also observed in the centrifuge tests [15], a distinct plateau cannot be identified in the V À w response (Fig. 12a). For w\0:1 m, the response is only slightly nonlinear becoming strongly nonlinear with the increase w, as the failure mechanism is being mobilized. After the failure mechanism is fully developed, K v tends to a residual (low, but nonzero) value. Consequently, the V À w response exhibits a progressive hardening regime, suggesting the mobilization of a local shear failure mechanism rather than a general one [61]. In the absence of a clear plateau, the definition of the vertical ultimate load V ult is not obvious. When a foundation exhibits such a hardening regime, what effectively limits its V ult is the accumulated settlement. Therefore, V ult is defined herein as the load necessary to produce settlement w ult = 0.3 m, which corresponds to a fairly large normalized settlement w ult =D f ¼ 15%. According to previous studies [3,36], this point is practically coincident with the initiation of the hardening regime after which K v is almost constant (Fig. 12b). Table 5 summarizes the V ult obtained from numerical simulations and centrifuge tests. The observed differences can be quantified as follows: The comparison is very satisfactory, with a maximum difference of ? 3.5% for the local scour condition. The numerical model successfully predicts the response of the caisson foundation under pure vertical loading and, most importantly, the (modest) effect of local scour. Local scour leads to a slight reduction of V num ult from 8.50 to 8.13 MN, due to the reduction of the effective overburden pressure in the soil below the foundation. The impact of general scour is much more pronounced in terms of both stiffness and bearing capacity (Fig. 12), leading to a 30% reduction of V num ult (from 8.50 to 5.95 MN).

Simulation of lateral pushover tests
The numerical predictions are compared to the results of the lateral pushover tests. A correction is applied to the experimental results to account for the effects of the radial gravity field on second-order additional moment DMðdÞ acting on the structure (see also [14]). The numerical model includes a rigid massless column representing the pier, the caisson foundation, and the deck and pier masses, modelled as vertical loads applied at the corresponding centres of mass of the two parts (Fig. 9).
The simulations include four steps: (i) definition of initial equilibrium through a geostatic step; (ii) ''wished in place'' installation of the caisson foundation; (iii) activation of vertical loads; and (iv) horizontal pushover loading under drained conditions.  An increasing horizontal displacement is imposed on the rigid pier at a height h act = 6.1 m to the base of the foundation, with the structure allowed to freely settle and rotate. The rocking response of the pier is calculated considering the in-plane rotation h and the settlement w, computed concerning the reference point (i.e. the centre of the base of the caisson, Fig. 9). The horizontal force is used to compute the overturning moment M: The developed 3D FE model of half of the domain is shown in Fig. 13. Consistent with the centrifuge tests, the distance between the base of the caisson and the bottom boundary of the model is 3D f ; the same applies to the distance of the caisson to model lateral boundaries. The latter has been shown to be sufficient to avoid boundary effects, in compliance with previous studies (e.g. [27,67]). To reduce the computational cost, a mesh size varying from 20 to 40 cm is adopted up to a distance of 1D f from the foundation, increasing to 80 cm further away. The soil and the caisson are modelled as previously described, while the pier is modelled as a discrete 2-node rigid beam (RB3D2) rigidly connected to the top of the caisson.

Results
The effects of scour on lateral pier-foundation response are analysed in terms of moment-rotation (M À h), settlementrotation (w À h), and secant rocking stiffness-rotation (K r;m À h) response. The numerical predictions (right) are compared to the experimental measurements (left) in Fig. 14 for the four hydraulic conditions examined: (a) before scour; (b) general scour; and (c) local scour around the pier, pushing the structure towards the upstream or downstream side.
The M À h response (Fig. 14a) is used to define the lateral pushover capacity of the structures. For a given rotation h, the actual resisting moment is theoretically given by the total moment M tot acting on the foundation, including the additional moment DMðdÞ due to P À d effects. However, the capacity of a foundation is usually defined as the maximum external load that the structure can withstand before suffering damage that could compromise the safety of the superstructure. Therefore, the moment capacity M ult is defined according to the external M À h response, thus implicitly considering second-order effects, as in previous studies (e.g. [19,47]). To allow for direct comparison among the different conditions considered, M ult is defined as the moment that develops for h ¼ 0:02 rad. The latter coincides approximately with the rotation after which the M À h curves start to decrease due to P À d effects.
The moment capacity obtained from the numerical simulations (M num ult ) and the centrifuge tests (M exp ult ) are compared in Table 6, along with the corresponding FS v computed through Eq. (4) according to the V ult reported in Table 5. The differences between numerical and experimental results are quantified as: The M À h response (Fig. 14a) and the effects of scour are very well predicted by the numerical model up to fairly large h. Their discrepancies range from -0.9% for the unscoured case to ? 11.4% for general scour. The effects of local scour are detrimental, leading to a reduction of M ult of about 25% when the structure is pushed towards the upstream direction. The effect of general scour is even more pronounced, leading to 43% reduction.
The numerical prediction is satisfactory also in terms of K r;m À h response (Fig. 14c). General scour leads to a significant decrease of K r;m for the entire range of rotations. The reduction is less pronounced for local scour, but still significant. The reduction of K r;m is more obvious for small rotations. For h ¼ 0:001 rad, K r;m is reduced from 115 MNm/rad before scouring, to 90 MNm/rad and 60 MNm/ rad in the presence of local and general scour, respectively. This result is qualitatively consistent with what has been observed from monitoring data of real scoured bridges (e.g. [17,34]).
Regarding the w À h response (Fig. 14b), the simulation slightly overpredicts the dilative response of dense sand, partially compensating for the plasticity-induced settlements, leading to moderate uplifting of the caisson for h [ 0:02 rad. However, it should be noted that the numerical prediction is overall consistent with the experimental results. Both in the numerical simulations and the centrifuge tests, the response is neither sinking nor uplifting dominated. As already observed in previous experimental studies (e.g. [3]), such behaviour is reasonable for a moderately loaded foundation (FS V % 5 Ä 10) on dense sand, where the effects of the two mechanisms (dilation and plasticity-induced settlement) are comparable. Under these conditions, the moderate effects of local scour on foundation bearing capacity (i.e. reduction of FS V from 8.0 to 7.7) do not affect the accumulation of settlement to a noticeable extent. In contrast, generalized scour leads to a moderate increase of the plasticity-induced settlement (Fig. 14b), resulting from the substantial decrease of FS v to 5.6. Although general scour is not very relevant on the vertical response of the bridge pier analysed herein, it could become substantial for heavily loaded structures, where a similarly pronounced decrease of FS v could lead to a transition to a sinking-dominated response.

Failure mechanisms
The total moment capacity of an embedded foundation is the result of two distinct resistance mechanisms: (a) base resistance (M base ) and (b) sides resistance (M sides ). In reality, the two mechanisms are coupled, but considering their effects separately allows for deeper insights into foundation response. The rotation of the structure leads to an increase of soil thrust (approaching passive conditions) along the pushing side of the foundation and a decrease along the opposite side (approaching active conditions). The corresponding sides' resistance M sides is then the result of passive-and active-type failure mechanisms [18,67]. Figure 15 compares the numerical failure mechanisms of the unscoured caisson to the cases of local (structure pushed towards the upstream direction) and general scour.
The results are presented in terms of deformed mesh with superimposed total displacement u tot (Fig. 15a) and void ratio e (Fig. 15b) contours for h ¼ 0:02 rad. The u tot contours (Fig. 15a) reveal an almost purely rotational scoop failure mechanism, as it would be expected for pure moment loading of a rocking-dominated system [67]. Below the edge of the foundation on the pushing side, soil yielding induced by the base resistance mechanism leads to dilation, and consequently to an increase of the void ratio e of the initially dense sand (Fig. 15b). This effect is more pronounced when the foundation is subjected to general scour, due to the more pronounced reduction of FS v ; the effect of local scour is much less pronounced. Active-and passive-type soil wedges can be observed along the sides of the caisson. While the active-type soil wedge is clearly visible, the passive-type wedge does not seem to be fully mobilized.
To quantify the relative importance of the two resisting mechanisms, the total resisting moment is subdivided into M base and M sides in Fig. 16. M base is computed by multiplying the normal nodal contact forces acting on the base of the caisson for the corresponding loading arm (measured at the reference point, R.P. in Fig. 9). M sides is then obtained as the difference between the total moment (including second-order effects) and M base .
For all the cases examined, M base is almost fully mobilized for h % 0:02 rad (as revealed by the attained plateau), indicating that the corresponding failure mechanism is fully developed at this stage. As expected from the vertical push tests (Table 5), local scour does not significantly affect this mechanism (Fig. 16a). This is not the case for general scour, which leads to a significant (14%) reduction of M base from 560 to 480 kNm. In stark contrast, M sides exhibits more of a hardening response with increasing h (Fig. 16b) for all cases examined. Interestingly, local and general scour have a markedly different effect on M sides , despite the fact that almost the same depth of soil is removed. An impressive 83% reduction is observed for general scour, with M sides (for h ¼ 0:02 rad) decreasing from 180 to 30 kNm. Local scour also leads to a significant but less pronounced reduction, ranging from 51% (M sides = 89 kNm) to 72% (M sides = 50 kNm) when the structure is pushed towards the downstream and upstream side, respectively. These significant differences with respect to the direction of loading, which are function of the asymmetry of the scour hole, are only applicable to M sides . For this specific caisson foundation, M base is dominant, and therefore the effect on total resistance is limited. This is not necessarily the case for deeper embedded foundations, where the contribution of M sides is usually dominant. The effects of local and general scour on M sides can be explained by looking at the stress fields in the soil elements along the sides of the foundation for h ¼ 0:02 rad (Fig. 17). Under general and local scour, soil erosion leads to a substantial decrease of soil pressures along the pushing side (Fig. 17a, b) and, consequently, of M sides (Fig. 16b). However, r 0 v (and consequently r 0 h ) is substantially higher for the local scour condition due to the presence of the soil slopes of the scour hole. This explains the differences observed between the M À h response obtained for the two hydraulic conditions (Fig. 16b).
The obtained earth pressure distributions are qualitatively consistent with the analytical solution proposed by Fröhlich [18]. The horizontal and vertical effective stresses follow a parabolic shape for all cases examined, with a maximum roughly at mid-height. Although the slope of the r 0 v distribution is almost equal to c 0 close to the top-left edge of the foundation, it is substantially higher along its pushing side (Fig. 17a). As the caisson rotation increases, the sides of the caisson are no longer vertical and shear stress occurs in the soil close to the footing. As a result, the initially horizontal contact pressure also begins to rotate with its vertical component leading to the observed increase of r 0 v . As a consequence, the distribution of r 0 h along the pushing side is also increased (Fig. 17b), explaining the hardening response of M sides observed in Fig. 16b. In contrast, r 0 h rapidly decreases to almost zero along the opposite side of the caisson. The distributions of the coefficient of lateral earth pressure K ¼ r 0 h r 0 v (computed locally for each integration point of the elements in contact to the foundation along the pushing direction) confirm these trends, which are the result of active-and passive-type failure mechanisms developing at the two sides of the foundation. While the active mechanism is fully mobilized for h ¼ 0:02 rad (K % 0:15 Ä 0:20), the passive mechanism along the pushing side appears to require further deformation to fully develop (K % 3).

Conclusions
This paper studied numerically the influence of foundation scour on the performance of a slender bridge pier supported on a cylindrical caisson foundation embedded in a layer of dense sand. The numerical simulations were conducted within the ABAQUS FE environment, employing the  Severn-Trent constitutive model [20,22,23] to simulate the nonlinear sand response. The model was initially calibrated on the basis of soil element tests, paying particular attention to the simulation of drained triaxial tests conducted on dense samples of Perth sand. The calibration process revealed the need to model the laboratory tests as boundary value problems due to stress and strain nonuniformities developing within the specimen. The FE model was then validated against the results of centrifuge tests conducted at ETH Zurich [15]. The numerical predictions were shown to be in good agreement with the experimental results. Specifically, the numerical models are able to reproduce both the vertical and the horizontal response of the foundation in terms of capacity, stiffness, and settlement accumulation. Most importantly, the comparison is satisfactory when looking at the effects of local and general scour. This is of particular relevance, given the independent calibration of the model based on soil element tests, without any ''fine-tuning'' attempts to better match the centrifuge test results. It may therefore be concluded that it is possible to infer the mechanical consequences of foundation scour through numerical simulations, provided that the geometry of the problem is properly taken into account. The latter can be estimated based on available numerical and experimental studies of the hydraulic process of local scour around bridge piers (e.g. [41,49]). Future mechanical studies should pay attention to the hydraulic conditions leading to foundation scour.
The numerical simulations have also provided new insights on the influence of foundation scour on the failure mechanisms developing in the soil around the caisson foundation. Generalized scour was shown to have a detrimental effect on vertical bearing capacity, leading to a reduction of FS v from 8.0 to 5.6, whereas the effect of local scour is almost negligible. As a consequence, the general scour of the riverbed reduces both the sidewall and the base resistance mechanisms. Conversely, local scour only affects the sides' resistance, exerted by the soil along the sides of the foundation. The latter is strongly asymmetric due to the morphology of the scour hole, which therefore plays an important role. For the specific foundation examined herein, such significant differences associated with the scour hole geometry, and therefore with the direction of loading, did not lead to any measurable effect in total moment resistance. This is because the examined cylindrical caisson is quite shallow (D f =d ¼ 1), and therefore the base resistance dominates the response. The effect of scour hole asymmetry can become significant for deeper and/or lightly loaded foundations, for which the sidewall contribution is expected to be more relevant.
Author contribution All authors contributed to the study conception and design. Material preparation, data collection, and numerical analyses were performed by Andrea Ciancimino. The first draft of the manuscript was written by Andrea Ciancimino, and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.
Funding Open access funding provided by Politecnico di Torino within the CRUI-CARE Agreement. No funding was received to assist with the preparation of this manuscript.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.
Data availability All data generated or analysed during the current study are available from the corresponding author upon reasonable request.