An overset-grid finite-difference algorithm for seismic wavefield propagations modelling in the polar coordinate system with a complex free-surface topography

SUMMARY

Nowadays, finite-difference method has been widely applied to simulate seismic wavefield propagation in a local seismic model with a complex topography by utilizing curvilinear grids and traction image method in Cartesian coordinates system. For a global seismic model in polar coordinate system, it is still a challenge for the conventional finite-difference method to deal with the grid singularity at the centre and the topographic free surface at the top. In this study, we develop a finite-difference method in polar coordinates for seismic wavefield propagation in the 2-D global model. In the proposed finite-difference method, the overset-grid algorithm is used to handle the grid singularity at the centre, and the curvilinear grid technique as well as the traction image method are applied to implement free-surface boundary condition on the complex topography. The proposed finite-difference method is validated in flat and topographic free-surface models by comparing synthetic waveforms with reference solutions. The seismic wavefield propagation in a realistic Mars profile is computed by the proposed finite-difference method. The proposed finite-difference method is an efficient and accurate method for seismic wave modelling in the polar coordinate system with a complex free-surface topography.

Numerical modelling, Planetary seismology, Computational seismology, Wave propagation

Subject

Seismology

Issue Section:

Seismology

1 INTRODUCTION

Seismic wavefield propagation simulations are important for full waveform inversions and global seismology studies. Many numerical methods have been applied to simulate seismic wavefield propagation, such as finite-element method (FEM; e.g. Lysmer & Drake 1972), pseudospectral method (PSM; e.g. Tessmer et al. 1992), spectral element method (SEM; e.g. Komatitsch et al. 2001) and Discontinuous Galerkin method (DG; e.g. Etienne et al. 2010; Pelties et al. 2012). Among these methods, the finite-difference method (FDM) is a popular algorithm for several reasons. The operator of FDM only involves several adjacent gridpoints, which is ideal for parallel computing. The character in parametrize is well suited for highly heterogeneous media and this inhomogeneity can be well mapped into grid. It does not require meshes to integrate interior interface topography which minimizes the effort in grid generation.

The complex free-surface topography can greatly affect seismic wavefield propagation. The key factors for FDM are the grid generation and free-surface boundary conditions. Compared with the numerical algorithm adopting unstructured grids, the conventional FDM of rectangular grids is not fit for complex topographies. Several novel FDM have been developed for free-surface topographies. Ohminato & Chouet (1997) approximate the surface topography by the stair-case approximation and implement the free-surface boundary conditions by a vacuum approximation. Oprsal & Zahradnik (1999) extend the vacuum method and the stair-case approximation to the second-order wave equation. In addition, Jih et al. (1988) approximate free-surface topography with segmented lines between gridpoints. Hestholm & Ruud (1994, 2000) adopt a vertical transformation grid technique in the staggered grid finite-differnce scheme to approximate the free-surface topography. The algorithms for correctly implementing free-surface boundary conditions have been greatly developed in recent years. In FDM, the normal derivatives of velocity components on free-surface gridpoints are usually computed by the horizontal derivatives of velocity components through the free-surface boundary condition (Zhang et al. 2012). One way to compute the normal derivatives of stress components on the free-surface gridpoints is the one-sided finite difference and it can maintain a second-order spatial accuracy in complex media. Another way is the stress image method which is usually for flat surface model (Robertsson 1996). Zhang & Chen (2006) use curvilinear grid to describe the free-surface topography and traction image method to implement the free-surface boundary condition.

For the seismic wavefield modelling in a global model, the earth's curvature should not be omitted. It has been approximated by grid deformation (Appelö & Petersson 2009) or level transformation (Li et al. 2014). It is natural in the 2-D polar coordinate system or 3-D spherical coordinate system (Toyokuni & Takenaka 2006) to consider the Earth's curvature. The staggered-grid finite-difference scheme is widely used in seismic wave modelling, where different components of the velocity and stress vectors are defined at different positions of the rectangular cell. Several studies have used this scheme for seismic wave modelling in spherical coordinates. However, the FDM of staggered grids is only suitable for solving elastic dynamic equations with a Cartesian mesh. In the polar coordinates, there is no efficient way to distribute the wavefield components in a staggered way while each component can only be distributed in its defined position. So, a collocated grid is used for simulation in this study.

The overset-grid algorithm (e.g. Benek et al. 1983) is a successful high-quality multiscale grid generation method in computational fluid dynamics. In this algorithm, the large-scale grid is first decomposed into small-scale subgrids in terms of physical, geometric and computing properties of the numerical model, and then subgrids are connected by overset-grids in the boundary areas. The overset-grid algorithm helps transform the grid generation problem of a complex geometry into a relatively easier-solved structured grid problem (Yang & Stern 2009), and makes sure that the subgrids are simpler, more flexible and independent, and higher-quality than the undecomposed large-scale grid.

There are lots of works that have introduced overset-grid algorithm into calculations in seismology. The overset-grid algorithm has been introduced to seismology for scalar wavefield propagation modelling (e.g. Beznosov & Appelö 2020) in the DG method. Chaljub et al. (2003) use the SEM to simulate seismic wave propagation in full-earth model, utilizing ‘cubed sphere’ transformation to generate a cube grid at the centre of earth and using Lagrange interpolation to complete the simulation of seismic waves. Komatitsch & Tromp (2002) use the SEM to simulate global propagation of seismic waves and places a cube in place of core to avoid grid singularity. Zang et al. (2021) extend the overset-grid algorithm in Cartesian coordinate system for elastic wavefield propagation modelling with complex topographies in the curvilinear grid finite-difference technique to simplify the grid generation and improve computing efficiency. In addition, apart from the overset-grid algorithm, some other algorithms are also used for global-scale seismic wave simulation. Nissen-Meyer et al. (2007) proposed a 2-D SEM for the axisymmetric Earth model to compute the response to a full 3-D moment tensor source. Wang et al. (2001) propose the PSM under polar coordinates and develop a scheme that uses extension of field variables in the radial direction, with which computation of the wavefield at the centre is avoided.

In this work, we first extend the curvilinear-grid finite-difference method (CGFDM, Zhang & Chen 2006) of Cartesian coordinate system into 2-D polar coordinate system to simulate seismic wavefield propagation in an annular or circle shape area with a complex topography. The overset-grid algorithm is utilized to solve the grid angularity in the centre areas and improves the parallel computing accuracy. The proposed overset-grid CGFDM is validated by comparing with reference solutions in global flat and topographic models. Topography effects in a realistic Mars profile are also simulated and demonstrated in this study.

2 METHOD

2.1 The first-order velocity-stress equation in the Polar coordinates

In seismology, the first-order velocity-stress equation is always utilized to simulate elastic wavefield propagation. It can be derived from the momentum equation and stress–strain relation,

\begin{array}{r} ρ \partial_{t} v = \nabla \cdot σ + f, \end{array}

(1)

\begin{array}{r} σ = C : ε \end{array}

(2)

where $ρ$ is the density, v is the velocity vector, $σ$ is the stress tensor, $ε$ is the strain tensor, C is the elastic tenor and f is the body-force source. In the elastic isotropic media, the elastic tensor is simplified as the Lamé constants $λ$ and $μ$ ⁠. In the polar coordinate system, the first-order velocity-stress equation can be expressed as

\begin{array}{r} ρ \partial_{t} v_{r} = \partial_{r} σ_{r r} + \frac{1}{r} (σ_{r r} + \partial_{θ} σ_{r θ} - σ_{θ θ}), \end{array}

(3)

\begin{array}{r} ρ \partial_{t} v_{θ} = \partial_{r} σ_{θ r} + \frac{1}{r} (σ_{r θ} + σ_{θ r} + \partial_{θ} σ_{θ θ}), \end{array}

(4)

\begin{array}{r} \partial_{t} σ_{r r} = (λ + 2 μ) \partial_{r} v_{r} + \frac{λ}{r} (v_{r} + \partial_{θ} v_{θ}), \end{array}

(5)

\begin{array}{r} \partial_{t} σ_{θ θ} = \frac{(λ + 2 μ)}{r} (v_{r} + \partial_{θ} v_{θ}) + λ \partial_{r} v_{r}, \end{array}

(6)

\begin{array}{r} \partial_{t} σ_{r θ} = μ \partial_{r} v_{θ} + \frac{μ}{r} (\partial_{θ} v_{r} - v_{θ}), \end{array}

(7)

where r is the polar diameter and $θ$ is the polar angle. In order to apply the curvilinear grid technique into eqs (3)–(7) for free-surface topographies, eqs (3)–(7) in the physical domain (r, $θ$ ⁠) should first map to the computing domain (⁠ $ξ$ ⁠, $η$ ⁠) by the following formulae (Fig. 1a),

\begin{array}{r} ξ = ξ (r, θ), \end{array}

(8)

\begin{array}{r} η = η (r, θ), \end{array}

(9)

to obtain the following equations,

\begin{array}{r} ρ \partial_{t} v_{r} = \partial_{ξ} σ_{r r} \partial_{r} ξ + \partial_{η} σ_{r r} \partial_{r} η + \frac{1}{r} (σ_{r r} + \partial_{ξ} σ_{r θ} \partial_{θ} ξ + \partial_{η} σ_{r θ} \partial_{θ} η - σ_{θ θ}), \end{array}

(10)

\begin{array}{r} ρ \partial_{t} v_{θ} = \partial_{ξ} σ_{r θ} \partial_{r} ξ + \partial_{η} σ_{r θ} \partial_{r} η + \frac{1}{r} (2 σ_{r θ} + \partial_{ξ} σ_{θ θ} \partial_{θ} ξ + \partial_{η} σ_{θ θ} \partial_{θ} η), \end{array}

(11)

\begin{array}{r} \partial_{t} σ_{r r} = (λ + 2 μ) (\partial_{ξ} v_{r} \partial_{r} ξ + \partial_{η} v_{r} \partial_{r} η) + \frac{λ}{r} (v_{r} + \partial_{ξ} v_{θ} \partial_{θ} ξ + \partial_{η} v_{θ} \partial_{θ} η), \end{array}

(12)

\begin{array}{r} \partial_{t} σ_{θ θ} = \frac{(λ + 2 μ)}{r} (v_{r} + \partial_{ξ} v_{θ} \partial_{θ} ξ + \partial_{η} v_{θ} \partial_{θ} η) + λ (\partial_{ξ} v_{r} \partial_{r} ξ + \partial_{η} v_{r} \partial_{r} η), \end{array}

(13)

\begin{array}{r} \partial_{t} σ_{r θ} = μ (\partial_{ξ} v_{θ} \partial_{r} ξ + \partial_{η} v_{θ} \partial_{r} η) + \frac{μ}{r} (\partial_{ξ} v_{r} \partial_{θ} ξ + \partial_{η} v_{r} \partial_{θ} η - v_{θ}) . \end{array}

(14)

Figure 1.

(a) Diagram of curvilinear-grid transformation, convert variables from physical space to computational space. (b) Diagram of traction image method in polar coordinate system.

Open in new tab Download slide

2.2 Implement the free-surface boundary condition by the Traction image method in the polar coordinates

The traction image method (Zhang & Chen 2006; Sun et al. 2016) is a high-accuracy algorithm to implement the free-surface boundary condition in the curvilinear grid (Fig. 1b). The free-surface boundary condition is

\begin{array}{r} T^{η} = n_{η} \cdot σ = 0, \end{array}

(15)

where $T^{η}$ is the traction, $n_{η}$ is the vertical unit normal vector at the free surface. Without generality, the free surface is assumed to be along the polar angel $θ$ direction, and the vertical unit normal vector $n_{η}$ can be expressed as

\begin{array}{r} n_{η} = \frac{\nabla η}{| \nabla η |} = \frac{\partial_{r} η {\hat{e}}_{r} + \frac{1}{r} \partial_{θ} η {\hat{e}}_{θ}}{| \partial_{r} η {\hat{e}}_{r} + \frac{1}{r} \partial_{θ} η {\hat{e}}_{θ} |} . \end{array}

(16)

The free-surface boundary condition in the curvilinear-grid polar coordinate system is obtained by substituting eq. (16) into eq. (15),

\begin{array}{r} \frac{σ_{r r} \partial_{r} η + \frac{1}{r} σ_{θ r} \partial_{θ} η}{| \partial_{r} η {\hat{e}}_{r} + \frac{1}{r} \partial_{θ} η {\hat{e}}_{θ} |} = 0, \end{array}

(17)

\begin{array}{r} \frac{σ_{r θ} \partial_{r} η + \frac{1}{r} σ_{θ θ} \partial_{θ} η}{| \partial_{r} η {\hat{e}}_{r} + \frac{1}{r} \partial_{θ} η {\hat{e}}_{θ} |} = 0. \end{array}

(18)

The relationship of the spatial derivatives of velocity components in the polar coordinate system can be obtained by substituting eqs (12)–(14) into equations of temporal derivatives of eqs (17) and (18), and it can be eventually expressed in the matrix form as

\begin{array}{r} A [\begin{array}{c} \partial_{η} v_{r} \\ \partial_{η} v_{θ} \end{array}] = B [\begin{array}{c} \partial_{ξ} v_{r} \\ \partial_{ξ} v_{θ} \end{array}] + C [\begin{array}{c} v_{r} \\ v_{θ} \end{array}], \end{array}

(19)

where

\begin{array}{r} A = [\begin{array}{c} \frac{μ}{r^{2}} \partial_{θ} η \partial_{θ} η + (λ + 2 μ) \partial_{r} η \partial_{r} η & \frac{λ}{r} \partial_{r} η \partial_{θ} η + \frac{μ}{r} \partial_{r} η \partial_{θ} η \\ \frac{λ}{r} \partial_{r} η \partial_{θ} η + \frac{μ}{r} \partial_{r} η \partial_{θ} η & μ \partial_{r} η \partial_{r} η + \frac{(λ + 2 μ)}{r^{2}} \partial_{θ} η \partial_{θ} η \end{array}], \end{array}

(20)

\begin{array}{r} B = - [\begin{array}{c} \frac{μ}{r^{2}} \partial_{θ} η \partial_{θ} ξ + (λ + 2 μ) \partial_{r} η \partial_{r} ξ & \frac{λ}{r} \partial_{r} η \partial_{θ} ξ + \frac{μ}{r} \partial_{r} ξ \partial_{θ} η \\ \frac{λ}{r} \partial_{r} ξ \partial_{θ} η + \frac{μ}{r} \partial_{r} η \partial_{θ} ξ & μ \partial_{r} η \partial_{r} ξ + \frac{(λ + 2 μ)}{r^{2}} \partial_{θ} ξ \partial_{θ} η \end{array}], \end{array}

(21)

\begin{array}{r} C = - [\begin{array}{c} \frac{λ}{r} \partial_{r} η & - \frac{μ}{r^{2}} \partial_{θ} η \\ \frac{(λ + 2 μ)}{r^{2}} \partial_{θ} η & - \frac{μ}{r} \partial_{r} η \end{array}] . \end{array}

(22)

In the traction image method, the traction values in the ghost points are antisymmetry about the free surface and they are expressed as follows,

\begin{array}{r} {{\vec{T}}_{{\hat{e}}_{r}} |}_{k_{η} = n + k} = - {{\vec{T}}_{{\hat{e}}_{r}} |}_{k_{η} = n - k}, \end{array}

(23)

\begin{array}{r} {{\vec{T}}_{{\hat{e}}_{θ}} |}_{k_{η} = n + k} = - {{\vec{T}}_{{\hat{e}}_{θ}} |}_{k_{η} = n - k}, \end{array}

(24)

where $k_{η}$ refers to the gridpoint in the $η$ direction in traction image method, ‘n’ refers to the gridpoint in the surface, and ‘k’ refers to the distance of the gridpoint in the traction image method from the surface.

Therefore, all the $η$ direction derivatives of the stress components in eqs (10)–(14) should be expressed through the traction force and its image. By calculating the $η$ direction derivatives for eqs (17) and (18), obtain

\begin{array}{r} \frac{1}{r^{2}} \partial_{θ} η σ_{r θ} \partial_{η} r + | T | \partial_{η} [\frac{\partial_{r} η σ_{r θ} + \frac{1}{r} \partial_{θ} η σ_{θ θ}}{| T |}] - \frac{\partial_{θ} η \partial_{θ} η \partial_{η} r}{r^{3} {| T |}^{2}} (\partial_{r} η σ_{r r} + \frac{1}{r} \partial_{θ} η σ_{θ r}), \end{array}

(25)

\begin{array}{r} \frac{1}{r^{2}} \partial_{θ} η σ_{θ θ} \partial_{η} r + | T | \partial_{η} [\frac{\partial_{r} η σ_{r θ} + \frac{1}{r} \partial_{θ} η σ_{θ θ}}{| T |}] - \frac{\partial_{θ} η \partial_{θ} η \partial_{η} r}{r^{3} {| T |}^{2}} (\partial_{r} η σ_{r θ} + \frac{1}{r} \partial_{θ} η σ_{θ θ}), \end{array}

(26)

where

\begin{array}{r} | T | = | \partial_{r} η {\hat{e}}_{r} + \frac{1}{r} \partial_{θ} η {\hat{e}}_{θ} | . \end{array}

(27)

So far, all the $η$ direction derivatives of velocity and stress components at the free surface in eqs (10)–(14) are expressed by the derivatives of stress component, traction component and traction component of $η$ direction. Free-surface boundary conditions involved in the free surface in eqs (10)–(14) are completed.

2.3 Overset-grid and data exchange

The general concept of overset-grid is to decompose computational domains with complex geometric boundaries into subregions with relatively simple geometric boundaries. After processing the free-surface conditions by the traction image method, the remaining question is how to generate a grid at the centre of the circle. A singularity will occur if we generate an ordinary polar coordinate grid at the centre of the circle and the spacing of the grid near the centre is really small which means the time interval is limited. In this work, an overset-grid is utilized to generate grid at the centre and interpolation is used to exchange the data.

In our overset-grid implementation, first the curvilinear grid is used to fit the surface topography of the circular study area. Then a Cartesian coordinate grid is generated to discrete the central part of the circular. The Cartesian coordinate grid and the curvilinear-grid arbitrarily overlap each other and the points of the two grids do not need to match within the overlapping region. The Cartesian coordinate grid not only avoids the problem that the polar coordinate system cannot be generated at the origin, but also greatly improves the calculation efficiency without calculating the deformation coefficient of the Cartesian coordinate grid.

In the process of seismic wave simulation using overset-grid, the key problem lies in the data exchange from one grid block to another. A variety of interpolation methods have been used for data exchange of overset-grid, including second-order linear interpolation (Delfs 2001; Khokhlov et al. 2019), higher order explicit interpolation (Sherer & Scott 2005) and Hermite interpolation (Guénanff 2004). Sherer & Scott (2001) make a full comparison of the different interpolation methods, recommending a generalized higher order explicit Lagrangian interpolation. Lagrangian interpolation and inverse distance weight (IDW) of different orders are used for testing.

The Lagrangian interpolation formula is as

\begin{array}{r} W_{target} (x_{0}, z_{0}) = \sum_{i = 1}^{N} \sum_{j = 1}^{N} L_{i j} W_{serve} (x_{i}, z_{j}), \end{array}

(28)

where $W_{target}$ is the value of point $(x_{0}, z_{0})$ and this value is provided by point $(x_{i}, z_{j})$ of $W_{serve}$ ⁠. N is the width of the interpolated template, and $L_{i j}$ is the Lagrangian polynomial. For the centred interpolation of templates with length N, $L_{i j}$ can be expressed as:

\begin{array}{r} L_{i j} = \prod_{l = 1, l i}^{N} \frac{x_{0} - x_{l}}{x_{i} - x_{l}} \prod_{k = 1, k j}^{N} \frac{z_{0} - z_{k}}{z_{j} - z_{k}} . \end{array}

(29)

The IDW interpolation formula is as

\begin{array}{r} Z (x_{0}, z_{0}) = \sum_{1}^{n} w_{i} * Z_{serve} (x_{i}, z_{i}), \end{array}

(30)

where $Z (x_{0}, z_{0})$ is the value of point $(x_{0}, z_{0})$ which is provided by point $(x_{i}, z_{i})$ of $Z_{serve}$ ⁠. $w_{i}$ is the weight of the value of the point $(x_{i}, z_{i})$ ⁠, which can be expressed as:

\begin{array}{r} w_{i} = \frac{1 / d_{i}}{\sum_{1}^{n} 1 / d_{i}}, \end{array}

(31)

where d_i is the distance of point $(x_{i}, z_{i})$ and point $(x_{0}, z_{0})$ ⁠. We will test the robustness and accuracy of the different interpolation methods in the following part.

The data exchange diagram between the Cartesian grid and the curvilinear grid is shown in Fig. 2. Taking the sixth-order interpolation as an example, the red point represents the receiving point in curvilinear and Cartesian grids that needs to be interpolated. The 6 × 6 interpolation points can be easily found (the points in light blue area). The same method can be used to interpolate data from the curvilinear to Cartesian grid. Considering each grid of Cartesian and the curvilinear grids requires three layers of ghost points on the surface when we use DRP (Dispersion Relation Preservation)/opt MacCormack differential format (Hixon 1997), these three-layer grids are also the points that need to obtain the data through interpolation for both Cartesian and curvilinear grids. It should also be noted that in order to ensure the stability of spatial difference, there must be at least three physical points in each grid overlapping each other to provide interpolation information in case of the difference format and interpolation format is three.

Figure 2.

Diagram of 2-D interpolation. (a) Diagram of data transfers from Cartesian grid to curvilinear- rid. The red point is the interpolation receiving point in the polar coordinate system that receives the interpolation information from black diamonds in the Cartesian grid. (b) Diagram of data transfers from curvilinear- rid to Cartesian grid. The red point belongs to the ghost point layer in the Cartesian grid. Black diamonds that belong to the curvilinear- rid provide the information used to interpolate.

Open in new tab Download slide

2.4 Parallel processing of the overset-grid program

Overset-grid algorithm allows overlapping between grid blocks when research areas involve complex geometry. It reduces unnecessary grids and avoids complex topological partitions, which compensates for the shortcomings of structural grids to some extent. However, compared with structured grids, the parallel processing of overset-grid is more complicated because the position relationship between the receiving point and the sending point in the interpolation process is hard to determine for their irregular position relationship. Overset-grid programs are always used to solve the simulation problem of regional-global scale, so the large computation amount requires us to process the program in parallel. FDM requires the data of each gridpoint at each time step and the operation of each gridpoint only needs the data of several surrounding points, which means it is suitable to use the parallel algorithm of regional decomposition type. So we used the Message Passing Interface method (MPI) to parallel the programs.

Unlike regular grid MPI, which only requires the data exchange at the edge of each MPI block with half-length layers of difference template, MPI of overset-grid has two types of data communication. One is similar to the regular grid, transferring the data in the half-length layers of difference template between the MPI blocks according to the topological relationship; the other communication is the transmission of interpolation data between the two subgrids. The second communication does not just involve the transmission of information between a fixed number of MPI blocks. For example, in the MPI model of four blocks, for a block, the first communication type only involves its own adjacent MPI blocks, while the second communication type may involve all four MPI blocks because of the geometric relationship in interpolation process. Therefore, how to produce the parallel of interpolation process is the core and difficulty.

The first step in parallel the interpolation process is establishment of MPI model. There are two common MPI models in overset-grid: each block handles only one type of subgrid; the other is to divide all the subgrids evenly to all blocks, where each block handles all types of subgrids. The first type of MPI model is more appropriate in problems with a very complex geometry or involves aerodynamics that a grid moves over time. Since our program deals with a spherical overset-grid, which is relatively simple and does not change with time, we chose the second type of MPI model (Fig. 3). This type also automatically meets the load balancing requirements.

Figure 3.

MPI communication model. The grey area indicates the area that needs to be interpolated and the serial number indicates the process number. The grid of the same process number operates on the same block.

Open in new tab Download slide

The next core problem is to find the send–recv relationship for the points that need to be interpolated in each process. It means the block number of the source and the location of the sending point in its own block. We use the adjacent lattice search method (Dong & Wang 2011) to determine the spatial location relationship. The process of adjacent lattice search method is similar to the template traversal method, which finds the set of adjacent points by comparing the distance between the goal point in the grid with its neighbour points. The whole MPI workflow is shown as follows:

gather grid of each block and obtain the whole grid information;
find the block number of the points that each block needs to receive and send and the coordinates of these points in the block;
at each time step, pack the information sent by each block and unpack the information each block receives.

2.5 Other numerical calculation and processing

The spatial difference scheme used in this algorithm is DRP/opt MacCormack difference scheme (Hixon 1997), which optimizes the dispersion and dissipation error of the scheme and greatly improves the grid resolution. The Runge–Kutta integral of order 4–6 is used in this work. Source terms can be added at a single point or spatially smoothed at several points following the pseudo-δ distribution (Herrmann 1979). In our experience, adding source terms at a single point always causes strong artificial reflections while spatially smoothing the source terms significantly reduces these reflections. Therefore, a smooth spatial distribution is used to add the source terms. If the source space is located in the curvilinear-grid area, it will be added to the curvilinear grid; otherwise, it will be added to the Cartesian grid. At the same time, it cannot be ignored that the source should not be located in overlapping areas. For circumstance where the source is located in the overlapping region, we can generate more curvilinear-grid layers and further move the overlapping region down to avoid the source locating in it.

For the three-layer gridpoints near the complex topography surface of the curvilinear grid, the difference template crosses the free surface, meaning that the $η$ direction derivative of velocity/stress component needs additional treatment. The processing method of $η$ direction derivative of stress component and velocity component follows the compact difference scheme (Hixon & Turkel 2000). It is used to calculate the $η$ direction derivative value of the velocity and stress component for the gridpoints of nearly complex topography surface with the fourth-order accuracy.

In addition, since the numerical simulation process of the curvilinear grid includes calculation of spatial derivatives whose number is about twice that of the Cartesian grid, the computation time required is also about twice of Cartesian grid. At the same time, with the increase of the number of curvilinear-grid layers, the grid distance in $θ$ direction will constantly decrease inevitably which reduces the computing efficiency. Therefore, we prefer to increase the area of the Cartesian grid area and reduce the area simulated by the curvilinear grid. In addition, compared with the non-overset-grid scheme of the same region model, the increased computation in overset grid is not obvious because only interpolation in three layers of each subgrid is added.

3 NUMERICAL VERIFICATION

In this section, we check the validity and accuracy of our method by comparing with solutions of the direct solution method (DSM) (Takeuchi et al. 1996) and DG (COMSOL Multiphysics 6.0). DSM gives exact waveforms for globally symmetric media without topography while DG solutions are reference solutions for the global model with a topography. To further quantitatively assess the agreement between two methods both in time and frequency domains, we calculate the time–frequency (TF) misfit (Kristekova et al. 2006) of the overset-grid FDM solutions and the reference solutions. The benefit of TF evaluation is that it can separate the envelope misfit (EM) and the phase misfit (PM), and provide the misfit information both in time and frequency domains simultaneously.

Some parameters of the two models are the same. The P-wave speed is 3 km s⁻¹ and the S-wave speed is 2 km s⁻¹, and the density is 1500 kg m⁻³. The source time function is Ricker wavelet with a central frequency of 0.96 Hz and time delay of 1 s. The curvilinear-grid spacing is no more than 0.05 km in the $θ$ -direction and 0.05 km in the r-direction, and the grid spacing of Cartesian grid is 0.05 km. The number of the curvilinear-gridpoints of the overset-grid FDM is 1800 × 300. The number of the Cartesian grid in the centre is 411 × 411.

3.1 Horizontal surface homogeneous medium model

To validate the accuracy of the overset-grid FDM, we first compare waveforms with the DSM in a spherical model of horizontal surface. The source and receivers are distributed on the same circle (Fig. 4). Receivers are evenly distributed on the surface and the source is 5 km depth under the surface. In the polar coordinate system, the source is a 2-D double-couple point source with moment tensor components $M_{θ θ}$ =1.0, $M_{r r}$ = −1.0, and it can be approximated to a 3-D point source in the spherical coordinate system by the transformation of Vidale et al. (1984),

\begin{array}{r} v_{3 - D} = \frac{1}{π \sqrt{D}} (\frac{1}{\sqrt{t}} * v_{2 - D}), \end{array}

(32)

where * is the convolution operation; $v_{2 - D}$ is the waveform of the 2-D point source and $v_{3 - D}$ is the conversion waveform of the point source. D is the distance between the source and the receiver. It is worth mentioning that generally the 3-D geometric diffusion effect can be correctly handled by 3-D and 2.5-D modelling, but not by 2-D modelling (Igel & Weber 1995; Toyokuni et al. 2005).

Figure 4.

Diagram of overset-grid horizontal surface model.

Open in new tab Download slide

The comparison of the converted overset-grid and DSM seismograms is shown in Fig. 5. The blue curve is the simulation result of the polar coordinate overset-grid FDM and the red curve is the result of DSM. The light blue curve is the difference of the two solutions. It is clear to see that the simulation results of the two algorithms are almost the same. Their PM and EM values are shown in Table 1. It can be seen that PM and EM values are lower than 0.4 and 0.7 per cent, respectively, meaning the excellent level of goodness of fit.

Figure 5.

Waveform comparison of the flat surface model.

Open in new tab Download slide

Table 1.

Open in new tab

EM and PM of $v_{θ}$ and $v_{r}$ components for all receivers, in the flat surface model with a double-couple point source, in percentage.

No.	EM (⁠ $v_{θ}$ ⁠)	EM (⁠ $v_{r}$ ⁠)	PM (⁠ $v_{θ}$ ⁠)	PM (⁠ $v_{r}$ ⁠)
Receiver 1	0.47	0.65	0.38	0.26
Receiver 2	0.49	0.43	0.38	0.38
Receiver 3	0.63	0.54	0.19	0.39
Receiver 4	0.38	0.42	0.17	0.26
Receiver 5	0.32	0.52	0.02	0.29
Receiver 6	0.57	0.47	0.38	0.36
Receiver 7	0.58	0.66	0.39	0.36
Receiver 8	0.61	0.41	0.38	0.28

No.	EM (⁠ $v_{θ}$ ⁠)	EM (⁠ $v_{r}$ ⁠)	PM (⁠ $v_{θ}$ ⁠)	PM (⁠ $v_{r}$ ⁠)
Receiver 1	0.47	0.65	0.38	0.26
Receiver 2	0.49	0.43	0.38	0.38
Receiver 3	0.63	0.54	0.19	0.39
Receiver 4	0.38	0.42	0.17	0.26
Receiver 5	0.32	0.52	0.02	0.29
Receiver 6	0.57	0.47	0.38	0.36
Receiver 7	0.58	0.66	0.39	0.36
Receiver 8	0.61	0.41	0.38	0.28

Table 1.

Open in new tab

EM and PM of $v_{θ}$ and $v_{r}$ components for all receivers, in the flat surface model with a double-couple point source, in percentage.

No.	EM (⁠ $v_{θ}$ ⁠)	EM (⁠ $v_{r}$ ⁠)	PM (⁠ $v_{θ}$ ⁠)	PM (⁠ $v_{r}$ ⁠)
Receiver 1	0.47	0.65	0.38	0.26
Receiver 2	0.49	0.43	0.38	0.38
Receiver 3	0.63	0.54	0.19	0.39
Receiver 4	0.38	0.42	0.17	0.26
Receiver 5	0.32	0.52	0.02	0.29
Receiver 6	0.57	0.47	0.38	0.36
Receiver 7	0.58	0.66	0.39	0.36
Receiver 8	0.61	0.41	0.38	0.28

No.	EM (⁠ $v_{θ}$ ⁠)	EM (⁠ $v_{r}$ ⁠)	PM (⁠ $v_{θ}$ ⁠)	PM (⁠ $v_{r}$ ⁠)
Receiver 1	0.47	0.65	0.38	0.26
Receiver 2	0.49	0.43	0.38	0.38
Receiver 3	0.63	0.54	0.19	0.39
Receiver 4	0.38	0.42	0.17	0.26
Receiver 5	0.32	0.52	0.02	0.29
Receiver 6	0.57	0.47	0.38	0.36
Receiver 7	0.58	0.66	0.39	0.36
Receiver 8	0.61	0.41	0.38	0.28

3.2 Model of overset-grid interpolation with complex topography

In this model, the computing area is a global area with a topography and the free surface is represented using the triangular function,

\begin{array}{r} y = 1000 * \sin (5 \cdot θ), \end{array}

(33)

where $θ$ is the radian of the surface gridpoint and y is the topography of the free surface. Two lines of receivers are set in this model. The first line is on the surface and the second line is along the polar radius in the computing area with interval of 500 m. The positions of the source and receivers are shown in Fig. 6. We use a single force point source; The time step is 0.01 s and the duration is 10 s. For overset-grid, the source is in a square grid.

Figure 6.

Diagram of overset-grid complex topography model.

Open in new tab Download slide

Synthetic waveforms of the two lines are shown in Fig. 7. EM and PM of receivers in two lines are shown in Tables 2 and 3, respectively. Their PM and EM values are lower than 0.3 and 0.8 per cent, respectively, which means an excellent level of agreement.

Waveform of the overset-grid comparing with DG algorithm. (a)Waveform of receivers at grid surface. (b) Waveform of receivers in a radius direction.

Open in new tab Download slide

Figure 7.

Waveform of the overset-grid comparing with DG algorithm. (a)Waveform of receivers at grid surface. (b) Waveform of receivers in a radius direction.

Table 2.

Open in new tab

EM and PM of $v_{x}$ and $v_{z}$ components for receivers at the surface, in the complex topography model with a single force point source, in percentage.

No.	EM (⁠ $v_{x}$ ⁠)	EM (⁠ $v_{z}$ ⁠)	PM (⁠ $v_{x}$ ⁠)	PM (⁠ $v_{z}$ ⁠)
Receiver 1	0.76	0.59	0.16	0.13
Receiver 2	0.63	0.58	0.13	0.15
Receiver 3	0.56	0.58	0.18	0.16
Receiver 4	0.58	0.61	0.14	0.27
Receiver 5	0.58	0.42	0.15	0.27
Receiver 6	0.48	0.49	0.23	0.19
Receiver 7	0.49	0.51	0.21	0.17
Receiver 8	0.49	0.50	0.17	0.21

No.	EM (⁠ $v_{x}$ ⁠)	EM (⁠ $v_{z}$ ⁠)	PM (⁠ $v_{x}$ ⁠)	PM (⁠ $v_{z}$ ⁠)
Receiver 1	0.76	0.59	0.16	0.13
Receiver 2	0.63	0.58	0.13	0.15
Receiver 3	0.56	0.58	0.18	0.16
Receiver 4	0.58	0.61	0.14	0.27
Receiver 5	0.58	0.42	0.15	0.27
Receiver 6	0.48	0.49	0.23	0.19
Receiver 7	0.49	0.51	0.21	0.17
Receiver 8	0.49	0.50	0.17	0.21

Table 2.

Open in new tab

EM and PM of $v_{x}$ and $v_{z}$ components for receivers at the surface, in the complex topography model with a single force point source, in percentage.

No.	EM (⁠ $v_{x}$ ⁠)	EM (⁠ $v_{z}$ ⁠)	PM (⁠ $v_{x}$ ⁠)	PM (⁠ $v_{z}$ ⁠)
Receiver 1	0.76	0.59	0.16	0.13
Receiver 2	0.63	0.58	0.13	0.15
Receiver 3	0.56	0.58	0.18	0.16
Receiver 4	0.58	0.61	0.14	0.27
Receiver 5	0.58	0.42	0.15	0.27
Receiver 6	0.48	0.49	0.23	0.19
Receiver 7	0.49	0.51	0.21	0.17
Receiver 8	0.49	0.50	0.17	0.21

No.	EM (⁠ $v_{x}$ ⁠)	EM (⁠ $v_{z}$ ⁠)	PM (⁠ $v_{x}$ ⁠)	PM (⁠ $v_{z}$ ⁠)
Receiver 1	0.76	0.59	0.16	0.13
Receiver 2	0.63	0.58	0.13	0.15
Receiver 3	0.56	0.58	0.18	0.16
Receiver 4	0.58	0.61	0.14	0.27
Receiver 5	0.58	0.42	0.15	0.27
Receiver 6	0.48	0.49	0.23	0.19
Receiver 7	0.49	0.51	0.21	0.17
Receiver 8	0.49	0.50	0.17	0.21

Table 3.

Open in new tab

EM and PM of $v_{x}$ and $v_{z}$ components for receivers in a radius, in the complex topography model with a single force point source, in percentage.

No.	EM (⁠ $v_{x}$ ⁠)	EM (⁠ $v_{z}$ ⁠)	PM (⁠ $v_{x}$ ⁠)	PM (⁠ $v_{z}$ ⁠)
Receiver 1	0.80	0.47	0.04	0.04
Receiver 2	0.68	0.77	0.00	0.03
Receiver 3	0.63	0.73	0.04	0.08
Receiver 4	0.61	0.59	0.18	0.11
Receiver 5	0.60	0.47	0.05	0.07
Receiver 6	0.58	0.57	0.13	0.05
Receiver 7	0.57	0.61	0.12	0.03
Receiver 8	0.61	0.59	0.04	0.12

No.	EM (⁠ $v_{x}$ ⁠)	EM (⁠ $v_{z}$ ⁠)	PM (⁠ $v_{x}$ ⁠)	PM (⁠ $v_{z}$ ⁠)
Receiver 1	0.80	0.47	0.04	0.04
Receiver 2	0.68	0.77	0.00	0.03
Receiver 3	0.63	0.73	0.04	0.08
Receiver 4	0.61	0.59	0.18	0.11
Receiver 5	0.60	0.47	0.05	0.07
Receiver 6	0.58	0.57	0.13	0.05
Receiver 7	0.57	0.61	0.12	0.03
Receiver 8	0.61	0.59	0.04	0.12

Table 3.

Open in new tab

EM and PM of $v_{x}$ and $v_{z}$ components for receivers in a radius, in the complex topography model with a single force point source, in percentage.

No.	EM (⁠ $v_{x}$ ⁠)	EM (⁠ $v_{z}$ ⁠)	PM (⁠ $v_{x}$ ⁠)	PM (⁠ $v_{z}$ ⁠)
Receiver 1	0.80	0.47	0.04	0.04
Receiver 2	0.68	0.77	0.00	0.03
Receiver 3	0.63	0.73	0.04	0.08
Receiver 4	0.61	0.59	0.18	0.11
Receiver 5	0.60	0.47	0.05	0.07
Receiver 6	0.58	0.57	0.13	0.05
Receiver 7	0.57	0.61	0.12	0.03
Receiver 8	0.61	0.59	0.04	0.12

No.	EM (⁠ $v_{x}$ ⁠)	EM (⁠ $v_{z}$ ⁠)	PM (⁠ $v_{x}$ ⁠)	PM (⁠ $v_{z}$ ⁠)
Receiver 1	0.80	0.47	0.04	0.04
Receiver 2	0.68	0.77	0.00	0.03
Receiver 3	0.63	0.73	0.04	0.08
Receiver 4	0.61	0.59	0.18	0.11
Receiver 5	0.60	0.47	0.05	0.07
Receiver 6	0.58	0.57	0.13	0.05
Receiver 7	0.57	0.61	0.12	0.03
Receiver 8	0.61	0.59	0.04	0.12

Fig. 8 also shows that the wave field snapshots (at 6, 8, 10 and 12 s) can well propagate through the junction of the overset-grid and are reflected from the free surface correctly where no spurious reflections are observed. In addition, the accuracy of Lagrange interpolation and IDW interpolation of different orders are also tested. For the surface receivers, average EM of the second- and sixth-order Lagrange interpolations are 0.57 and 0.49 per cent, respectively, and average EM of the second- and sixth-order IDW are 0.51 and 0.54 per cent, respectively. It is clear to see that Lagrange and IDW interpolations of second- and sixth-order have almost the same accuracy. In the following studies, the sixth-order Lagrange interpolation is used in the global-scale model.

Figure 8.

Wavefield of the Vz component of the overset-grid model at 6, 8, 10 and 12 s.

Open in new tab Download slide

The numerical consistency and stability for the overset-grid FDM is also tested. Receivers are uniformly in all directions on the flat surface of the curvilinear grid and the computing has 10 000 steps. Fig. 9 shows the waveform of Vz component. It can be seen that the amplitude phase of the waveform in each direction is highly consistent, and the algorithm is still stable after multiple interpolations and free-surface reflections, which proves the stability of this algorithm.

Figure 9.

Waveform in all directions on the surface of the overset-grid. The time step is 0.01 s and duration is 100 s.

Open in new tab Download slide

In general, the accuracy of the proposed finite-difference algorithm is validated by the flat surface global model and global model with complex topographies.

4 NUMERICAL SIMULATIONS OF THE FULL-MARS MODEL

Mars is an Earth-like planet in the solar system and it has the most similarity to Earth. The study of the internal structure and material composition of Mars is of great significance to understand the evolution of Mars, especially the generation and extinction mechanism of the magnetic field as well as the change of the habitable environment of Mars, which is an important topic of space exploration.

The study of seismology on Mars began in 1976, but Viking 1 and Viking 2 did not collect reliable seismic data (Anderson et al. 1977). Some studies using seismological methods to research Mars, such as Martian seismic activity, Martian quake localization, Martian shallow structure and crust-mantle structure (Zheng et al. 2014; Khan et al. 2016; Bissig et al. 2018; Knapmeyer-Endrun et al. 2018). On 2018 November 26, the US Insight (Insight) Mars probe landed near the Mars equator and installed a Mars seismometer. As of 2019 September 30, the seismometer has recorded 174 earthquakes, including more than 20 earthquakes of a 3–4 mag (Stahler et al. 2021). Seismic wave simulation of Mars can provide the reference of the acquired seismic data and the structure of Mars.

Sohl & Spohn (1997) proposes two models of the internal structure of Mars. He obtains the radial distribution of Mars mass, hydrostatic pressure, gravity, temperature and heat flow density under the conditions of hydrostatic balance and closed heat transfer. He also obtains the hierarchical density of Mars model and the structure model of Mars seismic velocity. The two models meet the constraint of the polar moment of inertia reaching the maximum possible value and the geochemical constraint of the SNC meteorite from Mars (Fig. 10a). Stahler et al. (2021) use three methods of seismology, geophysical and geodynamics modelling, picking up a variety of body waves in the seismic data and obtaining the vertical wave and transverse wave velocity with depth profile of Mars mantle and core. Combining the results of the three modelling methods and referring to the Mars gravity and satellite orbit measurements, Stahler gives the size of 1830 ± 40 km and the density of 5.7–6.3 g cm⁻³ of the Mars core. Considering the conclusions of the model proposed by Sohl and Stahler, Sohl's model B is choosen for simulation.

Figure 10.

(a) Sohl (1997). The P- and S-wave velocities, and density distribution of the full Mars model. The solid line represents model A, and the dotted line represents model B. (b) The topographic map of Mars. The line indicates the trace of the polar grid topography on the Martian surface.

Open in new tab Download slide

Considering that most of the Mars probes have historically landed between south latitude 30° and north latitude 30° and that the overall topography of Mars presents a strong dichotomization where the flat low land elevation fluctuates little in the Northern Hemisphere and the plateau in the Southern Hemisphere is much more pronounced. In addition, the overall topography of Mars is in the range of −8 to 22 km. The proposed overset-grid FDM supports the simulation of arbitrary complex topography, so the topographic data of the curvilinear grid in this work on the Martian surface is set between south latitude 30° and north latitude 30°, as shown in Fig. 10(b). Since the source mechanism of Mars earthquakes is not currently constrained by observational data, the double-couple source of $M_{θ θ}$ =1, $M_{r r}$ = −1 and $M_{θ r}$ = 0, which is equivalent to a dip fault with 45 inclination, is utilized in this simulation. The source is 30 km below the free surface. The source time function is a Ricker wavelet with a central frequency of 0.2 Hz and a time delay of 50 s. The radius of Mars is 3390 km and it is discretized 600 gridpoints in the r-direction by a grid interval of 2 km. The gridpoints in the $θ$ -direction is 10 700 and the grid interval is less than 2 km. The number of Cartesian gridpoints in the centre is 2207 × 2207 and the grid interval is 2 km. The gridpoints per wavelength in this simulation is about 8. The time step is 0.1 s and duration is 300 s. The inelastic attenuation factor Q in full-Mars model is set to 1 000 000 in the liquid core and 600 in other areas. In addition, the source is a 2-D line source in the simulation, and it can be approximated to a 3-D point source in the spherical coordinate system by eq. (32).

Fig. 11 shows synthetic seismic seismograms on the Martian surface. It is clear to see that due to the influence of the Mars core, the speed of P wave propagated to the nuclear mantle boundary will be significantly reduced and the direction of propagation is also changed. Therefore, there is a direct P-wave shadow area and the size of the shadow area is related to the radius of the Mars core. Direct P wave, direct S wave and multiple reflected wave PP wave of P wave can be observed on both components. Seismic phases of PcP, ScS and phase of multiple reflected wave SSS are also more prominent in the tangential component. Multiple reflected phase is generated between direct S and SS waves because of strong low-speed Mars crust. In addition, it can be seen that the surface waves of two components are relatively clear and obvious. With the propagation of the surface waves, the wave train is getting longer and longer and the dispersion phenomenon is more and more obvious. At the same time, by comparing waveforms of model with and without topography, it can be seen that topography affects both arrival time and amplitude of seismic wave. The amplitude of body waves is mainly influenced, while the phase velocity and frequency dispersion of surface waves are both affected.

Figure 11.

Theoretical seismograms of the radial and vertical components of Martian surface. Major phases are marked.

Open in new tab Download slide

Fig. 12 shows snapshots of wave field propagating in the full Mars model. It can be seen that at 200 s, the direct P and SV waves spread downward in the upper Martian mantle. The front of the wave is very clear and the P wave is followed by the pP and sP waves. The complex wavefront is generated in the reflection of the surface and the conversion of the direct P and SV waves due to the shallow and low-speed layer. In the Martian crust, slow surface waves are already beginning to occur. At 400 s, the direct SV wave reaches the boundary of mantle and core, while P and pP waves generated from the reflection of core–mantle boundary can be clearly seen. At 600 s, the reflected and converted P waves of the direct SV waves from the core–mantle boundary propagate upward. Due to the existence of the Martian core, the curved diff P wave can be seen near the core–mantle boundary and shadow area of direct P wave is generated. At 800 s, P wave has spread to the other side of the Martian core and the converted SV wave is generated. There are also various ‘Y’ -shaped multiple reflected waves in the wavefield. In the following time, the wavefield is strongly disturbed by waves from upper Martian mantle and various transformations occur at the core–mantle boundary. From these snapshots, the length of the wave trains produced by the conversion of interference in upper Martian crust with low velocity increases with the propagation distance.

Figure 12.

Continuous wavefield snapshots of P–SV waves in the full Mars model. The inner circle represents the Martian core.

Open in new tab Download slide

Combining theoretical seismograms with wavefield snapshots, it is clear that layers in upper mantle of low velocity are important to the generation of surface wave while conversion and reflection occurs as wave pass through the Martian core. Many areas of the surface of Mars are in the shadow zone. This provides an important reference for the future layout of seismometers. In addition, it can be seen that topography affects both arrival time and amplitude of seismic waves by comparing waveforms of model with and without topography, so the influence of the terrain should not be ignored in study of region areas.

5 SUMMARY

In this study, we develop a finite-difference numerical simulation method with arbitrary topography using an overset-grid program consisting of polar coordinate curvilinear and Cartesian grids, which can be applied to the regional-global simulation problem. The polar coordinate curvilinear grid is used to deal with complex topography and the Cartesian in the centre of the circle is used to solve the problem of the singularity at the centre of the circle. We exchange data between different grids using explicit Lagrange interpolation of different orders. Additionally, adjacent lattice search method is utilized to determine the geometric positional relationship between the grids when parallel the overset-grid program.

In addition, two examples and reference models are designed to compare the results, showing that our overset-grid algorithm is very accurate in wavefield propagation with free-surface conditions and grid interpolation. It also has a good stability in the long-time simulations. Finally, we simulate the wave field of full Mars model; the seismograms and wavefield snapshots show a variety of different seismic phases, indicating that our overset-grid algorithm of polar coordinate system deals with the simulation of region-global model well.

ACKNOWLEDGMENTS

This study is supported by National Key R&D Program of China (2022YFC3003504), Guangdong Provincial Key Laboratory of Geophysical High-resolution Imaging Technology (2022B1212010002), Shenzhen Science and Technology Program (grant no. KQTD20170810111725321) and Center for Computational Science and Engineering at Southern University of Science and Technology.

DATA AVAILABILITY

The data that support the findings of this study and the program are all available from the URL (https://github.com/qiuhk/Seismic-numerical-simulation).

REFERENCES

Anderson

D.L.

et al.

1977

Seismology on Mars

J. geophys. Res. (1896-1977)

4524

–

4546

Month:	Total Views:
September 2024	36
October 2024	32
November 2024	25
December 2024	13
January 2025	10
February 2025	14
March 2025	22
April 2025	56
May 2025	10

Article Contents

An overset-grid finite-difference algorithm for seismic wavefield propagations modelling in the polar coordinate system with a complex free-surface topography

SUMMARY

1 INTRODUCTION

2 METHOD

2.1 The first-order velocity-stress equation in the Polar coordinates

2.2 Implement the free-surface boundary condition by the Traction image method in the polar coordinates

2.3 Overset-grid and data exchange

2.4 Parallel processing of the overset-grid program

2.5 Other numerical calculation and processing

3 NUMERICAL VERIFICATION

3.1 Horizontal surface homogeneous medium model

3.2 Model of overset-grid interpolation with complex topography

4 NUMERICAL SIMULATIONS OF THE FULL-MARS MODEL

5 SUMMARY

ACKNOWLEDGMENTS

DATA AVAILABILITY

REFERENCES

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

This Feature Is Available To Subscribers Only