მთავარი
Computational Geosciences On an inverse source problem for enhanced oil recovery by wave motion maximization in reservoirs
On an inverse source problem for enhanced oil recovery by wave motion maximization in reservoirs
Karve, Pranav M., Kucukcoban, Sezgin, Kallivokas, Loukas F.როგორ მოგეწონათ ეს წიგნი?
როგორი ხარისხისაა ეს ფაილი?
ჩატვირთეთ, ხარისხის შესაფასებლად
როგორი ხარისხისაა ჩატვირთული ფაილი?
ტომი:
19
ენა:
english
ჟურნალი:
Computational Geosciences
DOI:
10.1007/s1059601494627
Date:
February, 2015
ფაილი:
PDF, 4,39 MB
თქვენი თეგები:
 გთხოვთ, პირველ რიგში შეხვიდეთ ანგარიშზე

დახმარება გჭირდებათ? გთხოვთ, გაეცნოთ ინსტრუქციას როგორ გადააგზავნოთ წიგნი Kindleზე
15 წუთის განმავლობაში ფაილი გადაიგზავნება თქვენს emailზე.
15 წუთის განმავლობაში ფაილი გადაიგზავნება თქვენს kindleზე.
შენიშვნა: აუცილებელია თითოეული წიგნის ვერიფიცირება, რომელსაც აგზავნით Kindleზე. შეამოწმეთ თქვენი საფოსტო ყუთი, მოსულია თუ არა Amazon Kindle Support წერილი დადასტურებით.
შენიშვნა: აუცილებელია თითოეული წიგნის ვერიფიცირება, რომელსაც აგზავნით Kindleზე. შეამოწმეთ თქვენი საფოსტო ყუთი, მოსულია თუ არა Amazon Kindle Support წერილი დადასტურებით.
ურთიერთდაკავშირებულ წიგნთა სია
0 comments
შეგიძლიათ დატოვოთ გამოხმაურება წიგნის შესახებ და გააზიაროთ თქვენი გამოცდილება. სხვა მკითხველისთვის საინტერესო იქნება თქვენი მოსაზრება წაკითხული წიგნების შესახებ. მიუხედავად იმისა მოგწონთ თუ არა წიგნი, მასზე გულწრფელი და დეტალური მსჯელობა, ადამიანებს მისცემს საშუალებას იპოვონ ახალი წიგნები, რომლებიც მათ დააინტერესებთ.
1


2


Comput Geosci DOI 10.1007/s1059601494627 ORIGINAL PAPER On an inverse source problem for enhanced oil recovery by wave motion maximization in reservoirs Pranav M. Karve · Sezgin Kucukcoban · Loukas F. Kallivokas Received: 10 July 2014 / Accepted: 11 December 2014 © Springer International Publishing Switzerland 2014 Abstract We discuss an optimization methodology for focusing wave energy to subterranean formations using strong motion actuators placed on the ground surface. The motivation stems from the desire to increase the mobility of otherwise entrapped oil. The goal is to arrive at the spatial and temporal description of surface sources that are capable of maximizing mobility in the target reservoir. The focusing problem is posed as an inverse source problem. The underlying wave propagation problems are abstracted in two spatial dimensions, and the semiinfinite extent of the physical domain is negotiated by a buffer of perfectlymatchedlayers (PMLs) placed at the domain’s truncation boundary. We discuss two possible numerical implementations: Their utility for deciding the tempospatial characteristics of optimal wave sources is shown via numerical experiments. Overall, the simulations demonstrate the inverse source method’s ability to simultaneously optimize load locations and time signals leading to the maximization of energy delivery to a target formation. Keywords Inverse source problem · Enhanced oil recovery · PDEconstrained optimization · Elastic wave energy focusing P. M. Karve · L. F. Kallivokas () Department of Civil, Architectural and Environmental Engineering, The University of Texas at Austin, 301 East Dean Keeton St. Stop C1747, Austin, TX 78712, USA email: loukas@mail.utexas.edu S. Kucukcoban Stress Engineering Services Inc., 13610 Westland East Blvd. Building 2, Houston, TX 77041, USA 1 Introduction Elastic wave stimulation of oil reservoirs has been suggested as a viable recourse for the purpose of enhanced oil recovery (EOR) [3, 4, 6, 10, 12, 16]. Seismic EOR, as it is somet; imes referred to, relies on the ability of propagating waves, generated by sources placed on or below the ground surface, to mobilize trapped oil particles. Laboratory and field experiments [4] suggest that such mobilization is possible either by dislodgement into the fluid flow of the oil particles adhering to the pore wall [12] or by the release of trapped oilblobs from the pores due to a Haines jumplike phenomenon [2, 3, 6]. After the droplets are dislodged, conventional methods can be used for recovery. In general, wavebased EOR methods are economically competitive and do not suffer from low sweep efficiency problems in heterogeneous reservoirs when compared to other EOR methods [10, 14]. From a technical perspective, uniform, or as uniform as possible, illumination of the reservoir by a sufficiently strong wave field is essential to the success of wavebased EOR methods. Radiation damping, intrinsic, and apparent attenuation present challenges in the field implementation of any wavebased EOR process. Although the dissipative character of the problem cannot be altered, intelligent choices for the locations and frequency content of the wave sources can help maximize reservoir shaking. For example, wave sources operating at certain amplification frequencies (characteristic of the reservoir and its surroundings) have been shown to produce stronger wave motion in the targeted formations [10] than blind sources operating at arbitrarily selected frequencies. It can also be argued that if the sources are situated and synchronized in a manner that promotes Comput Geosci constructive interference at the target formation, then sufficiently strong motion may result despite the attenuation. In this article, we discuss a systematic methodology for making a judicious choice about both the placement and the frequency content of wave sources used for focusing wave energy to subsurface formations. Frequency sweep is one possibility for determining the optimal (monochromatic) frequency of a ground motion actuator at which the oil mobility in the formation will be maximized. When given the location and strength of the wave sources, the frequency sweep relies on a linear elastodynamic or poroelastodynamic model for capturing the relevant physics and for determining the strength of the wave field within the reservoir over the driving frequency range of interest. Although traversing the frequency spectrum with fixed load locations is a computationally manageable task, an exhaustive search for the optimal load locations that could potentially maximize mobility is computationally prohibitive. This shortcoming of the frequency sweep approach could be overcome if the problem were formulated in terms of a search for the optimal source characteristics, both spatial and temporal, which, in turn, calls for casting the search as an inverse source problem. The inverse source problem approach for the maximization of motion in a target formation is akin to the inverse medium problem typically arising in exploration geophysics [1, 5, 7, 11, 18], where subsurface imaging relies on the minimization of the difference between experimentally collected and computationally obtained data. Similarly, the inverse source problem is formulated as a constrained Fig. 1 General framework of the inverse source process optimization problem, where minimization of an appropriate objective functional is tantamount to the maximization of wave motion in the reservoir. The search for optimal source characteristics is carried out by staying within the constraints of the governing physics and of equipment limitations. Schematically, the inversion procedure is depicted in Fig. 1. The procedure is initialized with the definition of the motion metric in the objective functional (L) and with initial guesses for the source characteristics, i.e., for the source location and the associated time signals (control variables). Following a constrained optimization approach, whereby the governing IBVP is sideimposed to L, the firstorder optimality conditions typically lead to a triad of forward, adjoint, and control problems (box A, Fig. 1). The wave source characteristics are updated, following, for example, a gradientbased scheme, which requires computation of the functional’s gradients in the space of the source characteristics. Thus, in general, the numerical solution of the inverse source problem performs two key operations: optimization (minimization) in the space of the control variables and discretization of the ensuing triad of BVPs. These operations can be performed in different order and would yield different gradients of the objective functional in question. Depending on the operational order, there arises either an optimizethendiscretize (OTD) or a discretizethenoptimize (DTO) approach; they are both shown in Fig. 2 [8, 15]. Either of the two approaches can be used to solve an inverse source problem aimed at the maximization of wave motion within a reservoir. For example, in Jeong et al. [9, Comput Geosci Fig. 2 Candidate algorithms for the core optimization scheme depicted in box A of Fig. 1 10], the OTD approach was used for optimizing the frequency content of surface sources used in wavebased EOR for reservoirs abstracted in one or two spatial dimensions, where, however, the load locations were treated as fixed. In this article, we adopt a similar twodimensional elastic setting for optimizing the wave source characteristics, but seek to optimize both the source frequency content (or time signal) and the source locations. The reservoir is represented by a target inclusion within an arbitrarily heterogeneous geological formation. We utilize hybrid perfectlymatchedlayers (PMLs) [13] to truncate the physical domain of interest. We discuss the formulation and implementation of both the OTD and the DTO approaches and conduct numerical experiments. We compare the performance of the frequency sweep approach with that of the inverse source problem approach and show that simultaneous optimization with respect to both source locations and time signals is important for a higher energy yield. We also discuss the effect of polarization of the applied surface tractions on the efficiency of energy focusing and report results attesting to the ability of the inverse source approach to prescribe wave sources that could improve the mobility of a reservoir’s bypassed oil. forward problem consists of finding the displacement field (velocity and acceleration) in the domain of interest, given suitable boundary and initial conditions, which include the applied surface tractions. We discuss next the elements of the forward problem. 2.1 Strong form We are concerned with elastic wave propagation in a twodimensional (2D) heterogeneous, elastic halfspace. We intend to use finite elements for the numerical solution of 2 The forward problem Our working hypothesis is that the material properties of the targeted geostructure and those of the surrounding formations are known a priori. For example, if a in Fig. 3 is the target reservoir, then the properties of both a and its surrounding formations reg are considered known. The Fig. 3 Problem definition Comput Geosci the BVPs. This choice, in turn, necessitates the truncation of the semiinfinite domain. We use hybrid perfectlymatchedlayers (PMLs)[13] to realize a physically faithful simulation of wave propagation within the computational domain. Figure 3 shows the truncated 2D elastic domain reg that has been augmented with PMLs (PML ). reg envelops the target inclusion a . ua (x, t) represents the displacement field in the target inclusion (a ), whereas ub (x, t) denotes the displacement in the rest of the domain (reg ∪ PML ). The governing equations in = a ∪reg ∪PML , for time t ∈ (0, T ] = J, are given as : div μa ∇ua + ∇uTa + div [{λa div ua } I] − ρa üa = 0, x ∈ a , (1) and interface conditions : = u− , x ∈ I , u+ b b μb ∇ub + ∇uTb + {λb div ub } I n+ I = ˜ e + ST ˜ p n− , x ∈ I , − ṠT I (5b) ua (x, t) = ub (x, t), (5c) σ a (x, t)T n− a = x ∈ a , −σ b (x, t)T n+ a , x ∈ a ; (5a) (5d) where σ a (x, t) = μa ∇ua + ∇uTa + λa (div ua )I, σ b (x, t) = μb ∇ub + ∇uTb + λb (div ub )I; (5e) (5f) div μb ∇ub + ∇uTb + div [{λb div ub } I] −ρb üb = 0, x ∈ reg , (2a) ˜ e + ST ˜ p − ρb (aüb + bu̇b + cub ) = 0, div ṠT x ∈ PML , 1 ˜p + ˜ p (∇ub )T D : (aS̈ + bṠ + cṠ) − ∇ub 2 ˜e + ˜ e (∇ u̇b )T = 0, x ∈ PML , +∇ u̇b (2b) (2c) where an overdot ( ˙ ) denotes a derivative with respect to time and a colon (:) represents tensor inner product. Equations 1 and 2 are the elastodynamics equations for a semiinfinite medium, augmented by PMLs [13]. S(x, t) is the stress history tensor, given by, t S11 (x, t) S12 (x, t) S(x, t) = = σ (x, τ )dτ, (3) S21 (x, t) S22 (x, t) 0 ˜ e and ˜ p are comwhere σ is the Cauchy stress tensor. ponents of the stretching tensor, and a, b, c are coefficients defining coordinate stretching in the PML region. Their detailed definitions are beyond the scope of this article and can be found in [13]. (λa , μa , ρa ) and (λb , μb , ρb ) are the Lamé parameters and mass density for the inclusion a and reg ∪ PML , respectively. D is the compliance tensor, so that the constitutive law takes the form D : Ṡ = 12 (∇ub + ∇uTb ). For t ∈ J, the governing equations are subjected to the following boundary conditions: PML ub (x, t) = 0, x ∈ fixed , μb ∇ub + ∇uTb + {λb div ub } I n = f(x, t), (4a) x ∈ load , (4b) + {λb div ub } I n = 0, x ∈ free , (4c) μb ∇ub +∇uTb PML ˜ e + ST ˜ p n = 0, x ∈ free ṠT ; and initial conditions : (4d) ua (x, 0) = 0, u̇a (x, 0) = 0, x ∈ a , ub (x, 0) = 0, u̇b (x, 0) = 0, x ∈ reg ∪ PML , (6b) S(x, 0) = 0, Ṡ(x, 0) = 0, x ∈ PML . (6a) (6c) The strong form of the forward problem can be stated as: given loads f(x, t), find ua ∈ H1 (a ) × J, ub ∈ H1 (reg ∪ PML ) × J, and S ∈ L2 (PML ) × J, so that they satisfy Eqs. 1 and 2 and conditions (4)–(6), where the pertinent function spaces for a scalar f , vector u, and tensor T are given by: L2 () = {f : f 2 d < ∞}, (7a) L2 () = {T : T ∈ (L2 ())2×2 }, H 1 () = {f : (f 2 + ∇f 2 )d < ∞, (7b) PML f (x) = 0 if x ∈ fixed }, 1 1 2 H () = {u : u ∈ (H ()) }. (7c) (7d) 2.2 Weak form Next, we cast the forward problem in its weak form in the Galerkin sense. We multiply (1) by a vector test function va , multiply (2a), (2b) by a vector test function vb , integrate by parts on their corresponding domains, and add them. Similarly, we multiply (2c) with tensor test function T and Comput Geosci integrate over PML . After some simplifications, using the boundary conditions (4) and interface conditions (5), we arrive at the following weak form: ∇va : μa (∇ua + ∇uTa ) + λa (div ua )I a 2.4 Time integration +va · ρa üa } d ∇vb : μb (∇ub + ∇uTb ) + λb (div ub )I + The time line is now discretized using a timestep t. We define the vector di = d, at time t = i t. The equation of motion of the spatiotemporally discretized system at time t = (i + 1) t can be written as, reg +vb · ρb üb } d ˜ e + ST ˜p ∇vb : ṠT + PML Md̈i+1 + Cḋi+1 + Kdi+1 = Fi+1 . +vb · ρb (aüb + bu̇b + cub )} d = vb load · f d, 1 ˜p T : D : (aS̈ + bṠ + cṠ) − ∇ub 2 PML ˜ p (∇ub )T + ∇ u̇b ˜e + ˜ e (∇ u̇b )T d = 0. + (8a) (8b) Numerical solution of the forward problem requires discretization in space and time. We introduce spatial approximations via shape functions (x) ∈ H1h () and (x) ∈ L2h (PML ). Thus, the trial and test functions can be expressed as, T va1 (x) (x)T ua1 (t) va = , u = , a vTa2 (x) (x)T ua2 (t) vTb1 (x) (x)T ub1 (t) vb = , u = , (9a) b (x)T ub2 (t) vTb2 (x) T = S = TT11 (x) TT12 (x) , TT21 (x) TT22 (x) (x)T S11 (t) (x)T S21 (t) (x)T S12 (t) (x)T S22 (t) (14) where d0 and ḋ0 are the prescribed displacement and velocity vectors. For any time t ∈ J, i.e., for i ≥ 0, we calculate the acceleration (d̈i+1 ) by solving, Meff d̈i+1 = Reff , where, Meff = M + Cγ ( t) + K( t)2 β = M + b4 C + b2 K, Reff = Fi+1 − Kdi − [C + ( t)K] ḋi − (1 − γ )( t)C + (0.5 − β)( t)2 K d̈i = Fi+1 − L0 di − L1 ḋi − L2 d̈i . (15) (16) (17) We, then, calculate the displacement and velocity at the (i + 1)th timestep using, di+1 = di + ( t)ḋi + (0.5 − β)( t)2 d̈i + β( t)2 d̈i+1 = di + b0 ḋi + b1 d̈i + b2 d̈i+1 , (18) ḋi+1 = ḋi + (1 − γ )( t)d̈i + γ ( t)d̈i+1 , (9b) where, henceforth, quantities with a tilde over the quantity symbol () denote vectors of nodal values of the subtended quantity. Introducing Eq. 9 into Eqs. 8a and 8b yields the following semidiscrete equation: Md̈ + Cḋ + Kd = F, (13) We employ Newmark’s time integration scheme to integrate Eq. 13 in time. At time t = 0: Md̈0 = F0 − Cḋ0 − Kd0 , 2.3 Semidiscrete form of unknown displacements (everywhere) and stress histories (PML only), and F is the force vector. The definitions of the element matrices used in our formulation are given in Appendix A. Further details of the global and element matrices can be found in [13]. (10) where, T reg reg d = ua1 ua2  ub1 ub2  uPML uPML S11 S22 S12 ,(11) b1 b2 T F = 0 0  F1 F2  0 0 0 0 0 . (12) We note that M, C, and K are the global mass, damping, and stiffness matrices, respectively, d is the vector = ḋi + b3 d̈i + b4 d̈i+1 . (19) Solution of the forward problem can be obtained by starting with i = 0 and marching in time (i ← i +1), using Eqs. 14– 19. In Eqs. 15 to 19, we have introduced constants b0 − b4 and matrices L0 − L2 . Their definitions are given below: b0 = t, b1 = ( 12 − β)( t)2 , b2 = β( t)2 , (20a) b3 = (1 − γ )( t), b4 = γ ( t), (20b) L0 = K, L1 = C + b0 K, L2 = b3 C + b1 K. (20c) 3 Load modeling In the inverse source formulation, the surface loads are treated as unknowns. The load descriptors are updated iteratively during the inversion process (e.g., ξ in Fig. 1). Comput Geosci This calls for parameterization of the spatiotemporal characteristics of the loads. Specifically, the tractions f(x, t) applied on load consist of contributions fi (x, t) from ns sources. The ith source consists of a spatial θi (x) and a temporal fi (t) component. θi is further decomposed into the x1 directional component θi1(x) and the x2 directional component θi2 (x). Thus, f(x, t) = ns fi (x, t) = i=1 ns θi1 (x) i=1 θi2 (x) fi (t). (21) In our numerical experiments, we use loads in either the x1 or x2 direction, and, thus, consistently with equipment limitations that allow load application in one direction only, either θi1 (x) = 0 or θi2 (x) = 0. We parameterize the (unknown) time signal using piecewise quadratic Lagrange polynomials τj (t) whose temporal nodal values are denoted by ξij . This allows us to express fi (t) as fi (t) = nf (22) ξij τj (t), j =1 where nf is the total number of Lagrange polynomials. The spatial variation of the ith load on load is captured by θiK , where subscript K describes the direction in which the load acts. For example, a constant pressure load applied vertically on part of the surface (x2 = 0) can be expressed as θi2 (x1 , 0) = H ηi − b2i − H ηi + b2i , (23) where H is the Heaviside function, ηi is the x1 coordinate of the load’s centerline, and bi is the ith load width. A vertical load varying like a Gaussian function in space about ηi is given by θi2 (x1 , 0) = −exp −(x1 −ηi )2 bi 4.1 Optimizethendiscretize (OTD) approach We use a standard PDEconstrainedoptimization approach to resolve the inverse source problem. The objective functional is cast in its continuous form. The constraints imposed by the governing PDEs and Neumanntype boundary conditions are incorporated by augmenting the objective functional with the PDEs and boundary conditions, multiplied by continuous Lagrange variables. Upon discretization of the firstorder optimality conditions, the resulting KKT (Karush–Kuhn–Tucker) system is solved numerically using a reducedspace approach. Maximization of the kinetic energy within the elastic inclusion (a ) is tantamount to minimization of the following objective functional: L = (24) subject to the governing Eqs. 1 and 2, boundary, interface, and initial conditions (4)–(6). We side impose the PDEs (1) to (2c), and the boundary conditions (4b), (4d) on L to form the Lagrangian A. Other boundary, initial, and interface conditions will be explicitly imposed. There results: A= 1 a T ρa [u̇a · u̇a ] dt d 0 a λua · μa div ∇ua + ∇uTa T + 0 +div [{λa div ua } I] − ρa üa ] dt d λub1 · μb div ∇ub + ∇uTb T + 0 reg + +div [{λb div ub } I] − ρb üb ] T ˜ e + ST ˜p λub2 · div ṠT 0 PML + −ρb (aüb + bu̇b + cub )] dt d T λS : D : (aS̈ + bṠ + cṠ) 0 PML 1 ˜p + ˜ p (∇ub )T ∇ub 2 ˜e + ˜ e (∇ u̇b )T dt d +∇ u̇b − Next, we discuss the formulation and solution of the inverse problem. The inverse source problem aims at maximizing a motion metric (e.g., the kinetic energy) in the target domain by seeking optimal time signals and locations for the surface tractions. In this section, we discuss the formulations for both the OTD and DTO approaches. (25) ρa [u̇a · u̇a ] dt d In our numerical experiments, we use (approximately) 2m wide loads. Thus, we set bi = 2 m in Eq. 23 and bi = 1.25 m in Eq. 24. 4 The inverse source problem , 0 a . 1 T T + load + PML free 0 λF · μb ∇ub + ∇uTb + {λb div ub } I] n − f(x, t)] dt d T ˜ e + ST ˜ p n dt d λb · ṠT (26) 0 Notice that λ = {λua , λub , λS , λF , λb } are the continuous Lagrange multipliers. We, now, seek to satisfy the firstorder optimality conditions. Comput Geosci 4.1.1 State problem interface conditions : The variation of the Lagrangian with respect to the Lagrange variables results in: δλua A = 0, i = 1, 2 i δλub A = 0, i = 1, 2 i (27) δλSij A = 0, ij = 11, 12, 21, 22 δλFi A = 0, i = 1, 2 δλbi A = 0, i = 1, 2. Equations 27 are satisfied when the (continuous) forward problem, i.e., Eqs. 1–6, is satisfied. The discrete form of the forward problem has been addressed in Section 2. (32b) λua (x, t) = λub (x, t), x ∈ a , (32c) λσa (x, t)T n− a = −λσb (x, t)T n+ a , x ∈ a ; (32a) (32d) where, λσa (x, t) = μa (∇λua + ∇λTua ) + λa (div λua )I, (32e) λσb (x, t) = μb (∇λub + ∇λTub ) + λb (div λub )I; (32f) and final value conditions : 4.1.2 Adjoint problem  strong form Variation of A with respect to the state variables yields: δuai A = 0, δubi A = 0, δSij A = 0, − λ+ u = λub , x ∈ I , b μb (∇λub + ∇λTub ) + λb (div λub )I n+ I = ˜ e + λS ˜ p n− , x ∈ I , − −λ̇S I i = 1, 2 i = 1, 2 ij = 11, 12, 21, 22. (28) Equations 28 give rise to the strong form of the adjoint problem, for time t ∈ (T , 0] = Ja , div μa ∇λua + ∇λTua + div λa div λua I −ρa λ̈ua = f ua , x ∈ a , (29) λua (x, T ) = 0, λ̇ua (x, T ) = −E u̇a , x ∈ a , (33a) λub (x, T ) = 0, λ̇ub (x, T ) = 0, x ∈ reg , (33b) λS (x, T ) = 0, λ̇S (x, T ) = 0, x ∈ PML . (33c) As it can be seen from the above, the adjoint problem is driven by the body forces that, in turn, depend on the value of the energy metric the forward problem computes over the target formation, while the time line is reversed. 4.1.3 Adjoint problem  weak form and div μb ∇λub + ∇λTub + div λb div λub I −ρb λ̈ub = 0, x ∈ reg , (30a) ˜ e +λS ˜ p −ρb aλ̈ub − bλ̇ub + cλub = 0, div −λ̇S x ∈ PML , (30b) T T ˜ e (∇ λ̇b ) − ˜ p (∇λb ) = 0, D : aλ̈S − bλ̇S + cλS + x ∈ PML , (30c) where f ua = E ρa üa , E = T a 0 λS = −2 ρa [u̇a · u̇a ] dt d 2 , sym λS , subject to the following boundary conditions: PML λub (x, t) = 0, x ∈ fixed , T μb (∇λub + ∇λub ) + λb (div λub )I n = 0, x ∈ load , μb (∇λub + ∇λTub ) + λb (div λub )I n = 0, x ∈ free , PML ˜ e + λS ˜ p n = 0, x ∈ free −λ̇S , PML ˜ e + λS ˜ p n = 0, x ∈ fixed −λ̇S , λF (x, t) = −λub (x, t), x ∈ load , λb (x, t) = −λub (x, t), x ∈ PML free ; (31a) (31b) (31c) We take an inner product of adjoint equations with appropriate test functions and integrate over the corresponding domains to obtain the following weak form of the adjoint problem: ∇va : μa ∇λua + ∇λTua + λa div λua I a +va · ρa λ̈ua + ρf λ̈w d + ∇vb : μb ∇λub + ∇λTub + λb div λub I reg +vb · ρb λ̈ub d ˜ e + λS ˜ p + vb · ρb aλ̈ub + ∇vb : −λ̇S PML −bλ̇ub + cλub d = − va · f ua d. (34a) a ˜ e (∇ λ̇b )T T : D : aλ̈S − bλ̇S + cλS + PML ˜ p (∇λb )T d = 0. − (34b) (31d) (31e) (31f) (31g) 4.1.4 Adjoint problem  semidiscrete form We introduce spatial approximations via shape functions (x) ∈ H1h () and (x) ∈ L2h (PML ). Thus, the trial and Comput Geosci test functions are given by T va1 (x) (x)T λua1 (t) va = , λua = , vTa2 (x) (x)T λua2 (t) vTb1 (x) (x)T λub1 (t) , λub = vb = , vTb2 (x) (x)T λub2 (t) T = λS = TT11 (x) TT12 (x) , TT21 (x) TT22 (x) (x)T λ (35a) (x)T λ S11 (t) (x)T λS21 (t) 4.1.6 Control problem(s) S12 (t) . (x)T λS22 (t) (35b) We define reg reg 1 2 PML PML λ = [λua1 λua2  λub λub  λub λub λS11 λS22 λS12 ]T , (36) 1 2 Fadj = [Fua1 Fua2  0 0  0 0 0 0 0]T , E ρa T u¨ ai (t) d. Fuai (t) = − (37) (38) Formulation of the control problem depends on the type of optimization to be performed. The inverse problem can be cast to obtain either the optimal time signal or the optimal position of the surface tractions, or both of the aforementioned load descriptors. Time signal optimization: In this case, the goal of the inverse problem is to find the optimal set of parameters ξij , as defined in Eq. 22. The control variable is given by: ξ = [ξ11 ξ12 · · · ξ1nf · · · ξ(ns )(nf −1) ξns nf ]. For a load applied in the xk direction, the gradient of the Lagrangian with respect to the control parameter ξmn is given by a Now, Eqs. 34a and 34b can be written in their semidiscrete form, using (35), as Madj λ̈ + Cadj λ̇ + Kadj λ = Fadj . ∇(ξmn ) A = load θmk (x)T T 0 λub k (t)τn (t)dt Madj = M, Cadj = −C, Kadj = K. (40) 4.1.5 Adjoint problem  time integration The adjoint problem is a final BVP, and it requires traversing the time line backwards. To this end, we introduce the following approximations: λi = λi+1 − λ̇i+1 t + ( t)2 λ̈i (0.5 − β) + λ̈i+1 β , λ̇i+1 = λ̇i+1 − t λ̈i (1 − γ ) + λ̈i+1 γ ( t), (41) (42) where the subscripts denote timestep. Equations 41 and 42 are used in the following system of equations: adj Madj λ̈i+1 + Cadj λ̇i+1 + Kadj λi+1 = Fi+1 . (43) Inserting Eqs. 41 and 42 into Eq. 43, after some simplifications, yields, Madj − t (1 − γ )Cadj + ( t)2 (0.5 − β)Kadj λ̈i adj = Fi − Cadj λ̇i+1 − tγ λ̈i+1 −Kadj λi+1 − t λ̇i+1 + ( t)2 β λ̈i+1 . (44) Similarly to the forward problem, we solve the adjoint problem using the stepbystep time integration procedure outlined by Eqs. 41–44. d. (46) (39) Comparison of Eqs. 8a and 8b with Eqs. 34a and 34b reveals that the system matrices for the forward and adjoint problems are related, as per: (45) Loadlocation optimization: Here, we seek the optimal load locations for maximizing the motion metric in the inclusion. The vector of control parameters, used in Eqs. 23 and 24, is given by η = η1 η2 · · · ηns . ∇(ηm ) A = ⎧ ⎫ nf ⎬ ∂θmk (x) T ⎨ T λub k (t) ξmj τj (t)dt d. ⎩ 0 ⎭ ∂ηm load j =1 (47) Detailed derivations of the control problems are given in Appendix B. The summary of the inversion procedure is discussed next. We solve the forward problem using an initial guess of the source characteristics to obtain the state variables. We use the state variables to form the body forces that drive the adjoint problem. Upon solution, the adjoint problem yields the values of the adjoint variable on the loaded boundary load . These values are used to compute the reduced gradient(s) in Eqs. 46 and/or 47. The procedure used for updating the vector of location parameters (η) is discussed next; the identical procedure is used for updating the vector of time parameters (ξ ). Let glk be the discretized reduced gradient for the load location vector obtained at the end of the kth inversion iteration. Thus, glk = (∇(η) A)k . (48) glk is used to compute a search direction slk in the space Comput Geosci of source characteristics. Here, we use the search direction given by the conjugate gradient method. slk = −glk , if k = 0, glk · glk slk = −glk + , otherwise. gl(k−1) · gl(k−1) (49) We, then, update the source characteristics by taking a step towards the search direction. The magnitude of the step is decided by the step length αlk , i.e., we set ηk+1 = η k + αlk slk . The iteration counter is updated (k ← k + 1), and the procedure is repeated until convergence is reached. Equation 51, which is tantamount to the discretized forward problem, will be used to formulate the inverse problem. Recall that the objective functional was given by Eq. 25 in its continuous form. The corresponding discrete form, modulo the mass matrix for compactness, is given by (compare with (50)): Ld = T 0 1 ρa 4.1.7 Discrete objective functional T T t [ 12 u˙ a0 u˙ a0 + 12 u˙ aN u˙ aN + N −1 T u˙ ai u˙ ai ] i=1 To complete the discrete problem, the objective functional (25) must also be discretized. Let L̂ denote the discrete L; then 1 L̂ = T , (50) T ˙ ua Ma u˙ a dt 0 where Ma is the mass matrix of the inclusion (Appendix A). 4.2 Discretizethenoptimize (DTO) approach Following the procedure outlined in [11, 15], we begin with the spatiotemporally discretized forward problem. The objective functional is cast in a discrete form and augmented by weighing the governing equations by discrete Lagrange multipliers. The resulting discrete Lagrangian is then subjected to the firstorder optimality conditions to arrive at the KKT system. Specifically, the timemarching scheme outlined in Eqs. 14 to 19 can be represented as a linear system of equations given by, Qu = f , 1 T ρa u˙ a u˙ a dt (51) 1 = , T ρa u B ua u (55) where B ua is a block diagonal matrix with tB i on its diagonals; B i are square matrices that are zero everywhere except on diagonals that correspond to the u̇a degrees of freedom. As before, we superimpose the governing equations (51), weighted by the discrete Lagrange multipliers p, on the objective functional to obtain the discrete Lagrangian, which is to be minimized. Thus, the constrained minimization problem can now be stated as, min Ad (u, p, f ) = Ld − p T (Qu − f ), f (56) where p = [λ̈0 λ̇0 λ0 λ̈1 λ̇1 λ1 · · · λ̈N λ̇N λN ]T , λ = [λua1 λua2  reg reg λub λub 1 2 λi = λ, at t = i t.  (57) PML PML λub λub λS11 λS22 λS12 ]T , (58) 1 2 (59) The firstorder optimality conditions can now be obtained by taking derivatives of Ad with respect to u, λ and the forceparameters ξ , η. where u = [d̈0 ḋ0 d0 d̈1 ḋ1 d1 · · · d̈N ḋN dN ]T , (52) f = [F0 ḋ0 d0 F1 0 0 · · · FN 0 0]T , (53) Q= ⎡ ⎤ M C K 0 00··· 0 0 0 0 00 ⎢ 0 I 0 0 00··· 0 0 0 0 0 0⎥ ⎢ ⎥ ⎢ 0 0 I 0 00··· 0 0 0 0 0 0⎥ ⎢ ⎥ ⎢ L 0 0 0 0 0⎥ ⎢ 2 L1 L0 Meff 0 0 · · · 0 ⎥ ⎢ −b I −I 0 −b I I 0 · · · 0 0 0 0 0 0⎥ ⎢ 3 ⎥ 4 ⎢ ⎥ 0 0 0 0 0⎥ . ⎢ −b1 I −b0 I −I −b2 I 0 I · · · 0 ⎢ . .. .. .. .. .. . . . .. .. .. .. .. ⎥ ⎢ . ⎥ ⎢ . . . . . . . .. . . . . .⎥ ⎢ ⎥ ⎢ 0 0 0 0 0 0 · · · L2 L1 L0 Meff 0 0 ⎥ ⎢ ⎥ ⎣ 0 0 0 0 0 0 · · · −b3 I −I 0 −b4 I I 0 ⎦ 0 0 0 0 0 0 · · · −b1 I −b0 I −I −b2 I 0 I (54) 4.2.1 State problem ∂ Ad = 0 =⇒ Qu = f , ∂p (60) which is the forward problem, given by Eq. 51. 4.2.2 Adjoint problem ∂ Ad −2B ua u = 0 =⇒ QT p = . ∂u ρa (uT B ua u)2 (61) Equation 61 represents the adjoint problem associated with the inverse problem of interest. Since the adjoint problem involves transpose of Q, we solve it by marching backwards in time. That is, from the last three rows of the matrix (61), Comput Geosci we get, (for i = N ), (update) λN = λuN , (62a) (update) λ̇N = λvN , (62b) (solve) MTeff λ̈N = λaN + b4 λ̇N + b2 λN . (62c) For (N − 1) ≤ i ≤ 1, we update λ̇i , λi and solve for λ̈i using the following, (update) λi = λui + λi+1 − LT0 λ̈i+1 , (63a) (update) λ̇i = λvi + b0 λi+1 + λ̇i+1 − LT1 λ̈i+1 , (63b) (solve) MTeff λ̈i = λai + b1 λi+1 + b3 λ̇i+1 −LT2 λ̈i+1 + b2 λi + b4 λ̇i . (63c) For i = 0, we get (solve)MT λ̈0 = λa0 + b1 λ1 + b3 λ̇1 − LT2 λ̈1 , (update)λ0 = (update)λ̇0 = λu0 + λ1 − LT0 λ̈1 − KT λ̈0 , λv0 + b0 λ1 + λ̇1 − LT1 λ̈1 − CT λ̈0 . provides the gradient of the Lagrangian with respect to the control parameter η. For any given loadlocation parameter ηm , we get T ∂Fi ∂ Ad = λ̈i . ∂ηm ∂ηm N For a load acting in the xp direction, we update each element, ηm , of the control parameter vector η using: nf N ∂θmp (x) ∂ Ad T = λ̈k,load ξmj τj (k t) d, ∂ηm ∂ηm load k=0 j =1 ⎞ ⎡⎛ nf N ⎣⎝ ξmj τj (k t)⎠ · = (64a) k=0 (64b) λ̈k,load (64c) In Eqs. 62–64, λui = λai = 0, 0 ≤ i ≤ N, (65a) v 2 ˙ T ˙ λi = −ρa t Ld [ua1 i ua2 i  0 0  0 0 0 0 0] , i = 0, N, (65b) λvi = −2ρa t L2d [u˙ a1 i u˙ a2 i  0 0  0 0 0 0 0]T , 0 < i < N. (65c) It can be seen from the above equations that the adjoint problem of the DTO approach is driven by the prescription of velocitylike adjoint variables (λvi ) at each time step. The value of λvi depends on the value of the objective functional as well as the value of the velocity vector for the target (u̇ai ). (70) i=0 T j =1 load ∂θmp (x) d . ∂ηm Effectively, Eqs. 68 and 71 of the DTO approach are used in lieu of Eqs. 46 and 47 of the OTD approach. Details of the derivations are given in Appendix B. As before, we start the inversion process with an initial guess of source characteristics. We solve the forward problem, compute the objective functional (Ld ), and obtain the velocity field in the target (u̇ai ). Ld and u̇ai are used to compute the velocitylike adjoint variables (λvi ) that drive the adjoint problem. Solution of the adjoint problem yields the values of the accelerationlike adjoint variables on the loaded boundary (λ̈k,load ). We, then, use Eqs. 68 and/or 71 4.2.3 Control problem(s) Time signal optimization: ∂ Ad ∂ξ = pT ∂f ∂ξ (66) provides the gradient of the Lagrangian with respect to the control parameter ξ . For any given nodalexcitation parameter ξmn , we get T ∂Fi ∂ Ad = λ̈i . ∂ξmn ∂ξmn N (67) i=0 For a load acting in the xp direction, we update each element, ξmn , of the control parameter vector ξ using: N ∂ Ad T = λ̈k,load θmp (x) τn (k t) d, ∂ξmn load k=0 = N τn (k t) T λ̈k,load k=0 load θmp (x) d. (68) Load location optimization: ∂ Ad ∂η = pT ∂f ∂η (69) (71) Fig. 4 Geological formation model Comput Geosci Table 1 Summary of numerical experiments Numerical experiment number Loading direction Spatial description of loads Optimization algorithm Optimization variable KEinc (J/m) 1 2 Vertical Horizontal Eq. 23 Eq. 24 OTD DTO Timesignal Timesignal Timesignal then loadlocation (sequentially) Timesignal and loadlocation (simultaneously) Timesignal Timesignal Timesignal Timesignal 1.63 3.22 3 3 4 4 Vertical Horizontal Horizontal Horizontal Eq. 23 Eq. 23 Eq. 23 Eq. 23 to compute the discrete reduced gradients. Equations 48 and 49 are then employed to update the source characteristics, and we follow the iterative inversion procedure until convergence is reached. 5 Numerical experiments In this section, we test the inversion algorithm by conducting numerical experiments. In our experiments, we use either horizontally or verticallypolarized surface loads. The maximum amplitude is set at 50 kN/m2 . Figure 4 shows the geometry as well as the P and Swave speeds for the geological formation model used to illustrate the capabilities of the methodology. The mass densities (ρa and ρb ) for Fig. 5 Frequency sweep results for uniform vertical load OTD OTD OTD DTO 3.85 4.75 1.63 3.10 3.09 3.10 all materials in the profile shown in Fig. 4 are 2200 kg/m3 . The profile has an elliptical target inclusion whose center is located at 260 m below the ground surface. The major axis of the inclusion is 30 m long, while the minor axis length is 15 m. In all experiments, we use a timestep t = 0.001 s. We define the following metrics to evaluate the performance of the proposed inversion method: if u(t) denotes the displacement at a computational node at time t, the timeaveraged kinetic energy (KETA ) at that node is defined as T KETA = 0 1 ρ [u̇(t) · u̇(t)] dt / T , 2 (72) where ρ is the mass density. Timeaveraged kinetic energy, further integrated over the target inclusion, is defined as Fig. 6 Frequency sweep results for uniform horizontal load Comput Geosci of the inversion technique with that of the frequency sweep approach. To this end, we conduct a frequency sweep for the geological formation model (Fig. 4) by applying first a uniform surface load, either horizontally or vertically polarized. The load has an amplitude of 2 kN/m2 everywhere along the line x2 = 0, and, temporally, it varies as a sine function of the driving frequency. We calculate the objective functionals (L̂ and Ld ) for frequencies ranging between 0.5 and 60 Hz. Figure 5 shows that for the verticallyacting load, there are many local minima at various driving frequencies. The global minimum occurs at 30.5 Hz. Notice that the value of the local minimum at 53 Hz is very close to that of the global minimum. For the horizontal load (Fig. 6), the objective functionals have a global minimum at 18.3 Hz. The figure also shows that both objective functionals exhibit a second strong local minimum at a driving frequency of 34.5 Hz. 5.1 Experiment 1—Source signal optimization In this experiment, we compare the results of the time signal optimization with those of the frequency sweep procedure. The performance is judged by comparing KEinc for the two approaches. To this end, we apply three vertical loads centered at (10,0)m, (−15, 0)m, and (0,0)m on the surface (x2 = 0) of the geological formation model. We use Eq. 23 to specify the spatial variability of the loads. If one were to choose time signals based on the results of the frequency sweep (Fig. 5), there are (at least) two possibilities: (i) monochromatic signals having 30.5 Hz as their dominant frequency and (ii) dichromatic signals having 30.5 and 53 Hz as dominant frequencies. To measure the energy delivery performance of these signals, we construct Fig. 7 Experiment 1—TS1 and TS2 signals KEinc . Thus, KEinc = = a T 0 T 0 1 ρa [u̇a (t) · u̇a (t)] dt d / T 2 1˙ ua (t)T Ma u˙ a (t)dt / T , 2 (73) where u˙ a (t) is the velocity vector corresponding to the computational nodes in the inclusion, and Ma is the mass matrix of the inclusion, as defined in Appendix A. We remark that the units of KETA are joule per cubic meter and those of KEinc are joule per meter. In the experiments, we use plots of KETA and values of KEinc to compare the degree of energy focusing achieved by various techniques. We report the results of four numerical experiments, which are used to highlight the performance of the inverse source approach; Table 1 shows a summary of the numerical experiments features, and the kinetic energy KEinc density resulting from each experiment. As discussed earlier, a frequency sweep is one possibility for deciding suitable temporal characteristics of the wave sources. In our experiments, we compare the performance Fig. 8 Frequency sweep results for three vertical loads located at (10,0)m, (−15, 0)m, and (0,0)m FFT of f1(t) 50 0 −50 0 0.2 0.4 Time [s] 0.6 0.8 1 FFT of f2(t) 50 0 −50 0 0.2 0.4 Time [s] 0.6 0.8 1 50 FFT of f3(t) 2 f3(t) [kN/m ] 2 f2(t) [kN/m ] 1 f (t) [kN/m2] Comput Geosci 0 −50 0 0.2 0.6 0.4 0.8 1 1 0.5 0 0 10 20 40 30 Frequency [Hz] 50 60 0 10 20 40 30 Frequency [Hz] 50 60 0 10 20 40 50 60 1 0.5 0 1 0.5 0 30 Frequency [Hz] Time [s] Fig. 9 Experiment 1–Initial guess FFT of f1(t) 50 0 −50 0 0.2 0.4 Time [s] 0.6 0.8 1 FFT of f2(t) 50 0 −50 0 0.2 0.4 Time [s] 0.6 0.8 1 50 0 −50 at the same three locations as before, the same (approximately) minima would be recovered, as the plot in Fig. 8 depicts. Of interest is how the inversion procedure will perform. Thus, next, we conduct the time signal optimization using the OTD approach (without changing the load locations). The initial guess for the time signals of the loads is shown in Fig. 9. It can be seen in Fig. 9b that the initial guesses have broad frequency support. The time signals obtained following the optimization are shown in Fig. 10, and as it can be seen in Fig. 10b for all three loads there are clearly two dominant frequencies at, approximately, 30 and 53.5 Hz. That is, the optimizer converged nicely to signals that, at a minimum, contain dominant contributions from the two “best” frequencies revealed by the frequency sweep, without any prior biasing of the inversion process. Notice also that the converged time signals have a rich content, with contributions from many frequencies other than those corresponding to the strongest minima seen in the frequency sweep. It can also be observed in Fig. 10b that the frequency spectra for FFT of f3(t) f3(t) [kN/m2] f2(t) [kN/m2] f1(t) [kN/m2] signals TS1 and TS2 using piecewise quadratic Lagrange polynomials (22), corresponding to the monochromatic and dichromatic signal, respectively. The time signals (TS1 and TS2) and their frequency spectra are shown in Fig. 7 measure KEinc for these two cases. Sources operating with TS1 as their driving signals are able to achieve KEinc = 1.06 J/m, whereas those with TS2 as the driving signals deliver KEinc of 1.39 J/m. In other words, operating three vertical loads at the combination of two sine functions, each operating at the two “best” frequencies obtained by the frequency sweep will result in an increase of the energy delivery of about 36 % over the energy delivered when operating the same three loads at the single best frequency resulting from the sweep. We remark that loading the entire width of the computational domain generates waves having, initially, plane wavefronts in the domain of interest. The waves emitted by the finitewidth sources are nearcylindrical in nature. This difference alone will not alter in any significant way the sweep’s results. For example, if three 2mwide verticallyacting loads operate 0 0.2 0.4 Time [s] 0.6 Fig. 10 Experiment 1— Converged time signals 0.8 1 1 0.5 0 0 10 20 30 40 Frequency [Hz] 50 60 0 10 20 30 40 Frequency [Hz] 50 60 0 10 20 30 40 Frequency [Hz] 50 60 1 0.5 0 1 0.5 0 Comput Geosci Fig. 11 Experiment 1—Timeaveraged kinetic energy 0 −50 0 0.2 0.4 Time [s] 0.6 0.8 1 0 0 0.2 0.4 Time [s] 0.6 0.8 1 50 0 −50 0.5 0 0 10 30 20 Frequency [Hz] 40 50 0 10 20 30 Frequency [Hz] 40 50 0 10 20 30 Frequency [Hz] 40 50 1 2 FFT of f (t) 50 −50 1 1 FFT of f (t) 50 This experiment highlights the key differences between the frequency sweep and time signal optimization procedures. The time signal optimization procedure was able to deliver about 17 % more kinetic energy (in the timeaveraged sense) than the best the frequency sweep information could do based on the dichromatic signal. Note that we conducted our experiment for a duration T = 1 s: in a field implementation of the seismic EOR, where the stimulation may be applied for days, 17 % improvement in the efficiency is significant. For the model formation used FFT of f (t) 3 3 f (t) [kN/m2] f (t) [kN/m2] 2 1 f (t) [kN/m2] the converged signals of the three loads are not identical. The difference can be attributed to the complex interacting wave patterns. Simply put, the inverse source problem is aware of these patterns by virtue of the forward (state) problem solution, and uses the information to steer the time signals to maximize energy focusing. Figure 11 depicts the distribution of timeaveraged kinetic energy (KETA ) in the computational domain before and after optimization; the KEinc delivered to the target due to the optimized time signals is 1.63 J/m. 0 0.2 Fig. 12 Experiment 2—Initial guess 0.4 Time [s] 0.6 0.8 1 0.5 0 1 0.5 0 0 −50 0 0.2 0.4 Time [s] 0.6 0.8 1 0 0 0.2 0.4 Time [s] 0.6 0.8 1 50 0 0 10 20 30 Frequency [Hz] 40 50 0 10 30 20 Frequency [Hz] 40 50 0 10 20 30 Frequency [Hz] 40 50 1 0.5 0 1 3 0 −50 0.5 2 FFT of f (t) 50 −50 1 1 FFT of f (t) 50 FFT of f (t) 3 f (t) [kN/m2] 2 f (t) [kN/m ] 2 1 f (t) [kN/m2] Comput Geosci 0 0.2 0.4 0.6 0.8 1 0.5 0 Time [s] Fig. 13 Experiment 2—Time signals after optimization (time signal only) in our experiment, the choice of time signals is not obvious from the frequency sweep results. A frequency sweep procedure, which is based on loading the entire surface of the (computational) domain, can become blind to the complex interference patterns created by waves emitted by a number of sources having finite widths. The degree of blindness depends, to a certain extent, on the heterogeneity of the geostructure. A frequency sweep procedure that involves using loads with smaller widths (e.g., 2 m) requires a combined frequencyandlocation sweep, which is computationally expensive. These observations indicate the utility of the proposed inverse source approach in deciding the time signals driving wave sources used for focusing vibrational energy. 5.2 Experiment 2–Source location and signal optimization In this experiment, we discuss the optimization of the wave sources’ spatial and temporal characteristics. The time signal and load location optimization can be performed either sequentially or simultaneously. In a sequential process, the time signals are optimized first, without changing the load Fig. 14 Experiment 2—Timeaveraged kinetic energy (sequential optimization) FFT of f1(t) 50 0 −50 0 0.2 0.4 Time [s] 0.6 0.8 1 FFT of f2(t) 50 0 −50 0 0.2 0.4 Time [s] 0.6 0.8 1 50 FFT of f3(t) f3(t) [kN/m2] 2 f2(t) [kN/m ] 2 f1(t) [kN/m ] Comput Geosci 0 −50 0 0.2 0.4 0.6 0.8 1 Time [s] 1 0.5 0 0 10 20 30 Frequency [Hz] 40 50 0 10 20 30 Frequency [Hz] 40 50 0 10 20 40 50 1 0.5 0 1 0.5 0 30 Frequency [Hz] Fig. 15 Experiment 2—Time signals after simultaneous optimization locations. The converged (optimized) signals are then used, while the source locations are optimized, until convergence is reached. By contrast, the simultaneous optimization process is initiated with guesses for both time signals and load locations, and both the spatial and temporal characteristics are updated during every inversion iteration until convergence is reached. In this experiment, we compare the performance of the optimization algorithm with that of the frequency sweep. Fig. 16 Experiment 2— Timeaveraged kinetic energy (simultaneous optimization) Sequential optimization: To test the performance of the sequential optimization procedure, we use three horizontal loads on the surface (x2 = 0) of the geological formation model. We use Eq. 24 to describe the spatial variability of the surface tractions. The loads are centered at (7,0)m, (−5, 0)m, and (0,0)m, i.e., η1 = 7m, η2 = −5m, and η3 = 0m in Eq. 24. The temporal optimization process is initiated with the time signals shown in Fig. 12. Their Fourier transforms (Fig. 12b) show the presence of a wide range Comput Geosci Fig. 17 Experiment 3—Frequency sweep results for vertical and horizontal loads of frequencies. We, next, perform the time signal optimization using the DTO procedure. The optimizer converges to the time signals shown in Fig. 13. The plots for the timeaveraged kinetic energy (KETA ) for the initially guessed and Fig. 18 Experiment 3—Time signals and frequency spectra for the converged time signals are compared in Fig. 14. KEinc for the converged time signals is 3.22 J/m. A closer look at the frequency spectra of the converged time signals (Fig. 13b) reveals that the optimizer converged to signals with dominant frequencies at 18.3 and 34.5 Hz. Recall that these frequencies correspond to the minima revealed by the frequency sweep (Fig. 6). Next, we carry out the optimization for the load locations. We use the converged time signals (Fig. 13) as fixed driving signals for the loads. We allow the optimizer to change the locations of loads in order to improve the energy delivery to the target inclusion. The loads move to the left and the converged locations are given by (−18.1, 0)m, (−26.5, 0)m, and (−25.2, 0)m. KEinc at the end of sequential optimization is 3.85 J/m. The plot of timeaveraged kinetic energy is shown in Fig. 14c. Note that the algorithm was able to deliver 20 % more energy merely by changing the locations of loads, i.e., by allowing the waves to constructively interfere in the target inclusion. In the sequential optimization process, fixing the frequency content of time signals while searching for optimal locations may hinder the desired focusing. We have Comput Geosci Fig. 19 Experiment 3—Timeaveraged kinetic energy observed that the optimal frequency content of the sources depends on their location. To investigate this effect, we perform simultaneous optimization, which retains the flexibility of temporal and spatial characteristics during the entire inversion process. Simultaneous optimization: We begin with an initial guess of time signals (Fig. 12), and load locations, which are the same as those used in the sequential optimization case, i.e., (7,0)m, (−5, 0)m, and (0,0)m). During each optimization iteration, we update the time signals, as well as the load locations, using the search direction for the respective Fig. 20 Experiment 4—Initial guess control variable. The converged time signals are shown in Fig. 15, and the timeaveraged kinetic energy is plotted in Fig. 16. Note that the converged frequency spectra for the loads in their optimized location (Fig. 15b) are different from those obtained based on time signal optimization for fixed load locations (Fig. 13b). The loads move to the left and the converged locations are (−19, 0)m, (−19.16, 0)m, and (−20.04, 0)m. KEinc for the converged time signals and locations is 4.75 J/m, which is about 50 % higher than that achieved by performing time signal optimization with constant load locations (Fig. 14b), and about 23 % higher than that achieved by the sequential time signal and load location optimization (Fig. 14c). Comput Geosci to the target formation by placing loads at advantageous locations. In the sequential optimization, the time signals are kept constant while performing load location optimization. Since the optimal frequency content of a load depends on its location, the sequential optimization cannot adequately compensate. The simultaneous optimization, on the other hand, is able to adjust both the temporal and spatial characteristics of loads and achieves (about 50 %) better energy delivery when compared to the time signal optimization. 5.3 Experiment 3—Polarization effect Fig. 21 Experiment 4—Growth of timeaveraged kinetic energy for the inversion iterations This experiment further highlights the importance of load location optimization: in the sequential optimization, we observed that about 20 % more energy can be delivered Fig. 22 Experiment 4—Time signals and frequency spectra Experiment 3 is designed to determine the effect of polarization of applied tractions on the wave energy focusing. Fundamental solutions for point loads on a (homogeneous) elastic halfspace [17] indicate that, for example, for a Poisson’s ratio of 0.2, a horizontal load radiates about 50 % of the energy via P and Sbody waves. On the other hand, a vertical point load imparts about 70 % of its energy in the form of Rayleigh surface waves. Thus, for a homogeneous, elastic halfspace, horizontal (point) loads are more Comput Geosci Fig. 23 Experiment 4—Timeaveraged kinetic energy efficient at delivering energy to deeply situated subterranean formations. To capture the behavior for a heterogeneous halfspace excited by loads having finite widths, we plot the frequency sweep of KEinc for horizontal and vertical loads (Fig. 17). Note that this is merely a different representation of the frequency sweeps reported earlier (Figs. 5 and 6). As seen in Fig. 17, the horizontally polarized load is able to deliver more vibrational energy to the target than the vertical load for a wide range of frequencies. To further illustrate the effects of loading direction, we compare the results of time signal optimization performed with three horizontal and three vertical loads. The locations of the loads are kept constant. We use the OTD approach to carry out the time signal optimization. The initial guess of time signals for both the horizontally and vertically polarized loads is shown in Fig. 9. The loads are centered at (10,0)m, (−15, 0)m, and (0,0)m. We use Eq. 23 to specify the spatial variability of the loads. The converged time signals and their frequency content are depicted in Fig. 18. Figure 19 shows the plots of KETA . KEinc for the vertical loads is 1.63 J/m (Fig. 19a), whereas that for the horizontal loads is 3.09 J/m (Fig. 19b), an almost twofold increase. Thus, our experiment indicates that horizontally polarized loads are able to deliver more vibrational energy. Though, both theory and the numerical experiment reported herein suggest that horizontal loads are preferable for focusing, we note that our model does not account for intrinsic and/or apparent attenuation. A judicious choice about the preferred direction of loading can be made after considering the Qfactors for P and S waves in the geostructure of interest. Table 2 “True” material properties for the geological formation model of Fig. 4 Target inclusion M1 M2 M3 M4 Cp (m/s) Cv (m/s) 762 953 1204 1394 1685 466 583 737 853 1032 Comput Geosci 5.4 Experiment 4—OTD versus DTO In this experiment, we use the same geological model to compare the performance of the OTD and DTO procedures. We carry out time signal optimization for three horizontal loads using both methods. The initial guess for the time signals is shown in Fig. 20. Figure 21 shows that the DTO approach is able to deliver slightly more kinetic energy to the target in the timeaveraged sense, but, in general, the differences are small. Nevertheless, in view of its ease of formulation and slightly better performance, we favor the DTO approach. In Fig. 22, we compare the converged time signals, and Fig. 23 shows the timeaveraged kinetic energy for the two procedures. Both procedures successfully recover the time signals, which exhibit a dominant frequency of, approximately, 18.3 Hz. 5.5 Experiment 5—The effect of formation property uncertainties on energy focusing The working hypothesis of our theoretical development and of the preceding numerical experiments is that the material properties of the geological formation (Fig. 4) are known with confidence. However, in practice, precise knowledge of the properties of the geostructure is elusive. Though a systematic treatment of the effect property uncertainties may have on energy focusing escapes the scope of this communication, in this section, we briefly outline methodologies that would allow the formal treatment of uncertainties and could lead to the quantification of their effects on focusing. The following simple numerical experiment highlights the issues: we first assume that the properties of the geological formation model depicted in Fig. 4 are inaccurate, and that the “true” properties are those given in Table 2, reflecting a 5 to 11 % change in the velocities compared to those in Fig. 4. Next, we are interested in quantifying the energy focusing when the “true” properties are used, while the formation is subjected to the optimal loads we obtained using the inaccurate wave velocity values. We are concerned with both the strength of the focusing, as well as with a potential shift to the focusing area. To this end, we use the spatiotemporally optimized loads obtained in Experiment 2 (Figs. 15 and 16) to excite the geological formation with the “true” properties as given in Table 2. The resulting distribution of KETA is shown in Fig. 24. The resulting value of KEinc was reduced from 4.75 J/m (Experiment 2) to 4.11 J/m, a reduction of about 13 %. It is clear that, while the material interfaces remained unaltered, a nonuniform increase in the properties of the formation led to a reduction of the strength of the focused energy, with no appreciable effects on the focusing area. The result is rather expected and warrants a formal treatment. We briefly outline two systematic approaches that could be used to quantify the effect of property uncertainties on the focusing. – Fig. 24 Timeaveraged kinetic energy for the geological formation with “true” material properties Sensitivity analysis: A sensitivity analysis of the elastodynamic system [19] can be carried out to quantify the dependence of KEinc on the material parameters. A firstorder sensitivity analysis gives rise to a discrete forward problem (13) cast in terms of the unknown ∂u sensitivity variable z, where z = ∂q and q1 is the mate1 rial parameter of interest (for example, q1 = λb1 –the first Lamé parameter of the top layer). The force vector for the sensitivity problem depends on the displacement, velocity, and acceleration fields (u, u̇, ü), and the derivatives of the global system matrices with respect ∂C ∂K to q1 ( ∂M ∂q1 , ∂q1 , ∂q1 ). Upon solution of the sensitivity problem, one recovers the time histories of the sensitivity variable and its time derivatives. The time histories can be used to calculate the derivative of KEinc with inc respect to q1 ( ∂KE ∂q1 ). Thus, the analysis determines the Comput Geosci – sensitivity of KEinc to material parameters of various layers in the formation. Reliability analysis: In this approach, the material properties are treated as random variables endowed with suitable probability distribution functions (PDFs). The deterministic inverse source problem can be solved for the mean values of the material properties to arrive at the spatiotemporally optimal loads (Fmean ). The timeaveraged kinetic energy of the inclusion (KEinc ) is also a random variable, and its value for the mean material properties and corresponding optimal loading (Fmean ) can be computed as KEmean inc . The probability of attaining a KEinc value less than or equal to a predefined threshold δKEmean inc (where 0 < δ < 1), i.e., P[KEinc − δKEmean ≤ 0], when Fmean is used as the surface inc excitation, can be evaluated using firstorder reliability methods (FORM). The reliability analysis requires computation of derivatives of KEinc with respect to inc various material parameters ( ∂KE ∂qi ), similarly to the sensitivity analysis. As explained above, it can be used to calculate the probability of failure in achieving a certain threshold KEinc , given the PDFs for the material parameters in the geostructure. 6 Conclusions We presented an inverse source approach for focusing energy in a target subterranean formation and discussed two possible numerical implementations, DTO and OTD. Through numerical experiments, we provided evidence of the inverse source method’s superiority over frequency sweeps for determining the optimal wave source locations and time signals. Moreover, the method’s ability to resolve the optimal spatiotemporal characteristics of the wave source was shown to result in the strongest energy focusing among all studied alternatives. Our numerical experiments also indicate that horizontally polarized loads may be preferable over vertical loads. Lastly, we observed no significant performance difference between the two numerical implementations we discussed (DTO and OTD). description and definitions of parameters αi (i = 1, 2), βi (i = 1, 2), etc. can be found in [13]. ∂ ∂T DMloc = DT d, DQloc = D d ij ∂xi ∂xj loc loc ∂ T Aij k = ij d Fi = T f i d ∂x k PML load 2μb + λb 4μ b (μb + λb ) PML λb k = 4μ (μ b b + λb ) PML 1 T = k d, 2μ b PML Nik = d, if i = 1, T d, if i = 2, if i = 3, Element matrices for the regular domain are: ρb Mreg 0 0 0 Mreg = , Creg = , 0 ρb Mreg 0 0 Kreg = reg reg reg reg μb Q21 + λb Q12 (2μb +λb ) Q11 +μb Q22 reg reg reg reg . μb Q12 + λb Q21 (2μb +λb ) Q22 + μb Q11 In the PML region, the element matrices can be computed as: ⎡ ⎤ ρb aMPML 0 0 0 0 ⎢ 0 ρb aMPML 0 0 0 ⎥ ⎢ ⎥ ⎢ MPML = ⎢ 0 0 −N1a N2a 0 ⎥ ⎥, ⎣ 0 0 N2a −N1a 0 ⎦ 0 0 0 0 −N3a ⎡ CPML Here, we present concise definitions of element matrices that form the global matrices in Eq. 10. A detailed T For a finite element within the target inclusion, the element mass, damping, and stiffness matrices are given by: ρa Ma 0 0 0 Ma = ,Ca = , 0 ρa Ma 0 0 μa Qa21 + λa Qa12 (2μa +λa ) Qa11 +μa Qa22 Ka= . μa Qa12 +λa Qa21 (2μa + λa )Qa22 + μa Qa11 Acknowledgments The authors’ work was partially supported by an Academic Alliance Excellence grant between the King Abdullah University of Science and Technology in Saudi Arabia (KAUST) and the University of Texas at Austin. The support is gratefully acknowledged. Appendix: A Element matrices k ⎤ ρb bMPML 0 Aα21 0 Aα12 ⎢ 0 ρb bMPML 0 Aα12 Aα21 ⎥ ⎢ ⎥ =⎢ 0 −N1b N2b 0 ⎥ ⎢ Aα21 ⎥, ⎣ 0 Aα12 N2b −N1b 0 ⎦ Aα12 Aα21 0 0 −N3b ⎡ KPML ⎤ ρb cMPML 0 Aβ21 0 Aβ12 ⎢ 0 ρb cMPML 0 Aβ12 Aβ21 ⎥ ⎢ ⎥ ⎢ = ⎢ Aβ21 0 −N1c N2c 0 ⎥ ⎥. ⎣ 0 Aβ12 N2c −N1c 0 ⎦ Aβ12 Aβ21 0 0 −N3c Comput Geosci Optimizethendiscretize approach where λ̈k,load contains the values of adjoint variable corresponding to the degrees of freedom represented by rows of vector on load . Timesignal optimization Loadlocation optimization We take the variation of the Lagrangian with respect to discrete scalar variables ξmn to obtain the gradient of A, i.e., For a given loadlocation parameter ηm , we get Appendix: B Control problem derivations δ(ξij ) A = ∇(ξij ) A = − load T 0 λF · ∂f dt d, or, ∂ξij ∇(ξij =ξmn ) A T ns ∂ λF 1 (x) · λF 2 (x) ∂ξmn i=1 load 0 ⎤ nf θi1 (x) ξij τj (t)dt d ⎦ θi2 (x) j =1 T λub1 (x) θ (x) = · m1 τn (t)dt d λub2 (x) θm2 (x) load 0 ) ( T λub1 (t) θm1 (x) · τn (t)dt d. = λub2 (t) 0 load θm2 (x) = − ∂ Ad ∂ηm = N * i=0 T ∂Fi ∂ηm . λ̈i For each loaded element, ns * Felem θi1 (x, ηi ) elem k,x 1 Fk = = · Felem i=1 load θi2 (x, ηi ) k,x2 nf * ∂Felem k ∂ηm N * k=0 ∂θm1 (x,ηm ) ∂ηm ∂θm2 (x,ηm ) ∂ηm = load T ∂Felem λ̈k ∂ηk m = N * k=0 nf * j =1 Loadlocation optimization Variation of A with respect to the control parameter η is given by δ(ηm ) A = ∇(ηm ) A =− T ∂f dt d ∂η m load 0 nf ∂θmk (x) T T = λub k (t) ξmj τj (t)dtd. ∂ηm 0 load λF · j =1 Discretizethenoptimize approach Timesignal optimization For any given nodalexcitation parameter ξmn , we get ∂ Ad ∂ξmn = N * i=0 T ∂Fi ∂ξmn . λ̈i We recall that for each element Felem k,x1 Felem k,x2 Felem = k ∂Felem k ∂ξmn N k=0 T λ̈k ∂Felem k ∂ξmn ns = i=1 load nf θi1 (x) ξij τj (k t) d, θi2 (x) j =1 θm1 (x) τ (k · t) d, θm2 (x) n N T θm1 (x) = λ̈k,load τn (k t) d, load θm2 (x) = load k=0 ξij τj (k t) d, j =1 T λ̈k,load nf ξmj τj (k t) d, j =1 load ∂θm1 (x,ηm ) ∂ηm ∂θm2 (x,ηm ) ∂ηm · ξmj τj (k t) d, where λ̈k,load contains the values of adjoint variable corresponding to the degrees of freedom represented by rows of vector on load . References 1. Bamberger, A., Chavent, G., Lailly, P.: About the stability of the inverse problem in 1D wave equations application to the interpretation of seismic profiles. Appl. Math. Optim. 5(1), 1–47 (1979) 2. Beresnev, I., Gaul, W., Vigil, R.D.: Direct porelevel observation of permeability increase in twophase flow by shaking. Geophys. Res. Lett. 38(20) (2011). doi:10.1029/2011GL048840 3. Beresnev, I.A., Deng, W.: Viscosity effects in vibratory mobilization of residual oil. Geophys. 4, N79–N85 (2010) 4. Beresnev, I.A., Johnson, P.A.: Elasticwave stimulation of oil production: a review of methods and results. Geophys. 59(6), 1000– 1017 (1994) 5. Bunks, C., Saleck, F., Zaleski, S., Chavent, G.: Multiscale seismic waveform inversion. Geophys. 60(5), 1457–1473 (1995) 6. Deng, W., Cardenas, M.B.: Dynamics and dislodgment from pore constrictions of a trapped nonwetting droplet stimulated by seismic waves. Water Res. Res. 49(7), 4206–4218 (2013) 7. Epanomeritakis, I., Akçelik, V., Ghattas, O., Bielak, J.: A NewtonCG method for largescale threedimensional elastic fullwaveform seismic inversion. Inverse Problems 24(3), 034,015 (2008) 8. Gunzburger, M.D.: Perspectives in flow control and optimization SIAM (2003) 9. Jeong, C., Kallivokas, L.F., Huh, C., Lake, L.W.: Optimization of sources for focusing wave energy in targeted formations. J. Geophys. Eng. 7(3), 242 (2010) Comput Geosci 10. Jeong, C., Kallivokas, L.F., Huh, C., Lake, L.W.: Maximization of oil mobility within a hydrocarbon reservoir for elastic wavebased enhanced oil recovery, In: SPE Annual Technical Conference and Exhibition, SPE 147150. Society of Petroleum Engineers (2011) 11. Kallivokas, L., Fathi, A., Kucukcoban, S., Stokoe(II), K., Bielak, J., Ghattas, O.: Site characterization using full waveform inversion. Soil Dyn. Earthquake Eng. 47, 62–82 (2013) 12. Kostrov, S.A., Wooden, B.O.: Mechanisms, field suitability, and case studies for enhancement of oil recovery and production using insitu seismic stimulation. In: 16th International Symposium on Nonlinear Acoustics (2002) 13. Kucukcoban, S., Kallivokas, L.: A symmetric hybrid formulation for transient wave simulations in pmltruncated heterogeneous media. Wave Motion 50(1), 57–79 (2013) 14. Lake, L.: Enhanced Oil Recovery Prentice Hall Incorporated (1989) 15. Petra, N., Stadler, G.: Model variational inverse problems governed by partial differential equations. Tech. Rep. 115, The Institute for Computational Engineering and Sciences, The University of Texas at Austin (2011) 16. Roberts, P.M., Majer, E.L., Lo, W.C., Sposito, G., Daley, T.M.: An integrated approach to seismic stimulation of oil reservoirs: Laboratory, field and theoretical results from doe/industry collaborations. In: 16th International Symposium on Nonlinear Acoustics (2002) 17. SánchezSesma, F.J., Weaver, R.L., Kawase, H., Matsushima, S., Luzón, F., Campillo, M.: Energy partitions among elastic waves for dynamic surface loads in a semiinfinite solid. Bulletin of the Seismological Society of America 101(4), 1704–1709 (2011) 18. Tarantola, A.: Inversion of seismic reflection data in the acoustic approximation. Geophys. 49(8), 1259–1266 (1984) 19. Tsikas, A.: Sensitivity analysis of elastodynamic systems. Master’s thesis, Department of Civil and Environmental Engineering, Carnegie Mellon University (1997)