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

Abstract

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 in the amplitude of the strains and settlements as the degree of V_s variability increases.

1 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, pore-water 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 build-up, 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). Shear-wave velocity V_s measurements are also commonly conducted to constrain the soil’s elastic properties. This parameter is defined aswhere 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 liquefaction-induced 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.

FIGURE 1

FIGURE 2

2 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 earthquake-induced 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, earthquake-induced 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 low-plasticity 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), 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).

FIGURE 3

3 Theory

3.1 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 time-domain 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., 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 cross-coherence CC) for a fixed distance r at frequency ω is related to the Rayleigh wave phase velocity c(ω) bywhere J₀ 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 cross-correlation and SPAC methods are equivalent (Chávez-García and Luzón, 2005; Tsai and Moschetti, 2010). With that in mind, the time-averaged CC is defined as (Ohori et al., 2002; Wapenaar et al., 2010; Xu et al., 2021)where 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).

3.2 Retrieving ground profiles from dispersion curves

The relationship between the dispersion curve and a one-dimensional 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-by-ratio 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 wavelength-defined 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

4.1 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 half-space 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 , representing the logarithm of the standard deviation of the velocity.

4.2 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).

FIGURE 4

4.3 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 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 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 bywhere 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.

FIGURE 5

4.4 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 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 free-field natural conditions that existed for our simulations.

5 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 , with ρ_b = 2.0[g/cm³] and being the density and shear-wave velocity of the underlying medium (i.e., the elastic homogeneous half-space) and A_elem = 10*10000 m². 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.

FIGURE 6

Effective stress conditions are enforced using stabilized-single-point 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.

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 post-liquefaction 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 stress-strain response is defined in the octahedral space in the following manner: the pressure-dependent small-strain shear modulus is defined bywhere G_max is the input shear modulus computed with Equation 1 at a reference effective confining stress , 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 seismic excitation. 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 in Saldaña Sotelo (2023); Saldaña et al., 2023 into our domain (Figure 6B). Table 2 summarizes the model parameters used for our analyses.

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 . 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: . 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.

6 Simulation results

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).

6.1 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).

FIGURE 7

FIGURE 8

For the zone surrounding representative node 2, we observe slightly smaller maximum r_u values (, 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 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.

6.2 Assessing the effect of V_s variability on key dynamic properties

In Figure 9, we plot, for different levels of , 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 levels are indicative of liquefaction behavior. Particularly, we can see that an increase in 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 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.

FIGURE 9

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 . Only a slight increase in the median maximum r_u values with increasing is observed, which is a contrary behavior to what is seen in western nodes. Anyhow, variations on 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 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 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 levels become apparent. In the same line, we see that all simulations (regardless of ) follow the pattern of permanent strain accumulation when high pore-water pressures are developing. This implies that liquefaction behavior occurs regardless of the simulation level. Additionally, the same tendency of greater 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 , but the increase is slight.

FIGURE 10

Black curves in Figure 10 show all model simulations’ 20-th and 80-th percentiles. For the strains and r_u time series, and 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 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.

7 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 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.

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. (2010, 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.

8 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 of 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 using Gaussian correlated random fields. Our simulation results demonstrate that increasing 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 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 stiffness-dependent. On the other hand, increasing did not significantly alter the deformation patterns of soil elements that exhibit non-liquefaction 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 . 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.

References

• 1

  Aki K. (1957). Space and time spectra of stationary stochastic waves, with special reference to microtremors. Bull. Earthq. Res. Inst.35, 415–456.

  • Google Scholar

• 2

  Aki K. (1980). Attenuation of shear-waves in the lithosphere for frequencies from 0.05 to 25 Hz. Phys. Earth Planet. Interiors21, 50–60. 10.1016/0031-9201(80)90019-9

  • CrossRef
  • Google Scholar

• 3

  Andrus R. D. Stokoe K. H. II (2000). Liquefaction resistance of soils from shear-wave velocity. J. Geotechnical Geoenvironmental Eng.126, 1015–1025. 10.1061/(asce)1090-0241(2000)126:11(1015)

  • CrossRef
  • Google Scholar

• 4

  Assimaki D. (2004). Topography effects in the 1999 Athens earthquake: Engineering issues in seismology. Ph.D. thesis, Massachusetts Institute of Technology.

  • Google Scholar

• 5

  Assimaki D. Ledezma C. Montalva G. A. Tassara A. Mylonakis G. Boroschek R. (2012). Site effects and damage patterns. Earthq. Spectra28, 55–74. 10.1193/1.4000029

  • CrossRef
  • Google Scholar

• 6

  Assimaki D. Pecker A. Popescu R. Prevost J. (2003). Effects of spatial variability of soil properties on surface ground motion. J. Earthq. Eng.7, 1–44. 10.1080/13632460309350472

  • CrossRef
  • Google Scholar

• 7

  Asten M. W. Henstridge J. (1984). Array estimators and the use of microseisms for reconnaissance of sedimentary basins. Geophysics49, 1828–1837. 10.1190/1.1441596

  • CrossRef
  • Google Scholar

• 8

  Baser T. Nawaz K. Chung A. Faysal S. Numanoglu O. A. (2023). Ground movement patterns and shallow foundation performance in Iskenderun coastline during the 2023 Kahramanmaras earthquake sequence. Earthq. Eng. Eng. Vib.22, 867–881. 10.1007/s11803-023-2205-9

  • CrossRef
  • Google Scholar

• 9

  Bassal P. C. Boulanger R. W. (2021). System response of an interlayered deposit with spatially preferential liquefaction manifestations. J. Geotechnical Geoenvironmental Eng.147, 05021013. 10.1061/(asce)gt.1943-5606.0002684

  • CrossRef
  • Google Scholar

• 10

  Bassal P. C. Boulanger R. W. (2023). System response of an interlayered deposit with a localized graben deformation in the Northridge earthquake. Soil Dyn. Earthq. Eng.165, 107668. 10.1016/j.soildyn.2022.107668

  • CrossRef
  • Google Scholar

• 11

  Bensen G. Ritzwoller M. Barmin M. Levshin A. L. Lin F. Moschetti M. et al (2007). Processing seismic ambient noise data to obtain reliable broad-band surface wave dispersion measurements. Geophys. J. Int.169, 1239–1260. 10.1111/j.1365-246x.2007.03374.x

  • CrossRef
  • Google Scholar

• 12

  Biot M. A. (1956). Theory of propagation of elastic waves in a fluid-saturated porous solid. II. Higher frequency range. J. Acoust. Soc. Am.28, 179–191. 10.1121/1.1908241

  • CrossRef
  • Google Scholar

• 13

  Bonnefoy-Claudet S. Cotton F. Bard P.-Y. (2006). The nature of noise wavefield and its applications for site effects studies: a literature review. Earth-Science Rev.79, 205–227. 10.1016/j.earscirev.2006.07.004

  • CrossRef
  • Google Scholar

• 14

  Boulanger R. W. Idriss I. (2014). CPT and SPT based liquefaction triggering procedures. Report No. UCD/CGM.-14 1.

  • Google Scholar

• 15

  Boulanger R. W. Meyers M. W. Mejia L. H. Idriss I. M. (1998). Behavior of a fine-grained soil during the Loma Prieta earthquake. Can. Geotechnical J.35, 146–158. 10.1139/t97-078

  • CrossRef
  • Google Scholar

• 16

  Boulanger R. W. Montgomery J. (2016). Nonlinear deformation analyses of an embankment dam on a spatially variable liquefiable deposit. Soil Dyn. Earthq. Eng.91, 222–233. 10.1016/j.soildyn.2016.07.027

  • CrossRef
  • Google Scholar

• 17

  Bray J. Arduino P. Ashford S. Asimaki D. Eldridge T. Frost J. et al (2010). Geo-Engineering reconnaissance of the 2010 Maule. Chile Earthq.1–347.

  • Google Scholar

• 18

  Bray J. Rollins K. Hutchinson T. Verdugo R. Ledezma C. Mylonakis G. et al (2012). Effects of ground failure on buildings, ports, and industrial facilities. Earthq. Spectra28, 97–118. 10.1193/1.4000034

  • CrossRef
  • Google Scholar

• 19

  Bray J. D. Luque R. (2017). Seismic performance of a building affected by moderate liquefaction during the Christchurch earthquake. Soil Dyn. Earthq. Eng.102, 99–111. 10.1016/j.soildyn.2017.08.011

  • CrossRef
  • Google Scholar

• 20

  Bray J. D. Sancio R. B. (2006). Assessment of the liquefaction susceptibility of fine-grained soils. J. Geotechnical Geoenvironmental Eng.132, 1165–1177. 10.1061/(asce)1090-0241(2006)132:9(1165)

  • CrossRef
  • Google Scholar

• 21

  Campos J. Hatzfeld D. Madariaga R. López G. Kausel E. Zollo A. et al (2002). A seismological study of the 1835 seismic gap in south-central Chile. Phys. Earth Planet. Interiors132, 177–195. 10.1016/s0031-9201(02)00051-1

  • CrossRef
  • Google Scholar

• 22

  Chávez-García F. J. Luzón F. (2005). On the correlation of seismic microtremors. J. Geophys. Res. Solid Earth110. 10.1029/2005jb003671

  • CrossRef
  • Google Scholar

• 23

  Chávez-García F. J. Rodríguez M. Stephenson W. R. (2005). An alternative approach to the SPAC analysis of microtremors: exploiting stationarity of noise. Bull. Seismol. Soc. Am.95, 277–293. 10.1785/0120030179

  • CrossRef
  • Google Scholar

• 24

  Cox B. R. Teague D. P. (2016). Layering ratios: a systematic approach to the inversion of surface wave data in the absence of a priori information. Geophys. J. Int.207, 422–438. 10.1093/gji/ggw282

  • CrossRef
  • Google Scholar

• 25

  Cubrinovski M. Rhodes A. Ntritsos N. Van Ballegooy S. (2019). System response of liquefiable deposits. Soil Dyn. Earthq. Eng.124, 212–229. 10.1016/j.soildyn.2018.05.013

  • CrossRef
  • Google Scholar

• 26

  Curtis A. Gerstoft P. Sato H. Snieder R. Wapenaar K. (2006). Seismic interferometry—turning noise into signal. Lead. Edge25, 1082–1092. 10.1190/1.2349814

  • CrossRef
  • Google Scholar

• 27

  Dashti S. Bray J. D. (2013). Numerical simulation of building response on liquefiable sand. J. Geotechnical Geoenvironmental Eng.139, 1235–1249. 10.1061/(asce)gt.1943-5606.0000853

  • CrossRef
  • Google Scholar

• 28

  de la Torre C. A. Bradley B. A. McGann C. R. (2022a). 2d Geotechnical site-response analysis including soil heterogeneity and wave scattering. Earthq. Spectra38, 1124–1147. 10.1177/87552930211056667

  • CrossRef
  • Google Scholar

• 29

  de la Torre C. A. Bradley B. A. Stewart J. P. McGann C. R. (2022b). Can modeling soil heterogeneity in 2D site response analyses improve predictions at vertical array sites?Earthq. Spectra, 87552930221105107. 10.1177/87552930221105107

  • CrossRef
  • Google Scholar

• 30

  Emergeo Working Group (2013). Liquefaction phenomena associated with the Emilia earthquake sequence of May June 2012. Available at: https://nhess.copernicus.org/articles/13/935/2013/ (Northern Italy: Natural Hazards and Earth System Sciences), 13 (4), 935–947. 10.5194/nhess-13-935-2013

  • CrossRef
  • Google Scholar

• 31

  Ekström G. (2014). Love and Rayleigh phase-velocity maps, 5–40 s, of the western and central USA from USArray data. Earth Planet. Sci. Lett.402, 42–49. 10.1016/j.epsl.2013.11.022

  • CrossRef
  • Google Scholar

• 32

  Ekström G. Abers G. A. Webb S. C. (2009). Determination of surface-wave phase velocities across USArray from noise and Aki’s spectral formulation. Geophys. Res. Lett.36. 10.1029/2009gl039131

  • CrossRef
  • Google Scholar

• 33

  El Haber E. Cornou C. Jongmans D. Abdelmassih D. Y. Lopez-Caballero F. Al-Bittar T. (2019). Influence of 2D heterogeneous elastic soil properties on surface ground motion spatial variability. Soil Dyn. Earthq. Eng.123, 75–90. 10.1016/j.soildyn.2019.04.014

  • CrossRef
  • Google Scholar

• 34

  Galli C. (1967). Geología urbana y suelo de fundación de concepción y talcahuano, Chile. Universidad de Concepción.

  • Google Scholar

• 35

  Garofalo F. Foti S. Hollender F. Bard P. Cornou C. Cox B. et al (2016). InterPACIFIC project: comparison of invasive and non-invasive methods for seismic site characterization. Part II: inter-comparison between surface-wave and borehole methods. Soil Dyn. Earthq. Eng.82, 241–254. 10.1016/j.soildyn.2015.12.009

  • CrossRef
  • Google Scholar

• 36

  Green R. A. Cubrinovski M. Cox B. Wood C. Wotherspoon L. Bradley B. et al (2014). Select liquefaction case histories from the 2010–2011 Canterbury earthquake sequence. Earthq. Spectra30, 131–153. 10.1193/030713eqs066m

  • CrossRef
  • Google Scholar

• 37

  Hallal M. M. Cox B. R. (2023). What spatial area influences seismic site response: insights gained from multiazimuthal 2D ground response analyses at the treasure island downhole array. J. Geotechnical Geoenvironmental Eng.149, 04022124. 10.1061/jggefk.gteng-11023

  • CrossRef
  • Google Scholar

• 38

  Harkrider D. G. (1964). Surface waves in multilayered elastic media I. Rayleigh and Love waves from buried sources in a multilayered elastic half-space. Bull. Seismol. Soc. Am.54, 627–679. 10.1785/bssa0540020627

  • CrossRef
  • Google Scholar

• 39

  Hashash Y. M. Park D. (2001). Non-linear one-dimensional seismic ground motion propagation in the Mississippi embayment. Eng. Geol.62, 185–206. 10.1016/s0013-7952(01)00061-8

  • CrossRef
  • Google Scholar

• 40

  Hashash Y. M. Park D. (2002). Viscous damping formulation and high frequency motion propagation in non-linear site response analysis. Soil Dyn. Earthq. Eng.22, 611–624. 10.1016/s0267-7261(02)00042-8

  • CrossRef
  • Google Scholar

• 41

  Hu J. (2023). Empirical relationships between earthquake magnitude and maximum distance based on the extended global liquefaction-induced damage cases. Acta Geotech.18, 2081–2095. 10.1007/s11440-022-01637-y

  • CrossRef
  • Google Scholar

• 42

  Huang D. Wang G. Du C. Jin F. (2021). Seismic amplification of soil ground with spatially varying shear wave velocity using 2D spectral element method. J. Earthq. Eng.25, 2834–2849. 10.1080/13632469.2019.1654946

  • CrossRef
  • Google Scholar

• 43

  Hutabarat D. Bray J. D. (2021). Effective stress analysis of liquefiable sites to estimate the severity of sediment ejecta. J. Geotechnical Geoenvironmental Eng.147, 04021024. 10.1061/(asce)gt.1943-5606.0002503

  • CrossRef
  • Google Scholar

• 44

  Inzunza D. A. Montalva G. A. Leyton F. Prieto G. Ruiz S. (2019). Shallow ambient-noise 3D tomography in the Concepción basin, Chile: implications for low-frequency ground motions. Bull. Seismol. Soc. Am.109, 75–86. 10.1785/0120180061

  • CrossRef
  • Google Scholar

• 45

  Ishihara K. (1993). Liquefaction and flow failure during earthquakes. Geotechnique43, 351–451. 10.1680/geot.1993.43.3.351

  • CrossRef
  • Google Scholar

• 46

  Joyner W. B. Chen A. T. (1975). Calculation of nonlinear ground response in earthquakes. Bull. Seismol. Soc. Am.65, 1315–1336. 10.1785/BSSA0650051315

  • CrossRef
  • Google Scholar

• 47

  Kayen R. Moss R. Thompson E. Seed R. Cetin K. Der Kiureghian A. et al (2013). Shear-wave velocity–based probabilistic and deterministic assessment of seismic soil liquefaction potential. J. Geotechnical Geoenvironmental Eng.139, 407–419. 10.1061/(asce)gt.1943-5606.0000743

  • CrossRef
  • Google Scholar

• 48

  Khosravifar A. Elgamal A. Lu J. Li J. (2018). A 3D model for earthquake-induced liquefaction triggering and post-liquefaction response. Soil Dyn. Earthq. Eng.110, 43–52. 10.1016/j.soildyn.2018.04.008

  • CrossRef
  • Google Scholar

• 49

  Kissling E. (1988). Geotomography with local earthquake data. Rev. Geophys.26, 659–698. 10.1029/rg026i004p00659

  • CrossRef
  • Google Scholar

• 50

  Lobkis O. I. Weaver R. L. (2001). On the emergence of the green’s function in the correlations of a diffuse field. J. Acoust. Soc. Am.110, 3011–3017. 10.1121/1.1417528

  • CrossRef
  • Google Scholar

• 51

  Luque R. Bray J. D. (2017). Dynamic analyses of two buildings founded on liquefiable soils during the Canterbury earthquake sequence. J. Geotechnical Geoenvironmental Eng.143, 04017067. 10.1061/(asce)gt.1943-5606.0001736

  • CrossRef
  • Google Scholar

• 52

  Luque R. Bray J. D. (2020). Dynamic soil-structure interaction analyses of two important structures affected by liquefaction during the Canterbury earthquake sequence. Soil Dyn. Earthq. Eng.133, 106026. 10.1016/j.soildyn.2019.106026

  • CrossRef
  • Google Scholar

• 53

  Lysmer J. Kuhlemeyer R. L. (1969). Finite dynamic model for infinite media. J. Eng. Mech. Div.95, 859–877. 10.1061/jmcea3.0001144

  • CrossRef
  • Google Scholar

• 54

  Ma Y. Clayton R. W. (2014). The crust and uppermost mantle structure of southern Peru from ambient noise and earthquake surface wave analysis. Earth Planet. Sci. Lett.395, 61–70. 10.1016/j.epsl.2014.03.013

  • CrossRef
  • Google Scholar

• 55

  McGann C. R. Arduino P. Mackenzie-Helnwein P. (2012). Stabilized single-point 4-node quadrilateral element for dynamic analysis of fluid saturated porous media. Acta Geotech.7, 297–311. 10.1007/s11440-012-0168-5

  • CrossRef
  • Google Scholar

• 56

  McKenna F. (2011). OpenSees: a framework for earthquake engineering simulation. Comput. Sci. Eng.13, 58–66. 10.1109/mcse.2011.66

  • CrossRef
  • Google Scholar

• 57

  Menke W. Jin G. (2015). Waveform fitting of cross spectra to determine phase velocity using Aki’s formula. Bull. Seismol. Soc. Am.105, 1619–1627. 10.1785/0120140245

  • CrossRef
  • Google Scholar

• 58

  Montalva G. Ruz F. Escribano D. Bastías N. Espinoza D. Paredes F. (2022). Chilean liquefaction case history database. Earthq. Spectra38, 2260–2280. 10.1177/87552930211070313

  • CrossRef
  • Google Scholar

• 59

  Montalva G. A. Bastías N. Rodriguez-Marek A. (2017). Ground-motion prediction equation for the Chilean subduction zone. Bull. Seismol. Soc. Am.107, 901–911. 10.1785/0120160221

  • CrossRef
  • Google Scholar

• 60

  Montalva G. A. Chávez-Garcia F. J. Tassara A. Jara Weisser D. M. (2016). Site effects and building damage characterization in Concepción after the Mw 8.8 Maule earthquake. Earthq. Spectra32, 1469–1488. 10.1193/101514eqs158m

  • CrossRef
  • Google Scholar

• 61

  Montgomery J. Boulanger R. W. (2017). Effects of spatial variability on liquefaction-induced settlement and lateral spreading. J. Geotechnical Geoenvironmental Eng.143, 04016086. 10.1061/(asce)gt.1943-5606.0001584

  • CrossRef
  • Google Scholar

• 62

  Moreno M. Rosenau M. Oncken O. (2010). 2010 Maule earthquake slip correlates with pre-seismic locking of Andean subduction zone. Nature467, 198–202. 10.1038/nature09349

  • CrossRef
  • Google Scholar

• 63

  National Academies of Sciences, E., Medicine (2016). State of the art and practice in the assessment of earthquake-induced soil liquefaction and its consequences.

  • Google Scholar

• 64

  Nelder J. A. Mead R. (1965). A simplex method for function minimization. Comput. J.7, 308–313. 10.1093/comjnl/7.4.308

  • CrossRef
  • Google Scholar

• 65

  Ohori M. Nobata A. Wakamatsu K. (2002). A comparison of ESAC and FK methods of estimating phase velocity using arbitrarily shaped microtremor arrays. Bull. Seismol. Soc. Am.92, 2323–2332. 10.1785/0119980109

  • CrossRef
  • Google Scholar

• 66

  Olivar-Castaño A. Pilz M. Pedreira D. Pulgar J. Díaz-González A. González-Cortina J. M. (2020). Regional crustal imaging by inversion of multimode Rayleigh wave dispersion curves measured from seismic noise: application to the Basque-cantabrian zone (N Spain). J. Geophys. Res. Solid Earth125, e2020JB019559. 10.1029/2020jb019559

  • CrossRef
  • Google Scholar

• 67

  Parolai S. Picozzi M. Richwalski S. Milkereit C. (2005). Joint inversion of phase velocity dispersion and H/V ratio curves from seismic noise recordings using a genetic algorithm, considering higher modes. Geophys. Res. Lett.32. 10.1029/2004gl021115

  • CrossRef
  • Google Scholar

• 68

  Parra-Colmenares E. J. (1996). Numerical modeling of liquefaction and lateral ground deformation including cyclic mobility and dilation response in soil systems. Rensselaer Polytechnic Institute.

  • Google Scholar

• 69

  Petracca M. Candeloro F. Camata G. (2017). STKO user manual. Pescara, Italy: ASDEA Software Technology, 551.

  • Google Scholar

• 70

  Picozzi M. Parolai S. Bindi D. Strollo A. (2009). Characterization of shallow geology by high-frequency seismic noise tomography. Geophys. J. Int.176, 164–174. 10.1111/j.1365-246x.2008.03966.x

  • CrossRef
  • Google Scholar

• 71

  Pilz M. Cotton F. (2019). Does the one-dimensional assumption hold for site response analysis? A study of seismic site responses and implication for ground motion assessment using KiK-net strong-motion data. Earthq. Spectra35, 883–905. 10.1193/050718eqs113m

  • CrossRef
  • Google Scholar

• 72

  Pilz M. Cotton F. Zhu C. (2021). How much are sites affected by two-and three-dimensional site effects? A study based on single-station earthquake records and implications for ground motion modelling. Geophys. J. Int.228, 1992–2004. 10.1093/gji/ggab454

  • CrossRef
  • Google Scholar

• 73

  Pilz M. Parolai S. Picozzi M. Bindi D. (2012). Three-dimensional shear wave velocity imaging by ambient seismic noise tomography. Geophys. J. Int.189, 501–512. 10.1111/j.1365-246x.2011.05340.x

  • CrossRef
  • Google Scholar

• 74

  Pilz M. Parolai S. Woith H. (2017). A 3-D algorithm based on the combined inversion of Rayleigh and Love waves for imaging and monitoring of shallow structures. Geophys. J. Int.209, 152–166. 10.1093/gji/ggx005

  • CrossRef
  • Google Scholar

• 75

  Poblete M. Dobry R. (1968). Modelo dinámico del suelo de concepción. Rev. IDIEM7, 12–18.

  • Google Scholar

• 76

  Popescu R. (1995). Stochastic variability of soil properties: data analysis, digital simulation, effects on system behavior. Princeton University.

  • Google Scholar

• 77

  Popescu R. Prévost J. H. Deodatis G. (1997). Effects of spatial variability on soil liquefaction: some design recommendations. Geotechnique47, 1019–1036. 10.1680/geot.1997.47.5.1019

  • CrossRef
  • Google Scholar

• 78

  Popescu R. Prevost J.-H. Deodatis G. (2005). 3D effects in seismic liquefaction of stochastically variable soil deposits. Geotechnique55, 21–31. 10.1680/geot.2005.55.1.21

  • CrossRef
  • Google Scholar

• 79

  Popescu R. Prevost J. H. Deodatis G. Chakrabortty P. (2006). Dynamics of nonlinear porous media with applications to soil liquefaction. Soil Dyn. Earthq. Eng.26, 648–665. 10.1016/j.soildyn.2006.01.015

  • CrossRef
  • Google Scholar

• 80

  Pretell R. Ziotopoulou K. Davis C. A. (2021). Liquefaction and cyclic softening at balboa boulevard during the 1994 northridge earthquake. J. Geotechnical Geoenvironmental Eng.147, 05020014. 10.1061/(asce)gt.1943-5606.0002417

  • CrossRef
  • Google Scholar

• 81

  Qiu Z. Prabhakaran A. Elgamal A. (2023). A three-dimensional multi-surface plasticity soil model for seismically-induced liquefaction and earthquake loading applications. Acta Geotech.18, 5123–5146. 10.1007/s11440-023-01941-1

  • CrossRef
  • Google Scholar

• 82

  Ramirez J. Barrero A. R. Chen L. Dashti S. Ghofrani A. Taiebat M. et al (2018). Site response in a layered liquefiable deposit: evaluation of different numerical tools and methodologies with centrifuge experimental results. J. Geotechnical Geoenvironmental Eng.144, 04018073. 10.1061/(asce)gt.1943-5606.0001947

  • CrossRef
  • Google Scholar

• 83

  Ritzwoller M. H. Lin F.-C. Shen W. (2011). Ambient noise tomography with a large seismic array. Comptes Rendus Geosci.343, 558–570. 10.1016/j.crte.2011.03.007

  • CrossRef
  • Google Scholar

• 84

  Ruegg J. Rudloff A. Vigny C. Madariaga R. De Chabalier J. Campos J. et al (2009). Interseismic strain accumulation measured by GPS in the seismic gap between Constitución and Concepción in Chile. Phys. Earth Planet. Interiors175, 78–85. 10.1016/j.pepi.2008.02.015

  • CrossRef
  • Google Scholar

• 85

  Saldaña H. Montalva G. Escribano D. Núñez-Jara S. Tiznado J. C. (2023). Two-dimensional nonlinear dynamic analysis of a liquefaction case history considering spatial variability, and long-duration megathrust earthquakes. Soil Dyn. Earthq. Eng. (in press).

  • Google Scholar

• 86

  Saldaña Sotelo H. (2023). Modelamiento numérico de asentamientos inducidos por licuación en la zona subductiva chilena. Master thesis, University of Concepcion.

  • Google Scholar

• 87

  Seed H. B. Idriss I. M. (1971). Simplified procedure for evaluating soil liquefaction potential. J. Soil Mech. Found. Div.97, 1249–1273. 10.1061/jsfeaq.0001662

  • CrossRef
  • Google Scholar

• 88

  Shapiro N. M. Campillo M. (2004). Emergence of broadband Rayleigh waves from correlations of the ambient seismic noise. Geophys. Res. Lett.31. 10.1029/2004gl019491

  • CrossRef
  • Google Scholar

• 89

  Shapiro N. M. Campillo M. Stehly L. Ritzwoller M. H. (2005). High-resolution surface-wave tomography from ambient seismic noise. Science307, 1615–1618. 10.1126/science.1108339

  • CrossRef
  • Google Scholar

• 90

  Taftsoglou M. Valkaniotis S. Papathanassiou G. Karantanellis E. (2023). Satellite imagery for rapid detection of liquefaction surface manifestations: the case study of Türkiye–Syria 2023 earthquakes. Remote Sens.15, 4190. 10.3390/rs15174190

  • CrossRef
  • Google Scholar

• 91

  Tao Y. Rathje E. (2020). Taxonomy for evaluating the site-specific applicability of one-dimensional ground response analysis. Soil Dyn. Earthq. Eng.128, 105865. 10.1016/j.soildyn.2019.105865

  • CrossRef
  • Google Scholar

• 92

  Thompson E. M. Baise L. G. Tanaka Y. Kayen R. E. (2012). A taxonomy of site response complexity. Soil Dyn. Earthq. Eng.41, 32–43. 10.1016/j.soildyn.2012.04.005

  • CrossRef
  • Google Scholar

• 93

  Tiznado J. C. Dashti S. Ledezma C. (2021). Probabilistic predictive model for liquefaction triggering in layered sites improved with dense granular columns. J. Geotechnical Geoenvironmental Eng.147, 04021100. 10.1061/(asce)gt.1943-5606.0002609

  • CrossRef
  • Google Scholar

• 94

  Tokimatsu K. Tamura S. Kojima H. (1992). Effects of multiple modes on Rayleigh wave dispersion characteristics. J. geotechnical Eng.118, 1529–1543. 10.1061/(asce)0733-9410(1992)118:10(1529)

  • CrossRef
  • Google Scholar

• 95

  Tokimatsu K. Uchida A. (1990). Correlation between liquefaction resistance and shear wave velocity. Soils Found.30, 33–42. 10.3208/sandf1972.30.2_33

  • CrossRef
  • Google Scholar

• 96

  Tsai V. C. Moschetti M. P. (2010). An explicit relationship between time-domain noise correlation and spatial autocorrelation (SPAC) results. Geophys. J. Int.182, 454–460. 10.1111/j.1365-246x.2010.04633.x

  • CrossRef
  • Google Scholar

• 97

  Van Ballegooy S. Malan P. Lacrosse V. Jacka M. Cubrinovski M. Bray J. et al (2014). Assessment of liquefaction-induced land damage for residential Christchurch. Earthq. Spectra30, 31–55. 10.1193/031813eqs070m

  • CrossRef
  • Google Scholar

• 98

  Vantassel J. P. Cox B. R. (2021). SWinvert: a workflow for performing rigorous 1-D surface wave inversions. Geophys. J. Int.224, 1141–1156. 10.1093/gji/ggaa426

  • CrossRef
  • Google Scholar

• 99

  Verdugo R. González J. (2015). Liquefaction-induced ground damages during the 2010 Chile earthquake. Soil Dyn. Earthq. Eng.79, 280–295. 10.1016/j.soildyn.2015.04.016

  • CrossRef
  • Google Scholar

• 100

  Vivallos J. Ramírez P. Fonseca A. (2010). Microzonificación sísmica de la ciudad de concepción, región del biobío. Serv. Nac. Geol. Minería. Carta Geol. Chile, Ser. Geol. Ambient.12.

  • Google Scholar

• 101

  Wapenaar K. Slob E. Snieder R. Curtis A. (2010). Tutorial on seismic interferometry: Part 2—underlying theory and new advances. Geophysics75, 75A211–75A227. 10.1190/1.3463440

  • CrossRef
  • Google Scholar

• 102

  Ward K. M. Porter R. C. Zandt G. Beck S. L. Wagner L. S. Minaya E. et al (2013). Ambient noise tomography across the central andes. Geophys. J. Int.194, 1559–1573. 10.1093/gji/ggt166

  • CrossRef
  • Google Scholar

• 103

  Weaver R. L. (1982). On diffuse waves in solid media. J. Acoust. Soc. Am.71, 1608–1609. 10.1121/1.387816

  • CrossRef
  • Google Scholar

• 104

  Welch P. (1967). The use of fast Fourier transform for the estimation of power spectra: a method based on time averaging over short, modified periodograms. IEEE Trans. audio electroacoustics15, 70–73. 10.1109/tau.1967.1161901

  • CrossRef
  • Google Scholar

• 105

  Xu Y. Lebedev S. Meier T. Bonadio R. Bean C. J. (2021). Optimized workflows for high-frequency seismic interferometry using dense arrays. Geophys. J. Int.227, 875–897. 10.1093/gji/ggab260

  • CrossRef
  • Google Scholar

• 106

  Yamanaka H. Ishida H. (1996). Application of genetic algorithms to an inversion of surface-wave dispersion data. Bull. Seismol. Soc. Am.86, 436–444. 10.1785/bssa0860020436

  • CrossRef
  • Google Scholar

• 107

  Yang Z. Elgamal A. Parra E. (2003). Computational model for cyclic mobility and associated shear deformation. J. Geotechnical Geoenvironmental Eng.129, 1119–1127. 10.1061/(asce)1090-0241(2003)129:12(1119)

  • CrossRef
  • Google Scholar

• 108

  Youd T. L. Idriss I. M. (2001). Liquefaction resistance of soils: summary report from the 1996 NCEER and 1998 NCEER/NSF workshops on evaluation of liquefaction resistance of soils. J. geotechnical geoenvironmental Eng.127, 297–313. 10.1061/(asce)1090-0241(2001)127:4(297)

  • CrossRef
  • Google Scholar

• 109

  Zhang G. Robertson P. Brachman R. W. (2002). Estimating liquefaction-induced ground settlements from CPT for level ground. Can. Geotechnical J.39, 1168–1180. 10.1139/t02-047

  • CrossRef
  • Google Scholar

• 110

  Zhang J. Wang Y. (2021). An ensemble method to improve prediction of earthquake-induced soil liquefaction: a multi-dataset study. Neural Comput. Appl.33, 1533–1546. 10.1007/s00521-020-05084-2

  • CrossRef
  • Google Scholar

• 111

  Zhang Y. Xie Y. Zhang Y. Qiu J. Wu S. (2021). The adoption of deep neural network (DNN) to the prediction of soil liquefaction based on shear wave velocity. Bull. Eng. Geol. Environ.80, 5053–5060. 10.1007/s10064-021-02250-1

  • CrossRef
  • Google Scholar

• 112

  Zhou Y.-G. Chen Y.-M. (2007). Laboratory investigation on assessing liquefaction resistance of sandy soils by shear wave velocity. J. Geotechnical Geoenvironmental Eng.133, 959–972. 10.1061/(asce)1090-0241(2007)133:8(959)

  • CrossRef
  • Google Scholar

• 113

  Zienkiewicz O. Shiomi T. (1984). Dynamic behaviour of saturated porous media; the generalized Biot formulation and its numerical solution. Int. J. Numer. Anal. methods geomechanics8, 71–96. 10.1002/nag.1610080106

  • CrossRef
  • Google Scholar