Conceptio › Archive › arXiv CS

arXiv CSopen access

Seismic Site Response Prediction from Sparse Observations Using Finite-Element-Pretrained Latent Dynamics

· arxiv_cs

arXiv CS · Papers · License: Open Access

Open Source ↗Direct PDF ↓

neural-networks

machine learning, deep learning, neural networks

arXiv:2609.20451v1 [cs.LG] 17 Sep 2026

S EISMIC S ITE R ESPONSE P REDICTION FROM S PARSE O BSERVATIONS U SING F INITE -E LEMENT-P RETRAINED L ATENT DYNAMICS

Yi Zhu Su Chen∗ Xiaojun Li State Key Laboratory of Bridge Safety and Resilience, Beijing University of Technology Beijing 100124, China

A BSTRACT Numerical site-response predictions often deviate from observations, yet correcting these discrepancies is difficult because records are limited in both sensor coverage and number of events. This study proposes the Transfer-Enabled Forced Latent Autoencoder for Response Equations (FLARET) to improve these predictions by learning and calibrating low-dimensional latent dynamics that connect the base acceleration input to acceleration outputs at multiple depths. FLARE-T learns a low-dimensional response manifold and input-driven dynamics from dense finite-element simulations. It then trains a sparse encoder to map simulated sensor responses into the learned coordinates and uses limited records to calibrate the dynamics within them. A short response window initializes each prediction, while the complete base motion drives the response. The framework was evaluated using a layered-soil centrifuge test and the Lotung field vertical array. Test-set results show that FLARE-T improved multi-depth acceleration histories and 5%-damped pseudoacceleration response spectra relative to the original finite-element models, reducing errors at every evaluated sensor for motions of different intensities and, at Lotung, for both horizontal components. Two Lotung source models with different constitutive parameters achieved comparable test-set accuracy, indicating reduced dependence on precise prior calibration. FLARE-T therefore provides a data-efficient means of combining dense numerical response information with limited field records to improve future site-response predictions. Keywords Seismic site response · Transfer learning · Machine learning · Reduced-order modeling

Introduction

Seismic site response analysis is routinely used in engineering practice to estimate earthquake motions at the ground surface and at selected depths within a soil profile, providing essential inputs for the seismic design and assessment of structures, foundations, and geotechnical systems (Seed et al., 1976; Borcherdt, 1994; Bazzurro and Cornell, 2004). For a prescribed input motion, the predicted response is governed jointly by the soil stratigraphy, shear-wave velocity profile, cyclic stress–strain behavior of the soil, and the amplitude, frequency content, and duration of the input motion (Rathje et al., 2010; Puzrin et al., 1997; Griffiths et al., 2016; Kim et al., 2016). In practice, these factors cannot be characterized exactly because subsurface investigations provide limited descriptions of the site, laboratory measurements may not fully represent in-situ soil behavior, and the motions selected as input represent only a finite range of possible earthquake loading (Rathje et al., 2010; Griffiths et al., 2016; Hallal et al., 2022). Additional uncertainty arises from the assumed initial state, constitutive relations, and boundary conditions used to represent the physical system (Chen, 1985; Kwok et al., 2007; Rodriguez-Marek et al., 2021). Consequently, a numerical model with a reasonable physical basis does not necessarily reproduce the response of a particular site. Comparisons with vertical-array recordings have shown that predicted acceleration time histories, amplification characteristics, and response spectra can differ appreciably from their observed counterparts, with the magnitude and frequency dependence of the discrepancies varying among sites and earthquake motions (Yee et al., 2013; Tao and Rathje, 2019; Zalachoris and Rathje, 2015; Stewart and Afshari, ∗

Corresponding author. E-mail address: [email protected]

2021; Zhu et al., 2022). Reducing this persistent discrepancy between numerical predictions and field observations therefore remains a central practical challenge in seismic site response analysis. Efforts to reconcile numerical predictions with observed site response have followed several complementary paths. Downhole array recordings have been used to infer in situ dynamic soil properties, calibrate constitutive parameters, and update numerical models against field measurements (Chang et al., 1996; Tsai and Hashash, 2008; Assimaki et al., 2011; Roten et al., 2014). More recently, advances in machine learning have enabled surrogate models to be trained on large ensembles of numerical results, providing efficient predictions of spectral amplification and acceleration time histories at multiple depths without repeated numerical analysis (Lee et al., 2023; Van Nguyen et al., 2024; Ilhan et al., 2025). Although these models can reproduce complex response patterns represented in the training simulations, agreement with numerical targets does not establish predictive accuracy against field observations because the simulated labels also reflect the assumptions and discrepancies of the underlying numerical model (Ilhan et al., 2025; Zhu et al., 2023). Other studies have trained machine-learning models directly on recorded ground motions, particularly surface and borehole records from the KiK-net network, to predict site amplification, response spectra, or acceleration time histories (Kim et al., 2020; Bergamo et al., 2021; Roten and Olsen, 2021; Zhu et al., 2023; Li et al., 2023). These observation-based models reduce their dependence on simulated response targets, but their development and transferability generally require datasets spanning numerous recording stations and earthquake events. Such data requirements are seldom satisfied in project-specific applications, where observations are sparse both spatially and in sample size: instruments are commonly installed at only a few depths, and the available dataset may contain only a limited number of usable events, particularly at strong shaking levels (Zhu et al., 2023, 2026a). Hybrid physics–data models and transfer-learning strategies have consequently been investigated as means of combining information derived from simulations, physical models, or large recording networks with limited observations from a target site (Tsai and Hashash, 2008; Zhang et al., 2025; Chen et al., 2025; Li et al., 2025). Against this background, the central problem addressed in this study is how to retain the site-specific response information contained in an existing numerical model while using limited field observations to correct its discrepancies and improve predictions for subsequent earthquake events. Conventional numerical site-response analysis computes ground motions by solving the equations of motion together with soil constitutive relations. In the present study, the soil profile is represented as a forced input–output system, with the acceleration at the reference depth as its input and the accelerations at selected locations above that depth as its outputs. These outputs are dynamically related because they arise from the propagation of the same input motion through the soil profile. System identification studies using vertical-array records have described site response through low-order dynamical models (Glaser and Baise, 2000), while reduced-order site models have reproduced downhole motions and surface response spectra using a limited number of generalized variables (Bantis et al., 2025). These findings motivate the assumption adopted here that the dominant acceleration patterns across depth can be represented on a low-dimensional response manifold parameterized by a small set of latent variables. Each latent state corresponds to a depth-dependent acceleration pattern, and its evolution under the base motion generates the response histories. The modeling task is therefore to learn the input-driven evolution of this low-dimensional state and the mapping that reconstructs the high-dimensional acceleration outputs. The Forced Latent Autoencoder for Response Equations (FLARE) proposed by Zhu et al. (2026b) provides a framework for learning this representation. It encodes high-dimensional response histories into a small set of latent variables, identifies sparse evolution equations that explicitly incorporate external forcing, and decodes the integrated latent trajectories into observable responses. To extend this capability to data-limited site response prediction, this study proposes the Transfer-Enabled Forced Latent Autoencoder for Response Equations (FLARE-T), a transfer framework that carries the response information learned from densely sampled numerical simulations into the sparse observation space and subsequently uses limited records to correct simulation-based predictions. FLARE-T is constructed through three sequential stages. First, densely sampled acceleration responses generated by a numerical site model are used to train a source FLARE model, in which the encoder, forced latent response equations, and decoder jointly characterize the propagation of input motions through the soil profile. Second, simulated responses extracted at the instrumented depths are used to train a sparse-response encoder that maps the available sensor measurements onto the latent state learned from the dense response field, thereby establishing a consistent latent representation between the dense numerical model and the sparse sensor configuration. Third, the transferred latent dynamics and observation-domain response reconstruction are calibrated against a limited set of recorded events to correct systematic discrepancies between simulated and measured responses while retaining the response structure acquired from the source model. The framework is evaluated using multi-depth acceleration records from geotechnical centrifuge tests (Afacan et al., 2014) and the field vertical array at the Lotung Large-Scale Seismic Test site (Borja et al., 1999). Events excluded from observation-based calibration are reserved to determine whether the corrections inferred from limited records remain predictive for subsequent events, with performance evaluated using acceleration time histories and 5%-damped response spectra at multiple depths. Two numerical source models with different degrees of prior parameter calibration are also examined in the field application to determine how the initial fidelity of the numerical model influences transfer performance. The results demonstrate that FLARE-T provides a systematic means of integrating physics-based

numerical response information with sparse observations, thereby improving the predictive fidelity of site response models under data-limited conditions.

Methodology

2.1 Problem formulation For a given soil profile, let u(t) ∈ R denote the absolute horizontal acceleration prescribed at the reference depth and let a(t) ∈ Rm denote the absolute acceleration outputs at m selected locations above that depth. The numerical dataset contains dense response vectors ad (t) ∈ Rmd , whereas the recorded dataset contains sparse response vectors as (t) ∈ Rms at the instrumented locations, with ms ≪ md . The numerical output locations and sensor depths establish the correspondence between these two response configurations. The reduced representation uses a latent state z(t) ∈ Rr , with r ≪ md , to describe the acceleration pattern across the dense output locations. A response mapping converts this state into the physical acceleration vector, defining the low-dimensional response manifold. The latent evolution equation advances z(t) under the prescribed input u(t), and the response mapping reconstructs the corresponding acceleration outputs. Predictions at the sensor locations are then extracted from the reconstructed dense response. Thus, the latent dimension r specifies the dimension of the learned dynamical state, while md and ms specify the spatial resolution of the numerical and observed outputs. Zhu et al. (2026b) introduced the Forced Latent Autoencoder for Response Equations (FLARE) to model forced systems from high-dimensional response data. FLARE combines an encoder that estimates a latent state from a response-history window, a latent dynamical model that evolves this state under a prescribed external input, and a decoder that reconstructs the physical response. For the densely sampled numerical response, the encoding and reconstruction operations are expressed as bd (t) = D(z(t)) , z(t) = ET (HL [ad ](t)) , a (1) where HL [ad ](t) contains L consecutive response samples ending at time t, ET is the source encoder, and D is the decoder. The response history provides the temporal information needed to estimate the dynamic state, which cannot generally be determined from an instantaneous acceleration vector. Building on FLARE, this study proposes the Transfer-Enabled Forced Latent Autoencoder for Response Equations (FLARE-T) to connect densely sampled numerical responses with sparse measured responses. FLARE-T first learns the latent response representation and its input-driven dynamics from numerical simulations, then establishes a sparse encoder that maps responses at the instrumented locations into the same latent space, and finally uses a limited number of recorded events to correct the coefficients of the transferred dynamics. The numerical model therefore provides the spatially resolved response information and the initial dynamical representation, while the recorded data account for systematic discrepancies between simulated and measured motions. The latent variables are treated as learned coordinates of the response and are not assigned direct equivalence to soil properties, constitutive parameters, or individual vibration modes. The three stages of FLARE-T are described in Section 2.2. After FLARE-T has been trained and adapted, prediction of a new event requires the prescribed input motion and a short initial window of measured response. The input u(t) is provided over the complete analysis interval [0, T ] and acts as the external forcing throughout the prediction. The sparse response as (t) is provided only over [0, Tini ] and is used to estimate z(Tini ) for initialization of the first-order latent dynamics. The prediction problem is written as bs (t) = P[u(0:T ), as (0:Tini )] (t), Tini < t ≤ T, a (2) where P denotes the trained and adapted FLARE-T model. The prescribed input remains available at every integration step, whereas the measured responses are used only for initialization and are not supplied after Tini . This separation between continuous seismic input and initial response information defines the prediction setting adopted in this study and provides the basis for the offline prediction procedure presented in Section 2.3. 2.2 FLARE-T architecture The FLARE-T architecture comprises the three sequential stages illustrated in Fig. 1. In the simulation-pretraining stage, the base excitations and densely sampled finite-element responses are used to train a source encoder ET , a latent response model, and a decoder D. In the sparse-distillation stage, simulated responses are extracted at locations corresponding to the available sensors, and a sparse encoder ES is trained to reproduce the latent states obtained from the dense response field. In the real-data-calibration stage, ES and D are retained, while the coefficients of the latent response model are calibrated using a limited number of recorded events. This three-stage procedure uses numerical simulations to establish the principal response representation, sparse simulated responses to accommodate the observation configuration, and measured records to correct the latent dynamics for subsequent site-response prediction. Stage 1 establishes the numerical source model from a collection of finite-element analyses conducted under different base motions. Each analysis is treated as a complete input–response trajectory, and the dense simulation dataset is denoted by NFE DFE = {(ui (t), ad,i (t))}i=1 , (3) 3

(a) Simulation pretraining Reconstructed FE reponses

Dense FE reponses

FE model

𝐙̈

𝒁𝑻

… =

𝑬𝑻

…

𝚵

"𝑻 𝒁

…

𝑫

ODE rollout

…

Causal window

𝚯(𝐙, 𝐙̇, 𝐔)

(fixed)

𝒖

Base excitation

(b) Sparse distillation Sparse FE reponses 𝒁𝑺

Causal window

𝑬𝑺

Reuse 𝐸!

Sparse recordings

Actual Site

Reuse 𝐷

Initialize

Latent matching

𝒁𝑻

Causal window

𝑬𝑺

𝐙̈

𝚯(𝐙, 𝐙̇, 𝐔)

… =

… (fixed)

𝒖

Observed-channel fit

𝚵

…

"𝑻 𝒁

ODE rollout

𝑫

update coeffcients

Fig. 1. Three-stage FLARE-T architecture: (a) simulation pretraining, (b) sparse distillation, and (c) real-data calibration.

where NFE is the number of simulated motions. Complete trajectories are assigned to the training and validation sets before response windows are extracted, ensuring that samples generated by the same base motion do not appear in different data subsets. The source model connects the dense response history, the input-driven latent dynamics, and the reconstructed response through z = ET (HL [ad ]) ,

ż = Θ(z, u)Ξ0 ,

bd = D(z), a

(4)

where Θ is a prescribed library containing a constant term, polynomial functions of the latent coordinates, and the base-acceleration terms, and Ξ0 is the coefficient matrix identified from the finite-element dataset. The input motion is excluded from ET and enters only through the latent response equation. The response-history window is flattened and processed by a multilayer-perceptron encoder, and the decoder maps the latent state back to the dense acceleration response. During rollout, the first-order latent equation is integrated from the state estimated at the beginning of a training segment, and the resulting latent trajectory b z is decoded without further response input. The source model is trained using LFE = λrec Lrec + λroll Lroll + λlat Llat + λeq Leq + λ1 ∥Ξ0 ∥1 , (5) with Lrec = MSE(D(z), ad ) , Llat = MSE(b z, z) ,

Lroll = MSE(D(b z), ad ) , Leq = MSE(ż, Θ(z, u)Ξ0 ) .

(6)

Here, MSE denotes the mean squared error over the samples and response components included in a training segment. The encoded state is treated as a fixed target when evaluating Llat , preventing this term from moving the encoder toward a degenerate representation. The reconstruction loss preserves the response information carried by the latent coordinates, the rollout losses constrain the integrated trajectory in the latent and response spaces, and the equation loss fits the local latent evolution. The ℓ1 penalty and sequential thresholding remove small coefficients and retain a sparse latent response equation. After convergence, ET , D, Ξ0 , and the retained candidate functions constitute the numerical source model. Stage 2 transfers the response representation learned from the dense finite-element output to the sparse observation configuration. This step is required because ET expects acceleration histories from all numerical output locations and therefore cannot be applied directly to the instrumented locations. Calibrating the model immediately against measured records would combine the effects of reduced spatial observation and simulation–measurement discrepancy. FLARE-T separates these effects by first establishing the dense-to-sparse connection within the numerical domain, 4

where synchronized dense and sparse responses are available for the same motion. A fixed observation matrix C maps the dense response to the sensor configuration: as,FE (t) = Cad (t).

(7)

Each row of C either selects a coincident numerical output or interpolates between adjacent output locations. Dense and sparse responses are normalized separately using statistics from their respective training simulations, and the same statistics are retained for the following stage. Synchronized response-history windows are then supplied to the fixed source encoder and the sparse encoder. The latter is trained by minimizing Ldist = MSE[ES (HL [as,FE ]) , ET (HL [ad ])] .

(8)

Because the two response windows originate from the same finite-element trajectory and terminate at the same time, the state produced by ET provides a direct target for ES . Only ES is updated during this stage; ET , D, Θ, and Ξ0 remain fixed. The target states are treated as constants during optimization, and the checkpoint is selected using the latent-matching error for the numerical validation trajectories. The resulting sparse encoder estimates states in the coordinate system of the numerical source model using only the channels available from the physical sensors. Stage 3 uses the measured events to correct the input-driven evolution of the transferred model. The sparse encoder ES , decoder D, observation matrix C, normalization parameters, and candidate library Θ are fixed, and only the latent-equation coefficients are calibrated. This separation is necessary because a latent representation does not have a unique coordinate system: changing the encoder can alter the latent coordinates without representing a corresponding physical change in the response. Simultaneous adjustment of the encoder, decoder, and dynamical coefficients using a limited number of measured events could therefore reduce the fitting error by changing the coordinate system rather than by correcting the response evolution. Fixing ES and D preserves the latent coordinates and the response manifold learned from the dense simulations, while fixing Θ preserves the functional basis of the dynamics. Calibration is consequently restricted to adjusting the vector field on this fixed manifold. The numerical simulations retain their information on the depth-dependent response pattern, and the measured events correct the evolution of that pattern under the imposed motion. The calibrated coefficients are written as Ξ = Ξ0 + ∆Ξ,

(9)

where ∆Ξ is initialized to zero and is the only trainable quantity in this stage. Calibration uses the candidate library established in Stage 1 and introduces no additional functional forms. For each measured event, the initial latent state is estimated from the available response window. The calibrated equation is then integrated under the recorded base motion, and the decoded dense response is mapped to the instrumented locations: bs = CD(z). z0 = ES (HL [as ](t0 )) , ż = Θ(z, u)Ξ, a (10) The response and input are processed using the normalization parameters retained from the preceding stages, and the decoded response is returned to physical units before comparison with the measurements. The measured response following initialization is used as the target of the complete rollout and is not supplied repeatedly to the encoder. The coefficient correction is determined by minimizing LR = Lobs + λ∆ R(∆Ξ),

(11)

where Lobs is the mean squared rollout error after each sensor residual has been normalized by the corresponding standard deviation from the sparse simulation data. The regularization term R is calculated separately for each latent equation as the squared norm of its coefficient correction divided by the squared norm of its initial coefficient vector, and the resulting values are averaged over the latent equations. This relative penalty limits unsupported departures from the numerical dynamics while allowing equations with different coefficient magnitudes to be corrected on a comparable basis. Sequential thresholding removes small calibrated coefficients to retain a sparse model. The final coefficients are selected according to the complete-rollout error for the measured validation events, while the held-out test events are excluded from coefficient estimation and model selection. 2.3 Offline Response Prediction After completing the three-stage model development, all components of FLARE-T are fixed for prediction, including the sparse encoder, calibrated latent dynamics, and decoder. For a new earthquake event, the horizontal input motion at the reference depth is prescribed over the complete analysis interval, while acceleration measurements at the instrumented locations are provided only within a short initial window. These initial measurements are used to establish the latent state at the start of prediction, after which the model is driven solely by the prescribed input motion. The objective of offline prediction is therefore to reproduce the subsequent acceleration responses at the selected locations through continuous latent-state evolution, without further response measurements or model-parameter updates. Let t0 denote the end of the initial observation window and the starting time of the subsequent prediction. The synchronized acceleration responses recorded at the instrumented locations over this window are processed using the 5

Initial state

𝐙#̈

𝐙𝟎 𝐙̇𝟎

𝑥! (𝑡 − 𝐿: 𝑡)

𝑥! (𝑡) 𝑥" (𝑡)

…

Latent rollout

𝑫

…

𝐙̈ = 𝚯(𝐙, 𝐙̇, 𝐮)𝚵!"#$

𝐙#

…

𝚯(𝐙, 𝐙̇, 𝐮) × 𝚵!"#$

…

𝑥# (𝑡 − 𝐿: 𝑡)

𝑬𝑺

… …

𝑥" (𝑡 − 𝐿: 𝑡)

Transferred latent dynamics

𝑥# (𝑡)

External input

Fig. 2. Offline prediction using an initial response window and the prescribed input motion. same normalization parameters adopted during model development and assembled into the response history HL [as ](t0 ). The fixed sparse encoder ES maps this response history to the initial latent state z0 , which provides the initial condition required by the first-order latent dynamics. This initialization is expressed as z0 = z(t0 ) = ES (HL [as ](t0 )) .

(12)

Starting from z0 , the calibrated latent equation is integrated over the prediction interval under the prescribed input motion. At each integration step, the current latent state and the corresponding input acceleration determine the rate of latent-state evolution, and the updated state is carried forward to the next step. The resulting latent trajectory is passed through the fixed decoder D to recover the densely distributed acceleration response, after which the observation matrix C extracts the responses at the instrumented locations. The reconstructed responses are finally transformed back to physical units using the normalization parameters retained from model development. This rollout produces continuous acceleration histories from the end of the initialization window to the end of the event: ż(t) = Θ(z(t), u(t)) Ξ, z(t0 ) = z0 , bd (t) = D(z(t)) , a bs (t) = Cb a ad (t).

t ∈ (t0 , T ],

(13)

Accordingly, the initial response window and the prescribed input motion serve distinct roles in the prediction. The former determines the latent state z0 at t0 , whereas the latter drives its evolution throughout the remaining interval (t0 , T ]. The measured responses after t0 are retained for comparison with the model output, and the corresponding predictions are generated as a continuous rollout of the calibrated FLARE-T model. In the following case studies, this procedure is applied to held-out events, and its predictive performance is evaluated by comparing the resulting multidepth acceleration histories and 5%-damped pseudo-acceleration response spectra with the corresponding measured responses.

Application

This section evaluates FLARE-T through two applications that provide distinct conditions for model development and validation. The first application uses a layered-soil centrifuge test conducted under controlled shaking, with acceleration responses recorded at multiple depths, to examine whether the response information learned from numerical simulations can be transferred to the experimental sensor configuration and corrected using a limited number of physical records (Afacan et al., 2014). The second application uses the Lotung downhole array, where earthquake records from two horizontal components allow the methodology to be evaluated under field conditions and natural seismic excitation (Elgamal et al., 1995; Zeghal et al., 1995). The two applications therefore cover model-scale and field-scale soil profiles, different material and response characteristics, and different levels of observational availability. Predictive performance is assessed using held-out records through comparisons of multidepth acceleration histories and 5%-damped pseudoacceleration response spectra. The final part of this section further examines the influence of the initial numerical model by comparing FLARE-T models developed from two sets of source-model parameters for the Lotung site. 3.1 Centrifuge Test of a Layered Soil Profile 3.1.1 Test Configuration and Numerical Model The physical benchmark was selected from the soft-clay centrifuge testing program conducted using the 9-m-radius geotechnical centrifuge at the University of California, Davis, and was identified as AHA02 in the experimental database (Afacan et al., 2014). The test was performed at a centrifugal acceleration of 57.2g, for which the 49.7-cm-high soil model represented a 28.43-m-deep prototype profile. The profile comprised a 5.66-m-thick upper layer of dense Monterey sand over seven layers of San Francisco Bay Mud, with adjacent clay layers separated by 0.572-m-thick Monterey sand drainage layers. The soil model was constructed in a hinged-plate container that permitted shear deformation of the soil column, while a rigid base connected to the servo-hydraulic shaking table applied horizontal 6

Hinged-plate container

Accelerometer

depth (m)

A15 0.00 m

Upper Monterey sand (5.66 m) A14 3.37 m

A13 6.46 m

49.7 cm model = 28.43 m prototype

Bay Mud 7 (3.26 m)

A12 7.12 m

A11 10.58 m

Bay Mud 6 (3.03 m)

A10 11.10 m

A9 14.01 m

Bay Mud 5 (3.20 m)

A8 15.22 m A7 17.39 m

Bay Mud 4 (2.86 m)

A6 18.88 m

Bay Mud 3 (3.15 m)

A5 21.68 m A4 22.19 m

Bay Mud 2 (1.09 m)

A3 23.97 m A2 26.25 m

Bay Mud 1 (2.75 m) A40

A41

Rigid base / shake table Horizontal input motion Bay Mud (MH)

Monterey sand

Fig. 3. Soil profile, accelerometer locations, and base excitation of the centrifuge model.

input motions. Simultaneous acceleration measurements along the depth of the model provide a controlled physical dataset for evaluating multidepth site-response prediction. Two accelerometers mounted on the rigid base, A40 and A41, recorded the imposed horizontal motion, and their average was adopted as the reference-depth input. The response vector contained the absolute accelerations recorded by the 14 accelerometers A2–A15 in the central vertical array. A2 was located 26.25 m below the prototype ground surface, A15 was located at the surface, and the remaining sensors sampled the intervening soil layers, as illustrated in Fig. 3. The input and observed response were therefore defined as u(t) =

aA40 (t) + aA41 (t) , 2

as (t) = [aA2 (t), aA3 (t), . . . , aA15 (t)] .

(14)

The numerical model reproduced the prototype layer thicknesses as a 28.4284-m-high and 1-m-wide plane-strain soil column in Abaqus/Standard (Dassault Systèmes Simulia Corp., 2020). The column was discretized using 200 four-node reduced-integration plane-strain elements (CPE4R), with one element across the width and 200 elements distributed over the depth while preserving the experimental layer boundaries. The opposing nodes at each elevation were constrained to have equal horizontal displacement, and vertical displacement was restrained to maintain the shear-column kinematics; the ground surface remained traction free. A geostatic step established the initial stress state before each 30-s implicit dynamic analysis. The measured or simulated input motion was imposed as horizontal acceleration at the base, using an initial time increment of 0.001 s and a maximum increment of 0.005 s. Absolute horizontal accelerations were recovered at the centers of all 200 elements to form the dense numerical response. Responses at A2–A14 were obtained by linear interpolation between the adjacent numerical output locations, whereas the A15 response was taken directly from the surface nodes. Numerical and experimental trajectories were aligned to a common time origin. Dense and sparse simulated responses were normalized separately using statistics from their respective training sets, and these normalization parameters were retained when the experimental records were introduced. The Bay Mud layers were represented using the Modified Cam–Clay model with logarithmic porous elasticity and exponential isotropic hardening (Roscoe and Burland, 1968; Dassault Systèmes Simulia Corp., 2020). Let p′ denote the mean effective stress, q the deviatoric stress, p′c the preconsolidation pressure, and εpv the plastic volumetric strain, with 7

Table 1. Material and constitutive parameters of the centrifuge finite-element model. (a) Layer-specific properties Layer

Depth (m)

Model

ρ (kg m−3 )

Vs (m s−1 )

OCR

p′c0 (kPa)

Upper Monterey sand Bay Mud 7 Monterey sand 6 Bay Mud 6 Monterey sand 5 Bay Mud 5 Monterey sand 4 Bay Mud 4 Monterey sand 3 Bay Mud 3 Monterey sand 2 Bay Mud 2 Monterey sand 1 Bay Mud 1

0.00–5.66 5.66–8.92 8.92–9.50 9.50–12.53 12.53–13.10 13.10–16.30 16.30–16.87 16.87–19.73 19.73–20.31 20.31–23.45 23.45–24.02 24.02–25.11 25.11–25.68 25.68–28.43

MC MCC MC MCC MC MCC MC MCC MC MCC MC MCC MC MCC

2029 1651 2029 1662 2029 1672 2029 1733 2029 1733 2029 1733 2029 1733

130.91 78.85 169.67 89.30 181.62 102.60 192.11 128.25 200.97 138.70 209.36 144.40 213.05 148.20

– 1.658–1.663 – 1.607–1.611 – 1.570–1.573 – 1.354–1.356 – 1.351–1.353 – 1.350–1.351 – 1.348–1.349

– 1.27 – 1.28 – 1.15 – 3.31 – 2.76 – 2.47 – 2.27

0.426 0.563 0.426 0.566 0.426 0.536 0.426 0.910 0.426 0.831 0.426 0.786 0.426 0.753

– 89.56 – 117.28 – 142.86 – 435.73 – 435.73 – 435.73 – 435.73

(b) Constitutive parameters Material Bay Mud Monterey sand

c′ (kPa)

ϕ′ (◦ )

ψ (◦ )

0.18675 –

0.01737 –

1.20 –

1.00 –

– 0.30

– 1.00

– 35

– 5

compression taken as positive. The yield surface and isotropic hardening relation are written as fMCC = q 2 + M 2 p′ (p′ − p′c ) = 0, dp′c dεpv = . ′ pc λ−κ

(15)

The Monterey sand layers were represented using the Mohr–Coulomb model with linear elasticity and nonassociated plastic flow (Dassault Systèmes Simulia Corp., 2020). In terms of the major and minor principal effective stresses σ1′ and σ3′ , the yield criterion is expressed as fMC = (σ1′ − σ3′ ) − (σ1′ + σ3′ ) sin ϕ′ − 2c′ cos ϕ′ = 0,

(16)

where c′ and ϕ′ are the effective cohesion and friction angle, respectively. A dilation angle ψ = 5◦ was adopted for the plastic potential. The stiffness assigned to each layer was determined from its mass density ρi and shear-wave velocity Vs,i . The corresponding shear modulus was used directly in the logarithmic elastic formulation for Bay Mud, while the Young’s modulus of Monterey sand was obtained using a Poisson’s ratio of 0.30: 2 Gi = ρi Vs,i ,

Ei = 2 (1 + νi ) Gi .

(17)

The initial vertical effective stress increased geostatically from zero at the ground surface to approximately 218.7 kPa at the base. The initial horizontal effective stress in each layer was assigned using the corresponding at-rest earth-pressure coefficient: ′ ′ σh0,i = K0,i σv0,i . (18) The shear-wave velocities used in the numerical model were 95% of the interpreted experimental values, and the initial preconsolidation pressures of the three upper Bay Mud layers were increased by 15% during preliminary numericalmodel calibration. Rayleigh coefficients α = 1.1919 s−1 and βR = 2.8812 × 10−3 s were applied to provide 10% damping at 1.0479 and 10 Hz. The resulting material and initial-state parameters are summarized in Table 1. 3.1.2 Dense-Response Model Training and Validation The numerical dataset was generated using 100 synthetic base motions represented by nonstationary, band-limited random waveforms with varied frequency content and strong-motion duration. The target peak ground accelerations (PGAs) were stratified over the range of 0.025–0.55g, which covers the approximately 0.03–0.54g range of the centrifuge records used in this study. Following baseline correction, each waveform was multiplied by a single scale 8

factor to attain its prescribed PGA; the resulting motions had PGAs between 0.0265 and 0.5397g. Each motion was applied to the base of the finite-element model for a 30-s analysis, and the absolute horizontal accelerations at the 200 element-center locations were retained as the dense response field. Motions 001–080 were assigned to training and motions 081–100 to validation. This division was made at the level of complete input motions, so all response windows and sequence segments derived from a given motion remained within the same dataset. For the first-stage model, each response history contained 600 time samples with an integration time step of ∆t = 0.05 s. The encoder used a response window of L = 3 consecutive samples, and the 200-dimensional response field was represented by r = 8 latent variables. The candidate library contained a constant term, all linear and quadratic terms of the latent variables, and a standalone linear term for the base acceleration; sinusoidal terms and products between the base acceleration and latent variables were excluded. For optimization, each complete response history was divided into ten contiguous 60-sample segments, resulting in 800 training sequences and 200 validation sequences, while model selection was performed using the complete validation motions. The retained model minimized the validation score formed from the rollout, reconstruction, and latent-rollout errors with relative weights of 1, 0.25, and 0.05, respectively. Prediction accuracy at the jth response location was measured using the normalized root-mean-square error (NRMSE), NRMSEj =

RMS[b aj (t) − aj (t)] × 100%, RMS[aj (t)]

(19)

where aj (t) and b aj (t) are the finite-element response and model prediction, respectively. Each box in Fig. 4 represents the distribution of NRMSEj over the 200 response depths for one complete validation motion. Table 2. Data and model settings for the centrifuge application. Data source Training

Stage 1 Dense FE responses Cases 001–080

Stage 2 Dense–sparse FE pairs Cases 001–080

Validation Test Input Target Channels L r

Cases 081–100 – HL [ad ], u ad 200 3 8

Cases 081–100 – HL [as,FE ] z 14 3 8

Stage 3 Centrifuge records S006–S008, S010–S012, S015–S017, S019–S021 S013, S018 S009, S014, S022 HL [as ], u as 14 3 8

Figure 4(a) and Fig. 4(b) compare the finite-element response field and the corresponding full-record rollout for validation motion 081. The model reproduces the timing and depth-dependent propagation of the principal response bands, together with their amplitude distribution and decay after the strong-motion interval. The remaining discrepancies are concentrated near the largest response pulses and vary gradually with depth, without altering the overall response pattern. For the 20 validation motions in Fig. 4(c), the median depthwise NRMSE ranges from 3.9% to 15.8%, and the median over all motion-depth combinations is 9.6%. The relatively compact interquartile ranges for most motions indicate that the prediction accuracy is generally maintained throughout the profile, providing the dense numerical response model required for the subsequent transfer to the experimental sensor configuration. 3.1.3 Sparse-Response Transfer and Validation The source encoder developed in Stage 1 estimates the latent state from acceleration histories at all 200 finite-element output locations, whereas the centrifuge measurements are available only at accelerometers A2–A15. The source encoder therefore cannot be applied directly to the experimental records. To reconcile these different observation configurations before introducing measured data, acceleration responses at A2–A15 were extracted from each finiteelement simulation to form a sparse dataset consistent with the centrifuge instrumentation. These sparse responses were used to train the sparse encoder ES to recover the latent states produced by the source encoder ET . This stage establishes a mapping from the instrumented responses to the latent coordinates learned from the dense numerical response field. For each simulated motion, a response-history window of length L = 3 was formed from the 14 acceleration channels corresponding to A2–A15 and supplied to ES . The synchronized dense-response window was supplied to the fixed source encoder to provide the target latent state. Cases 001–080 were used for training and cases 081–100 were retained for validation, following the complete-motion allocation summarized in Table 2. Only the parameters of ES were updated by minimizing the latent-state matching loss in Eq. (8); the source encoder ET , decoder D, candidate 9

(a) FE simulation: case 081

(b) FLARE prediction: case 081

Acceleration, a (m s−2)

10 15 20

2 1 0 −1 −2

Depthwise NRMSE (%)

−3 −4

100

099

098

097

0 096

095

094

093

Time (s)

092

091

090

089

088

087

086

Time (s)

085

084

081

083

082

Depth below surface (m)

Validation case

Fig. 4. Stage-1 validation: (a) finite-element response field for case 081; (b) corresponding model rollout; and (c) depthwise NRMSE distributions for cases 081–100. (b) A9 (NRMSE = 13.7%)

(g) Validation NRMSE by output

FE target FLARE decoder

0 −1

A2 A3

−1

−2

A4 A5 A6

−2 −3

Pseudo-spectral acceleration, Sa (m s−2)

Time (s)

(d) A2: 5%-damped spectrum

(e) A9: 5%-damped spectrum

(f) A15: 5%-damped spectrum 10

A7 A8 A9 A10 A11

8 4

Decoder output

Acceleration, a (m s−2)

(a) A2 (NRMSE = 6.1%) 2

A12

A13

A14 2 A15

0 10−1

100

Period, T (s)

0 10−1

100

Period, T (s)

100

Period, T (s)

NRMSE across validation cases (%)

Fig. 5. Stage-2 validation: (a)–(c) reconstructed acceleration time histories at A2, A9, and A15; (d)–(f) corresponding 5%-damped pseudoacceleration response spectra; and (g) NRMSE distributions at A2–A15.

library Θ, and coefficient matrix Ξ0 remained fixed. The final sparse encoder was selected according to the lowest latent-state matching error over the validation motions. Figure 5 evaluates the transferred representation for a representative validation motion and for the complete validation set. Accelerometers A2, A9, and A15 represent response locations near the base, at an intermediate depth, and at the ground surface, respectively. As shown in Figs. 5(a)–5(c), the responses reconstructed through the sparse encoder and fixed decoder reproduce the principal phases, amplitudes, and duration of the finite-element targets. The agreement is closest at A2 and decreases toward A15, with NRMSE values of 6.1%, 13.7%, and 23.0% for the three illustrated responses. The corresponding 5%-damped pseudoacceleration response spectra in Figs. 5(d)–5(f) preserve the principal spectral peaks, although local amplitude differences become more evident at A9 and A15. In Fig. 5(g), each boxplot represents the NRMSE distribution for one accelerometer across validation cases 081–100. The distributions remain centered near 10% at most locations, while several upper sensors exhibit a wider range of errors among the validation motions. These results show that the sparse encoder can locate the response state within the existing latent coordinates using only the channels available from the centrifuge instrumentation. 3.1.4

Record-Based Calibration and Held-Out Prediction

Stage 3 used the measured centrifuge motions according to the allocation summarized in Table 2. Records S006–S008, S010–S012, S015–S017, and S019–S021 were used to calibrate the latent-dynamics coefficients, while S013 and S018 were used exclusively for model selection based on the validation rollout error. Records S009, S014, and S022 were excluded from both coefficient calibration and model selection and were retained for final testing. Throughout this stage, the sparse encoder ES , decoder D, candidate library Θ, and response mapping C remained fixed. Only the coefficient matrix governing the latent-state evolution was adjusted using the calibration records. 10

(b) Stage 3

−2

SINDy coefficient, ξij

SINDy library term

(a) Stage 1 1 z1 z2 z3 z4 z5 z6 z7 z8 z12 z1 z2 z1 z3 z1 z4 z1 z5 z1 z6 z1 z7 z1 z8 z22 z2 z3 z2 z4 z2 z5 z2 z6 z2 z7 z2 z8 z32 z3 z4 z3 z5 z3 z6 z3 z7 z3 z8 z42 z4 z5 z4 z6 z4 z7 z4 z8 z52 z5 z6 z5 z7 z5 z8 z62 z6 z7 z6 z8 z72 z7 z8 z82 u

−4

−6

z1̇

z2̇

z3̇

z4̇

z5̇

z6̇

z7̇

z8̇

z1̇

z2̇

Latent-state equation

z3̇

z4̇

z5̇

z6̇

Latent-state equation

z7̇

z8̇

z1̇

z2̇

z3̇

z4̇

z5̇

z6̇

z7̇

z8̇

−8

Latent-state equation

Fig. 6. Latent-dynamics coefficients: (a) Stage 1, (b) Stage 3, and (c) their difference. Because Stage 2 did not modify the latent-dynamics coefficients, the Stage 1 coefficient matrix in Fig. 6(a) also represents the model immediately before record-based calibration. The calibrated coefficients in Fig. 6(b) retain the dominant positive and negative patterns of the source model, while the difference in Fig. 6(c) shows localized changes of smaller magnitude. The cosine similarity between the two coefficient matrices is 0.98, confirming that calibration preserved their overall structure while adjusting selected terms. Since the sparse encoder and decoder remained fixed, these coefficient changes represent corrections to the latent-state evolution within the response coordinates established from the finite-element simulations. Each held-out motion was initialized using the same three-sample response window from A2–A15. The latent equation was subsequently integrated over the remaining analysis interval under the complete measured base motion, and no additional response measurements were supplied during prediction. The resulting FLARE-T responses were compared with the centrifuge measurements and with the finite-element responses obtained without record-based correction. For a prediction apj (t) at sensor j, the time-dependent normalized absolute error was defined as apj (t) − aj (t) × 100%, (20) RMS[aj (t)] where aj (t) is the measured response and apj (t) denotes either the FLARE-T or finite-element result. The pseudoacceleration response spectra were calculated with 5% damping over periods from 0.10 to 5.00 s. In Figs. 7(g)– 9(g), each boxplot contains the values of ej (t) over all time samples of one test motion at the indicated sensor. For the lower-amplitude motion S009, FLARE-T closely follows the measured phase and principal acceleration cycles at A2, A9, and A15 and reproduces the decay of the response following the main shaking interval. The finiteelement result shows increasingly evident phase and amplitude differences toward A9 and A15, whereas FLARE-T remains aligned with the measured oscillations. The principal spectral peaks at A9 and A15 are also represented more closely by FLARE-T, although the shorter-period peak at A2 remains underestimated. The sensorwise distributions in Fig. 7(g) show lower median errors and narrower interquartile ranges for FLARE-T than for the finite-element model throughout A2–A15. For the higher-amplitude motion S014, both predictions reproduce the principal acceleration pulse near the base, but their differences become more pronounced at the intermediate and surface locations. FLARE-T more closely captures the phase, peak sequence, and subsequent decay at A9 and A15. Its response spectra reproduce the dominant A2 and A9 peaks more accurately and substantially reduce the excessive spectral amplification predicted by the finite-element model at A15. Consistent reductions in the median and spread of the time-dependent error are observed at every sensor in Fig. 8(g), indicating that the improvement is maintained as the shaking amplitude increases. Motion S022 produces the largest response amplitudes among the three held-out events. FLARE-T continues to reproduce the timing of the principal wave groups and the decay of the measured motions more closely than the ej (t) =

(a) A2

(b) A9 Measured FLARE-T FE

0.2

0.4 0.2

FLARE-T FE

0.4

0.2

0.1 0.0

0.0 −0.1

−0.2

−0.4

−0.6

−0.4

−0.3

0.0

−0.2

Pseudo-spectral acceleration, Sa (m s−2)

(g) Error distribution by sensor

0.6

A7 0

Time (s)

(d) A2

(e) A9

(f) A15

0.8

2.5

1.4 1.2

0.6

Sensor

Acceleration, a (m s−2)

0.3

0.8 0.6

2.0

A11

1.5

A12 A13

1.0

A14

0.4

0.2

0.5

0.2 0.0 −1 10

A15

0.0 −1 10

100

A9 A10

1.0

0.4

0.0 −1 10

100

Period, T (s)

100

Period, T (s)

Normalized absolute error (%)

Period, T (s)

Fig. 7. Held-out motion S009: (a)–(c) acceleration time histories; (d)–(f) 5%-damped pseudoacceleration response spectra; and (g) sensorwise error distributions. (a) A2

(b) A9 Measured FLARE-T FE

(g) Error distribution by sensor

−1

−2

−3

Pseudo-spectral acceleration, Sa (m s−2)

−3

Time (s)

(d) A2

(e) A9

(f) A15

A8 A9

A11

A10

A12

6 A13

3 2

4 A14

1 0 10−1

A7 0

Time (s)

FLARE-T FE

−2

−3

Sensor

Acceleration, a (m s−2)

A15 100

Period, T (s)

0 10−1

100

Period, T (s)

100

Period, T (s)

Normalized absolute error (%)

Fig. 8. Held-out motion S014: (a)–(c) acceleration time histories; (d)–(f) 5%-damped pseudoacceleration response spectra; and (g) sensorwise error distributions.

finite-element model at the three representative locations. The agreement is strongest at A2, while residual amplitude differences remain at A9 and A15 during the most intense portion of the response. FLARE-T provides a closer representation of the measured spectral shape and reduces the pronounced short-period overprediction of the finiteelement model at the surface, although differences in individual spectral peaks remain. Figure 9(g) nevertheless shows consistently lower error distributions for FLARE-T at all sensors, demonstrating that the record-based correction remains effective for this stronger event. Across the three held-out motions, FLARE-T provides more consistent multi-depth acceleration histories than the finite-element model and generally improves the amplitudes and locations of the principal response-spectrum peaks. The reductions in the sensorwise error distributions are maintained over the range of shaking amplitudes represented by the test events, with the largest improvements occurring at locations where discrepancies in the original finite-element responses accumulate during upward propagation. The remaining spectral differences for the strongest motion indicate that the correction does not reproduce every local peak, but the overall results show that the information obtained from the calibration records remains predictive for events excluded from model development. The following application evaluates the same framework using field vertical-array observations and separately modeled horizontal components. 12

(a) A2

(b) A9 Measured FLARE-T FE

−2 −2

−1

−2

−3 0

Time (s)

(d) A2

(e) A9

(f) A15 12

17.5

15.0 12.5 10.0

15.0

12.5

10.0

7.5

5.0

2.5

0.0 −1 10

A7 0

Time (s)

17.5

100

Period, T (s)

FLARE-T FE

−4

Pseudo-spectral acceleration, Sa (m s−2)

(g) Error distribution by sensor

0.0 −1 10

A8 A9 A10 A11 A12

A13

A14

A15

0 10−1

100

Sensor

Acceleration, a (m s−2)

100

Period, T (s)

100

Normalized absolute error (%)

Period, T (s)

Fig. 9. Held-out motion S022: (a)–(c) acceleration time histories; (d)–(f) 5%-damped pseudoacceleration response spectra; and (g) sensorwise error distributions.

3.2 3.2.1

Field Vertical-Array Records at the Lotung Site Site Conditions and Numerical Model

The Lotung Large-Scale Seismic Test (LSST) site in northeastern Taiwan was established for field investigations of soil–structure interaction and was instrumented with free-field surface and downhole accelerometers. Because these instruments recorded the same earthquakes at several elevations within the soil deposit, the resulting vertical-array data have been widely used to identify site dynamic properties, infer depth-dependent stiffness and nonlinear soil behavior, and evaluate numerical site-response models (Elgamal et al., 1995; Zeghal et al., 1995; Chang et al., 1996; Glaser and Baise, 2000; Borja et al., 1999; Lee et al., 2006). The present application uses the free-field records to evaluate whether FLARE-T can predict depth-dependent motions for earthquake events excluded from record-based calibration. The modeled profile extends from the ground surface to a depth of 17 m and predominantly consists of layered silty sand and sandy silt deposits. The groundwater table is located 0.6 m below the ground surface. The horizontal acceleration recorded by DHB17 at a depth of 17 m is prescribed as the input motion, while the records from DHB11 at 11 m, DHB6 at 6 m, and FA1-5 at the ground surface define the response quantities considered in this study. These stations provide synchronized acceleration histories for multiple earthquake events in both the east–west and north–south directions, allowing the prediction performance to be examined for separate horizontal components under field conditions. The soil profile and the accelerometer locations used in the analysis are shown in Fig. 10. The numerical model was developed in Abaqus/Standard as a 17-m-deep and 0.085-m-wide plane-strain soil column. The profile was discretized using 200 four-node reduced-integration elements (CPE4R), giving 201 nodal elevations and 402 nodes. Nodes at the same elevation on the two lateral boundaries were constrained to undergo identical motion, the ground surface was traction free, and vertical displacement was restrained at the base. The prescribed horizontal acceleration was applied uniformly along the base. The initial effective-stress field was established under self-weight using K0 = 0.50, with the groundwater condition introduced at a depth of 0.6 m. Each dynamic analysis covered 15 s, with an initial time increment of 0.001 s, a maximum increment of 0.005 s, and an output interval of 0.005 s. Absolute accelerations were obtained at 200 elevations above the input boundary to form the dense numerical response field; responses corresponding to DHB11, DHB6, and FA1-5 were subsequently extracted to form the sparse simulated dataset. The soil layers were represented using the Dafalias–Manzari SANISAND model, which describes pressuredependent stiffness, state-dependent hardening and dilatancy, and the effect of fabric evolution under cyclic loading (Dafalias and Manzari, 2004). The elastic shear and bulk moduli are expressed as (2.97 − e)2 1+e 2(1 + ν) K= G, 3(1 − 2ν) G = G0 pa

p̄ pa

1/2 , (21)

17 m field profile

Water table 0.6 m

Layer / depth

Layer parameters

L01 0.00–1.21 m

Vs = 115 m/s γ = 16.2 kN/m³

L02 1.21–2.73 m

Vs = 126 m/s γ = 16.5 kN/m³

L03 2.73–4.24 m

Vs = 137 m/s γ = 20.4 kN/m³

L04 4.24–5.76 m

Vs = 148 m/s γ = 18.6 kN/m³

L05 5.76–7.27 m

Vs = 159 m/s γ = 17.8 kN/m³

L06 7.27–8.48 m

Vs = 167 m/s γ = 17.9 kN/m³

L07 8.48–10.90 m

Vs = 179 m/s γ = 19.3 kN/m³

L08 10.90–12.70 m

Vs = 193 m/s γ = 18.6 kN/m³

L09 12.70–14.20 m

Vs = 201 m/s γ = 19.7 kN/m³

L10 14.20–15.50 m

Vs = 210 m/s γ = 22.5 kN/m³

L11 15.50–17.00 m

Vs = 218 m/s γ = 19.3 kN/m³

Station depth FA1–5 0 m

DHB6 6 m

DHB11 11 m

DHB17 17 m

DHB17 horizontal input motion Unsaturated soil

Saturated soil

Fig. 10. Lotung soil profile and free-field accelerometers used for model input and response prediction. where e is the current void ratio, pa is atmospheric pressure, ν is Poisson’s ratio, and p̄ is the mean effective stress used in the constitutive integration. A layer-specific value of G0 was determined so that the initial shear modulus Gmax = ρVs2 reproduced the shear-wave velocity assigned to each layer in Fig. 10. The dependence of the constitutive response on the current density and confinement is introduced through the state parameter Ψ = e − ec ,

ec = e0 − λc

p̄ pa

ξ ,

(22)

where ec is the critical-state void ratio at the current mean effective stress. Plastic loading is governed by the yield condition r 2 f = ∥s − p̄α∥ − mp̄ = 0, (23) 3 where s is the deviatoric effective-stress tensor, α is the back-stress-ratio tensor, and m controls the opening of the yield surface. The evolution of α and the fabric tensor governs the hardening, dilatancy, and response following load reversal. Below the groundwater table, the excess pore-pressure response was evaluated within the constitutive implementation using a water bulk modulus of 2.2 GPa. The parameter set used for the numerical model is summarized in Table 3. 3.2.2 Dense-Response Model Training and Validation A dataset of 100 nonstationary horizontal input motions was used to generate the numerical responses for source-model training. The motions had peak ground accelerations ranging from 0.0075 to 0.112g and comprised 60 weak motions with PGA < 0.02g, 20 moderate motions with 0.02g ≤ PGA < 0.06g, and 20 strong motions with PGA ≥ 0.06g. Each motion had a duration of 15 s and was applied at the base of the finite-element model described in Section 3.2.1. The simulated motions represent a scalar horizontal component and were not assigned east–west or north–south labels; consequently, the same numerical dataset was used to establish the common source model for both horizontal components. Cases 001–080 were assigned to training and Cases 081–100 were reserved for validation, with each complete input motion and its corresponding response field retained within a single data partition. 14

Table 3. SANISAND parameters used for the Lotung finite-element model. Parameter pa (kPa) λc Mc m ν h0 nb nd cz Kw (GPa)

Value 100 0.019 1.25 0.0085 0.05 6.00 1.10 3.50 600 2.20

Parameter e0 ξ Me G0 eini ch A0 zmax pt /pa K0

Value 0.934 0.70 0.882 266.99–430.82 0.86 0.968 0.704 4.00 0.05 0.50

Table 4. Data and model settings for the Lotung application. Data source

Stage 1 Dense FE responses

Stage 2 Dense–sparse FE pairs

Training Validation Test Input Target Channels L r Components

Cases 001–080 Cases 081–100 – HL [ad ], u ad 100 3 14 Common to EW/NS

Cases 001–080 Cases 081–100 – HL [as,FE ] z 3 5 14 Common to EW/NS

Stage 3 Vertical-array records LSST05–06, LSST08–10, LSST12–13, LSST18 LSST11, LSST15–16 LSST07, LSST14, LSST17 HL [as ], u as 3 5 14 EW/NS separate

Of the 200 numerical output locations, 100 approximately uniformly spaced response channels were retained for Stage 1. The source encoder used a response-window length of L = 3, and the depth-distributed acceleration response was represented by r = 14 latent variables. The candidate library contained constant, linear, and quadratic functions of the latent state, while the base acceleration entered the latent dynamics as a separate linear forcing term. The synchronized input and response histories were supplied to the model at a time interval of 0.05 s. After allocation of the complete motions to the training and validation sets, each motion was divided into five contiguous sequences for optimization. Model selection was based on the error obtained by rolling the latent dynamics through the complete validation motions and reconstructing all 100 response channels. Prediction errors were evaluated using the NRMSE defined in Eq. (19). Table 4 summarizes the data allocation and principal model settings for all three stages of the Lotung application. Figures 11(a) and 11(b) compare the finite-element response field and the complete rollout prediction for validation Case 081. The prediction reproduces the timing and depth coherence of the principal acceleration cycles, including the concentration of stronger response during the main shaking interval and the subsequent decay in amplitude. The major positive and negative acceleration bands are recovered throughout the modeled depth, although localized differences remain in the amplitudes of several strong cycles. Figure 11(c) extends the evaluation to Cases 081–100. Each boxplot contains the NRMSE values obtained at the 100 response depths for one validation motion. The median depthwise errors are approximately 10–26%, while differences in the widths of the distributions show that the variation of prediction accuracy with depth depends on the input motion. These results confirm that the source model retains the principal temporal and depth-dependent characteristics of numerical site response for motions excluded from training. 3.2.3 Sparse-Response Transfer and Validation The Stage 1 encoder requires acceleration histories from 100 depth-distributed response channels, whereas the Lotung vertical array provides responses only at DHB11, DHB6, and FA1-5 within the modeled profile. Responses at these three locations were therefore extracted from each finite-element response field to form a sparse simulated dataset with the same observation configuration as the field records. This paired dataset was used to train a sparse encoder that estimates the latent coordinates established in Stage 1 from the three available response histories. The sparse encoder received acceleration windows from DHB11, DHB6, and FA1-5, using a window length of L = 5 and the previously established latent dimension of r = 14. Cases 001–080 were used for training, and Cases 081–100 were used for validation, consistent with the simulation partition adopted in Stage 1. Only the parameters of the sparse encoder were updated during this stage; the source encoder, decoder, latent dynamics, candidate library, 15

(b) FLARE prediction: case 081

0.2

20 10

−0.2

100

099

098

15.0

097

12.5

096

10.0

095

7.5

Time (s)

094

5.0

093

2.5

092

0.0

091

15.0

090

12.5

089

10.0

088

7.5

Time (s)

087

5.0

081

2.5

086

−0.1

085

0.0

084

8 10

0.1

083

Depthwise NRMSE (%)

Acceleration, a (m s−2)

Depth below surface (m)

0.0

082

(a) FE simulation: case 081

Validation case

Fig. 11. Stage 1 validation: (a) finite-element response field for Case 081, (b) complete rollout prediction, and (c) depthwise NRMSE for Cases 081–100. (b) DHB6_a1 (NRMSE = 20.3%)

0.10

0.3

0.1 0.0

0.0

0.00

(g) Validation NRMSE by output

0.2 0.1

0.05

0.2

DHB11_a1

FE target FLARE decoder

−0.1 −0.05

−0.1

−0.10

−0.2

0.0

2.5

5.0

7.5

10.0

12.5

15.0

0.0

(d) DHB11_a1: 5%-damped spectrum 0.6

5.0

7.5

10.0

12.5

15.0

−0.3 0.0

2.5

5.0

7.5

10.0

12.5

Time (s)

(e) DHB6_a1: 5%-damped spectrum

(f) FA1_5_a1: 5%-damped spectrum

1.25

0.5

15.0

1.5

1.00 0.4 0.75

0.3

0.50

0.2

0.5 0.25

0.1 0.0 −1 10

1.0 FA1_5_a1

Pseudo-spectral acceleration, Sa (m s−2)

Time (s)

2.5

DHB6_a1

−0.2

Decoder output

Acceleration, a (m s−2)

(a) DHB11_a1 (NRMSE = 19.1%)

100

Period, T (s)

0.00 −1 10

0.0 −1 10

100

Period, T (s)

100

Period, T (s)

NRMSE across validation cases (%)

Fig. 12. Stage 2 validation: (a)–(c) sparse-response reconstructions, (d)–(f) corresponding 5%-damped pseudoacceleration spectra, and (g) NRMSE distributions. response mapping, and normalization parameters remained fixed. Training employed the latent-state matching loss defined in Eq. (8), and the model with the lowest latent-state matching error over the validation cases was retained. The corresponding data allocation and model settings are summarized in Table 4. Figures 12(a)–12(c) compare the finite-element responses and their reconstructions for representative validation Case 081 at DHB11, DHB6, and FA1-5. The reconstructed histories follow the phase, principal amplitudes, and decay of the target motions at all three locations, with NRMSE values of 19.1, 20.3, and 20.0%, respectively. The corresponding 5%-damped pseudo-acceleration spectra in Figs. 12(d)–12(f) reproduce the period of the dominant spectral peak and the overall spectral shape, although the peak spectral amplitudes are underestimated to varying degrees. Figure 12(g) summarizes the results for Cases 081–100; each boxplot represents the NRMSE distribution of one output sensor over the 20 validation motions. DHB11 exhibits a lower median error, while the error distributions for DHB6 and FA1-5 remain similar to each other and within the same overall range. The consistent reconstruction of motions at the three elevations supports the use of a single sparse encoder to place the available sensor responses in the latent coordinates established by the dense-response model. 3.2.4 Record-Based Calibration and Held-Out Prediction The field records were divided by earthquake event, with LSST05, LSST06, LSST08, LSST09, LSST10, LSST12, LSST13, and LSST18 used for calibration; LSST11, LSST15, and LSST16 used for model selection; and LSST07, LSST14, and LSST17 retained for testing. The EW and NS components of each event were assigned to the same subset. Both directions used the source model and sparse encoder obtained from Stages 1 and 2, while the latent-dynamics coefficients were calibrated separately for each component. Only these coefficients were updated, and the test events were excluded from both coefficient adjustment and model selection. Figures 13(a) and 13(b) compare the Stage 1 and Stage 3 coefficient matrices for the EW component, while Figs. 13(c) and 13(d) provide the corresponding comparison for the NS component. The locations and signs of the 16

Latent-state equation

z1̇ z2̇ z3̇ z4̇ z5̇ z6̇ z7̇ z8̇ z9̇ ̇ z10 ̇ z11 ̇ z12 ̇ z13 ̇ z14

2.5

(b) EW - Stage 3

0.0

−2.5

SINDy coefficient, ξij

Latent-state equation

z1̇ z2̇ z3̇ z4̇ z5̇ z6̇ z7̇ z8̇ z9̇ ̇ z10 ̇ z11 ̇ z12 ̇ z13 ̇ z14

5.0

−5.0

2.5

(d) NS - Stage 3

0.0

−2.5

SINDy coefficient, ξij

Latent-state equation Latent-state equation

z1̇ z2̇ z3̇ z4̇ z5̇ z6̇ z7̇ z8̇ z9̇ ̇ z10 ̇ z11 ̇ z12 ̇ z13 ̇ z14

(a) EW - Stage 1

−5.0

1 z1 z2 z3 z4 z5 z6 z7 z8 z9 z10 z11 z12 z13 z14 2 z1 zz1 z1 z2 z1 z3 z1 z4 5 z1 z z1 z6 z1 z7 z1 8 z1 zz9 z1 z10 z1 z11 z1 z12 z1 z13 14 2 z2 zz2 z2 z3 z2 z4 z2 z5 z2 z6 z2 z7 z2 8 z2 zz9 z2 z10 z2 z11 z2 z12 z2 z13 14 2 z3 zz3 z3 z4 z3 z5 z3 z6 z3 z7 z3 8 z3 zz9 z3 z10 z3 z11 z3 z12 z3 z13 14 2 z4 zz4 z4 z5 z4 z6 7 z4 z z4 8 z4 zz9 z4 z10 z4 z11 z4 z12 z4 z13 14 2 z5 zz5 z5 z6 z5 z7 z5 z8 9 z5 z z5 z10 z5 z11 z5 z12 z5 z13 14 2 z6 zz6 z6 z7 z6 8 z6 zz9 z6 z10 z6 z11 z6 z12 z6 z13 14 2 z7 zz7 z7 8 z7 zz9 z7 z10 11 z7 z z7 z12 z7 z13 14 2 z8 z8 z8 zz9 z8 z10 z8 z11 z8 z12 z8 z13 14 2 z9 zz9 z9 z10 z9 z11 z9 z12 z9 z13 14 2 z10 zz10 z10 z11 z10 12 z10 zz13 14 2 z11 zz11 z11 z12 z11 z13 14 2 z12 zz12 z12 z13 14 2 z13 zz13 14 z142 u

z1̇ z2̇ z3̇ z4̇ z5̇ z6̇ z7̇ z8̇ z9̇ ̇ z10 ̇ z11 ̇ z12 ̇ z13 ̇ z14

SINDy library term

Fig. 13. Latent-dynamics coefficients for (a) EW at Stage 1, (b) EW at Stage 3, (c) NS at Stage 1, and (d) NS at Stage 3.

dominant coefficients remain largely unchanged after calibration, indicating that both direction-specific models retain the principal latent dynamics learned from the numerical simulations. The measured records primarily modify the magnitudes of a limited subset of coefficients, and the modified entries differ between EW and NS. These differences arise from the independent calibration of the two measured components within the same latent coordinates and do not represent direction-dependent soil material parameters. Each held-out event was initialized using the same five-sample response window from DHB11, DHB6, and FA1-5. The complete DHB17 acceleration history was then supplied as the forcing input to predict the subsequent responses at the three upper sensors. FLARE-T predictions were compared with the measured records and the corresponding finiteelement results before record-based correction. Timewise errors were evaluated using Eq. (20), and the 5%-damped pseudo-acceleration spectra were evaluated over periods from 0.10 to 5.00 s. In Figs. 14–16, panels (g) and (n) show the EW and NS error distributions, respectively; each boxplot contains the timewise error samples for one sensor and one event component. For LSST07, FLARE-T closely follows the phase and amplitudes of the principal EW and NS response cycles at DHB11 and DHB6 and substantially reduces the excessive late-time oscillations produced by the finite-element model. The improvement extends to FA1-5, where the timing of the principal surface pulses and the subsequent decay are reproduced more consistently in both directions. The predicted spectra also recover the dominant spectral bands at the three elevations more accurately than the finite-element results, particularly around the main intermediate-period peaks. The remaining surface discrepancies are concentrated near the largest acceleration pulses and the amplitudes of the dominant spectral peaks. Consistently lower error distributions in Figs. 14(g) and 14(n) show that the improvement is present in both horizontal components. LSST14 has smaller response amplitudes and a greater concentration of energy at short periods. FLARE-T reproduces the onset and decay of the measured oscillations in both directions, whereas the finite-element responses retain excessive high-frequency motion after the principal wave group. The short-period spectral peaks are also brought closer to the measured amplitudes at DHB11, DHB6, and FA1-5, although moderate underestimation remains near several peak ordinates. For every sensor, the FLARE-T error distributions lie below their finite-element counterparts in both Figs. 15(g) and 15(n), demonstrating that the improvement is not confined to one horizontal component. 17

(a) DHB11_a1 (EW)

(b) DHB6_a1 (EW) Measured FLARE-T FE

1.0

(g) EW error distribution

1.5

FLARE-T FE

1.0

DHB11_a1

0.5 0.0

0.0

−0.5

−0.5 −1

−1.0 2.5

5.0

7.5

10.0

12.5

15.0

0.0

2.5

5.0

7.5

10.0

Time (s)

(d) DHB11_a1 (EW)

(e) DHB6_a1 (EW)

12.5

15.0

0.0

2.5

5.0

7.5

10.0

12.5

15.0

Time (s) (f) FA1_5_a1 (EW)

Sensor

−1.0 0.0

DHB6_a1

a (m s−2)

0.5

0 10−1

100

Measured FLARE-T FE

100

T (s) (i) DHB6_a1 (NS)

Normalized absolute error (%)

T (s) (j) FA1_5_a1 (NS)

(n) NS error distribution

FLARE-T FE

1.0 1

0.5

a (m s−2)

0 10−1

100

T (s) (h) DHB11_a1 (NS) 1.0

FA1_5_a1

3 2

DHB11_a1

Sa (m s−2)

3 2

0.0 0

0.0

2.5

5.0

7.5

10.0

12.5

15.0

−1.0 0.0

Time (s)

5.0

7.5

10.0

12.5

−1 0.0

15.0

Time (s)

(k) DHB11_a1 (NS)

2.5

5.0

7.5

10.0

Time (s)

(l) DHB6_a1 (NS)

(m) FA1_5_a1 (NS)

12.5

15.0

Sensor

−0.5

DHB6_a1

−0.5

4 3

1.5

FA1_5_a1

Sa (m s−2)

2.0

1.0

1 1

0.5 0.0 −1 10

100

T (s)

0 10−1

100

T (s)

100

T (s)

100

Normalized absolute error (%)

Fig. 14. LSST07 predictions: (a)–(g) EW and (h)–(n) NS time histories, response spectra, and error distributions.

LSST17 contains multiple prominent wave groups distributed over a longer portion of the record. FLARE-T follows their phase evolution and depth-dependent amplitudes in both components, including the slower response cycles that dominate the latter part of the motion. The predicted spectra generally reproduce the locations and amplitudes of the principal peaks, while the finite-element results exhibit larger discrepancies in the short-period range and at the ground surface. The reductions in timewise error are comparable between EW and NS and occur at all three elevations, with slightly larger residual surface errors in the EW component. Local differences remain around individual peaks and intervals containing rapid changes in frequency content, but they do not obscure the principal response characteristics. Across the three held-out events, FLARE-T provides more consistent predictions of response phase, amplitude evolution, and the principal spectral peaks than the original finite-element model. The timewise error distributions are reduced at all three sensor elevations for both horizontal components, covering motions with different amplitudes, frequency contents, and durations. These results demonstrate that a common numerical source model can be transferred to sparse field observations and subsequently calibrated into direction-specific predictive models. The following subsection examines how the initial parameterization of the numerical source model influences this transfer process. 3.3

Influence of Prior Numerical-Model Calibration

A physics-based numerical source model is fundamental to FLARE-T because its simulated response fields provide the spatial information required to establish the latent response manifold and its input-driven dynamics. In engineering applications, however, the extent to which constitutive parameters can be calibrated before prediction depends on the available field observations. To examine the influence of this prior calibration, two numerical source models are considered for the Lotung site. The models share the same soil profile, shear-wave velocity distribution, SANISAND formulation, and numerical configuration, but use different values for selected constitutive parameters. This controlled comparison evaluates whether the record-based calibration stage can maintain consistent predictive performance for held-out events under different degrees of prior parameter calibration. The investigation therefore concerns the sensitivity of FLARE-T to specific parameter values within a common and physically representative numerical framework, for which the source model remains the basis of response learning and transfer. 18

0.15

0.2

0.10

0.1

0.05 0.0

FLARE-T FE

0.0

0.00 −0.1

−0.05 −0.1 2.5

5.0

7.5

10.0

12.5

15.0

0.0

2.5

5.0

7.5

10.0

12.5

15.0

0.0

2.5

5.0

7.5

10.0

Time (s)

(d) DHB11_a1 (EW)

(e) DHB6_a1 (EW)

(f) FA1_5_a1 (EW)

1.0

0.6

1.5

0.4

1.0

0.2

0.5

12.5

15.0

Sensor

−0.2

−0.10 0.0

(g) EW error distribution

DHB6_a1

0.1

a (m s−2)

(b) DHB6_a1 (EW) Measured FLARE-T FE

DHB11_a1

(a) DHB11_a1 (EW)

0.6 0.4

FA1_5_a1

Sa (m s−2)

0.8

0.2 0.0 −1 10

100

0.2

100

120

(n) NS error distribution FLARE-T FE

0.2 0.1 0.0 −0.2

−0.2 5.0

7.5

10.0

12.5

15.0

−0.3 0.0

Time (s) (k) DHB11_a1 (NS)

−0.4 2.5

5.0

7.5

10.0

12.5

15.0

0.0

2.5

5.0

7.5

10.0

Time (s)

(l) DHB6_a1 (NS)

(m) FA1_5_a1 (NS)

12.5

15.0

Sensor

2.5

DHB6_a1

a (m s−2)

Normalized absolute error (%)

0.2

−0.1

−0.2 0.0

T (s) (j) FA1_5_a1 (NS)

0.0

0.0 −0.1

100

T (s) (i) DHB6_a1 (NS) Measured FLARE-T FE

0.1

0.0 −1 10

100

T (s) (h) DHB11_a1 (NS)

DHB11_a1

0.0 −1 10

1.0 1.5

0.8

0.75

0.6

0.50

FA1_5_a1

Sa (m s−2)

1.00

1.0

0.4

0.25

0.2

0.00 −1 10

0.0 −1 10

100

T (s)

0.5 0.0 −1 10

100

T (s)

100

T (s)

Normalized absolute error (%)

Fig. 15. LSST14 predictions: (a)–(g) EW and (h)–(n) NS time histories, response spectra, and error distributions.

Table 5. SANISAND parameters used in the two Lotung source models. Parameter Baseline source model Precalibrated source model m 0.0100 0.0085 h0 7.05 6.00 eini 0.88 0.86

The baseline source model adopted the initial SANISAND parameter set, with m = 0.01, h0 = 7.05, and eini = 0.88, whereas the precalibrated source model used m = 0.0085, h0 = 6.00, and eini = 0.86, as summarized in Table 5. These parameters affect the onset and accumulation of plastic deformation, the rate of hardening, and the initial soil state; their modification therefore produces meaningful differences in the simulated cyclic response. Because the elastic stiffness in SANISAND also depends on the void ratio, the layer-specific G0 values were recalculated after changing eini so that the shear-wave velocity assigned to each layer remained equal to the profile presented in Fig. 10. This treatment allows the comparison to focus on differences in nonlinear constitutive response rather than differences in the initial wave-propagation velocity. The remaining SANISAND parameters and the numerical configuration were retained, and an independent simulated-response dataset was generated from each source model before applying the same three-stage FLARE-T procedure. Figure 17 compares the normalized training and validation MSE during third-stage calibration for the baseline and precalibrated source models, with the EW and NS results shown in Figs. 17(a) and 17(b), respectively. For both horizontal components, the training errors decrease to comparable levels, while the validation errors remain within the same overall range despite differences in their convergence paths. The changes in loss at epoch 200 coincide with the first scheduled sequential thresholding operation, during which coefficients below the prescribed threshold are removed and the optimizer is reinitialized before calibration continues; the subsequent changes occur for the same reason. The bounded loss histories following these operations show that both source-model parameterizations can be calibrated using the available field records, although the initial parameterization affects the detailed optimization path and the 19

(a) DHB11_a1 (EW)

0.0

−0.1

0.0

0.0 −0.2

−0.1

−0.2

−0.2 2.5

5.0

7.5

10.0

12.5

15.0

0.0

FLARE-T FE

0.2 DHB11_a1

0.1

0.0

(g) EW error distribution

−0.4 2.5

5.0

7.5

10.0

12.5

15.0

0.0

2.5

5.0

7.5

10.0

Time (s)

(d) DHB11_a1 (EW)

(e) DHB6_a1 (EW)

(f) FA1_5_a1 (EW)

12.5

15.0

Sensor

0.1

0.4

0.2

DHB6_a1

a (m s−2)

(b) DHB6_a1 (EW) Measured FLARE-T FE

0.2

1.25

0.75

0.6

1.0

0.50

0.4

0.5 0.25

0.2

0.00 −1 10

100

(n) NS error distribution FLARE-T FE

0.2

0.1

0.0

−0.1

−0.2

100

0.0

2.5

5.0

7.5

10.0

12.5

15.0

0.0

2.5

Time (s)

5.0

7.5

10.0

12.5

15.0

0.0

2.5

Time (s)

(k) DHB11_a1 (NS)

10.0

12.5

15.0

(m) FA1_5_a1 (NS) 1.25

0.8

0.6

7.5

Time (s)

(l) DHB6_a1 (NS)

1.0

0.8

5.0

1.00

0.6

0.75

0.4

0.2

0.50

0.2

0.0 −1 10

0.25

0.0 −1 10

100

Sensor

−0.2

DHB6_a1

a (m s−2)

Normalized absolute error (%)

T (s) (j) FA1_5_a1 (NS)

0.2

0.1

Sa (m s−2)

100

T (s) (i) DHB6_a1 (NS) Measured FLARE-T FE

0.2

0.0 −1 10

100

T (s) (h) DHB11_a1 (NS)

DHB11_a1

0.0 −1 10

0.0

FA1_5_a1

1.5

1.00

0.8

FA1_5_a1

Sa (m s−2)

1.0

0.00 −1 10

100

T (s)

100

T (s)

Normalized absolute error (%)

T (s)

Fig. 16. LSST17 predictions: (a)–(g) EW and (h)–(n) NS time histories, response spectra, and error distributions.

(a) EW 100

(b) NS 100

Normalized MSE

Baseline FLARE-T training Baseline FLARE-T validation Calibrated FLARE-T training Calibrated FLARE-T validation

10−1

200

400

600

800

1000

Epoch

200

400

600

800

1000

Epoch

Fig. 17. Third-stage optimization histories for the two source models: (a) EW and (b) NS.

validation error attained during training. For each source model and horizontal component, the model corresponding to the lowest validation error was retained for subsequent evaluation. Figure 18 compares the held-out predictions obtained from the baseline and precalibrated source models. For each model, the third-stage coefficients selected from the validation events were used to predict LSST07, LSST14, and LSST17, none of which contributed to coefficient calibration or model selection. The latent state was initialized from the same five-sample response window, after which the complete DHB17 motion was applied as the external input. The predicted and measured accelerations at DHB11, DHB6, and FA1-5 were pooled over the prediction interval, so that each point in Fig. 18 represents one measured–predicted acceleration pair at a particular sensor and time step. The 20

agreement was quantified using the coefficient of determination, N X

(b a i − ai )

R2 = 1 − i=1 N X

(24)

(ai − a)

i=1

where N is the total number of pooled samples and a is their mean measured acceleration. The predictions from both source models are concentrated around the 1:1 line, with R2 ranging from 0.814 to 0.912 across the three events and two horizontal components. The clearest effect of prior numerical-model calibration occurs for LSST07-EW, for which R2 increases from 0.858 to 0.907. Only small changes are observed for LSST07-NS and both components of LSST17, while the differences for LSST14 are also minor, with the baseline-source predictions giving slightly higher R2 . These results show that prior calibration can improve prediction for individual events, while third-stage calibration using field records produces broadly consistent held-out predictions from both numerical source models within the parameter range considered. 1.6

(a) LSST07-EW

(b) LSST14-EW

Predicted acceleration (m s−2)

Baseline FLARE-T: R = 0.858

0.10

Calibrated FLARE-T: R 2 = 0.907

0.8

Baseline FLARE-T: R = 0.834

Baseline FLARE-T: R 2 = 0.905

Calibrated FLARE-T: R 2 = 0.833

Calibrated FLARE-T: R 2 = 0.907

0.2

0.05

0.0

0.00

0.0

−0.05

−0.8

−0.2 Baseline FLARE-T Calibrated FLARE-T 1:1 line

−0.10 −1.6 −1.6

−0.8

0.0

0.8

1.6

(d) LSST07-NS

Baseline FLARE-T: R 2 = 0.911

Predicted acceleration (m s−2)

−0.10

0.12

Calibrated FLARE-T: R 2 = 0.912

−0.05

0.00

0.05

0.10

−0.2

(e) LSST14-NS

(f) LSST17-NS

Baseline FLARE-T: R 2 = 0.819

Baseline FLARE-T: R 2 = 0.899

Calibrated FLARE-T: R 2 = 0.814

Calibrated FLARE-T: R 2 = 0.905

0.00

0.0

−0.06

−1

0.2

0.06

0.0

−0.2

−0.12 −2

−2

−1

Measured acceleration (m s−2)

−0.12

−0.06

0.00

0.06

Measured acceleration (m s−2)

0.12

−0.2

0.0

0.2

Measured acceleration (m s−2)

Fig. 18. Measured and predicted accelerations for the held-out events: (a) LSST07-EW; (b) LSST14-EW; (c) LSST17EW; (d) LSST07-NS; (e) LSST14-NS; and (f) LSST17-NS. Taken together, the results of this chapter establish the complementary roles of numerical simulation and physical observations within FLARE-T. Across the centrifuge and Lotung applications, the dense numerical responses provided the depth-dependent response patterns and their input-driven evolution, the sparse-response transfer connected this information to the available sensor configuration, and the recorded motions corrected the transferred dynamics for prediction of held-out events. The resulting models reproduced multi-depth acceleration histories and the principal features of the corresponding 5%-damped pseudoacceleration spectra under both controlled experimental and field conditions. The comparison of the two Lotung source models further showed that the numerical model remains essential for providing a physically based response representation, whereas the final predictions depend less strongly on the precise prior calibration of the selected constitutive parameters. Better prior calibration can reduce the discrepancy presented to the transfer procedure and improve predictions for particular events, as observed for LSST07-EW, but 21

comparable performance for the remaining events and components indicates that the record-based correction can compensate for a substantial portion of the initial parameter-related error. FLARE-T therefore provides a practical means of retaining the physical information supplied by a site-response model while reducing the level of prior parameter calibration required to obtain reliable predictions for subsequent earthquake events.

Conclusion

This study proposed the Transfer-Enabled Forced Latent Autoencoder for Response Equations (FLARE-T) to improve site-response predictions by correcting discrepancies between numerical simulations and measured motions using limited site-specific records. The soil profile is represented as an input–output system, with base acceleration as its input and accelerations at multiple depths as its outputs. FLARE-T transfers the depth-dependent response representation learned from dense numerical simulations to the available sensor configuration and subsequently calibrates the latent dynamics using recorded motions. The method was evaluated using a layered-soil centrifuge test and the Lotung field vertical array, with predictive performance assessed on separate test sets in terms of multi-depth acceleration histories and 5%-damped pseudoacceleration response spectra. The test-set results showed that FLARE-T consistently provided more accurate site-response predictions than the original finite-element models in both applications. For the centrifuge test, the improvement was maintained across all instrumented depths and over the range of shaking intensities considered. For the Lotung site, lower prediction errors were obtained at all evaluated depths, for both horizontal components, and for motions with different amplitudes, frequency contents, and durations. In both cases, FLARE-T more accurately reproduced the response phase, amplitude evolution, decay, and principal spectral characteristics. The agreement across depths and test events further supports the ability of each calibrated model to reproduce coordinated responses at multiple locations through a shared lowdimensional latent state driven by the prescribed base motion. These results demonstrate that the response information transferred from numerical simulations can be effectively corrected using limited records to improve the prediction of test-set motions under both experimental and field conditions. The comparison between the two Lotung source models further clarifies the respective roles of numerical simulation and field observation in FLARE-T. Despite their different parameter settings, the two Lotung source models produced broadly comparable predictions for the test set. The generally similar performance indicates that FLARE-T reduces the dependence of the final prediction on the precise prior calibration of the selected constitutive parameters, which is particularly valuable when the available site data are insufficient for exhaustive numerical-model calibration. This result does not diminish the importance of the initial numerical model. Its essential role is to provide a physically reasonable coordinate system for the site response by describing the principal patterns of wave propagation and the coordinated variation of acceleration across depth. This role is directly connected to the low-dimensional responsemanifold assumption underlying FLARE-T: the soil profile is represented as a forced input–output system, in which base acceleration drives a low-dimensional latent state that describes the coordinated acceleration response across depth. From this perspective, the comparable predictions obtained from the two parameterizations suggest that both source models captured response patterns suitable for prediction, despite differences in their simulated evolution. Limited records could therefore correct much of the prediction discrepancy within the response representations supplied by these models. This distinction between representing the dominant response patterns and predicting their evolution explains how FLARE-T can reduce the demand for precise prior calibration while retaining the physical information provided by the numerical model. By addressing the sparse instrumentation and small event samples commonly encountered in engineering practice, the framework provides a practical route for converting existing numerical site models into observation-corrected predictive models for subsequent earthquake events.

Data and Code Availability Data and code generated or analysed during this study are available from the corresponding author upon reasonable request.

Declarations The authors have no conflicts of interest to declare that are relevant to the content of this article.

Acknowledgements This work is under the support of National Natural Science Foundation of China (Grant Numbers U2539204 and 52192675). 22

References Afacan, K. B., Brandenberg, S. J., and Stewart, J. P. (2014). “Centrifuge modeling studies of site response in soft clay over wide strain range.” Journal of Geotechnical and Geoenvironmental Engineering, 140(2), 04013003. Assimaki, D., Li, W., and Kalos, A. (2011). “A wavelet-based seismogram inversion algorithm for the in-situ characterization of nonlinear soil behavior.” Pure and Applied Geophysics, 168(10), 1669–1691. Bantis, J., Miranda, E., and Heresi, P. (2025). “Simplified site response analysis for regional seismic risk assessments.” Soil Dynamics and Earthquake Engineering, 188, 109022. Bazzurro, P. and Cornell, C. A. (2004). “Nonlinear soil-site effects in probabilistic seismic-hazard analysis.” Bulletin of the Seismological Society of America, 94(6), 2110–2123. Bergamo, P., Hammer, C., and Fäh, D. (2021). “On the relation between empirical amplification and proxies measured at swiss and japanese stations: Systematic regression analysis and neural network prediction of amplification.” Bulletin of the Seismological Society of America, 111(1), 101–120. Borcherdt, R. D. (1994). “Estimates of site-dependent response spectra for design: Methodology and justification.” Earthquake Spectra, 10(4), 617–653. Borja, R. I., Chao, H.-Y., Montáns, F. J., and Lin, C.-H. (1999). “Nonlinear ground response at lotung LSST site.” Journal of Geotechnical and Geoenvironmental Engineering, 125(3), 187–197. Chang, C.-Y., Mok, C. M., and Tang, H.-T. (1996). “Inference of dynamic shear modulus from lotung downhole data.” Journal of Geotechnical Engineering, 122(8), 657–665. Chen, A. T. F. (1985). “Transmitting boundaries and seismic response.” Journal of Geotechnical Engineering, 111(2), 174–180. Chen, S., Hu, X., Jiang, W., Wang, S., Chen, X., and Li, X. (2025). “Data-physical fusion deep learning for site seismic response using KiK-net records.” Earthquake Engineering & Structural Dynamics, 54(3), 993–1008. Dafalias, Y. F. and Manzari, M. T. (2004). “Simple plasticity sand model accounting for fabric change effects.” Journal of Engineering Mechanics, 130(6), 622–634. Dassault Systèmes Simulia Corp. (2020). Abaqus 2020 Analysis User’s Guide. Dassault Systèmes Simulia Corp., Johnston, RI. Elgamal, A.-W., Zeghal, M., Tang, H. T., and Stepp, J. C. (1995). “Lotung downhole array. i: Evaluation of site dynamic properties.” Journal of Geotechnical Engineering, 121(4), 350–362. Glaser, S. D. and Baise, L. G. (2000). “System identification estimation of soil properties at the lotung site.” Soil Dynamics and Earthquake Engineering, 19(7), 521–531. Griffiths, S. C., Cox, B. R., Rathje, E. M., and Teague, D. P. (2016). “Mapping dispersion misfit and uncertainty in Vs profiles to variability in site response estimates.” Journal of Geotechnical and Geoenvironmental Engineering, 142(11), 04016062. Hallal, M. M., Cox, B. R., and Vantassel, J. P. (2022). “Comparison of state-of-the-art approaches used to account for spatial variability in 1D ground response analyses.” Journal of Geotechnical and Geoenvironmental Engineering, 148(5), 04022019. Ilhan, O., Hashash, Y. M. A., Stewart, J. P., Rathje, E. M., Nikolaou, S., and Campbell, K. W. (2025). “Artificial neural network-based simulated site amplification models for central and eastern north america.” Earthquake Spectra, 41(4), 3190–3212. Kim, B., Hashash, Y. M. A., Stewart, J. P., Rathje, E. M., Harmon, J. A., Musgrove, M. I., Campbell, K. W., and Silva, W. J. (2016). “Relative differences between nonlinear and equivalent-linear 1-D site response analyses.” Earthquake Spectra, 32(3), 1845–1865. Kim, S., Hwang, Y., Seo, H., and Kim, B. (2020). “Ground motion amplification models for japan using machine learning techniques.” Soil Dynamics and Earthquake Engineering, 132, 106095. Kwok, A. O. L., Stewart, J. P., Hashash, Y. M. A., Matasovic, N., Pyke, R., Wang, Z., and Yang, Z. (2007). “Use of exact solutions of wave propagation problems to guide implementation of nonlinear seismic ground response analysis procedures.” Journal of Geotechnical and Geoenvironmental Engineering, 133(11), 1385–1398. Lee, C.-P., Tsai, Y.-B., and Wen, K.-L. (2006). “Analysis of nonlinear site response using the LSST downhole accelerometer array data.” Soil Dynamics and Earthquake Engineering, 26(5), 435–460. Lee, Y.-G., Kim, S.-J., Achmet, Z., Kwon, O.-S., Park, D., and Di Sarno, L. (2023). “Site amplification prediction model of shallow bedrock sites based on machine learning models.” Soil Dynamics and Earthquake Engineering, 166, 107772. 23

Li, L., Jin, F., Huang, D., He, C., and Ma, F. (2025). “Generalized deep neural network for seismic site response prediction with transfer learning.” Engineering Applications of Artificial Intelligence, 158, 111546. Li, L., Jin, F., Huang, D., and Wang, G. (2023). “Soil seismic response modeling of KiK-net downhole array sites with CNN and LSTM networks.” Engineering Applications of Artificial Intelligence, 121, 105990. Puzrin, A., Frydman, S., and Talesnick, M. (1997). “Effect of degradation on seismic response of israeli continental slope.” Journal of Geotechnical and Geoenvironmental Engineering, 123(2), 85–93. Rathje, E. M., Kottke, A. R., and Trent, W. L. (2010). “Influence of input motion and site property variabilities on seismic site response analysis.” Journal of Geotechnical and Geoenvironmental Engineering, 136(4), 607–619. Rodriguez-Marek, A., Bommer, J. J., Youngs, R. R., Crespo, M. J., Stafford, P. J., and Bahrampouri, M. (2021). “Capturing epistemic uncertainty in site response.” Earthquake Spectra, 37(2), 921–936. Roscoe, K. H. and Burland, J. B. (1968). “On the generalized stress–strain behaviour of wet clay.” Engineering Plasticity, J. Heyman and F. A. Leckie, eds., Cambridge University Press, Cambridge, UK, 535–609. Roten, D., Fäh, D., and Bonilla, L. F. (2014). “Quantification of cyclic mobility parameters in liquefiable soils from inversion of vertical array records.” Bulletin of the Seismological Society of America, 104(6), 3115–3138. Roten, D. and Olsen, K. B. (2021). “Estimation of site amplification from geotechnical array data using neural networks.” Bulletin of the Seismological Society of America, 111(4), 1784–1794. Seed, H. B., Ugas, C., and Lysmer, J. (1976). “Site-dependent spectra for earthquake-resistant design.” Bulletin of the Seismological Society of America, 66(1), 221–243. Stewart, J. P. and Afshari, K. (2021). “Epistemic uncertainty in site response as derived from one-dimensional ground response analyses.” Journal of Geotechnical and Geoenvironmental Engineering, 147(1), 04020146. Tao, Y. and Rathje, E. M. (2019). “Insights into modeling small-strain site response derived from downhole array data.” Journal of Geotechnical and Geoenvironmental Engineering, 145(7), 04019023. Tsai, C. C. and Hashash, Y. M. A. (2008). “A novel framework integrating downhole array data and site response analysis to extract dynamic soil behavior.” Soil Dynamics and Earthquake Engineering, 28(3), 181–197. Van Nguyen, D., Choo, Y., and Kim, D. (2024). “Deep learning application for nonlinear seismic ground response prediction based on centrifuge test and numerical analysis.” Soil Dynamics and Earthquake Engineering, 182, 108733. Yee, E., Stewart, J. P., and Tokimatsu, K. (2013). “Elastic and large-strain nonlinear seismic site response from analysis of vertical array recordings.” Journal of Geotechnical and Geoenvironmental Engineering, 139(10), 1789–1801. Zalachoris, G. and Rathje, E. M. (2015). “Evaluation of one-dimensional site response techniques using borehole arrays.” Journal of Geotechnical and Geoenvironmental Engineering, 141(12), 04015053. Zeghal, M., Elgamal, A.-W., Tang, H. T., and Stepp, J. C. (1995). “Lotung downhole array. ii: Evaluation of soil nonlinear properties.” Journal of Geotechnical Engineering, 121(4), 363–378. Zhang, H., Zheng, K., and Miao, Y. (2025). “Combining physical model with neural networks for earthquake site response prediction.” Soil Dynamics and Earthquake Engineering, 189, 109116. Zhu, C., Cotton, F., Kawase, H., Haendel, A., Pilz, M., and Nakano, K. (2022). “How well can we predict earthquake site response so far? site-specific approaches.” Earthquake Spectra, 38(2), 1047–1075. Zhu, C., Cotton, F., Kawase, H., and Nakano, K. (2023). “How well can we predict earthquake site response so far? machine learning versus physics-based modeling.” Earthquake Spectra, 39(1), 478–504. Zhu, C., Kwak, D.-Y., Pilz, M., Haendel, A., Xie, J., Kawase, H., and Cotton, F. (2026a). “OpenAmp: An open-source site amplification database to unlock the potential of machine learning.” Earthquake Spectra, 42(1), e70054. Zhu, Y., Chen, S., Li, X., and Du, X. (2026b). “Discovering latent response laws in forced physical systems.” arXiv preprint arXiv:2607.09801.

Record · ID 978444 · SHA-256 e175442499feab1d
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.