This paper presents a full-waveform inversion method for reconstructing the temporal and spatial distribution of unknown, incoherent dynamic traction in a heterogeneous, bounded solid domain from sparse, sur/uniFB01cial responses. This work considers SH wave motions in a two-dimensional (/two.pnumD) domain. The partial-di/uniFB00erential-equation (PDE)-constrained optimization framework is employed to search a set of control parameters, by which a mis/uniFB01t between measured responses at sensors on the top surface in- duced by targeted traction and their computed counterparts induced by estimated traction is minimized. To mitigate the solution multiplicity of the presented inverse problem, we employ the Tikhonov (TN) reg- ularization on the estimated traction function. We present the mathematical modeling and numerical implementation of both optimize-then-discretize (OTD) and discretize-then-optimize (DTO) approaches. The /uniFB01nite element method (FEM) is employed to obtain the numerical solutions of state and adjoint problems. Newton’s method is utilized for estimating an optimal step length in combination with the conjugate-gradient scheme, calculating a desired search direction, throughout a minimization process. Numerical results present that the complexity of a material pro/uniFB01le in a domain increases the error between reconstructed traction and its target. Second, the OTD and DTO approaches lead to the same inversion result. Third, when the sampling rate of the measurement is equal to the timestep for dis- cretizing estimated traction, the ratio of the size of measurement data to the number of the control parameters can be as small as /one.pnum:/one.pnum/two.pnum in the presented work. Fourth, it is acceptable to tackle the presented inverse modeling of dynamic traction without the TN regularization. Fifth, the inversion performance is more compromised when the noise of a larger level is added to the measurement data, and using the TN regularization does not improve the inversion performance when noise is added to the measure- ment. Sixth, our minimizer su/uniFB00ers from solution multiplicity less when it identi/uniFB01es dynamic traction of lower frequency content than that of higher frequency content. The wave responses in a computational domain, induced by targeted traction and its reconstructed one, are in excellent agreement with each other. Thus, if the presented dynamic-input inversion algorithm is extended in realistic /three.pnumD settings, it could reconstruct seismic input motions in a truncated domain and, then, replay the wave responses in a computational domain. /one.pnum Introduction There is a need for estimating incident seismic wave/uniFB01elds in a soil-structure system from limited seismic measurement data because, by using the identi/uniFB01ed seismic inputs, engineers can reconstruct (i.e., replay) responses within structures and soils during an earthquake event. There have been two dominant, conven- tional methods for the purpose mentioned above: the one is deconvolution and the other is the inversion of a seismic source pro/uniFB01le at a fault in a very large domain. The deconvolution algorithm has been used for the identi/uniFB01cation of an incoming seismic wave signal into a soil column by using vibrational measurement on the ground surface. For instance, Mejia and Dawson (/two.pnum/zero.pnum/zero.pnum/six.pnum) have presented the deconvolution to compute a seismic input signal by using the SHAKE program (Schnabel, /one.pnum/nine.pnum/seven.pnum/two.pnum), which solves the /one.pnumD seismic wave propagation problem in a domain of a semi-in/uniFB01nite ex- tent. Recently, there have been studies on the deconvolution of both vertical and horizontal components of sur/uniFB01cial measurement data to identify the vertical and horizontal input wave motions (Poul and Zerva, 1 /two.pnum/zero.pnum/one.pnum/eight.pnuma,b). We note that, although the deconvolution has been widely used in geotechnical earthquake engineering, it is e/uniFB00ective on individual soil columns only when soil strati/uniFB01cation is horizontally uniform, and incoming seismic waves vertically propagate. Namely, when the soil property is arbitrarily heteroge- neous (not horizontally layered), and incoming seismic waves, consisted of P , S, and/or surface waves, are highly incoherent (not vertically propagating), the deconvolution cannot e/uniFB00ectively reconstruct the incom- ing waves. On the other hand, there have been studies for inverting for seismic-source parameters (however simple or complicated an adopted seismic source model may be) at a hypocenter. For instance, Akcelik et al. (/two.pnum/zero.pnum/zero.pnum/two.pnum) presented a method to invert for a simpli/uniFB01ed seismic source time signal in a large /three.pnumD domain that includes a source at a fault. This method requires forward and inverse wave simulations of a very large domain. Upon characterizing the seismic source via inversion, then attention is typically turned on how to propagate the motion from the source to the surface, where the real interest is. However, in the large-scale source-inversion problem, the material properties of a large domain could be poorly characterized (on the other hand, the material properties could be better characterized by virtue of active wave source-based geotechnical characterization method (Fathi et al., /two.pnum/zero.pnum/one.pnum/six.pnum)). The limitations of the two methods, mentioned above, necessitate developing an alternative method that can identify arbitrary, incoherent incoming seismic waves in a truncated /two.pnumD or /three.pnumD domain by using sparse seismic measurement. Such a potential method could serve as an alternative to the deconvolution, while bypassing all the complexities associated with the inversion of the source at the hypocenter and the subse- quent propagation steps. Recently, Jeong and Seylabi (/two.pnum/zero.pnum/one.pnum/eight.pnum) presented prototype research that can recon- struct a seismic input signal propagating into a /one.pnumD truncated, heterogeneous, undamped solid system by using the partial di/uniFB00erential equation (PDE)-constrained optimization method. Lloyd and Jeong (/two.pnum/zero.pnum/one.pnum/eight.pnum) also show that the PDE-constrained optimization can reconstruct the discretized parameters of moving vibra- tional body forces in both space and time in a /one.pnumD heterogeneous, linear, elastic, undamped solid by using the sparse measurement of wave motions. These works were cast into a minimization problem, where a mis/uniFB01t between a measured response(s) at a sensor(s) induced by a targeted wave source pro/uniFB01le and a com- puted wave solution(s) induced by an estimated source pro/uniFB01le is minimized, and the PDE-constrained op- timization scheme analytically evaluates the gradient of a mis/uniFB01t with respect to control parameters, which parameterize an estimated dynamic input function. Because of such an analytical nature, its computa- tional e/uniFB03ciency of computing the gradient of a mis/uniFB01t with respect to control parameters does not depend on the number of them. Thus, it can update a large set of control parameters very e/uniFB03ciently. In addition, the PDE-constrained optimization can accommodate regularization to address the solution multiplicity of full-waveform inversion problems, typically caused by the sparsity of measurement, and to stabilize the convergence by penalizing an undesired aspect of an estimated pro/uniFB01le while enhancing a selected feature (e.g., smoothness) of targeted pro/uniFB01les. There had been a wide range of studies on elastodynamic inverse problems—e.g., full-waveform material inversion, inverse material design, full-waveform inverse-scattering, and source inversion—based on the PDE-constrained optimization as shown in the following literature review. Kang and Kallivokas (/two.pnum/zero.pnum/one.pnum/zero.pnuma) and Kang and Kallivokas (/two.pnum/zero.pnum/one.pnum/one.pnum) examined the numerical algorithms to image the distributions of the scalar-wave speeds in one-dimensional and two-dimensional solids that are surrounded by Perfectly-Matched-Layers (PML), where waves are forced to decay and they are prevented from re/uniFB02ecting o/uniFB00 the surrounding bound- aries (Kang and Kallivokas, /two.pnum/zero.pnum/one.pnum/zero.pnumb; Kucukcoban and Kallivokas, /two.pnum/zero.pnum/one.pnum/one.pnum, /two.pnum/zero.pnum/one.pnum/three.pnum; Fathi et al., /two.pnum/zero.pnum/one.pnum/five.pnumb; Poul and Zerva, /two.pnum/zero.pnum/one.pnum/eight.pnumc). Pakravan et al. (/two.pnum/zero.pnum/one.pnum/six.pnum) devised a new methodology to probe the elastic and attenuating parameters of two-dimensional viscoelastic layered solids. Kallivokas et al. (/two.pnum/zero.pnum/one.pnum/three.pnum), Fathi et al. (/two.pnum/zero.pnum/one.pnum/five.pnuma), Fathi et al. (/two.pnum/zero.pnum/one.pnum/six.pnum), and Kucukcoban et al. (/two.pnum/zero.pnum/one.pnum/nine.pnum) studied algorithms to invert for the Lamé parameters in two-dimensional and three-dimensional PML-truncated solid domains, and these computational studies made signi/uniFB01cant ad- vancement in geotechnical site characterization using dynamic tests. Tran and McVay (/two.pnum/zero.pnum/one.pnum/two.pnum) investigated the Gauss-Newton-based full waveform inversion approach to estimate the elastic modulus pro/uniFB01le in a two-dimensional domain. Mashayekh et al. (/two.pnum/zero.pnum/one.pnum/eight.pnum) investigated a new methodology to estimate the me- chanical properties of layered elastic or viscoelastic media by taking into account the dispersion relation of the layered medium. Tromp et al. (/two.pnum/zero.pnum/zero.pnum/eight.pnum) and Zhu et al. (/two.pnum/zero.pnum/one.pnum/seven.pnum) investigated the adjoint-tomography geo- physical inversion using the spectral element wave modeling in global- or regional-scale domains by using earthquake waves emitted from a seismic source of a known location and a known signal. Recently, Goh and Kallivokas (/two.pnum/zero.pnum/one.pnum/nine.pnum) investigated a new inverse metamaterial design method by minimizing the distance between the target and the computed group velocity pro/uniFB01les via a dispersion-constrained optimization method, and the method can be used for designing metamaterials, such as user-de/uniFB01ned omnidirectional band gaps in an elastic medium. In addition, it had been shown that strong discontinuities within solids, /two.pnum such as the boundaries of voids, can be identi/uniFB01ed by using inverse modelings. Namely, Guzina et al. (/two.pnum/zero.pnum/zero.pnum/three.pnum) and Jeong et al. (/two.pnum/zero.pnum/zero.pnum/nine.pnum) studied inverse scattering algorithms using the PDE-constrained optimization, asso- ciated with the boundary element method (BEM) wave solver, taking advantages of the moving boundary concept and the total derivative (Petryk and Mroz, /one.pnum/nine.pnum/eight.pnum/six.pnum). Nguyen-Tuan et al. (/two.pnum/zero.pnum/one.pnum/nine.pnum) also made recent progress in the inverse scattering algorithm, using the PDE-constrained optimization associated with the level-set based extended /uniFB01nite element method (XFEM) solver, so as to identify the geometry of voids in a static-hydro-mechanical system. Both BEM and XFEM wave solvers can model the boundaries of the strong discontinuities and update their geometries without cumbersome remeshing during an inversion process as opposed to a conventional /uniFB01nite element method (FEM) wave solver, which should remesh a domain to update the boundaries’ geometries (Jung et al., /two.pnum/zero.pnum/one.pnum/three.pnum). On the other hand, Aquino et al. (/two.pnum/zero.pnum/one.pnum/nine.pnum) devised a novel algorithm to detect debonded interfaces (i.e., interface cracks or incomplete weld bonds) in compos- ite solids by using steady-state vibrational tests and a density function that characterizes the bonding at the interfaces in composite solids. Besides, the following studies have investigated the methods to identify dynamic input functions. Hasanov and Baysal (/two.pnum/zero.pnum/one.pnum/four.pnum) studied an algorithm to detect the time-independent spatial load distributions of a dynamic source on a cantilever beam. Binder et al. (/two.pnum/zero.pnum/one.pnum/five.pnum) also recover virtual, stationary wave sources at possible locations of structural anomalies using the adjoint equation approach. Walsh et al. (/two.pnum/zero.pnum/one.pnum/three.pnum) reported the inverse problems for the identi/uniFB01cation of dynamic sources in acoustics and elastodynamics, employing a DTO approach. That is, the discrete form of the forward wave equation at each time step is imposed into a Lagrangian, and the adjoint equation and the gradient of the Lagrangian with respect to the source parameters are derived in discrete forms. In particular, Walsh et al. (/two.pnum/zero.pnum/one.pnum/three.pnum) suggested that the DTO approach is more suitable than the OTD counterpart when the nonlinearity is considered in the forward problem because the discrete, linearized forward equation of every time step can be individu- ally side-imposed into the Lagrangian. The PDE-constrained optimization has been also used for identifying optimal, non-moving wave source pro/uniFB01les that can focus wave energy to speci/uniFB01c areas in solids (Tadi et al., /one.pnum/nine.pnum/nine.pnum/six.pnum; Jeong et al., /two.pnum/zero.pnum/one.pnum/zero.pnum, /two.pnum/zero.pnum/one.pnum/five.pnum; Karve et al., /two.pnum/zero.pnum/one.pnum/five.pnum; Karve and Kallivokas, /two.pnum/zero.pnum/one.pnum/five.pnum; Jeong and Kallivokas, /two.pnum/zero.pnum/one.pnum/six.pnum; Karve et al., /two.pnum/zero.pnum/one.pnum/six.pnum). Due to the aforementioned robustness and scalability of the PDE-constrained optimization for inverse prob- lems, it is worth continuing to investigate it so as to identify incoherent seismic input motions in a multi- dimensional domain from sparse seismic measurement data. This research reconstructs the spatial and temporal distributions of incoherent dynamic traction on a boundary of a heterogeneous, bounded, un- damped solid system of anti-plane motions by using the PDE-constrained optimization method. Our numerical experiments show that the presented dynamic-input identi/uniFB01cation approach can success- fully identify unknown targeted traction without knowing any information about the target in heteroge- neous domains. Our parametric studies investigate the performance of the presented inverse modeling with respect to the complexity of material heterogeneity in a domain, the number of sensors, the regu- larization intensity factor, the optimization modeling type (i.e., OTD versus DTO), the noise level in mea- surement data, and the traction signal type (i.e., a high-frequency Ricker wavelet versus a low-frequency realistic seismic signal). /two.pnum Problem De/uniFB01nition This study is aimed at reconstructing the spatial and temporal distributions of dynamic traction on a bound- ary of an undamped solid by using measured wave responses at sparsely-distributed sensors on the top surface of the solid (see Fig. /one.pnum). The geometries and the material properties of the solid are assumed to be known in advance. /two.pnum./one.pnum The governing equation The governing equation for the SH wave propagation in the undamped solid domain is (we omit to show the spatial and/or temporal dependency of the variables): r(Gru)@2u @t2= 0;on J; (/one.pnum) whereu=u(x;y;t )denotes the displacement /uniFB01eld in the anti-plane ( z) direction of the wave motion of a solid particle (i.e., SH wave motion); x,y, andtdenote horizontal and vertical coordinates and time; G(x;y)and(x;y)denote the shear modulus and the mass density of the solid; denotes the domain, /three.pnum Sensors measuring vibrations SH waveF(x,t): dynamic shear traction in the anti-plane directionyx(0,0)Figure /one.pnum: Problem setting. and J= (0;T]is the time interval of interest. The solid is subject to a traction-free condition on the top surface ( t) and dynamic shear stress on the bottom surface ( b): G@u @y(x;0;t) = 0;0xL; (/two.pnum) G@u @y(x;D;t ) =F(x;t);0xL; (/three.pnum) whereDis they-coordinate of b, andF(x;t)denotes the dynamic shear stress applied on b. The solid is constrained by /uniFB01xed boundary conditions on the left ( l) and right ( r) boundaries: u(0;y;t) = 0; D<y< 0; (/four.pnum) u(L;y;t ) = 0; D<y< 0: (/five.pnum) whereLis thex-coordinate of r. The governing wave physics is also subject to zero initial-value condi- tions: u(x;y;0) = 0; (/six.pnum) @u @t(x;y;0) = 0: (/seven.pnum) We note that this work considers a /two.pnumD bounded domain as a prototype for the seismic-input inversion prob- lem in the multi-dimensional setting. Continuing this work, we will investigate the seismic-input inversion in a /two.pnumD//three.pnumD unbounded (truncated) domain. /two.pnum./two.pnum Parameterization of an estimated dynamic traction function We discretize an estimated dynamic traction function, F(x;t), over space and time as: F(x;t) =NxX k=1NtX j=1k(x) j(t)Fkj; (/eight.pnum) where k(x)denotes the k-th component of a vector of global basis functions used for the spatial dis- cretization of F(x;t); j(t)denotes the j-th component of a vector of global basis functions used for the /four.pnum temporal discretization of F(x;t);Fkjdenotes the discretized value of F(x;t)at each discrete location xk and timetj; andNxandNtdenote the numbers of discretization points over space and time, respectively. Although the shape functions to construct k(x)and j(t)can be of any low order, linear shape func- tions are used for both of them in the presented inverse modeling. The sizes of the temporal and spatial discretization are set to be, respectively, the time-step size ( t) of the forward time integration and the element size ( x) of a mesh on bfor a forward wave solver. The presented inverse modeling is aimed at reconstructing the set of control parameters Fkj, of which corresponding wave responses in the domain are consistent with the measurement on t. /three.pnum Inverse Modeling—the optimize-then-discretize (OTD) approach This section presents the OTD modeling for identifying the temporal and spatial distributions of unknown tractionF(x;t)based on measured wave responses on the top surface of the solid. First, this section presents the mathematical modeling of deriving the /uniFB01rst-order optimality conditions in a continuous form. Second, this section shows the discrete forms of the state and adjoint equations and the gradient of the objective functional with respect to the control parameters. /three.pnum./one.pnum The objective functional We cast the presented inverse problem into a minimization problem, where we seek the values of control parameters (i.e., Fkjin (/eight.pnum) for all kandj) that correspond to a minimum (either global or local one) of an objective functional: L=ZT 0NsX i=1(umiui)2dt+RTN; (/nine.pnum) whereumidenotes the displacement /uniFB01eld of the measured wave response at the i-th sensor induced by a targetedF(x;t);uidenotes the computed counterpart due to an estimated F(x;t), which is constructed by estimated control parameters; and Nsdenotes the number of sensors. In this computational study, umiis synthetically created by using a forward wave solver with a pseudo target of F(x;t). The /uniFB01rst term of (/nine.pnum) is a mis/uniFB01t between umiandui. We speculate that the mis/uniFB01t functional is quadratic and convex (please see the Appendix C) so that the /uniFB01rst-order optimality conditions of the Lagrangian functional (/one.pnum/one.pnum) would be su/uniFB03cient for obtaining the inversion solution. However, due to the sparsity of measurement data, we hypothesize that the considered inverse problem would su/uniFB00er from the solution multiplicity, and the regularization could be e/uniFB00ective for improving the convergence of the inversion solution to a targeted pro/uniFB01le. To test the hypothesis, we employ the second term of (/nine.pnum), RTN, which denotes the Tikhonov (TN) regularization term: RTN=R 2ZT 0Z b@F(x;t) @x2 +@F(x;t) @t2 ddt; (/one.pnum/zero.pnum) whereRis the regularization factor, which adjusts the amount of penalty on the derivative of F(x;t). By minimizing the regularization term RTNalong with the mis/uniFB01t, we attempt to minimize the discontinuity ofF(x;t)and smooth it while mitigating the solution multiplicity of the presented inverse problem. It is well known that, when the material inversion is performed, the TN regularization on a material pro/uniFB01le overly smooths the discontinuity at the interface of layered media. Thus, the TN regularization is suited for identifying a smooth material pro/uniFB01le while the total variation (TV) is used for enhancing the discontinuity of the material pro/uniFB01le. Meanwhile, for the seismic input inversion, a typical time signal of a seismic input motion should be a smooth function because the high-frequency content of a discontinuous signal cannot be retained along the propagation path from a seismic source to a near-surface domain due to attenuation. Thus, this work used the TN regularization under a hypothesis that using the TN regularization can improve the uniqueness of the inversion solution while smoothing the estimated traction function. This hypothesis is tested in Example /four.pnum shown in the later section ‘Numerical Experiments’. /five.pnum /three.pnum./two.pnum Lagrangian functional By imposing the governing equation (/one.pnum) and the Neumann boundary condition (/three.pnum) onto the side of the objective functional via Lagrange multipliers, a Lagrangian functional is built as: A=ZT 0NsX i=1(umiui)2dt+RTN +ZT 0Z  r(Gru)@2u @t2 d dt +ZT 0Z bF G@u @yF(x;t) ddt; (/one.pnum/one.pnum) where=(x;y;t )andF=F(x;t)are the Lagrange multipliers. Note that the boundary condition and initial conditions are implicitly imposed in (/one.pnum/one.pnum): they are not shown in (/one.pnum/one.pnum) but used for the derivation of the adjoint and control equations. The /uniFB01rst-order optimality conditions of the Lagrangian functional lead to state, adjoint, and control equations. The satisfaction of these equations leads to an optimal solution, corresponding to the minimal value of the objective functional. /three.pnum./three.pnum The /uniFB01rst-order optimality conditions The /uniFB01rst-order optimality conditions of the Lagrangian functional Arequire the vanishing variations of A with respect to the state variable u(x;y;t ), the Lagrange variables (x;y;t )andF(x;t), and the control parameter=Fkj. Such vanishing conditions lead to a triad of state, adjoint, and control equations (Lions, /one.pnum/nine.pnum/seven.pnum/one.pnum): ;FA= 0 : The /uniFB01rst condition (state problem) ; (/one.pnum/two.pnum) uA= 0 : The second condition (adjoint problem) ; (/one.pnum/three.pnum) A= 0 : The third condition (control problem) : (/one.pnum/four.pnum) For the /uniFB01rst condition, the variation of Awith respect to the Lagrange variables (x;y;t )andF(x;t) vanishes when the state problem—the original governing wave equation (/one.pnum) and its associated boundary and initial-value conditions—is satis/uniFB01ed. Our inverse modeling procedure automatically satis/uniFB01es it by nu- merically solving the state problem for estimated control parameters. As the second condition, the variation of Awith respect to the state variable u(x;y;t )should vanish. Such a vanishing variational condition leads to the following adjoint equation (see the derivation of the adjoint problem in Appendix): r(Gr)@2 @t2= 2(umu)NsX i=1(xxi;yyi);on [0;T); (/one.pnum/five.pnum) where (xxi;yyi)is the Dirac delta. In the adjoint PDE, the di/uniFB00erence between umanduserves as the time signal of a point wave source at each location of a sensor. It is noteworthy that the strong form of the adjoint PDE is derived from the weak-form like equation (/seven.pnum/seven.pnum). Although the strong form of the adjoint PDE is weakly satis/uniFB01ed, the FEM solution of the adjoint PDE fully satis/uniFB01es the weak-form like equation (/seven.pnum/seven.pnum) and, thus, satis/uniFB01es the second condition of the /uniFB01rst-order optimality condition. On the other hand, if the Lagrangian functional is built by imposing the discrete form of the state PDE, the aforementioned issue does not arise because the corresponding discrete form of the adjoint problem fully satis/uniFB01es the second condition of the /uniFB01rst-order optimality condition. To study the latter aspect, this paper also presents the discretize-then-optimize (DTO) counterpart in the later section ‘Inverse Modeling—the discretize–then- optimize (DTO) approach’, and our numerical experiments tests its inversion performance in Example /two.pnum in the later section ‘Numerical Experiments’. The adjoint PDE is also subject to the following boundary conditions: (0;y;t) =(L;y;t ) = 0; D<y< 0; @ @y(x;0;t) =@ @y(x;D;t ) = 0;0xL; (/one.pnum/six.pnum) /six.pnum and the following /uniFB01nal-value conditions: (x;y;T ) = 0; @ @t(x;y;T ) = 0: (/one.pnum/seven.pnum) The third condition states that the variation of Awith respect to a scalar-valued control parameter =Fkj should vanish. The vanishing variation condition leads to the following control equation (see the derivation of the control problem in Appendix): A=@A @ =ZT 0Z b@ @(F(x;t))ddt+@RTN @ =ZT 0Z bk(x) j(t)ddtRZT 0Z b@2F(x;t) @x2+@2F(x;t) @t2 k(x) j(t)ddt | {z } the closed form of@A @=@L @= 0:(/one.pnum/eight.pnum) Note thatA=@A @is the derivative of Awith respect to a control parameter =Fkj. Since the side- imposed terms inAvanish,@A @is equivalent to@L @, which constitutes a gradient vector rL, whereis a vector of all the control parameters. The control equation (/one.pnum/eight.pnum) implies that rLat any estimated values ofcan be evaluated in a semi-analytical manner by using its closed form once the solutions of state and adjoint problems are computed. /three.pnum./four.pnum Finite element solution of the state problem To /uniFB01ndu(t)2Ufor allv2V, we cast the weak form of the state problem as: Z rv(Gru)d +Z v@2u @t2d =Z bvF(x;t)d; (/one.pnum/nine.pnum) wherevdenotes a test function. The function spaces for a scalar-valued uandvare de/uniFB01ned as: U=fu:u2H1; u l;r= 0g V=fv:v2H1; v l;r= 0g (/two.pnum/zero.pnum) To resolve the weak form numerically, we use the standard /uniFB01nite-element approximation. We approximate the test and trial functions, respectively, as: v(x;y)'vT (x;y); u(x;y;t )' (x;y)Tu(t); (/two.pnum/one.pnum) where (x;y)denotes a vector of global basis functions constructed by shape functions of each /uniFB01nite element mesh in the domain, and u(t)denotes a vector of nodal solutions of the state problem. Then, (/one.pnum/nine.pnum) reduces to the following discrete form: Mu(t) +Ku(t) =F(t); (/two.pnum/two.pnum) where ()denotes the second-order derivative of its subtended variable with respect to t;Mdenotes a global mass matrix; Kdenotes a global sti/uniFB00ness matrix; F(t)denotes a global force vector. They are de/uniFB01ned as: K=Z G@ @x@ T @x+@ @y@ T @y d ; M=Z  Td ; F(t) =Z b F(x;t)d: (/two.pnum/three.pnum) /seven.pnum /three.pnum./five.pnum Finite element solution of the adjoint problem To /uniFB01nd(t)2for allv2V, the weak form of the adjoint equation (/one.pnum/five.pnum) is obtained as: Z rv(Gr)d +Z v@2 @t2d =NsX i=12v(xi;yi)(um(xi;yi;t)u(xi;yi;t)): (/two.pnum/four.pnum) The function space for a scalar-valued is de/uniFB01ned as:  =f:2H1;  l;r= 0g (/two.pnum/five.pnum) The test and trial functions are approximated as follows: v(x;y)'vT (x;y); (x;y;t )' (x;y)T(t); (/two.pnum/six.pnum) where(t)is a vector of the nodal adjoint solution. The weak form of the adjoint problem changes to the following time-dependent discrete form: M(t) +K(t) =Fadj(t); (/two.pnum/seven.pnum) where Fadj(t)is de/uniFB01ned as: Fadj(t) = 2NsX i=1( (xi;yi)u(xi;yi;t) (xi;yi)um(xi;yi;t)): (/two.pnum/eight.pnum) We note that the speci/uniFB01c forms of the matrices KandMin the discrete form of the adjoint problem in (/two.pnum/seven.pnum) are identical to those for the state problem in (/two.pnum/three.pnum). /three.pnum./six.pnum Time integration We solve the time-dependent discrete form of the state problem in (/two.pnum/two.pnum) by using the implicit Newmark time integration (i.e., average-acceleration scheme). We omit to show the detail of the forward time integration procedure of the state problem. On the other hand, for the time-dependent discrete form of the adjoint problem in (/two.pnum/seven.pnum), the time inte- gration begins from the /uniFB01nal time t=Tand ends at the initial time t= 0. That is, the backward time integration begins with the /uniFB01nal-value conditions, (T) = 0 ,_(T) = 0 , and (T) =M1Fadj(T): (/two.pnum/nine.pnum) For the backward time-marching procedure from t=T, the following approximations are used: _n=_n+1t 2nt 2n+1; (/three.pnum/zero.pnum) n=4 (t)2(n+1n)4 t_nn+1: (/three.pnum/one.pnum) Substituting (/three.pnum/zero.pnum) into (/three.pnum/one.pnum) yields: n=4 (t)2(n+1n) +4 t_n+1n+1: (/three.pnum/two.pnum) By inserting (/three.pnum/two.pnum) into the discrete form of the adjoint equation (/two.pnum/seven.pnum), we obtain the following equation that will be used for computing nat everyn-th time step: n= K+4 (t)2M1 Fadj(tn) +M4 (t)2n+14 t_n+1+n+1 : (/three.pnum/three.pnum) Once we obtain nfrom (/three.pnum/three.pnum), we obtain nfrom (/three.pnum/two.pnum) and, in turn, _nfrom (/three.pnum/zero.pnum). /eight.pnum /three.pnum./seven.pnum The discrete form of the gradient The gradient of the objective functional Lwith respect to a scalar variable can now be numerically com- puted as: r(=Fkj)L=xt(xk;D;tj)Rxt@2F(x;t) @x2+@2F(x;t) @t2 atxk;tj: (/three.pnum/four.pnum) Here, this works uses the aforementioned FEM solution of the adjoint problem for evaluating the gradient in (/three.pnum/four.pnum) under the OTD approach. /four.pnum Inverse Modeling—the discretize–then-optimize (DTO) approach This section presents the inverse modeling based on the DTO approach. The Lagrangian functional is built by imposing the discrete form of the state problem, using the discrete adjoint variable, into the objective functional in the discrete form. The /uniFB01rst-order optimality conditions are derived in the discrete form. The time-integration implementation of the state and adjoint problems are already embedded in their discrete forms. /four.pnum./one.pnum The discrete objective functional The discrete-form counterpart of the objective functional (/nine.pnum) is given by: ^L= (^um^u)TB(^um^u) +R 2^FTR^F; (/three.pnum/five.pnum) where ^u= [u0_u0u0u1_u1u1:::uN_uNuN]Tcorresponds to the space-time discretization of u(x;y;t )for (x;y)2 andt2[0;T], induced by an estimated F(x;t)(N is the number of time steps, and uiare the spatial degrees of freedom at the i-th time step); ^umis the space-time discretization of um(x;t)induced by a targeted F(x;t); and Bis a block diagonal matrix, determined as B= tBon the diagonal, where B is a square matrix that is zero everywhere except on the diagonals that correspond to a degree of freedom for which measured data are available; Ris the regularization factor; ^F= [0 0 F0F10 0:::FN0 0]Tis a global force vector corresponding to all the time steps—that is, the discrete control parameter Fkjare populated in ^F; and Ris the matrix corresponding to the discretization scheme used for the regularization terms de/uniFB01ned as: R=ZT 0Z b@w(x;0;t) @x@wT(x;0;t) @x+@w(x;0;t) @t@wT(x;0;t) @t ddt; (/three.pnum/six.pnum) wherew(x;y;t )denotes a vector of global basis functions, in both space and time, constructed by shape functions of each /uniFB01nite element mesh in the domain and the shape functions over the time. That is, an estimated traction function F(x;t)can be discretized as: F(x;t) =wT(x;y= 0;t)^F: (/three.pnum/seven.pnum) /four.pnum./two.pnum The discrete Lagrangian functional The Lagrangian functional corresponding to (/three.pnum/five.pnum) is built by imposing the discrete form of the state problem using the discrete adjoint variable: ^A= (^um^u)TB(^um^u) +R 2^FTR^F+^T(Q^u^F); (/three.pnum/eight.pnum) where ^= [0_001_11:::N_NN]Tis the discrete (space-time) Lagrange multiplier that en- forces the discrete forward problem as a constraint; and Qis the discrete forward operator de/uniFB01ned as: /nine.pnum Q=2 6666666666666664I 0 0 0 0 0 ::: 0 0 0 0 0 0 0 I 0 0 0 0 ::: 0 0 0 0 0 0 K 0M 0 0 0 ::: 0 0 0 0 0 0 L1L2L3Ke 0 0::: 0 0 0 0 0 0 a1I I 0a1I I 0::: 0 0 0 0 0 0 a0Ia2I Ia0I0I::: 0 0 0 0 0 0 ....................................... 0 0 0 0 0 0 ::: L1L2L3Ke 0 0 0 0 0 0 0 0 ::: a 1I I 0a1I I 0 0 0 0 0 0 0 ::: a 0Ia2I Ia0I0I3 7777777777777775; (/three.pnum/nine.pnum) where: a0=4 (t)2; a1=2 t; a2=4 t; (/four.pnum/zero.pnum) Ke =a0M+K; (/four.pnum/one.pnum) L1=a0M;L2=a2M;L3=M: (/four.pnum/two.pnum) /four.pnum./three.pnum The /uniFB01rst-order optimality condition in the DTO modeling The discrete optimality conditions of (/three.pnum/five.pnum) require that the variations of ^Awith respect to ^,^uand^Fvanish. The /uniFB01rst condition, taking the variation with respect to ^, recovers the discrete form of the state equation: @^A @^=Q^u^F= 0: (/four.pnum/three.pnum) For the second condition, the variation of ^Awith respect to ^ushould vanish: @^A @^u=QT^+ 2B(^um^u) = 0: (/four.pnum/four.pnum) Equation (/four.pnum/four.pnum) represents the discrete adjoint equation. Since it involves the transpose of Q, we solve it by marching backwards in time. For example, from the last two rows of (/three.pnum/nine.pnum), we obtain the /uniFB01nal condi- tions: N= 0; (/four.pnum/five.pnum) _N= 0; (/four.pnum/six.pnum) respectively; and the third row from the bottom yields: Ke/uniFB00TN= 2tB(uNumN) +a1_N+a0N; (/four.pnum/seven.pnum) which can be solved for N. For time steps n=N1; N2;:::; 1;we /uniFB01rst update nand_nas the following: n=MTn+1n+1; (/four.pnum/eight.pnum) _n=a2MTn+1_n+1a2n+1; (/four.pnum/nine.pnum) and, then, solve the following: Ke/uniFB00Tn=2tB(unumn) +a1_n+a0n +a0MTn+1a1_n+1a0n+1; (/five.pnum/zero.pnum) Finally, the /uniFB01rst three rows of (/one.pnum/five.pnum) result in the following equations. First, we solve: MT0=MT11; (/five.pnum/one.pnum) and, then, update _0and0as the following: _0=a2MT1_1a21; (/five.pnum/two.pnum) 0=KT0+a0MT1a1_1a01+ tB(u0um0): (/five.pnum/three.pnum) /one.pnum/zero.pnum We note that the backward adjoint time integration in the DTO approach di/uniFB00ers from that shown in its counterpart of the OTD approach. The third condition states that the variation of Awith respect to ^Fshould vanish: @^A @^F=RR^F^= 0; (/five.pnum/four.pnum) which represents the discrete control equation and implies that@^A @=@^L @is the component of the vector: RR^F^; (/five.pnum/five.pnum) at its row corresponding to =Fkj. /four.pnum./four.pnum Implementation of the regularization term in the gradient Figure /two.pnum:Fkjis surrounded by four elements in the space in terms of xandt: the horizontal and vertical axes represent, respectively, the xcoordinate and time t. This subsection shows the detail of evaluating the regularization term, RR^F, in (/five.pnum/five.pnum). The non-zero contri- bution of an element (shown in Fig. /two.pnum) in the space in terms of xandttoRmatrix is: Re= xt2 6664(1 31 (x)2+1 31 (t)2) (1 31 (x)2+1 61 (t)2) (1 61 (x)21 61 (t)2) (1 61 (x)21 31 (t)2) (1 31 (x)2+1 61 (t)2) (1 31 (x)2+1 31 (t)2) (1 61 (x)21 31 (t)2) (1 61 (x)21 61 (t)2 (1 61 (x)21 61 (t)2) (1 61 (x)21 31 (t)2) (1 31 (x)2+1 31 (t)2) (1 31 (x)2+1 61 (t)2) (1 61 (x)21 31 (t)2) (1 61 (x)21 61 (t)2) (1 31 (x)2+1 61 (t)2) (1 31 (x)2+1 31 (t)2)3 7775: (/five.pnum/six.pnum) Accordingly, the four elements surrounding Fkj(see Fig. /two.pnum) contribute to the 99submatrix (i.e., R/four.pnumE) of R. Then, the component of RR^Fcorresponding to Fkjcan be computed as: @RTN @Fkj=R(R/four.pnumE /five.pnum-th row )F/four.pnumE; (/five.pnum/seven.pnum) where the /five.pnum-th row of R/four.pnumE(i.e., R/four.pnumE /five.pnum-th row ) is: R/four.pnumE /five.pnum-th row = xt2 666666666666664(1 6(x)21 6(t)2) (1 3(x)22 3(t)2) (1 6(x)21 6(t)2) (2 3(x)2+1 3(t)2) (4 3(x)2+4 3(t)2) (2 3(x)2+1 3(t)2) (1 6(x)21 6(t)2) (1 3(x)22 3(t)2) (1 6(x)21 6(t)2):3 777777777777775T ; (/five.pnum/eight.pnum) /one.pnum/one.pnum andF/four.pnumEis de/uniFB01ned as: F/four.pnumE=2 6666666666664F(k1)(j1) F(k1)j F(k1)(j+1) Fk(j1) Fkj Fk(j+1) F(k+1)(j1) F(k+1)j F(k+1)(j+1)3 7777777777775: (/five.pnum/nine.pnum) Thus, under the DTO approach, (/five.pnum/seven.pnum) can be implemented as: @RTN @FkjDTO=Rxt 1 6F(k1)(j1)2F(k1)j+F(k1)(j+1) (x)2 +4 6Fk(j1)2Fkj+Fk(j+1) (x)2 +1 6F(k+1)(j1)2F(k+1)j+F(k+1)(j+1) (x)2 +1 6F(k1)(j1)2Fk(j1)+F(k+)(j1) (t)2 +4 6F(k1)j2Fkj+F(k+1)j (t)2 +1 6F(k1)(j+1)2Fk(j+1)+F(k+1)(j+1) (t)2 ;(/six.pnum/zero.pnum) while, its counterpart in the OTD approach can be implemented as: @RTN @FkjOTD=RxtFk(j1)2Fkj+Fk(j+1) (x)2 +F(k1)j2Fkj+F(k+1)j (t)2 ;(/six.pnum/one.pnum) which corresponds to (/three.pnum/four.pnum). /five.pnum Numerical Implementation of the Inversion Process By utilizing the semi-analytically evaluated gradient vector rL, this work iteratively updates a set of esti- mated control parameters by using the gradient-based minimization scheme as follows: (a) First, we compute synthetic measured data umat sensors by using pseudo-target traction F(x;t). (b) Then,u(x;t)is obtained by using estimated F(x;t)that is constructed by estimated control param- eters. (c) The adjoint problem is, then, solved by using the solution of the state problem in the previous step. (d) The gradient of the objective functional, rL, is evaluated. (e) Finally, the gradient-based minimization scheme updates the estimated control parameters via the conjugate-gradient method and the Newton’s method. The conjugate-gradient method determines the best search direction, and the Newton’s method determines an optimal step length. The numerical optimizer repeats the above steps (b) to (e) and iteratively solve for the control parameters that satisfy the vanishing control equation. A set of these steps is counted as an inversion iteration. The detailed procedure of the numerical optimizer is summarized in Algorithm /one.pnum. /five.pnum./one.pnum Conjugate gradient In every s-th inversion iteration, the gradient vector is computed as g(s)=rL(s). By using g(s), the search direction d(s)is computed by using the conjugate-gradient scheme (Fletcher and Reeves, /one.pnum/nine.pnum/five.pnum/one.pnum): d(s)=g(s) (s= 0and everym(e.g.,m= /five.pnum)); d(s)=g(s)+g(s):g(s) g(s1):g(s1)d(s1) (s1): (/six.pnum/two.pnum) Theg(s)is reset to be equal to g(s)at everyminversion iteration in order to eliminate the progressively- accumulated error in the search direction (Kang and Kallivokas, /two.pnum/zero.pnum/one.pnum/zero.pnuma). In the presented numerical exper- iments, we used m= /five.pnum. /one.pnum/two.pnum Algorithm /one.pnum Minimization Algorithm /one.pnum: Set TOL= 1010 /two.pnum: Build MandKmatrices for forward and adjoint wave solvers. /three.pnum: Compute umby using a forward wave solver with target traction F(x; t). /four.pnum: Set an iteration index s= 1 and initial control parameters ((s=1) = 0) and computeL(s=1). /five.pnum:while (L(s+1)>TOLL(1)ands < 103)do /six.pnum: Compute u(x; y; t )by solving the discrete form of the forward problem using (s). /seven.pnum: Solve the discrete form of the adjoint problem. /eight.pnum: Compute the components of a gradient vector, g(s)=rL(s). /nine.pnum: Compute an optimal search direction d(s)by using the conjugate-gradient scheme. /one.pnum/zero.pnum: Compute an optimal step length h(s)by using the Newton’s method. /one.pnum/one.pnum: Update estimated control parameters as (s+1) =(s)+d(s)h(s)and computeL(s+1). /one.pnum/two.pnum: s s+ 1 /one.pnum/three.pnum:end while /five.pnum./two.pnum Adaptively-calculated regularization factor The value of the regularization factor Rin the above gradient (/one.pnum/eight.pnum) determines the extent, to which the penalty is imposed on the oscillation of the spatial and temporal variation of F(x;t). IfRis too large, the estimated traction pro/uniFB01le may remain too smooth. If Ris too small, the numerical optimizer could su/uniFB00er from the solution multiplicity. To determine the optimal value of Rduring the inversion process, this work adopts the regularization factor continuation scheme (Kang and Kallivokas, /two.pnum/zero.pnum/one.pnum/zero.pnuma). To this end, we decomposeLintoLmandRTN=RLR, which are de/uniFB01ned as: Lm=ZT 0NsX i=1(umiui)2dt; LR=1 2ZT 0Z b@F(x;t) @x2 +@F(x;t) @t2 ddt: (/six.pnum/three.pnum) Accordingly,r(=Fkj)Lin (/one.pnum/eight.pnum) can be decomposed into the following two in the OTD approach: r(=Fkj)Lm=ZT 0Z bk(x) j(t)ddt; (/six.pnum/four.pnum) Rr(=Fkj)LR=RZT 0Z b@2F(x;t) @x2+@2F(x;t) @t2 k(x) j(t)ddt: (/six.pnum/five.pnum) or the following two in the DTO approach: r(=Fkj)^LmandRr(=Fkj)^LR, which are the components of the vectors, respectively, ^andRR^Fin (/five.pnum/five.pnum), at their rows corresponding to =Fkj. Here, Kang and Kallivokas (/two.pnum/zero.pnum/one.pnum/zero.pnuma) suggested imposing the following inequality: RjrLRj<jrLmj;orR<jrLmj jrLRj: (/six.pnum/six.pnum) Thus, in each inversion iteration, Rcan be set as: R=IRjrLmj jrLRj; (/six.pnum/seven.pnum) whereIRdenotes the regularization intensity factor. Since Kang and Kallivokas (/two.pnum/zero.pnum/one.pnum/zero.pnuma) heuristically found that the value of IRshould be 0<IR<0:5for the material inversion, we tested the performance of the presented inversion with respect to IRin the numerical experiments as well. /five.pnum./three.pnum Updating the estimated control parameters Starting with an initial guess for the control parameter vector , which is comprised of Fkj, the estimate ofcan be updated iteratively as: (s+1)=(s)+h(s)d(s); (/six.pnum/eight.pnum) whereh(s)is a scalar-valued step size and d(s)is a search-direction vector computed by the conjugate- gradient scheme, and sis the iteration index in the optimization process. For each iteration, s, the Newton’s /one.pnum/three.pnum method (Lloyd and Jeong, /two.pnum/zero.pnum/one.pnum/eight.pnum) is used to determine the optimal step size, h(s). Namely, in each s-th iteration in the numerical optimization process, there are up to /four.pnum sub-iterations ( ris the sub-iteration index). We begin from an initially estimated hs(r=1), and, for the next sub-iteration of (r+ 1) up torof /four.pnum, we updatehs(r+1)as the following: hs(r+1)=hsrL0(hsr)(s+1) L00(hsr)(s+1) =hsrL(hsr+)(s+1)L(hsr)(s+1) 2 L(hsr+)(s+1)2L(hsr)(s+1)+L(hsr)(s+1) 2 (/six.pnum/nine.pnum) whereL(hsr)(s+1)is the objective functional given the updated in (/six.pnum/eight.pnum) using hsr. For each sub-iteration, r, values for the /uniFB01rst and second derivatives of the objective functional with respect to hs—i.e.,L0(hsr)(s+1) andL00(hsr)(s+1)—are determined numerically via central di/uniFB00erence approximations. In the very right hand side term of (/six.pnum/nine.pnum), L(hsr)(s+1)isL(s+1)evaluated when (s+1)=s+ (hsr)ds. /six.pnum Numerical Experiments This section shows a set of numerical examples, investigating the performance of the presented inverse modeling with respect to various factors. In all the examples, we consider a square-shaped solid domain, of which extent is /six.pnum/zero.pnum m /six.pnum/zero.pnum m. To avoid an inverse crime, we compute the synthetic measured data um using an element size set as /zero.pnum./five.pnum m, while an element size of /one.pnum m is used in the FEM solvers for obtaining state and adjoint solutions in the computational domain. The same time step of /zero.pnum./zero.pnum/zero.pnum/one.pnum s is used in the forward and inversion problems. The solid has a uniform mass density of /one.pnum/zero.pnum/zero.pnum/zero.pnum kg/m3. Two targeted dynamic tractions, F1(x;t)andF2(x;t), are considered in this section. The targeted F1(x;t), of which amplitude changes over space and time, is shown in Fig. /six.pnum(a). Its peak amplitude is /two.pnum/zero.pnum/zero.pnum N/m2, and the total observation duration T1is /one.pnum./five.pnum s in the inversion simulation of F1(x;t). The time-dependent value of the targeted F1(x;t)at any speci/uniFB01c value of xis a Ricker wavelet (for instance, see F1(30;t)in Fig. /three.pnum(a)) with its central frequency of /one.pnum/zero.pnum Hz (see the frequency contents of F1(30;t)in Fig. /three.pnum(b)). The targeted F2(x;t)is shown in Fig. /one.pnum/nine.pnum(a), and its amplitude changes over space and time. Fig. /three.pnum(c) shows the time-dependent value of the targeted F2(30;t)and Fig. /three.pnum(d) shows its frequency contents. The signal is a part of a recorded ground motion signal during the /one.pnum/nine.pnum/nine.pnum/four.pnum Northridge earthquake from the Paci/uniFB01c Earth- quake Engineering Research Center (PEER) ground motion database (PEER, /two.pnum/zero.pnum/zero.pnum/zero.pnum). The total observation durationT2is /six.pnum s for the inversion simulation to identify the targeted F2(30;t). For the inversion process, we discretize estimated F1(x;t)by using /nine.pnum/one.pnum,/five.pnum/zero.pnum/zero.pnum control parameters—i.e., /six.pnum/one.pnum (over space)/one.pnum/five.pnum/zero.pnum/zero.pnum (over time)—and F2(x;t)by using /three.pnum/six.pnum/six.pnum,/zero.pnum/zero.pnum/zero.pnum control parameters—i.e., /six.pnum/one.pnum (over space) /six.pnum/zero.pnum/zero.pnum/zero.pnum (over time). The temporal and spatial intervals for the discretization are /zero.pnum./zero.pnum/zero.pnum/one.pnum s and /one.pnum m, respec- tively, for both forces. Initially-estimated values of all the control parameters are zero. Our numerical opti- mizer iteratively updates the values of all the control parameters by using the presented inverse modeling procedure. In what follows, six examples of numerical experiments are presented. The /uniFB01rst example is focused on the performance of inverting for F1(x;t)with respect to the complexity of the material pro/uniFB01le in the domain. To this end, we consider the following material pro/uniFB01les: •Material pro/uniFB01le /one.pnum: A homogeneous solid with the wave speed of Vs= /two.pnum/five.pnum/zero.pnum m/s, •Material pro/uniFB01le /two.pnum: A /two.pnum-layered solid with /one.pnum inclusion as in Fig. /four.pnum(a) with wave speeds of Vs1= /three.pnum/zero.pnum/zero.pnum m/s,Vs2= /three.pnum/five.pnum/zero.pnum m/s, and Vs3= /two.pnum/zero.pnum/zero.pnum m/s, •Material pro/uniFB01le /three.pnum: A /three.pnum-layered solid with /three.pnum inclusions as in Fig. /four.pnum(b) with wave speeds of Vs1= /four.pnum/zero.pnum/zero.pnum m/s,Vs2= /four.pnum/five.pnum/zero.pnum m/s,Vs3= /three.pnum/zero.pnum/zero.pnum m/s,Vs4= /three.pnum/five.pnum/zero.pnum m/s,Vs5= /two.pnum/zero.pnum/zero.pnum m/s, and Vs6= /two.pnum/five.pnum/zero.pnum m/s. The second example compares the performance of identifying F1(x;t)by using the OTD method versus the DTO method. The third example tests the inversion performance of identifying F1(x;t)with respect to the number of sensors on the top surface of the domain. The fourth example examines the convergence of the estimated F1(x;t)into the target with respect to the regularization intensity factor, IR. The /uniFB01fth example is focused on the performance of inverting for F1(x;t)with respect to the noise level. Lastly, the sixth example shows the capability of our inverse modeling to reconstruct F2(x;t), which is a realistic seismic signal as opposed to F1(x;t). Each example considers the result data of multiple cases, of which /one.pnum/four.pnum 00.511.5-1000100200 020400200040006000800010000 0246-2024 0510050010001500(a)(b) (c)(d)Figure /three.pnum: (a) The time signal of F1(x= 30;t); (b) the amplitude of Fourier Transform of F1(x= 30;t); (c) the time signal of F2(x= 30;t); and (d) the FFT of F2(x= 30;t). (a)(b) Figure /four.pnum: Heterogeneous solids: (a) Material pro/uniFB01le /two.pnum; and (b) Material pro/uniFB01le /three.pnum. /one.pnum/five.pnum input parameters are summarized in Table /one.pnum. For the sake of assessing the accuracy to reconstruct F(x;t) in the numerical results, the following error norm between estimated F(x;t)and its target is used: E=RT 0R jF(x;t)targetF(x;t)estimatej2d dt RT 0R jF(x;t)targetj2d dt100[%]: (/seven.pnum/zero.pnum) Table /one.pnum: Summary of all cases. Case number Material pro/uniFB01le Force pro/uniFB01le Approach # of sensors IR Noise level [%] /one.pnum /one.pnum /one.pnum OTD /three.pnum/zero.pnum /zero.pnum./five.pnum /zero.pnum /two.pnum /two.pnum /one.pnum OTD /three.pnum/zero.pnum /zero.pnum./five.pnum /zero.pnum /three.pnum /three.pnum /one.pnum OTD /three.pnum/zero.pnum /zero.pnum./five.pnum /zero.pnum /four.pnum /three.pnum /one.pnum DTO /three.pnum/zero.pnum /zero.pnum./five.pnum /zero.pnum /five.pnum /three.pnum /one.pnum OTD /one.pnum/five.pnum /zero.pnum./five.pnum /zero.pnum /six.pnum /three.pnum /one.pnum OTD /one.pnum/zero.pnum /zero.pnum./five.pnum /zero.pnum /seven.pnum /three.pnum /one.pnum OTD /five.pnum /zero.pnum./five.pnum /zero.pnum /eight.pnum /three.pnum /one.pnum OTD /one.pnum/five.pnum /one.pnum./zero.pnum /zero.pnum /nine.pnum /three.pnum /one.pnum OTD /one.pnum/five.pnum /zero.pnum./one.pnum /zero.pnum /one.pnum/zero.pnum /three.pnum /one.pnum OTD /one.pnum/five.pnum /zero.pnum./zero.pnum/one.pnum /zero.pnum /one.pnum/one.pnum /three.pnum /one.pnum OTD /one.pnum/five.pnum /zero.pnum./zero.pnum /zero.pnum /one.pnum/two.pnum /three.pnum /one.pnum OTD /one.pnum/five.pnum /zero.pnum./zero.pnum /one.pnum /one.pnum/three.pnum /three.pnum /one.pnum OTD /one.pnum/five.pnum /zero.pnum./zero.pnum /two.pnum /one.pnum/four.pnum /three.pnum /one.pnum OTD /one.pnum/five.pnum /zero.pnum./zero.pnum /three.pnum /one.pnum/five.pnum /three.pnum /one.pnum OTD /one.pnum/five.pnum /zero.pnum./zero.pnum/one.pnum /two.pnum /one.pnum/six.pnum /three.pnum /one.pnum OTD /one.pnum/five.pnum /zero.pnum./one.pnum /two.pnum /one.pnum/seven.pnum /three.pnum /one.pnum OTD /one.pnum/five.pnum /zero.pnum./five.pnum /two.pnum /one.pnum/eight.pnum /three.pnum /two.pnum OTD /one.pnum/five.pnum /zero.pnum./zero.pnum /zero.pnum Before our parametric studies on the inversion performance, we verify the theoretical derivation of the adjoint and control equations and the numerical implementation. That is, the gradients obtained by our semi-analytical OTD and DTO approaches are compared with that from the /uniFB01nite di/uniFB00erence (FD) approach. In this veri/uniFB01cation, we used the targeted F1(x;t), a heterogeneous domain of the material pro/uniFB01le /three.pnum, shown in Fig. /four.pnum(b), only one sensor in the center of t,IRof /zero.pnum, and the noise level of /zero.pnum%. In order to reduce the computational cost of the FD approach, the spatial distribution of F1(x;t)is set to be uniform, and, thus, only the temporal variation of F1(x;t)is considered. Fig. /five.pnum shows excellent agreement among the normalized gradients (i.e., rL=jrLj) that are calculated by using the OTD approach, the DTO approach, and the FD approach, respectively, at the /uniFB01rst inversion iteration. Thus, our theoretical derivation and numerical implementation of the presented inverse modeling are trustable. /six.pnum./one.pnum Example /one.pnum: Investigating the inversion performance with respect to the material pro/uniFB01le complexity. In this example, we test the performance of the presented inverse modeling with respect to the complexity of the material pro/uniFB01le by using the results of Cases /one.pnum to /three.pnum. We used the OTD approach, /three.pnum/zero.pnum sensors, IRof /zero.pnum./five.pnum, and the noise level of /zero.pnum% in these Cases. Fig. /six.pnum(b,c,d) show reconstructed dynamic traction for Cases /one.pnum, /two.pnum, and /three.pnum, respectively. Figs. /seven.pnum and /eight.pnum show that the values ofLandE, in general, decrease as the number of iterations increases. Fig. /seven.pnum depicts that Lshows a sawtooth behavior over iterations. We suggest that it occurs because of the penalty that is imposed by the regularization term on the derivative of F1(x;t). Namely, as seen in the later Examples /four.pnum and /six.pnum, when IRis equal to zero, the sawtooth behavior of Ldoes not occur. Besides, Fig. /eight.pnum shows that the more complex the material pro/uniFB01le is, the higher terminal value of Eis obtained. Fig. /nine.pnum shows that the wave responses of umdue to the targeted F1(x;t)match those of udue to the reconstructed one at two sensors in Cases /one.pnum-/three.pnum. It implies that our numerical minimizer is very e/uniFB00ective in minimizing the mis/uniFB01t: there is a very small di/uniFB00erence, of the scale of 1013to1012, among the terminal values of Lin Fig. /seven.pnum for Cases /one.pnum to /three.pnum. Even though umatchumat the end of the inversion simulation, our optimizer results /one.pnum/six.pnum 050010001500-2-1.5-1-0.500.51Figure /five.pnum: Comparison between the gradients generated by the OTD approach, DTO approach and the FD approximation. (a)(b) (c)(d) Figure /six.pnum: (a) Target and (b-d) Reconstructed F1(x;t)for Cases /one.pnum-/three.pnum at the /one.pnum/zero.pnum/zero.pnum/zero.pnumth iteration. Horizontal and vertical axis in the contour plot represent, respectively, the numbering of the discretized points over space (x) and time (t) of the distribution of F1(x;t). /one.pnum/seven.pnum 0200400600800100010-1310-1210-1110-1010-910-810-7Figure /seven.pnum: Example /one.pnum - Objective function, L, versus iterations with respect to the material pro/uniFB01les. 020040060080010000102030405060708090100 Figure /eight.pnum: Example /one.pnum - Error, E, versus iterations with respect to the material pro/uniFB01les. /one.pnum/eight.pnum in a higher terminal value of Efor a more complex material pro/uniFB01le. To explain this aspect, we note that Lloyd and Jeong (/two.pnum/zero.pnum/one.pnum/eight.pnum) reported that a more heterogeneous background solid leads to a larger terminal value of an error between a targeted moving, body force-typed wave source function and its estimated counterpart in a /one.pnumD solid setting. It is because, the more heterogeneous the material property of a domain is, the waves reverberate more inside the domain. That is, because of the reverberation, as seen in Fig. /six.pnum, the reconstructed traction has a stronger noise-like behavior, leading to a higher /uniFB01nal value of Ewhen a more heterogeneous material pro/uniFB01le is used. Thus, in a more heterogeneous solid setting, the inversion solver is less likely to converge towards a targeted traction function. Fig. /one.pnum/zero.pnum shows the snapshots of the wave responses in the entire computational domain, induced by (a) the targetedF1(x;t)and (b) its reconstructed counterpart in Case /three.pnum. Both responses are in excellent agreement with each other even though the terminal value of Ein Case /three.pnum is the largest among the Cases /one.pnum-/three.pnum. Thus, the presented dynamic-input inversion algorithm could be further developed for reconstructing seismic input motions and, then, “replaying" the corresponding wave responses in a truncated computational domain during seismic events. 00.511.5-20210-500.511.5-101210-5 00.511.5-20210-500.511.5-20210-5 00.511.5-20210-500.511.5-20210-5(a)(b)(c)(d)(e)(f) Figure /nine.pnum:umandu, at sensors placed on the top surface. (a-b) measured at x= /two.pnum/zero.pnum m and x= /four.pnum/zero.pnum m, respectively, for Case /one.pnum. (c-d) measured at x= /two.pnum/zero.pnum m andx= /four.pnum/zero.pnum m, respectively, for Case /two.pnum. (e-f) measured atx= /one.pnum/zero.pnum m andx= /two.pnum/zero.pnum m, respectively, for Case /three.pnum. /six.pnum./two.pnum Example /two.pnum: Comparison between the inversion performances of the OTD and DTO approaches This example compares the performance of the OTD approach with that of the DTO counterpart by using the results of Cases /three.pnum and /four.pnum, in both of which we reconstruct the targeted traction F1(x;t)by using the material pro/uniFB01le /three.pnum,IRof /zero.pnum./five.pnum, and the noise level of /zero.pnum%. Fig. /one.pnum/one.pnum shows the targeted F1(x;t)and its reconstructed counterparts obtained by using the OTD and DTO approaches, in Cases /three.pnum and /four.pnum, respectively. We could not visually noice the di/uniFB00erence between the two reconstructed ones in Fig. /one.pnum/one.pnum. Thus, we suggest that there is no di/uniFB00erence between the /uniFB01nal reconstructed traction function obtained by using the OTD approach and that by using the DTO approach in this /two.pnumD SH wave work. Fig. /one.pnum/two.pnum also shows that the terminal values of L andE, when using the OTD approach, are in excellent agreement with those using the DTO counterpart. However, a sawtooth behavior of Loccurs more when we use the DTO approach than the OTD approach. As mentioned earlier, the sawtooth behavior occurs because of the regularization. Namely, as shown in (/six.pnum/zero.pnum) and (/six.pnum/one.pnum), the gradient of the regularization term in the OTD approach is implemented di/uniFB00erently from the DTO counterpart so that their corresponding behaviors of Ldi/uniFB00er from each other. /one.pnum/nine.pnum 01020304050600102030405060 -4-3-2-10123410-6 (a) 01020304050600102030405060 -4-3-2-10123410-6 (b)Figure /one.pnum/zero.pnum: Wave responses, u(x;y;t ), at /one.pnum./five.pnum seconds due to (a) target traction and (b) reconstructed traction in Case /three.pnum. /two.pnum/zero.pnum 01020304050600102030405060 -4-3-2-10123410-601020304050600102030405060 -4-3-2-10123410-6 (a) (b)01020304050600102030405060 -4-3-2-10123410-6 (c)Figure /one.pnum/one.pnum: (a) Targeted F1(x;t)and (b,c) Reconstructed F1(x;t)for the OTD and DTO approaches, respec- tively, at the /one.pnum/zero.pnum/zero.pnum/zero.pnumth iteration. /six.pnum./three.pnum Example /three.pnum: Investigating the inversion performance with respect to the number of sensors In this example, we study the inversion performance with respect to the number of sensors on the top surface. We considered Cases /three.pnum, /five.pnum, /six.pnum, and /seven.pnum, which use /three.pnum/zero.pnum, /one.pnum/five.pnum, /one.pnum/zero.pnum, and /five.pnum sensors, respectively, and the material pro/uniFB01le /three.pnum, IRof /zero.pnum./five.pnum, and the noise level of /zero.pnum%. Fig. /one.pnum/three.pnum shows that, although increasing the number of sensors decreases the terminal value of Eof the inversion, using /five.pnum, /one.pnum/zero.pnum, /one.pnum/five.pnum, and /three.pnum/zero.pnum sensors gives rise to the terminal values of Eof the almost same magnitude, i.e., /two.pnum/five.pnum to /two.pnum/six.pnum% with a slight di/uniFB00erence of only up to /one.pnum% from each other. Thus, we suggest that when the sampling rate of the measurement is equal to the timestep for discretizing F1(x;t), the ratio of the size of measurement data to the number of the control parameters can be as small as /one.pnum:/one.pnum/two.pnum—e.g., (/five.pnum sensors /one.pnum/five.pnum/zero.pnum/zero.pnum timesteps):(/six.pnum/one.pnum nodes /one.pnum/five.pnum/zero.pnum/zero.pnum timesteps)—in the presented example. /six.pnum./four.pnum Example /four.pnum: Investigating the inversion performance with respect to IR This example studies the accuracy of the inversion with respect to the regularization intensity factor, IR. We used /one.pnum/five.pnum sensors on the top surface, the material pro/uniFB01le /three.pnum, and the noise level of /zero.pnum% for all the cases in this example. In addition to the Case /five.pnum in Example /three.pnum, where IRof /zero.pnum./five.pnum is used, we considered Cases /eight.pnum-/one.pnum/one.pnum, which useIRof /one.pnum./zero.pnum, /zero.pnum./one.pnum, /zero.pnum./zero.pnum/one.pnum, and /zero.pnum./zero.pnum, respectively. As can be seen in Fig. /one.pnum/four.pnum, the terminal value of Lfor a large value of IR(e.g.,IR= /one.pnum./zero.pnum) is larger than those for smaller values of IR. We also note that, when IRis /zero.pnum./zero.pnum, the sawtooth behavior of Ldoes not occur in Fig. /one.pnum/four.pnum. Furthermore, Fig. /one.pnum/five.pnum shows that the terminal value of EforIRof /one.pnum./zero.pnum is /five.pnum/four.pnum./nine.pnum % and is about twice higher than those (about /two.pnum/five.pnum %) for 0IR0:5. The inversion performance, in terms of Eover the iterations, is relatively unstable for IRof /one.pnum./zero.pnum, while it is stable for 0IR0:5. That is, using IRof /one.pnum./zero.pnum increasesEafter the /one.pnum/eight.pnum/zero.pnum-th iteration. Fig. /one.pnum/six.pnum compares the /uniFB01nal estimated traction function, when IRis /one.pnum./zero.pnum, with that for IRof /zero.pnum./zero.pnum. Fig. /one.pnum/six.pnum shows that the noise-like behavior is more severe in the reconstructed traction for IRof /one.pnum./zero.pnum than that for IRof /zero.pnum./zero.pnum. We also note that whenIRof /zero.pnum./zero.pnum was used,Eis the smallest (/two.pnum/four.pnum./six.pnum/five.pnum%) in this example as shown in Fig. /one.pnum/five.pnum. We discuss why using IRof /zero.pnum./zero.pnum leads to as a small terminal value of EasIRgreater than /zero.pnum./zero.pnum in the fol- /two.pnum/one.pnum 0200400600800100010-1310-1210-1110-1010-910-810-7(a) Mis/uniFB01t, L, versus iterations 020040060080010002030405060708090100 (b) Error, E, versus iterations Figure /one.pnum/two.pnum: Example /two.pnum - LandEwith respect to the inverse approach. /two.pnum/two.pnum 020040060080010002030405060708090100(a) Overall scale 02004006008001000242526272829303132 (b) Close-up scale Figure /one.pnum/three.pnum: Example /three.pnum - Error, E, versus iterations with respect to the number of sensors. /two.pnum/three.pnum lowing. As discussed in the previous /one.pnumD work of the seismic-input inversion (Jeong and Seylabi, /two.pnum/zero.pnum/one.pnum/eight.pnum), the proposed dynamic-traction inversion is naturally equipped with the smoothing e/uniFB00ect even without using the TN regularization. Namely, the FEM solver naturally /uniFB01lters out high frequencies of the estimated trac- tion function and smooths its temporal variation even without the regularization method. Similarly, it /uniFB01lters out the high-wavelength content of the spatial variation of an estimated traction function. Because of such inherent low-pass /uniFB01ltering of the FEM solver, as shown in our presented numerical results, even when TN regularization is not used, the inversion solver smooths an estimated traction function. 0200400600800100010-1310-1210-1110-1010-910-8 Figure /one.pnum/four.pnum: Example /four.pnum - Objective functional, L, versus iterations with respect to IR. /six.pnum./five.pnum Example /five.pnum: Investigating the inversion performance with respect to the noise level In this example, we focus on examining the performance of inverting F1(x;t)with respect to the noise level of random noise that is added to umprior to the inversion. We used the OTD approach, the material pro/uniFB01le /three.pnum, /one.pnum/five.pnum sensors, IRof /zero.pnum./zero.pnum, and examined the noise level of /zero.pnum%, /one.pnum%, /two.pnum%, and /three.pnum%, which correspond to Cases /one.pnum/one.pnum-/one.pnum/four.pnum, respectively. Fig. /one.pnum/seven.pnum shows that, the larger the noise level is, the larger the terminal values of LandEare obtained when no regularization is used. Fig. /one.pnum/eight.pnum shows the inversion performance with respect to IR(Cases /one.pnum/five.pnum-/one.pnum/seven.pnum in addition to the Case /one.pnum/three.pnum) when we use the material pro/uniFB01le /three.pnum, /one.pnum/five.pnum sensors, and the noise level of /two.pnum%. We note that IRof /zero.pnum./zero.pnum/one.pnum, /zero.pnum./one.pnum, and /zero.pnum,/five.pnum did not make any di/uniFB00erence in the terminal value of Ecompare to that for IRof /zero.pnum./zero.pnum. Thus, we suggest that using the TN regularization does not improve the inversion performance in the presented work when um contains noises. /six.pnum./six.pnum Example /six.pnum: Examining the feasibility of the presented inverse modeling to recon- struct a realistic seismic signal F2(x; t) In this example, we focus on examining the feasibility of inverting for a realistic seismic signal F2(x;t). We used the OTD approach, the material pro/uniFB01le /three.pnum, /one.pnum/five.pnum sensors, IRof /zero.pnum./zero.pnum, and the noise level of /zero.pnum%. Fig. /one.pnum/nine.pnum shows the excellent agreement between the targeted and reconstructed dynamic tractions, F2(x;t), in Case /one.pnum/eight.pnum. Fig. /two.pnum/zero.pnum and /two.pnum/one.pnum show the values of LandE, respectively, over iterations, and Ldecreases without the sawtooth behavior because IRof /zero.pnum./zero.pnum is used. Overall, F2(x;t)has much lower frequency content thanF1(x;t)so that the wave responses induced by F2(x;t)are less complex than those by F1(x;t) (see the wave responses in Fig. /two.pnum/two.pnum and compare them with those in Fig. /one.pnum/zero.pnum). Thus, our minimizer su/uniFB00ers from solution multiplicity less severely when it identi/uniFB01es F2(x;t)thanF1(x;t). Accordingly, the terminal /two.pnum/four.pnum 020040060080010000102030405060708090100(a) Overall scale 02004006008001000242526272829303132 (b) Close-up scale Figure /one.pnum/five.pnum: Example /four.pnum - Error, E, versus iterations with respect to IR. /two.pnum/five.pnum (a)(b)Figure /one.pnum/six.pnum: Example /four.pnum - Estimated traction F1(x;t)for (a) Case /eight.pnum, where IR= 1:0, and (b) Case /one.pnum/one.pnum, where IR= 0:0 value ofE, /three.pnum./eight.pnum/six.pnum%, for reconstructing F2(x;t)in Case /one.pnum/eight.pnum is much smaller than its counterpart, /two.pnum/nine.pnum./seven.pnum/six.pnum%, of reconstructing F1(x;t)in Case /one.pnum/one.pnum. Fig. /two.pnum/two.pnum shows the snapshots of the wave responses in the entire computational domain, induced by (a) the targeted F2(x;t)and (b) its reconstructed counterpart in Case /one.pnum/eight.pnum. In general, both responses are in excellent agreement with each other. /seven.pnum Conclusion We present the mathematical modeling and numerical implementation of a new inversion process for iden- tifying the spatial and temporal distributions of dynamic traction applied on a boundary of a /two.pnumD solid do- main with SH scalar wave motions. We tackle the inverse problem by using a gradient-based minimization scheme. The gradient of an objective functional is evaluated semi-analytically by using the adjoint solution. We present both OTD and DTO methods, each of which resolves the adjoint problem di/uniFB00erently from each other. Numerical results show the following /uniFB01ndings of the performance of this new inversion method. First, the complexity of the material pro/uniFB01le in a domain increases the error between the reconstructed traction and its target. The more heterogeneous the material property of a domain is, the waves reverberate more inside the domain. Because of the reverberation, the reconstructed traction has a stronger noise-like behavior, leading to a higher terminal value of Ewhen a more heterogeneous material pro/uniFB01le is used. Thus, in a more heterogeneous solid setting, the inversion solver is less likely to converge towards a targeted traction function. Second, the OTD and DTO methods lead to the same inversion results, but the distribution of L over the iterations shows a more signi/uniFB01cant sawtooth behavior when we use the DTO method than the OTD method. Third, when the sampling rate of the measurement is equal to the timestep for discretizing F(x;t), the ratio of the size of measurement data to the number of the control parameters can be as small as /one.pnum:/one.pnum/two.pnum in the presented work. Fourth, the regularization intensity IRshould not be too large: for instance,IRis recommended to be smaller than and equal to /zero.pnum./five.pnum. We also note that the terminal value ofE, whenIRof /zero.pnum is used, is as small as those when using /zero.pnum./zero.pnum/one.pnum IR/zero.pnum./five.pnum. Thus, it is acceptable to tackle the presented inverse modeling of dynamic tractions without the regularization. Fifth, the terminal values ofLandEincrease as the noise of a larger level is added to um, and using the TN regularization does not improve the inversion performance when noise is added to um. Sixth, our minimizer su/uniFB00ers from solution multiplicity less when it identi/uniFB01es dynamic traction of lower frequency content than that of higher frequency content. As shown in the numerical results, the wave responses in the entire computational domain, induced by the targeted traction and the reconstructed one, are in excellent agreement with each other in the presented highly-reverberating domain. Thus, if the presented dynamic-input inversion algorithm is extended into realistic /three.pnumD settings (see below for the details of the extension), it could allow engineers to reconstruct incident seismic motions and, then, to replay the wave responses in a /three.pnumD truncated domain. /two.pnum/six.pnum 0200400600800100010-1310-1210-1110-1010-9(a) Mis/uniFB01t, L, versus iterations 020040060080010002030405060708090100 (b) Error, E, versus iterations Figure /one.pnum/seven.pnum: Example /five.pnum - LandEwith respect to the noise level. /two.pnum/seven.pnum 0200400600800100011.21.41.61.822.22.42.62.8310-12(a) Mis/uniFB01t, L, versus iterations 020040060080010002030405060708090100 (b) Error, E, versus iterations Figure /one.pnum/eight.pnum: Example /five.pnum - LandEfor /two.pnum% of noise with respect of IR. /two.pnum/eight.pnum (a)(b)Figure /one.pnum/nine.pnum: Example /six.pnum - (a) Target and (b) Reconstructed F2(x;t)for Case /one.pnum/eight.pnum at the /six.pnum/zero.pnum/zero.pnum/zero.pnum-th iteration. 010002000300040005000600010-1710-1610-1510-1410-1310-12 Figure /two.pnum/zero.pnum: Example /six.pnum - Objective function, L, versus iterations. /two.pnum/nine.pnum 01000200030004000500060000102030405060708090100Figure /two.pnum/one.pnum: Example /six.pnum - Error, E, versus iterations. /seven.pnum./one.pnum Future extensions The presented domain does not fully represent a realistic one, which should be truncated by a wave- absorbing boundary condition (WABC) and subject to remote seismic excitation. Thus, we will extend this work as follows. A large extent of the computational domain will be truncated by using a WABC. Then, one can invert for seismic input motions in a truncated /two.pnumD//three.pnumD domain in the following two possible meth- ods. First, an estimated incident seismic wave can be modeled as an equivalent traction function on a WABC, e.g., dashpot (Lysmer and Kuhlemeyer, /one.pnum/nine.pnum/six.pnum/nine.pnum). The estimated traction function will be discretized over space and time, and all the discretized values will be control parameters to be identi/uniFB01ed. For instance, an x-directional traction function, Fx(x;y;t ), on a face that is perpendicular to the z-axis of a WABC will be discretized over space ( xandy) and time (t). Then, by using the presented traction-inversion approach, we can reconstruct the traction on the surfaces of a WABC. Second, the Domain Reduction Method (DRM) will be featured in a forward wave solver. Bielak and Chris- tiano (/one.pnum/nine.pnum/eight.pnum/four.pnum) and Bielak et al. (/two.pnum/zero.pnum/zero.pnum/three.pnum) had developed the DRM, by which free-/uniFB01eld wave motions are applied, as a dynamic input, along a /uniFB01ctitious boundary (also known as a DRM boundary) enclosed by the WABC. The DRM has been widely used for modeling wave behaviors of truncated solid domains subject to remote seismic excitations (Paolucci and Pitilakis, /two.pnum/zero.pnum/zero.pnum/seven.pnum; Tripe et al., /two.pnum/zero.pnum/one.pnum/three.pnum; Jeremić et al., /two.pnum/zero.pnum/one.pnum/three.pnum; Rahnema et al., /two.pnum/zero.pnum/one.pnum/six.pnum; Poursartip et al., /two.pnum/zero.pnum/one.pnum/seven.pnum; Zhang et al., /two.pnum/zero.pnum/one.pnum/nine.pnum). Thus, the extension of the presented method will be aiming at reconstructing free-/uniFB01eld seismic input motions at a DRM boundary. That is, we will spatially and temporally discretize estimated incident seismic wave/uniFB01eld functions at the two boundary surfaces of a single-element DRM bu/uniFB00er layer of the domain. For instance, we will discretize an x-component incident-wave/uniFB01eld func- tion at the horizontal boundary of a DRM bu/uniFB00er layer, i.e., u0 bx(x;y;t )oru0 ex(x;y;t )—the subscripts b andedenote the boundaries of the DRM layer neighboring the interior and exterior domains, respectively. Next, we will reconstruct the spatial and temporal distributions of the estimated incident seismic wave/uniFB01eld functions. Data Availability Statement Some or all data, models, or code generated or used during the study are available from the corresponding author by request. /three.pnum/zero.pnum 01020304050600102030405060 -1.5-1-0.500.511.510-7 (a) 01020304050600102030405060 -1.5-1-0.500.511.510-7 (b)Figure /two.pnum/two.pnum: Wave responses, u(x;y;t ), at /six.pnum./zero.pnum seconds due to (a) target load and (b) reconstructed load in Case /one.pnum/eight.pnum. /three.pnum/one.pnum •MATLAB code (.m format) of the presented inverse modeling that contains the optimization solver and the forward and adjoint wave solvers. •MATLAB datasets (.mat format) of the presented numerical results (Cases /one.pnum to /one.pnum/eight.pnum). Acknowledgment This material is based upon work supported by the National Science Foundation under Award No. CMMI- /one.pnum/eight.pnum/five.pnum/five.pnum/four.pnum/zero.pnum/six.pnum and CMMI-/two.pnum/zero.pnum/four.pnum/four.pnum/eight.pnum/eight.pnum/seven.pnum. Any opinions, /uniFB01ndings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily re/uniFB02ect the views of the National Science Foundation. /three.pnum/two.pnum A Derivation of the adjoint problem The variation ofAwith respect to ushould vanish as: uA=ZT 0uNsX i=1(umiui)2dt |{z} a+ZT 0Z r(Gru)d dt |{z} bZT 0Z @2u @t2d dt |{z} c +ZT 0Z bFG@u @yddt |{z} dZT 0Z bFF(x;t)ddt |{z} =0= 0: (/seven.pnum/one.pnum) Partain (/seven.pnum/one.pnum) can be written as: a=ZT 0uNsX i=1(umiui)2dt=ZT 0NsX i=12(umiui)uidt =ZT 0Z 2(umu)uNsX i=1(xxi;yyi)d dt; (/seven.pnum/two.pnum) where (xxi;yyi)is the Dirac delta function. Integrating part bin (/seven.pnum/one.pnum) by parts over space twice leads to: b=ZT 0Z r(Gru)d dt =ZT 0Z r(Gru)d dtZT 0Z (rGru)d dt =ZT 0Z r(Gru)d dt+ZT 0Z ur(Gr)d dtZT 0Z r(uGr)d dt: (/seven.pnum/three.pnum) Due to the Divergence Theorem, (/seven.pnum/three.pnum) becomes: b=ZT 0Z G@u @nddt+ZT 0Z ur(Gr)d dtZT 0Z uG@ @nddt =ZT 02 664Z tG@u @n|{z} =0d +Z bG@u @nd +Z l,rG@u @nd3 775dt +ZT 0Z ur(Gr)d dtZT 02 4Z t,buG@ @nd +Z l,ru|{z} =0G@ @nd3 5dt =ZT 0Z bG@u @yddt+ZT 0Z l,rG@u @nddt +ZT 0Z ur(Gr)d dtZT 0Z t,buG@ @nddt: (/seven.pnum/four.pnum) where@() @ndenotes a directional derivative of a variable in the direction of an outward unit normal vector non. Integrating cin (/seven.pnum/one.pnum) by parts over time twice leads to: c=ZT 0Z @2u @t2d dt =Z  @u @tT 0d +ZT 0Z @ @t@u @td dt =Z  @u @tT 0d +Z @ @tuT 0d ZT 0Z @2 @t2ud dt; (/seven.pnum/five.pnum) /three.pnum/three.pnum Because@u @t(t= 0) andu(t= 0) vanish, (/seven.pnum/five.pnum) becomes: c=Z  @u @t Td +Z  @ @tu Td ZT 0Z @2 @t2ud dt: (/seven.pnum/six.pnum) Due to (/seven.pnum/two.pnum), (/seven.pnum/four.pnum), and (/seven.pnum/six.pnum), (/seven.pnum/one.pnum) can be written as: uA=ZT 0Z b(F)G@u @yddt +ZT 0Z l,rG@u @nddtZT 0Z t,buG@ @nddt Z  @u @t Td +Z  @ @tu Td +ZT 0Z u" 2(umu)NsX i=1(xxi;yyi) +r(Gr)@2 @t2# d dt= 0: (/seven.pnum/seven.pnum) uA= /zero.pnum in (/seven.pnum/seven.pnum) is satis/uniFB01ed when we satisfy the adjoint problem in (/one.pnum/five.pnum) to (/one.pnum/seven.pnum). B Derivation of the control equation The variation ofAwith respect to should vanish as: A=@A @=ZT 0NsX i=1 2(umiui)@ui @ dt | {z } e+ZT 0Z @ @ r(Gru)@2u @t2 d dt | {z } =0 +ZT 0Z b@F @ G@u @yF(x;t) ddt | {z } =0+ZT 0Z r Gr@u @ d dt | {z } f ZT 0Z @2 @t2@u @ d dt |{z} g +ZT 0Z bFG@ @y@u @ ddtZT 0Z bF@ @(F(x;t))ddt +@RTN @= 0: (/seven.pnum/eight.pnum) Partein (/seven.pnum/eight.pnum) can be written as: e=ZT 0Z 2(umu)@u @NsX i=1(xxi;yyi)d dt: (/seven.pnum/nine.pnum) Integratingfin (/seven.pnum/eight.pnum) by parts twice over space leads to: f=ZT 0Z bG@ @y@u @ ddt+ZT 0Z l,rG@ @n@u @ ddt +ZT 0Z @u @ r(Gr)d dtZT 0Z t,b@u @ G@ @nddt: (/eight.pnum/zero.pnum) By integrating gin (/seven.pnum/eight.pnum) by parts over time twice and knowing that@ @t@u @(t= 0) and@u @(t= 0) vanish, we rewritegin (/seven.pnum/eight.pnum) as: g=Z  @ @t@u @ Td +Z  @ @t@u @ Td ZT 0Z @2 @t2@u @ d dt: (/eight.pnum/one.pnum) /three.pnum/four.pnum Due to (/seven.pnum/nine.pnum), (/eight.pnum/zero.pnum), and (/eight.pnum/one.pnum), (/seven.pnum/eight.pnum) becomes: A=ZT 0Z b(F)|{z} =0G@ @y@u @ ddt +ZT 0Z l,r|{z} =0G@ @n@u @ ddtZT 0Z t,b@u @ G@ @n|{z} =0ddt Z 2 4 |{z} =0@ @t@u @3 5 Td +Z 2 664@ @t|{z} =0@u @3 775 Td +ZT 0Z @u @" 2(umu)NsX i=1(xxi;yyi) +r(Gr)@2 @t2# | {z } =0d dt ZT 0Z bF@ @(F(x;t))ddt+@RTN @= 0: (/eight.pnum/two.pnum) Eq. (/eight.pnum/two.pnum) reduces to: A=ZT 0Z b@ @(F(x;t))ddt+@RTN @= 0: (/eight.pnum/three.pnum) The /uniFB01rst term of (/eight.pnum/three.pnum) is: ZT 0Z b@ @(F(x;t))ddt=ZT 0Z bk(x) j(t)ddt |{z } the /uniFB01rst term of (/one.pnum/eight.pnum); (/eight.pnum/four.pnum) In addition, the second term in (/eight.pnum/three.pnum) is: @RTN @=@ @R 2ZT 0Z b@F(x;t) @x2 +@F(x;t) @t2 ddt =RZT 0Z b@F @x@^F @x +@F @t@^F @t dddt = RZT 0@F @x^FdtL 0|{z} =0RZT 0Z b@2F @x2^Fddt + RZ b@F @t^FdT 0|{z} =0RZT 0Z b@2F @t2^Fddt =RZT 0Z b@2F(x;t) @x2+@2F(x;t) @t2 k(x) j(t)ddt | {z } the second term of (/one.pnum/eight.pnum); (/eight.pnum/five.pnum) where ^Fis de/uniFB01ned as: ^F=@F(x;t) @= k(x) j(t); (/eight.pnum/six.pnum) and we enforce that: @F @x= 0;atx= 0;L; @F @t= 0;att= 0;T: (/eight.pnum/seven.pnum) /three.pnum/five.pnum C On the quadratic and convex objective functional In this section, we prove that the objective functional is quadratic and convex. From the compact form of the state problem, ^ucan be de/uniFB01ned as: ^u=Q1^F: (/eight.pnum/eight.pnum) Due to (/eight.pnum/eight.pnum), the discrete objective functional ( ^L) can be written as: ^L= (^um^u)TB(^um^u) = (^umQ1^F)TB(^umQ1^F) = (^uT m^FT(Q1)T)B(^umQ1^F) =^uT mB^um|{z} a^FT(Q1)TB^um|{z} b^uT mB Q1^F|{z} c+^FT(Q1)TB Q1^F|{z} d; (/eight.pnum/nine.pnum) where the part ais a constant, and the parts bandcare linear functions so that their Hessians vanish. The partdin (/eight.pnum/nine.pnum) is a quadratic function. Therefore, proving that the part din (/eight.pnum/nine.pnum) is convex will show that ^L is convex. Sensor measuring vibrations F(x,t): dynamic shear traction in the anti-plane directionyx(0,0) Fixed boundaryFixed boundary Nodes Figure /two.pnum/three.pnum: Problem setting of Appendix C To this end, we consider the following simple example—a homogeneous square-shaped solid domain, of which extent is /six.pnum/zero.pnum m /six.pnum/zero.pnum m with its shear wave speed of /four.pnum/zero.pnum/zero.pnum m/s and mass density of /one.pnum/zero.pnum/zero.pnum/zero.pnum kg/m3. The shear stress is applied on the bottom surface, and one sensor is placed in the middle of top surface (see Fig. /two.pnum/three.pnum). The solid is constrained by /uniFB01xed boundary conditions on the left and right boundaries, and the element size is /three.pnum/zero.pnum m. Therefore, in this case, there are only three degrees of freedom in the discretized domain. We consider /two.pnum time steps, where t= 0:1s. In such a case, the matrices M,K, andQare the followings: M=2 4199980 99990 0 99990 399960 99990 0 99990 1999803 5;K=2 421334400053344000 0 53344000 426688000 53344000 053344000 2133440003 5; (/nine.pnum/zero.pnum) and Q=2 6666664I 0 0 0 0 0 0 I 0 0 0 0 K 0M 0 0 0 L1L2L3Ke 0 0 a1I I 0a1I I 0 a0Ia2I Ia0I0I3 7777775; (/nine.pnum/one.pnum) /three.pnum/six.pnum where: a0=4 (t)2; a1=2 t; a2=4 t; (/nine.pnum/two.pnum) Ke =a0M+K; (/nine.pnum/three.pnum) L1=a0M;L2=a2M;L3=M: (/nine.pnum/four.pnum) The vector ^Fis built as: ^F=2 66666640 0 F0 F1 0 03 7777775;F0=2 4F30;0 0 03 5;F1=2 4F30;1 0 03 5; (/nine.pnum/five.pnum) whereF30;0andF30;1are the components of a force vector at the node of x=/three.pnum/zero.pnum m and at the timesteps oft=/zero.pnum s andt=/zero.pnum./one.pnum s, respectively. In addition, the block diagonal matrix Bis de/uniFB01ned as: B= t2 6666664B10 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 B20 0 0 0 0 0 0 0 0 0 0 0 0 03 7777775;B1=B2=2 40 0 0 0 0 0 0 0 13 5: (/nine.pnum/six.pnum) Therefore, in this case, due to (/nine.pnum/zero.pnum)-(/nine.pnum/six.pnum), the part din (/eight.pnum/nine.pnum) can be written as: d=^FT(Q1)TB Q1^F = 1:251024F2 30;0+ 1:251024F2 30;1+ 2:51024F30;0F30;1: (/nine.pnum/seven.pnum) As shown in (/nine.pnum/seven.pnum), the part din (/eight.pnum/nine.pnum) is a quadratic function in terms of F30;0andF30;1. We note that the part din (/eight.pnum/nine.pnum) is convex if the eigenvalues of its Hessian matrix are non-negative (greater than or equal to /zero.pnum). Here, the Hessian matrix of the part din (/eight.pnum/nine.pnum) is: H=2 4@2d @F2 30;0@2d @F30;0@F30;1 @2d @F30;1@F30;0@2d @F2 30;13 5=2:510242:51024 2:510242:51024 ; (/nine.pnum/eight.pnum) of which eigenvalues are obtained as: |I*- H|= 2:510242:51024 2:510242:51024 = 0: (/nine.pnum/nine.pnum) Namely, the eigenvalues are: 1= 0;  2= 51024: (/one.pnum/zero.pnum/zero.pnum) Both eigenvalues are non-negative so that the part din (/eight.pnum/nine.pnum) is convex, and, therefore, the objective func- tional ^Lis also convex. /three.pnum/seven.pnum