Spatial variability of shear wave velocity: implications for the liquefaction response of a case study from the 2010 Maule Mw 8.8 Earthquake, Chile

the effect of

Assessing the potential and extent of earthquake-induced liquefaction is paramount for seismic hazard assessment, for the large ground deformations it causes can result in severe damage to infrastructure and pose a threat to human lives, as evidenced by many contemporary and historical case studies in various tectonic settings.In that regard, numerical modeling of case studies, using state-of-the-art soil constitutive models and numerical frameworks, has proven to be a tailored methodology for liquefaction assessment.Indeed, these simulations allow for the dynamic response of liquefiable soils in terms of effective stresses, large strains, and ground displacements to be captured in a consistent manner with experimental and in-situ observations.Additionally, the impact of soil properties spatial variability in liquefaction response can be assessed, because the system response to waves propagating are naturally incorporated within the model.Considering that, we highlight that the effect of shear-wave velocity V s spatial variability has not been thoroughly assessed.In a case study in Metropolitan Concepción, Chile, our research addresses the influence of V s spatial variability on the dynamic response to liquefaction.At the study site, the 2010 Maule M w 8.8 megathrust Earthquake triggered liquefaction-induced damage in the form of ground cracking, soil ejecta, and building settlements.Using simulated 2D V s profiles generated from real 1D profiles retrieved with ambient noise methods, along with a PressureDependentMultiYield03 sand constitutive model, we studied the effect of V s spatial variability on pore pressure generation, vertical settlements, and shear and volumetric strains by performing effective stress site response analyses.Our findings indicate that increased V s variability reduces the median settlements and strains for soil units that exhibit liquefaction-like responses.On the other hand, no significant changes in the dynamic response are observed in soil units that exhibit non-liquefaction behavior, implying that the triggering of liquefaction is not influenced by spatial variability in V s .We infer that when liquefaction-like behavior is triggered, an increase of the damping at the shallowest part of the soil domain might be the explanation for the decrease

Introduction
Earthquake-induced liquefaction occurs when granular soils, frequently loose to medium-dense saturated sands, exhibit a significant loss of strength and stiffness and start behaving as a liquid due to the seismic generation of excess pore-water pressure (National Academies of Sciences, E. Medicine, 2016).The large ground deformations and damages resulting from liquefaction have been well-documented in recent case histories.Notable examples include the 2010 Maule Chile earthquake M w 8.8 (Bray et al., 2012;Verdugo and González, 2015), the 2010-2011 Canterbury earthquake sequence in New Zealand (maximum magnitude M w 7. 1 Green et al., 2014;Van Ballegooy et al., 2014), the 2012 Emilia Earthquake sequence in Northern Italy (maximum magnitude M w 6. 1, Emergeo Working Group, 2013), and the recent Kahramanmaras earthquake sequence (maximum magnitude M w 7.8, Baser et al., 2023;Taftsoglou et al., 2023).Although these earthquakes varied in magnitude, tectonic setting, and resultant ground motion intensities, they produced moderate to severe liquefaction-induced damage in the form of sand boils, building and ground settlements, ground cracking, soil ejecta, and lateral spreading.Furthermore, these effects occurred within bounds defined by an empirical relationship between magnitude and maximum fault distance at which liquefaction was reported (Hu, 2023).
To comprehend the mechanisms behind the observed ground deformations, researchers conduct backanalyses of liquefaction case studies.This involves collecting geotechnical and geophysical data and then utilizing empirical or analytical methods to establish a connection between predictions and observations.This is a challenging task due to the complex and nonlinear nature of liquefaction, which involves significant deformations and fluid flow.Additionally, in-situ data is often sparse and may not accurately reflect the spatial variability in the site conditions, the input ground motion can only be approximated to a limited extent, and predictor models may not be tailored for reproducing the diverse range of surficial manifestations of liquefaction.The most widely known methodology for achieving this goal, the simplified procedure (Seed and Idriss, 1971;Youd and Idriss, 2001;Boulanger and Idriss, 2014), comprises a collection of semi-empirical relationships that proxy the seismic load and soil strength for a given layer with its geomechanical properties derived from in-situ measurements.However, the assumption of free-field conditions (i.e., the absence of a structural load) and 1D vertical stratification of isolated layers falls short in accurately reproducing the dynamic response of the media, causing discrepancies between observations and predicted outcomes as evidenced in recent case studies (Dashti and Bray, 2013;Luque and Bray, 2017;Cubrinovski et al., 2019;Hutabarat and Bray, 2021).
On the other hand, Effective Stress Site Response Analyses (ESSRAs) have been demonstrated to be a more suitable alternative than the simplified procedure for evaluating liquefaction response.In this approach, the seismic load is integrated as equivalent forces applied at the boundaries of a defined soil domain, and the soil media is modeled as a two-phase material (i.e., solid and fluid).By using a proper constitutive model, the wave propagation problem is solved under undrained conditions to reproduce the nonlinear elastoplastic stress-strain response, porewater pressure generation, and post-shaking dissipation generated by an earthquake (Popescu et al., 2006).ESSRAs of case histories have been able to successfully reproduce the pore-pressure buildup, flow patterns, and large shear and volumetric strains that lead to ground deformations consistent with surficial evidence (Bray and Luque, 2017;Luque and Bray, 2017;Bassal and Boulanger, 2021;Pretell et al., 2021;Qiu et al., 2023;Saldaña et al., 2023;Saldaña Sotelo, 2023).Furthermore, parametric studies have been conducted to evaluate the influence of soil properties and ground motion intensity variability on liquefaction response (Popescu et al., 1997;2005;Boulanger and Montgomery, 2016;Montgomery and Boulanger, 2017).
Under the assumption that the constitutive model, the earthquake loading, and the numerical framework to perform the simulations are well defined, the soil domain must be sufficiently characterized to evaluate the liquefaction response of documented case histories.Certainly, research indicates that identifying the strength, extent, and behavior type of geotechnical units is pivotal for understanding the observed deformation patterns (Cubrinovski et al., 2019;Luque and Bray, 2020;Bassal and Boulanger, 2023).Geotechnical vertical soundings, such as the Standard Penetration Test (SPT) and the Cone Penetration Test (CPT), are commonly carried out for the characterization of the soil structure because the measurements provided by these tests correlate to essential parameters that control the liquefaction potential of the soil, such as the relative density (D r ), soil behavior type (I c ), and hydraulic conductivity (k).Shearwave velocity V s measurements are also commonly conducted to constrain the soil's elastic properties.This parameter is defined as where G max is the shear modulus at small strains and ρ is the density of the media.Intuitively, V s is a measurement of the stiffness of the soil and its ability to undergo shear strains.This parameter can be determined through invasive methods, and non-invasive methods, based on the propagation of surface waves generated by passive or active sources.It has been shown that, despite the higher resolution of the invasive methods, the precision of both methods is comparable (Garofalo et al., 2016).The use of shear-wave velocity as a proxy for soil response to liquefaction has a long-standing history: V s depends on the overburden effective stress and the void ratio of the soil (Kayen et al., 2013).Liquefaction potential is highly sensitive to the void ratio, as soils with a high void ratio have the ability to store more pore-water pressure.With that in mind, in the context of the simplified procedure, experimental and in-situ data have been used to derive deterministic and probabilistic empirical relationships between V s and the cyclic resistance and stress ratio of soils (Tokimatsu and Uchida, 1990;Andrus and Stokoe II, 2000;Zhou and Chen, 2007;Kayen et al., 2013).In recent years, data-driven methods based on artificial intelligence have been developed for predicting liquefaction (Zhang et al., 2021;Zhang and Wang, 2021).The performance of these methods has significantly improved when coupling geotechnical data (CPT and SPT) with V s measurements.In the context of site response analyses, V s is an essential parameter for wave propagation problems because it governs seismic wave amplitude and constrains the stress-strain response of soil units.
Of particular interest for this research, the influence of V s spatial variability on ground motion has been highlighted by various studies.In general, the V s structure of a site is three-dimensional due to the presence of basins, topographic irregularities and inherent soil variability.Whether the site response may be accurately represented as 2D or 1D must be assessed from case to case (Thompson et al., 2012;Pilz and Cotton, 2019;Tao and Rathje, 2020;Pilz et al., 2021;Hallal and Cox, 2023).At shallow depths (< 50 m), linear elastic and viscoelastic 2D site response studies have shown that increasing the degree of spatial variability in V s correlates with an increased ground motion variability, especially at frequencies higher than the resonant frequency of the soil (El Haber et al., 2019;Huang et al., 2021;de la Torre et al., 2022a).This can be attributed to the frequency-dependent wave attenuation (Aki, 1980;Assimaki, 2004).Nonlinear site response analyses have also been conducted to assess the effect of V s spatial variability using a total stress formulation (Assimaki et al., 2003).In this case, the hysteretic loops of the soil do not follow the same path during loading and unloading, implying energy dissipation and shear modulus degradation, which in turn leads to increased motion damping when the soil is subjected to higher strains (Hashash and Park, 2001;Hashash and Park, 2002).Keeping this in mind, to the best of our knowledge, there is a limited exploration of the influence of V s spatial variability in liquefaction response through nonlinear effective stress analyses.
This research aims to assess the effect of V s spatial variability in liquefaction triggering, stress-strain response, and liquefactioninduced settlements.To do so, we consider a residential area in Metropolitan Concepción, Chile, where a variable response to liquefaction was reported in the frame of the 2010 M w 8.8 Maule Earthquake (Figures 1, 2).This site has already been thoroughly studied and geotechnically characterized (Bray et al., 2012;Verdugo and González, 2015;Montalva et al., 2022;Saldaña et al., 2023;Saldaña Sotelo, 2023).Using seismic arrays deployed within the site, we retrieve several V s profiles from ambient noise data using the SPatial AutoCorrelation (SPAC) technique (Aki, 1957;Chávez-García et al., 2005;Tsai and Moschetti, 2010).This allows us to identify key soil profiles consistent with the site characteristics.We then generate 2D heterogeneous V s model realizations from the soil profiles, using correlated random fields at different levels of variance.This process is similar to approaches used in non-liquefaction site response analyses (Assimaki et al., 2003;El Haber et al., 2019;de la Torre et al., 2022a;b).Finally, we perform ESSRA, using each model realization as an input, to assess the impact of different levels of variance on the dynamic response of the soil in terms of excess pore water pressure, generated shear and volumetric strains, and vertical settlements.

Background and case study
The megathrust 2010 M w 8.8 Maule Earthquake that struck South-Central Chile ruptured a mature seismic gap that was quiescent since the 1835 M w 8.5 Concepción Earthquake (Campos et al., 2002;Ruegg et al., 2009;Moreno et al., 2010).The rupture propagated bilaterally north and south of the epicenter, covering an along-strike length of around 500 km.The reported damages spanned from the city of Valparaíso (33.0°S) to Temuco (38.7°S) (Bray et al., 2010;Assimaki et al., 2012).Figure 1A illustrates the main metropolitan areas affected by earthquakeinduced damages.The city of Concepción (Figure 1B), located in South-Central Chile, experienced particularly strong ground motions.This behavior was attributed to site and basin effects that amplified the seismic load within the city (Assimaki et al., 2012;Montalva et al., 2016;Inzunza et al., 2019).Particularly, earthquakeinduced liquefaction was reported throughout the city with varying levels of severity, including sediment ejecta, lateral spreading, flow failure, excessive settlements, surface cracks, and structural damage (Bray et al., 2010;Bray et al., 2012;Verdugo and González, 2015;Montalva et al., 2022).
Concepción is located on Pleistocene-Holocene sediments in a fluvial basin created by a horst-graben structure associated with NE-oriented normal faults that cut and displaced Upper Cretaceous to Pliocene sedimentary rocks on top of Upper Paleozoic granitic bodies (Galli, 1967;Vivallos et al., 2010).The basin infill is primarily made up of fluvial sand terraces that were formed due to the Bío-Bío River's meandering after crossing the Coastal Cordillera toward its mouth at the Pacific Ocean.As a result of this process, metric banks of carefully selected sand grains were deposited, creating a relatively uniform unit (Montalva et al., 2016).Borehole, gravimetric, and tomographic studies have shown that the depth to bedrock is variable between outcropping rock up to 160 m, with an average basin thickness of 100 m (Poblete and Dobry, 1968;Inzunza et al., 2019).Geotechnically speaking (Montalva et al., 2016), the basin consists mostly of layers of poorly graded (SP) and silty sands (SM) that were uniformly deposited with non-plastic fines (approximately 22% fine content).This information was gathered from various boreholes conducted across the city.Additionally, layers of lowplasticity or no-plasticity silts (ML) ranging from 1 to 4 m thick can be found at different depths throughout the basin, indicating that they were deposited during a period of low velocity in the Bío-Bío River.
Figure 2A shows a satellite view of our case study area.At the time of the 2010 earthquake, there were four eight-story buildings on the site (towers A, B, C, and D); currently, there are six of them.Surficial evidence of liquefaction was reported on the site in the form of sediment ejecta (brown shapes in Figures 2A, C, D),  (Montalva et al., 2022).The PGAs were computed with Montalva et al. (2017) Ground Motion Prediction Equation.The dashed light-green circle encapsulates our study site, Los Presidentes.The yellow cross depicts the location of strong-motion station CCP, from which the earthquake ground motion was deconvolved to use as input for the site response analysis.ground cracking (orange lines in Figures 2A, E), and building settlements Bray et al., 2010, Bray et al., 2012).A Light Detection and Ranging (LiDAR) image (Figure 2B, Robert Kayen, Personal Communication, and in-situ reports (Bray et al., 2010;Bray et al., 2012) revealed that tower A experienced settlements ranging from 6.9 to 34.5 cm, while the ground surrounding the building settled roughly 22 cm.Tower B, on the other hand, experienced settlements ranging from 7.6 to 10.7 cm, while the surrounding ground settled about 18 cm.These large settlements caused extensive structural damage, leading to the demolition of both towers in 2013; they were subsequently reconstructed between 2016 and 2018.The other two towers, C and D, only experienced minimal settlements.Outside of the studied area, no evidence of liquefaction was reported.In this study, we have focused our region of interest to the area between buildings A and B (Figure 3A).Six geotechnical tests were conducted in this region -three SPT, and three CPT boreholes (cyan stars in Figure 3A).A comprehensive geotechnical characterization of the site, derived from the data obtained through the boreholes, is described in Saldaña Sotelo (2023); Saldaña et al. (2023).

Retrieving dispersion curves from ambient noise data
In the subsequent section, we introduce the seismic interferometry technique in a broad way.First, we explain how through cross-correlations in the time-domain, a representation of the Green Function can be obtained.Following this, we explain the SPAC technique and the assumptions for which it is valid.Finally, we establish a connection between the timedomain cross-correlation method and the frequency-domain SPAC method.
The noise wavefield refers to the combination of seismic waves caused by ambient vibrations of both natural and anthropogenic origin, such as tides, oceanic waves, meteorological phenomena, heavy machinery, and vehicles (Asten and Henstridge, 1984;Bonnefoy-Claudet et al., 2006).This wavefield is primarily composed of surface waves (Chávez-García and Luzón, 2005;Bonnefoy-Claudet et al., 2006).The dispersive nature of the surface waves allows us to illuminate the Earth's internal structure through the seismic interferometry technique (Curtis et al., 2006).Under the assumption that the noise wavefield is diffuse (i.e., waves with uncorrelated random amplitudes and phases propagate in all directions, Weaver, 1982;Lobkis and Weaver, 2001), the Green's Function and the group or phase velocity dispersion curve between two receivers can be calculated by averaging the correlograms over time (Shapiro and Campillo, 2004;Wapenaar et al., 2010).In spite of the strong assumptions made about the wavefield, the widespread success of ambient seismic noise methodologies in recent years has demonstrated their effectiveness in a variety of applications, ranging from soil structure characterization to lower crust/upper mantle tomographic models (e.g., Shapiro et al., 2005;Ritzwoller et al., 2011;Ward et al., 2013;Ma and Clayton, 2014;Inzunza et al., 2019).The bandwidth of the retrieved dispersion curves is primarily constrained by instrumental characteristics, the interstation distance between the receivers, and the deployment time of the seismic array.The number of stacked correlograms, whose duration is frequency-dependent, increases the signal-to-noise ratio (Bensen et al., 2007).In fact, for shallow imaging (i.e., <100 meters), recording for a few hours should suffice (e.g., Picozzi et al., 2009;Pilz et al., 2012).On the other hand, researchers have stacked years of data to obtain dispersion curves for deeper explorations (e.g., Ward et al., 2013;Ma and Clayton, 2014).
The aforementioned approach is based on time-domain cross-correlation.However, a similar method, called SPatial AutoCorrelation (SPAC, Aki (1957), states that if the noise wavefield is stochastic and stationary in space and time, the azimuthal average of the cross-correlation in the frequency domain (i.e., the crosscoherence CC) for a fixed distance r at frequency ω is related to the Rayleigh wave phase velocity c(ω) by where J 0 represents the zero-th Bessel function of the first kind, and A is an amplitude factor that considers attenuation and normalization errors in the CC (e.g., Menke and Jin, 2015).It is noteworthy that the azimuthal average can be replaced by the time average if we consider that the noise wavefield is diffuse (Aki, 1957;Chávez-García et al., 2005).This implies that the crosscorrelation and SPAC methods are equivalent (Chávez-García and Luzón, 2005; Tsai and Moschetti, 2010).With that in mind, the timeaveraged CC is defined as (Ohori et al., 2002;Wapenaar et al., 2010;Xu et al., 2021) where R (S ij ) is the real part of the cross power spectral density (PSD) between station pair i and j, S ii and S jj are the individual PSDs, and the ⟨⋅⟩ operator denotes averaging over time segments.
For each pair of simultaneous recordings, the CC is calculated with Equation 3, and substituted into Equation 2. It is clear from these equations that the relationship between the data and the phase velocity is nonlinear.Still, the dispersion curve c(ω) can be retrieved through the zero-crossings (Ekström et al., 2009;Ekström, 2014) or by fitting a Bessel function to the CC waveform (e.g., Menke and Jin, 2015;Pilz et al., 2017;Olivar-Castaño et al., 2020).

Retrieving ground profiles from dispersion curves
The relationship between the dispersion curve and a onedimensional ground profile is also nonlinear.Furthermore, the dispersion curve includes contributions from both fundamental and higher modes, which cannot be trivially separated.The curve that is composed of this superposition of modes is known as the apparent dispersion curve.This curve is related to each mode through their respective medium response functions (Harkrider, 1964;Tokimatsu et al., 1992;Ohori et al., 2002).Therefore, the nonlinear inverse problem relating the apparent dispersion curve and the ground profile can be solved using adequate methods, such as the simplex downhill method (Nelder and Mead, 1965) or a genetic algorithm (Yamanaka and Ishida, 1996;Parolai et al., 2005).
The genetic algorithm employed in this research is a nonlinear, partially probabilistic approach that aims to explore the entire parameter space in order to find the optimal solution (i.e., the space within the predefined lower and upper bounds for density and P-S wave velocity).Since each inversion is defined by a random seed, all resulting solutions are unique, even if they start with the same parameterization.Consequently, there are two significant issues regarding the inverse problem solution and parameterization: the non-uniqueness of the solution and how well it is constrained by the available data.The first issue implies that many different ground profiles may fit the data reasonably well, making it difficult to determine a single 'true' solution.The second issue highlights the limitations of dispersion curve data and how this data constrains each layer's thickness and the maximum depth exploration.A layering-byratio scheme was designed to address these issues Cox and Teague (2016); Vantassel and Cox (2021).This approach allows for the exploration of various parameterizations, with different numbers of layers, in order to determine the optimal number of layers that suitably fit the data, preventing over or under-parameterization.Moreover, it ensures that no layers begin below the wavelengthdefined spatial resolution.The best-constrained solution can be selected by minimizing the misfit between the data and the prediction, and a measurement of the inter-parameterization uncertainty can be obtained from the N-lowest misfit solutions (e.g., N = 10, 100).
4 Characterization of the site velocity structure

Data and methods
In the frame of this study, we conducted five geophysical surveys from November 2021 to December 2022.Three-component high-frequency seismic Tromino ® instruments, with a sampling rate of 512 Hz, were deployed at the study site to record the ambient seismic wavefield at a total of 50 station positions (magenta triangles in Figure 3A).The instruments were configured to record synchronously for a duration of approximately 40 min at each position.For each pair of simultaneous recordings (black lines in Figure 3A), the CC was calculated with Equation 3 using the Welch Method (Welch, 1967).This method involves computing periodograms from detrended segments of 60 s in length, with a 50% overlap between segments.Additionally, a 5% cosine taper was applied to both ends of each segment, followed by a one-bit normalization filter.This procedure reduces spectral leakage and increases the signal-to-noise ratio (Bensen et al., 2007).Finally, the CC was smoothed using a moving-average filter and substituted into Equation 2. The lower and upper phase velocity limits for the solution were set to 100 [m/s] and 400 [m/s], respectively, based on results from previous surveys (Montalva et al., 2022;Saldaña et al., 2023;Saldaña Sotelo, 2023).We defined the minimum resolvable wavelength, λ min , as half the interstation distance (r/2).Due to the different interstation distances, we manually defined the frequency band for the 198 CC waveforms produced.Furthermore, we conducted a visual inspection of the resulting dispersion curves, obtained using the methodology of Menke and Jin (2015), discarding any curves that fit the CC waveforms poorly and outliers corresponding to geologically-implausible interpretations (i.e., they do not align with the fluvial sand-filled Concepcion basin).A subset of 123 curves, presented in Figure 3B, aligned with the aforementioned criteria.The minimum and maximum velocities are 116 and 320 [m/s], respectively, while the frequency range spans from 2.85 to 25 Hz.The majority of dispersion curves fall within the 6-16 Hz frequency band, and are representative of very soft, shallow conditions.Subsequently, each dispersion curve is transformed from linear-frequency to the log-wavelength domain.Given that the variation of phase velocity as a function of wavelength is less than as a function of frequency, this transformation reduces the gap between points in the dispersion data without needlessly increasing the number of samples (Vantassel and Cox, 2021).
We inverted ground profiles for each dispersion curve using a generic algorithm (Yamanaka and Ishida, 1996;Parolai et al., 2005) with parameterizations consisting of 3, 4, and 5 layers for each curve, with the deepest layer representing the half-space base.The V s lower and upper bounds of each layer (regardless of their depth) were set to depend on the lower and upper values of the dispersion curve, ranging from approximately 100 to 350 m/s for the upper sedimentary layers and up to 450 m/s for the halfspace layer.As the generic algorithm is part-probabilistic, each dispersion curve was inverted using 12 seeds, resulting in a total of 36 inversions for each curve considering the different numbers of layers used.Then, through visual inspection of the ground profiles and their goodness-of-fit to the data, we determined the number of layers most adequate for each curve and identified the representative ground profile as the one exhibiting the lowest misfit with respect to the dispersion curve.Furthermore, from the 12 solutions, we estimated the inter-model uncertainty by calculating σ ln,v s , representing the logarithm of the standard deviation of the velocity.

Inverted ground profiles
Figure 4 displays three representative ground profile solutions, their uncertainties, and the dispersion curve fit obtained through the aforementioned procedure.Subsequently, we narrow our focus on the first 30 m of the V s profiles, as the deeper structure is poorly resolved by the frequency range of our phase velocity data.It can be seen that the thickness of the upper layers decreases as the interstation distance decreases because the wavelengths constrained by the dispersion curves become shorter.Additionally, the uncertainties (blue dashed lines in Figure 4A) tend to increase near the depths corresponding to layer interfaces.This happens because the trade-off between layer thickness and V s cannot be resolved using only dispersion data.Nevertheless, the ground profile uncertainties are within the expected range for surface wave surveys (Vantassel and Cox (2021) and references therein).

Merging scheme
We used a simple merging scheme to quantify the spatial variability of V s from the ground profiles, reminiscent of the approximative tomographic inversion proposed in Kissling (1988).First, we divided the site into a rectangular grid of 35 m in the x direction and 42 m in the y direction (Figure 5).Assuming that the surface waves travel along straight lines linking each pair of receivers, the velocity profile corresponding to each grid cell is obtained by weighting all the ground profiles corresponding to pairs of stations whose interstation path passes through it.The weighting scheme takes into account the length of each ray within the grid cell and the σ ln,v s of each profile.In other words, the ground profile associated with the ray with the longest length within a specific grid cell carries greater weight than rays crossing shorter distances in the same grid cell.Additionally, ground profiles with lower values of σ ln,v s for a given depth also have greater weight.Then, for a given cell that is crossed by n rays of cell length l and logarithmic standard deviation σ z at a given depth, the shear-wave velocity at depth z, V z , is given by where the subindex s of the shear-wave velocity term V s is omitted here for clarity.Equation 4 ensures that rays covering more distance within a cell carry more weight, that short interstation distance models have more weight at shallow depths, and that large interstation distance models have more weight at greater depths.

Resulting velocity structure and representative profiles
Figures 5A-C display the average velocity structure of the site at depths of 0-10 [m], 10-20 [m], and 20-30 [m].Generally, lower velocities are observed in the shallow structure (Figure 5A), with increasing velocities at greater depths (Figures 5B, C).Notably, the grid cells near the towers show relatively high velocities at all depths with respect to the middle and southernmost cells.This can be attributed to the fact that the southernmost part of the study area is located near public streets, which were likely not consolidated by heavy machinery before the building reconstruction.On the other hand, the cells in the middle are in an area that was likely consolidated and refilled before the reconstruction of the towers, but to a lesser extent than the cells adjacent to the structures.Moreover, the weight of the buildings themselves may also affect Frontiers in Earth Science 08 frontiersin.orgThe velocity structure of the site calculated by the weighting scheme.the velocity structure of the adjacent cells.Figure 5E shows the weighted ground profiles for all the grids (grey profiles).The black profile represents cells directly adjacent to the buildings.The blue profile represents cells in the ground area between the buildings.The red profile represents cells south of the buildings.As severe liquefaction was observed near buildings A and B in 2010, the red profile (soils not consolidated by heavy machinery prior to reconstruction) is considered the most representative of the freefield natural conditions that existed for our simulations.

Numerical model setup
We perform ESSRA using the Open System for Earthquake Engineering Simulation (OpenSEES) finite-element framework (McKenna, 2011) and the Scientific Toolbox for OpenSEES (Petracca et al., 2017).The wave propagation over fluid-saturated porous media is solved in a 2D domain under plane-strain conditions.We consider the cross-section T-T ′ depicted in Figure 3A, a domain of 140 m in length in the x-direction and 50 m in depth in the z-direction, as shown in Figure 6A.In our modeling, we did not account for the influence of building's load.To enforce free-field conditions at the lateral boundaries of the model, we model free-field columns of 10 m length in the x-direction and 10,000 m size in the out-of-plane-direction, and tie the displacement degree-of-freedom to enforce periodic boundary conditions, in a very similar manner to what is performed.The nodes at the base of the model were fixed in the z-direction, under the assumption that the domain is underlain by an elastic homogeneous half-space, and they were forced to move horizontally in the same direction by tying them to the lower-left node.The ground motion was input at this lower-left node as a force time history with a Lysmer-Kuhlemeyer dashpot (Lysmer and Kuhlemeyer, 1969;Joyner and Chen, 1975), and propagated to the constrained base nodes and across the domain.The LK dashpot method requires to input a coefficient c = ρ b V s b A elem , with ρ b = 2.0[g/cm 3 ] and V s b = 339[m/s] being the density and shear-wave velocity of the underlying medium (i.e., the elastic homogeneous half-space) and A elem = 10 * 10000m 2 .As neither borehole nor outcrop bedrock seismic recordings were available near our study site, the ground motion at the base was estimated by applying simple linear deamplification, to a depth of 50 m, of the ground motion recorded at station CCP (see Figure 1 for station location; Figures 5F, G for the station ground motion).As mentioned earlier, our profiles are well-constrained only until a depth of 30 m.To address this issue, we created a joint profile by taking the median of all of our non-weighted ground profiles until 30 m, and using the profile derived for the Concepción basin from Inzunza et al. (2019) below this depth.The merged profile is shown in Table 1.
Effective stress conditions are enforced using stabilized-singlepoint quadrilateral u − p elements (SSPQUADUP, McGann et al., 2012), which are based on the u − p formulation.This formulation assumes that saturated soils are a continuum composed of a solid and a fluid phase, with the displacement of the solid phase u and the pore-fluid pressure p of the fluid phase being the main variables.These underlying assumptions are valid for most earthquake and soil dynamics problems (Biot, 1956;Zienkiewicz and Shiomi, 1984).SSPQUADUP elements include an additional pressure degree of freedom.As the groundwater table at the site varies from roughly 0.5-1.5 m (Saldaña et al., 2023;Saldaña Sotelo, 2023), we fixed the pore-water pressure degree-of-freedom at the surface level.Benchmark finite element model realization of the T-T ′ cross-section defined in Figure 3. (A), top panel: V s random field obtained using σ ln,v s = 0.075.Free-field columns at the lateral boundaries are depicted in thick black rectangles.The red arrows start and end at the main and constraint nodes, respectively.Input ground motion is applied at the green node in the lower-left corner.Base nodes (blue) are constrained to the green node.The recorder nodes at a depth of 3 m are depicted as magenta triangles and are further analyzed in Figure 9. (B), lower panel: D r structure of the cross-section (Saldaña et al., 2023;Saldaña Sotelo, 2023).
To model the elastoplastic behavior of the sands in the shallower 30 m, we used the PressureDependentMultiYield03 model (Khosravifar et al., 2018).This multi-yield surface model was originally designed to capture the cyclic mobility and postliquefaction accumulation of shear strain on sands and has been updated to account for the influence of the number of loading cycles on liquefaction triggering.A thorough description of the model formulation can be found in Parra-Colmenares (1996); Yang et al. (2003); Khosravifar et al. (2018).We highlight from the abovementioned literature that (1) the model assumes that elastic and plastic deformations occur simultaneously in the soil, and the elastic behavior is linear and isotropic, while the plastic behavior is nonlinear and anisotropic, and (2) the soil nonlinear shear stressstrain response is defined in the octahedral space in the following manner: the pressure-dependent small-strain shear modulus G p max is defined by where G max is the input shear modulus computed with Equation 1 at a reference effective confining stress p ′ r , d is the pressure-dependent coefficient set to 0.5, and p′ is the effective confining stress that varies during the earthquake loading.Thus, shear modulus reduction curves are computed using a hyperbolic relationship given by Where τ oct is the octahedral shear stress, γ oct is the deviatoric strain in the octahedral space, and γ r is a parameter that constrains the shape of the backbone curve.Equation 6 describes the shape of the hysteretic loops for a given effective stress, small-strain shear modulus at a given pressure (given by Equation.5), and  6B).Table 2 summarizes the model parameters used for our analyses.To simulate V s heterogeneity and wave attenuation, we followed the method of de la Torre et al. (2022a); de la Torre et al. (2022b).We generated perturbations of the 1D V s profile shown in Figure 5E, using spatially anisotropic correlated Gaussian random fields with varying σ ln,v s values, to simulate different levels of heterogeneity.The parameters that we need to define to achieve this are the following: • The logarithm of the standard deviation of the velocity σ ln,v s : Here we set five different values: 0.075, 0.125, 0.175, 0.225, and 0.275, to simulate conditions with very low to high variance.
Values mentioned in de la Torre et al. (2022a), and references therein range from 0.1 to 0.37.The σ ln,v s values obtained for our profiles in Figure 4 vary roughly from 0.04 to 0.23, consistent with the aforementioned literature and the values chosen.We run 20 simulations per level of variance, resulting in a total of 100 simulations.• The horizontal r H and vertical r V correlation lengths were set to 15 and 2 m.This is within the bounds of values previously used for generating heterogeneous V s random fields.(Popescu, 1995;Assimaki et al., 2003;de la Torre et al., 2022a).An exponential function was used to simulate the spatial correlation.
Additionally, a zero-variance zone of 25 m length was established at the lateral boundaries, where the 1D deterministic soil profile is valid.Figure 6 shows the benchmark V s realization for σ ln,v s = 0.075.Lastly, a single homogeneous elastic layer of V s = 339 m/s was used for the subdomain within 30 and 50 m depth.We implemented the procedure of Ramirez et al. (2018); Tiznado et al. (2021) to determine the vertical size of the elements efficiently.For small-strain problems, the maximum vertical element size, h max , depends on the V s value of the element and the maximum frequency of the input motion: h max = V s 4 f max .As the shear modulus degrades in nonlinear problems, this implies that V s will decrease as the simulations progress.Therefore, to ensure that our elements do not exceed their theoretical maximum size, we further divided h max by a factor of 4 to account for this phenomenon and guarantee that a sufficient number of elements cover one wavelength at all times.

Simulation results
We now describe and analyze the simulation results from five simulation outputs: • The shear strain γ.
• The volumetric strain ϵ v , defined as the mean of the vertical and horizontal strains.• The excess pore-water pressure ratio r u , defined as the difference between the current and initial pore-water pressure divided by the overburden effective stress.• The vertical settlements u z .
• The horizontal acceleration a x Threshold values previously defined in the literature for liquefaction triggering are between 80% and 100% for r u , and between 3% and 5% for the shear strains (Ishihara, 1993;Boulanger et al., 1998;Bray and Sancio, 2006).

Benchmark simulation
Figure 7 shows the maximum shear strain γ, maximum volumetric strain ϵ v , excess pore water pressure ratio r u , maximum horizontal acceleration a x , and vertical settlement u z for the benchmark model realization shown in Figure 6.We focus on three representative nodes, which are depicted as white circles 1, 2, and 3 in Figure 7.For the zone surrounding representative node 1, liquefaction of a loose sand layer (see Figure 6B) occurs at depths between 7 and 8 m due to the accumulation of large strains (Figures 7A, B) and r u values exceeding 100% being located in this part (Figure 7C).The highly nonlinear hysteretic loops of the stress-strain curve and the r u time series are also consistent with liquefaction phenomena (Figure 8A).Moreover, a steep gradient of horizontal accelerations is observed below the node (Figure 7D).
For the zone surrounding representative node 2, we observe slightly smaller maximum r u values (∼90 − 100%, Figure 7C) and maximum shear strains in the range of 2-4% (Figure 7A).These r u and shear strain values are close to the previously stated threshold values (Ishihara, 1993;Boulanger et al., 1998;Bray and Sancio, 2006).The stress-strain response recorded is nonlinear (Figure 8B), but to a considerably lesser extent than node 1.Furthermore, the maximum volumetric strains are negative (indicating contraction) near node 2, but the surface values are positive (indicating dilation) (Figure 7A).The area around node 2 is the one that exhibits the largest vertical settlements, reaching roughly 50 mm (Figure 7E).
As for the zone surrounding representative node 3, no liquefaction occurs.Indeed, r u values do not exceed 50% (Figure 7C), the stress-strain behavior is close to linear (Figure 8C), no significant settlements are observed (Figure 7E), and the Mass density ρ [g/cm 3 ] 1.9 1.9 1.9 1.9  maximum accelerations are significantly higher than in the zone around node 2 (Figure 7D), implying that this zone did not experience significant soil softening.As shown in Figure 6B, node 3 is embedded in a zone of dense sands, which are less susceptible to liquefaction compared to the layers in which nodes 1 and 2 are situated.

Assessing the effect of V s variability on key dynamic properties
In Figure 9, we plot, for different levels of σ ln,v s , the median maximum γ, ϵ v , r u , and u z of all simulations as a function of the distance along the cross-section at a depth of 2 m.The cross-section's nine nodes-five positioned to the west and four to the east-display distinct behaviors.Consequently, we conduct separate analyses for the two cases.For the western nodes, the high r u (80-100%) values for all σ ln,v s levels are indicative of liquefaction behavior.Particularly, we can see that an increase in σ ln,v s leads to a decrease in the median maximum γ and ϵ v .In a similar manner, a slight reduction in the median maximum u z is observed at the nodes x = 46 [m] and x = 54 [m] as σ ln,v s drops.This tendency is not observed for the rest of the western nodes, and the median maximum u z values are very similar.Additionally, no significant changes in the median maximum r u values are observed.
For the eastern nodes, the low r u (45-50%) values and the developed small strains and settlements indicate that no liquefaction occurred.Furthermore, we appreciate very slight variabilities in the soil response for different levels of σ ln,v s .Only a slight increase in the median maximum r u values with increasing σ ln,v s is observed, which is a contrary behavior to what is seen in western nodes.Anyhow, variations on σ ln,v s did not signify a change in the eastern nodes from non-liquefaction to liquefaction behavior.Although the values are slightly different, a very similar behavior for the western and eastern nodes is observed at depths of 4 and 5 m (Supplementary Figures S1, S2).
Figure 10 shows the time series of the aforementioned variables for all the simulations (grey lines) and the median time series at different levels of σ ln,v s for the western node x = 54 [m] (see Figure 9 for node location).This is the node that exhibited the largest median Maximum shear strain (A), volumetric strain (B), excess pore-water pressure ratio (C), horizontal acceleration (D), and settlements (E) for the benchmark model realization shown in Figure 6.The magenta numbered circles represent the locations of the three nodes depicted in Figure 8. maximum shear strain and settlements.Overall, we see that all the median time series closely follow each other in the first 20 s, when r u values are not large.After the 20 s, and even more markedly at around 45 s (i.e., when r u stabilizes and stops increasing), differences between the σ ln,v s levels become apparent.In the same line, we see that all simulations (regardless of σ ln,v s ) follow the pattern of permanent strain accumulation when high pore-water pressures are developing.This implies that liquefaction behavior occurs regardless of the simulation σ ln,v s level.Additionally, the same tendency of greater σ ln,v s leading to smaller strains and settlements observed in Figure 9 is noted for the median time series.It can be appreciated as well that starting from t > 40s, the higher median r u values happen for σ ln,v s = 0.075, but the increase is slight.
Black curves in Figure 10 show all model simulations' 20-th and 80-th percentiles.For the strains and r u time series, σ ln,v s = 0.075 and σ ln,v s = 0.275 curves follow closely the 80-th and 20-th percentile curves, respectively.On the other hand, we see from the vertical displacement time series that the percentile curves are farther off from the median curves than for the strain and r u time series.Actually, the settlement u z (i.e., the last value of the time series) for the 20-th and 80-th percentiles considerably varies from 41 to 53 mm.We can also see that the volumetric and vertical displacements are still increasing at the last time step.This could indicate that flow is still occurring after the strong motion has ended, and hence larger settlements could be obtained by increasing the simulation time.A similar pattern is observed for the time series at depths of 4 and 5 m (Supplementary Figures S3, S4).
For comparison, we also show the time series for the eastern node x = 86 [m] (Supplementary Figure S5).From here, we mention two noteworthy points: The first is that shear strain, pore-water pressure ratio, and vertical displacement median time series for each σ ln,v s are much more stacked together than in the previous case.Additionally, we appreciate more variability in the volumetric strains, although the maximum magnitude of the strains accumulated is smaller than 0.05.With that in mind, even though the 20-th and 80-th of volumetric strains and vertical displacements are farther off from the median curves, their values are very small.The second is that no liquefaction behavior is observed even for the most extreme outliers simulations (i.e., the ones that exhibit the largest r u ).This is in agreement with what we had shown for the eastern nodes in Figure 9.

Discussion
Our findings indicate distinct responses between western and eastern nodes due to the varying levels of V s spatial variability.While strains and settlements tend to decrease as variability increases for the western nodes, the response of the eastern nodes remains relatively unchanged.In Figure 10, we see that the western node was subjected to large plastic deformations across all σ ln,v s levels, while Supplementary Figure S5 demonstrates minimal deformations for all simulations at the eastern node.From this, and recalling from Figure 6 that the shallow western part of the domain corresponds to soil units with low relative density, and the shallow eastern part of the domain corresponds to units with high relative density, it can be inferred that V s spatial variability may not alter the soil's intrinsic response to liquefaction.This response is likely primarily controlled by the soil's mechanical properties and the earthquake loading.In other words, our results support the notion that V s spatial variability does not play a primary role in liquefaction triggering, but rather influences the liquefaction response when liquefaction is already ongoing.
We propose a physical explanation for the decrease in liquefaction-induced settlements and strains as σ ln,v s increases.As seen in Figure 10A, the median curves diverge as nonlinear behavior increases (i.e., plastic deformation develops during the simulation).We can infer that the more nonlinear the behavior of the soil is, the increased importance of σ ln,v s becomes apparent.Certainly, it is known that for a given shear-wave velocity value, the soil damping increases with larger shear strains (Hashash and Park, 2001;2002).Therefore, if variability increases, the overall soil stiffness decreases, intensifying the damping and reproducing the observed deformation patterns.It should be noted that this pattern depends on the median V s values and the confining pressure of the soil units.If V s is high and the soil unit is deep, the damping will be lower due to the increased shear strength of the soil, reducing the importance of σ ln,v s .See, for instance, Supplementary Figure S6, which displays the shear-strain time series of node number 1 from Figure 7, where significant liquefaction was observed between 7 and 8 m depth.Even though larger shear strains close to 10% are appreciated, the median time series do not exhibit a clear relationship between σ ln,v s and the developed shear strains.In summary, we propose that the influence of σ ln,v s on the liquefaction response becomes relevant in the following circumstances: 1.The ground motion intensity is strong enough to trigger nonlinear behavior and significant r u values.2. The mechanical properties of the soil allow liquefaction to occur.3. The median V s of the soil at a given depth is sufficient to allow significant damping when subjected to large (e.g., >1) shear strains for a given ground motion intensity.
In terms of liquefaction-induced settlements, one of the most detrimental liquefaction hazards, our results show that the effect of V s spatial variability on the median computed settlements was consistent but slight.However, the settlements' 20-th and 80-th percentiles are 41 and 53 mm, respectively.This implies that V s variability may significantly influence the computed settlements.Therefore, it is crucial to constrain the velocity structure of a site accurately in order to gain a better insight into the liquefaction response of a case study.With this in mind, we identified three distinct V s profiles consistent with the site's geophysical and geotechnical characteristics within our small study site.For bigger soil domains (e.g., Pretell et al., 2021;Qiu et al., 2023, we can expect to see even more diversity on V s signatures. We note that the calculated maximum shear strains near the surface are not larger than 3%, which is not considered to be exceedingly large.Saldaña Sotelo (2023); Saldaña et al. (2023) conducted ESSRAs at the same study site and obtained maximum surface shear strains near the surface up to 10%.Additionally, simulations of Luque and Bray (2017) on two buildings -that exhibited liquefaction-induced settlements during the Canterbury Earthquake Sequence of 2010-2011 -reached maximum shear strains beneath the buildings of around 8%.We suggest that the effect of V s heterogeneities on near-surface deformations may be more significant in these cases compared to our study.
Even though it was not our primary research goal, we compared the computed settlements at the representative western and eastern nodes with the ground settlements reported in Bray et al. (2010Bray et al. ( , 2012) ) (Table 3).Our average settlements for both nodes are  significantly lower than the reported values, especially for the eastern node, where no liquefaction behavior was captured in our simulations.We suggest motives for that discrepancy.Recent studies have shown that including structures in the simulations may greatly increase the vertical settlements in the proximity (Bray and Luque, 2017;Saldaña et al., 2023;Saldaña Sotelo, 2023).At our study site, Saldaña Sotelo (2023); Saldaña et al. (2023) calculated the liquefaction-induced ground settlements from the available CPT data using Zhang et al. (2002) simplified methodology, obtaining values from 10 to 37 mm.These values also do not match the observations, but are more in alignment with what we obtained.Vertical settlements computed with empirical relationships exhibit mismatches with observations because they (1) assume free-field conditions, (2) only consider volumetric-induced settlements, and (3) cannot reproduce the complex system response as numerical modeling does.While our modeling considers volumetric and shear-induced settlements and reproduces the system response, it was conceived under a free-field assumption.We therefore suggest that not accounting for the buildings' load is the principal reason our simulations underestimate the settlements.Other reasons, such as how the seismic load is incorporated into the domain and the constitutive model employed, might also be relevant, but to a lesser extent.

Conclusion
We performed effective stress analyses (ESSRAs) at the Los Presidentes site, which experienced liquefaction damage during the 2010 M w 8.8 Maule Earthquake, in order to assess the effects Frontiers in Earth Science 16 frontiersin.orgof V s heterogeneities on liquefaction response at the site.Using seismic interferometry methods, we identified three distinct velocity profiles consistent with the geological setting of Metropolitan Concepcion, and representative of the site conditions.Using the profile that most accurately represented free-field conditions, we generated 2D V s heterogeneous model realizations at various levels of σ ln,v s using Gaussian correlated random fields.Our simulation results demonstrate that increasing σ ln,v s in near-surface soil elements, which exhibit nonlinear and liquefaction-like behavior due to earthquake-induced pore-water pressure buildup, leads to a decrease in median maximum shear, volumetric strains, and vertical settlements.We argue that when nonlinear behavior begins, an increase in damping at the element level is responsible for this reduction in the computed measurements.Furthermore, the influence of σ ln,v s on the liquefaction response is more pronounced at shallow depths, where both V s values and the confining pressure values are lower because the damping of the soil is stiffnessdependent.On the other hand, increasing σ ln,v s did not significantly alter the deformation patterns of soil elements that exhibit nonliquefaction behavior.This implies that V s influence on liquefaction triggering is minimal and that the soil's mechanical properties and the intensity of the seismic motion are the primary factors in determining whether liquefaction is triggered or not.
In the frame of liquefaction hazard assessment, the maximum computed settlements at the 20-th and 80-th percentiles were 41 and 53 mm, respectively.This highlights a significant difference considering that our simulations were conducted without accounting for the building's load, and that the computed maximum shear strains were relatively small (<3%).We infer that for higher amplitude ground motions and softer soil conditions, V s variability may play an important role in the final settlements and deformations.
While our findings are within the frame of a case study of liquefaction due to a megathrust earthquake, this analysis can be extended to other tectonic and geotechnical environments.In this regard, we encourage researchers to thoroughly characterize the velocity structure of a site when conducting ESSRAs for past or expected liquefaction case studies as V s heterogeneities, which are expected to be found in the scale of liquefaction problems, can significantly alter the dynamic response of liquefiable soils, and become a key factor for assessment.
10.3389/feart.2024.1354058their affiliated organizations, or those of the publisher, the editors and the reviewers.Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

FIGURE 1 (
FIGURE 1 (A), upper panel: Shaded relief map of Central Chile.Red dots indicate major cities subjected to damages in the frame of the Maule 2010 Earthquake (epicenter depicted as a yellow star).The yellow rectangle represents Metropolitan Concepción, zoomed in the lower panel (B).Lower-right inset represents Central Chile (blue rectangle) within South America.(B) Satellite View of Metropolitan Concepción.Peak ground acceleration (PGA) color-coded circles show sites where surficial evidence of liquefaction was reported(Montalva et al., 2022).The PGAs were computed withMontalva et al. (2017) Ground Motion Prediction Equation.The dashed light-green circle encapsulates our study site, Los Presidentes.The yellow cross depicts the location of strong-motion station CCP, from which the earthquake ground motion was deconvolved to use as input for the site response analysis.

FIGURE 2
FIGURE 2 Summary of the damages at the Los Presidentes site.(A), upper left: Site boundary, damages in the form of ground cracks and boil ejecta.Only buildings A, B, C, and D are labeled because the northernmost ones were not built at the time of the earthquake.(B), lower left: A LiDAR image with the liquefaction-induced building settlements measured at the site (Robert Kayen, personal communication).The color palette indicates the distance between the laser sensor and the scanned surface.In this representation, red indicates larger distances, while blue corresponds to smaller distances.(C,D), top and middle right: Evidence of soil ejecta at the northeastern corner and inside the residential complex, respectively.(E), bottom right: Evidence of ground cracking between buildings A and C. Figure was extracted from Saldaña Sotelo (2023), and references therein.

FIGURE 3 (
FIGURE 3 (A), left: Schematic view of the study site.Magenta triangles depict the positions of seismic sensors, black lines connecting them show the station pairs from which dispersion curves shown in (B) were retrieved.The dashed red line represents the T-T ′ cross-section employed in the numerical model.(B), right: Distribution of all Rayleigh Wave dispersion curves across the study site.Cyan stars illustrate the location of the geotechnical borehole tests.
FIGURE 4 (A) Top panels: The ground profiles inverted from the dispersion curves displayed in the lower panels.Grey profiles show the lowest misfit solution for each seed.The red profile is the lowest-misfit solution for all seeds.The dashed blue line represents the logarithm of the standard deviation of V s .(B) Lower panels: Dispersion curve fit (dotted black line) by each ground profile.Colors represent the same as at the top panels.

FIGURE 5
FIGURE 5 (A-C), leftmost top panels: Average V s structure at depths of (0-10 [m]),(10-20  [m]) and(20-30 [m]), respectively.In (A-C), the X and Y coordinates are within the coordinate system of Figure3A.(D), top right: Distribution of ground profiles averaged on each grid to obtain the weighted ground profile.(E), lower left: The weighted ground profiles of all grids.Three representative profiles of different ground conditions are depicted with red, blue, and black lines.The red velocity profile, which represents free-field conditions, is used for the ESSRAs.(F, G), rightmost lower panels: Black and orange ground motion and their Fourier spectrums depict the CCP EW station recording and its deconvolution to 50 m depth, respectively.

FIGURE 8
FIGURE 8 Recorder acceleration, excess pore-water pressure, and stress-strain time series of the numbered nodes of Figure 7. (A), top: Node 1. (B), middle: Node 2. (C), bottom: Node 3. The time series' black, red, and blue color codes represent the first 15 s of motion, from 15 to 75 s, and from 75 to the end of motion.

FIGURE 9
FIGURE 9 Median of the Max.Shear strain (A), top left; max.Volumetric strain (B), top right; max.PWP ratio (C), bottom left; and vertical settlements (D), bottom right, at different levels of σ ln,v s across the T-T ′ cross-section.Nodes are located at a depth of 2 m.Thick black vertical line divides the western from eastern nodes.Thick red vertical lines depict the location of the western and eastern nodes analyzed in Figure 10 and Supplementary Figure S5.Horizontal black lines near the x-axis represent the locations of towers A and B relative to the cross-section.
FIGURE 10All simulations of shear strain (A), volumetric strain (B), vertical displacement (C), and excess pore-water pressure ratio (D) time series computed for the node at the western node at position x = 54 and depth of 2 m, depicted in Figure9.Color-coded time series represent the median time series at different levels of σ ln,v s .Black curves represent all simulations' 20-th and 80-th percentile time series.

TABLE 1
Saldaña Sotelo (2023)obtained by merging the median profile from this work and the uppermost two layers of the profile derived inInzunza et al. (2019).Most of the model input parameters apart from G max and the Bulk Modulus at a reference pressure B r , which can be computed from G max , have already been calibrated by the developers for different relative densities.Therefore, we mapped the D r values from the geostatistical model described inSaldaña Sotelo (2023);Saldaña et al., 2023 into our domain (Figure

TABLE 2 Pressure
Dependent Multi Yield 03 model parameters employed in this research.