• ISSN 2096-8957
  • CN 10-1702/P

大数据时代的地球动力学数值计算方法回顾与展望

张贝 Gabriele Morra David A.Yuen HenryM. Tufo MatthewG. Knepley 陈石

引用本文:
Citation:

大数据时代的地球动力学数值计算方法回顾与展望

    作者简介: 张贝(1985-),男,助理研究员,主要从事计算地球动力学与数值模拟. E-mail:rular099@qq.com.
  • 中图分类号: P315

Numerical methods of geodynamics in big data era: Review and outlook

  • CLC number: P315

  • 摘要: 随着大数据时代的到来,计算地球动力学数值计算方法体系更加完善. 本文系统地回顾了传统数值模拟方法在计算地球动力学领域的应用进展,包括:有限差分法、有限单元法、谱方法和谱元法;并对近年来一些新发展的算法和应用前景进行了综述,如:不连续Galerkin法、小波方法和格子玻尔兹曼方法等. 本综述有助于读者以整体视角了解地球动力学数值计算方法的发展脉络,并对大数据时代下研究适应日益丰富的数据和新算法提供有益参考.
  • 图 1  有限元法模拟的Central Andes俯冲过程(修改自Salomon, 2018).图中有限元法可以使用非规则网格,灵活设置计算模型,其中展示模型计算50000年之后的结果,分别是(a)位移(km);(b)压力(MPa);(c)热流(mW/m2);(d)温度(K)

    Figure 1.  Central Andes subduction process simulated by finite element method (modified from Salomon, 2018). The finite element method can use irregular grid to flexibly set the calculation model, in which the results calculated by the model after 50000 years are shown, respectively: (a) displacement (km); (b) Pressure (MPa); (c) Heat flow (mW/m2); (d) Temperature (K)

    图 2  DGM模拟声波—弹性波传播(修改自Hecht-Nielsen, 1992

    Figure 2.  DGM simulates acoustic wave-elastic wave propagation (modified from Hecht-Nielsen, 1992)

    图 3  小波方法表示的二维图像(Mallat, 1999). (a)原始图像f(N),分辨率N=256×256;(b)原始图像用正交小波分解的系数矩阵<f,Ψk.*?>=>j,n>,k=1,2,3,4,小波尺度因子为2j,大于阈值T的系数用黑点表示;(c)使用尺度因子最大的3组小波系数(N/16个)重建图像;(d)使用最大的N/16个系数重建图像

    Figure 3.  The wavelet method represents a two-dimensional image (Mallat, 1999). (a) Original image f(N), resolution N=256×256; (b) the original image using orthogonal wavelet decomposition coefficient matrix <f,Ψk.*?>=>j,n>, k = 1, 2, 3, 4, wavelet scale factor is 2j, is greater than the threshold value of T coefficient with black spots; (c) Three sets of wavelet coefficients (N/16) with the largest scale factor were used to reconstruct the image; (d) Reconstruct the image using the maximum N/16 coefficients

    图 4  火星表面重力和地形的小波展示(修改自Kido and Yuen, 2003).(a)原始重力;(b)原始地形;(c)不同尺度小波重建的重力;(d)不同尺度小波重建的地形. l.*?>=>w=8、16、32、64对应的尺度分别是424 km、212 km、106 km、53 km

    Figure 4.  A wavelet display for the gravity and topography of Mars (modified from Kido and Yuen, 2003). (a) Original gravity; (b) Original topography; (c) The gravity reconstruction by wavelets of different scales; (d) The topography reconstruction by wavelets of different scales. The corresponding scales of l.*?>=>w=8, 16, 32 and 64 are 424 km, 212 km, 106 km and 53 km respectively.

  • [1]

    Abboud T, Barbier D. Hi-BoX: A generic library of fast solvers for boundary element methods[C]// 2016 IEEE Conference on Antenna Measurements Applications (CAMA). Syracuse NY USA, 2016: 1-4.
    [2]

    Asgari A, Moresi L. Multiscale particle-in-cell method: from fluid to solid mechanics. In: Advanced Methods for Practical Applications in Fluid Mechanics[M]. InTech, Croatia, 2012: 185-208.
    [3]

    Audet P, Mareschal J -C. Wavelet analysis of the coherence between Bouguer gravity and topography: application to the elastic thickness anisotropy in the Canadian Shield[J]. Geophysical Journal International, 2007, 168(1): 287-298. doi: 10.1111/j.1365-246X.2006.03231.x
    [4]

    Balachandar S, Yuen D A. Three-dimensional fully spectral numerical method for mantle convection with depth-dependent properties[J]. Journal of computational physics, 1994, 113(1): 62-74. doi: 10.1006/jcph.1994.1118
    [5]

    Balay S, Abhyankar S, Adams M, et al. PETSc Users Manual Revision 3.7[M]. Argonne National Laboratory ANL, Lemont, IL, 2016.
    [6]

    Barnes J, Hut P. A hierarchical O(NlogN) force-calculation algorithm[J]. Nature, 1986, 324(6096): 446-449. doi: 10.1038/324446a0
    [7]

    Bauer S, Huber M, Ghelichkhan S, et al. Large-scale simulation of mantle convection based on a new matrix-free approach[J]. Journal of computational science, 2019, 31: 60-76. doi: 10.1016/j.jocs.2018.12.006
    [8]

    Benedetti I, Aliabadi M H, Davì G. 45. A fast 3D dual boundary element method based on hierarchical matrices[J]. International Journal of Solids and Structures, 2008(7-8): 2355-2376.
    [9]

    Bhatnagar P L, Gross E P, Krook M. A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems[J]. Physics Review, American Physical Society, 1954, 94(3): 511-525.
    [10]

    Bozdağ E, Peter D, Lefebvre M, et al. Global adjoint tomography: first-generation model[J]. Geophysical Journal International, 2016, 207(3): 1739-1766. doi: 10.1093/gji/ggw356
    [11]

    Braeck S, Podladchikov Y Y, Medvedev S. Spontaneous dissipation of elastic energy by self-localizing thermal runaway[J]. Physical review. E, Statistical, nonlinear, and soft matter physics, 2009, 80(4 Pt 2): 046105.
    [12]

    Brogi F, Ripepe M, Bonadonna C. Lattice Boltzmann modeling to explain volcano acoustic source[J]. Scientific reports, 2018, 8(1): 9537. doi: 10.1038/s41598-018-27387-0
    [13]

    Brune P R, Knepley M G, Smith B F, et al. Composing Scalable Nonlinear Algebraic Solvers[J]. SIAM Review, Society for Industrial and Applied Mathematics, 2015, 57(4): 535-565.
    [14]

    Buffett B A, Becker T W. Bending stress and dissipation in subducted lithosphere[J]. Journal of Geophysical Research: Solid Earth, 2012, 117: B05413.
    [15]

    Buffett B, Rowley D. Plate bending at subduction zones: Consequences for the direction of plate motions[J]. Earth and Planetary Science Letters, 2006, 245(1-2): 359-364. doi: 10.1016/j.jpgl.2006.03.011
    [16]

    Burstedde C, Ghattas O, Gurnis M, et al. Scalable adaptive mantle convection simulation on petascale supercomputers[C]// SC'08: Proceedings of the 2008 ACM/IEEE Conference on Supercomputing, Austin Texas, 2008: 1-15.
    [17]

    Burstedde C, Stadler G, Alisic L, et al. Large-scale adaptive mantle convection simulation[J]. Geophysical Journal International, 2013, 192(3): 889-906. doi: 10.1093/gji/ggs070
    [18]

    Cai X-C, Keyes D E. Nonlinearly preconditioned inexact Newton algorithms[J]. SIAM Journal of Scientific Computing, Society for Industrial and Applied Mathematics, 2002, 24(1): 183-200. doi: 10.1137/S106482750037620X
    [19]

    Capitanio F A, Morra G, Goes S. Dynamics of plate bending at the trench and slab-plate coupling[J]. Geochemistry, Geophysics, Geosystems, 2009, 10(4): 1-15.
    [20]

    Cazeaux P, Zahm O. A fast boundary element method for the solution of periodic many-inclusion problems via hierarchical matrix techniques[J]. ESAIM: Proceedings and Surveys, 2015, 48: 156-168. doi: 10.1051/proc/201448006
    [21]

    Chaillat S, Bonnet M. Recent advances on the fast multipole accelerated boundary element method for 3D time-harmonic elastodynamics[J]. Wave Motion, 2013, 50(7): 1090-1104. doi: 10.1016/j.wavemoti.2013.03.008
    [22]

    Chaillat S, Semblat J-F, Bonnet M. A preconditioned 3-D multi-region fast multipole solver for seismic wave propagation in complex geometries[J]. Communications in computational physics, 2012, 11(2): 594-609. doi: 10.4208/cicp.231209.030111s
    [23]

    Chen Y, Li Y, Valocchi A J, et al. Lattice Boltzmann simulations of liquid CO.*?>=>2 displacing water in a 2D heterogeneous micromodel at reservoir pressure conditions[J]. Journal of Contaminant Hydrology, 2018, 212: 14-27. doi: 10.1016/j.jconhyd.2017.09.005
    [24]

    Choblet G, Čadek O, Couturier F, et al. OEDIPUS: a new tool to study the dynamics of planetary interiors[J]. Geophysical Journal International, 2007, 170(1): 9-30. doi: 10.1111/j.1365-246X.2007.03419.x
    [25]

    Chorin A J. Numerical study of slightly viscous flow[J]. Journal of fluid mechanics, 1973, 57(4): 785-796. doi: 10.1017/S0022112073002016
    [26]

    Chorin A J. A numerical method for solving incompressible viscous flow problems[J]. Journal of computational physics, 1997, 135(2): 118-125. doi: 10.1006/jcph.1997.5716
    [27]

    Christensen U. Convection with pressure- and temperature-dependent non-Newtonian rheology[J]. Geophysical Journal International, 1984, 77(2): 343-384. doi: 10.1111/j.1365-246X.1984.tb01939.x
    [28]

    Cockburn B, Karniadakis G E, Shu C-W. Discontinuous Galerkin Methods: Theory, Computation and Applications[M]. Berlin: Springer Science & Business Media, 2012.
    [29]

    Colella P, Woodward P R. The Piecewise Parabolic Method (PPM) for gas-dynamical simulations[J]. Journal of computational physics, 1984, 54(1): 174-201. doi: 10.1016/0021-9991(84)90143-8
    [30]

    Conrad C P, Hager B H. Effects of plate bending and fault strength at subduction zones on plate dynamics[J]. Journal of Geophysical Research: Solid Earth, 1999, 104(B8): 17551-17571. doi: 10.1029/1999JB900149
    [31]

    Dabrowski M, Krotkiewski M, Schmid D W. MILAMIN: MATLAB-based finite element method solver for large problems[J]. Geochemistry, Geophysics, Geosystems, 2008, 9(4): 1-24.
    [32]

    Dangla P, Semblat J-F, Xiao H, et al. A simple and efficient regularization method for 3D BEM: Application to frequency-domain elastodynamics[J]. Bulletin of the Seismological Society of America, 2005, 95(5): 1916-1927. doi: 10.1785/0120050012
    [33]

    Davies D R, Wilson C R, Kramer S C. Fluidity: A fully unstructured anisotropic adaptive mesh computational modeling framework for geodynamics[J]. Geochemistry, Geophysics, Geosystems, 2011, 12(6): Q06001.
    [34]

    Di Ilio G, Chiappini D, Ubertini S, et al. Hybrid lattice Boltzmann method on overlapping grids[J]. Physical Review E, 2017, 95(1-1): 013309.
    [35]

    Dupont T, Hoffman J, Johnson C, et al. The Fenics Project[M]. Chalmers Finite Element Centre, Chalmers University of Technology, 2003.
    [36]

    Duretz T, Räss L, Podladchikov Y Y, et al. Resolving thermomechanical coupling in two and three dimensions: spontaneous strain localization owing to shear heating[J]. Geophysical Journal International, 2019, 216(1): 365-379. doi: 10.1093/gji/ggy434
    [37]

    Duretz T, Souche A, Borst R, et al. The benefits of using a consistent tangent operator for viscoelastoplastic computations in geodynamics[J]. Geochemistry, Geophysics, Geosystems, 2018, 171: 374.
    [38]

    Erlebacher G, Yuen D A, Dubuffet F. Case study: Visualization and analysis of high Rayleigh number-3D convection in the Earth’s mantle[C]// IEEE Visualization, 2002: 493-496.
    [39]

    Fang H-R, Saad Y. Two classes of multisecant methods for nonlinear acceleration[J]. Numerical Linear Algebra with Applications, 2009, 16(3): 197-221. doi: 10.1002/nla.617
    [40]

    Feichtinger C, Donath S, Köstler H, et al. WaLBerla: HPC software design for computational engineering simulations[J]. Journal of Computational Science, 2011, 2(2): 105-112. doi: 10.1016/j.jocs.2011.01.004
    [41]

    Fichtner A. Full Seismic Waveform Modelling and Inversion[M]. Berlin: Springer Science & Business Media, 2010.
    [42]

    Fischer P F, Kruse G W, Loth F. Spectral elementmethods for transitional flows in complex geometries[J]. Journal of scientific computing, 2002, 17(1): 81-98.
    [43]

    Fornberg B. A Practical Guide to Pseudospectral Methods[M]. Cambridge: Cambridge University Press, 1998.
    [44]

    Fowler A C. A mathematical model of magma transport in the asthenosphere[J]. Geophysical and Astrophysical Fluid Dynamics, 1985, 33(1-4): 63-96. doi: 10.1080/03091928508245423
    [45]

    Frankel S P. Convergence rates of iterative treatments of partial differential equations[J]. Mathematical Tables and Other Aids to Computation, 1950, 4(30): 65-75. doi: 10.2307/2002770
    [46]

    Fraters M, Thieulot C, van den Berg A, et al. The Geodynamic World Builder: a solution for complex initial conditions in numerical modeling[J]. Solid Earth Discussions, 2019, 10(5): 1785-1807. doi: 10.5194/se-10-1785-2019
    [47]

    Frisch U, Hasslacher B, Pomeau Y. Lattice-gas automata for the Navier-Stokes equation[J]. Physical review letters, 1986, 56(14): 1505-1508. doi: 10.1103/PhysRevLett.56.1505
    [48]

    Gerardi G, Ribe N M. Boundary element modeling of two-plate interaction at subduction zones: Scaling laws and application to the aleutian subduction zone[J]. Journal of Geophysical Research, 2018, 123(6): 5227-5248.
    [49]

    Gerardi G, Ribe N M, Tackley P J. Plate bending, energetics of subduction and modeling of mantle convection: A boundary element approach[J]. Earth and Planetary Science Letters, 2019, 515: 47-57. doi: 10.1016/j.jpgl.2019.03.010
    [50]

    Gerya T. Introduction to Numerical Geodynamic Modelling[M]. Cambridge: Cambridge University Press, 2019.
    [51]

    Gharti H N, Langer L, Tromp J. Spectral-infinite-element simulations of earthquake-induced gravity perturbations[J]. Geophysical Journal Internationala, 2019a, 217(1): 451-468. doi: 10.1093/gji/ggz028
    [52]

    Gharti H N, Langer L, Tromp J. Spectral-infinite-element simulations of coseismic and post-earthquake deformation[J]. Geophysical Journal Internationalb, 2019b, 216(2): 1364-1393. doi: 10.1093/gji/ggy495
    [53]

    Gharti H N, Tromp J. Spectral-infinite-element simulations of magnetic anomalies[J]. Geophysical Journal International, 2019, 217(3): 1656-1667. doi: 10.1093/gji/ggz107
    [54]

    Glatzmaier G A. Numerical simulations of mantle convection: Time-dependent, three-dimensional, compressible, spherical shell[J]. Geophysical and Astrophysical Fluid Dynamics, 1988, 43(2): 223-264. doi: 10.1080/03091928808213626
    [55]

    Glatzmaier G A. Introduction to Modeling Convection in Planets and Stars: Magnetic Field, Density Stratification, Rotation[M]. New Jersey: Princeton University Press, 2013.
    [56]

    Gmeiner B, Rüde U, Stengel H, et al. Performance and scalability of hierarchical hybrid multigrid solvers for Stokes systems[J]. SIAM Journal of Scientific Computing, 2015, 37(2): C143-C168. doi: 10.1137/130941353
    [57]

    Griffiths D V, Smith I M. Numerical Methods for Engineers, Second Edition[M]. Boca Raton: CRC Press, 2006.
    [58]

    Guo Z, Shi B, Zheng C. A coupled lattice BGK model for the Boussinesq equations[J]. International Journal for Numerical, 2002, 39(4): 325-342.
    [59]

    Guzina B B, Pak R Y S. On the analysis of wave motions in a multi-layered solid[J]. The Quarterly Journal of Mechanics and Applied Mathematics, 2001, 54(1): 13-37. doi: 10.1093/qjmam/54.1.13
    [60]

    Harder H, Hansen U. A finite-volume solution method for thermal convection and dynamo problems in spherical shells[J]. Geophysical Journal International, 2005, 161(2): 522-532. doi: 10.1111/j.1365-246X.2005.02560.x
    [61]

    He X, Chen S, Doolen G D. A novel thermal model for the lattice boltzmann method in incompressible limit[J]. Journal of computational physics, 1998, 146(1): 282-300. doi: 10.1006/jcph.1998.6057
    [62]

    Hecht-Nielsen R. Theory of the backpropagation neural network[J]. Neural Networks, 1992: 65-93.
    [63]

    Hennenfent G, Herrmann F J. Seismic Denoising with nonuniformly sampled curvelets[J]. Computing in Science & Engineering, 2006, 8(3): 16-25.
    [64]

    Hennenfent G, Herrmann F J. Simply denoise: Wavefield reconstruction via jittered undersampling[J]. Geophysics, 2008, 73(3): V19-V28. doi: 10.1190/1.2841038
    [65]

    Herman F, Seward D, Valla P G, et al. Worldwide acceleration of mountain erosion under a cooling climate[J]. Nature, 2013, 504(7480): 423-426. doi: 10.1038/nature12877
    [66]

    Hesthaven J S, Warburton T. Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications[M]. Berlin: Springer Science & Business Media, 2007.
    [67]

    Heuveline V, Latt J. THE OPENLB PROJECT: AN OPEN SOURCE AND OBJECT ORIENTED IMPLEMENTATION OF LATTICE BOLTZMANN METHODS[J]. International Journal of Modern Physics C, 2007, 18(04): 627-634. doi: 10.1142/S0129183107010875
    [68]

    Higuera F J, Jiménez J. Boltzmann Approach to Lattice gas simulations[J]. Europhysics Letters, 1989, 9(7): 663-668. doi: 10.1209/0295-5075/9/7/009
    [69]

    Huang H, Sukop M C, Lu X-Y. Multiphase lattice Boltzmann methods: Theory and application[M]. Wiley-Blackwell, 2015.
    [70]

    Hwang F-N, Cai X-C. Improving robustness and parallel scalability of Newton method through nonlinear preconditioning[C]// Domain Decomposition Methods in Science and Engineering. Springer Berlin Heidelberg, 2005: 201-208.
    [71]

    Hwang F-N, Su Y-C, Cai X-C. A parallel adaptive nonlinear elimination preconditioned inexact Newton method for transonic full potential equation[J]. Computers & fluids, 2015, 110: 96-107.
    [72]

    Igel H. Computational Seismology: A Practical Introduction[M]. Oxford: Oxford University Press, 2017.
    [73]

    Ismail-Zadeh A, Tackley P. Computational Methods for Geodynamics[M]. Cambridge: Cambridge University Press, 2010.
    [74]

    Jadamec M A, Billen M I. Reconciling surface plate motions with rapid three-dimensional mantle flow around a slab edge[J]. Nature, 2010, 465(7296): 338-341. doi: 10.1038/nature09053
    [75]

    Jain C, Rozel A B, Tackley P J, et al. Growing primordial continental crust self-consistently in global mantle convection models[J]. Gondwana Research, 2019, 73: 96-122. doi: 10.1016/j.gr.2019.03.015
    [76]

    Jarvis G T, Mckenzie D P. Convection in a compressible fluid with infinite Prandtl number[J]. Journal of fluid mechanics, 1980, 96(3): 515-583. doi: 10.1017/S002211208000225X
    [77]

    Jing G, Zheng Y, Wang L. Consensus of multiagent systems with distance-dependent communication networks[J]. IEEE Transactions on Neural Networks and Learning Systems, 2017, 28(11): 2712-2726. doi: 10.1109/TNNLS.2016.2598355
    [78]

    John V, Knobloch P, Novo J. Finite elements for scalar convection-dominated equations and incompressible flow problems: a never ending story?[J]. Computing and Visualization in Science, 2018, 19(5): 47-63.
    [79]

    Jolliffe I. Principal component analysis[M]. Berlin Heidelberg: Springer, 2011.
    [80]

    Kageyama A, Sato T. “Yin-Yang grid”: An overset grid in spherical geometry[J]. Geochemistry, Geophysics, Geosystems, 2004, 5(9): Q09005.
    [81]

    Kameyama M, Kageyama A, Sato T. Multigrid iterative algorithm using pseudo-compressibility for three-dimensional mantle convection with strongly variable viscosity[J]. Journal of computational physics, 2005, 206(1): 162-181. doi: 10.1016/j.jcp.2004.11.030
    [82]

    Kameyama M, Kageyama A, Sato T. Multigrid-based simulation code for mantle convection in spherical shell using Yin–Yang grid[J]. Physics of the Earth and Planetary Interiors, 2008, 171(1): 19-32.
    [83]

    Kang Q, Lichtner P C, Zhang D. Lattice Boltzmann pore-scale model for multicomponent reactive transport in porous media[J]. Journal of Geophysical Research: Solid Earth, 2006, 111: B05203.
    [84]

    Kido M, Yuen D A. Mapping the admittance and correlation functions on the Martian surface using a spherical wavelet fiido MJ][J]. Physics of the Earth and Planetary Interiors, 2003, 135(1): 1-16. doi: 10.1016/S0031-9201(02)00176-0
    [85]

    Kiss D, Podladchikov Y, Duretz T, et al. Spontaneous generation of ductile shear zones by thermal softening: Localization criterion, 1D to 3D modelling and application to the lithosphere[J]. Earth and planetary science letters, 2019, 519: 284-296. doi: 10.1016/j.jpgl.2019.05.026
    [86]

    Klawonn A, Lanser M, Rheinbach O, et al. Nonlinear FETI-DP and BDDC methods: Aunified framework and parallel results[J]. SIAM Journal of Scientific Computing, 2017, 39(6): C417-C451. doi: 10.1137/16M1102495
    [87]

    Kohl N, Thönnes D, Drzisga D, et al. The HyTeG finite-element software framework for scalable multigrid solvers[J]. International Journal of Parallel, 2019, 34(5): 477-496.
    [88]

    Komatitsch D, Vilotte J-P. The spectral element method: An efficient tool to simulate the seismic response of 2D and 3D geological structures[J]. Bulletin of the Seismological Society of America, 1998, 88(2): 368-392.
    [89]

    Kronbichler M, Heister T, Bangerth W. High accuracy mantle convection simulation through modern numerical methods[J]. Geophysics Journal International, 2012, 191: 12-29. doi: 10.1111/j.1365-246X.2012.05609.x
    [90]

    Krüger T, Kusumaatmaja H, Kuzmin A, et al. The lattice Boltzmann method[J]. Springer International Publishing, 2017, 10(978-3): 4-15.
    [91]

    Lagrava D, Malaspinas O, Latt J, et al. Advances in multi-domain lattice Boltzmann grid refinement[J]. Journal of Computational Physics, 2012, 231(14): 4808-4822. doi: 10.1016/j.jcp.2012.03.015
    [92]

    Lanzkron P J, Rose D J, Wilkes J T. An analysis of approximate nonlinear elimination[J]. SIAM Journal of Scientific Computing, 1996, 17(2): 538-559. doi: 10.1137/S106482759325154X
    [93]

    Larsen T B, Yuen D A, Moser J, et al. A high-order finite-difference method applied to large Rayleigh number mantle convection[J]. Geophysical and Astrophysical Fluid Dynamics, 1997, 84(1-2): 53-83. doi: 10.1080/03091929708208973
    [94]

    Leclaire S, Parmigiani A, Malaspinas O, et al. Generalized three-dimensional lattice Boltzmann color-gradient method for immiscible two-phase pore-scale imbibition and drainage in porous media[J]. Physical Review E, 2017, 95(3-1): 033306.
    [95]

    Leng W, Zhong S. Implementation and application of adaptive mesh refinement for thermochemical mantle convection studies[J]. Geochemistry, Geophysics, Geosystems, 2011, 12(4): Q04006.
    [96]

    Li Z H, Di Leo J F, Ribe N M. Subduction-induced mantle flow, finite strain, and seismic anisotropy: Numerical modeling[J]. Journal of geophysical research, 2014, 119(6): 5052-5076.
    [97]

    Li Z-H, Ribe N M. Dynamics of free subduction from 3-D boundary element modeling[J]. Journal of Geophysical Research, 2012, 117: B06408.
    [98]

    Liu S, Zhong Z, Takbiri-Borujeni A, et al. A case study on homogeneous and heterogeneous reservoir porous media reconstruction by using generative adversarial networks[J]. Energy Procedia, 2019, 158: 6164-6169. doi: 10.1016/j.egypro.2019.01.493
    [99]

    Luo H, Luo L Q, Nourgaliev R, et al. A reconstructed discontinuous Galerkin method for the compressible Navier–Stokes equations on arbitrary grids[J]. Journal of Computational Physics, 2010, 229(19): 6961-6978.
    [100]

    Luporini F, Louboutin M, Lange M, et al. Architecture and performance of Devito, a system for automated stencil computation[J]. ACM Transactions on Mathematical Software, 2020, 46(1): 1-28.
    [101]

    Lynch D R. Numerical Partial Differential Equations for Environmental Scientists and Engineers: A First Practical Course[M]. Berlin: Springer Science & Business Media, 2004.
    [102]

    Machetel P, Yuen D A. The onset of time-dependent convection in spherical shells as a clue to chaotic convection in the Earth's mantle[J]. Geophysical Research Letters, 1986, 13(13): 1470-1473. doi: 10.1029/GL013i013p01470
    [103]

    MacNeice P, Olson K M, Mobarry C, et al. PARAMESH: A parallel adaptive mesh refinement community toolkit[J]. Computer Physics Communications, 2000, 126(3): 330-354. doi: 10.1016/S0010-4655(99)00501-9
    [104]

    Mai P M, Schorlemmer D, Page M, et al. The earthquake-source inversion validation (SIV) project[J]. Seismological Research Letters, 2016, 87(3): 690-708. doi: 10.1785/0220150231
    [105]

    Malevsky A V. Spline-characteristic method for simulation of convective turbulence[J]. Journal of Computational Physics, 1996, 123(2): 466-475. doi: 10.1006/jcph.1996.0037
    [106]

    Malevsky A V, Yuen D A. Plume structures in the hard-turbulent regime of three-dimensional infinite Prandtl number convection[J]. Geophysical Research Letters, 1993, 20(5): 383-386. doi: 10.1029/93GL00293
    [107]

    Mallat S. A Wavelet Tour of Signal Processing[M]. Elsevier, 1999.
    [108]

    Manga M, Stone H A. Low Reynolds number motion of bubbles, drops and rigid spheres through fluid-fluid interfaces[J]. Journal of Fluid Mechanics, 1995, 287: 279-298. doi: 10.1017/S0022112095000954
    [109]

    Margrave G F, Lamoureux M P. Numerical methods of exploration seismology: with algorithms in MATLAB[M]. Cambridge: Cambridge University Press, 2019.
    [110]

    May D A, Moresi L. Preconditioned iterative methods for Stokes flow problems arising in computational geodynamics[J]. Physics of the Earth and Planetary Interiors, 2008, 171(1): 33-47.
    [111]

    May D A, Schellart W P, Moresi L. Overview of adaptive finite element analysis in computational geodynamics[J]. Journal of Geodynamics, Elsevier, 2013, 70: 1-20. doi: 10.1016/j.jog.2013.04.002
    [112]

    McKENZIE D A N. The Generation and Compaction of Partially Molten Rock[J]. Journal of Petrology, 1984, 25(3): 713-765. doi: 10.1093/petrology/25.3.713
    [113]

    McKenzie D P, Roberts J M, Weiss N O. Convection in the Earth's mantle: towards a numerical simulation[J]. Journal of fluid mechanics, 1974, 62(3): 465-538. doi: 10.1017/S0022112074000784
    [114]

    Middleton R S, William Carey J, Currier R P, et al. Shale gas and non-aqueous fracturing fluids: Opportunities and challenges for supercritical CO.*?>=>2[J]. Applied Energy, 2015, 147: 500-509. doi: 10.1016/j.apenergy.2015.03.023
    [115]

    Mora P, Morra G, Yuen D A. A concise python implementation of the lattice Boltzmann method on HPC for geo-fluid flow[J]. Geophysical Journal International, 2020, 220(1): 682-702. doi: 10.1093/gji/ggz423
    [116]

    Mora P, Yuen D A. Simulation of plume dynamics by the Lattice Boltzmann Method[J]. Geophysical Journal International, 2017, 210(3): 1932-1937. doi: 10.1093/gji/ggx279
    [117]

    Mora P, Yuen D A. Comparison of convection for Reynolds and Arrhenius temperature dependent viscosities[J]. Fluid Mechanics Research International Journal, 2018, 2(3): 99-107.
    [118]

    Moresi L, Dufour F, Mühlhaus H-B. A Lagrangian integration point finite element method for large deformation modeling of viscoelastic geomaterials[J]. Journal of computational physics, 2003, 184(2): 476-497. doi: 10.1016/S0021-9991(02)00031-1
    [119]

    Moresi L-N, Solomatov V S. Numerical investigation of 2D convection with extremely large viscosity variations[J]. Physics of Fluids, 1995, 7(9): 2154-2162. doi: 10.1063/1.868465
    [120]

    Morra G. Pythonic Geodynamics[M].Springer, 2018.
    [121]

    Morra G, Chatelain P, Tackley P, et al. Large scale three-dimensional boundary element simulation of subduction[J]. Computational Science, 2007: 1122-1129.
    [122]

    Morra G, Chatelain P, Tackley P, et al. Earth curvature effects on subduction morphology: Modeling subduction in a spherical setting[J]. Acta Geotechnica, 2009, 4(2): 95-105. doi: 10.1007/s11440-008-0060-5
    [123]

    Morra G, Regenauer-Lieb K. A coupled solid–fluid method for modelling subduction[J]. Philosophical Magazine, 2006, 86(21-22): 3307-3323. doi: 10.1080/14786430500256359
    [124]

    Morra G, Yuen D A, Lee S-M, et al. Source of the cenozoic volcanism in central Asia. In: AGU Geophysical Monograph(ed) Subduction Dynamics: From Mantle Flow to Mega Disasters[M]. John Wiley & Sons, Inc., 2015, 211: 97-113.
    [125]

    Ogawa M, Schubert G, Zebib A. Numerical simulations of three-dimensional thermal convection in a fluid with strongly temperature-dependent viscosity[J]. Journal of Fluid Mechanics, 1991, 233: 299-328. doi: 10.1017/S0022112091000496
    [126]

    Omlin S, Räss L, Podladchikov Y Y. Simulation of three-dimensional viscoelastic deformation coupled to porous fluid flow[J]. Tectonophysics, 2018, 746: 695-701. doi: 10.1016/j.tecto.2017.08.012
    [127]

    Oosterlee C W, Lorenz F J G. Multigrid methods for the Stokes system[J]. Computing in Science Engineering, 2006, 8(6): 34-43. doi: 10.1109/MCSE.2006.115
    [128]

    Oosterlee C W, Washio T. Krylov Subspace Acceleration of Nonlinear Multigrid with Application to Recirculating Flows[J]. SIAM Journal of Scientific Computing, 2000, 21(5): 1670-1690. doi: 10.1137/S1064827598338093
    [129]

    Parmentier E M, Turcotte D L, Torrance K E. Numerical experiments on the structure of mantle plumes[J]. Journal of Geophysical Research, 1975, 80(32): 4417-4424. doi: 10.1029/JB080i032p04417
    [130]

    Parmigiani A, Degruyter W, Leclaire S, et al. The mechanics of shallow magma reservoir outgassing[J]. Geochemistry, Geophysics, Geosystems, 2017, 18(8): 2887-2905. doi: 10.1002/2017GC006912
    [131]

    Parmigiani A, Faroughi S, Huber C, et al. Bubble accumulation and its role in the evolution of magma reservoirs in the upper crust[J]. Nature, 2016, 532(7600): 492-495. doi: 10.1038/nature17401
    [132]

    Patankar S. Numerical Heat Transfer and Fluid Flow[M]. Taylor & Francis, 2018.
    [133]

    Phillips P J. Finite Element Methods in Linear Poroelasticity: Theoretical and Computational Results[D]. Austin, TX, USA: University of Texas at Austin, 2005.
    [134]

    Phillips T N, Davies J H, Oldham D. Towards global SEM mantle convection simulations on polyhedral-based grids[J]. Journal of Computational and Applied Mathematics, 2019, 348: 48-57. doi: 10.1016/j.cam.2018.08.042
    [135]

    Popov A, Kaus B. LaMEM (Lithosphere and Mantle Evolution Model): advancing a staggered-grid finite difference version of the code[A]. ui.adsabs.harvard.edu, 2013: EGU2013-7761.
    [136]

    Pozrikidis C. A Practical Guide to Boundary Element Methods with the Software Library BEMLIB[M]. CRC Press, 2002.
    [137]

    Pradana B, Fauzi U. Three-dimensional modeling of multiphase fluid flow with linear temperature gradients in porous media using lattice Boltzmann method Rothman-Keller[J]. Journal of Physics, 2019, 1245: 012026.
    [138]

    Prieto G A, Florez M, Barrett S A, et al. Seismic evidence for thermal runaway during intermediate-depth earthquake rupture[J]. Geophysical Research Letters, 2013, 40(23): 6064-6068. doi: 10.1002/2013GL058109
    [139]

    Qian Y H, D'Humières D, Lallemand P. Lattice BGK models for Navier-Stokes equation[J]. EPL, IOP Publishing, 1992, 17(6): 479-484.
    [140]

    Räss L, Duretz T, Podladchikov Y Y, et al. M2Di: Concise and efficient MATLAB 2-DS tokes solvers using the Finite Difference Method[J]. Geochemistry, Geophysics, Geosystems, 2017, 18(2): 755-768. doi: 10.1002/2016GC006727
    [141]

    Räss L, Simon N S C, Podladchikov Y Y. Spontaneous formation of fluid escape pipes from subsurface reservoirs[J]. Scientific Reports, 2018, 8(1): 11116. doi: 10.1038/s41598-018-29485-5
    [142]

    Räss L, Duretz T, Podladchikov Y Y. Resolving hydromechanical coupling in two and three dimensions: spontaneous channelling of porous fluids owing to decompaction weakening[J]. Geophysical Journal International, 2019, 218(3): 1591-1616. doi: 10.1093/gji/ggz239
    [143]

    Ratcliff J T, Schubert G, Zebib A. Steady tetrahedral and cubic patterns of spherical shell convection with temperature-dependent viscosity[J]. Journal of Geophysical Research, 1996, 101(B11): 25473-25484. doi: 10.1029/96JB02097
    [144]

    Regenauer-Lieb K, Yuen D A. Rapid conversion of elastic energy into plastic shear heating during incipient necking of the lithosphere[J]. Geophysical Research Letters, 1998, 25(14): 2737-2740. doi: 10.1029/98GL02056
    [145]

    Regenauer-Lieb K, Yuen D A, Branlund J. The initiation of subduction: criticality by addition of water?[J]. Science, 2001, 294(5542): 578-580. doi: 10.1126/science.1063891
    [146]

    Reis T. The Lattice Boltzmann Method for Fluid Flows[D]. University of Wales Aberystwyth, 2003.
    [147]

    Richtmyer R D, Morton K W. Difference Methods for Initial-value Problems[J]. Interscience, New York, 1967.
    [148]

    Rozel A, Golabek G J, Näf R, et al. Formation of ridges in a stable lithosphere in mantle convection models with a viscoplastic rheology[J]. Geophysical research letters, 2015, 42(12): 4770-4777. doi: 10.1002/2015GL063483
    [149]

    Saad Y, Schultz M H. Conjugate gradient-like algorithms for solving nonsymmetric linear systems[J]. Mathematics of Computation, 1985, 44(170): 417-424. doi: 10.1090/S0025-5718-1985-0777273-9
    [150]

    Salomon C. Finite element modelling of the geodynamic processes of the Central Andes subduction zone: A reference model[J]. Geodesy and Geodynamics, 2018, 9(3): 246-251. doi: 10.1016/j.geog.2017.11.007
    [151]

    Salvadori A. Analytical integrations in 3D BEM for elliptic problems: Evaluation and implementation[J]. International Journal for Numerical Methods in Engineering, 2010, 84(5): 505-542. doi: 10.1002/nme.2906
    [152]

    Schaefer L, Ikeda M, Bao J. The Lattice Boltzmann Equation Method for Complex Flows[C]// ASME 2012 10th International Conference on Nanochannels, Microchannels, and Minichannels. American Society of Mechanical Engineers Digital Collection, 2013: 687-699.
    [153]

    Schmalholz S M, Duretz T. Impact of grain size evolution on necking in calcite layers deforming by combined diffusion and dislocation creep[J]. Journal of Structural Geology, 2017, 103: 37-56. doi: 10.1016/j.jsg.2017.08.007
    [154]

    Schmieschek S, Shamardin L, Frijters S, et al. LB3D: A parallel implementation of the Lattice-Boltzmann method for simulation of interacting amphiphilic fluids[J]. Computer Physics Communications, 2017, 217: 149-161. doi: 10.1016/j.cpc.2017.03.013
    [155]

    Schornbaum F, Rüde U. Massively Parallel Algorithms for the Lattice Boltzmann Method on NonUniform Grids[J]. SIAM Journal on Scientific Computing, 2016, 38(2): C96-C126. doi: 10.1137/15M1035240
    [156]

    Scott D R, Stevenson D J. Magma ascent by porous flow[J]. Journal of Geophysical Research, 1986, 91(B9): 9283. doi: 10.1029/JB091iB09p09283
    [157]

    Semblat J. F. Modeling seismic wave propagation and amplification in 1D/2D/3D linear and nonlinear unbounded media[J]. International Journal of Geomechanics, 2011, 11(6): 440-448. doi: 10.1061/(ASCE)GM.1943-5622.0000023
    [158]

    Shahnas M H, Yuen D A, Pysklywec R N. Inverse problems in geodynamics using machine learning llgorithms[J]. Journal of Geophysical Research, 2018, 123(1): 296-310.
    [159]

    Shan X, Chen H. Lattice Boltzmann model for simulating flows with multiple phases and components[J]. Physical Review E, 1993, 47(3): 1815-1819. doi: 10.1103/PhysRevE.47.1815
    [160]

    Shan X. Simulation of Rayleigh-Bénard convection using a lattice Boltzmann method[J]. Physical Review E, 1997, 55(3): 2780-2788. doi: 10.1103/PhysRevE.55.2780
    [161]

    Shearer P M. Introduction to Seismology[M]. Cambridge: Cambridge University Press, 2019.
    [162]

    Shukla K, Chan J, Maarten V, et al. A weight-adjusted discontinuous Galerkin method for the poroelastic wave equation: Penalty fluxes and micro-heterogeneities[J]. Journal of Computational Physics, 2020, 43: 109061.
    [163]

    Simo J C, Hughes T J R. Computational Inelasticity[M]. New York: Springer, 2000.
    [164]

    Sioshansi R, Conejo A J. Optimization in Engineering: Models and Algorithms[M]. Springer, 2019.
    [165]

    Sizova E, Gerya T, Stüwe K, et al. Generation of felsic crust in the Archean: a geodynamic modeling perspective[J]. Precambrian Research, 2015, 271: 198-224. doi: 10.1016/j.precamres.2015.10.005
    [166]

    Skarbek R M, Rempel A W. Dehydration-induced porosity waves and episodic tremor and slip[J]. Geochemistry, Geophysics, Geosystems, 2016, 17(2): 442-469. doi: 10.1002/2015GC006155
    [167]

    So B D, Yuen D A, Lee S M. An efficient implicit-explicit adaptive time stepping scheme for multiple-time scale problems in shear zone development[J]. Geochemistry, Geophysics, Geosystems, Wiley Online Library, 2013, 14(9): 3462-3478. doi: 10.1002/ggge.20216
    [168]

    Spiegelman M, Cox K G, McKenzie D P, et al. Physics of melt extraction: theory, implications and applications[J]. Philosophical Transactions of the Royal Society of London, 1993, 342(1663): 23-41.
    [169]

    Stadler G, Gurnis M, Burstedde C, et al. The dynamics of plate tectonics and mantle flow: from local to global scales[J]. Science, 2010, 329(5995): 1033-1038. doi: 10.1126/science.1191223
    [170]

    Stein S, Wysession M. An Introduction to Seismology, Earthquakes, and Earth Structure[M]. John Wiley & Sons, 2009.
    [171]

    Stemmer K, Harder H, Hansen U. A new method to simulate convection with strongly temperature- and pressure-dependent viscosity in a spherical shell: Applications to the Earth's mantle[J]. Physics of the Earth and Planetary Interiors, 2006, 157(3): 223-249.
    [172]

    Succi S, Succi S. The Lattice Boltzmann Equation: For Complex States of Flowing Matter[M]. Oxford: Oxford University Press, 2018.
    [173]

    Sulsky D, Chen Z, Schreyer H L. A particle method for history-dependent materials[J]. Computer Methods in Applied Mechanics and Engineering, 1994, 118(1): 179-196.
    [174]

    Tackley P J. Effects of strongly temperature-dependent viscosity on time-dependent, three-dimensional models of mantle convection[J]. Geophysical Research Letters, 1993, 20(20): 2187-2190. doi: 10.1029/93GL02317
    [175]

    Tackley P J. Self-consistent generation of tectonic plates in time-dependent, three-dimensional mantle convection simulations[J]. Geochemistry, Geophysics, Geosystems, 2000, 1(8): 2000GC000043.
    [176]

    Tackley P J. Modelling compressible mantle convection with large viscosity contrasts in a three-dimensional spherical shell using the yin-yang grid[J]. Physics of the Earth and Planetary Interiors, 2008, 171(1): 7-18.
    [177]

    Thielmann M, Rozel A, Kaus B J P, et al. Intermediate-depth earthquake generation and shear zone formation caused by grain size reduction and shear heating[J]. Geology, 2015, 43(9): 791-794. doi: 10.1130/G36864.1
    [178]

    Thompson T B, Meade B J. Boundary element methods for earthquake modeling with realistic 3D geometries[J]. EarthArXiv, 2019a.
    [179]

    Thompson T B, Meade B J. Earthquake cycle modeling of the Cascadia subduction zone[J]. EarthArXiv, 2019b.
    [180]

    Tornberg A-K, Greengard L. A fast multipole method for the three-dimensional Stokes equations[J]. Journal of Computational Physics, 2008, 227(3): 1613-1619. doi: 10.1016/j.jcp.2007.06.029
    [181]

    Torrance K E, Turcotte D L. Thermal convection with large viscosity variations[J]. Journal of Fluid Mechanics, 1971, 47(1): 113-125. doi: 10.1017/S002211207100096X
    [182]

    Tsuchiya J, Umemoto K. First-Principles Determination of the Dissociation Phase Boundary of Phase H MgSiO.*?>=>4H.*?>=>2[J]. Geophysical Research Letters, 2019, 46(13): 7333-7336. doi: 10.1029/2019GL083472
    [183]

    Tsuboi S, Ando K, Miyoshi T, et al. A 1.8 trillion degrees-of-freedom, 1.24 petaflops global seismic wave simulation on the K computer[J]. The International Journal of High Performance Computing Applications, 2016, 30(4): 411-422. doi: 10.1177/1094342016632596
    [184]

    Turcotte D L, Schubert G. Geodynamics[M]. Cambridge: Cambridge University Press, 2002.
    [185]

    Turcotte D L, Torrance K E, Hsui A T. Convection in the Earth's mantle[J]. Methods in Computational Physics, 1973, 13: 431-454.
    [186]

    Udias A, Buforn E. Principles of seismology[M]. Cambridge: Cambridge University Press, 2017.
    [187]

    Umemoto K, Wentzcovitch R M, Allen P B. Dissociation of MgSiO.*?>=>3 in the cores of gas giants and terrestrial exoplanets[J]. Science, 2006, 311(5763): 983-986. doi: 10.1126/science.1120865
    [188]

    van den Berg A, Segal G, Yuen D A. SEPRAN: A versatile finite-element package for a wide variety of problems in geosciences[J]. Journal of Earth Science, 2015, 26(1): 89-95. doi: 10.1007/s12583-015-0508-0
    [189]

    van den Berg A P, Yuen D A, Umemoto K, et al. Mass-dependent dynamics of terrestrial exoplanets using ab initio mineral properties[J]. Icarus, 2019, 317: 412-426. doi: 10.1016/j.icarus.2018.08.016
    [190]

    Van Kan J, Segal A, Vermolen F. Numerical Methods in Scientific Computing. Department of Applied Mathematics[M]. Netherlands: Delft University of Technology, 2008.
    [191]

    Wang J, Wang D, Lallemand P, et al. Lattice Boltzmann simulations of thermal convective flows in two dimensions[J]. Computers & Mathematics with Applications, 2013, 65(2): 262-286.
    [192]

    Xie C, Raeini A Q, Wang Y, et al. An improved pore-network model including viscous coupling effects using direct simulation by the lattice Boltzmann method[J]. Advances in Water Resources, 2017, 100: 26-34. doi: 10.1016/j.advwatres.2016.11.017
    [193]

    Yarushina V M, Podladchikov Y Y. (De) compaction of porous viscoelastoplastic media: Model formulation[J]. Journal of Geophysical Research, 2015, 120(6): 4146-4170.
    [194]

    Yarushina Viktoriya M., Podladchikov Yuri Y., Minakov Alexander, et al. On the mechanisms of stress-triggered seismic events during fluid injection[J]. Poromechanics VI, 2017: 795-800.
    [195]

    Yavneh I. Why multigrid methods are so efficient[J]. Computing in Science & Engineering, aip.scitation.org, 2006, 8(6): 12-22.
    [196]

    Ye R, de Hoop M V, Petrovitch C L, et al. A discontinuous Galerkin method with a modified penalty flux for the propagation and scattering of acousto-elastic waves[J]. Geophysical Journal International, 2016, 205(2): 1267-1289. doi: 10.1093/gji/ggw070
    [197]

    Yoshida M, Kageyama A. Application of the Yin-Yang grid to a thermal convection of a Boussinesq fluid with infinite Prandtl number in a three-dimensional spherical shell[J]. Geophysical Research Letters, 2004, 31(12): L12609.
    [198]

    Zhang S, Yuen D A. The influences of lower mantle viscosity stratification on 3D spherical-shell mantle convection[J]. Earth and Planetary Science Letters, 1995, 132(1): 157-166.
    [199]

    Zhang S, Yuen D A. Intense local toroidal motion generated by variable viscosity compressible convection in 3-D spherical-shell[J]. Geophysical Research Lettersa, 1996a, 23(22): 3135-3138. doi: 10.1029/96GL03029
    [200]

    Zhang S, Yuen D A. Various influences on plumes and dynamics in time-dependent, compressible mantle convection in 3-D spherical shell[J]. Physics of the Earth and Planetary Interiorsb, 1996b, 94(3): 241-267.
    [201]

    Zhao B, MacMinn C W, Primkulov B K, et al. Comprehensive comparison of pore-scale models for multiphase flow in porous media[C]// Proceedings of the National Academy of Sciences of the United States of America. National Acad Sciences, 2019, 116(28): 13799-13806.
    [202]

    Zhong S, Yuen D A, Moresi L N, et al. Numerical methods for mantle convection[J]. Treatise on Geophysics, 2007, 7: 227-252. doi: 10.1016/B978-044452748-6/00118-8
    [203]

    Zhong S, Zuber M T, Moresi L, et al. Role of temperature-dependent viscosity and surface plates in spherical shell models of mantle convection[J]. Journal of Geophysical Research, 2000, 105(B5): 11063-11082. doi: 10.1029/2000JB900003
    [204]

    Zhu B, Yan H, Zhong Y, et al. Relativistic HPIC-LBM and its application in large temporal-spatial turbulent magnetic reconnection. Part I. model development and validation[J]. Applied Mathematical Modelling, 2020a, 78: 932-967. doi: 10.1016/j.apm.2019.09.043
    [205]

    Zhu B, Yan H, Zhong Y, et al. Relativistic HPIC-LBM and its application in large temporal-spatial turbulent magnetic reconnection. Part II. Role of turbulence in the flux rope interaction[J]. Applied Mathematical Modelling, 2020b, 78: 968-988. doi: 10.1016/j.apm.2019.05.027
    [206]

    Zienkiewicz O C, Taylor R L, Sherwin S J, et al. On discontinuous Galerkin methods[J]. International Journal for Numerical Methods in Engineering, 2003, 58(8): 1119-1148. doi: 10.1002/nme.884
    [207]

    Zobin V M. Introduction to Volcanic Seismology[M]. Elsevier, 2012.
    [208]

    Zou Q, He X. On pressure and velocity boundary conditions for the lattice Boltzmann BGK model[J]. Physics of Fluids, 1997, 9(6): 1591-1598. doi: 10.1063/1.869307
  • 加载中
图(4)
计量
  • 文章访问数:  1529
  • HTML全文浏览量:  761
  • PDF下载量:  54
  • 被引次数: 0
出版历程
  • 收稿日期:  2020-06-25
  • 网络出版日期:  2020-09-15
  • 刊出日期:  2021-01-01

大数据时代的地球动力学数值计算方法回顾与展望

    作者简介: 张贝(1985-),男,助理研究员,主要从事计算地球动力学与数值模拟. E-mail:rular099@qq.com

摘要: 随着大数据时代的到来,计算地球动力学数值计算方法体系更加完善. 本文系统地回顾了传统数值模拟方法在计算地球动力学领域的应用进展,包括:有限差分法、有限单元法、谱方法和谱元法;并对近年来一些新发展的算法和应用前景进行了综述,如:不连续Galerkin法、小波方法和格子玻尔兹曼方法等. 本综述有助于读者以整体视角了解地球动力学数值计算方法的发展脉络,并对大数据时代下研究适应日益丰富的数据和新算法提供有益参考.

English Abstract

    • 随着大数据时代的到来,计算地球动力学数值计算方法体系日益完善. 本文系统性地综述了数值计算技术方法在地球科学领域的研究新进展. 2016年Yuen和Knepley在《Treatise of Geophysics》(Zhong et al., 2007)中曾系统总结过在地幔对流研究中使用的数值计算技术. 但在这之后计算领域又出现了诸多新发展,特别是大数据概念提出后,新的数值计算方法不断发展,我们所处的时代也正在经历一场数据革命,而地球科学家正是这场革命的经历者,因此我们衷心希望可以通过本文进一步对这些地球动力学数值计算方法进行全新的总结,为青年科技人员提供更多参考.

      如今,人们经常听到“机器学习”这个词. 但是,机器学习到底是什么?机器学习工具对地球物理学有什么作用呢?简而言之,机器学习(Machine Learning, ML)是一种让计算机以自主的方式去解决给定任务的方法. 尽管计算机没有智能,但它们可以通过主成分分析、聚类分析等方法提取并记住数据的参数化特征,在这一点上远胜人脑. 在构建一个机器学习算法之前,我们首先需要将期望解决的问题进行定义,即希望计算机执行什么任务. 比如,如何找出美国加利福尼亚一个月内发生的地震次数?现有的大多数ML工具都需要大量数据来进行培训和测试,收集尽可能多的数据在机器学习中十分重要,可以是真实的数据(比如照片或地震波形),也可以是仿真产生的数据. 接下来是找到数据中有效的特征. 它们可以是P波和S波的初动,通过地震台站的前10 min波形或模拟裂纹形成/结晶过程的数据. 许多机器学习算法,例如卷积神经网络或主成分分析(Jolliffe, 2011)都可以在训练时自动找到这些特征. 在训练过程中,人们使用各种最优化技术来优化参数(Sioshansi and Conejo, 2019). 训练通常和测试同时进行,机器学习领域有一些专门设计的技术来提高训练的效率,例如神经网络中的反向传播(Hecht-Nielsen, 1992; Jing et al., 2017).

      Zhong等(2007)的论文中仍然强调传统的有限元和有限体积方法,这是因为两位作者(Zhong和Moresi)以及世界上大多数地球动力学家,他们毕生都致力于这些技术. 但是在解决“超级地球”的问题时,可压缩性和绝热加热比地球要强得多(Umemoto et al., 2006; Shahnas et al., 2018; Tsuchiya and Umemoto, 2019; van den Berg et al., 2019),超级计算机的功能无法满足地壳和地幔动力学模型复杂性的要求(Herman et al., 2013; Sizova et al., 2015; Jain et al., 2019),因此在求解高黏流体动量方程时去设计不需要存大矩阵的数值方法十分重要. 本文将回顾一些非传统技术,例如:边界元素、径向基函数、格子玻尔兹曼方法(Mora et al., 2020; Succi and Succi, 2018)和伪瞬态方法(Omlin et al., 2018; Räss et al., 2019)、拟压缩法(Frankel, 1950; Chorin and Others, 1973; Chorin, 1997; Kameyama et al., 2005). 本文认为无矩阵方法将为地球动力学家提供更多的可能,并在高分辨率模型的计算中更有效地节约内存和性能开销. 这对于实现诸如在非线性热—化学—力学相互作用方面的突破性研究,无疑将是十分必要的.

    • 相对于地震科学和地震波传播中的建模方法,固体地球物理现象的数值模拟技术的发展(Stein and Wysession, 2009; Fichtner, 2010; Zobin, 2012; Igel, 2017; Udias and Buforn, 2017; Margrave and Lamoureux, 2019; Shearer, 2019)相对缓慢,新的地球动力学教科书数量也相对较少(Ismail-Zadeh and Tackley, 2010; Glatzmaier, 2013; Gabriele Morra, 2018; Gerya, 2019). 也许是因为地球动力学中要解决的大规模多尺度非线性问题过于复杂,与数据驱动型科学(例如:金融或生物医学领域)相比,近年来新出版的地球动力学教科书相对较少. 但我们仍然认为,为了跟上科学技术发展的脚步,在地球动力学研究中使用更多数据驱动方法十分必要.

      经典的地球动力学教科书(Turcotte and Schubert, 2002)仍然有其教学价值,人们可以将数值结果与解析解(如:考虑剪切生热和热对流情况的剪切流中稳态边界层问题)进行比较,从而调节诸如瑞利数之类的控制参数等. 下面我们从离散化数值计算方法分类、复杂问题的数值计算策略和非线性问题对策等三方面进行综述,并尝试在介绍每个领域方法研究进展的同时,也介绍每个研究领域算法的实现程序或软件,这样有助于帮助读者更加深入理解算法原理和加快解决实际问题.

    • 在地幔对流的数值模拟研究领域,有限元(finite element)、有限差分(finite difference)、有限体积(finite volume)和谱方法(spectral method)有着悠久的历史,这些工作最早可以追溯到半世纪以前,如McKenzie等(1974)以及Machetel和Yuen(1986),而且现阶段在主流学术期刊中仍然大量使用这些方法(Tackley, 2000; Kageyama and Sato, 2004; Harder and Hansen, 2005; Stemmer et al., 2006; John et al., 2018; Phillips et al., 2019).

      一些地球动力学模型采用"从头算"(ab initio)的方法来计算现在的地球或行星内部状态,这类模型不需要真实的约束条件,Tackley等(2000)对于地幔的研究就是使用这类模型的典型案例. 另外一些方法是尝试构造真实的约束条件来开展研究,当然这是极具挑战性的任务. 例如:最新的开源代码库(Geodynamic World Builder)就是用于在直角坐标或球坐标中为地球动力学模型设置合适的初始条件而开发的(Fraters et al., 2019).

      在讨论不同空间离散化的方法之前,我们先从整体角度来回顾一下. 一般来讲,无论使用FD、FV、FE还是 SM方法,任何计算都需要在当前的软硬件基础上来实现,而当前硬件和软件基础架构所能提供的内存相比计算能力来讲更有限,这意味着大规模的计算时间主要由内存带宽和吞吐量决定,因此优化内存使用的算法对提高计算效率更有效,例如:无矩阵方法相比使用大型稀疏矩阵的传统方法更加高效.

    • 该方法建立在微积分的基础上,使用局部泰勒级数展开,得到具有指定精度的计算式. 已经有很多讲解有限差分方法的优秀书籍,例如:Lynch(2004)的有限差分入门著作;Griffiths和Smith(2006)的书中则有许多Fortran 95实现的例子. 有限差分方法具有容易理解和容易编程的优点,自1950年代以来,就一直是数值计算的最主要工具之一,Richtmyer和Morton(1967)为气象学和应用武器研究编写的开创性代码就是很好的证明.

      几乎所有开拓性的地幔对流程序(Torrance and Turcotte, 1971; Turcotte et al., 1973; McKenzie et al., 1974; Parmentier et al., 1975; Jarvis and Mckenzie, 1980)都是采用有限差分方法来设计的. 而为了研究非线性问题,人们使用更高阶精度的有限差分公式,例如:一阶扩展到二阶(Torrance and Turcotte, 1971; McKenzie et al., 1974),用于精确求解二维变黏滞系数问题的四阶公式(Christensen, 1984),Malevsky(1996)开发了一种基于三维样条的地幔对流程序,该程序能够完成高达108瑞利数的模拟计算(Malevsky and Yuen, 1993),如果使用现今的计算机硬件,容许的瑞利数还可以再提高几个数量级.

    • 在传统的有限元方法中,网格区域内的解用一组局部基函数的线性组合(形函数)来描述,这些局部基函数的值仅在一个单元内是非零的. 该方法的目标是找到描述解的最佳线性组合. 一般来说,形函数是一个多项式,所有网格单元互相毗连. 最近还发展出了一种称为无单元伽辽金(Galerkin)的方法,该方法的核函数有相互重叠的部分(Zienkiewicz et al., 2003). 其思路与径向基函数(Radial Basis Functions, RBF)方法类似,结果就是增加了该方法的复杂性,但也比传统的有限元方法更灵活.

      在解决不可压缩介质中的斯托克斯流问题时,为了避免出现虚假的解则建议选择速度阶数高于压力阶数的形状函数(Simo and Hughes, 2000). 单元形状的选择也会影响解的质量,许多有限元建模人员已经观察到了力学模拟中使用三角形网格会出现“网格锁定”的现象. 有限元中的积分通常使用在节点上的线性插值(高斯积分)来计算. 更高级一点的方法是使用物质点上的积分(Sulsky et al., 1994; Moresi et al., 2003),该方法使用起来更加灵活,而且可以考虑介质的非均匀性(例如对于地球内部介质). 而有限元方法的主要优点之一是适用于复杂的几何形状,设置边界条件可以更便捷. 最近,有限元法用于创建俯冲带的参考模型(图1)(Salomon, 2018),可实现热—力耦合和一般的黏塑性(Duretz et al., 2018)程序.

      图  1  有限元法模拟的Central Andes俯冲过程(修改自Salomon, 2018).图中有限元法可以使用非规则网格,灵活设置计算模型,其中展示模型计算50000年之后的结果,分别是(a)位移(km);(b)压力(MPa);(c)热流(mW/m2);(d)温度(K)

      Figure 1.  Central Andes subduction process simulated by finite element method (modified from Salomon, 2018). The finite element method can use irregular grid to flexibly set the calculation model, in which the results calculated by the model after 50000 years are shown, respectively: (a) displacement (km); (b) Pressure (MPa); (c) Heat flow (mW/m2); (d) Temperature (K)

      用于模拟物理上具有不连续的问题的非连续伽辽金法(Discontinuous Galerkin Method, DGM),可看作是传统有限元方法的一种推广(Hesthaven and Warburton, 2007). DGM方法使用不连续的近似解,因此很适合模拟具有不连续的问题,如:核幔边界(Core-Mantle Boundary, CMB)模拟或者一些涉及到部分熔融体的壳内过程模拟,同时不会像有限元或有限体积方法那样会由于吉布斯(Gibbs)现象而引发数值算法不稳定(Cockburn et al., 2012).

      近年来,DGM方法有很多新的研究方向,其中一些对地球物理学问题研究具有重要意义,尤其是对求解Navier-Stokes方程. 例如:重构的不连续伽辽金(RDG)方法,最初是为可压缩欧拉方程设计的,应用到求解Navier-Stokes方程,则是从不连续DGM解出发,使用最小二乘法得到光滑的黏滞性和热通量,从而实现对它们的离散化(Luo et al., 2010).

      最近,一个全新的DGM应用是用于弹性声波的传播和散射建模(图2)(Ye et al., 2016). 在过去的几年里,De Hoop团队在构建该方法的软件包方面做出了重要贡献,使得DGM技术能够在大规模并行计算平台上用于求解三维偏微分方程. 该软件可以在考虑旋转以及高度非球对称的全球模型中计算地球动力学问题,在计算黏弹性正态模态问题时,即使模型介质具有三维横向不均匀以及横向变化的固液交界面,也能得到精确的解(Liu et al., 2019; Shukla et al., 2020). 该程序包也可以处理地球动力学中的非线性特征值问题,例如涉及旋转以及科里奥利力的问题,由于在计算过程中需要使用复数,因此处理这类问题会消耗更多的系统内存.

      图  2  DGM模拟声波—弹性波传播(修改自Hecht-Nielsen, 1992

      Figure 2.  DGM simulates acoustic wave-elastic wave propagation (modified from Hecht-Nielsen, 1992)

    • 有限体积方法与有限元和有限差分方法有许多共同的特点. 在1980年Patankar曾指出,有限体积方法可以看成有限元方法中加权残差的一个特例,只要让加权函数在单元内部是常数,单元外部为零. 与有限差分相比较,有限差分基于泰勒级数,有限体积法则是将方程在控制体积或控制单元上积分,并在单元边界处逼近微分算子(Patankar, 2018).

      至少从1990年代早期开始,有限体积方法就广泛用于地幔对流的数值模拟(Ogawa et al., 1991). Tackley(1993)对二维/三维笛卡儿模型实现了一种标准的多重网格有限体积方法. Ratcliff等(1996)开发了一种用于地幔对流的三维球壳有限体积实现. Harder和Hansen(2005)Stemmer等(2006)以及Choblet等(2007)也编写了有限体积代码计算球壳内的地幔对流,他们使用立方—球网格剖分. Kageyama和Sato(2004)在剖分网格时使用了一种类似棒球的拓扑结构,称为“阴阳”网格,随后在Kageyama和Sato(2004)Yoshida和Kageyama(2004)以及Kameyama等(2008)的工作中进一步完善了该技术,其在许多地幔对流的球形模型中得到了应用(Tackley, 2008). 上述这些地幔对流研究都使用了二阶精度,也有一些研究使用更高阶的公式,并已用于三维直角坐标的地幔对流模拟(Larsen et al., 1997).

    • 谱方法是一种经典的数值方法,主要特点是分析简单、计算准确. 它基于Sturm-Liuoville型二阶微分方程引出的正交本征函数展开式的概念,适用于正交曲线坐标系. 对于地幔对流问题,这意味着我们可以通过笛卡尔几何的傅里叶展开和三维球壳的球谐函数展开来表达水平方向上的数值变化. 这方面有很多优秀的数值计算教材(Glatzmaier, 2013),适用于初学者来了解和学习地球和恒星内部的对流问题.

      Balachandar和Yuen(1994)使用谱变换的方法实现了一套在三维笛卡尔坐标下计算地幔对流的程序,该程序使用Chebyshev函数展开求解竖直方向的两点边值问题,在水平方向上则使用快速傅里叶变换. 借助Krylov子空间迭代技术(Saad and Schultz, 1985),谱方法也被扩展到求解可变黏度的动量方程,可容许103的黏度变化.

      最经典的三维球壳对流代码是由Machetel和Yuen(1986)开发的,通过在径向使用有限差分,沿经向和纬向使用球面谐波分解. Glatzmaier(1988, 2013)则在径向采用Chebyshev多项式展开,在水平向使用基于快速傅里叶变换的球面谐波分解. Chebyshev多项式在热边界层附近可以达到非常高的精度;也可以选择在内部边界层上进行展开,适用于处理包含670 km深度的不连续面的问题.

      Zhang和Yuen(1995, 1996)基于高阶FD方法(Fornberg, 1998)设计了一种适用于三维球形几何的算法,同样是沿圆周方向做球面谐波展开. 他们还开发了一种求解动量方程的迭代技术,用于可压缩对流中横向黏度对比高达200的可变黏度问题(Zhang and Yuen, 1996). 但是谱方法的一个重要缺点是难以解决吉布斯现象的问题,这可以通过黏度随球谐阶数衰减的谱来监测.

      Erlebacher等(2002)也发展了一种谱方法与有限差分结合的方法,主要是在竖直方向上使用Fornberg高阶格式,在水平方向上则使用快速傅里叶变换,并在本世纪初利用MPI并行计算技术实现了高达108瑞利数问题的模拟,网格计算规模达到约600×600×600,如此规模的计算在当时并不多见. 如今的算法比当时有了很大的发展,网格计算规模可以达到1 500×1 500×1 500,瑞利数可以达到1012. Erlebacher等(2002)对这些多尺度特征热湍流的模拟结果进行了细致的可视化工作.

      普林斯顿大学的Jeroen Tromp为谱元法在地球动力学领域的应用做出了很多贡献(Fischer et al., 2002). Tromp团队为科研领域开发了一套SpecFem软件包,可用于研究三维地震波传播、具有横向变化黏度的冰后期回弹和地震波在含有孔隙介质的地壳中的传播等问题. 该软件支持并行计算,通过Github仓库可以很容易获取. 关于谱元法的基本原理,可以参考Komatitsch和Vilotte(1998)的开创性论文. 2016年,在三维地球模型上实现了前所未有的1.2 s地震周期精度的1.8万亿自由度谱元法模拟工作,在82 134个节点的超级计算机上完成了模拟实验(Tsuboi et al., 2016). 伴随着数据同化的趋势,该方法已被用于创建第一代全球联合层析成像(Bozdağ et al., 2016).

      谱元法与无限元(IE)的联合是谱元法的一个新发展方向. 这种单元与边界元有一些类似,都是建立在无界区域的泛函空间上,但是它们比边界元更容易使用,因为不涉及奇异格林函数的积分(Gharti et al., 2019a, 2019b; Gharti and Tromp, 2019).

    • 在极端条件下,FD、FV和FE方法求解与动量方程相关的椭圆方程时,需要较大的内存开销. 比如:Prandtl数趋于无穷大的算例. 当材料参数(如:黏滞系数)变化非常剧烈的时候,可以使用直接法或直接迭代法完成计算,但这类方案都需要借助大型稀疏矩阵来完成,耗费大量的内存和带宽资源. 相反,如果使用无矩阵方法,则可以避免使用大型稀疏矩阵,节约计算资源.

      无矩阵方法的一个典型例子是HyTeG有限元软件框架,使用了无矩阵的多重网格方法,在地球物理学方面有应用(Kohl et al., 2019). 该软件采用低阶多项式逼近变系数偏微分方程的代理元素矩阵. 使用无矩阵算子,最大可以使用47 250个计算核心来模拟大尺度地幔对流问题,应用离散算子后其并行效率高达93%,计算模型自由度高达1万亿,模拟全球网格分辨率可以达到1.5 km(Bauer et al., 2019).

      科研人员还尝试开发无矩阵方法来解决地球物理学中偏微分方程的问题,如Devito项目,由于使用了无矩阵方法,所以无需显式地构建Jacobian矩阵和Hessian矩阵,而是使用对向量的操作代替. 虽然该项目主要关注地震波问题,但同样的框架也适用于求解流体的偏微分方程(Luporini et al., 2020).

      此外,另一种伪暂态方法也使用了一种快速收敛的迭代无矩阵方法,内存占用非常小. 该方法使用改进的物理伪压缩性和物理激励的瞬态项来收敛残差,在5 120个GPU上实现了三维两相流模拟,总体收敛速率达到N4/3Räss et al., 2018; Räss et al., 2019). 无矩阵方法的另一个重要应用是解双曲方程问题,例如:基于分段抛物线法(Piecewise Parabolic Method, PPM)(Colella and Woodward, 1984). 这种方法极少在地球物理学中应用,而在天体物理学中通常用于模拟恒星对流和发电机、超音速流和冲击波,因此在许多地球物理学问题中也具有应用前景,如高压矿物物理学中的冲击波模拟问题. Paul Woodward在研究数年时间尺度上恒星核燃烧的问题时采用了无矩阵方法,实现了4 0003规模的网格计算.

    • 在各种工程问题中,边界单元方法(Boundary Element Method, BEM)的应用已经具有悠久的历史. 特别是在声学方面,许多应用都是专门为地震学而研发(Guzina and Pak, 2001; Dangla et al., 2005; Semblat, 2011; Chaillat and Bonnet, 2013). BEM的大多数用途是用单个边界单元层分隔两个计算区域. 在每个边界元上,积分运算可以采用解析公式(Salvadori, 2010)或数值积分(Pozrikidis, 2002)来计算. 最近很多对多圈层地质问题领域的新应用(Chaillat et al., 2012)表明,如果能适当加速,边界单元方法有望成为地球物理学中的主要数值工具之一.

      在地球动力学领域,有的研究将BEM与FEM联合使用(Morra and Regenauer-Lieb, 2006),也有单独使用的情况,这些研究包括俯冲板片的动力学建模(Morra et al., 2007, 2009; Li and Ribe, 2012; Li et al., 2014; Gerardi and Ribe, 2018;),全球尺度的板块构造(Morra et al., 2009)和俯冲带系统上的能量估计问题(Gerardi et al., 2019),这是地球动力学领域的一个基础性科学问题(Conrad and Hager, 1999; Buffett and Rowley, 2006; Capitanio et al., 2009; Buffett and Becker, 2012).

      在火山学研究方面,Manga和Stone(1995)率先将BEM方法引入,后续有很多工作跟进,同时,快速多极加速方法(Fast Multipole Acceleration, FMA)也被用来做大规模非均匀系统的建模,模拟了在地幔中的板片下插或者底辟构造形成的羽流上升问题(Morra et al., 2015). 在地震学研究中,BEM方法常用于模拟断层滑动引起的应力/应变过程(Mai et al., 2016). 最近,Thompson和Meade(2019a)还新开发了一种使用连续线性位移和滑移基函数的无奇异伽辽金边界元法,以弥补线性边界单元边界上奇异点积分所产生的误差,并用于研究Cascadia俯冲带上的周期性地震活动(Thompson and Meade, 2019b).

      BEM中的线性系统是适定的,但其关联的矩阵是稠密的. 在过去10年中,多种无矩阵方法被用于解决边界单元问题,包括:快速多极加速(FMA)方法(Barnes and Hut, 1986; Tornberg and Greengard, 2008)、多级矩阵方法(Benedetti et al., 2008; Cazeaux and Zahm, 2015)等. 第一种方法具有处理多尺度问题的潜在优势,因为它与三维非结构化表面网格兼容,并且网格分辨率可以在需要加密的区域自适应加密. 构建好的方程组可以使用GMRES迭代算法求解(Saad and Schultz, 1985). Hi-BoX就是一个边界单元法软件包,其使用了目前最先进的快速直接解法和迭代解法,并包含H矩阵和FMA方法以及将二者结合用于大规模并行异构体系下的解决方案(Abboud and Barbier, 2016).

    • 格子玻尔兹曼方法(Lattice Boltzmann Method, LBM)处理的问题尺度介于分子和连续介质之间,这种方法最初被设计用于模拟计算流体力学(Computational Fluid Dynamics, CFD)问题,现在被广泛应用于多种问题(Succi and Succi, 2018). 与其他CFD方法不同的是,该方法不是直接求解Navier-Stokes方程,而是通过计算粒子在固定网格上的流动和碰撞(弛豫)来求解,其中用到统计力学中平衡态扰动解导出的函数. 这种方法有许多优点,其中之一就是软件可以写得非常简洁、有效(Mora et al., 2020). 诸如汽液共存、液滴、多孔介质复杂环境中的流体等许多系统都可以采用LBM方法模拟.

      LBM方法起源于统计物理学中的一种叫格子气体自动机(Lattice Gas Automata,LGA)方法,该方法用于模拟粒子在离散格点上的运动和碰撞,以重现分子在气体中的运动规律. LGA在宏观极限下符合Navier-Stokes方程(Frisch et al., 1986),该方法计算量非常大. 而在LBM方法中,仅计算晶格中参与运动和碰撞的粒子的分布. 换句话说,LBM模拟的是经典玻尔兹曼方程的一阶近似.

      Bhatnagar、Gross和Krook (BGK)(Bhatnagar et al., 1954)提出了一种通过松弛来计算粒子碰撞的有效方法,并被用于高效地实现LBM模拟(Higuera and Jiménez, 1989; Qian et al., 1992; Reis, 2003),以及施加精确的压力和速度边界条件(Zou and He, 1997). 后续对LBM方法应用和研究得到了飞速的发展,有更多的相关文献可以参考(Huang et al., 2015b; Krüger et al., 2017).

      使用LBM方法的很多研究是模拟热对流问题(Shan, 1997; He et al., 1998; Guo et al., 2002; Wang et al., 2013),包括多相流方法(Shan and Chen, 1993),到目前为止仍然算是一个非常活跃的研究领域(Huang et al., 2015; Di Ilio et al., 2017; Xie et al., 2017). 由于LBM方法能很好地在高性能计算集群(High Performance Computing, HPC)上部署,该方法在许多科学和工程领域内得到广泛应用. 最近的大规模并行多重物理过程模拟软件WaLBerla(Feichtinger et al., 2011) 能在千万亿次HPC上计算1012节点规模的问题,其万亿规模的非均匀网格在Schornbaum和Rüde(2016)的文章中有详细描述. 广泛使用的Palabos(Lagrava et al., 2012)、OpenLB(Heuveline and Latt, 2007)、 LB3D(Heuveline and Latt, 2007; Schmieschek et al., 2017)系统,通过HPC测试都显示具有较好的计算性能.

      尽管LBM方法在热对流方面取得了很多研究进展,但最近才被应用于研究地幔对流类地球动力学问题(Mora and Yuen, 2017, 2018). 该方法的Python简版可以参考(Mora et al., 2020). 对于地球物理学中的其它研究方向,也有一些研究进展,特别是多孔介质内流体的模拟问题(Leclaire et al., 2017; Zhao et al., 2019),具体地说如反应运移(Kang et al., 2006)、页岩气(Middleton et al., 2015)、复杂流(Schaefer et al., 2013)、多相流(Pradana and Fauzi, 2019)和水文学(Chen et al., 2018)等的模拟. 值得注意的是,最近在火山声学(Brogi et al., 2018)、岩浆脱气(Parmigiani et al., 2016; Parmigiani et al., 2017)等方面也出现了一些应用案例. Zhu等(2020a, 2020b)使用4 0003格点规模的LBM方法求解MHD方程,模拟了太阳内部的磁流体运动.

    • 对于线性或非线性斯托克斯流(Stokes flow),上述方法(FD、FV、FE)中的刚度矩阵是对称正定的,这使得很多方法可以用来解决整个线性问题. 多重网格法(MG)是地幔对流模拟中最常用的方法. 该方法早期的实现可以追溯到Tackley(1993)、CITCOM (Moresi and Solomatov, 1995)及其后续的软件包CitcomS/CitcomCU(Zhong et al., 2000). 另外,阴阳网格方法也是一个重要进展,该方法让多重网格法在球坐标中的实现成为可能(Yoshida and Kageyama, 2004; Tackley, 2008).

      多重网格法的核心思想是在一组多重嵌套的不同粗细的网格上来解同一个问题,首先在较粗尺度的网格上快速获得近似解,此时的残差在较粗尺度上达到极小,然后再降低细尺度上的残差,这一步通过在细尺度的网格上的局部求解器来实现(也称为“平滑器”). 该方法的一个特点是,随着问题规模的增加,迭代的次数不会显著增长. 逐点平滑的策略比较常见(如高斯—赛德尔迭代),但对于复杂的方程组,更复杂的平滑方法是十分必要的. Yavneh(2006)Oosterlee和Gaspar-Lorenz(2006)Gmeiner等(2015)对多重网格法的许多细节进行了综述. LaMeM最近尝试将广义多重网格方法移植到用于岩石圈—地幔问题的有限差分方程(Popov and Kaus, 2013).

      MILAMIN(Dabrowski et al., 2008)是一套紧凑的多重网格代码,仅适用于拉普拉斯算子,然后使用Schur补的方法将其用于求解Stokes方程(May and Moresi, 2008). 因为此处是应用于较简单的拉普拉斯算子,所以使用的是逐点平滑. 同样的方法可以用于有限差分求解Stokes方程(Räss et al., 2017).

    • 对地球系统进行建模涉及到多尺度物理,这也是地幔对流的一个重要特征. 地幔对流的特征波长可以很长,例如:太平洋板块的尺度为10 000 km. 然而,地幔对流也从根本上受薄的热边界层以及被称为海洋岩石圈的边界层控制,在高瑞利数的情形下会导致小规模的上升流(100 m). 板块边界过程和热化学对流过程发生的尺度可能更小. 此外,材料属性受到近微观属性的影响(如:晶粒大小),除了少数尝试外,在地幔对流研究中还无法很好地考虑这类问题(Rozel et al., 2015; Schmalholz and Duretz, 2017).

      大多数现有的地幔对流编码(如:CitcomS和STAG3D)都适用于基本一致的网格,迄今为止最大的计算规模使用了数亿个未知数(Jadamec and Billen, 2010). 还有一些更复杂的代码,如:Underworld (Asgari and Moresi, 2012)和SEPRAN(Van Kan et al., 2008; van den Berg et al., 2015). SEPRAN不仅适用于地球动力学,还可用于解决工业问题,例如具有惯性项和复杂流变性(如塑性)的车辆碰撞. 这两套代码均使用Python语言编写,因此,相比严格用C或FORTRAN编写的代码更方便学习使用. 最近的代码如ASPECT(Kronbichler et al., 2012)和Rhea(Burstedde et al., 2008; Burstedde et al., 2013)使用了AMR技术(Davies et al., 2011; Leng and Zhong, 2011; May et al., 2013). 它们可以运行在大规模并行计算机上(103~105的核),在三维全球模型中其局部分辨率甚至能达到数千米. 因此,它们为解决一些计算上具有挑战性的地球动力学瞬态问题提供了可能(Stadler et al., 2010).

      AMR最适合处理具有局部非均匀性的问题,而不擅长捕捉整个区域的小尺度特征. ASPECT和Rhea不能处理地球动力学中的所有问题,例如孔隙介质中的波动问题. 其他更灵活的有限元软件包如芝加哥大学R. Scott开创的Fenics(Dupont et al., 2003),正在被许多大学使用,哥伦比亚大学的M. Spiegelman和卡内基研究所的C. Wilson用它处理多物理场孔隙波问题.

      二维热—化学对流AMR(Leng and Zhong, 2011)为数值分析带来了极大的拓扑和数学复杂性. 它需要高效的网格细化、粗化、重分区算法,还需要高效的并行方法在多处理器之间做好负载均衡. 为了应对这些挑战,Rhea使用了一种称为八叉树森林的技术,计算区域被划分为具有八叉树结构的子域(Burstedde et al., 2008; Burstedde et al., 2013). 八叉树是一种逻辑结构,其中每个节点要么是没有后代节点的叶节点,要么是有后代节点的中间节点,这种逻辑结构在计算机科学中早已有所应用,例如:PARAMESH(MacNeice et al., 2000). 该结构是递归的、层次化的,可以有效地实现FE分析和并行计算,就像在Rhea中一样. ASPECT是在Deal.II数值软件包基础之上开发的. 它支持高阶时空离散化方案,并使用八叉树森林进行并行网格剖分.

      处理千米分辨率的含时间项地幔对流问题的一个挑战是数值稳定需要较小的时间步长. Rhea和ASPECT在求解能量方程时都使用显格式来离散化时间项. 到目前为止,ASPECT尚未用于计算具有复杂物理机制的时间依赖性问题,例如:需要计算数百万个时间步的晶粒尺寸减小问题,这类问题的求解难度很大. 相反,ASPECT的很多应用都少于10 000个时间步,也没有用到它最复杂的功能.

      AMR算法也有局限性,它可能会导致某种类型的“过度聚焦”,在解决高度非线性且有局部形变的问题时,AMR算法可能在某一点处不断优化网格而忽略其它区域的网格加密. 另外,除了复杂性和性能问题,目前尚未有严格的算法来解决何时何处需要重新优化网格,这仍需要建模者自己把握.

    • 由于岩石具有孔隙性,地下流体流动需要涉及一个以上的相,化学反应过程十分复杂且经常发生,但这些在研究壳内近地表地质问题中具有重要意义. 这方面模拟需要的基本方程详见Phillips(2005)文献中的阐述. 在推导描述多相流体问题的偏微分方程时,为了保证方程的普适性,包含热—力学—化学耦合作用十分重要(Yarushina and Podladchikov, 2015; Räss et al., 2019).

      地壳内的许多地质过程涉及到广泛分布的流体的排出,例如:在沉积过程中滞留的流体或部分熔融产生的流体. 孔隙波在涉及多相流的许多问题中出现,有时甚至超过三相问题. 在这些问题中,渗透率是一个很重要的动力学参数,它由不同密度的岩石与流体的流变性和动力学不稳定性协同控制.

      需要注意的是在孔隙波方程中我们将孔隙度作为自由变量,而渗透率作为达西定律中的非线性输运参数,则看作孔隙度的函数. 它们比黏性流体的斯托克斯方程更复杂、非线性更强. 渗透率是孔隙度的高度非线性函数,就像黏滞系数是岩石流变学中温度和粒度的高度非线性函数一样. 孔隙波不同于声波、弹性波和孔隙弹性波. 孔隙波以其特有的物理性质在地壳中传播,传播速度从人类游泳的速度到高速公路上行驶的汽车的速度范围变化. 因此,它比声波、弹性波和孔弹性波慢许多个数量级,因为浅层地壳岩石的孔隙度随时间的变化而变化,这是由于岩石基质内的流体运动造成的. 这种相互作用可能会产生一种水文机制,在这种机制中,流动是由充满流体的孔隙度和类波运动的自传播域来完成(McKENZIE, 1984; Fowler, 1985; Scott and Stevenson, 1986; Spiegelman et al., 1993).

      对于两相或多相的压实,人们可以用以下几种流变学来建立多种类型的方程: 黏性、弹性、黏弹塑性. 它们代表由多孔岩石基质和黏性较低的间隙流体组成的系统. 为了求解它们,计算成本随着流变性质的复杂度而增加,因此,必须开发新的计算方法,例如:使用无矩阵伪瞬态方法(Omlin et al., 2018; Räss et al., 2019). 由Ludovic Räss和Sam Omlin设计的孔隙波三维模拟可以被认为是2020年计算地球动力学领域的壮举,它使用了洛桑大学300个GPU组成的集群. 到目前为止,这些十分困难的模拟工作还没有在其他地方实现.

      脱水引起的孔隙波也与偶发的地震和断层滑动有关. Skarbek和Rempel(2016)建模了一个一维的俯冲板片,让它俯冲到慢滑动和震颤区域,结果表明周期性的慢滑动和震颤现象可以通过孔隙波在板片界面上的迁移来解释. Yarushina等(2017)也基于多孔岩石的传播变形,讨论了流体注入在远离作业场地的地方诱发微震活动的机理.

    • 当在柔性变形的本构方程中同时考虑弹性和黏性时,变形引起的剪切生热可能导致热失稳. 弹性提供了驱动热失稳的势能. 为了捕获剪切生热失稳的较短时间尺度和收缩区域,必须使用高分辨率网格和特殊的时间步长(So et al., 2013),包括使用了隐式和显式的自动化方案. 可以证明,在非绝热条件下,黏度对温度的简单Arrhenius依赖性就足以产生热失稳. 此外,非绝热导致应变和温度分布的极端局部化,可能导致自然界中经常观察到的窄剪切带(Braeck et al., 2009).

      从韧性区域地震的产生(Thielmann et al., 2015)到俯冲带的形成(Regenauer-Lieb and Yuen, 1998; Regenauer-Lieb et al., 2001),地球物理学领域的许多应用都与剪切加热和热失稳过程有一定的联系,但这一过程的诸多细节和普适判据仍然存在争议(Prieto et al., 2013).

      Kiss等(2019)的最新研究通过热软化机制阐明了韧性剪切带产生的机理. 他们通过比例分析找到了局部化发生的判据,利用从一维到三维的数值模拟来支持其研究结果,并将其应用于岩石圈相关研究(Duretz et al., 2019). 时间步长的判据主要需要考虑与动量方程无关的方程,即热方程. 在这类问题中,源于简单流体力学的Courant时间步长准则必须根据此处的物理过程重新定义,例如从温度方程,或孔隙波中出现的孔—黏—弹性流变学中弹性强度的Deborah数准则.

    • 随着考虑更真实、非线性更强的流变行为,如:橄榄岩的非牛顿剪切减薄,或地幔结构的塑性变形,传统的基于简单牛顿迭代的非线性求解方法变得十分低效. 而非线性迭代的难以收敛,且不能通过线性预条件子来改善,因此研究人员不得不使用耗时但稳定的Picard迭代. 最近人们发明了一类非线性预条件子(nonlinearity PreConditioner, NPC)(Brune et al., 2015),特别适合多物理问题. 与线性预条件子相似,NPC也可以自然地分为“左”或“右”预条件子(Hwang and Cai, 2005).

      左非线性预条件子引入一些小规模的非线性问题,这些问题比原问题能更高效的求解,同时也能很好地收敛到原问题的解. 该方法的区域分解版本,如ASPIN(Cai and Keyes, 2002)解决每个子域上的小型非线性问题,并将它们相加或相乘地组合起来. 小型非线性问题也可以使用Krylov子空间生成,如非线性GMRES(Oosterlee and Washio, 2000)或非线性Broyden(Fang and Saad, 2009)算法,其中小型非线性最小二乘问题在子空间中解决,并用于更新原问题的解.

      另一方面,右非线性预条件子则像右线性预条件子的作用方式一样,直接改变近似解所在的空间,通过这种方式改善非线性. 这里最简单的形式是使用一个辅助求解器(如:Picard)来引导Newton迭代. 更复杂的版本,如非线性消除法(Lanzkron et al., 1996; Hwang et al., 2015),在每次迭代中精确地求解非线性问题的子集,这种方法类似于线性块消除或Schur补运算.

      非线性预条件框架现在已经有可用的版本(Balay et al., 2016),允许用户组合左和右非线性预条件子,以加快问题的收敛. 由于所有的非线性预条件子都是近似的非线性求解器,这相当于非线性迭代组合的一般理论(Brune et al., 2015),用户可以挑选特定的求解器来处理具体的非线性问题. 特别是对于多物理问题,一个领域的非线性可以被其他领域主导的Jacobian方向掩盖,导致收敛速度停滞(Klawonn et al., 2017).

      一些迭代求解器还可以同时减少线性和非线性残差. 如果解不依赖于一致的切线或不符合完整的牛顿框架,那么线性问题就不需要在远离平衡的地方精确求解. 对非线性项采用单层循环松弛的伪瞬态求解器已经被用来模拟剪切加热(Duretz et al., 2019)和水—力学耦合(Räss et al., 2019)引起的应变局部化.

    • 小波方法的出现最早可以追溯到1980年代. 地球物理学家Alex Grossman和Jean Morlet最初应用该方法捕捉地震波序列中的尖峰. 完整的小波发展史可以参考《A Wavelet Tour of Signal Processing》(图3)(Mallat, 1999).

      图  3  小波方法表示的二维图像(Mallat, 1999). (a)原始图像f(N),分辨率N=256×256;(b)原始图像用正交小波分解的系数矩阵<f,Ψk.*?>=>j,n>,k=1,2,3,4,小波尺度因子为2j,大于阈值T的系数用黑点表示;(c)使用尺度因子最大的3组小波系数(N/16个)重建图像;(d)使用最大的N/16个系数重建图像

      Figure 3.  The wavelet method represents a two-dimensional image (Mallat, 1999). (a) Original image f(N), resolution N=256×256; (b) the original image using orthogonal wavelet decomposition coefficient matrix <f,Ψk.*?>=>j,n>, k = 1, 2, 3, 4, wavelet scale factor is 2j, is greater than the threshold value of T coefficient with black spots; (c) Three sets of wavelet coefficients (N/16) with the largest scale factor were used to reconstruct the image; (d) Reconstruct the image using the maximum N/16 coefficients

      小波方法被用于求解非线性微分方程,该方法的优势在于可以用很少的参数来捕捉信号中的多尺度特征. Felix Hermann还应用小波方法开发了很多在地震机器学习方法的应用(https://slim.gatech.edu),例如:Jittered Under-sumpling(Hennenfent and Herrmann, 2008)应用Curvelets处理地震波形(Hennenfent and Herrmann, 2006). 值得注意的是,虽然这些使用小波的机器学习应用程序都需要大量的计算资源,但是由于小波具有紧凑的多尺度特性,它们在CNN和RNN神经网络上非常高效. 在一些研究的例子中,我们可以看到小波方法用来捕捉火星岩石圈导纳函数的横向变化(图4)(Kido and Yuen, 2003)和地球表面的地形特征(Audet and Mareschal, 2007).

      图  4  火星表面重力和地形的小波展示(修改自Kido and Yuen, 2003).(a)原始重力;(b)原始地形;(c)不同尺度小波重建的重力;(d)不同尺度小波重建的地形. l.*?>=>w=8、16、32、64对应的尺度分别是424 km、212 km、106 km、53 km

      Figure 4.  A wavelet display for the gravity and topography of Mars (modified from Kido and Yuen, 2003). (a) Original gravity; (b) Original topography; (c) The gravity reconstruction by wavelets of different scales; (d) The topography reconstruction by wavelets of different scales. The corresponding scales of l.*?>=>w=8, 16, 32 and 64 are 424 km, 212 km, 106 km and 53 km respectively.

    • 本文系统地回顾了经典的数值计算方法,包括:空间离散化方法、谱方法、无矩阵方法、小波方法和一些非线性过程的模拟技术,与Zhong等(2007)的综述相比,本文更加深入讨论了一些技术细节,如:孔隙波和剪切生热的一些细节,同时,涉及了更多种类的数值方法介绍. 文中同时也兼顾介绍了一些新的数值计算技术发展与应用进展,如:格子波尔兹曼方法以及它在地幔对流(Mora and Yuen, 2018; Mora et al., 2020)和磁流体动力学(Zhu et al., 2020a, 2020b)中的应用.

      本文旨在从整体视角给出经典数值计算方法的发展脉络,使读者了解这些发展的演化过程,并未对每种方法的技术细节进行详细的剖析. 我们建议感兴趣的读者们能进一步通过查阅本文的参考文献来了解更多技术细节,通过这些文献可以更好地理解相关研究的发展趋势以及为快速地掌握相关方法原理提供有益的帮助. 本文内容能为从事地球动力学建模的科研人员提供更多的信息,全面了解数值计算技术的新进展.

参考文献 (208)

目录

    /

    返回文章
    返回