This paper introduces a new inversion method for the reconstruction of complex, incoherent SH in- cident wave/uniFB01eld in a domain that is truncated by a wave-absorbing boundary condition (WABC), using a partial di/uniFB00erential equation (PDE)-constrained optimization method. In numerical examples, dynamic traction at the WABC mimics seismic incident wave/uniFB01eld. Estimated traction is discretized over space and time, and the discretized values are reconstructed by using seismic motion data that are sparsely made by sensors on the top surface of a domain and a vertical array. The discretize-then-optimize (DTO) approach is used in the mathematical modeling and numerical implementation, and the /uniFB01nite element method (FEM) is applied to solve state and adjoint problems. The numerical results show that incident, inclined plane waves, cannot be fully reconstructed if using only the top surface sensors. In order to im- prove the inversion performance, a vertical array of sensors on the side boundary of a domain should be included. Second, a su/uniFB03ciently large number of sensors must be employed to improve the algorithm’s in- version performance. Third, the minimizer su/uniFB00ers less from solution multiplicity when it identi/uniFB01es lower frequency traction (e.g., a realistic seismic signal). Fourth, the larger value of the inversion error in the reconstructed traction does not necessarily translate to an error of the same order of magnitude in the corresponding reconstructed wave responses in the computational domain. Fifth, our presented inver- sion algorithm’s accuracy is not compromised by the material complexity of a background domain. Lastly, the error in the reconstructed traction and the error in the corresponding wave responses grow when the noise of a larger level is added to the measurement data, but not in the same proportion. By extend- ing the presented method into realistic /three.pnumD settings, this algorithm can indicate where large amplitudes of stress waves (i.e., weak points) occur in built environments and soils in a domain of interest during seismic events. /one.pnum Introduction The ability to properly reconstruct response of infrastructure systems to actual earthquake loads, which are di/uniFB03cult to identify, would allow for improvement in the inspection of existing infrastructures that are hit by earthquakes as well as the design of new infrastructures. In order to properly reconstruct the re- sponse of infrastructure objects during earthquakes, the so-called earthquake soil-structure interaction (ESSI) systems, engineers need knowledge about earthquake wave/uniFB01elds that are exciting ESSI systems. In other words, engineers, analysts need to be able to identify incident seismic wave/uniFB01eld in a domain ad- jacent, within a few hundred meters, to the ESSI system. By using such reconstructed, realistic seismic wave/uniFB01elds, engineers can assess the impact of earthquakes on ESSI systems. It is important to note that the reconstructed incident seismic wave/uniFB01eld will contain a realistic combination of body and surface waves with broad frequencies and complex wave incoherence. To reconstruct seismic wave/uniFB01eld in a given ESSI domain, engineers could perform the following types of analysis. 1 The deconvolution algorithm has been used for the identi/uniFB01cation of an incoming seismic wave signal into a one-dimensional (/one.pnumD) soil column by using the transfer function of a soil column and sur/uniFB01cial measure- ment data (Ju, /two.pnum/zero.pnum/one.pnum/three.pnum; Mejia and Dawson, /two.pnum/zero.pnum/zero.pnum/six.pnum; Poul and Zerva, /two.pnum/zero.pnum/one.pnum/eight.pnuma,b). Namely, Mejia and Dawson (/two.pnum/zero.pnum/zero.pnum/six.pnum) showed the procedure to back-calculate an incident wave signal that induces seismic wave responses in a /one.pnumD soil column by using the deconvolution method and the SHAKE code (Schnabel, /one.pnum/nine.pnum/seven.pnum/two.pnum) as a /one.pnumD seismic wave solver. We note that, although the deconvolution has been widely used in geotechnical earthquake engineering, it is e/uniFB00ective on an individual soil column only when soil strati/uniFB01cation is horizontally uniform, and when /one.pnum component (/one.pnumC) incoming seismic waves propagate vertically. However, it is well known that seismic waves typically consist of (a) /three.pnum component (/three.pnumC) body waves, which consist of primary, compres- sional, P waves, and secondary, shear, S waves, and (b) /three.pnumC surface waves, which consist of Rayleigh waves, Love waves, and other types of irregular surface waves. Because incoming seismic waves, consisting of body and surface waves, are in general, highly incoherent (not vertically propagating), the deconvolution cannot reconstruct these /three.pnumC incoming waves. In addition, one could identify seismic-source pro/uniFB01les at ahypocenter in a regional-scale domain. For in- stance, 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 simula- tions in a very large domain (e.g., an area of /one.pnum/zero.pnum/zero.pnum km-by-/one.pnum/zero.pnum/zero.pnum km with a depth of /one.pnum/zero.pnum km). Upon characterizing the seismic source via inversion in a low-frequency scale (e.g., f <1Hz), seismologists could compute the low-frequency wave motions propagated from the source to a near-surface, where the real interest is. However, it is infeasible to accurately characterize the material properties of such a large domain. Thus, the reconstructed wave responses from this approach could be highly inaccurate. Moreover, the computa- tional cost of such a large-scale source-inversion problem is typically too large to identify the seismic source pro/uniFB01le with higher frequency components. In other words, the discretization spacing of /uniFB01nite element or /uniFB01nite di/uniFB00erence model in such a regional-scale source-inversion simulation cannot be su/uniFB03ciently small to make it e/uniFB03cient and accurate. Thus, in a practical sense, such methodology cannot reconstruct seismic wave/uniFB01eld with enough accuracy and with a frequency range of interest. As an alternative to the aforementioned conventional methods (i.e., the deconvolution in a /one.pnumD soil column and the source-characterization method in a large-scale domain), we present a new numerical method to identify arbitrary, incoherent incoming seismic waves in a multi-dimensional truncated domain by using sparse seismic measurement. To /uniFB01ll this algorithmic gap, we are investigating a new numerical method, where we could discretize the spatial and temporal distribution of an estimated incident wave and cal- culate the discretized parameters of the estimated wave that can minimize the discrepancy between the measured data at sensor locations and their computed counterparts induced by the estimated wave. Since the number of discretized parameters is quite large (e.g., millions), the PDE-constrained optimization-based inversion framework could be favorable in this problem because it allows semi-analytical calculation of a gradient of an objective functional with respect to a large set of control parameters. The PDE-constrained optimization framework has advanced the inverse problems associated with elasto- dynamic wave motions, and the associated applications span from material and scatterer characteriza- tion to dynamic-input identi/uniFB01cation. For examples, the geotechnical site characterization methodologies have been investigated in a soil domain that is truncated by Perfectly-Matched-Layers (PML) Fathi et al. (/two.pnum/zero.pnum/one.pnum/five.pnumb); 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); Mashayekh et al. (/two.pnum/zero.pnum/one.pnum/eight.pnum); Poul and Zerva (/two.pnum/zero.pnum/one.pnum/eight.pnumc) by using state-adjoint-control equation approaches Fathi et al. (/two.pnum/zero.pnum/one.pnum/five.pnuma, /two.pnum/zero.pnum/one.pnum/six.pnum); Kalli- vokas et al. (/two.pnum/zero.pnum/one.pnum/three.pnum); Kang and Kallivokas (/two.pnum/zero.pnum/one.pnum/zero.pnuma, /two.pnum/zero.pnum/one.pnum/one.pnum); Kucukcoban et al. (/two.pnum/zero.pnum/one.pnum/nine.pnum). The Gauss-Newton-based full-waveform inversion method, which is based on a gradient vector and a Hessian matrix, had also been studied for characterizing the material pro/uniFB01les in a truncated two-dimensional solid domain Pakravan et al. (/two.pnum/zero.pnum/one.pnum/six.pnum); Tran and McVay (/two.pnum/zero.pnum/one.pnum/two.pnum) for the application of the site characterization. Askan et al. (/two.pnum/zero.pnum/zero.pnum/seven.pnum, /two.pnum/zero.pnum/one.pnum/zero.pnum) and Askan and Bielak (/two.pnum/zero.pnum/zero.pnum/eight.pnum) described anelastic full-waveform seismic inversion algorithms for estimating seismic velocity and attenuation properties of large heterogeneous basins from surface observations of ground motions. Seidl et al. (/two.pnum/zero.pnum/two.pnum/zero.pnum, /two.pnum/zero.pnum/one.pnum/nine.pnum) used the PDE-constrained optimization to estimate both shear modulus and traction boundary condition distributions in a joint-inversion problem for the biomechanical imaging (BMI) application. Hamdaoui et al. (/two.pnum/zero.pnum/one.pnum/nine.pnum) utilized an adjoint method to identify the shear modu- lus and the loss factor of a viscoelastic sandwich structure. Lamert et al. (/two.pnum/zero.pnum/one.pnum/eight.pnum) investigated a method that utilizes elastodynamic full-waveform inversion to localize disturbances (e.g., concrete remnants) ahead of a tunnel. Reuber et al. (/two.pnum/zero.pnum/two.pnum/zero.pnum) presented an adjoint formulation that uses principal stress directions as data in the mis/uniFB01t function to identify rheological parameters (e.g., density, viscosity) of lithosphere and mantle. Wang and Brigham (/two.pnum/zero.pnum/one.pnum/seven.pnum) presented an inverse material characterization algorithm that combines direct inversion and optimization-based inversion to estimate the spatial distribution of elastic modulus in /two.pnum A targeted dynamic traction mimics a plane incident wave that impinges the domain with an incident angle of 45 degree. WABCWABCWABC Figure /one.pnum: The problem con/uniFB01guration, showing an exemplary domain with a homogeneous material of its wave speed Vs1and inclusions with their wave speeds Vs2andVs3. plate structures from incomplete response /uniFB01eld measurements. Kraft et al. (/two.pnum/zero.pnum/one.pnum/six.pnum) studied an adjoint state approach for the identi/uniFB01cation of vehicle parameters (e.g., sti/uniFB00ness of the suspension) from on-track mea- surements for calibration of railway vehicle models. Lloyd and Jeong (/two.pnum/zero.pnum/one.pnum/eight.pnum) has shown the adjoint-control equation approach to identify the pro/uniFB01les of moving wave sources of an arbitrary number in a heteroge- neous /one.pnumD solid. The discretize-then-optimize (DTO) approach—deriving the adjoint and control equations in discrete forms based on a discrete objective functional—has also advanced the /uniFB01eld of the identi/uniFB01cation of dynamic wave sources in solids Walsh et al. (/two.pnum/zero.pnum/one.pnum/three.pnum). An inverse-scattering problem—the /uniFB01eld of detect- ing complex geometries of wave-scattering objects Guzina et al. (/two.pnum/zero.pnum/zero.pnum/three.pnum); Jeong et al. (/two.pnum/zero.pnum/zero.pnum/nine.pnum); Nguyen-Tuan et al. (/two.pnum/zero.pnum/one.pnum/nine.pnum)—has also been studied by the PDE-constrained optimization framework. To investigate the feasibility to identify arbitrary, incoherent incoming seismic waves in a truncated do- main by using sparse seismic measurement by using the PDE-constrained optimization framework, Jeong and Seylabi (/two.pnum/zero.pnum/one.pnum/eight.pnum) studied a full-waveform source-inversion method to identify an incoming seismic wave in a /one.pnumD semi-in/uniFB01nite solid. Guidio and Jeong (/two.pnum/zero.pnum/two.pnum/zero.pnum) also presented the mathematical modeling and numer- ical implementation of a new inverse-source method to identify the spatial and temporal distributions of targeted traction on the boundary of a /two.pnumD heterogeneous, bounded, undamped solid by using the sparse measurement of anti-plane wave motions without any prior information about targeted traction. The work showed that the DTO-based inversion leads to the same performance as the optimize-then-discretize (OTD) counterpart, but the former is more straightforward in terms of the mathematical modeling than the lat- ter. The work also shows that using the Tikhonov (TN) regularization does not improve the accuracy of the inversion. Continuing the aforementioned papers Guidio and Jeong (/two.pnum/zero.pnum/two.pnum/zero.pnum); Jeong and Seylabi (/two.pnum/zero.pnum/one.pnum/eight.pnum), the presented paper discusses a numerical method that can fully reconstruct the comprehensive pro/uniFB01les (i.e., both spa- tial and temporal distributions) of complex, incoherent seismic incident wave/uniFB01elds propagating into a /two.pnumD truncated domain of anti-plane shear wave motions. /two.pnum Problem De/uniFB01nition This study explores a new full-waveform inversion method for reconstructing the spatial and temporal dis- tributions of dynamic traction in a /two.pnumD reduced domain from data of wave responses measured at sparsely- located sensors. Here, the semi-extent of the domain is truncated by using a wave-absorbing boundary condition (WABC). Sensors are located horizontally on the top surface and vertically at the right-hand side of a domain. The problem con/uniFB01guration is shown in Fig. /one.pnum. /three.pnum /two.pnum./one.pnum The governing equation The governing equation for the anti-plane shear wave propagation in an undamped solid domain is (for brevity, the spatial and/or temporal dependencies of variables are omitted): 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, and J= [0;T]is the time interval of interest. The solid is subject to a traction-free condition on the top surface ( t): @u @y(x;0;t) = 0; 0xL: (/two.pnum) The absorbing boundary conditions are presented on the bottom ( b), left ( l) and right ( r) boundaries, respectively as: @u @y(x;D;t ) =f(x;D;t ) G1 Vs@u @t; 0xL; (/three.pnum) @u @x(0;y;t) =f(0;y;t) G1 Vs@u @t; D y0; (/four.pnum) @u @x(L;y;t ) =1 Vs@u @t; D y0; (/five.pnum) whereDis they-coordinate of b;Lis thex-coordinate of r; andVs(x;y)is the shear wave speed of the medium. We note that f(x;y;t )denotes the dynamic traction of shear stress applied in the z-direction on the boundaries, and f(x;y;t )is used for modeling the incoming seismic wave/uniFB01eld arriving at the WABC boundary/one.pnum: f(x;y;t ) =2G Vs@uI @t; (/six.pnum) whereuIdenotes the incident seismic wave/uniFB01eld. We consider that uIimpinges the domain at land rwith an inclined angle. The zero initial-value conditions are presented as the system is initially at rest: u(x;y;0) = 0;@u @t(x;y;0) = 0: (/seven.pnum) /three.pnum Forward Wave Modeling This section presents the /uniFB01nite element modeling for obtaining the numerical solutions of the governing equation (/one.pnum). /three.pnum./one.pnum Finite element method to solve the state problem From the governing equation (/one.pnum), the weak form of the state problem is cast in order to /uniFB01nd u(t)2Ufor allv2V: Z rv(Gru)d +Z rvG Vs@u @td +Z lvG Vs@u @td +Z bvG Vs@u @td +Z v@2u @t2d =Z bvf(x;y;t )d +Z lvf(x;y;t )d; (/eight.pnum) wherevdenotes a test function. The function spaces, UandV, are de/uniFB01ned as: U=fu:u2H1g;V=fv:v2H1g (/nine.pnum) /one.pnumWe note that Bielak’s Domain Reduction Method (DRM) Bielak et al. (/two.pnum/zero.pnum/zero.pnum/three.pnum) can be used for modeling the wave responses in a truncated domain induced by an incident wave/uniFB01eld, impinging a WABC of a domain. In this work, we are testing our inverse algorithm in a truncated domain by using the boundary conditions (/three.pnum-/six.pnum) as an alternative to the DRM. /four.pnum Next, the test and trial function are approximated by the standard /uniFB01nite-element approximation as follows: u(x;y;t )' (x;y)Tu(t); v(x;y)'vT (x;y); (/one.pnum/zero.pnum) where (x;y)denotes a vector of global basis functions that is made of shape functions in the local co- ordinate of each element in the domain; vdenotes a vector of nodal solution of the test function; and u(t)denotes the nodal solution vector of the state problem at time t. Then, the weak form (/eight.pnum) reduces to: Mu(t) +C_u(t) +Ku(t) =F(t); (/one.pnum/one.pnum) where ( _ ) and(  ) denotes, respectively, the /uniFB01rst and second-order derivative of its subtended variable with respect to t;Mdenotes a global mass matrix; Cdenotes a global dumping matrix; Kdenotes a global sti/uniFB00ness matrix; and Fdenotes a global force vector. The vectors and matrices are de/uniFB01ned as: K=Z G@ @x@ T @x+@ @y@ T @y d ; C=Z rG Vs Td +Z lG Vs Td +Z bG Vs Td; M=Z  Td ; F=Z b fd +Z l fd: (/one.pnum/two.pnum) /three.pnum./two.pnum Time integration of the state problem We solve the time-dependent discrete form (/one.pnum/one.pnum) in every discrete j-th time step as: Muj+C_uj+Kuj=Fj: (/one.pnum/three.pnum) Applying the initial value-conditions (/seven.pnum) onto (/one.pnum/three.pnum), the solution vector at the initial time step is obtained by solving the following: Mu0=F0; (/one.pnum/four.pnum) By using the implicit Newmark time integration, this work solves the system of equations for each time step. The solution vector of a next time step is related to that of a previous time step as: _uj=2 t(ujuj1)_uj1; (/one.pnum/five.pnum) and uj=4 (t)2(ujuj1)4 t_uj1uj1; (/one.pnum/six.pnum) where tdenotes the size of a time step. By plugging (/one.pnum/five.pnum) and (/one.pnum/six.pnum) into (/one.pnum/three.pnum), it turns into the follow- ing: 4 (t)2M+2 tC+K uj=Fj+M4 (t)2uj1+4 t_uj1+uj1 +C2 tuj1+_uj1 :(/one.pnum/seven.pnum) By using (/one.pnum/seven.pnum), this work solves for uj. Then, the values of _ujandujcan be updated by using ujand, respectively, (/one.pnum/five.pnum) and (/one.pnum/six.pnum). /three.pnum./three.pnum Compact form of the state problem Equations (/one.pnum/four.pnum) to (/one.pnum/seven.pnum) can be cast in the following compact form: Q^u=^F; (/one.pnum/eight.pnum) where the vector ^ucorresponds to the space-time discretization of u(x;y;t )for(x;y)2 andt2[0;T], induced byf(x;y;t ); and ^Fis a global force vector corresponding to all the time steps. These vectors are built as: /five.pnum ^u=2 6666666664u0 _u0 u0 ... uN _uN uN3 7777777775;^F=2 66666666640 0 F0 ... FN 0 03 7777777775; (/one.pnum/nine.pnum) where the total number of time steps is de/uniFB01ned as N. Finally, the matrix Qis the discrete forward operator de/uniFB01ned as: 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 C M 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; (/two.pnum/zero.pnum) where: a0=4 (t)2; a1=2 t; a2=4 t; Ke =a0M+a1C+K;L1=a0Ma1C;L2=a2MC;L3=M: (/two.pnum/one.pnum) In (/two.pnum/zero.pnum), the /one.pnumst and /two.pnumnd rows correspond to the initial conditions u0= 0 and_u0= 0; the /three.pnumrd row is equivalent to (/one.pnum/four.pnum); the /four.pnum, /seven.pnum, /one.pnum/zero.pnum, ...-th rows correspond to (/one.pnum/seven.pnum); the /five.pnum, /eight.pnum, /one.pnum/one.pnum, ...-th rows are (/one.pnum/five.pnum); and the /six.pnum, /nine.pnum, /one.pnum/two.pnum, ...-th rows are (/one.pnum/six.pnum). Instead of the time marching scheme, by solving the matrix system of (/one.pnum/eight.pnum) "all at once", we could obtain the forward wave solutions. However, to save the computational cost, we opt for using the time marching scheme for solving (/one.pnum/eight.pnum). /four.pnum Inverse Modeling This section presents mathematical modeling for identifying the temporal and spatial distributions of un- known dynamic traction f(x;y;t )—mimicking unknown incoming seismic wave/uniFB01eld per (/six.pnum)—based on measured wave responses while the properties of the system, such as the extents and the material prop- erties of the solid, are assumed to be known. Speci/uniFB01cally, the components of global force vectors corre- sponding to targeted traction will be identi/uniFB01ed. To this end, the presented work studies the discretize-then-optimize (DTO) approach because, in contrast to the OTD counterpart, the DTO method involves quite a compact mathematical modeling procedure to side- impose the governing wave equation to an objective functional and to derive the derivatives of a Lagrangian functional with respect to adjoint and control variables. Namely, an objective functional is formulated in a discrete form, and, then, the discrete optimality conditions are derived and implemented. /four.pnum./one.pnum Control parameters Under this inversion method, we identify Pkj—a component of ^Fcorresponding to kandtj—as control parameters, where kis thek-th discrete node on the boundaries landb, andtjis thej-th time step. The discrete numbering of kincreases counter-clockwise from the top-left corner of the domain. Fig. /two.pnum shows the discrete numbering kinPkjwith a reduced number of discrete nodes. We note that the components of^Fthat do not correspond to discrete nodes on landbvanish, and they will not be inverted for. /six.pnum NodesFigure /two.pnum: Discrete numbering kinPkjcorresponding the discrete node kon the boundaries landb. /four.pnum./two.pnum Discrete Lagrangian functional As a minimization problem, we seek the values of the control parameters (i.e., Pkjfor allkandj) that correspond to a minimum of a discrete objective functional, which is de/uniFB01ned as: ^L=1 2(^u^um)TB(^u^um); (/two.pnum/two.pnum) where ^ucorresponds to the space-time discretization of a forward simulation variable induced by a set of estimatedPkj;^umis the space-time discretization of um(x;y;t )induced by a set of Pkjof to targeted force vector. In this computational study, we synthetically create the measured displacement /uniFB01eld, um(x;y;t ), by using our FEM solver with a pseudo target force vector. In (/two.pnum/two.pnum), Bis de/uniFB01ned as tB, where Bis a square matrix, of which components are all /zero.pnum except for those of the diagonal, having values of all /one.pnum, if they correspond to the degrees of freedom of sensor locations. In this work, the measured data is obtained from sparsely-distributed sensors. To tackle the minimization problem, the following Lagrangian functional ^Ais built by imposing the discrete form of the forward wave problem, using a Lagrange multiplier vector ^, onto ^L: ^A=1 2(^u^um)TB(^u^um)^T(Q^u^F); (/two.pnum/three.pnum) where ^is the discrete (space-time) Lagrange multiplier: ^=h T 0;_T 0;T 0;::;T N;_T N;T N;iT : (/two.pnum/four.pnum) We show that the objective functional (/two.pnum/two.pnum) is quadratic and convex (please see the Appendix A) so that the /uniFB01rst-order optimality conditions of the Lagrangian functional (/two.pnum/three.pnum) could be su/uniFB03cient for obtaining an inver- sion solution. Thus, we do not need to include any regularization in the presented work: our previous work Guidio and Jeong (/two.pnum/zero.pnum/two.pnum/zero.pnum), considering a /two.pnumD bounded domain, also showed that the Tikhonov (TN) regular- ization does not improve the performance of inverting dynamic traction on a boundary of /two.pnumD bounded solid domain. However, as shown in the previous work, because of the sparsity of the measurement data, the inversion solution may not be unique in the presented work unless a su/uniFB03ciently large number of sensors measure wave motions: We discuss this aspect in Example /two.pnum in the presented work. /four.pnum./three.pnum The /uniFB01rst-order optimality conditions To identify unknown, target control parameters, the following optimality conditions should be satis/uniFB01ed: @^A @^= 0 The /uniFB01rst condition (state problem) ; (/two.pnum/five.pnum) @^A @^u= 0 The second condition (adjoint problem) ; (/two.pnum/six.pnum) @^A @^F= 0 The third condition (control problem) : (/two.pnum/seven.pnum) /seven.pnum The /uniFB01rst condition, (/two.pnum/five.pnum), the vanishing derivative of ^Awith respect to ^, will be automatically satis/uniFB01ed when we solve the discrete forward problem (/one.pnum/eight.pnum), which is also referred to as the discrete state equation: @^A @^=Q^u+^F= 0: (/two.pnum/eight.pnum) As the second condition, (/two.pnum/six.pnum), the derivate of ^Awith respect to ^ushould vanish, and it will be automatically satis/uniFB01ed when we solve the following discrete adjoint problem: @^A @^u=@^L @^u|{z} a@(^TQ^u) @^u|{z} b= 0: (/two.pnum/nine.pnum) The partain (/two.pnum/nine.pnum) can be written as (see the derivation in Appendix B): a=B(^u^um): (/three.pnum/zero.pnum) The partbin (/two.pnum/nine.pnum) can be written as: b=@ @^u(^TQ^u)|{z} scalar=@ @^u(^TQ^u)T =@ @^u((Q^u)T^) =@ @^u(^uTQT^) =QT^: (/three.pnum/one.pnum) Due to (/three.pnum/zero.pnum) and (/three.pnum/one.pnum), (/two.pnum/nine.pnum) can be written as: @^A @^u=QT^+B(^u^um) = 0|{z} adjoint equation: (/three.pnum/two.pnum) Because (/three.pnum/two.pnum) includes transpose of Q, we solve it by marching backwards in time. Namely, from the last two rows of QT, we obtain the /uniFB01nal conditions: N= 0;_N= 0; (/three.pnum/three.pnum) respectively; and the third row from the bottom of QTresults in the following equation: Ke/uniFB00TN= tB(uNumN) +a1_N+a0N; (/three.pnum/four.pnum) which can be solved for N. For time steps j=N1;N2;:::; 1;we /uniFB01rst update the followings: j=MTj+1j+1; (/three.pnum/five.pnum) _j= (a2MT+CT)j+1_j+1a2j+1: (/three.pnum/six.pnum) Then, solve the following: Ke/uniFB00Tj=tB(ujumj) +a1_j+a0j+ (a0MT+a1CT)j+1a1_j+1a0j+1: (/three.pnum/seven.pnum) Eventually, the /uniFB01rst three rows of QTlead to the followings: MT0=MT11; (/three.pnum/eight.pnum) _0=CT0+ (a2MT+CT)1_1a21; (/three.pnum/nine.pnum) 0=KT0+ (a0MT+a1CT)1a1_1a01+ tB(u0um0): (/four.pnum/zero.pnum) As the third condition, (/two.pnum/seven.pnum), the vanishing derivative of ^Awith respect to ^Fis tantamount to the following control equation: @^A @^F=^= 0; (/four.pnum/one.pnum) /eight.pnum which implies that the inversion solution of ^Fshould lead to the vanishing ^. In other words, if the inversion solution of ^Fmatches its target, ^uwill match ^umat sensor locations. Thus, the driving force term, B(^u ^um)in (/three.pnum/two.pnum) vanishes so that ^will vanish per (/three.pnum/two.pnum). Eq. (/four.pnum/one.pnum) also implies that@^A @(=Pkj)=@^L @(=Pkj)is the component of the vector ^corresponding to the global node numbering and the time step of =Pkj. In other words, (/four.pnum/one.pnum) becomes: r(=Pkj)^A=r(=Pkj)^L=(( k;tj)whenj= /one.pnum, ( k;tj)whenj/two.pnum,(/four.pnum/two.pnum) where ( k;tj)or( k;tj)denotes the component of ^corresponding to the degree of freedom and the timestep ofPkj. Eq. (/four.pnum/two.pnum) implies that this works uses the solution of the adjoint problem for evaluating the gradient of ^Lwith respect to in (/four.pnum/two.pnum) in a semi-analytical manner. In addition, it also implies that searching a set of the control parameters, , that leads to the vanishing ^is equivalent to searching leading to the zero-valuedr^L. This process is equivalent to the gradient-based minimization. The gradient-based minimization process, shown in the Algorithm /one.pnum. /five.pnum Numerical Implementation of the Inversion Process By utilizing the semi-analytically evaluated, reduced/two.pnumgradient vectorr^L, this work iteratively updates a set of estimated control parameters by using the gradient-based minimization scheme as follows: (a) First, we use pseudo-traction f(x;y;t ), which corresponds to an inclined incident wave, to create a pseudo targeted force vector ^Fand generate synthetic measured data ^umat sensors. (b) Then, ^uis obtained by using estimated ^Fthat is constructed by estimated control parameters (i.e., Pkj). (c) Next, the adjoint problem is solved (i.e., computing ^) by using ^uthat is calculated in the previous step. (d) By using the adjoint solution, we compute the gradient of the objective functional, r^L. (e) The conjugate-gradient method Guidio and Jeong (/two.pnum/zero.pnum/two.pnum/zero.pnum) determines the best search direction, and an optimal step length is calculated by Newton’s method Guidio and Jeong (/two.pnum/zero.pnum/two.pnum/zero.pnum). (f) Lastly, the gradient-based minimization scheme updates the approximate control parameters using the search direction and optimal step length. The numerical optimizer repeats the steps (b) to (f) above, and a set of these steps is counted as an inver- sion iteration. Throughout inversion iterations, our numerical optimizer iteratively calculate the values of estimated control parameters that lead to the vanishing gradient of the objective functional. The compre- hensive method for the numerical optimizer is outlined in the Algorithm /one.pnum. /five.pnum./one.pnum Updating the estimated control parameters Starting from the initially-guessed (zero-valued) control parameter vector consisting of Pkj, the estima- tion ofcan be modi/uniFB01ed iteratively as follows: (s+1)=(s)+h(s)d(s); (/four.pnum/three.pnum) whereh(s)is a scalar-valued step size and d(s)is a search-direction vector computed by the conjugate- gradient scheme Guidio and Jeong (/two.pnum/zero.pnum/two.pnum/zero.pnum), and sis the iteration index in the optimization process. For each iteration,s, the Newton’s method Guidio and Jeong (/two.pnum/zero.pnum/two.pnum/zero.pnum); Lloyd and Jeong (/two.pnum/zero.pnum/one.pnum/eight.pnum) is used to determine the optimal step size, h(s). /six.pnum Numerical Experiments In this section, we investigate the performance of the presented inversion method in a domain truncated by using a WABC in consideration of various factors. /two.pnumThe full gradient ish @^L @^(= 0);@^L @^u(= 0);@^L @^i . /nine.pnum Algorithm /one.pnum Inversion Algorithm /one.pnum: Set TOL= 107and maximum iteration number = /one.pnum/zero.pnum/zero.pnum/zero.pnum. /two.pnum: Build M,C, andKmatrices. /three.pnum: Create synthetic measurement data ^umfrom a forward wave solver using a targeted force vector ^Fcorresponding to targeted tractionf(x;y;t ). /four.pnum: Set an iteration index s= 1 and initial control parameters ((s=1)= 0) and compute ^L(s=1). /five.pnum:while (^L(s+1)>TOL^L(1)ands< (maximum iteration number)) do /six.pnum: Solve the discrete forward problem using (s)and compute ^u. /seven.pnum: Obtain ^by solving the discrete 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 Newton’s method. /one.pnum/one.pnum: Update estimated control parameters as (s+1)=(s)+d(s)h(s)and compute ^L(s+1). /one.pnum/two.pnum:s s+ 1. /one.pnum/three.pnum: if(s== (maximum iteration number) and^L(s10)^L(s) ^L(s10)1001%)then /one.pnum/four.pnum: maximum iteration number = maximum iteration number +/one.pnum/zero.pnum/zero.pnum/zero.pnum. /one.pnum/five.pnum: end if . /one.pnum/six.pnum:end while . Figure /three.pnum: Material pro/uniFB01le /two.pnum: a /two.pnum-layered solid with two inclusions. The considered domain size is set to be /two.pnum/zero.pnum/zero.pnum m /six.pnum/zero.pnum m, and the mass density is uniform as /three.pnum/zero.pnum/zero.pnum/zero.pnum kg/m3in the domain. Two material pro/uniFB01les, which di/uniFB00er from each other in terms of wave speeds, are considered in the presented numerical experiments. The material pro/uniFB01le /one.pnum is of a homogeneous background solid domain with two inclusions as in Fig. /one.pnum with shear wave speeds of Vs1= /three.pnum/five.pnum/three.pnum./five.pnum/five.pnum m/s, Vs2= /eight.pnum/zero.pnum/zero.pnum m/s, and Vs3= /one.pnum/zero.pnum/zero.pnum/zero.pnum m/s. The material pro/uniFB01le /two.pnum is a /two.pnum-layered background solid with the same two inclusions, as shown in Fig. /three.pnum, with shear wave speeds of Vs1= /three.pnum/five.pnum/three.pnum./five.pnum/five.pnum m/s, Vs2= /two.pnum/two.pnum/three.pnum./six.pnum m/s, Vs3= /eight.pnum/zero.pnum/zero.pnum m/s, and Vs4= /one.pnum/zero.pnum/zero.pnum/zero.pnum m/s. The synthetic measured data ^umis computed by using an element size of /zero.pnum./five.pnum m for a targeted force vector. To avoid an inverse crime/three.pnum, we use an element size of /one.pnum m for evaluating state solution ^uand the adjoint solution ^for an estimated force vector. t of /zero.pnum./zero.pnum/zero.pnum/one.pnum s is used for computing all ^um,^u, and ^. For a homogeneous background domain (i.e., the material pro/uniFB01le /one.pnum), four targeted force vectors ^F1;2;3;4 (corresponding to four di/uniFB00erent targeted dynamic tractions) are considered in the presented numerical experiments. The dynamic tractions are applied on the left and bottom boundaries of a homogeneous background domain, and they mimic incident plane wave/uniFB01eld propagating with an incident angle 45. The time-dependent values of targeted ^F1;2;3corresponding to a node at a speci/uniFB01c location x= /one.pnum/zero.pnum/zero.pnum m and y = -/six.pnum/zero.pnum m—Ricker wavelets with di/uniFB00erent frequency contents—are shown in Fig. /four.pnum(a-f). The signals corre- sponding to ^F1;2;3have the central frequency of /five.pnum, /one.pnum/zero.pnum, /one.pnum/five.pnum Hz, respectively. The total observation duration is /one.pnum./five.pnum s in the inversion simulation of identifying targeted ^F1;2;3. On the other hand, the time-dependent signal corresponding to ^F4at the same location and its frequency content are shown in Fig. /four.pnum(g,h). This signal is a portion of a registered ground motion signal during the /one.pnum/nine.pnum/nine.pnum/four.pnum Northridge earthquake from the Paci/uniFB01c Earthquake Engineering Research Center (PEER) ground motion database PEER (/two.pnum/zero.pnum/zero.pnum/zero.pnum). The total observation duration is /six.pnum s for the inversion simulation to identify the targeted ^F4. On the other hand, ^F5is a force vector due to dynamic traction on the left and bottom boundaries of a two- layered background domain (i.e., the material pro/uniFB01le /two.pnum), and this signal mimics incident plane wave/uniFB01eld of an incident angle 45in the bottom layer. The signal corresponding to ^F5atx= /one.pnum/zero.pnum/zero.pnum m and y= -/six.pnum/zero.pnum m is a Ricker wavelet with a central frequency of /five.pnum Hz. /three.pnumWe should replicate a situation where the wave solver for computing ^uis only approximating its true physics that governs ^um. /one.pnum/zero.pnum (a)(b) (c)(d) (e)(f) (g)(h)Figure /four.pnum: (a,c,e,g) The time signals of ^F1;2;3;4corresponding to a node at x= /one.pnum/zero.pnum/zero.pnum m and y= -/six.pnum/zero.pnum m, respec- tively; (b,d,f,h) the amplitude of Fourier Transform of each signal. /one.pnum/one.pnum Table /one.pnum: Summary of all cases. Case Material Target Vertical # of Noise Terminal Error number pro/uniFB01le force pro/uniFB01le array sensors level error E Eu /one.pnum /one.pnum ^F1N /five.pnum/zero.pnum /zero.pnum% /seven.pnum./six.pnum/zero.pnum% /zero.pnum./eight.pnum/nine.pnum% /two.pnum /one.pnum ^F1Y /six.pnum/five.pnum /zero.pnum% /one.pnum./four.pnum/five.pnum% /zero.pnum./one.pnum/six.pnum% /three.pnum /one.pnum ^F1Y /three.pnum/two.pnum /zero.pnum% /one.pnum./six.pnum/zero.pnum% /zero.pnum./one.pnum/eight.pnum% /four.pnum /one.pnum ^F1Y /one.pnum/six.pnum /zero.pnum% /one.pnum/eight.pnum./six.pnum/one.pnum% /eight.pnum./eight.pnum/four.pnum% /five.pnum /one.pnum ^F2Y /six.pnum/five.pnum /zero.pnum% /two.pnum./eight.pnum/two.pnum% /zero.pnum./six.pnum/two.pnum% /six.pnum /one.pnum ^F3Y /six.pnum/five.pnum /zero.pnum% /one.pnum/zero.pnum./nine.pnum/zero.pnum% /three.pnum./four.pnum/five.pnum% /seven.pnum /one.pnum ^F4Y /six.pnum/five.pnum /zero.pnum% /zero.pnum./three.pnum/three.pnum% /zero.pnum./zero.pnum/zero.pnum% /eight.pnum /two.pnum ^F5Y /six.pnum/five.pnum /zero.pnum% /one.pnum./five.pnum/six.pnum% /zero.pnum./one.pnum/six.pnum% /nine.pnum /two.pnum ^F5Y /six.pnum/five.pnum /one.pnum% /two.pnum./eight.pnum/eight.pnum% /zero.pnum./one.pnum/eight.pnum% /one.pnum/zero.pnum /two.pnum ^F5Y /six.pnum/five.pnum /two.pnum% /seven.pnum./one.pnum/zero.pnum% /zero.pnum./two.pnum/zero.pnum% /one.pnum/one.pnum /two.pnum ^F5Y /six.pnum/five.pnum /five.pnum% /four.pnum/one.pnum./zero.pnum/two.pnum% /zero.pnum./four.pnum/zero.pnum% /one.pnum/two.pnum /two.pnum ^F5Y /six.pnum/five.pnum /one.pnum/zero.pnum% /one.pnum/one.pnum/two.pnum./four.pnum/zero.pnum% /zero.pnum./eight.pnum/six.pnum% For the inversion process, the temporal and spatial intervals for the discretization are /zero.pnum./zero.pnum/zero.pnum/one.pnum s and /one.pnum m, respectively, for computing wave responses induced by all the estimated ^Fin all the cases. Cases /one.pnum-/six.pnum and /eight.pnum-/one.pnum/one.pnum, which consider ^F1;2;3;5, use /three.pnum/nine.pnum/one.pnum,/seven.pnum/nine.pnum/one.pnum control parameters (i.e., /two.pnum/six.pnum/one.pnum (over space) /one.pnum/five.pnum/zero.pnum/one.pnum (over time)), and Case /seven.pnum, which consider ^F4, uses /one.pnum,/five.pnum/six.pnum/six.pnum,/two.pnum/six.pnum/one.pnum control parameters (i.e, /two.pnum/six.pnum/one.pnum (over space) /six.pnum/zero.pnum/zero.pnum/one.pnum (over time)). The control parameters are set as zero at the initial inversion iteration. In what follows, six examples of numerical investigations are presented: •The /uniFB01rst example is aimed at comparing the performance of inverting Pkj, corresponding to ^F1, between Case /one.pnum—utilizing only the sensors on the top surface (t)of the domain—versus Case /two.pnum— using a vertical array of sensors on the right boundary (r)with those on t. •The second example tests the inversion performance of identifying ^F1with respect to the number of sensors on randtin Cases /two.pnum to /four.pnum. •The third example is focused on analyzing the performance of inverse modeling to reconstruct tar- geted force vectors, of which time-dependent signals contain di/uniFB00erent central frequencies: Cases /two.pnum, /five.pnum, and /six.pnum use ^F1,^F2, and ^F3, respectively, corresponding to tractions with central frequencies of /five.pnum, /one.pnum/zero.pnum, and /two.pnum/zero.pnum Hz. •The fourth example investigates the inversion performance of identifying ^F4, which corresponds to a realistic seismic signal, in Case /seven.pnum. Namely, the performance in Case /seven.pnum will be compared with that in Case /two.pnum. •The /uniFB01fth example examines the convergence of the estimated control parameters into targeted ^F5 by using the material pro/uniFB01le /two.pnum (the /two.pnum-layered background solid with two inclusions) in Case /eight.pnum. •Lastly, the sixth example studies the e/uniFB00ect of noise level in measurement data on the inversion per- formance of identifying ^F5. In addition to Case /eight.pnum, where a noise level of /zero.pnum% is used, Case /nine.pnum-/one.pnum/two.pnum are considered and correspond to noise levels of /one.pnum%, /two.pnum%, /five.pnum%, and /one.pnum/zero.pnum%, respectively. The summary of all the cases presented in this study is shown in Table /one.pnum. For appraising the e/uniFB00ectiveness to reconstruct ^Fin the numerical results, the following error norm be- tween estimated force vector and its target is used: E=j^Ftarget^Festimatej2 j^Ftargetj2100[%]; (/four.pnum/four.pnum) where we compare the components of ^Ftarget and^Festimate corresponding to the nodes on the boundaries where tractions are applied./four.pnumFor the sake of assessing the accuracy to reconstruct the wave responses /four.pnumMore speci/uniFB01cally, they are compared with each other at the nodes of the domain that is discretized for simulating the wave responses for ^Festimate . Since a twice smaller element size is used for computing the wave responses for ^Ftarget , the components of ^Ftarget at those nodes are multiplied twice when they are compared with ^Festimate in (/four.pnum/four.pnum). /one.pnum/two.pnum in the entire computational domain, the numerical results show the following error norm between the wave response solution vector at the j-th time step, uj, induced by /uniFB01nal ^Festimate and its counterpart, umj, induced by ^Ftarget for all the time steps/five.pnum: Eu=NX j=1jumjujj2 jumjj2100[%]: (/four.pnum/five.pnum) /six.pnum./one.pnum The exemplary forward wave responses Fig. /five.pnum(a-d) show the wave responses in the whole computational domain induced by targeted ^F1in Case /one.pnum att= /zero.pnum./three.pnum s, /zero.pnum./four.pnum s, /zero.pnum./five.pnum s, and /zero.pnum./six.pnum s, respectively. Fig. /six.pnum presents the wave responses caused by targeted ^F5 in Case /eight.pnum, in which material pro/uniFB01le /two.pnum is employed, at the same time steps. The plots at t= /zero.pnum./three.pnum s in Fig. /five.pnum(a) and Fig. /six.pnum(a) show that the dynamic input traction at the WABC mimics seismic incident, inclined plane waves. Fig. /six.pnum(a,b) shows the change in the incident angle of the wave caused by the di/uniFB00erence of the shear wave speeds in the two-layered medium. This behavior complies with Snell’s Law. Fig. /five.pnum(b) and Fig. /six.pnum(b) display the waves passing through the /uniFB01rst inclusion, which is located at /six.pnum/zero.pnum m from the left boundary. Fig. /five.pnum(c) and Fig. /six.pnum(c) show the wave motions passing through the second inclusion, which is /one.pnum/two.pnum/zero.pnum m from the left boundary. Finally, Fig. /five.pnum(d) and Fig. /six.pnum(d) exhibit the waves arriving at the right boundary. Fig. /seven.pnum shows the comparison between the components of targeted ^F1and those of a vector2G V@u @t—per (/six.pnum)—on l,b, and rin a homogeneous domain induced by the same ^F1. As shown in Fig. /seven.pnum(b), they are in excellent agreement with each other except for those that are attributed to the waves re/uniFB02ected from t, in turn, arriving at landb. Therefore, our underlying rationale of (/six.pnum) is justi/uniFB01ed. /six.pnum./two.pnum The veri/uniFB01cation of the inverse modeling The theoretical derivation of the adjoint and control equations and the numerical implementation are veri/uniFB01ed prior to the parametric studies. The gradient obtained by our semi-analytical DTO approach is compared with that from the /uniFB01nite di/uniFB00erence (FD) procedure, which uses the following equation: @^L @FD=^L(+h)^L() h; (/four.pnum/six.pnum) whereindicates a control parameter. In this veri/uniFB01cation, we used the targeted ^F1, material pro/uniFB01le /one.pnum, and only one sensor in the center of t. Fig. /eight.pnum reveals an outstanding agreement between the normalized gradi- ent (i.e.,r^L=jr^Lj) that is calculated by using the presented DTO approach and that by the FD approach. Both of them are computed at the /uniFB01rst inversion iteration for control parameters corresponding to a node atx= /one.pnum/zero.pnum/zero.pnum m and y= -/six.pnum/zero.pnum m. Therefore, the theoretical derivation and mathematical implementation are trustable. /six.pnum./three.pnum Example /one.pnum: Investigating the inversion performance with or without a vertical ar- ray of sensors In this example, we compare the performance of the presented inverse modeling of identifying a targeted set ofPkjof^F1in Case /one.pnum, where we use the sensors only on the top surface t, with that in Case /two.pnum, where we add a vertical (e.g., at borehole) array of sensors on the boundary on the right side r. We used the material pro/uniFB01le /one.pnum and /five.pnum/zero.pnum sensors on tin both Cases. In Case /two.pnum, a vertical array of /one.pnum/five.pnum sensors are added on r. Fig. /nine.pnum shows the targeted set of Pkjand its reconstructed counterpart obtained by using sensors on only tversus both tandr, in Case /one.pnum and /two.pnum, respectively. When Case /one.pnum is tested, as shown in Fig. /nine.pnum(b), the targeted control parameters on the bottom boundary b(especially, those at the right-bottom portion ofb) cannot be fully recovered. However, by adding the vertical array of sensors on r, they are fully reconstructed (Fig. /nine.pnum(c)). /five.pnumAgain, the wave responses ujandumjare compared with each other at the nodes of the domain that is discretized for simulating the wave responses for ^Festimate . /one.pnum/three.pnum (a) (b) (c) (d)Figure /five.pnum: (a-d)The wave motions of u(x;y;t )[m], att= /zero.pnum./three.pnum s, /zero.pnum./four.pnum s, /zero.pnum./five.pnum s, and /zero.pnum./six.pnum s, respectively, in Case /one.pnum. /one.pnum/four.pnum (a) (b) (c) (d)Figure /six.pnum: (a-d) The wave motions of u(x;y;t )[m], att= /zero.pnum./three.pnum s, /zero.pnum./four.pnum s, /zero.pnum./five.pnum s, and /zero.pnum./six.pnum s, respectively, in Case /eight.pnum. /one.pnum/five.pnum ReflectionReflectionSame as the applied tractionFigure /seven.pnum: (a) A set of Pkjcorresponding to targeted traction and (b)2G V@u @tof the wave responses induced by the targeted traction in N/m2in Case /one.pnum. 02004006008001000120014001600-1-0.8-0.6-0.4-0.200.2 Figure /eight.pnum: Comparasion between the gradients generated by the presented inverse modeling and the FD approximation. /one.pnum/six.pnum (a) (b)(c)not accurately reconstructedFigure /nine.pnum: Example /one.pnum: (a) A set of Pkjcorresponding to targeted traction and (b,c) The sets of reconstructed Pkjfor Case /one.pnum and /two.pnum at the /one.pnum/zero.pnum/zero.pnum/zero.pnumth iteration. The dominant frequency of /five.pnum Hz is used. The horizontal axis in the contour plot represents the nodes on the absorbing boundaries of the domain that is discretized for simulating the wave responses for ^Festimate . Moreover, since the element size used to calculate the wave responses for ^Ftarget is twice smaller, the targeted force components are twice multiplied at those nodes. We repeat this in the other contour plots where the targeted set of Pkjand its reconstructed counterpart are compared. /one.pnum/seven.pnum 010020030040050060070080090010000102030405060708090100Figure /one.pnum/zero.pnum: Example /one.pnum - After /one.pnum/zero.pnum/zero.pnum/zero.pnum iterations, the /uniFB01nal value of Eis /seven.pnum./six.pnum/zero.pnum% in Case /one.pnum, and that in Case /two.pnum is /one.pnum./four.pnum/five.pnum%, which shows a reduction of /six.pnum./one.pnum/five.pnum% in E. Fig. /one.pnum/zero.pnum shows that, by using the vertical array of sensors in Case /two.pnum, the terminal value of Eof the inversion becomes much smaller than that in Case /one.pnum. After /one.pnum/zero.pnum/zero.pnum/zero.pnum iterations, the /uniFB01nal value of Eis /seven.pnum./six.pnum/zero.pnum% in Case /one.pnum, but that in Case /two.pnum is /one.pnum./four.pnum/five.pnum%, which shows a reduction of /six.pnum./one.pnum/five.pnum% in E. Thus, we suggest that using not only the sensors on the top surface but also a vertical array of sensors improves the inversion performance in the presented problem setting compared to the case using only the sensors on the top surface. As shown in Table /one.pnum, even though the terminal error Eis /seven.pnum./six.pnum/zero.pnum% in Case /one.pnum, the error Eubetween the wave responses induced by targeted and /uniFB01nal estimated force vector is only /zero.pnum./eight.pnum/nine.pnum% in Case /one.pnum. It is because the wave induced by a particular set of Pkjat the right bottom corner passes rimmediately after its onset. /six.pnum./four.pnum Example /two.pnum: Investigating the inversion performance with respect to the number of sensors In this example, we test the performance of identifying a targeted set of Pkjof^F1with respect to the number of sensors on tandr. We considered Cases /two.pnum, /three.pnum, and /four.pnum, which use /six.pnum/five.pnum sensors (/five.pnum/zero.pnum on tand /one.pnum/five.pnum onr), /three.pnum/two.pnum sensors ( /two.pnum/five.pnum on tand /seven.pnum on r), and /one.pnum/six.pnum sensors (/one.pnum/two.pnum on tand /four.pnum on r), respectively, and material pro/uniFB01le /one.pnum. Fig. /one.pnum/one.pnum shows the targeted set of Pkjcorresponding to ^F1and its counterparts obtained in Cases /two.pnum, /three.pnum, and /four.pnum. Fig. /one.pnum/one.pnum shows that when Cases /two.pnum and /three.pnum are tested (/six.pnum/four.pnum and /three.pnum/two.pnum sensors are used), the inversion algorithm can accurately reconstruct the targeted control parameters. However, in Case /four.pnum, where /one.pnum/six.pnum sensors are used, the control parameters cannot be accurately recovered. As shown in Fig. /one.pnum/two.pnum, while the terminal values of the errors Ein Case /two.pnum and /three.pnum, where /six.pnum/five.pnum and /three.pnum/two.pnum are used, are /one.pnum./four.pnum/five.pnum% and /one.pnum./six.pnum%, the error Eis increased to /one.pnum/eight.pnum./six.pnum/one.pnum% in Case /four.pnum, where /one.pnum/six.pnum sensors are utilized. Table /one.pnum shows that the error Eubetween the wave responses induced by the target and its /uniFB01nal-estimated counterpart also grows as the number of sensors is decreased. For instance, the error Eugrows from /zero.pnum./one.pnum/six.pnum% in Case /two.pnum to /eight.pnum./eight.pnum/four.pnum% in Case /four.pnum. Fig. /one.pnum/three.pnum shows the wave responses’ snapshots in the entire computational domain obtained in Case /four.pnum due to targeted and reconstructed control parameters. Although the two wave responses look overall alike in Fig. /one.pnum/three.pnum, the value of Euin Case /four.pnum is the highest among all the cases /one.pnum to /one.pnum/two.pnum. Thus, we suggest that the desirable ratio of the size of measurement data to the number of control pa- rameters can be as small as /one.pnum:/eight.pnum—e.g., (/three.pnum/two.pnum sensors /one.pnum/five.pnum/zero.pnum/zero.pnum timesteps):(/two.pnum/six.pnum/one.pnum nodes /one.pnum/five.pnum/zero.pnum/zero.pnum timesteps)—in the presented example. /one.pnum/eight.pnum (a) (c)(b) (d)Figure /one.pnum/one.pnum: Example /two.pnum: (a) The set of Pkjcorresponding to targeted traction and (b,c,d) the sets of recon- structedPkjfor Case /two.pnum, /three.pnum, and /four.pnum, respectively, at the /one.pnum/zero.pnum/zero.pnum/zero.pnumth iteration. 010020030040050060070080090010000102030405060708090100 Figure /one.pnum/two.pnum: Example /two.pnum - After /one.pnum/zero.pnum/zero.pnum/zero.pnum iterations, the /uniFB01nal value of Eis /one.pnum./four.pnum/five.pnum% in Case /two.pnum, /one.pnum./six.pnum/zero.pnum% in Case /three.pnum, and /one.pnum/eight.pnum./six.pnum/one.pnum% in Case /four.pnum. /one.pnum/nine.pnum (a) (c) (e) (g)(b) (d) (f) (h)Figure /one.pnum/three.pnum: Example /two.pnum - Case /four.pnum: Wave responses, u(x;y;t )[m], in the domain induced by ( left) the targeted set ofPkjand ( right ) its reconstructed counterpart at (a-b) /zero.pnum./three.pnum s, (c-d) /zero.pnum./four.pnum s, (e-f) /zero.pnum./five.pnum s, and (g-h) /zero.pnum./six.pnum s. /two.pnum/zero.pnum 02004006008001000120014001600180020000102030405060708090100Figure /one.pnum/four.pnum: Example /three.pnum - The terminal value of E, /one.pnum./four.pnum/five.pnum%, for reconstructing Pkjin Case /two.pnum is smaller than its counterparts, /two.pnum./eight.pnum/two.pnum% in Case /five.pnum and /one.pnum/zero.pnum./nine.pnum/zero.pnum% in Case /six.pnum. /six.pnum./five.pnum Example /three.pnum: Investigating the inversion performance with respect to the dominant frequency of the time signal of targeted traction This example studies the accuracy of the inversion with respect to the dominant frequency of the time signal of targeted traction. We used /five.pnum/zero.pnum sensors on the top surface t, /one.pnum/five.pnum sensors on the vertical array along r, and material pro/uniFB01le /one.pnum. In addition to the Case /two.pnum in Example /one.pnum, where the central frequency of targeted traction is /five.pnum Hz, we examined Case /five.pnum and /six.pnum, which consider targeted traction of central frequencies of /one.pnum/zero.pnum and /two.pnum/zero.pnum Hz, respectively. Fig. /one.pnum/four.pnum shows that, the higher frequency content the targeted traction contains, the higher terminal value ofEis obtained. Namely, the terminal value of E, /one.pnum./four.pnum/five.pnum%, for reconstructing ^F1in Case /two.pnum is smaller than its counterparts /two.pnum./eight.pnum/two.pnum% for identifying ^F2in Case /five.pnum and /one.pnum/zero.pnum./nine.pnum/zero.pnum% for reconstructing ^F3in Case /six.pnum. Thus, we suggest that the inversion modeling presented in this study su/uniFB00ers from solution multiplicity more when it identi/uniFB01es traction of a higher dominant frequency. Due to the traction time signals’ higher frequencies in Case /five.pnum and /six.pnum, more inversion iterations are used for reaching convergence of inversion solution in these cases than Case /one.pnum (i.e., /one.pnum/zero.pnum/zero.pnum/zero.pnum iterations in Case /one.pnum vs /two.pnum/zero.pnum/zero.pnum/zero.pnum iterations in Cases /five.pnum and /six.pnum). Fig. /one.pnum/five.pnum(a,c,e) show sets of Pkjcorresponding to targeted traction for Cases /two.pnum, /five.pnum, and /six.pnum, respectively, while Fig. /one.pnum/five.pnum(b,d,f) show their reconstructed counterparts. Figs. /one.pnum/six.pnum to /one.pnum/eight.pnum also show the snapshots of the wave re- sponses in the entire computational domain induced by the targeted ^F1;2;3and their reconstructed coun- terparts, respectively. All the /uniFB01gures present agreements between wave responses induced by the targeted and reconstructed force vectors. As shown in Table /one.pnum, Eu, the error between wave responses induced by ^Ftarget and its /uniFB01nal-estimated counterpart, is smaller than the terminal value of error Efor all Cases. For instance, in Case /six.pnum, the terminal value of Eis /one.pnum/zero.pnum./nine.pnum/zero.pnum%, butEuThus, we suggest that the larger terminal value ofEin the reconstructed force vector of a higher central frequency does not necessarily lead to Eu of the same order of magnitude due to inherent low-pass /uniFB01ltering of a typical FEM wave solver. /six.pnum./six.pnum Example /four.pnum: Examining the feasibility of the presented inverse modeling to recon- struct targeted traction of a realistic seismic signal In this example, we study the performance of identifying a targeted set of Pkjof^F4, which is created by a realistic seismic signal (Case /seven.pnum). In Case /seven.pnum, we use the material pro/uniFB01le /one.pnum with /five.pnum/zero.pnum sensors on tand /one.pnum/five.pnum sensors onr. Fig. /one.pnum/nine.pnum shows the excellent agreement between the targeted and reconstructed Pkjin Case /seven.pnum. The inversion performance in Case /seven.pnum is also compared with that in Case /two.pnum. Cases /two.pnum and /seven.pnum use the same material and sensor con/uniFB01gurations except di/uniFB00erent force vectors ^F1and^F4, respectively. Fig. /two.pnum/zero.pnum shows that the terminal value ofE, /zero.pnum./three.pnum/three.pnum%, in Case /seven.pnum is smaller than that in Case /two.pnum (/one.pnum./four.pnum/five.pnum%). We note that, even though the signal to generate ^F4is more complex than the other corresponding to ^F1, our inversion solver identi/uniFB01es ^F4more /two.pnum/one.pnum Figure /one.pnum/five.pnum: Example /three.pnum: (a) Target and (b) Reconstructed Pkjin Case /two.pnum; (c) Target and (d) Reconstructed Pkj in Case /five.pnum; and (e) Target and (f) Reconstructed Pkjin Case /six.pnum, for Case /two.pnum at the /one.pnum/zero.pnum/zero.pnum/zero.pnumth iteration, and Case /five.pnum-/six.pnum at the /two.pnum/zero.pnum/zero.pnum/zero.pnumth iteration. /two.pnum/two.pnum (a) (c) (e) (g)(b) (d) (f) (h)Figure /one.pnum/six.pnum: Example /three.pnum - Case /two.pnum: Wave responses, u(x;y;t )[m], in the domain induced by ( left) the targeted set ofPkjand ( right ) its reconstructed counterpart at (a-b) /zero.pnum./three.pnum s, (c-d) /zero.pnum./four.pnum s, (e-f) /zero.pnum./five.pnum s, and (g-h) /zero.pnum./six.pnum s. /two.pnum/three.pnum (a) (c) (e) (g)(b) (d) (f) (h)Figure /one.pnum/seven.pnum: Example /three.pnum - Case /five.pnum: Wave responses, u(x;y;t )[m], in the domain induced by ( left) the targeted set ofPkjand ( right ) its reconstructed counterpart at (a-b) /zero.pnum./three.pnum s, (c-d) /zero.pnum./four.pnum s, (e-f) /zero.pnum./five.pnum s, and (g-h) /zero.pnum./six.pnum s. /two.pnum/four.pnum (a) (c) (e) (g)(b) (d) (f) (h)Figure /one.pnum/eight.pnum: Example /three.pnum - Case /six.pnum: Wave responses, u(x;y;t )[m], in the domain induced by ( left) the targeted set ofPkjand ( right ) its reconstructed counterpart at (a-b) /zero.pnum./three.pnum s, (c-d) /zero.pnum./four.pnum s, (e-f) /zero.pnum./five.pnum s, and (g-h) /zero.pnum./six.pnum s. /two.pnum/five.pnum Figure /one.pnum/nine.pnum: Example /four.pnum: (a) A target set of Pkjand (b) its reconstructed counterpart corresponding to ^F4for Case /seven.pnum at the /one.pnum/zero.pnum/zero.pnum/zero.pnumth iteration. 010020030040050060070080090010000102030405060708090100 Figure /two.pnum/zero.pnum: Example /four.pnum - Eover iterations in Case /two.pnum and Case /seven.pnum. accurately than ^F1because the signal in ^F4has a lower dominant frequency than ^F1. This behavior is consistent with our /uniFB01nding from Example /three.pnum. Table /one.pnum shows that Eubetween the wave responses induced by ^F4and its /uniFB01nal-estimated counterpart is zero. In addition, Figs. /two.pnum/one.pnum and /two.pnum/two.pnum show the agreement between the wave responses induced by the target and those induced by its reconstructed counterpart in Case /seven.pnum. /six.pnum./seven.pnum Example /five.pnum: Investigating the inversion performance in a /two.pnum-layered background domain This example focuses on examining the performance of inverting a targeted set of Pkjof^F5by considering a two-layered background domain with inclusions. The performance is compared with that using the single- layered, homogeneous background domain counterpart. To this end, we used the material pro/uniFB01le /two.pnum, a vertical array of /one.pnum/five.pnum sensors on r, and /five.pnum/zero.pnum sensors on tin Case /eight.pnum. It is compared with Case /two.pnum, in which the material pro/uniFB01le /one.pnum is used with ^F1. The traction signals of both ^F5and^F1have the same central frequency of /five.pnum Hz. The sensor con/uniFB01guration in Case /eight.pnum is same as Case /two.pnum. Fig. /two.pnum/three.pnum shows the excellent agreement between the targeted Pkjand its reconstructed counterpart. Fig. /two.pnum/four.pnum shows the value of Eover iterations in Case /two.pnum and /eight.pnum for material pro/uniFB01les /one.pnum and /two.pnum, respectively. Despite using a more complex background domain in Case /eight.pnum, the di/uniFB00erence between the terminal values of Ein Cases /two.pnum /two.pnum/six.pnum (a) (c) (e) (g)(b) (d) (f) (h)Figure /two.pnum/one.pnum: Example /four.pnum - Case /seven.pnum: Wave responses, u(x;y;t )[m], in the domain induced by ( left) the targeted set ofPkjand ( right ) its reconstructed counterpart at (a-b) /one.pnum./zero.pnum s, (c-d) /two.pnum./five.pnum s, (e-f) /three.pnum./zero.pnum s, and (g-h) /three.pnum./five.pnum s. /two.pnum/seven.pnum (a) (c) (e) (g)(b) (d) (f) (h)Figure /two.pnum/two.pnum: Example /four.pnum - Case /seven.pnum: Wave responses, u(x;y;t )[m], in the domain induced by ( left) the targeted set ofPkjand ( right ) its reconstructed counterpart at (a-b) /three.pnum./seven.pnum/five.pnum s, (c-d) /four.pnum./five.pnum s, (e-f) /five.pnum./two.pnum/five.pnum s, and (g-h) /six.pnum./zero.pnum s. /two.pnum/eight.pnum (a)(b)Figure /two.pnum/three.pnum: Example /five.pnum: (a) A target set of Pkjand (b) its reconstructed counterpart corresponding to ^F5for Case /eight.pnum at the /one.pnum/zero.pnum/zero.pnum/zero.pnumth iteration. 010020030040050060070080090010000102030405060708090100 Figure /two.pnum/four.pnum: Example /five.pnum - The /uniFB01nal value of E, in Case /eight.pnum, is /one.pnum./five.pnum/six.pnum%, and that in Case /two.pnum is /one.pnum./four.pnum/five.pnum%. and /eight.pnum is only /zero.pnum./one.pnum/one.pnum%. Table /one.pnum shows thatEuis only /zero.pnum./one.pnum/six.pnum% in Case /eight.pnum, and it is as low as Case /two.pnum. Fig. /two.pnum/five.pnum presents the wave responses in the domain induced by the targeted ^F5and its /uniFB01nal-reconstructed counterpart in Case /eight.pnum at t= /zero.pnum./three.pnum s, /zero.pnum./four.pnum s, /zero.pnum./five.pnum s, and /zero.pnum./six.pnum s. In general, the responses are in excellent agreement with each other. Thus, we suggest that the accuracy of our presented inversion algorithm is not compromised by the material complexity of a background domain. /six.pnum./eight.pnum Example /six.pnum: Examining the e/uniFB00ect of noisy data on the inversion performance In this example, we compare the accuracy of the inversion with respect to the noise level of random noise that is added to measured signals um(x;y;t )at sensor locations prior to the inversion. We considered the material pro/uniFB01le /two.pnum with /five.pnum/zero.pnum sensors on tand /one.pnum/five.pnum sensors on r, and examined the inversion performance of identifying a targeted set of Pkjof^F5in Cases /eight.pnum-/one.pnum/two.pnum, which consider the noise level of /zero.pnum%, /one.pnum%, /two.pnum%, /five.pnum%, and /one.pnum/zero.pnum%, respectively. Fig. /two.pnum/six.pnum shows the targeted Pkjand its reconstructed counterparts obtained in Cases /eight.pnum-/one.pnum/two.pnum. Fig. /two.pnum/six.pnum shows a more signi/uniFB01cant di/uniFB00erence between the target and estimated Pkjfor a larger noise level. Fig. /two.pnum/seven.pnum reveals that the terminal value of Eincreases as the noise level increases. For instance, while the terminal value of Ein Case /eight.pnum, of a noise level /zero.pnum%, is /one.pnum./five.pnum/six.pnum%, the /uniFB01nal value of Ein Case /one.pnum/two.pnum, where /one.pnum/zero.pnum% of noise level is studied, is /one.pnum/one.pnum/two.pnum./four.pnum/zero.pnum%. Although the terminal value of Eexponentially grows as the noise level grows, Table /one.pnum shows that Euin- creases in a much smaller ratio: Euis less than /one.pnum% for all the Cases in this example. For instance, while the terminal value ofEin Case /one.pnum/two.pnum is /one.pnum/one.pnum/two.pnum./four.pnum/zero.pnum%, Euis only /zero.pnum./eight.pnum/six.pnum%. The snapshots of the wave responses in the entire computational domain induced by the targeted and estimated force vector in Case /one.pnum/two.pnum are shown in Fig. /two.pnum/eight.pnum, and they are in good agreement with each other. It is because the noise-like error in the recon- /two.pnum/nine.pnum (a) (c) (e) (g)(b) (d) (f) (h)Figure /two.pnum/five.pnum: Example /five.pnum - Case /eight.pnum: Wave responses, u(x;y;t )[m], in the domain induced by ( left) the targeted set ofPkjand ( right ) its reconstructed counterpart at (a-b) /zero.pnum./three.pnum s, (c-d) /zero.pnum./four.pnum s, (e-f) /zero.pnum./five.pnum s, and (g-h) /zero.pnum./six.pnum s. /three.pnum/zero.pnum (a) (c) (e)(b) (d) (f)Figure /two.pnum/six.pnum: Example /six.pnum: (a) A set of Pkjcorresponding to targeted traction and (b-f) reconstructed sets of Pkj in Cases /eight.pnum-/one.pnum/two.pnum, corresponding to noise level of /zero.pnum%, /one.pnum%, /two.pnum%, /five.pnum%, and /one.pnum/zero.pnum%, respectively, after /one.pnum/zero.pnum/zero.pnum/zero.pnum iterations. structedPkjdoes not a/uniFB00ect the wave responses due to the low-pass /uniFB01ltering of the FEM solver. Thus, we suggest that the terminal values of EandEuincrease as the noise level increases, but not in a similar ratio to each other. Furthermore, as discussed previously in Example /three.pnum, a large value of Edoes not necessarily mean a value of Euof the same order of magnitude. /seven.pnum Conclusion We present a robust mathematical and computational modeling that identi/uniFB01es complex, incoherent seis- mic incident wave/uniFB01eld in a /two.pnumD solid, truncated by a WABC. The PDE-constrained optimization method is utilized to tackle the inverse problem. The DTO approach is employed to solve the adjoint problem and to semi-analytically compute the gradient of the discrete objective functional in the inverse modeling proce- dure. Numerical results depict the performance of this inversion approach in a domain enclosed by WABC. First, targeted traction—mimicking incident, inclined plane waves—can not be fully reconstructed if we use only the sensors on the top surface. It is because wave motions induced by a part of dynamic traction at the bottom boundary do not reach the top surface of the considered truncated domain. Thus, to improve the inversion performance, a vertical array of sensors on the side boundary of a domain should be adopted. Second, to optimize the algorithm’s inversion e/uniFB03ciency, a su/uniFB03ciently large number of sensors should be used. In the presented example, a desirable ratio of measurement data size to the number of the control parameters can be as small as /one.pnum:/eight.pnum. Third, due to the intrinsic low-pass /uniFB01ltering of the FEM wave solver, a /three.pnum/one.pnum 01002003004005006007008009001000020406080100120Figure /two.pnum/seven.pnum: Example /six.pnum - The terminal value of Ecomparing Cases /eight.pnum-/one.pnum/two.pnum with respect to the noise level. (a) (c) (e) (g)(b) (d) (f) (h) Figure /two.pnum/eight.pnum: Example /six.pnum - Case /one.pnum/two.pnum: Wave responses, u(x;y;t )[m], in the domain induced by ( left) the targeted set ofPkjand ( right ) its reconstructed counterpart at (a-b) /zero.pnum./three.pnum s, (c-d) /zero.pnum./four.pnum s, (e-f) /zero.pnum./five.pnum s, and (g-h) /zero.pnum./six.pnum s. /three.pnum/two.pnum large value of the inversion error in the reconstructed traction does not necessarily lead to an error of the same order of magnitude in the corresponding reconstructed wave responses in the entire computational domain. Fourth, the minimizer su/uniFB00ers less from solution multiplicity when it identi/uniFB01es lower dominant frequency traction (e.g., a real earthquake signal). Fifth, the accuracy of our presented inversion algorithm is not compromised by the material complexity of a background domain. Lastly, the terminal values of E and errorEugrow as the noise level that is added to umis increased, but not in a similar ratio to each other. The wave responses in the entire computational domain, caused by the targeted traction and the restored one, are in excellent agreement. Therefore, by extending the presented approach to practical /three.pnumD elasto- dynamic settings, this algorithm can determine where signi/uniFB01cant amplitudes of stress waves (or structural failures) occur in built-in environments of interest during seismic events. The presented proposed method could also be extended in the /three.pnumD setting such that it could /uniFB02exibly ac- commodate measurement types (e.g., full translational /three.pnumC motions, /six.pnumC (translational /three.pnumC with rotational /three.pnumC motions), vertical-only /one.pnumC motions, etc.). Thus, our method can fully take advantage of big data of seismic measurement data from modern infrastructures and built environments equipped with ubiquitous sen- sors. Data Availability Some or all data, models, or code generated or used during the study are available from the corresponding author by request. •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. Acknowledgment This material is based upon work supported by the National Science Foundation under Award 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. /eight.pnum Nomenclature /three.pnum/three.pnum Symbol Comment x;y;z Horizontal, vertical, and anti-plane directions, respectively t Time T Final time of the total observation duration u(x;y;t ) Displacement /uniFB01eld in the anti-plane ( z) direction of the wave motion G(x;y) Shear modulus Vs(x;y) Shear wave speed (x;y) Mass density Domain J Total observation duration t;b;l;rTop, bottom, left, and right boundaries of the domain D;L y -coordinate of bandx-coordinate of r f(x;y;t ) Dynamic traction uIIncident seismic wave/uniFB01eld v Test function U;V Function spaces H1Hilbert space (x;y) Vector of global basis function v Vector of nodal solution of test function u(t) Solution vector of the state problem at t _u(t);u(t) First- and second-order derivatives of solution vector with respect to t M;K;C Global mass, sti/uniFB00ness, and dumping matrices F Global force vector j Thej-th time step; the j-th component in Pkj t Size of time step Q Discrete forward operator ^u Space-time discretization of u(x;y;t )for all the time steps ^F Global force vector corresponding to all the time steps N Total number of time steps k Thek-th discrete node on the absorbing boundary k Numbering of the discrete node k; the numbering of the k-th component in Pkj Pkj Component of ^Fcorresponding to kandtj ^L Discrete objective functional um(x;y;t ) Wave signal measured by a sensor at a location of xandyinduced by a targeted traction ^um Space-time discretization of um(x;y;t )for all the time steps B Square matrix that is zero everywhere except on the diagonals that correspond to a degree of freedom B Block diagonal matrix ^A Discrete Lagrangian functional ^ Lagrange multiplier vector for all the time steps  A set of control parameters, =Pkj g Gradient vector d Search direction vector h Scalar-value step size s Thes-th inversion iteration E Error norm between estimated traction and its target EuError norm between wave responses induced by targeted traction and its counterpart /three.pnum/four.pnum Sensor NodesWave-absorbing boundaryFigure /two.pnum/nine.pnum: Problem setting of Appendix A A On proof about the objective functional being quadratic and convex in terms of control parameters In this section, we prove that the objective functional is quadratic and convex in terms of control parame- ters, which are housed in ^F. From the compact form of the state problem, ^ucan be de/uniFB01ned as: ^u=Q1^F: (/four.pnum/seven.pnum) Due to (/four.pnum/seven.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; (/four.pnum/eight.pnum) where the part ais a constant, and the parts bandcare linear functions, in terms of ^F, so that their Hessians with respect to the control parameters vanish. The part din (/four.pnum/eight.pnum) is a quadratic function in terms of^F. Therefore, proving that the part din (/four.pnum/eight.pnum) is convex will show that ^Lis convex. To this end, we consider the following simple example—a homogeneous square-shaped solid domain, of which extent is /nine.pnum/zero.pnum m /six.pnum/zero.pnum m with its shear wave speed of /three.pnum/five.pnum/three.pnum./five.pnum/five.pnum m/s and mass density of /three.pnum/zero.pnum/zero.pnum/zero.pnum kg/m3. The shear-traction nodal forces are applied on one node in the bottom surface ( x= /three.pnum/zero.pnum m and y= -/six.pnum/zero.pnum m) and one node in the left boundary ( x= /zero.pnum m andy= /three.pnum/zero.pnum m), and one sensor is placed at x= 30 m andy= /zero.pnum m on the top surface (see Fig. /two.pnum/nine.pnum). The solid is truncated by wave-absorbing boundary conditions on the left, bottom, and right boundaries, and the element size is /three.pnum/zero.pnum m. Therefore, in this case, there are twelve degrees of freedom in the discretized domain. We are using only two shear-traction nodal forces in this simple example, but we can extend it to a case with more shear-traction nodal forces (e.g., six dynamic traction nodal forces—three in the left and three in the bottom surface). We consider /two.pnum time steps, where t= 0:1s. Then, the matrices M,C, andKare calculated by using (/one.pnum/two.pnum), and Qis obtained as the following: Q=2 6666664I 0 0 0 0 0 0 I 0 0 0 0 K C M 0 0 0 L1L2L3Ke 0 0 a1I I 0a1I I 0 a0Ia2I Ia0I0I3 7777775; (/four.pnum/nine.pnum) /three.pnum/five.pnum where: a0=4 (t)2; a1=2 t; a2=4 t; Ke =a0M+a1C+K;L1=a0Ma1C;L2=a2MC;L3=M: (/five.pnum/zero.pnum) The vector ^Fis built as: ^F=2 66666640 0 F0 F1 0 03 7777775;F0= [0P4;00 0P2;00 0 0 0 0 0 0]T;F1= [0P4;10 0P2;10 0 0 0 0 0 0]T; (/five.pnum/one.pnum) whereP2;0andP2;1are the components of a force vector at the node of k= 2and at the timesteps of t=/zero.pnum s andt=/zero.pnum./one.pnum s, respectively. Similarly, P4;0andP4;1are the components of a force vector at the node of k= 4and at the timesteps of t=/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; (/five.pnum/two.pnum) where B1andB2are a /one.pnum/two.pnumx/one.pnum/two.pnum square matrix that is zero everywhere except for the /one.pnum/zero.pnum-th diagonal compo- nent with its value of /one.pnum that corresponds to the degree of freedom where the sensor is placed. Based on (/four.pnum/nine.pnum) to (/five.pnum/two.pnum), the part din (/four.pnum/eight.pnum) can be written as: d= 1:471021P2 2;0+ 1:311021P2;0P4;0+ 2:941021P2;0P2;1+ 1:311021P2;0P4;1 + 2:921022P2 4;0+ 1:311021P4;0P2;1+ 5:841022P4;0P4;1 + 1:471021P2 2;1+ 1:311021P2;1P4;1 + 2:921022P2 4;1: (/five.pnum/three.pnum) As shown in (/five.pnum/three.pnum), the part din (/four.pnum/eight.pnum) is a quadratic function in terms of P2;0,P4;0,P2;1, andP4;1. We note that the part din (/four.pnum/eight.pnum) is convex if the eigenvalues of its Hessian matrix are greater than or equal to /zero.pnum. Here, the Hessian matrix of the part din (/four.pnum/eight.pnum) is: H=2 666664@2d @P2 2;0@2d @P2;0@P4;0@2d @P2;0@P2;1@2d @P2;0@P4;1 @2d @P4;0@P2;0@2d @P2 4;0@2d @P4;0@P2;1@2d @P4;0@P4;1 @2d @P2;1@P2;0@2d @P2;1@P4;0@2d @P2 2;1@2d @P2;1@P4;1 @2d @P4;1@P2;0@2d @P4;1@P4;0@2d @P4;1@P2;1@2d @P2 4;13 777775 =2 6642:9410211:3110212:9410211:311021 1:3110215:8410221:3110215:841022 2:9410211:3110212:9410211:311021 1:3110215:8410221:3110215:8410223 775; (/five.pnum/four.pnum) of which eigenvalues are obtained as: |I*- H|= 2:9410211:3110212:9410211:311021 1:3110215:8410221:3110215:841022 2:9410211:3110212:9410211:311021 1:3110215:8410221:3110215:841022 = 0:(/five.pnum/five.pnum) Namely, the eigenvalues are: 1;2= 0;  3= 7:051021;  4= 4:881025: (/five.pnum/six.pnum) The eigenvalues are non-negative so that the part din (/four.pnum/eight.pnum) is convex, and, therefore, the objective func- tional ^Lis convex. /three.pnum/six.pnum B Derivation of the /uniFB01rst term of the adjoint problem In this section, we show the derivation of the part ain (/three.pnum/two.pnum). a=@^L @^u=@1 2(^u^um)TB(^u^um) @^u(/five.pnum/seven.pnum) To this end, we compute the derivative using the product rule: a=1 2@(^u^um)T @^uB(^u^um) +1 2(^u^um)TB@(^u^um) @^u: (/five.pnum/eight.pnum) Next, we use the following rule for vector di/uniFB00erentiation: @(rTw) @r=@(wTr) @r=w; (/five.pnum/nine.pnum) where wis treated as constant in the derivative. Therefore, by denoting w1=B(^u^um)and(w2)T= (^u^um)TB, and introducing them in (/five.pnum/eight.pnum), part acan be written as: a=1 2@(^u^um)T @^uw1+1 2wT 2@(^u^um) @^u =1 2@f(^u^um)Tw1g @^u+1 2@fwT 2(^u^um)g @^u =1 2w1+1 2w2 =1 2B(^u^um) +1 2BT(^u^um) =1 2(B+BT) (^u^um): (/six.pnum/zero.pnum) The block diagonal matrix Bis symmetric, then: a=B(^u^um) (/six.pnum/one.pnum) /three.pnum/seven.pnum