Anderson localization and reentrant delocalization of tensorial elastic waves in two-dimensional fractured media
AAnderson localization and reentrant delocalization of tensorial elastic waves intwo-dimensional fractured media
Qinghua Lei ∗ and Didier Sornette
1, 2, 3, 4 Department of Earth Sciences, ETH Zurich, Zurich, Switzerland Department of Management, Technology and Economics, ETH Zurich, Zurich, Switzerland Department of Physics, ETH Zurich, Zurich, Switzerland Institute of Risk Analysis, Prediction and Management, Academy for Advanced Interdisciplinary Studies,Southern University of Science and Technology, Shenzhen, China (Dated: February 4, 2021)We study two-dimensional tensorial elastic wave transport in densely fractured media and docu-ment transitions from propagation to diffusion and to localization/delocalization. For large fracturestiffness, waves are propagative at the scale of the system. For small stiffness, multiple scatteringprevails, such that waves are diffusive in disconnected fracture networks, and localized in connectedones with a strong multifractality of the intensity field. A reentrant delocalization is found in well-connected networks due to energy leakage via evanescent waves and cascades of mode conversion.
Understanding elastic wave transport in crustal rocksis crucial for earthquake prediction and geophysical ex-ploration [1]. Fractures, as an important type of crustalheterogeneities [2, 3], can strongly affect the wavefieldbehavior [4–6]. When a wave packet enters an inhomo-geneous system, it can exhibit various regimes such aspropagation, diffusion, and localization [7–9]. The waveis propagative at scales less than the mean free pathlength, beyond which it becomes diffusive if the disorderof the medium is not too strong [10]. If the disorder isstrong enough so that the Ioffe-Regel criterion [11] is met,the wavefield enters a regime without effective transport,known as Anderson localization [12], as a result of verystrong interferences between multiply scattered waves bythe disordered medium acting as effective random mirrorcavities [13, 14]. This transition was originally predictedfor electron wave functions to explain the metal-insulatortransition in disordered crystals, but was later realized asa generic phenomenon for all wave types [15, 16]. Ander-son localization has been experimentally demonstratedfor microwaves [17], light waves [18, 19], acoustic waves[20, 21], and matter waves [22]. It is speculated that An-derson localization of elastic waves may also happen inthe Earth’s crust [16], as implied by the field evidenceof weak localization [23] and coda localization [24, 25].However, the conditions under which Anderson localiza-tion occurs in fractured crustal rocks remain unknown.In this Letter, we study the transport and localizationof tensorial elastic waves in fractured media based on two-dimensional (2D) numerical simulations. We focus on theheterogeneity effect induced by fractures and limit ourscope to pure elastic problems. We capture complex wavephenomena in fracture systems arising collectively fromthe transmission, reflection, refraction, and diffraction bynumerous individual fractures. We show a full spectrumof wave energy transport transitioning from propagationto diffusion and to localization/delocalization. ∗ [email protected] Model approach and setup. —The numerical model usesthe finite element method to solve the 2D time-domainelastodynamic equation that governs tensorial waves inlinear elastic solids under a plane-strain condition [26].Fractures are represented as discrete linear elastic in-terfaces of vanishing thickness and non-welded contactbased on the displacement discontinuity method [27],such that across each fracture, stress is continuous butdisplacement is discontinuous controlled by the fracturestiffness [4]. More details of our model including the gov-erning equations are in the Supplemental Material [28].We generate synthetic fracture networks embedded inan otherwise isotropic and homogeneous matrix within adomain Ω of size 2 L × L . The matrix has a density ρ ,and P and S wave velocities V and V (cid:48) , respectively, with V /V (cid:48) = 1 .
73. The location and orientation of fracturesare purely random, whereas their length is constant as2 l = L/
5. The fracture network connectivity is quanti-fied by the percolation parameter p = N l /L [29], where N is the number of fractures. Fracture system transformsfrom a disconnected to a connected state as p exceeds thepercolation threshold p c ≈ p values, 2, 6, and 10, corresponding to a non-connected,critically-connected, and well-connected network, respec-tively. The network of a higher p is constructed by addingfractures into that of a lower p . For each p value, 100 re-alizations are generated for ensemble averaging.The fracture networks are placed in a cell surroundedby absorbing layers [30] that suppress unwanted reflec-tions at the model boundary so as to mimic an un-bounded space. We define a Cartesian coordinate sys-tem with the origin at the domain center and two axes(x and y) orthogonal to the boundary. A point sourceis located at the origin to excite purely circular P wavesby applying a force-type sinusoid signal with the wave-length λ = L/
10. The displacement amplitude of theincident wave at the source is A ( A /L → k is explored with the normalized stiffness a r X i v : . [ phy s i c s . g e o - ph ] F e b ˜ k = k/ ( ωZ ) = 0 . .
01, 0 .
1, 1, 10, and 100, where ω = 2 πV /λ is the angular frequency and Z = ρV isthe seismic impedance. The normal and shear stiffnesscomponents are identical. We normalize the time t as˜ t = ( t − t ) / ( L/V − t ), where t = λ/ (2 V ) is the half-period of the wave emitted from the source. We computeup to ˜ t = 1 and 10 for the one-cycle and continuous ex-citation simulations, respectively. More details of ourmodel including material properties, numerical settings,a sensitivity analysis of spatial/temporal discretization,and a demonstration of the statistical convergence are inthe Supplemental Material [28]. Excitation of a one-cycle wave. —Fig. 1a shows thewavefield spatial distribution by plotting the normal-ized amplitude ˜ A (local amplitude A normalized by A ) at ˜ t = 1, for which no wave has escaped outof the domain. The wavefield evolution is quantifiedby the normalized gyradius [18, 31] (Fig. 1b), ˜ σ = L [ (cid:82)(cid:82) Ω ˜ Ir dxdy/ (cid:82)(cid:82) Ω ˜ Idxdy ] / , where ˜ I = ˜ A is the nor-malized wave intensity, and r is the distance between thewave position ( x , y ) and the source (0, 0). The temporalscaling of ˜ σ is expected to obey ˜ σ ∼ ˜ t ν with ν taking 1,0 .
5, and 0 for ballistic propagation, normal diffusion, andAnderson localization, respectively.For ˜ k = 100, fractures have no visible influence on thewavefield, whose circular wavefront is fully preserved ir-respective of p . For ˜ k = 10, the wavefront still keepsa circular shape, although some scattered waves are leftbehind it, which are more noticeable as p increases. Forboth ˜ k = 100 and 10, ˜ σ scales linearly with ˜ t ( ν = 1);some deviation from linearity exists in the early stage(˜ t < .
05) since the source excitation has not completedits full-period yet. For ˜ k = 1, the wavefront becomes con-siderably distorted with extensive scattered waves emerg-ing in its interior, such that the wavefield is character-ized by the fast propagation of frontal waves and slowdiffusion of coda waves. Increasing p tends to cause areduction of the wavefront speed and an enhancementof coda generation. Consequently, ˜ σ scales sublinearlywith ˜ t ( ν < k = 0 .
1, the wavefront is essen-tially destroyed and significant wave energy is backscat-tered/backreflected, phenomena which become more pro-nounced as p increases; the scaling of ˜ σ with ˜ t is com-patible with ν = 0 . k = 0 .
01, waves can still escape from the disconnectednetwork of p = 2 through the gaps between fractures,but are trapped near the source by the networks of p = 6and 10; wave diffusion becomes anomalous ( ν < . k = 0 . p = 10), qualified by a saturation of ˜ σ ( ν →
0) atthe late stage (˜ t > . p = 2 and 6) still exhibit abnormaldiffusion (0 < ν < . ν with the reduction of ˜ k (Fig. 2). Excitation of a continuous wave. —Fig. 3a shows thespatial distribution of ˜ A at ˜ t = 10, while Fig. 3b plotsthe temporal evolution of ˜ σ . For ˜ k = 100 and 10, thewavefield exhibits a ripple-like pattern consisting of a sus-tained wave train radially propagating outwards; ˜ σ scaleslinearly with ˜ t for ˜ t ≤ ν ≈ k = 1, the wavefield becomesscattered, more noticeable as p increases. However, thespatial footprints of waves are still relatively homoge-neous due to the dominant propagating part of the waveenergy (0 . < ν < k = 0 .
1, multiple scatteringprevails, resulting in highly-diffusive wavefields ( ν ≈ . p . For ˜ k = 0 .
01 and 0 . p = 2 is similar to that for ˜ k = 0 . ν ≈ . p = 6 becomes more subdiffusive orweakly localized ( ν < . p = 10 becomes more localized around the source whilealso globally more diffusive with ν ≈ .
5. Such a reen-trant phenomenon of delocalization is attributed to thedominance of evanescent waves [1], which leak wave en-ergy out of the closed “cavity” formed by interconnectedfractures. Since the properties of the matrix at the twosides of a fracture are the same, in order to generateevanescent waves, the only possible mechanism is that aS wave impinges on the fracture at an incidence angleexceeding the critical angle θ c = arcsin( V (cid:48) /V ) = 35 . ◦ ,associated with a decreased transmission of converted Pwaves [4]. This prerequisite can be easily found whenstrong multiple scattering occurs, promoting an equipar-tition of P and S waves [32] (although the source onlyexcites P waves) and distributing them into arbitrary di-rections [7]. The generated evanescent waves then propa-gate along fractures as interface waves [33, 34] at a speedasymptotically approaching that of Rayleigh waves forlow ˜ k .To further elucidate the emergence and consequence ofevanescent waves, we calculate the transmitted wave en-ergy at any given point within Ω as ˜ ψ = VL (cid:82) t ˜ Idt , whichis the normalized temporal integration of the wave in-tensity. Fig. 4 shows the spatial distribution of ˜ ψ atdifferent times ˜ t in a well-connected network of p = 10and ˜ k = 0 . t = 0). However, some waves (evanescent waves) man-age to escape the cavity and travel along fractures asinterface waves, marked by a fingering pattern in the ˜ ψ field (˜ t = 0 . p = 2 p = 10 p = 6(a) log à = 1 = 1 = 1 = 0.001= 0.01= 0.1 xy -2 -1 -2 -1 -2 -1 -2 -1 -2 -1 -2 -1 FIG. 1. (a) Wavefields illustrated by the normalized amplitude ˜ A at time ˜ t = 1 after exciting a circular, one-cycle P wavefrom the center of fractured media with different percolation parameters p and normalized fracture stiffnesses ˜ k . (b) Temporalvariation of the normalized gyradius ˜ σ ; shaded area indicates the standard deviation over 100 realizations; insets show ˜ σ versus˜ t in linear scales. Propagation-diffusion transitionDiffusion-localization transition Ballistic propagation ( ν = 1) Anderson localization ( ν = 0) Normal diffusion ( ν = 0.5) -3 -2 -1 FIG. 2. Dependence of the exponent ν as a function of thenormalized fracture stiffness ˜ k for fractured media with differ-ent percolation parameters p . Markers and error bars indicatethe mean and standard deviation over 100 realizations. (˜ t = 0 . k = 0 . (cid:104) ˜ I (cid:105) of the entire wavefield. The energy grows continu-ously in the connected networks of p = 6 and 10 dueto wave localization, while that of p = 2 tends to satu-rate because of outward diffusion. The probability den-sity function P ( ˜ I ) (Fig. 5b) shows that, as p increases,the intensity distribution changes from a lognormal-liketo a power law-like behavior. P ( ˜ I ) for p = 10 has amarginally heavier tail than that for p = 6 (bottom leftinset), suggesting a slightly stronger localization in theformer system. A higher proportion of intermediate in-tensity (0 . ≤ ˜ I <
1) is noticed for p = 10 compared to p = 6 (upper right inset) due to wave delocalization.We also perform a multifractal analysis to characterizethe spatial heterogeneity and hierarchy of the wavefield p = 2 p = 10 p = 6(a) log à = 1 = 1 = 1 = 0.001= 0.01= 0.1 xy -2 -1 -2 -1 -2 -1 -2 -1 -2 -1 -2 -1 FIG. 3. (a) Wavefields illustrated by the normalized amplitude ˜ A at time ˜ t = 10 after exciting a circular, continuous P wavefrom the center of fractured media with different percolation parameters p and normalized fracture stiffnesses ˜ k . (b) Temporalvariation of the normalized gyradius ˜ σ ; shaded area indicates the standard deviation over 100 realizations; insets show ˜ σ versus˜ t in linear scales. log xy = = = FIG. 4. Transmitted wave energy ˜ ψ at different times ˜ t inthe near-field ( − . ≤ x/L ≤ . − . ≤ y/L ≤ . p = 10 andnormalized fracture stiffness ˜ k = 0 . in the near-field of the source ( − . ≤ x/L ≤ . − . ≤ y/L ≤ . D q of order q (Fig. 5c) isdefined from the moment M q = (cid:80) ni =1 Θ qi ∼ ε ( q − D q [2],where the spatial domain is partitioned into n boxes ofsize ε , Θ is the sum of local intensities within the i thbox normalized by the total intensity. The multifractalspectrum [2] (Fig. 5d) is determined by computing theH¨older exponent α = (cid:80) ni =1 Θ qi ln Θ i / ( M q ln ε ), which is -10 -5 -9 -9 -7 -5 -1 -4 -2 -5 0 5135 FIG. 5. (a) Spatially-averaged wave intensity (cid:104) ˜ I (cid:105) versus time˜ t ; (b) probability density function of normalized intensity P ( ˜ I ); (c) generalized fractal dimension D q of order q and (d)singularity spectrum f α versus the H¨older exponent α for thenear-field intensity at ˜ t = 10. Insets give zoomed-in views. related to the singularity spectrum f α via the Legendretransform f α = αq − D q ( q − p increases from 2 to6 and then to 10, the high energy content of the wavesbecomes more and more localized, as revealed by the en-hanced nonlinearity of D q for q > f α in the small α region (Fig. 5d). Thelow-to-intermediate energy content also becomes more lo-calized when p increases from 2 to 6, associated with asteeper D q for q < f α for large α (Fig. 5d). However, as p further increases from 6 to10, D q for q < f α at large α narrows, im-plying that low-to-intermediate energy content is morediffuse and homogeneously partitioned over space in thewell-connected network, due to wave delocalization.To conclude, we have shown transitions of elastic wavetransport in fractured media from propagation to diffu- sion and to localization/delocalization, governed by frac-ture stiffness and network connectivity. Waves propagateballistically when the stiffness is large. As the stiffnessdecreases, multiple scattering occurs, leading to diffusivetransport in disconnected fracture networks and local-ization in connected ones. We have documented a reen-trant delocalization in well-connected networks due toevanescent waves and cascades of mode conversion. Wehave demonstrated that wave localization is qualified bya saturated growth of wavefield gyradius and a strongmultifractality of intensity field, whereas delocalizationdestroys these characteristics. This research opens thedoor to understanding complex wave phenomena in frac-tured media and has important implications for manygeophysical problems. [1] S. Stein and M. Wysession, An Introduction to Seismol-ogy, Earthquakes, and Earth Structure (Blackwell Pub-lishing Ltd, Oxford, 2003).[2] D. Sornette,
Critical Phenomena in Natural Sciences -Chaos, Fractals, Selforganization and Disorder: Con-cepts and Tools , 2nd ed. (Springer, Berlin, 2006).[3] Q. Lei and X. Wang, Geophys. Res. Lett. , 1551 (2016).[4] L. J. Pyrak-Nolte, L. R. Myer, and N. G. W. Cook, J.Geophys. Res. , 11345 (1990).[5] S. Vlastos, E. Liu, I. G. Main, and X.-Y. Li, Geophys. J.Int. , 649 (2003).[6] H. Hamzehpour, M. Asgari, and M. Sahimi, Phys. Rev.E , 063305 (2016).[7] A. Ishimaru, Wave Propagation and Scattering in Ran-dom Media (Academic Press, New York, 1978).[8] P. Sheng,
Introduction to Wave Scattering, Localizationand Mesoscopic Phenomena (Springer, Berlin, 2006).[9] D. Sornette, J. Stat. Phys. , 669 (1989).[10] D. Sornette, Acustica , 199 (1989).[11] A. F. Ioffe and A. R. Regel, Prog. Semicond. , 237(1960).[12] P. W. Anderson, Phys. Rev. , 1492 (1958).[13] D. Sornette, Acustica , 251 (1989).[14] M. C. W. van Rossum and T. M. Nieuwenhuizen, Rev.Mod. Phys. , 313 (1999).[15] P. W. Anderson, Philos. Mag. B , 505 (1985).[16] A. Lagendijk, B. V. Tiggelen, and D. S. Wiersma, Phys.Today , 24 (2009).[17] R. Dalichaouch, J. P. Armstrong, S. Schultz, P. M. Platz-man, and S. L. McCall, Nature , 53 (1991).[18] T. Sperling, W. B¨uhrer, C. M. Aegerter, and G. Maret,Nat. Photonics , 48 (2013).[19] T. Schwartz, G. Bartal, S. Fishman, and M. Segev, Na-ture , 52 (2007).[20] H. Hu, A. Strybulevych, J. H. Page, S. E. Skipetrov, andB. A. V. Tiggelen, Nat. Phys. , 945 (2008).[21] D. Sornette, Acustica , 15 (1989). [22] J. Billy, V. Josse, Z. Zuo, A. Bernard, B. Hambrecht,P. Lugan, D. Cl´ement, L. Sanchez-Palencia, P. Bouyer,and A. Aspect, Nature , 891 (2008).[23] E. Larose, L. Margerin, B. A. van Tiggelen, andM. Campillo, Phys. Rev. Lett. , 048501 (2004).[24] C. Friedrich and U. Wegler, Geophys. Res. Lett. , 1(2004).[25] K. Aki and V. Ferrazzini, J. Geophys. Res. , 16617(2000).[26] Q. Lei and D. Sornette (2020), J. Geophys. Res.Solid Earth (under review), e-print: 10.1002/es-soar.10504593.1.[27] M. Schoenberg, J. Acoust. Soc. Am. , 1516 (1980).[28] See Supplemental Material at [URL will be inserted bypublisher], for more details about the governing equa-tions, model setup, model validation, sensitivity analysis,and statistical convergence.[29] I. Balberg, C. H. Anderson, S. Alexander, and N. Wag-ner, Phys. Rev. B , 3933 (1984).[30] J. R. Pettit, A. Walker, P. Cawley, and M. J. S. Lowe,Ultrasonics , 1868 (2014).[31] P. Sebbah, D. Sornette, and C. Vanneste, Phys. Rev. B , 12506 (1993).[32] R. Hennino, N. Tregoures, N. M. Shapiro, L. Margerin,M. Campillo, B. A. van Tiggelen, and R. L. Weaver,Phys. Rev. Lett. , 3447 (2001).[33] L. J. Pyrak-Nolte and N. G. W. Cook, Geophys. Res.Lett. , 1107 (1987).[34] L. J. Pyrak-Nolte, J. Xu, and G. M. Haley, Phys. Rev.Lett. , 3650 (1992).[35] J. C. de Bremaecker, Geophysics , 253 (1958).[36] C. Castellani, C. D. Castro, and L. Peliti, J. Phys. A.Math. Gen. , L1099 (1986).[37] M. Schreiber and H. Grussbach, Phys. Rev. Lett. , 607(1991).[38] F. Evers and A. D. Mirlin, Rev. Mod. Phys.80