Subdomain Chebyshev Spectral Method for 2 D and 3 D Numerical Differentiations in a Curved Coordinate System

A new numerical approach, called the “subdomain Chebyshev spectral method” is presented for calculation of the spatial derivatives in a curved coordinate system, which may be employed for numerical solutions of partial differential equations defined in a 2D or 3D geological model. The new approach refers to a “strong version” against the “weak version” of the subspace spectral method based on the variational principle or Galerkin’s weighting scheme. We incorporate local nonlinear transformations and global spline interpolations in a curved coordinate system and make the discrete grid exactly matches geometry of the model so that it is achieved to convert the global domain into subdomains and apply Chebyshev points to locally sampling physical quantities and globally computing the spatial derivatives. This new approach not only remains exponential convergence of the standard spectral method in subdomains, but also yields a sparse assembled matrix when applied for the global domain simulations. We conducted 2D and 3D synthetic experiments and compared accuracies of the numerical differentiations with traditional finite difference approaches. The results show that as the points of differentiation vector are larger than five, the subdomain Chebyshev spectral method significantly improve the accuracies of the finite difference approaches.


Introduction
With development of computer technology, computer modeling that imitates the underground earth interior becomes now a powerful tool to geological, environmental and geophysical problems, e.g.numerical simulations of tectonic dynamics [1] [2], underground water flow or geothermal system [3]- [5], seismic wave propagation [6]- [9] and geo-electromagnetic fields [10] [11].Mathematically, computer modeling seeks approximate solutions of partial differential equations (PDE) that govern geological evolution or geophysical fields in the Earth, subject to a Dirichlet or Neumann boundary conditions.The simplest method of computer modeling is finitedifference approach to the partial derivatives of PDE, e.g.various high-order finite-difference formulae [12] and staggered grid schemes [13].Due to simple principle and easy implementation, they become now widely used for many geological and geophysical problems [14] [15].However, these methods often employ regular grids, and consequently suffer from complex geometries of geological models, e.g.various free-surface topography, such as hill, trench and folded rocks, which are often encountered in practice.In these cases, the finite-difference approach has to refine the free-surface topography and subsurface interfaces by stepwise curves or stepwise blocks [16], which cannot exactly match the geometries of geological models.To overcome the disadvantage, many researchers prefer the so-called "weak" solutions of PDE, i.e. finite-element method [17] or subspace spectral element method [18] [19], which is based on the variational principle or Galerkin's weighting scheme.Alternatively, one has to define a set of coordinate transformations and adapt the finite-difference approach to a curved coordinate system so as to match the free-surface topography and subsurface interfaces of geological models [20]- [22].
Spectral method, e.g.Chebyshev spectral method is superior to the traditional finite-difference approach in numerical differentiations because of exponential convergence [23].Due to its high accuracy, easy computer programming, and applicability to complex model geometries, the spectral method becomes one of the most popular accurate solvers of PDE [23]- [25].However, the most disadvantageous aspect of the spectral methods is expensive consumption of computer memory and CPU time, because it employs global collocation points in calculations of the spatial derivatives, resulting in a fully filling-in differentiation matrix [23].Therefore, it suffers the disadvantage in computer modeling for a large 3D geological or geophysical model [24] [25].
In this paper, we present a subdomain Chebyshev spectral method to calculate spatial derivatives in a curved coordinate system, which is similar to subspace spectral element method (weak form) [18] that divides the global domain into subdomains (elements) exactly matching the geometry of the model, and then Chebyshev points are applied to discretization of physical quantities and calculation of the spatial derivatives in each subdomain.The Chebyshev collocation points for the spatial differentiations only involve the local physical samples instead of the global values; therefore, the method saves huge computer memory for a 3D computer modeling.Such discretization results in a sparse assembled matrix like the "weak" version of subspace spectral element method [18] [19] and Gaussian quadrature grid approach [9].We conduct 2D and 3D synthetic experiments and compare the results with the finite-difference approaches.The results show that the new method is superior to the traditional finite-difference method in numerical differentiations and applicable to complex geological models having arbitrary free-surface topography and multiple subsurface interfaces.

Strong Solution
In general, computer modeling is to solve the following PDE [26]: subject to some Dirichlet or Neumann boundary conditions.Here ( ) ⋅ L is often a second-order linear differential operator depending on model vector m and partial derivatives ( ) represent location and physical field in a 2D or a 3D model domain Ω.The right-hand side vector ( ) s x is considered as a source of the physical field F defined in Ω.The vectors F and m are different physical quantities according to geological or geophysical problems, e.g. in tectonic dynamics, F may consist of pressure field p, deviatoric stress σ ′ and velocity field v, and m comprises the density ρ and viscosity η : [2]; in seismic wave modeling, F is the ground displacement vector u, or composed of speed v and stress tensor σ , and m includes density ρ and elastic moduli ijkl c :

{ }
, ijkl c ρ = m [8] [9]; in geo-electromagnetics, F becomes either electric field E or magnetic field H, and m involves electric permittivity ε , magnetic permeability µ and conductivity σ : [11].All these quantities are functions of the spatial coordinates ( ) . In order to match free-surface topography, e.g.hill, trench and ridge, one may employ the global coordinate transformations [20]- [22]: mapping the computational coordinates ( ) , , ξ ξ ξ into the physical coordinates ( ) , , x x x .In general, we assume F be a component of the vector F, and according to the chain law, we have the first and second derivatives: Here, summation convention of the repeated integer subscripts k and l has been applied.Accordingly, the discrete versions of Equations ( 3) and ( 4) can be written in the following forms where Here, ) and

D
are the first and second differentiation vectors and have the following generalised forms:

{ }
Substituting Equations ( 5) and ( 6) into (1), and then applying the PDE to every point { } , 1, 2, , , which have 3N unknowns of ( ) , , and gives a linear system.Solving the linear system subject to Dirichlet or Neumann boundary conditions, one obtains the so-called "strong" numerical solutions:

D
are determined by the neighbour points of p x , e.g. using a fourth-order finite-difference approach, the assembled matrix of Equation ( 10) becomes sparse and saves much computer memory.Then, the linear system can be efficiently solved by iterative or parallel solvers for a large 3D geological model [27] [28].Therefore, the crucial step of geological or geophysical modeling is to find accurate and efficient differentiation vectors
It should be noticed that the differentiation vectors

Coordinate Transformation
In order to exactly match free-surface topography and subsurface interfaces of the Earth, we divide the model domain Ω into non-overlapping subdomains, e.g. in 3D case, we have [ ] ( ) ( ) , , which maps the local coordinates ( x y z represent the size and central point of the subdomain ijk Ω , respectively, e.g. in the z-direction, we have and both are functions of x and y.Therefore, Equation (11) gives non-zero components of the Jacobian matrix: , , and non-zero components of the Hessian matrix: Equations ( 13) and (14) indicate that the Jacobian and Hessian matrices depend on the derivatives ( ) where ( ) ijk a αβ are coefficients defined in the subdomain ijk Ω and determined from the known samples ( ) z x y can be expressed by Equation (15), and these expressions are differentiable eve- rywhere in the model domain.Therefore, Equation (15) guarantees the partial derivatives ( ) for computing the Jacobian and Hessian matrices involved in Equations ( 13) and ( 14).

Subdomain Chebyshev Spectral Method
In the transformed subdomain to three coordinates ( ) ξ η ς , where N ξ , N η and N ς are numbers of the points in the three directions.Geo- metrically, these points are visualized as the projections on [ ] 1,1 − of equispaced points on the upper half of the unit circle.Figure 1 gives the two subdomains of the 1D case.One can easily extend it to 2D and 3D cases.
Substituting the Chebyshev points into equation (11), one obtains the global coordinates ( ) that demonstrate a one-to-one mapping between the two coordinate systems.All these Chebyshev points form a grid that discretizes the model parameters m and physical field vector F in Equation ( 1).This grid may be called "subdomain Chebyshev differentiation grid", which is similar to the "Gaussian quadrature grid" we used to sample the physical quantities in the elements [8] [9].Basing on the subdomain Chebyshev differentiation grid, we construct the Lagrange interpolation polynomials with the Chebyshev points in the subdomain, so that we have the first differentiation vectors and second differentiation vectors , , ′′ ⋅ = are their first and second derivatives.These differentiation vectors have exponential convergences of numerical differentiations in the subdomains and avoid the Runge phenomenonthe errors increase exponentially in the equispaced case [23].
Substitutions of Equations ( 13) and ( 17) for (7), we obtain general forms of the first differentiation vectors with respect to the global coordinates: , .
Similarly, substituting Equations ( 13), ( 14) and ( 19) into ( 8 .We omitted them here for saving space.However, the following simple method may be used to calculate these second differentiation vectors, e.g. according to Equation (5), we have following identity: and obtain , Here, we used ( ) ( )  21) is avoidance of computing the second derivatives of coordinate transformations, so that computation of the second differentiation vectors becomes much simpler than using Equation (6).
It should be noticed that Equations ( 19) and ( 21) are valid for the points inside the subdomains.At a boundary point of the subdomains, they may yield different values of the partial derivatives due to multiple subdomains' sharing.For the "strong" solution, the partial derivatives at the boundary points must be unique, for which we pick up half of Chebyshev points from the connecting subdomains (see points in the "boundary domain" in Figure 1) and construct the Lagrange interpolation polynomials given by Equation ( 17) and (18).Consequently, Equations ( 19) and ( 21) become well defined (uniqueness) for the partial derivatives at the boundary points.Apparently, the differentiation vectors at the boundary points differ from the ones at the points inside the subdomains due to the different distributions of Chebyshev points in the boundary domains and the ordinary subdomains (see Figure 1).This difference does not break down Equations ( 19) and (21), and on the contrary, it makes them applicable to any point in the whole domain.We call these differentiation schemes the "subdomain Chebyshev spectral method", because it applies the Chebyshev points to the subdomain rather than the globe domain presented in standard spectral method [23].Particularly, if we replace the Chebyshev points with equispaced points in the subdomains, Equations ( 19) and ( 21) become the central finite-difference approaches.It shows that the subdomain Chebyshev spectral method may regress to the finite-difference approaches.

Numerical Experiments
To demonstrate the capability of the new approach, we designed a 2D and a 3D geological model (see Figure 2).The former is a valley and the latter represents a ridge, both have two interfaces under the ground.With the given model geometries, we used subdomain Chebyshev points to discretise the model domain and obtained the subdomain Chebyshev differentiation grids.Figure 2 shows the grids of the two models, from which one can see 10 × 8 and 10 × 10 × 10 bent subdomains of the 2D and 3D geological models, respectively.The grids exactly match free-surface topography and subsurface interfaces.
In order to examine the accuracy of the subdomain Chebyshev spectral method, we employed the following testing functions as the physical quantities defined in the two models: , , sin 2π cos 2π sin 2π Here, L x , L y , and L z are lengths of the model domain in three coordinate directions.From the two functions, one can analytically calculate the first and second partial derivatives, which are used to check the accuracies of the subdomain Chebyshev spectral method.Figure 3 gives five partial derivatives ( ) , , , , .This is because of the linear transformation from ζ to z and non-linear relation- ship between ζ and x in Equation (11), which defines the changes of free-surface topography and subsurface interfaces in the horizontal direction.However, one can find that all absolutely relative errors are less than 0.3%.These results prove that the subdomain Chebyshev spectral method can provide accurate spatial derivatives in a 2D geological model.
In order to compare the subdomain Chebyshev spectral method with the traditional finite-difference approach, we repeated the above numerical experiments with finite-difference method and compared the average and maximum absolutely relative errors of the two methods.Meanwhile, five algorithms have been implemented three finite-difference approaches (FDM0, FDM1, FDM2) and two subdomain Chebyshev spectral methods (SCSM1, SCSM2).Here, the number "0" denotes the algorithm applying Lagrange interpolation functions to the free-surface topography ( ) x and subsurface interfaces ( ) k z x instead of spline interpolation functions.The integers "1" and "2" stand for the algorithms using Equation ( 8) and (21), respectively, for computing the second differentiation vectors, i.e. xx D , xz D and zz D .As mentioned above, the difference between the two algorithms is inclusion or exclusion of the high-order derivatives of the free-surface topography and subsurface interfaces, i.e.

( )
( ) 4 shows the average absolutely relative errors of five partial derivatives computed by the subdomain-Chebyshev spectral methods and finite-difference approaches.From these results, four facts are observed: firstly,

D
given by Equations ( 7) and ( 8 4).This is because of linear transformation from ζ to z (see Equation ( 11)); thirdly, the two algorithms of FDM2 and SCSM2 performs "0" denotes the scheme that applies Lagrange interpolations to free-surface topography and subsurface interfaces."1" and "2" stand for the algorithms using and not using high-order derivatives in coordinate transformations, respectively.better convergences than the algorithms of FDM1 and SCSM1 in computations of the second derivatives, because the former does not require high-order derivatives ( )  and actually reduces the order of smooth- ness in the coordinate transformations; finally, increasing the points in the subdomains (larger than five), the subdomain Chebyshev spectral methods significantly improve accuracies of the finite-difference approaches, and show achievement of almost one-order higher accuracies than the finite-difference methods.
In the 3D case (see Figure 2), we once again applied the subdomain Chebyshev spectral method to compute nine partial derivatives ( ) and their absolutely relative errors.Figure 5 shows an example in which seven points were also employed for the subdomain Chebyshev differentiation vector, and 10 × 10 points were applied to 2D cubic spline interpolations defining the free-surface topography  shows that most of the errors appear at the boundaries of the subdomains and places of the rough free-surface and subsurface interfaces; the maximum absolutely relative errors are less than 0.7% for the second derivatives.These performances show that the subdomain Chebyshev spectral method is an accurate tool for spatial derivatives in such a 3D geological model.
Figure 6 gives overall comparisons of the subdomain Chebyshev spectral method with the finite-difference approaches in the 3D case.The average absolutely relative errors of nine spatial derivatives are plotted against points of the differentiation vectors.These results exhibit similar performances to the 2D case, e.g. the convergence curves of z f ∂ and zz f ∂ have slight difference between two methods due to the linear transformation Here, the integers "1" and "2" denote the algorithms involving and excluding the high-order derivatives of spline interpolations.
from ζ to z.Other curves demonstrate that, as the points of differentiation vectors are larger than 5, the sub- domain Chebyshev spectral method significantly improve accuracies of the finite-difference approaches.Therefore, through the 3D experiments, we once again show that the subdomain Chebyshev spectral method is superior to the finite-difference approximation in numerical differentiation with curved coordinate system.

Conclusions
We have demonstrated a subdomain Chebyshev spectral method to numerical differentiations in a curved coordinate system, which may be employed to seek a "strong" solution of PED defined in ageological or a geophysical model.The new method employs a set of non-linear local coordinate transformations, global spline interpolations and subdomain Chebyshev points to make the discretization grid automatically match free-surface topography and subsurface interfaces, and guarantees the numerical solution converging to the "strong" solution of partial differential equations.The new method possesses advantages from the traditional finite-difference approaches, i.e. simple discretization procedure and sparse assembled matrix, but it has stronger capability to deal with complex geometries of geological models and better accuracies in numerical differentiations.Methodologically speaking, the subdomain Chebyshev spectral method may replace the finite-difference approaches in geological or geophysical computer modelling and enables us to easily hand free-surface topography and significantly improve accuracy of modelling results.2D and 3D synthetic experiments demonstrate the accuracies and capabilities of the subdomain Chebyshev spectral method.As points of the differentiation vectors are large than 5, the subdomain Chebyshev spectral method can reach one-order higher accuracy in numerical differentiations than the traditional finite-difference approaches in curved coordinates systems.Two algorithms of the subdomain Chebyshev spectral method (SCSM1, SCSM2) are both valid for numerical differentiations in PDE, but the algorithm SCSM2 is more efficient in computations than SCSM1, due to exclusion of computing the high-order derivatives of the coordinate transformations.
the Jacobian and Hessian matrices respectively, both can be analytically calculated from the given Equation (2); the number of points in the k ξ -direction and l ξ -directions, e.g. in Cartesian system, if Equation (2) is i become the standard finite-difference or spectral differentiation vectors[23].According to Equation (9), the discrete forms of the field quantity F are obtained, i.e.
are the free-surface topography and subsurface interfaces if they are present.If there is not any physical interface in the model domain, are pure mathematical boundaries of the subdomains.In each subdomain, we write Equation (2) into the following forms: of the subdomains.Spline interpolation theory has shown that arbitrary scattered or regularly-distributed free-surface topography ( ) ), we obtain the second differentiation vectors with respect to the global coordinates { } D D D .These substitutions are more or less tedious and involve computations of the second derivatives their components so as to clearly indicate the relationship between xy D and { } , x y D D .Calculating yx D F , gives the same result as Equation (21).This shows that the second differentiation vector xy D given by Equation (21) satisfies Clairaut's theorem xy yx = D D , and can be obtained by multiplication of two first differentiation vectors x D and y D .The advantage of Equation (

Figure 1 .
Figure 1.Illustration of the 1D subdomain Chebyshev differentiation scheme.In each subdomain, five Chebyshev points are applied.At a boundary point, half Chebyshev points are taken from the neighbour subdomains (see the points in "boundary domain").

Figure 2 .
Figure 2. 2D and 3D geological models used for numerical experiments.The subdomain Chebyshev differentiation grids are plotted in the backgrounds and applied to compute the partial derivatives.
column) and their absolutely relative errors (right column) of the 2D model, computed by the presented method.In these experiments, we applied seven Chebyshev points ( ) and vertical directions in each subdomain, and 11 points for the 1D spline interpolation functions matching the free-surface topography and subsurface interfaces.The left column shows that x f ∂ and xz f ∂ are odd functions of x; z f ∂ , xx f ∂ and zz f ∂ are even functions of x; and xx f ∂ equals zz f ∂ , all of which are predicted from exact derivatives of Equation (22).The right column demonstrates that accuracies of the derivatives with respect to the vertical direction, i.e. z f ∂ and zz f ∂ are much higher than the derivatives with respect to horizontal direction, i.e. x f ∂ , xx f ∂ and xz f ∂

Figure 3 .
Figure 3. First and second partial derivatives and their absolutely relative errors computed by subdomain Chebyshev spectral method for the 2D geological model shown in Figure 2. Seven Chebyshev points were applied to each subdomain.
); secondly, the convergence curves of the derivatives z f ∂ and zz f ∂ computed by the two methods are slightly different (see z f ∂ -diagram and zz f ∂ -diagram in Figure

Figure 4 .
Figure 4. Convergence curves of the 2D geological model.FDM0, FDM1 and FDM3 are three finite-difference methods.SCSM1 and SCSM2 are two subdomain Chebyshev spectral methods."0"denotes the scheme that applies Lagrange interpolations to free-surface topography and subsurface interfaces."1" and "2" stand for the algorithms using and not using high-order derivatives in coordinate transformations, respectively.
. In this figure, we give two panels: (a) and (b), presenting the computed derivatives and absolutely relative errors, respectively; both have three sections at y = −20, 0, +20 km showing the patterns of the derivatives and distributions of the absolutely relative errors in the model domain.From these results, one can see that patterns of the numerical results are completely consistent with the analytic derivatives, e.g.

Figure 5 (
a) show that xz f ∂ is an even function of x and y and that yz f ∂ is an odd function of x and y; meanwhile, xx f ∂ equals yy f ∂ .Figure 5(b)

FDM1Figure 6 .
Figure 6.Convergence curves of four numerical differentiation schemes: SCSM1 and SCSM2 are two algorithms of subdomain Chebyshev spectral method, FDM1 and FDM2 are two finite-difference approaches.Here, the integers "1" and "2" denote the algorithms involving and excluding the high-order derivatives of spline interpolations.
) must be smooth (differentiable) in the domain.In the next section, we present a set of the local coordinate transformations which are applicable for arbitrary free-surface topography and make the numerical differentiation vectors ), are applicable only if the Jacobian ( )