Shallow strip foundations subjected to earthquake-induced soil liquefaction: Validation, modelling uncertainties, and boundary effects Soil Dynamics and Earthquake Engineering

Despite recent advancements in predicting the response of shallow strip foundations during earthquake-induced liquefaction, significant modelling – related uncertainties remain, which are the focus of this paper. The problem is analysed through coupled hydromechanical analyses, employing an advanced constitutive model. The model is calibrated based only on the initial void ratio, and then validated against 6 centrifuge model tests, conducted at the University of Cambridge. Through a strict validation procedure, based on pore pressures, settlement and rotation time histories, as well as deformation mechanisms, the strengths and weaknesses of the numerical model are identified. It is shown that final settlement and rotation can be predicted with adequate accuracy, but more work is needed to achieve accurate predictions of settlement rate, maximum rotation, and pore pressures in the vicinity of the foundation. The numerical model is then used to investigate key modelling uncertainties. After revealing the sensitivity to initial soil density and to parasitic vertical acceleration, the effects of the centrifuge model container and of the distance of lateral model boundaries ( L ) are parametrically investigated. Boundary effects are minimized with a laminar container, where a normalized boundary distance L / D L ≥ 1 is shown to be adequate for all liquefiable layer depths ( D L ) examined. The use of a rigid container is proven problematic, as it always imposes an unrealistic wave propagation pattern. The use of Duxseal inclusions offers a major advantage, allowing accurate reproduction of foundation settlement even with L / D L ≥ 1, a key conclusion for the design of centrifuge tests .


Introduction
Recent contributions have led to deeper understanding of the key mechanisms controlling the response of shallow strip foundations during earthquake-induced liquefaction [1][2][3]. Despite such advancements, the accurate prediction of liquefaction-induced settlement and rotation remains challenging, due to a number of persisting uncertainties. Besides the inherent uncertainties in defining the input seismic motion and soil parameters, the problem is sensitive to a number of modelling-related uncertainties.
Numerical modelling of foundation-structure systems during earthquake-induced liquefaction requires nonlinear deformation analyses (NDAs), employing advanced fully-coupled constitutive models that can capture the nonlinear stress-strain response of liquefiable soil (e.g., Refs. [1,2,4]. Such advanced constitutive models typically contain a large number of model parameters, requiring extensive soil element testing for proper calibration. Centrifuge modelling can then be used as a benchmark for model validation.
However, validation against a single experimental result is not sufficient. Different geometrical and loading characteristics may lead to the development of different stress-paths, resulting to significantly different deformation mechanisms. The ability of a constitutive model to successfully reproduce one specific centrifuge model test cannot guarantee its success when modelling a different case, with different geometry or loading conditions. Hence, before conducting parametric studies, thorough validation against centrifuge test results of varied boundary value problems is essential.
Centrifuge model tests offer valuable insight to the studied problem and are widely accepted as validation benchmarks. However, physical modelling also has certain limitations and is not "immune" to modelling uncertainties. For instance, centrifuge-mounted shakers tend to produce parasitic input motions, which are typically ignored. Moreover, the realised material parameters (e.g., soil density ρ and saturation degree S r ) may deviate from the targeted values, a source of error that is often not considered. The soil properties in a centrifuge model are prone to spatial variability, even within the same soil layer [5,6]. Many researchers have attempted to quantify the extent of the soil's spatial non-uniformity in centrifuge model testing. Li et al. [7], developed an electrical resistivity needle probe, capable of assessing with high-resolution the spatial variability of soil electrical resistivity. They reported a coefficient of variation (COV) in terms of porosity in the range of 1-4% for loose sand, and 1-6% for dense sand. Similarly, Bolton et al. [8], examined the spatial non-uniformity of centrifuge models through in-flight cone penetration testing. They observed that the COV in terms of tip resistance, void ratio and normalized cone resistance varies between 2 and 15% in the vertical direction.
Compared to the laterally infinite soil layer of the prototype, the boundary effects created by the rigid and smooth sidewalls of a model container lead to strain and stress dissimilarities and generation of Pwaves [9]. In order to alleviate these issues, flexible-boundary containers (e.g., flexible laminar container -LC, equivalent shear beam -ESB) have been used in centrifuge modelling [9,10]. Lee et al. [11] investigated the boundary effects introduced by a laminar container on the acceleration and pore pressure generation of a saturated sand deposit. They concluded that the acceleration is not affected by the boundaries at a distance of one-twentieth of the model length from the sidewalls. The required distance for the pore pressure was 5 times larger.
Before evaluating the efficacy of numerical modelling techniques (and proceeding with parametric studies), the effect of such modelling uncertainties should be identified and quantified, which is the main scope of this paper. Initially, 6 different centrifuge model tests, conducted at the University of Cambridge [12] with variable geometries and foundation bearing pressures, are numerically analysed. The analyses are conducted employing the fully-coupled constitutive model PM4Sand [13], as implemented in the finite difference (FD) code FLAC (Itasca, 2016). The constitutive model is thoroughly calibrated based on soil element tests, conducted at the ETH Zurich (ETHZ) geotechnical laboratory [14].
The validation is not restricted to pore pressure build-up and foundation settlement and rotation, but extends to deformation mechanisms, captured in the centrifuge tests through image analysis [12]. Modelling and calibration are identical for all 6 centrifuge tests. The calibration is strictly based on the void ratio (e), without any effort to "fine-tune" model parameters to better match the centrifuge tests. This is a first key contribution of the paper, which allows in-depth assessment of sensibly calibrated numerical simulations versus a range of centrifuge experiments.
The second key contribution is the use of the thoroughly and consistently (but not perfectly) validated numerical model to investigate the effect of physical modelling uncertainties and boundary effects. The sensitivity to uncertainties related to soil density and input motion is examined, and the effect of boundary conditions imposed by different types of centrifuge containers is identified. The quantification of the influence of lateral boundaries in function of their type and distance to the structure is considered of importance to centrifuge modellers, allowing for optimized design of future experiments, but also to numerical modellers simulating such experiments for validation purposes.

Centrifuge model testing
The carefully calibrated constitutive model is employed to numerically simulate 6 centrifuge tests, that modelled the seismic response of a plane strain structure with shallow strip foundation resting on liquefiable sand. By comparing the numerical predictions to the centrifuge test results, the reliability of the constitutive model is assessed, offering insights on the sensitivity of the numerical analysis to key model parameters and setting the base for future parametric studies. The combined numerical and experimental study offers deeper understanding of the mechanisms leading to foundation settlements. Moreover, the extended validation against 6 experiments with different geometry, density, stress fields, and shear loading offers an opportunity for more thorough assessment of the strengths and weaknesses of numerical predictions.

Experimental setup
The centrifuge model tests were conducted on the Turner Beam Centrifuge of the Schofield Centre at the University of Cambridge [12]. Fig. 1a presents a typical cross-section of a centrifuge model test. In all tests, a rigid structure was supported on a shallow foundation of width B = 4.65 m (in prototype scale), having an aspect ratio h cm /B = 0.58 (where h cm is the height of the center of mass). The models were prepared by air pluviation of Hostun HN31 sand in a rigid container (equipped with a Perspex window) and tested at 50 g. As outlined Table 1, 6 centrifuge tests are considered, with different liquefiable layer depths (D L ), and different bearing pressures (q).
The fluid viscosity was scaled using a de-aired aqueous hydroxypropyl methylcellulose solution, prepared according to Adamidis & Madabhushi [15] to be 50 times more viscous than water, in order to render dynamic and pressure dissipation time scaling compatible. Images captured through the Perspex window using a high frame rate camera were used for Particle Image Velocimetry analysis [16], from which displacement fields were produced. Nevertheless, the use of a rigid container inevitably introduced boundary effects, which were alleviated using a material called Duxseal (Parker Hannifin Corporation, 2013) at the lateral model boundaries. Its effectiveness has been examined by Steedman & Madabhushi (1991), who concluded that it is able to absorb about 65% of incident wave energy for P and S waves. As illustrated in Fig. 1a, the instrumentation included miniature pore pressure transducers (PPT), piezoelectric accelerometers (ACC), micro-electromechanical system accelerometers (MEMS), and linear variable displacement transducers (LVDT). More details on the experimental setup can be found in Adamidis [17].

Numerical modelling of centrifuge tests
Fully-coupled, effective stress, time history numerical analyses are performed in the FD code FLAC 8.0 (Itasca 2016) employing the PM4sand (Version 3.1) constitutive model (Boulanger & Ziotopoulou [13]. Fig. 1b depicts the FD mesh, model geometry, interface and water table location. Duxseal is modelled as linear elastic, with its permeability set 100 times lower than that of the adjacent soil. Following Popescu & Prevost [18], the Duxseal Young's modulus is set to E = 800 kPa, its Poisson's ratio ν = 0.46, and its mass density ρ = 1.65 Mg/m 3 . The boundary conditions imposed by the rigid container are accounted for by applying the seismic excitation at the bottom and edge nodes of the FD model, using a velocity time history that corresponds to the acceleration trace that was measured at the base of the centrifuge container. The bottom nodes are fixed in the vertical direction, while the side nodes are free. The input motions for the 6 centrifuge model tests are presented in Fig. 1c. In all cases considered, the width and height of the structure are 4.65 m and 5.36 m, respectively. They are modelled using linear elastic beam elements of mass density ρ, Young's modulus E, cross section area A, and moment of inertia I. A fine assembly of beam elements is adopted in this paper, which leads to a more accurate prediction of the structural rotation than the coarser beam assembly that was used for the numerical modelling of the OA06 centrifuge experiment in Kassas et al. [14]. The properties of the beam elements are adjusted to match the bearing pressure of the respective centrifuge test (Table 1). In the centrifuge experiments, the fixed-base frequency of the structure was selected to be higher than that of the input motion, promoting a rigid-body type of response [12]. This is maintained in the numerical analyses, setting a large enough stiffness to the structure to promote rigid-body response. A frictional soil-foundation interface is introduced, allowing for sliding and separation. The foundation width is discretized using 20 elements, while the friction angle at the interface is set to δ = 28.8 ο , following measurements using the direct shear apparatus [17].
Following the recommendations of Boulanger & Ziotopoulou [13], Rayleigh damping is set to 0.5% at a center frequency of 1 Hz, which is the dominant frequency of the input motion. The PostShake Flag parameter of PM4sand is enabled after strong shaking (20 s) to improve modeling of post-liquefaction reconsolidation strains, as suggested by Ziotopoulou [19]. Large deformations and second order effects are accounted for, by automatically updating mesh node coordinates after each analysis step.

Constitutive model calibration
The seismic response of fully saturated Hostun (HN31) sand is simulated using the plane-strain, critical stress compatible, bounding surface plasticity model PM4Sand, Version 3.1 [13]. The model has been developed for geotechnical earthquake engineering applications, following the framework of the DM04 model [20]. PM4sand has been shown to be capable of predicting soil behavior under cyclic loading at different relative densities, under different confining pressures and stress paths. A detailed description of model formulation and its modifications can be found in Ziotopoulou & Boulanger [21] and Ziotopoulou & Boulanger [19].
Accurate modelling of cyclic soil response is key for the successful simulation of centrifuge model tests. PM4Sand comprises 27 input parameters, of which 6 are considered primary, 2 are defined as "flag" parameters, and 21 as secondary parameters. The primary parameters include the apparent relative density (D r ), the small-strain shear modulus coefficient (G o ) that controls the small strain shear modulus (G max ), and the contraction rate parameter (h p0 ), which can be adjusted to obtain the cyclic resistance ratio (CRR). Preset values can be used for the secondary parameters, or they can be calibrated using experimental data. The calibration procedure followed here is based on an extensive element testing database, with the full process reported in Kassas et al. [14]. The key calibration parameter for each centrifuge model test is the initial void ratio (e). Table 2 summarizes the calibration methodology for Hostun (HN31) sand, considering the void ratio of the centrifuge model test that is to be numerically simulated. The small-strain shear modulus coefficient (G o ) is calculated for the desired void ratio using the analytical expression of Azeiteiro [22] for the small strain shear modulus of Hostun (HN31) sand (G max ), validated against bender element tests conducted at ETHZ [14]: The contraction rate parameter (h po ) ( Table 3) is adjusted to obtain the target cyclic resistance ratio (CRR). Bolton [23] critical state constants are estimated as Q = 8.4 and R = 0.78 based on triaxial test results. The permeability (k) is calculated using the analytical expression where: C CK is an empirical coefficient; and S the ratio of surface to soil particles' volume. The cyclic stress resistance curve for each centrifuge model test is calculated as a function of the initial void ratio (e) and number of cycles to liquefaction (N): This equation was produced from constant volume simple shear tests of Hostun HN31 sand samples with initial vertical effective stress σ ′ vo = 100 kPa [14]. Fig. 2 presents the cyclic stress ratio versus number of cycles required for liquefaction for the void ratio corresponding to each of the six simulated centrifuge model tests, as produced using Eq. (3).

Assessment of numerical simulations
The numerical prediction is compared to the centrifuge test results in terms of settlement and rotation of the structures, pore pressure buildup, as well as the deformation mechanisms that develop within the liquefiable layer. Both the calibration methodology and the simulation technique are exactly the same for all cases examined. The numerical methodology is assessed for three different liquefiable layer depths (D L = 0.5, 1.0, 2.5 ×B) and two levels of bearing pressure(q = 50, 100 kPa).

Pore water pressure build-up
Figs. 3 and 4 present the comparison between the computed pore water pressure time histories and those measured in the centrifuge experiments, at the free-field and below the structure, respectively. The comparison is presented for all 6 centrifuge tests, corresponding to two levels of bearing pressure (q = 50, 100 kPa) and three liquefiable layer depths = 0.5, 1.0, 2.5 × B). For all cases examined, the numerical analysis predicts adequately well the developed pore-water pressures at the free-field (Fig. 3). Both the initial build-up and the maximum values of pore water pressures are predicted satisfactorily, irrespective of D L . In the experimental time-histories, there is more significant co-seismic excess pore water pressure dissipation than what is predicted by the numerical analysis.
The numerical predictions of pore water pressure build-up under the structure are less accurate (Fig. 4). The numerical analysis captures adequately well the pore water pressures directly below the footing (P3) only for the deep (D L /B = 2.5) layer. At this location, the initial (for t < 5s) drop of pore water pressure is due to dilation and is wellpredicted for both bearing pressures examined (q = 50kPa, OA04; and q = 100kPa, OA05). The discrepancies between analysis and experiment are more pronounced for deeper locations (P2 and P1). With the progression of seismic shaking (t ≥ 5s), the excess pore water pressure is partly dissipated in the experiments, a phenomenon that is not captured by the numerical analysis. The numerical predictions are less accurate for intermediate (D L /B = 1.0, OA06 and OA07) and shallow (D L /B = 0.5, OA08 and OA10) layers, both in terms of initial pore water pressure build-up and to its subsequent, co-seismic dissipation.
The inferior prediction of the initial abrupt increase of pore water pressures is attributed to the lower excess pore pressures predicted by the constitutive model during the initial load cycles [14]. As discussed later on, this discrepancy is also reflected in the reduced accumulation of settlement during the initial load cycles (Figs. 5 and 9). The discrepancies with respect to co-seismic pore water pressure dissipation (observed for all cases examined) are attributed to the inherent inability of the constitutive model to accurately predict plastic volumetric strains [21]. Another potential source of error is the modelling of the potentially complicated fluid flow at the soil-foundation interface. In the numerical simulations, the potential for fluid flow along the interface and the consequent localized material transfer were not taken into account.     response to the experimentally measured. With the exception of the lightly-loaded (q = 50kPa) structure lying on a shallow (D L /B = 0.5) layer (test OA08) and the heavily-loaded (q = 100kPa) structure lying on the intermediate deep (D L /B = 1.0)layer (test OA07), the final settlement is predicted sufficiently well. With the exception of the shallow layers (D L /B = 0.5), the numerical predictions are reasonably conservative, over-estimating (not to a large extent) the final foundation settlement. Looking closer at the evolution of settlements during shaking, the discrepancies are more pronounced. The analysis significantly underpredicts the settlement during the first load cycle, while a slight overprediction is observed for later load cycles. The under-prediction during the first cycles of loading can be linked to the inability of the constitutive model to capture the abrupt drop of initial effective stresses [14]. A related phenomenon is the sudden reorientation of sand grains during initiation of shear loading, which is difficult to be reproduced in the framework of continuum mechanics. The potential spatial variability of soil properties in the centrifuge experiment and the imperfect contact between the footing and the ground surface are not accounted for in the numerical model, where a uniform and perfectly symmetric soil layer is assumed. The over prediction of co-seismic settlements is attributed to the lower cyclic stress resistance predicted by the model after the first ten load cycles, compared to the corresponding soil element tests. The latter is related to the calibration of the constitutive model, in which it was not possible to capture the targeted shear stress resistance over a wide range of loading cycles. This led to over-prediction of the cyclic shear resistance before the first ten load cycles, and to its under-prediction in the following load cycles [14].

Settlement-rotation (wθ) response
The evolution of rotation is also predicted with non-negligible discrepancies. With the exception of test OA04, the rotational oscillation amplitude was experimentally observed to increase with subsequent load cycles. Although the numerical analysis reproduced the trend for test OA04, higher levels of rotational oscillation than those experi-mentally recorded were predicted for all deep and intermediate layers examined. This discrepancy was reversed for the shallow layers (D L /B = 0.5), where the experimentally recorded rotational amplitude was higher than the numerical prediction. The discrepancies for the deep and intermediate layers can be attributed to the previously discussed underprediction of excess pore water pressures below the structure in the beginning of shaking. In the experiments, the larger excess pore water pressures developing below the foundation led to more pronounced softening of the soil beneath, limiting the acceleration transmitted to the structure, consequently leading to lower magnitude of rotation at the beginning of the motion. The higher amplitude of rocking oscillations that is predicted later can be attributed to the difference in the deformation mechanism between physical and numerical modelling, as shown in Kassas et al. [14].
Furthermore, the imperfect soil-foundation interface and the progressive embedment of the foundation, which was observed in the centrifuge experiments, but cannot be accurately modelled in the numerical simulation may lead to some discrepancies. At this point, it should be pointed out that in some centrifuge experiments (OA07, OA10) the structure had accumulated a small amount of rotation before the seismic shaking was imposed. This initial tilt is not included in the numerical simulations.
The numerically predicted final settlements are compared to the experimentally measured for all cases examined in Fig. 6 (left), along with settlement estimations according to simplified methodologies. Table 4 presents the required parameters for the aforementioned simplified prediction methods. Given the very thorough calibration process and the higher level of sophistication, it is not surprising that the numerical analysis offers the best results. However, certain simplified predictions perform better than others. The methods of Tokimatsu & Seed [24] and, often used in practice, consistently under-predict (except one case) the measured settlement. This is also no surprise, as these methods were developed for prediction of free-field settlements, capturing the volumetric contraction of soil. They cannot accurately predict structural settlement, which is primarily due to shear mechanisms [12,25].
The method of Bertalot et al. [26], developed for the estimation of an upper-bound of liquefaction-induced settlement of structures on shallow foundations, performed much better, indeed offering an upper bound estimate of settlement for all cases examined.
In Fig. 6 (right), the numerically predicted residual foundation rotations are compared to the experimentally measured for all cases examined. Additionally, the rotation estimation according to the simplified methodology of Adamidis & Madabhushi [3] is presented, which offers predictions of maximum rotation and conservative estimations of residual rotation for D L /B ≤ 1. The sophisticated numerical simulation does not seem to outperform the simplified method. Fig. 8 and compare the numerical predictions to the experimental results in terms of contours of total shear strain at the end of shaking, for the centrifuge tests with bearing pressure q = 50 and 100 kPa respectively. Overall, the numerical simulation successfully predicts the concentration of shear strains below the foundation edges and in the adjacent soil towards the free field. With the exception of experiment OA06 (Fig. 7), the shearing observed directly below the foundation is also predicted reasonably well. Particularly for experiment OA04, the numerical simulation captures the limited amplitude of shearing bellow the foundation (Fig. 7), due to the reduced amplitude of rocking oscillation, that is also predicted adequately well (Fig. 5).

Deformation mechanisms
The reliability of the numerically calculated deformation mechanisms depends on the depth of the liquefiable layer. For intermediate (D L /B = 1.0) and shallow (D L /B = 0.5) layers, significant shear strains are predicted deep below the foundation, close to the base of the layer, which were not observed in the centrifuge model tests. This is related to the numerical prediction of increased excess pore water pressures in these locations, leading to more significant soil softening and larger shear strains. On the contrary, in the case of the deep layer (D L /B = 2.5), where the pore water pressure generation below the foundation is captured more adequately, the calculated mechanism is similar to the one experimentally observed.
In the deep layer (D L /B = 2.5), boundary effects are significant, as evidenced by the shear strains numerically calculated at the model edges, in the immediate vicinity of the Duxseal inclusions. This behavior could not be directly observed in the centrifuge experiments due to the limited area monitored by the camera, which did not encompass the entire model. However, Adamidis & Madabhushi [12] derived the same conclusion, based on the magnitude of displacements that reached the boundaries of the monitored soil.
Particularly interesting is the depth of the deformation mechanism predicted in the case of the deep layer (D L /B = 2.5). The mechanism computed for experiment OA05 is more shallow than for experiment OA04. Adamidis & Madabhushi [12] also reported a shallower mechanism for the heavily-loaded structure (OA05), though it seems that the numerical prediction is even shallower, as implied by the pore water pressure prediction at P1 (Fig. 4). The choice of a rigid centrifuge container, which despite the Duxseal inclusions did not allow the lower part of the soil layer to shear sufficiently, in combination with the localization of strains due to the higher amplitude of rocking oscillation in experiment OA05 (Fig. 5), are responsible for the prediction of a more shallow mechanism.
Despite the aforementioned discrepancies, the numerical simulations reproduce successfully some of the key aspects of the problem. Both the final rotation and settlement are predicted with reasonable accuracy. Unsurprisingly, the numerical simulation outperforms all of the examined simplified methods for settlement prediction (Fig. 6). However, the rate of settlement is not predicted with great success (Figs. 3 and 10), and the maximum rotation is over-predicted for most of the cases examined. The excess pore water pressures are accurately predicted in the free field (Fig. 3), but less successfully below the foundation (Fig. 4), especially for shallow layers. Finally, the deformation mechanisms are reasonably reproduced and the boundary effects in the case of deep layers are rightly identified (Fig. 7 and Fig. 8 ).

Physical modelling uncertainties
Model validation against centrifuge tests requires careful and consistent calibration. The previously discussed validation against 6 centrifuge tests was performed consistently, without any "fine-tuning" of calibration parameters to improve predictions. The validation is not perfect, and the key factors leading to the observed discrepancies were identified and discussed. However, in addition to the limitations of any constitutive model, there are also physical modelling uncertainties that may affect the results of such a validation process. These uncertainties are related to material properties, input motion characteristics, and boundary conditions. This section explores the sensitivity of the results to such uncertainties.

Sand relative density
The element tests for the characterization of Hostun HN31 sand and the constitutive model calibration were conducted at ETHZ [14], with the sand acquired from the same provider as for the University of Cambridge centrifuge model tests [17]. Fig. 9 illustrates the limiting void ratios (e min , e max ) of the studied sand, according to 5 different published studies [22,[31][32][33]36], and compares them to the ones measured at ETHZ according to ASTM D4253-00 [28] and ASTM D4254-00 [29]. Regarding e min , it was found that a denser state could be achieved by air pluviation through sieves.
As shown in Fig. 9, the reported limiting void ratios exhibit a considerable range. For example, in the case of experiment OA05 (e = 0.838), differences in estimation of e min and e max lead to a range of estimated relative density D r from 38% to 56% -a difference of 19%. This is critical in numerical modelling with the PM4sand constitutive model, where relative density (D r ), rather than void ratio (e), is used as an input parameter. Furthermore, considering that an error in the estimation of D r is also possible, in particular in the centrifuge model tests but also in the element tests, a sensitivity analysis based on the initial relative density of each model test is conducted, assuming two cases with initial D r ranging ± 4 %. Fig. 10 shows that the prediction of settlement (w) is very sensitive to the assumed relative density D r , with the ± 4% range examined leading to significant deviations. Interestingly, the settlement response of test OA06 appears to be more sensitive than that of the remaining tests. This can be attributed to the fact that the cyclic strength resistance predicted by PM4sand is not a linear function of the initial D r . In this specific test (e = 0.788), the initial D r is higher in comparison to other tests, and hence the sensitivity to the predicted soil cyclic shear strength is also higher.

Duxseal boundaries
Another significant source of uncertainty are the properties of the Duxseal material. In the absence of more detailed characterization, the Duxseal material was modelled as linear elastic, as proposed Popescu & Prevost [18]. However, the assumption of linear elastic response is not necessarily accurate. Especially at larger depths, where the developing lateral stresses transmitted to Duxseal are higher, its response may very well be nonlinear, leading to non-negligible boundary effects. To explore the uncertainty in the properties of Duxseal as well as the effect of potential nonlinearity, a sensitivity analysis is conducted, varying the elasticity modulus of the Duxseal material (E DUX ) from half (400 kPa) to double (1600 kPa) of the initial estimate. (800 kPa). Fig. 11a summarizes the results of the sensitivity analysis in terms of the final co-seismic settlement. As one might expect, the response is more significantly affected for the deeper D L /B = 2.5layers (OA04 and OA05), where larger lateral stresses are transmitted to the Duxseal boundaries. The sensitivity to E DUX practically disappears for the shallow D L /B = 0.5 layers (OA08 and OA10). In addition to the smaller lateral stresses, in these cases the deformation mechanisms did not reach the Duxseal boundaries (Fig. 8 and ). Interestingly, the results indicate the presence of boundary effects, albeit limited, also for the intermediate D L /B = 1.0 layers (OA06 and OA07).
It is particularly important that the stiffness of the Duxseal inclusions affects the developed deformation mechanism, as shown in Fig. 11b. The softer material allows the soil to shear at greater depth, leading to a deeper mechanism, more similar to the one observed in the corresponding centrifuge experiment (OA05).

Input motion
In addition to the targeted horizontal input motion, (centrifuge) seismic shakers also tend to produce a parasitic vertical motion, which is  the subject of this section. The amplitude of the vertical acceleration is significantly lower than the horizontal, and is thus often neglected in numerical analysis. In order to quantify the sensitivity of the results to such parasitic vertical acceleration, the numerical analyses are repeated taking it into account.
In the centrifuge experiments, the horizontal and vertical acceleration of the model container was recorded only at its bottom left corner (Fig. 1a). Based on this recording, two cases are considered: (a) Fig. 10. Sensitivity of predicted settlement (w) to the assumed relative density (D r ): comparison of numerical predictions to experimental measurements for the 6 centrifuge model tests. Horizontal (V H ) and vertical (V V ) input velocity applied throughout the rigid box; and (b) Horizontal velocity (V H ) applied throughout the rigid box and angular velocity (ω) around the midpoint of the base of the rigid box, resulting in the recorded vertical motion at the bottom left corner. The latter case aims at simulating a parasitic yawing motion introduced by the shaker. Such a response is considered as it has been observed in certain centrifuge shakers [30].
In Fig. 12a, the effect of the vertical velocity component is presented for the two shallow layer (D L /B = 0.5) tests OA08 (q = 50kPa) and OA10 (q = 100kPa). Fig. 12b depicts the velocity time histories, as computed from the recorded accelerations. For both tests, the addition of the rotational velocity (ω) does not have any significant influence, leading to slightly reduced final settlement (w). On the other hand, the scenario of uniform parasitic vertical velocity leads to larger settlements for both cases examined, despite the relatively low amplitude of the vertical velocity (V V ) in comparison to the horizontal one(V H ). Especially for the lightly loaded foundation (q = 50kPa) of test OA08, the addition of the parasitic vertical velocity leads to a substantial increase of settlement (w). A similar sensitivity was observed in Fig. 10, where the settlement of this test also increased significantly with a small reduction of relative density D R (by 4 %), unlike the more heavily loaded (but of same depth) test OA10. It may be inferred that in the simulation of experiment OA08, a relatively small increase in loading or decrease in strength is sufficient to produce significantly larger settlements, bringing the numerical prediction much closer to the experimentally measured value.

The effect of centrifuge model container
Despite its unavoidable limitations, the (thoroughly, but not perfectly) validated numerical model can be employed to extend the knowledge obtained from the centrifuge model tests, and to explore the effects of physical modeling uncertainties. Such insights may allow optimization of future centrifuge tests.
This section investigates one of the key uncertainties in physical modelling: the effects of the type of centrifuge model container. As previously discussed, instead of a flexible laminar container (LC) or an equivalent shear beam (ESB), a rigid container was used to allow for recording of deformation fields and the associated deformation mechanisms through PIV analysis. However, a rigid container unavoidably introduces boundary effects [9]. Duxseal material was added at the edges of the container, in an effort to absorb some of the incident wave energy. However, especially for the deeper layers tested (OA04 and OA05), the deformation mechanisms appear to have reached the model boundaries, revealing the importance of their lateral distance (which is explored in the next section).
In order to investigate the effect of the model container type, the numerical model is employed to simulate the response of an idealizedadequately longlaminar container, where free-field response (i.e., away from the structure) is not affected by the boundaries. To achieve this condition in simulation, the lateral model boundaries are now moved further away from the structure: 4.0B in simulation, instead of 1.8B in the centrifuge model tests. The idealized laminar container is simulated by connecting opposite lateral boundary nodes of the same depth. Fig. 13 compares the settlements of the 6 centrifuge model tests using the original rigid container with Duxseal inclusions (termed Duxseal boundary from now on), to those obtained using the idealized laminar container (LC). In all cases examined, the shape of the settlement curve is not significantly affected. With the exception of test OA05, the Duxseal boundary only leads to a slight increase of the final settlement compared to the idealized LC, which is considered here as a benchmark. This indicates that boundary effects associated with the use of a rigid container with Duxseal inclusions lead to conservative results in terms of accumulated settlement. The influence of the Duxseal inclusions and the distance of the lateral boundaries to the foundation are examined in more detail in the following section.
To gain more insight on boundary effects, the deformation mechanisms are examined in more detail. Fig. 14 presents the total computed shear strain at the end of shaking for the heavily loaded structures, comparing the rigid container with the Duxseal boundaries to the idealized laminar box. The deformation mechanisms do not appear to be affected by the boundary modification, either for the shallow (OA10) or the intermediate layer (OA07), as the Duxseal inclusions were far enough in both cases. In stark contrast, the deformation mechanism is significantly different for the deep layer (OA05). In this case, the shearing zone below the foundation clearly reaches the Duxseal inclusions, propagating along the Duxsealsoil interface towards the ground surface, developing a preferential lower energy mechanism. In the case of the idealized laminar container (LC), a deeper deformation mechanism develops, with increased shear strains below the edges of the foundation, all the way to the base of the liquefiable layer. This deeper mechanism cannot develop in the case of the rigid container with Duxseal inclusions, as the shear band is interrupted by the lateral model boundaries. In the case of the deep layer experiments (OA04 and OA05), the interference of lateral model boundaries and the resulting interception of the deeper deformation mechanism, is also reflected in the pore-water pressure build-up. Fig. 15 compares the pore-water pressures below the foundation (left) and in the free-field (right) for test OA05. Interestingly, boundary effects are practically negligible in the free field (P2, P4, P14) and close to the foundation (P3, P12). The only visible difference is observed at the deepest point of measurement below the foundation  (P1). In the case of the Duxseal inclusions, the soil in the free field has already been fully liquefied (r u = 1). This is not the case for the soil bellow the foundation, due to the increased vertical load imposed by the superstructure and the more shallow failure mechanism. In the case of the idealized laminar container, the development of the deeper deformation mechanism leads to higher levels of shearing in the soil below the foundation (P1), resulting in increased excess pore-water pressure generation.

Parametric study: container type and distance of lateral boundaries
In order to obtain deeper understanding of boundary effects, in function of centrifuge model container and distance of the lateral boundaries (L) to the foundation, a parametric study is conducted. In all cases examined, the container itself is not simulated, but imposed by appropriate boundary conditions. In the case of the laminar container, the horizontal input seismic motion is applied directly to the base nodes of the numerical model, and the previously discussed node connectivity is employed to model the laminates in a simplified manner (i.e., ignoring the small but potentially non-negligible friction between successive laminates, and their inertia). For the other two rigid container types, the input motion is applied concurrently at the base and at the lateral nodes of the numerical model. The Duxseal inclusions are modelled as previously described.
Three response variables are investigated: (i) the final co-seismic settlement (w); (ii) the maximum structural acceleration (a); and (iii) the maximum co-seismic rotation (ϑ). Each response variable is normalized by its corresponding value for the idealized laminar container (L = 20D L ), considered here as the best approximation of realityin a numerical analysis, this would be the boundary of choice, if computational effort was not an issue. In centrifuge modelling, centrifuge capacity and space limitations render the use of such a long LC impractical. The normalized response variables are plotted versus the distance of the lateral boundaries (L), normalized by the liquefiable layer depth (D L ) corresponding to each analysis. Fig. 16a compares the final co-seismic settlement (w /w ILC ) normalized by that of the idealized LC, for the three model containers and three liquefiable layer depths (D L ), in function of the distance of the lateral boundaries (L/D L ) normalized by the liquefiable layer depth. The laminar container offers the best response, as w converges for L ≥ 1 D L for all liquefiable layer depths examined. This is an important conclusion for centrifuge model test campaigns, allowing the examination of larger foundations at lower g levels, without introducing boundary effects in the settlement response.

Settlement (w)
In the case of a rigid boundary with Duxseal inclusions, a distance of L ≥ 6 D L is required for the settlement not to be affected by the model boundaries. Although the low shear stiffness of the Duxseal material allows some soil deformation, it cannot compete with the effectiveness of a laminar container. In the case of the shallow and the intermediate layer, the fact that w/w ILC > 1 for L/D L = 1is attributed to the Duxseal being located within the developing deformation mechanism. It is worth noting, that the distance to the lateral boundaries was shorter than needed for convergence of settlement for all of the examined centrifuge model tests (Table 1). Nevertheless, thanks to the Duxseal inclusions, the deviations of the normalized settlement are not expected to have exceeded 10%.
With a fully rigid container (without Duxseal), the normalized settlement w/w ILC converges to 1 (i.e. it is not affected by boundary effects anymore) for a boundary distance of L ≥ 6D L . However, the initial deviation of w/w ILC from 1 (before convergence) is significantly higher Fig. 15. The effect of centrifuge model container in terms of pore pressure build-up for the deep layer test OA05 (q = 100 kPa). Comparison of computed pore water pressures below the foundation (left) and in the free-field (right) for the original rigid container with Duxseal inclusions, to an idealized-adequately long-laminar container (considered here as the benchmark).
compared to the case with Duxseal inclusions. The decrease of L/ D L (i.e., the boundaries getting closer to the foundation) below 6, leads to a significant decrease of the normalized settlement w/ w ILC . In the absence of the deformable Duxseal inclusions, the rigid container restricts soil shearing, even hindering the occurrence of liquefaction. Evidently, the closer the boundaries to the point of interest, the bigger their effect. Nevertheless, it is worth noting that the settlement can be accurately modelled in a centrifuge even with a rigid container, provided that the lateral boundaries are kept far enough from the foundation (L ≥ 6D L for the cases examined).

Maximum acceleration (a)
Regarding the maximum structural acceleration (a), convergence is observed for a similar distance L/D L ≥ 6 for all container types examined (Fig. 16b). The results are similar for the laminar container and the one with Duxseal inclusions, with maximum deviations of up to 20% for L/D L < 6 (i.e., before convergence takes place). The deviations are more significant in the case of the rigid container (without Duxseal). Interestingly, these large deviations are only observed for the intermediate layer (D L /B = 1), indicating resonance as the main cause. Indeed, the first natural frequency (F n,1 ) of the D L /B = 1 soil column is 6.8 Hz (without considering liquefaction), not far from the fixed based frequency of the modelled structure (7.0 Hz). The first natural frequencies of the other two soil columns with D L /B = 0.5 and 2.5 are 11.4 Hz and 3.4 Hz, respectively. Fig. 16c illustrates the evolution of maximum normalized rotation (ϑ/ϑ ILC ) in function of the normalized distance of the lateral boundaries (L/ D L ). As shown in Fig. 5, the numerical model is not able to fully predict the rotational response of the centrifuge model tests. This discrepancy is most likely due to the limitations of the constitutive model and the FD code [14]. Despite this limitation, the influence of the boundary conditions is still worth investigating numerically, as it can offer insight for the design of future centrifuge model tests. In the case of the laminar container, the maximum rotation ϑ/ϑ ILC starts fluctuating around a value of 1 forL/D L ≥ 6. Below that threshold, ϑ/ϑ ILC decreases below 1, as the boundaries get closer. A similar value of for convergence (L/D L ≥ 6) is required for the rigid box with Duxseal inclusions. However, deviations from 1 are larger than for the LC, for small distances of the lateral boundaries.

Maximum rotation (ϑ)
The effects of model boundaries on rotation are more pronounced for a fully rigid container. The scatter in maximum rotation is significantly higher and no clear threshold for convergence to ϑ / ϑ ILC = 1 can be detected. For example, in the case of the intermediate D L / B = 1 layer, the normalized maximum rotation ϑ / ϑ ILC ≈ 1.6 for L/ D L = 4, reducing to 0.9 for L/D L = 5, and increasing again to 1.2 for L/ D L = 6. Fig. 17 compares the settlement-rotation (w − θ) response for the three model containers, progressively increasing L/D L from 4 to 6. Already at L/ D L = 4, the rotation is practically unaffected by the boundaries of the laminar container. Though not equally effective, the Duxseal container exhibits minimal boundary effects on the w − θ response, also for L/ D L = 4. In the fully rigid container, the rotational response is significantly (and erratically) influenced by the model boundaries for all L/ D L values examined, with no obvious trend for convergence.
In an effort to gain deeper insight on the response of the three model containers, Fig. 18 compares the developing wave fields at the ground surface for D L /B = 1, in function of model container type and varying L/ D L from 4 to 5. For each model container type, the wave field is constructed by plotting together the acceleration time histories of all the nodes at the ground surface. To allow direct comparison, both L/ D L = 4 and 5 refer to the same area of soil (to the right of the structure), with L/D L = 5 plot flipped. The vertical axis depicts the normalized distance x/D L , measured from the center of the foundation. The empty area from L/D L = − 0.5 to 0.5 corresponds to the nodes under the foundation, which are not plotted for better comparison.
In the case of the laminar container (Fig. 18a), a shear wave is detected at the edge of the foundation for t = 6.46 s, arriving to the boundary with a time lag. At this time, the pore-water pressure at the far-field is higher than below the foundation. Hence, the effective shear modulus (G) and the corresponding shear wave velocity (V s ) at the farfield are lower than below the foundation, leading to the observed time lag. Of particular importance is that the shear wave arrives at the same time at the edge of the foundation for L/D L = 4 and 5. For t < 6.46 s, the time lag between waves arriving at the foundation and at the far-filed is much less pronounced, as the pore water pressures are still limited. Fig. 18c depicts a very different response for the rigid container. In this case, the shear wave is first detected at the model boundary, being transmitted towards the foundation and arriving there with a time lag. Due to this inverted response (compared to the laminar box), the arrival time at the edge of the foundation depends on the distance of the lateral boundaries (L). Indeed, the shear wave arrives at the foundation at t = 6.30 s for L/D L = 4, and a bit later, at t = 6.43s for L/D L = 5. This difference in wave arrival may be the culprit for the larger scatter in maximum structural rotation observed for the rigid box in Figs. 16 and 17. Another key implication for rigid containers, is that the acceleration at the ground surface is not reduced after the occurrence of liquefaction, as the container walls keep transmitting accelerations at the lateral edges of the model. Fig. 18b confirms the improvement in rigid container response offered by Duxseal inclusions. In this case, the acceleration wave arrives at the edge of the foundation through the soil column below the structure. The arrival is the same for both L/D L = 4 and 5, similarly to what was observed for the laminar container. Afterwards, the acceleration wave is transmitted away from the foundation, towards the model boundaries. However, the Duxseal does not manage to completely filter out the acceleration transmitted from the rigid container, as indicated by the arrival of acceleration waves close to the edges of the model. Despite this limitation, the addition of Duxseal undoubtedly significantly improves the performance of a rigid container.

Conclusions
This paper has studied the seismic response of shallow strip foundations during earthquake-induced liquefaction. Coupled hydromechanical analyses were conducted employing the FD code FLAC2d, modelling nonlinear soil response with PM4Sand. The model was calibrated according to Kassas et al. [14], based only on the initial void ratio of the soil and was then validated against 6 centrifuge model tests of different geometry and foundation bearing pressure. The centrifuge tests used herein for validation were previously conducted at the University of Cambridge [12].
Through a strict validation procedure, based on pore pressures, foundation settlement and rotation time histories, as well as deformation mechanisms (extracted in the centrifuge tests through image analysis), the strengths and weaknesses of the numerical modelling technique were identified. The final settlement and rotation were systematically predicted with adequate accuracy, but more work seems to be necessary to achieve accurate prediction of settlement rate, maximum co-seismic rotation, and the evolution of pore water pressures in the vicinity of the foundation.
The thoroughly (yet not perfectly) validated numerical model was then used to investigate key centrifuge modelling uncertainties. The results were shown to be highly sensitive to the initial density of the liquefiable soil layer, with small changes leading to significant fluctuations in the predicted settlement. The uncertainties regarding the stiffness of Duxseal inclusions at the edges of the physical model were found to influence the results only for deep layers. However, the sensitivity to Duxseal's properties highlights the need for more detailed modelling of this material. Furthermore, the response of Duxseal needs to be taken into consideration for the interpretation of a centrifuge experiment's result. The parasitic vertical component of the input motion (simulated by the centrifuge-mounted shaker) was also shown to affect the response for deep layers.
Further investigation focused on the type of centrifuge model container and the distance of lateral boundaries. The rigid container with Duxseal inclusions that was used in the examined centrifuged model tests was compared to an idealized laminar container (ILC), with the lateral boundaries at an adequately large distance allowing the reproduction of free-field conditions. Using the ILC as benchmark, it was shown that the use of the rigid container with Duxseal inculsions led to a slight over-estimation of settlement. Moreover, in deep layers, the proximity of the Duxseal boundaries led to a more shallow deformation mechanism than the one numerically predicted for the ILC benchmark.
A parametric study followed, which allowed more general conclusions to be drawn on boundary effects. Three types of containers were comparatively assessed: a laminar container (LC), a rigid container, and a rigid container with Duxseal inclusions. For all container types, the effect of the distance (L) of the lateral boundaries normalized by the liquefiable layer depth (D L ) on settlement and structural acceleration prediction was parametrically investigated. As might be expected, boundary effects are minimized for a laminar container. A normalized boundary distance L/D L ≥ 1 was shown to be adequate for all liquefiable layer depths examinedan important conclusion for the design of centrifuge tests.
The use of a rigid container was shown to be problematic, as it always imposes (through its rigid lateral boundaries) an unrealistic wave propagation pattern. Nevertheless, the settlement can still be reproduced with reasonable accuracy, provided that the lateral boundaries are at an adequately large distance of L/D L ≥ 6. Evidently, this poses a severe practical restriction for centrifuge testing, leading to a requirement for a very long rigid container, in order to achieve a realistic response. The use of Duxseal inclusions offers a major advantage, allowing adequate reproduction of foundation settlement even with L/ D L ≥ 1 -a key conclusion for the design of centrifuge tests. However, a larger distance of L/D L ≥ 6 is still needed to correctly reproduce structural acceleration, though potential errors are much smaller than for a fully rigid container.