Fatemeh Pourahmadian

Bojan B Guzina
email: guzin001@umn.edu

Houssem Haddar

A synoptic approach to the seismic sensing of heterogeneous fractures: from geometric reconstruction to interfacial characterization

Keywords: inverse scattering, elastic waves, fractures, heterogeneous contact condition, specific stiffness, hydraulic fractures

A non-iterative waveform sensing approach is proposed toward (i) geometric reconstruction of penetrable fractures, and (ii) quantitative identification of their heterogeneous contact condition by seismic i.e. elastic waves. To this end, the fracture support Γ (which may be non-planar and unconnected) is first recovered without prior knowledge of the interfacial condition by way of the recently established approaches to noniterative waveform tomography of heterogeneous fractures, e.g. the methods of generalized linear sampling and topological sensitivity. Given suitable approximation Γ of the fracture geometry, the jump in the displacement field across Γ i.e. the fracture opening displacement (FOD) profile is computed from remote sensory data via a regularized inversion of the boundary integral representation mapping the FOD to remote observations of the scattered field. Thus obtained FOD is then used as input for solving the traction boundary integral equation on Γ for the unknown (linearized) contact parameters. In this study, linear and possibly dissipative interactions between the two faces of a fracture are parameterized in terms of a symmetric, complex-valued matrix K(ξ) collecting the normal, shear, and mixedmode coefficients of specific stiffness. To facilitate the high-fidelity inversion for K(ξ), a 3-step regularization algorithm is devised to minimize the errors stemming from the inexact geometric reconstruction and FOD recovery. The performance of the inverse solution is illustrated by a set of numerical experiments where a cylindrical fracture, endowed with two example patterns of specific stiffness coefficients, is illuminated by plane waves and reconstructed in terms of its geometry and heterogeneous (dissipative) contact condition.

Introduction

Geometric and interfacial properties of fractures in rock and like quasi-brittle solids (e.g. concrete and composites) are the subject of critical importance to a wide spectrum of scientific and technological facets of our society including energy production from natural gas and geothermal resources [START_REF] Baird | Monitoring increases in fracture connectivity during hydraulic stimulations from temporal variations in shear wave splitting polarization[END_REF][START_REF] Verdon | Measurement of the normal/tangential fracture compliance ratio (z N /z T ) during hydraulic fracture stimulation using s-wave splitting data[END_REF][START_REF] Taron | Coupled mechanical and chemical processes in engineered geothermal reservoirs with dynamic permeability[END_REF], seismology [START_REF] Mclaskey | Fault healing promotes high-frequency earthquakes in laboratory experiments and on natural faults[END_REF], non-destructive evaluation (NDE) [START_REF] Doqueta | Experimental and numerical study of crack healing in a nuclear glass[END_REF], hydrogeology [START_REF] Cook | Natural joints in rock: Mechanical, hydraulic and seismic behaviour and properties under normal stress[END_REF], and environmental protection [START_REF] Place | wet fault or dry fault? a laboratory approach to remotely monitor the hydro-mechanical state of a discontinuity using controlled-source seismics[END_REF]. One particular quantity embodying the fracture's linearized contact law is the so-called specific stiffness matrix K, relating the surface traction to the jump in displacement across the interface. In practical terms, the spatial heterogeneity of the contact parameters reflected in K(ξ)due to e.g. variable distribution of normal stress [START_REF] Martel | Mechanics of slip and fracture along small faults and simple strike-slip fault zones in granitic rock[END_REF], may be responsible for progressive failure along discontinuities that may occur well before the frictional resistance of the entire interface is surpassed [START_REF] Mutlu | Slip initiation on frictional fractures[END_REF]. This may lead to a catastrophic failure of dams, tunnels and slopes [START_REF] Gu | Progressive failure of lower San Fernando dam[END_REF][START_REF] Eberhardt | Numerical analysis of initiation and progressive failure in natural rock slopes -the 1991 Randa rockslide[END_REF][START_REF] Bakun-Mazor | Dynamic sliding of tetrahedral wedge: the role of interface friction[END_REF], particularly when fractures show slip-weakening behavior [e.g. 25] while the underlying design is based upon averaged contact properties. It is further shown in [START_REF] Hedayat | Precursors to the shear failure of rock discontinuities[END_REF] that the onset of slip along an interface can be identified via temporal variation of the fracture's specific stiffness in shear direction. Accordingly, real-time monitoring of K(ξ)'s evolution may not only serve as an early indicator of the interfacial instability and failure, but may also help understand the mechanism of shallow earthquakes [START_REF] Jaeger | Fundamentals of Rock Mechanics[END_REF]. Active seismic sensing of fractures' contact condition has also come under a spotlight in energy production from unconventional resources [START_REF] Baird | Monitoring increases in fracture connectivity during hydraulic stimulations from temporal variations in shear wave splitting polarization[END_REF][START_REF] Verdon | Measurement of the normal/tangential fracture compliance ratio (z N /z T ) during hydraulic fracture stimulation using s-wave splitting data[END_REF], owing to the strong correlation between the hydraulic conductivity of a fracture network and the spatiotmporal variations of K [START_REF] Pyrak-Nolte | Approaching a universal scaling relationship between fracture stiffness and fluid flow[END_REF]. For sensing purposes, one should bear in mind that the fracture's response to given activation is driven not only by its contact condition K, but also by its geometry which is not limited to the planar condition [e.g. [START_REF] Brace | A note on brittle crack growth in compression[END_REF][START_REF] Barton | The shear strength of rock and rock joints[END_REF][START_REF] Thomas | The geometry of echelon fractures in rock: implications from laboratory and numerical experiments[END_REF]. Thus, the objective of this research is to establish a robust framework for the seismic waveform sensing of heterogeneous fractures, that is capable of resolving both their geometry and interfacial characteristics without iterations. Traditionally, seismic waves have been used in the context of acoustic emission (AE) [START_REF] Lockner | The role of acoustic emission in the study of rock fracture[END_REF][START_REF] Shah | Damage mechanisms in stressed rock from acoustic emission[END_REF] to monitor the progression of evolving fractures via the detection of underlying microseismic events, whose energy -captured by the receivers -is used to track the failure process. Such "passive" sensing approach is, however, ineffective when trying to either image in situ fractures or to assess the fracture's interfacial condition. Another approach, motivating this research, is the concept of active seismic sensing applied to fracture identification. This approach, where the discontinuity is "illuminated" by an external seismic source, carries the potential of simultaneous fracture imaging and characterization thanks to the sensitivity of scattered waves to the interfacial condition.

In general the relationship between the wavefield scattered by an obstacle and its geometry and mechanical characteristics is nonlinear, which invites two overt solution strategies: (i) linearization via e.g. Born approximation and ray theory [START_REF] Bleistein | Mathematics of Multidimensional Seismic Imaging, Migration, and Inversion[END_REF], or ii) pursuit of the nonlinear minimization approach [START_REF] Virieux | An overview of full-waveform inversion in exploration geophysics[END_REF]. Over the past two decades, however, a number of sampling methods have emerged that consider the nonlinear nature of the inverse problem in an iteration-free way. In the context of extended scatterers, examples of such paradigm include the linear sampling method (LSM) [START_REF] Colton | A simple method for solving inverse scattering problems in the resonance region[END_REF][START_REF] Cakoni | The linear sampling method for cracks[END_REF][START_REF] Guzina | A linear sampling approach to inverse elastic scattering in piecewise-homogeneous domains[END_REF], the factorization method (FM) [START_REF] Kirsch | The Factorization Method for Inverse Problems[END_REF][START_REF] Boukari | The factorization method applied to cracks with impedance boundary conditions[END_REF], the generalized linear sampling method (GLSM) [START_REF] Audibert | Qualitative Methods for Heterogeneous Media[END_REF][START_REF] Pourahmadian | Generalized linear sampling method for elastic-wave sensing of heterogeneous fractures[END_REF][START_REF] Cakoni | Inverse Scattering Theory and Transmission Eigenvalues[END_REF], the concept of topological sensitivity (TS) [START_REF] Guzina | Topological derivative for the inverse scattering of elastic waves[END_REF][START_REF] Chikichev | Generalized topological derivative for the navier equation and inverse scattering in time domain[END_REF][START_REF] Bellis | Crack identification by 3D time-domain elastic or acoustic topological sensitivity[END_REF] and the subspace migration technique [START_REF] Park | Multi-frequency subspace migration for imaging of perfectly conducting, arc-like cracks in full-and limited-view inverse scattering problems[END_REF]. Among these, the TS and GLSM approaches have been recently adapted to permit elastic-wave sensing of heterogeneous fractures [START_REF] Pourahmadian | On the elastic-wave imaging and characterization of fractures with specific stiffness[END_REF][START_REF] Pourahmadian | Generalized linear sampling method for elastic-wave sensing of heterogeneous fractures[END_REF].

As indicated earlier, there is a mounting interest in medical diagnosis, target recognition, seismology, and energy production [START_REF] Cakoni | A Qualitative Approach to Inverse Scattering Theory[END_REF][START_REF] Minato | Imaging and characterization of a subhorizontal non-welded interface from point source elastic scattering response[END_REF][START_REF] Verdon | Measurement of the normal/tangential fracture compliance ratio (z N /z T ) during hydraulic fracture stimulation using s-wave splitting data[END_REF][START_REF] Place | wet fault or dry fault? a laboratory approach to remotely monitor the hydro-mechanical state of a discontinuity using controlled-source seismics[END_REF] to develop hybrid sensing schemes that also reveal the boundary condition of hidden anomalies. In this vein, it is noteworthy that some anomaly-indicator functionals -such as those featured by the FM [START_REF] Kirsch | The Factorization Method for Inverse Problems[END_REF], (G)LSM [START_REF] Cakoni | The determination of anisotropic surface impedance in electromagnetic scattering[END_REF][START_REF] Pourahmadian | Generalized linear sampling method for elastic-wave sensing of heterogeneous fractures[END_REF] and TS [START_REF] Pourahmadian | On the elastic-wave imaging and characterization of fractures with specific stiffness[END_REF] are largely insensitive to the (unknown) boundary condition of a hidden anomaly. For instance, [START_REF] Cakoni | A Qualitative Approach to Inverse Scattering Theory[END_REF][START_REF] Ben Hassen | Application of the linear sampling method to identify cracks with impedance boundary conditions[END_REF] demonstrate that the LSM is successful in reconstructing electromagnetic obstacles and cracks regardless of their boundary condition. Recent advancements on the recovery of boundary (or interfacial) conditions are, on the other hand, mostly optimizationbased as proposed in the context of acoustic and electromagnetic inverse scattering. Hitherto, a variational method is proposed in [START_REF] Colton | The determination of the surface impedance of a partially coated obstacle from far-field data[END_REF] within the framework of the LSM to determine the essential supremum of the electrical impedance at the boundary of partially coated obstacles. More recently, [START_REF] Cakoni | A Qualitative Approach to Inverse Scattering Theory[END_REF] combined the LSM with iterative algorithms as a tool to expose the surface properties of obstacles from acoustic and electromagnetic data. In elastodynamics, [START_REF] Minato | Imaging and characterization of a subhorizontal non-welded interface from point source elastic scattering response[END_REF] proposed a Fourier-based approach employing the reverse-time migration and wavefield extrapolation to retrieve the heterogeneous compliance of a planar interface under the premise of high frequency and absence of evanescent waves along the interface.

In light of the above developments, the twin aim of (non-iterative) seismic imaging and interfacial characterization of heterogeneous fractures is pursued sequentially via a 3-step approach where: (1) the shape reconstruction of penetrable fractures is effected, without prior knowledge of their contact condition, via recently established GLSM [START_REF] Pourahmadian | Generalized linear sampling method for elastic-wave sensing of heterogeneous fractures[END_REF] and TS [START_REF] Pourahmadian | On the elastic-wave imaging and characterization of fractures with specific stiffness[END_REF] approaches to elastic waveform tomography of discontinuity surfaces; (2) given the reconstructed fracture geometry Γ, the fracture opening displacement (FOD) profile is recovered via a double-layer potential representation mapping the FOD to the scattered field observations, and (3) the specific stiffness profile (as given by its normal, shear, and mixed-mode components) -is resolved from the knowledge of FOD and Γ by making use of the traction boundary integral equation written for Γ. To help construct a robust inverse solution, a three-step regularization scheme is also devised to minimize the error due to: (i) compactness of the double-layer potential map used to recover the FOD; (ii) inexact geometric reconstruction of the fracture surface, and (iii) presence (if any) of areas on Γ with near-zero FOD values. The performance of the inverse solution is illustrated by a set of numerical experiments assuming seismic illumination in the resonance region.

Preliminaries

With reference to Fig. 1, consider a heterogeneous (possibly unconnected) fracture Γ ⊂ R 3 embedded in a homogeneous, isotropic elastic solid endowed with mass density ρ and Lamé parameters µ and λ. In the spirit of the linear slip model [START_REF] Schoenberg | Elastic interface waves along a fracture[END_REF] used to describe the seismic response of fractures in rock, the contact condition over Γ is given by a 3 × 3 symmetric matrix, K(ξ), of specific stiffness coefficients -synthesizing the spatially-varying nature of its rough and possibly multi-phase interface. Assuming time-harmonic seismic illumination, K is taken to be complex-valued with (K) 0 in order to allow for energy dissipation at the interface and to ensure the well-posedness of the forward scattering problem [START_REF] Pourahmadian | Generalized linear sampling method for elastic-wave sensing of heterogeneous fractures[END_REF].

Let Ω denote the unit sphere centered at the origin. For a given triplet of vectors d ∈ Ω and q p , q s ∈ R 3 such that q p d and q s ⊥ d, consider the case when the fracture is illuminated by a combination of compressional and shear plane waves

u i (ξ) = q p e ikpd•ξ + q s e iksd•ξ (1)
propagating in direction d, where k p and k s = k p (λ+2µ)/µ denote the respective wave numbers. The interaction of Γ with the incident field u i gives rise to the scattered field ũ ∈ H 1 loc (R 3 \Γ) 3 which satisfies

∇•(C : ∇ũ) + ρω 2 ũ = 0 in R 3 \Γ, n • C : ∇ũ = K(ξ) ũ -t i on Γ. (2)
Here, ω 2 = k 2 s µ/ρ is the frequency of excitation; ũ = [ũ +ũ-] is the jump in ũ across Γ, hereon referred to as the fracture opening displacement (FOD);

C = λ I 2 ⊗ I 2 + 2µ I 4 (3)
is the fourth-order isotropic elasticity tensor; I m (m = 2, 4) denotes the mth-order symmetric identity tensor; t i = n • C : ∇u i is the incident-field traction vector, and n := n -is the unit normal on Γ. For clarity it is noted that, owing fo the continuity of the incident field u i , the jump in the total field u = u i + ũ across Γ equals that of the scattered field, namely

+ n + n K(⇠)
u = ũ on Γ. (4)
On uniquely decomposing ũ into an irrotational part and a solenoidal part as ũ = ũp + ũs where ũp = 1

k 2 s -k 2 p (∆ + k 2 s )ũ, ũs = 1 k 2 p -k 2 s (∆ + k 2 p )ũ, (5)
( ξ) + O(r -2 ) as r = |ξ| → ∞, (7)
where ξ = ξ/r is the unit direction of observation, while ũ∞ p ∈ L 2 (Ω) 3 and ũ∞ s ∈ L 2 (Ω) 3 denote respectively the far-field patterns of ũp and ũs that admit [START_REF] Pourahmadian | Generalized linear sampling method for elastic-wave sensing of heterogeneous fractures[END_REF] the integral representation

ũ∞ p ( ξ) = -ik p ξ Γ λ ũ •n + 2µ n• ξ ũ • ξ e -ikp ξ•x dS x , ũ∞ s ( ξ) = -ik s ξ × Γ µ ũ × ξ (n• ξ ) + µ n× ξ ( ũ • ξ) e -iks ξ•x dS x . (8)
In this setting, we define the far-field pattern of ũ by ũ∞ := ũ∞ p ⊕ ũ∞ s .

For the purposes of seismic fracture sensing, ũ∞ is monitored over a subset Ω obs ⊆ Ω of the unit sphere.

Remark 1. The featured framework of elastodynamic scattering in R 3 is introduced for the reasons of convenience only. In general, the ensuing approach to quantitative reconstruction of the specific stiffness profile K(ξ), once the fracture geometry has been resolved, applies equally to situations when the reference elastic domain D ⊂ R 3 is semi-infinite, bounded, or heterogeneous (e.g. layered half-space). In this vein, the analysis also caters for sensing configurations entailing near-field observations of the scattered field ũ over a limited-aperture surface S obs ⊂ D.

Elastic-wave sensing of heterogeneous fractures

In what follows, the model problem in Section 2 is used as a basis to describe the inverse solution schematically shown in Fig. 2, where the sensory waveform data (in this case the far-field pattern ũ∞ ) provide an input to the 3-step approach for geometric reconstruction and interfacial characterization of subsurface discontinuities in an elastic solid, e.g. natural or hydraulic fractures in rock. In this framework,

• The fracture geometry Γ is reconstructed, irrespective of its contact condition, via a suitable approach to non-iterative seismic waveform tomography [START_REF] Pourahmadian | On the elastic-wave imaging and characterization of fractures with specific stiffness[END_REF][START_REF] Pourahmadian | Generalized linear sampling method for elastic-wave sensing of heterogeneous fractures[END_REF];

• The FOD profile, ȗ , over the reconstructed fracture support Γ is recovered from the germane boundary integral representation (double-layer potential) of the scattered field, and

• The specific stiffness profile K(ξ) -as given by its normal, shear, and mixed components -is recovered from the knowledge of ȗ and Γ.

These basic steps are elucidated in the sequel.

Seismic experiment Specific stiffness FOD profile ➠ , ˘ ➠ JȗK K(⇠) K(⇠) u i

Geometric reconstruction

Non-iterative elastic waveform tomography

➠ ˘ ũ , ˘ ũ JȗK on ˘ ũ Figure 2
. Three-step approach to non-iterative geometric reconstruction and interfacial characterization of heterogeneous fractures by elastic waves.

Geometric fracture reconstruction

The essential first ingredient toward comprehensive sensing of heterogeneous fractures is an ability to decipher the observed elastic waveforms toward geometric reconstruction of germane discontinuity surfaces. By building on the earlier works in scalar inverse scattering [START_REF] Boukari | The factorization method applied to cracks with impedance boundary conditions[END_REF][START_REF] Audibert | A generalized formulation of the linear sampling method with exact characterization of targets in terms of farfield measurements[END_REF], electromagnetism [START_REF] Park | Topological derivative strategy for one-step iteration imaging of arbitrary shaped thin, curve-like electromagnetic inclusions[END_REF], and elastic-wave sensing of impenetrable (tractionfree) discontinuities [START_REF] Bellis | Crack identification by 3D time-domain elastic or acoustic topological sensitivity[END_REF][START_REF] Bellis | Qualitative identification of cracks using 3D transient elastodynamic topological derivative: formulation and FE Comp[END_REF], recent research efforts have shown a path toward non-iterative elastic waveform tomography of penetrable fractures (irrespective of their contact condition) via the TS approach [START_REF] Pourahmadian | On the elastic-wave imaging and characterization of fractures with specific stiffness[END_REF] -as corroborated by high-frequency simulations and laboratory observations [START_REF] Tokmashev | Experimental validation of the topological sensitivity approach to elastic-wave imaging[END_REF][START_REF] Guzina | Why the high-frequency inverse scattering by topological sensitivity may work[END_REF], and the GLSM paradigm [START_REF] Pourahmadian | Generalized linear sampling method for elastic-wave sensing of heterogeneous fractures[END_REF]. To provide an example and to help establish an explicit setting for the discussion, the GLSM approach to elastic-wave imaging of heterogeneous fractures is summarized next.

3.1.1. Generalized Linear Sampling Method. For a given vector density g = g p ⊕ g s ∈ L 2 (Ω) 3 comprised of its compressional (g p ) and shear (g s ) wave components mimicking the decomposition of q in (1), consider the elastic Herglotz wave function [START_REF] Dassios | Elastic Herglotz functions[END_REF]

u i g (ξ) := Ω g p (d)e ikpd•ξ dS d ⊕ Ω g s (d)e iksd•ξ dS d , ξ ∈ R 3 , (10)
and define the far-field operator F :

L 2 (Ω) 3 → L 2 (Ω) 3 by F (g) := ũ∞ g , (11)
where ũ∞ g is the far-field pattern (9) of ũ ∈ H 1 loc (R 3 \Γ) 3 solving ( 2) and ( 6) with data

t i = n • C : ∇u i g on Γ.
In this setting, the idea of the GLSM is to construct a nearby solution of the (ill-posed) far-field equation

F g Φ ∞ L , (12)
where Φ ∞ L is the far-field pattern of a test radiating field, generated by a small trial fracture L ⊂ R 3 with prescribed FOD a ∈ H1/2 (L) 3 (see [START_REF] Pourahmadian | Generalized linear sampling method for elastic-wave sensing of heterogeneous fractures[END_REF] for details). Without loss of generality, L can be taken as a vanishing penny-shaped fracture at z ∈ R 3 with normal n ∈ Ω and constant (mode I) FOD profile a ∝ n, in which case (8) yields

Φ ∞ L ( ξ) = -ik p ξ λ + 2µ(n • ξ) 2 e -ikp ξ•z ⊕ 2iµk s ξ × (n × ξ)(n • ξ) e -iks ξ•z . ( 13
)
With reference to [START_REF] Boukari | The factorization method applied to cracks with impedance boundary conditions[END_REF], a heterogeneous fracture Γ with arbitrary (linear) contact condition can then be reconstructed [START_REF] Pourahmadian | Generalized linear sampling method for elastic-wave sensing of heterogeneous fractures[END_REF] from the support of the GLSM characteristic function

C L(z, n) = (F ) 1 2 g L α,δ 2 + δ g L α,δ 2 -1/2 , z ∈ R 3 , n ∈ Ω ( 14
)
where

F := 1 2 |(F + F * )| + 1 2i (F -F * ) (15)
is a self-adjoint surrogate for F [START_REF] Kirsch | The Factorization Method for Inverse Problems[END_REF], and g L α,δ minimizes the penalized least-squares functional

J δ α (g; Φ ∞ L ) : = F g -Φ ∞ L 2 + α (F ) 1 2 g 2 + δ g 2 , g ∈ L 2 (Ω) 3 (16)
that features an absolute measure, δ > 0, of perturbation in the data (given by the far-field operator F ) and a penalization parameter α = α(δ) > 0 as in [START_REF] Audibert | A generalized formulation of the linear sampling method with exact characterization of targets in terms of farfield measurements[END_REF].

Remark 2. For generality, it should be noted that geometric fracture reconstruction can also be accomplished by other application-specific inversion schemes, e.g. microseismic imaging -which deploys travel time inversion or migration of seismic events generated during the fracturing process [START_REF] Maxwell | Microseismic Imaging of Hydraulic Fracturing[END_REF]. The featured TS and GLSM techniques are, however, part of a more general waveform inversion platform which is (i) computationally efficient due to its non-iterative nature; (ii) self-contained for it uses the common set of sensory data for both geometric reconstruction and interfacial characterization, and (iii) robust by providing significant flexibility in terms of the sensing configuration and requiring no a priori knowledge on the fracture geometry nor its contact condition.

Inversion of the FOD profile

Given a geometric reconstruction of the fracture surface Γ -obtained as described in Section 3.1, a double-layer (elastodynamic) potential representation of the scattered field [START_REF] Bonnet | Boundary Integral Equation Methods for Solids and Fluids[END_REF] serves as a map M : H1/2 ( Γ) 3 → L 2 (Ω obs ) 3 relating the sought FOD, ȗ , to the sensory data ũ∞ , namely

M ȗ ( ξ) = ũ∞ ( ξ), M ȗ ( ξ) := Γ Σ ∞ ( ξ, x) : ȗ (x) ⊗ n(x) dS x , ξ ∈ Ω obs (17)
where n = nis the unit normal on Γ, and Σ ∞ is the far-field pattern as |ξ| → ∞ of the elastodynamic fundamental stress tensor Σ(ξ, x) [START_REF] Achenbach | Reciprocity in Elastodynamics[END_REF]. In terms of dyadic notation and Einstein summation convention, the latter quantity can be written as

Σ(ξ, x) = Σ ij (ξ, x) e ⊗ e i ⊗ e j , ξ, x ∈ R 3 , ξ = x, i, j, ∈ {1, 2, 3} (18)
where Σ ij = Σ ji signify the components of the Cauchy stress tensor at ξ due to point force e (i.e. the unit vector in the th coordinate direction) acting at x. For completeness, it is noted (roughly speaking) that the Sobolev space H1/2 ( Γ) denotes the space of functions in H 1/2 ( Γ) whose extension by zero to a larger Lipshitz surface, Λ ⊃ Γ, gives functions that belong to H 1/2 (Λ) [START_REF] Pourahmadian | Generalized linear sampling method for elastic-wave sensing of heterogeneous fractures[END_REF]. The space H1/2 ( Γ) is also the dual of H -1/2 ( Γ).

Lemma 3.1. Linear operator M: H1/2 ( Γ) 3 → L 2 (Ω obs ) 3 introduced in ( 17) is compact and injective (see also Lemma 5.2 in [START_REF] Pourahmadian | Generalized linear sampling method for elastic-wave sensing of heterogeneous fractures[END_REF]).

Proof. Thanks to the explicit asymptotic expansion of Σ as |ξ| → ∞ [e.g. 1] and [START_REF] Bellis | Crack identification by 3D time-domain elastic or acoustic topological sensitivity[END_REF], one finds that

Σ ,∞ ij ( ξ, x) = ik p (2µ ξi ξj + λδ ij ) ξ e -ikp ξ•x ⊕ ik s µ(δ i ξj + δ j ξi -2 ξi ξj ξ ) e -iks ξ•x , i, j, ∈ {1, 2, 3} (19)
which are holomorphic over Ω×R 3 . As a result, integral operator M has a smooth kernel and is therefore compact from H1/2 ( Γ) 3 to L 2 (Ω obs ) 3 . To establish the injectivity of M , one may observe as in [START_REF] Pourahmadian | Generalized linear sampling method for elastic-wave sensing of heterogeneous fractures[END_REF] that the vanishing far-field pattern [START_REF] Bellis | Qualitative identification of cracks using 3D transient elastodynamic topological derivative: formulation and FE Comp[END_REF] requires that ũ = 0 in R 3 \ Γ by the Rellich Lemma, and consequently that ȗ = 0 thanks to the unique continuation principle and the fundamental property of double-layer potentials by which ȗ = ũ on Γ.

Remark 3. In the sequel, it is assumed that M: H1/2 ( Γ) 3 → L 2 (Ω obs ) 3 has a dense range, a hypothesis that is supported by Lemma 5.3 in [START_REF] Pourahmadian | Generalized linear sampling method for elastic-wave sensing of heterogeneous fractures[END_REF] giving sufficient conditions on the fracture geometry Γ and excitation frequency ω for the range denseness of (17).

Remark 4. In physical terms, the compactness of M is reflected by the presence of interface waves [START_REF] Pyrak-Nolte | Elastic interface waves along a fracture[END_REF] on Γ. In theory these local waves, propagating along the surfaces of discontinuity in an elastic medium, are characterized by an exponential decay [START_REF] Stoneley | Elastic waves at the surface of separation of two solids[END_REF][START_REF] Pyrak-Nolte | Elastic interface waves along a fracture[END_REF] with normal distance to the interface. Despite the apparent leakage of such surface-wave energy into the exterior due to the curvature of Γ (if any) and the interaction of interfacial waves with the fracture edge ∂ Γ, these local FOD modes may have only a marginal fingerprint in the remote sensory data that warrants a custom treatment in the inversion scheme.

Remark 5. For a generic sensing configuration entailing (i) semi-infinite or bounded reference domain D ⊂ R 3 and (ii) near-field measurements of the scattered field, ũ∞ , Ω obs , and Σ ∞ in ( 17) are suitably replaced by ũ, S obs ⊂ D, and Σ D -the germane elastodynamic Green's function. Assuming a physical separation between the fracture surface and the sensing grid, the relevant double-layer potential map M D can likewise be shown (on a case-by-case basis) to be compact thanks to the smoothness of Σ D .

Discretization. Taking Ω obs as a union of discrete sensing directions -as may be the case in a physical experiment, and allowing the reconstructed fracture surface Γ to be arbitrarily complex, a discrete version of M ȗ may be obtained via suitable discretization of Γ and FOD in terms of surface i.e. boundary elements [START_REF] Bonnet | Boundary Integral Equation Methods for Solids and Fluids[END_REF][START_REF] Pourahmadian | On the elastic-wave imaging and characterization of fractures with specific stiffness[END_REF]. As a result, ( 17) can be recast via the collocation method as

M ȗ = ũ∞ , ( 20
)
where M is a 3N obs × 3N nds coefficient matrix; ȗ = ȗ is a vector sampling the FOD at N nds geometric nodes over Γ (see [START_REF] Pourahmadian | Generalized linear sampling method for elastic-wave sensing of heterogeneous fractures[END_REF]Fig. 11]), and ũ∞ collects the far-field pattern (9) of the scattered field in N obs sensing directions over Ω obs .

Solution. In this setting, the idea is to have N nds < N obs and to compute the FOD profile on Γ from the sensory data ũ∞ by solving an overdetermined linear system [START_REF] Colton | A simple method for solving inverse scattering problems in the resonance region[END_REF]. As a result, the spatial resolution of the FOD reconstruction will be inherently limited (in addition to other factors) by the number of observation directions on S obs . Also note that (20) can be solved anew for each available direction d p ∈ Ω (p = 1, . . . P ) of the incident plane wavefield (1), where M is invariant i.e. source-independent.

Regularization. A critical point in recovering the FOD -which subsequently affects the inversion of the specific stiffness -is that ( 17) is ill-posed due to the compactness of M .

In the context of discrete statement [START_REF] Colton | A simple method for solving inverse scattering problems in the resonance region[END_REF], M may accordingly contain unacceptably small singular values, whose number will (for a given fracture geometry and contact characteristics) depend on the properties of the incident field u i . To provide a physical interpretation of such behavior, consider the singular value decomposition (SVD) of M as

M = U M Λ M V * M , (21)
where U M (resp. V M ) collects the left (resp. right) eigenvectors of M, and Λ M contains its singular values. In this setting, the right eigenvectors V M,q (q = 1, Q) corresponding to the Q smallest (unacceptably small) singular values of M, Λ q (which reflect the compactness of M ), can be affiliated with the interface wave modes [START_REF] Pyrak-Nolte | Elastic interface waves along a fracture[END_REF] on Γ, see Remark 4. In the context of ( 20), the ill-posedness of ( 17) can be dealt with in a standard way via e.g. truncated SVD or Tikhonov regularization aided by the Morozov discrepancy principle [START_REF] Kirsch | An Introduction to the Mathematical Theory of Inverse Problems[END_REF](which takes the penalty parameter to be commensurate with the level of noise in the data). It will be shown later, however, that the approximating effect of such regularization can be mitigated assuming the availability of sensory data (ũ ∞ ) for multiple incident fields -which is a customary premise for most inverse scattering solutions.

Inversion of the specific stiffness profile

The last step in the proposed inverse scheme is to substitute the identified FOD, ȗ , into the fracture's boundary condition on Γ, and to solve thus-obtained equation for the specific stiffness profile K(ξ). In this vein, the true contact condition (2) on Γ is recast over the reconstructed fracture surface Γ as

K(ξ) ȗ = t + ti , ξ ∈ Γ, (22)
where ti = n•C :∇u i is the incident-field traction on Γ, while t is the scattered-field traction on Γ expressed in terms of ȗ = ȗ by invoking the traction boundary integral equation (TBIE) [START_REF] Pourahmadian | On the elastic-wave imaging and characterization of fractures with specific stiffness[END_REF][START_REF] Bonnet | Boundary Integral Equation Methods for Solids and Fluids[END_REF] as a map T : H1/2 ( Γ) → H -1/2 ( Γ) such that

t(ξ) = T ȗ := -n(ξ) • C : - Γ Σ(ξ, x) : D x ȗ (x) dS x + ρω 2 n(ξ) • C : Γ U (ξ, x) • ȗ ⊗ n (x) dS x , ξ ∈ Γ, (23)
where U = U i (ξ, x) e i ⊗ e and Σ = Σ ij (ξ, x) e i ⊗ e j ⊗ e denote respectively the elastodynamic displacement and stress fundamental solution at ξ ∈ R 3 due to point force acting at x ∈ R 3 (see [46, Appendix A], [START_REF] Guzina | On the analysis of wave motions in a multi-layered solid[END_REF][START_REF] Bonnet | Boundary Integral Equation Methods for Solids and Fluids[END_REF]); -signifies the Cauchy-principal-value integral, and D x is the tangential differential operator given by

D x (f ) = D kl (f m ) e l ⊗ e m ⊗ e k , D kl (f m ) = nk f m,l -nl f m,k , (24)
with nk = nk (x) and f m,k = ∂f m /∂x k in the global coordinate frame {e 1 , e 2 , e 3 }.

3.3.1. Indirect solution approach. Given ȗ on Γ, the scattered-field traction t(ξ) in ( 22) can be computed without any reference to K(ξ). This is accomplished by imposing ȗ = ȗ as a Dirichlet boundary condition on Γ, and solving the resulting (exterior) boundary value problem in the reference domain (R 3 in the present study, or D ⊃ Γ in a more general case). With the knowledge of t and ti on Γ at hand, K(ξ) can then in principle be solved from [START_REF] Dassios | Elastic Herglotz functions[END_REF]. In particular:

• If K is assumed to be diagonal in the fracture's local coordinates, namely K := diag(κ n , κs 1 , κs 2 ), one arrives via ( 22) at the uncoupled system κp (ξ) ȗp (ξ) = tp (ξ) + ti p (ξ), p ∈ {n, s 1 , s 2 }, ξ ∈ Γ for the diagonal entries of K(ξ) (no summation implied), which is solvable from the sensory data for a single incident field -provided that ȗ (ξ) does not vanish at ξ.

• In situations where K is taken to be fully populated -implying dilatant contact behavior (i.e. coupling between the normal and shear resistance), one obtains the coupled system Kpq (ξ) ȗq (ξ) = tp (ξ) + ti p (ξ), p, q ∈ {n, s 1 , s 2 }, ξ ∈ Γ which makes use of the Einstein summation convention. In this case there are (recalling the symmetry of K) six unknowns at given ξ ∈ Γ, requiring the availability of sensory data for at least two incident fields.

This approach has an advantage that it does not require the Green's function for the reference domain, which is convenient when dealing with fracture sensing inside finite bodies, or using numerical solution techniques such as finite element or finite difference methods.

Direct solution approach.

One apparent drawback of the foregoing scheme is that it requires an additional forward solution for each incident field, as required to compute the scattered-field traction t(ξ). In situations where the Green's function for the reference domain is available (as in the present case), this problem can be partly circumvented by combining ( 22) and [START_REF] Doqueta | Experimental and numerical study of crack healing in a nuclear glass[END_REF] as

K(ξ) ȗ = T ȗ + ti , ξ ∈ Γ, (25)
where the right-hand side can be computed as follows.

Discretization. By deploying the collocation method as examined in Section 3.2, a discretized version of the integro-differential operator T ȗ in [START_REF] Doqueta | Experimental and numerical study of crack healing in a nuclear glass[END_REF] can be computed as T ȗ , where T is a 3N col × 3N nds coefficient matrix; N col is the number of points on Γ where (25) is collocated, and ȗ is a vector sampling the FOD at N nds geometric nodes over Γ as before.

As examined in the sequel, N col -carrying the parametrization of K(ξ) -can be either smaller or larger than N nds . To meet the computational C 1 smoothness requirement on FOD in (23) (due to appearance of the tangential differential operator D x ) while allowing ȗ on Γ to be parametrized via standard C 0 (boundary element) interpolation, the collocation points for solving [START_REF] French | Displacement and dynamic weakening processes in smectite-rich gouge from the central deforming zone of the San Andreas fault[END_REF] are conveniently placed in the interior of the boundary elements used to represent the geometry and kinematics of the reconstructed fracture surface (see [START_REF] Pourahmadian | On the elastic-wave imaging and characterization of fractures with specific stiffness[END_REF], Appendix B for details). However, since the collocation points belong to the fracture surface, the first integral on the right-hand-side of ( 23) remains singular and must be properly regularized [START_REF] Bonnet | Boundary Integral Equation Methods for Solids and Fluids[END_REF][START_REF] Pourahmadian | On the elastic-wave imaging and characterization of fractures with specific stiffness[END_REF]. In this setting, a discretized statement of the contact condition ( 25)

reads K ȗ = T ȗ + ti , ( 26
)
where K is a 3N col × 3N nds block-diagonal matrix of specific stiffness coefficients, and ti is the free-field traction vector sampled at N col collocation points over Γ.

Solution.

To solve [START_REF] Gu | Progressive failure of lower San Fernando dam[END_REF] for the entries of K, one may consider the following situations.

• Assuming K := diag(κ n , κs 1 , κs 2 ), K in ( 26) becomes diagonal and the specific stiffness profile can be evaluated directly at collocation points. In particular, the left-hand side of ( 26) can be recast as

K ȗ := Ȃ ȗ k, ( 27
)
where Ȃ ȗ is a 3N col × 3N col diagonal coefficient matrix given by the reconstructed FOD profile, and k is a 3N col -vector sampling the (normal and shear) specific stiffness coefficients, namely {κ n , κs 1 , κs 2 } at N col collocation nodes over Γ. Assuming m internal collocation points per boundary element in the evaluation of ( 26), one accordingly has N col = mN el N nds , where N el denotes the number of boundary elements. In situations when Ȃ ȗ is free of zero or near-zero diagonal elements (the contrary case will be discussed shortly), the substitution of ( 27) into [START_REF] Gu | Progressive failure of lower San Fernando dam[END_REF] immediately yields a stable solution for the specific stiffness profile, k, on Γ.

• In situations where K is taken as fully populated whereby the number of unknowns per contact point is six, the specific stiffness distribution is parametrized using the previously established geometric interpolation -so that the unknowns are evaluated at N nds geometric nodes. In this setting, one has

K ȗ := B ȗ k, ( 28
)
where B ȗ is a 3N col × 6N nds coefficient matrix given by the recovered FOD profile, and k is a 6N nds -vector collecting the six entries of the symmetric stiffness matrix K sampled at N nds geometric nodes over Γ. Accordingly, by taking m internal collocation points per element so that N col = mN el 2N nds , one obtains the sufficient number of equations to solve [START_REF] Gu | Progressive failure of lower San Fernando dam[END_REF] for k. Similar to the previous case, this yields a stable solution provided that B ȗ is free of zero or near-zero singular values, which physically implies that no geometric nodes be affiliated with permanent contact (zero FOD) on Γ. For completeness, this and other likely sources of ill-posedness jeopardizing the evaluation K(ξ) are addressed next.

Regularization

To help construct a robust solution to [START_REF] Gu | Progressive failure of lower San Fernando dam[END_REF], a three-step regularization scheme is devised to suppress the error and instabilities due to: (1) compactness of the double-layer potential M ȗ in [START_REF] Cakoni | Inverse Scattering Theory and Transmission Eigenvalues[END_REF], physically stemming from the presence of interface waves on Γ;

(2) tangential differentiation of the recovered FOD along inexact fracture surface Γ, and (3) presence (if any) of areas on Γ with near-zero FOD values. These three steps of regularization are elucidated in the sequel. To aid the discussion, it is convenient to recall [START_REF] Cook | Natural joints in rock: Mechanical, hydraulic and seismic behaviour and properties under normal stress[END_REF] and to similarly introduce the SVD of matrix T discretizing the integro-differential operator T in [START_REF] Doqueta | Experimental and numerical study of crack healing in a nuclear glass[END_REF] as

T = U T T V * T , (29)
where U T (resp. V T ) collects the left (resp. right) eigenvectors of T, and Λ T contains its singular values.

3.4.1. Suppression of interface waves on Γ. Recalling [START_REF] Colton | A simple method for solving inverse scattering problems in the resonance region[END_REF] as the basis for FOD recovery and assuming the availability of sensory data for multiple incident fields, (p = 1, . . . P ), the idea is to synthetically recombine the available scattered-field data ũ∞,p in order to minimize the participation of interface waves in the associated FOD profile ȗ . To do so, recall that M is excitation-independent and let U M,q (q = 1, 2, . . . Q) denote the left eigenvectors affiliated with the Q smallest singular values of M that are suppressed by the regularization of [START_REF] Colton | A simple method for solving inverse scattering problems in the resonance region[END_REF]. On conveniently rewriting the latter equation for the pth incident field as Λ M V * M ȗ = U * M ũ∞,p , the idea is to solve

U * M,q P p=1 g p ũ∞,p = 0, q = 1, Q, (30)
for the (synthetic) source density g = (g 1 , g 2 , . . . , g P ), designed to minimize the participation of (unrecoverable) interface waves in the FOD on Γ -without prior knowledge of K(ξ).

In this way, the solution error due to regularization of ( 20) via e.g. truncated SVD is suppressed, while maintaining the stability of the solution by deploying the sensory data ũ∞ = P p=1 g p ũ∞,p as the basis for solving [START_REF] Colton | A simple method for solving inverse scattering problems in the resonance region[END_REF]. In this setting, one proceeds with solving (26) for the specific stiffness profile k by setting ti = P p=1 g p ti,p where tf p i,p is the free-field traction on Γ due to pth incident field.

Remark 6. Assuming P < Q in solving [START_REF] Guzina | Why the high-frequency inverse scattering by topological sensitivity may work[END_REF], the " strength" g p of each incident field can be computed via least squares. When P > Q, on the other hand, one can take Q out of P available sets of sensory data to construct an even-determined problem.

3.4.2. Regularization of T. In general the linear operator T : H1/2 ( Γ) → H -1/2 ( Γ) introduced in [START_REF] Doqueta | Experimental and numerical study of crack healing in a nuclear glass[END_REF], that is discretized as T, is constructed on the basis of the recovered fracture support Γ and entails tangential differentiation according to [START_REF] Eberhardt | Numerical analysis of initiation and progressive failure in natural rock slopes -the 1991 Randa rockslide[END_REF]. As a result, T may contain large singular values due to the approximate nature of Γ and its roughness, leading to the amplification of small errors in ȗ while computing the right-hand-side of [START_REF] Gu | Progressive failure of lower San Fernando dam[END_REF]. To ensure a stable solution in terms of k, consider the minimal subspace of R 3N nds , spanned by the first N right eigenvectors V T,n (n = 1, 2, . . . N ) of T, that contributes to the construction of ȗ below a designated error δ (a measure of error in FOD reconstruction), namely ȗ -

N n=1 ( ȗ •V T,n )V T,n < δ. (31)
In this setting, the first term on the right-hand side in ( 26) can be regularized as

T ȗ N n=1 T ȗ •U T,n U T,n , (32)
where U T,n (n = 1, 2, . . . N ) are the left eigenvectors of T.

3.4.3. Treatment of vanishing FOD on Γ. Solving [START_REF] Gu | Progressive failure of lower San Fernando dam[END_REF] for the interfacial stiffness profile k, one may find that the coefficient matrix Ȃ ȗ in [START_REF] Guzina | Topological derivative for the inverse scattering of elastic waves[END_REF] or B ȗ in ( 28) is singular. To interpret this situation, one may recall [START_REF] Guzina | Topological derivative for the inverse scattering of elastic waves[END_REF] and observe that the near-zero eigenvalues of Ȃ ȗ are affiliated with those collocation points ξ ∈ Γ where ȗ (ξ) → 0. One such example can be constructed by considering a penny-shaped fracture with locally orthotropic K(ξ), subjected to a pair of antiplane shear incident waves (1) whose directions of incidence are symmetric with respect to the fracture plane, and whose polarization is parallel to the fracture plane; in this case the FOD can be shown to vanish along the entire fracture. To ensure a stable solution in situations where the reconstructed FOD vanishes over a subset of germane collocation points, (26) can be solved via suitable regularization e.g. by invoking Tikhonov regularization or truncated SVD. In this case, however, additional incident fields may be needed to quantitatively resolve K(ξ) at those locations.

Numerical implementation and results

In light of the existing elastodynamic simulations [START_REF] Pourahmadian | On the elastic-wave imaging and characterization of fractures with specific stiffness[END_REF][START_REF] Pourahmadian | Generalized linear sampling method for elastic-wave sensing of heterogeneous fractures[END_REF] -specifically focused on the geometric reconstruction of heterogeneous fractures via TS and GLSM, this section examines the effectiveness of the proposed 3-step approach (see Fig. 2) with a particular emphasis on the inversion of the specific stiffness profile. The testing configuration is illustrated in Fig. 3, where a cylindrical fracture Γ of width L = 0.7, arclength = 0.55, and radius R = 0.35 is embedded in a linear, isotropic and homogeneous elastic solid with shear and compressional wave speeds c s = 1 and c p = 2.08, respectively. As shown in Fig. 4, the fracture is endowed with two sets of (diagonal and orthotropic) interfacial stiffness profiles K(ξ), namely: (1) a piecewise-constant "zebra" distribution in both shear and normal directions, and (2) a "cheetah" pattern defined in the fracture's local coordinates. By making reference to the local orthonormal basis (n, e 1 , e 2 ) on Γ, one may conveniently write the specific stiffness matrix for both configurations as

K(ξ) = κ n (ξ)n ⊗ n + κ s (ξ) e 1 ⊗ e 1 + e 2 ⊗ e 2 }, ξ ∈ Γ. ( 33
)
As can be seen from Fig. 4, the pair (κ n , κ s ) is assumed to be complex-valued, signifying energy dissipation (due to e.g. friction) at the interface. The fracture is illuminated by a set of incident plane waves (1), taking k s such that the ratio between the shear wavelength and the arclength of Γ is λ s / = 0.7. Thus-induced scattered field is then measured in terms of the far-field pattern, ũ∞ , given by ( 7)-( 9). The spatial density of sensory data, for both illumination and sensing purposes, is given by N θ ×N φ = 25×12 directions given by the polar (θ j , j = 1, . . . N θ ) and azimuthal (φ k , k = 1, . . . N φ ) angle values spanning the unit sphere Ω. In what follows, all forward simulations are performed by an elastodynamic boundary integral equation code [START_REF] Pak | Seismic soil-tructure interaction analysis by direct boundary element methods[END_REF][START_REF] Pourahmadian | On the elastic-wave imaging and characterization of fractures with specific stiffness[END_REF] and polluted by 5% random noise in terms of the observed far-field pattern. Geometric reconstruction. For each interface scenario the fracture support Γ is recovered, by way of the GLSM framework described in Section 3.1.1, from the knowledge of the far-field operator [START_REF] Bonnet | Boundary Integral Equation Methods for Solids and Fluids[END_REF] as shown in Fig. 5 (see [START_REF] Pourahmadian | Generalized linear sampling method for elastic-wave sensing of heterogeneous fractures[END_REF] for details).

FOD recovery. Given Γ, one may construct a "double-layer potential" operator M in (20) -mapping the sought FOD, ȗ , to the sensory data ũ∞ by way of [START_REF] Cakoni | Inverse Scattering Theory and Transmission Eigenvalues[END_REF]. As shown by Lemma 3.1,[START_REF] Colton | A simple method for solving inverse scattering problems in the resonance region[END_REF] represents an ill-posed problem due to the compactness of M which is reflected by the appearance of interfacial waves [START_REF] Pyrak-Nolte | Elastic interface waves along a fracture[END_REF] on Γ -which are known to decay exponentially away from the fracture surface. This is illustrated in Fig. 6 where the "true" FOD over Γ -obtained by simulating the forward problem (2) due to a single incident S-wave, is decomposed into two components, namely: (i) the retrievable part ȗ computed by solving (20) over Γ via Tikhonov regularization endowed with the Morozov discrepancy principle (5% random noise), and (ii) the residual part uȗ , comprised of interfacial waves, obtained by subtracting the reconstructed displacement jump from the "true" FOD. Note that M used to recover FOD in Fig. 6 is intentionally constructed on the basis of the "true" fracture geometry Γ, so that the computed residual part is not polluted by errors due to geometric fracture reconstruction. With such remark, Fig. 7 compares the "true" FOD over Γ, u -induced by the interaction of a single incident P-wave with the "zebra" interface on Γ, with the recovered displacement jumps, ȗ according to [START_REF] Colton | A simple method for solving inverse scattering problems in the resonance region[END_REF], over Γ and Γ. For a robust inversion of the interfacial stiffness, however, one should reconstruct the FOD profile ˘ ˘ zebra cheetah from ũ∞ by making use of the first regularization step in Section 3.3 -where the observed wavefield data from multiple incident fields are synthetically recombined to minimize, and possibly eliminate, the participation of interface waves in the FOD on Γ. With reference to [START_REF] Guzina | Why the high-frequency inverse scattering by topological sensitivity may work[END_REF], this is accomplished by selecting Q as the number of singular values of M that are smaller than 15% of its largest singular value. The result is is shown in Fig. 8, which compares the "true" FOD over Γ -due to recombined incident fields -with the corresponding reconstruction ȗ over Γ. Reconstruction of the specific stiffness. By virtue of (33) and the developments in Section 3.3, one finds the regularized statement of (26) to read

Ȃ ȗ k = N n=1 T ȗ •U T,n U T,n + ti . (34)
Here, Ȃ ȗ a diagonal coefficient matrix given by the reconstructed FOD at collocation points on Γ; ti collects the incident-field tractions on Γ, and U T,n (n = 1, 2, . . . N ) are the left eigenvectors of T, truncated as in [START_REF] Hedayat | Precursors to the shear failure of rock discontinuities[END_REF] with δ = 0.001 to minimize the effect of imperfect geometric reconstruction. As examined in Section 3.4.3, ( 34) is solved via Tikhonov regularization and Morozov discrepancy principle applied to the noise level of 5%. In this setting, Fig. 9 shows the reconstructed zebra pattern over Γ -assuming prior knowledge of the "true" fracture geometry Γ, while Fig. 10 (resp. Fig. 11) compares the "true" zebra (resp. cheetah) K-distributions over Γ with their recovered K-counterparts on Γ. As can be seen from both displays, the fidelity of specific stiffness reconstruction is rather remarkable given (i) no prior information about the fracture geometry nor its contact condition, and (ii) multiple steps of regularization used to accomplish the sequential geometric reconstruction and interfacial characterization.

Acknowledgments

The second author kindly acknowledges the support provided by the National Science Foundation (Grant #1536110) and the University of Minnesota Supercomputing Institute.

Summary

In this study, a three-step inverse solution is proposed for the geometric reconstruction and interfacial characterization of heterogeneous fractures by seismic waves. In the approach the fracture surface is (a) allowed to be unconnected, and (b) endowed with spatially- FOD recovery assuming prior knowledge of the "true" fracture geometry Γ: shear component of ȗ (along the width of Γ), zebra interface scenario, single incident wave. The "true" FOD, u (top row) consists of (i) the trace ȗ of the propagating waves (middle row) -whose non-trivial fingerprint in the farfield data allows for their robust reconstruction, and (ii) the residual part uȗ (bottom row), given by the trace of local modes i.e. interface waves -whose vanishing signature in the far-field data prevents their recovery. dependent (linear) contact condition via a 3 × 3 complex-valued matrix of specific stiffness coefficients, which is capable of representing the effects of dilatancy, energy dissipation, and anisotropy in the fracture's contact law. In the first stage of the sensing scheme, the fracture surface is reconstructed with the aid of recently established techniques for elastic waveform tomography of heterogeneous fractures (namely the methods of topological sensitivity and generalized linear sampling) which operate irrespective of the fracture's (unknown) contact condition. With such result in place, the fracture's opening displacement (FOD) is, for given sensory data, recovered over the reconstructed fracture surface Γ via a double-layer potential representation mapping the FOD to the scattered field in the exterior domain. In the third and last step of the proposed scheme the specific stiffness profile (as given by its normal, shear, and mixed-mode components) -is recovered from the knowledge of FOD and Γ by making use of the traction boundary integral equation written for Γ. To help construct a robust inverse solution, a three-step regularization scheme is also devised to minimize the error due to: (i) compactness of the double-layer potential map used to recover the FOD, physically manifested by the presence of interfacial waves on Γ; (ii) inexact geometric reconstruction, and (iii) possible presence of areas on Γ with near-zero FOD values. The numerical results, assuming no prior knowledge of the fracture geometry nor its contact condition, suggest a remarkable performance of the staggered inverse solution and its potential toward quantitative recovery of both the elastic and dissipative contact characteristics of heterogeneous fractures.

dFigure 1 .

 1 Figure 1. Fracture Γ ⊂ R 3 endowed with a heterogeneous distribution of specific interfacial stiffness K(ξ) is illuminated by a set of plane (compressional and shear) waves propagation in direction d.

Figure 3 .

 3 Figure 3. Elastodynamic sensing configuration featuring a heterogeneous cylindrical fracture.

Figure 4 .

 4 Figure 4. Two interface scenarios for the cylindrical fracture: zebra (left) and cheetah (right) patterns representing the distribution of complex-valued specific stiffness coefficients k s and k n in (33).

Figure 5 .

 5 Figure 5. Geometric reconstruction by way of GLSM: "true" geometry Γ vs. the reconstructed fracture support Γ, obtained for the zebra and cheetah interface scenarios.

  Figure 6.

Figure 7 .Figure 8 .Figure 9 .

 789 Figure 7. FOD recovery from the far-field data collected for a single incident P-wave, normal component of ȗ , zebra interface scenario: "true" FOD on Γ as given by the forward model (left); recovered FOD over the true fracture geometry Γ (middle); recovered FOD over the reconstructed fracture geometry Γ (right).

Figure 10 .Figure 11 .

 1011 Figure 10. Recovered distribution of the specific stiffness, over the reconstructed fracture surface Γ, from the far-field patterns collected for multiple incident waves: "true" zebra pattern K(ξ) over Γ (left) vs. the recovered distribution K(ξ) over Γ (right).