Development of an Estimation Method for the Seismic Motion Reproducibility of a Three-dimensional Ground Structure Model by combining Surface-observed Seismic Motion and Three-dimensional Seismic Motion Analysis

Structures on or within the ground can experience substantial impacts from the ground structure during seismic events [1][2]. For instance, soft ground atop the hard ground, coupled with the characteristics of seismic motion and the ground structure, can lead to localized amplification of seismic motion, resulting in damage. Enhancing the reliability of ground structure assessments is required to enable accurate evaluation of structural behavior during seismic events and mitigate earthquake-induced damage effectively.

Many methods for estimating ground structures have been proposed and employed in practical applications. For instance, methods directly assessing in situ geotechnical properties, such as borehole testing, yield reliable outcomes. Nonetheless, due to the limited number of measurement points, interpolation becomes necessary to estimate the spatial extent of the ground structure. In the construction of 3D ground structure models, interpolation methods based on geotechnical engineering knowledge, such as inverse distance weighting methods [3], curvature minimization principles [4], and Kriging methods [5], are employed. These interpolations rely on various assumptions, resulting in multiple 3D ground structure models, necessitating additional evaluation to determine the reliability of each model.

Such evaluation can utilize the observed ground motion of small amplitudes. For instance, the reliability of a candidate ground structure model can be assessed by comparing the phase velocity derived from long-term microtremor observations [6] or by comparing the Green’s functions obtained via seismic interferometry [7]. Alternatively, the quality of 3D ground structure models can be evaluated by directly utilizing observed seismic ground motions of small amplitudes, which typically exhibit stronger signal strengths. For example, [8] demonstrates the feasibility of evaluating a 3D ground structure model using observed seismic ground motions of small amplitudes. However, the method in [8] assumes that seismic motions input to the 3D ground structure model are observed at underground observation points, making it challenging to apply to sites where underground observation points cannot be installed.

Building upon the above insights, we propose a method that selects a credible ground structure model from many generated 3D ground structure models using small amplitude seismic motions solely observed at the ground surface. We illustrate the effectiveness of the method through numerical experiments by conducting $10^{2}$–$10^{3}$ cases of fast 3D seismic wave propagation analysis on GPUs, and show that we can select a credible ground structure model based on the amount of seismic motion information observed at the ground surface. Additionally, we show that the selected 3D ground structure model can be utilized to evaluate ground motion with sufficient accuracy.

We introduce a method for extracting a credible 3D ground structure model from a pool of candidate 3D ground structure models using small amplitude seismic motions observed at the ground surface. The general setting, depicted on the left of Fig. 1, involves inputting seismic waves to a 3D ground structure comprising soft and hard soil, with resulting seismic ground motions observed at the designated points on the surface. Notably, the 3D ground structure and the input seismic waves are typically unknown.

Next, we elucidate the numerical analysis model for the target system, as depicted on the right side of Fig. 1. This study uses the cartesian coordinate system $x_{1},x_{2},x_{3}$ for simplicity. The ground enclosed by the dashed line in the left panel of Fig. 1 is designated as $V(\mathbf{m})$, defined by a set of model parameters $\mathbf{m}$. An external force, represented by the input velocity wave $f_{i}(t)$, is applied to the bottom of $V(\mathbf{m})$, and its dynamic response is analyzed to yield the time-history response $obs^{k}_{i}(t)$ at each observation point $k$ on the ground surface. To compute the small amplitude seismic response of the ground to an earthquake, we model $V(\mathbf{m})$ as linearly elastic and solve the governing equations for a linear dynamic elastic body

Here, $c_{ijkl}$, $u_{j}$, and $\rho$ represent the elasticity tensor, displacement in the $j$-th direction, and density, respectively. The notation $(\ddot{~{}})$ denotes second-order derivatives in time, and $d_{i}$ indicates the differentiation in the $i$-th direction. Note that stress-free, non-reflective, and semi-infinite absorbing boundary conditions are enforced on the model’s top, bottom, and side surfaces. We aim to develop a method that can evaluate $V(\mathbf{m})$ and $f_{i}(t)$ that is consistent with $obs^{k}_{i}(t)$.

We describe the parametrization employed for this evaluation. First, we suppose that the input velocity wave can be described as follows:

We utilize unit impulse waves $p(t)$, for instance, trigonometric functions, to represent the input velocity wave. Here, $c_{ij}$ represents unknown scalar values assessed through error minimization, and $\Delta t$ denotes the time-stepping stride for discretization. While arbitrary $p(t)$ and $\Delta t$ are permissible, we configure them to facilitate the reconstruction of the input wave within the target frequency range. Moreover, we denote the Green’s function of the $i$-th directional velocity response at observation point $k$ when $p(t)$ is input in the $j$-th direction at the bottom of $V(\mathbf{m})$, as $G^{k}_{ij}(\mathbf{m},t)$. The $i$-th directional velocity at observation point $k$ for $f_{j}(t)$ can thus be expressed as

Utilizing the observed ground velocity in the $i$-th direction at the $k$-th observation point ($obs^{k}_{i}(t)$), the error can be assessed as follows:

With the parameter settings outlined above, we compute $c_{ij}$ that minimizes the error from the observed seismic motion $obs^{k}_{i}(t)$ for each candidate $l$-th 3D ground structure model characterized by model parameters $\mathbf{m}_{l}$. First, we compute $G^{k}_{ij}(\mathbf{m}_{l},t)$ for the model parameters $\mathbf{m}_{l}$. Subsequently, by considering the stationary condition of $ERR$, which is a second-order function of $c_{ij}$, the coefficients $c_{ij}$ can be determined by solving:

Here, $\mathbf{A}$, $\mathbf{b}$ represent a constant matrix and a constant vector, respectively, while $\mathbf{c}$ denotes an unknown vector with components $c_{ij}$.

By applying singular value decomposition to $\mathbf{A}$ and discarding negligible singular values, we construct a pseudoinverse matrix and compute $\mathbf{c}$ and $ERR$. Note that the dimension of $\mathbf{A}$ is small because this computation can be performed independently for each trial, and thus the cost of this computation can be kept small enough. That is, if the seismic ground motions contain sufficient information (i.e., the number of observation points and the number of events is adequate) and are suitably constrained by the model-specific time-history Green’s functions, it is possible to assess the credibility of various 3D ground structure models by comparing the error $ERR$ obtained through attempts to reproduce observed seismic ground motions using $\mathbf{m}_{l}$.

We conducted a numerical experiment to demonstrate the efficacy of the proposed method in selecting a credible 3D ground structure model from among geotechnically estimated ground models, using seismic ground motions observed at the ground surface. Here, a credible ground model can reproduce the seismic ground motions observed at the surface. Specifically, the method described in Section 2 estimates the seismic ground motion that best fits the observed seismic ground motion using a given 3D ground structure model, albeit with introduced errors. This analysis in estimating seismic ground motion that best fits the observed ground motion for a given 3D ground structure model is conducted each time an earthquake is observed. Discarding ground models that inadequately reproduce seismic ground motions can identify a credible ground model capable of reproducing observed seismic ground motions at a specific site.

First, we set a reference model to serve as the correct solution for the numerical experiment. The reference model comprises a two-layer structure comprising a sedimentary layer and bedrock, which imitates a real ground structure. It spans 600, 600, and 100 m in the $x_{1}$, $x_{2}$, and $x_{3}$ directions, respectively (see Fig. 2). The ground surface is flat at an elevation of 0 m, with the thickness of the first layer specified as depicted in Fig. 7a). The physical properties of each layer are given in Fig. 2. This ground structure emulates valley floor lowlands formed by sediment accumulation from river meandering. During earthquakes, local ground amplification occurs, emphasizing the importance of accurately assessing the ground structure for evaluating seismic ground motion and making effective earthquake mitigation strategies. Assuming that the reference model is situated at KiK-net [9] station IBRH19, we simulate a scenario where earthquakes listed in Table 1 occur sequentially, and seismic motions are observed. To simulate this scenario, actual ground motions observed at underground stations of IBRH19 were inputted as forces from the bottom of the reference ground model, and pseudo-observed seismic ground motions were recorded at the ground surface. The analysis was performed up to 2.5 Hz, which is the frequency range considered to have a significant impact on structural damage. Semi-infinite absorbing boundary conditions are applied to the sides and bottom of the numerical model.

Next, a set of candidate ground structure models is generated using geotechnical engineering methods. While many methods generating ground structure models exist, we assume that the physical properties of the first and second layers, along with the boundary location between them, are determined from a simple ground survey conducted at the 120 survey points illustrated in Fig. 7a) (i.e., $\mathbf{m}$ represents model parameters describing the boundary shape between the first and second layers). The geometry of the boundary between the first and second layers within the target area of 600 $\times$ 600 m is estimated based on information regarding the boundary layer location at these 120 survey points, creating a set of candidate ground structure models. Although many methods for estimating the layer boundary geometry under these conditions exist, we employ an inverse distance weighting method. This method assigns weights to observed values obtained near the evaluation point based on the inverse of the distance, yielding an interpolated result. Utilizing the inverse distance weighting method, the layer thickness $z$ at any point in the domain is evaluated as follows:

where $\bar{z}_{i}$ denotes the layer thickness at the $i$-th measurement point located at $(\bar{x}_{1}^{i},\bar{x}_{2}^{i})$, $M$ represents the number of measurement points nearest to the target position $(x_{1},x_{2})$ utilized for interpolation, and $q$ is a positive constant parameter. The values of $M$ and $q$ are arbitrary. In this context, we designate $M=i~{}(i=1,2,...,50)$ and $q=0.1i~{}(i=1,2,....,40)$ to generate a total of $50\times 40=2000$ 3D ground structure models. Subsequently, from the pool of candidate 3D ground structure models, we select 236 models that exhibit substantial variations in the layer boundary geometry, based on a 3D ground structure model generated using commonly used parameters (i.e., $M=20$ and $q=2.0$). Note that Laplace smoothing was applied five times to mitigate steep slopes in certain areas.

Utilizing the method described in Section 2 and the pseudo-observed seismic ground motions described above, we attempt to identify credible 3D ground structure models by reproducing the observed seismic ground motions for each of the 236 obtained models. First, we consider a scenario where only one observation point is available on the ground surface at $(x_{1},x_{2})=(300,300)$ m. The analysis conducted for each candidate model is the same as the analysis using the reference model (i.e., the target frequency is set up to 2.5 Hz), with semi-infinite absorbing boundary conditions applied to the sides and bottom of the model (the same applies to the other analysis cases in this Section). Figure 3a) illustrates the errors obtained for each ground structure model across each earthquake event. Although the error varies depending on the model, it remains very small regardless of the ground structure model in the case of a single observation point (from the definition of $ERR$ in Eq. (2), $ERR=0.05$ means that the waveform is matched within an error of about 5% of its amplitude). This indicates that the waveform constraints are insufficient when the number of observation points is limited. In other words, despite the complexity of the time history Green’s function in a 3D medium, observed waveforms can be accurately reconstructed regardless of the model used by imposing errors on the estimated input waveform. Indeed, the estimated input waveform obtained using model000083 (i.e., model number 83 of the 236 candidate models), which exhibits a small error, substantially differs from the true waveform in amplitude and phase characteristics (Fig. 4). Conversely, increasing the number of observation points becomes imperative to leverage model-specific time-history constraints in Green’s functions for accurately reconstructing both observed and input waveforms. This constraint can be leveraged to diminish the ability to reproduce the observed waveforms depending on the ground structure model, thereby facilitating the selection of the appropriate ground structure model.

Next we employ nine observation points on the ground surface $(x_{1},x_{2})=(100+200i,100+200j)$ m, where $i,j=0,1,2$. Figure 5a) illustrates the error obtained for each model across each earthquake event. Compared with the single observation point case, the information derived from the pseudo-observed earthquake ground motion adequately constrains the ground model. Consequently, models exhibiting small error levels are consistently identified irrespective of the earthquake event. The comparison with the historical maximum error depicted in Fig. 5b) reveals that the credible ground model remains consistently selected even as the number of experienced earthquakes and the quantity of information increases. Furthermore, we augment the information by increasing the number of surface observation points to 25, located at $(x_{1},x_{2})=(100+100i,100+100j)$ m, where $i,j=0,1,....,4$. The error estimation is shown in Fig. 6a). As in the 9-point scenario, models with small errors are systematically chosen. The comparison with the historical maximum error exhibited in Fig. 6b) underscores the continued stable selection of the credible ground model even with the heightened number of experienced earthquakes and augmented information. Comparing the 9-point case with the 25-point case, it is evident that the increase in information facilitates the stable and systematic extraction of ground models characterized by smaller errors.

Finally, we checked the extracted results. Figures 7a,b, and c) display the layer thickness distribution of the reference ground model, model000018 (with a minimum error of 0.388 for the case with 25 observation points), and model000228 (with error of 0.781 for the case with 25 observation points). Model000018, with the smallest error, closely resembles the reference ground model, indicating the extraction of a plausible 3D ground structure model through this method. While no model among the candidate ground models exactly matches the reference model, model000018 was selected due to its minimal difference compared to the target wavelength used for analysis. Conversely, model000228, with a large error, substantially deviates from the reference model. Next, Fig. 8 illustrates the $x_{1}$ component of the input seismic motion estimated by model000018 and model000228 using pseudo-observed seismic ground motion at 25 observation points during earthquake event #1. The input waveform estimated using model000018, evaluated to have a small error, closely resembles the true input wave, while the waveform estimated using model000228 differs substantially from the true input wave. This underscores that sufficient observed seismic ground motions enable the Green’s function constraint of the time history in the 3D medium obtained from a 3D wavefield analysis, ensuring estimation performance for both observed and input seismic motions when using a plausible model. Finally, Fig. 9 depicts the ground motion distribution using the reference model, model000018, and model000228. The responses of the reference model and model000018 closely align regardless of the location from the observation points, whereas the response of model000228 differs markedly. As shown in Fig. 10, model000018, selected as the most credible 3D ground structure model in this study, can estimate not only the time series response at the observation point ($(x_{1},x_{2})=(300,300)$ m) but also the time series response at a point which is located away from the observation points ($(x_{1},x_{2})=(350,350)$ m). This indicates that a consistent ground model that matches the observed seismic ground motion has been extracted. These results underscore the significance of selecting a plausible 3D ground structure model, as demonstrated in this study, for mitigating earthquake damage, since even 3D ground structure models estimated by geotechnical engineering methods exhibit substantial differences in their performance in reproducing seismic ground motions.

As described above, it is evident that a 3D ground structure model demonstrating high performance in reproducing observed seismic motions can be extracted by using seismic motion observations at the ground surface. Achieving this necessitates multiple 3D seismic response analyses. In this study, finite element models are generated for the 236 ground models with varying layer thicknesses, and for each model, the Green’s function response of Eq. (1) to the unit impulse wave $p(t)$ in each of the $x_{1}$, $x_{2}$, and $x_{3}$ directions are computed, resulting in a total of 236$\times$3 = 708 instances of 3D wave field analysis. An automatic method for generating high-quality finite element models is required to facilitate such a large number of analyses. Moreover, a fast implicit time integration-based 3D seismic analysis method capable of performing stable calculations, even on finite element models with small local elements that are essential to faithfully model the complex ground geometry, is required. Consequently, in this study, finite element models comprising second-order tetrahedral elements are automatically generated using a mesh generator [10] employing octree background cells (the finite element models utilized in this study average 2,320,401 degrees of freedom). The Newmark-$\beta$ method, a type of implicit time integration, is employed to compute the seismic ground motion over 16,000 time steps with a time increment of $dt=0.001$ s to ensure stability. While implicit time integration is suitable for stable computation, it entails solving solutions of large sparse matrix equations, resulting in substantial computational costs. To address this issue, we use a fast solver from [11], based on the conjugate gradient method with variable preconditioning [12]. Consequently, even a 3D seismic analysis with 2.32 million degrees of freedom (516748 second-order tetrahedral elements with minimum element size of 5 m) and 16,000 time steps can be executed in approximately 270 s on a computing environment with an NVIDIA A100 PCIe 40 GB GPU [13]. The entire evaluation can be completed stably within a short duration of approximately 1.5 h using 44 A100 GPUs. Due to time constraints, conducting 3D wave propagation analysis for numerous ground models of this scale would have been impractical in the past. However, leveraging a fast 3D wave propagation method, as demonstrated here, enables such analyses to be conducted within feasible timeframes.

This study introduced a method to extract a 3D ground structure model that can reproduce observed seismic ground motions from a pool of candidate ground structure models using surface-observed earthquake ground motions and 3D seismic ground motion analysis. Through numerical experiments, we illustrated its effectiveness. Even when ground models are generated using methods based on geotechnical engineering aspects, their performance in reproducing observed ground motion varies substantially. A comparison of simulation results between models exhibiting high and low performance underscored substantial differences in seismic damage estimation during earthquakes, highlighting the efficacy of our method. While our numerical experiments focused solely on treating the layer boundary geometry as a model parameter, there are no constraints on parameterization in this context. Both the material properties and geometry can serve as model parameters for extracting a plausible 3D ground structure model. The above findings indicate that our proposed method is expected to contribute to enhancing the reliability of 3D ground structure models generated by supplementing borehole data and other databases, which is expected to become more common in the future. There is also possibility in enhancing the method by utilizing further HPC and AI for efficiently generating better candidate models around the ground models that are chosen to be credible and searching for a ground model with better performance in reproducing observed ground motions.