3-D velocity models — transformation from general to natural splines

The functions describing material parameters and structural interfaces in velocity models are frequently represented by splines. The general cubic splines differ from the natural cubic splines by the boundary conditions at the outermost gridpoints. The general cubic splines have a general curvature at the outermost gridpoints used for interpolation, whereas the natural splines have a zero normal curvature at the outermost gridpoints. It is thus very useful to employ a simple algorithm for the transformation between the general and natural splines. The transformation from the natural to general (bi–) (tri–) cubic splines is straightforward, because the natural splines represent a special case of the general splines. This paper is devoted to the algorithm of transformation from the general to natural (bi–) (tri–) cubic splines. We present the formulae necessary for the transformation together with their derivation. We illustrate the presented formulae on the example of fitting a 1–D quadratic function by natural cubic splines, and on the example of a velocity model of a layered structure with two 3–D structural interfaces.


INTRODUCTION
Numerical modelling of the propagation of seismic waves requires a mathematical description of the medium through which the waves propagate. The medium is described by means of the space distribution of some material parameters, for example, the values of seismic wave velocities, or the values of elastic parameters, the values of density, loss factors, etc. We refer to this space distribution as the velocity model. For some numerical methods, it is sufficient to specify the values of the material parameters only at discrete points in the space, usually at the gridpoints of a sufficiently dense grid. Other numerical methods require the values of the material parameters to be known at each point of the model volume, and each numerical method usually has specific requirements in respect of the material parameters.
For example, for ray tracing, the values of the material parameters must be defined everywhere in the velocity model, and, moreover, they should be continuous up to the second-order partial derivatives. The discontinuities of the material parameters may occur only across so-called structural interfaces. The structural interfaces must also be continuous, and should be described by one or several smooth surfaces, continuous again up to the second-order partial derivatives.
The smooth surfaces in the velocity models prepared for the application of ray methods are often described by splines. The authors of this paper use the natural bicubic or tricubic splines with tension by Cline (1974a,b;1981). The general cubic splines differ from the natural cubic splines by the boundary conditions at the outermost gridpoints. The general cubic splines have a general curvature at the outermost gridpoints used for the interpolation, whereas the natural splines have a zero normal curvature at the outermost gridpoints. Analogously for mixed derivatives f 1122 at the corner gridpoints in 2-D, f 1122 , f 1133 f 2233 at the edge gridpoints in 3-D, and f 112233 at the tip gridpoints in 3-D. The number of parameters thus equals the number of gridpoints for the natural splines, whereas the number of parameters of the general splines also includes the normal curvatures at the outermost gridpoints and the above mentioned mixed derivatives, i.e. it corresponds to the number of gridpoints of the grid increased by 2 additional gridpoints in each direction. It is thus very useful to have a simple algorithm for the transformation between the general and natural splines.
The transformation from the natural to general (bi-) (tri-) cubic splines is straightforward: all normal curvatures at the outermost gridpoints and the above mentioned mixed derivatives, which are implicitly zeros in the natural splines specification, are explicitly specified as zeros in the general splines description. We did not know the algorithm for the transformation from the general to natural (bi-) (tri-) cubic splines, and we thus propose it in Section 2. In Section 3, we illustrate the proposed algorithm by fitting a 1-D quadratic function by natural cubic splines. Section 4 then illustrates the proposed formulae using the example of a velocity model of a layered structure with two 3-D structural interfaces.

TRANSFORMATION FROM THE GENERAL TO NATURAL SPLINES
Assume that function f is specified by the general (bi-) (tri-) cubic splines, whose coefficients are given at the gridpoints of a rectangular, not necessarily regular, 1-D, 2-D or 3-D grid. If we wish to specify the function by natural splines, we follow the procedure described below.
First we generate the values of the function and its partial derivatives at the gridpoints of the original 1-D, 2-D or 3-D grid. We thus generate 3, 9 or 27 values at each gridpoint of the original grid depending on the dimensionality of the grid: (1) We then extend the original 1-D, 2-D or 3-D grid by 2 additional gridpoints in the direction of each grid axis, situated outside of the original grid, see Fig. 1. The distances of the new gridpoints from the original grid may be arbitrary. We denote them h 1− , h 1+ , h 2− , h 2+ , h 3− , h 3+ , and introduce up to 6 vectors Derivatives (1) are continuous at the outermost gridpoints of the original grid and are thus known, whereas the homogeneous third-order derivatives and mixed derivatives containing homogeneous third-order derivatives may be discontinues there and are to be determined.
If x is an outermost gridpoint, the values at the neighbouring additional gridpoint x+h 1 , where h 1 stands for h 1− or h 1+ , are given by cubic function where the values of f (x), f 1 (x) and f 11 (x) are given. We now determine the value of f 111 (x). For the natural splines, the second-order derivative at x+h 1 should vanish, Stud. Geophys. Geod., 63 (2019) Relation (3) with condition (5) then reads If x is an outermost gridpoint, the values at the neighbouring additional gridpoint x + h 1 , where h 1 stands for h 1− or h 1+ , may be evaluated using Eq. (6), and analogously for x+h 2 or x+h 3 . If x is a corner gridpoint in 2-D or an edge gridpoint in 3-D, the values at point x+h 1 +h 2 outside the original grid are given by bicubic function where the values of and f 1122 (x) are given. We now determine the values of other partial derivatives. For the natural splines, the second-order derivative should vanish for given h 1 and all values of h 2 within the new outermost grid cell, which yields conditions Analogously, the second-order derivative should vanish for given h 2 and all values of h 1 within the new outermost grid cell, which yields conditions Equation (7) with conditions (9) reads Equation (12) with conditions (11) finally yields which may also be expanded as . If x is a corner gridpoint in 2-D or an edge gridpoint in 3-D, the values at the neighbouring additional corner or edge gridpoint x+h 1 +h 2 may be calculated using relation (14), or analogously for x+h 1 +h 3 or x+h 2 +h 3 . If x is a tip gridpoint in 3-D, the values at point x + h 1 + h 2 + h 3 outside the original grid are given by tricubic function where the values of f ( , f 12233 (x) and f 112233 (x) are given. We now determine the values of other partial derivatives. For the natural splines, the second-order derivative should vanish for given h 1 and all values of h 2 and h 3 within the new outermost grid cell, which yields conditions , j = 0, 1, 2, 3 , k = 0, 1, 2, 3 . (17) Analogously, the second-order derivative should vanish for given h 2 and all values of h 1 and h 3 within the new outermost grid cell, which yields conditions and the second-order derivative should vanish for given h 3 and all values of h 1 and h 2 within the new outermost grid cell, which yields conditions Equation (22) with conditions (19) then reads and Eq. (23) with conditions (21) finally yields which may also be expanded as . If x is a tip gridpoint in 3-D, the values at the neighbouring additional tip gridpoint x+h 1 +h 2 +h 3 may be calculated using relation (25).

SIMPLE EXAMPLE OF FITTING A 1-D QUADRATIC FUNCTION BY NATURAL CUBIC SPLINES
Assume that we wish to specify quadratic function f (x) = x 2 in the interval x ∈ −1, 1 .
For the general cubic splines, we may specify values f = 1, 0, 1 at gridpoints x = −1, 0, 1, with the corresponding boundary conditions f 11 = 2, 2 for the secondorder derivative at outermost gridpoints x = −1, 1, see the black curve in Fig. 2. Note that we could obtain an equal quadratic function if we specified values f = 4, 1, 0, 1, 4 at gridpoints x = −2, −1, 0, 1, 2, with the corresponding boundary conditions f 11 = 2, 2 for the second-order derivative at outermost gridpoints x = −2, 2, see the yellow extension of the black curve in Fig. 2. x -2 -1 0 1 2 f 1 2 3 4 Fig. 2. Fitting 1-D quadratic function f (x) = x 2 by cubic splines. Black: The quadratic function in interval −1, 1 may exactly be fitted by the general cubic splines with respective boundary conditions. Yellow: The quadratic function may be extended to interval −2, 2 and may again be exactly fitted by the general cubic splines with respective boundary conditions. Blue: The quadratic function fitted by natural cubic splines at points x = −1, 0, 1. Green: The quadratic function fitted by natural cubic splines at points x = −2, −1, 0, 1, 2. Red: The natural cubic splines exactly fitting the quadratic function in interval −1, 1 according to Eq. (6).
In conclusion, if we specify quadratic function f (x) = x 2 at gridpoints x = −1, 0, 1 and interpolate the values using the natural cubic splines, we obtain the blue curve in Fig. 2, which considerably differs from the correct black one. If we specify the values at gridpoints x = −2, −1, 0, 1, 2 according to Eq. (6) and interpolate the values using the natural cubic splines, we obtain exactly the black curve f (x) = x 2 in the interval x ∈ −1, 1 .

EXAMPLE OF A VELOCITY MODEL WITH STRUCTURAL INTERFACES
Velocity model L7, provided by the Institut Français du Pétrole (IFP), is a velocity model of a real geological structure. Velocity model L7 is composed of layers separated by curved interfaces, see Fig. 3. The Earth surface is flat. The presented version contains only the upper two of the several structural interfaces.
The structural interfaces were interpolated using the general bicubic splines in IFP, whereas the software package MODEL (Bucha and Bulant, 2017 ) used by the authors of this paper considers the natural (bi-) (tri-) cubic splines or analogous splines with tension by Cline (1974a,b;1981). We thus had to transform the given general bicubic splines to the natural bicubic splines using the procedure described in Section 2.
The velocity model is limited by a rectangular box in Cartesian coordinates. Coordinate x 1 extends from 3.130 to 11.10134 length units, probably kilometres, coordinate x 2 extends from 1.301 to 10.701. Coordinate x 3 points downwards and extends from -0.001 to 3.000 length units. The top of the box has been selected 0.001 length units above the Earth surface situated at elevation x 3 = 0.000.
The layers in the velocity model are limited by 3 surfaces numbered 1, 2, 3 from top to bottom. The velocity model is composed of the free space above the flat Earth surface (surface No. 1), and of 3 layers numbered 1, 2, 3 from top to bottom. The surfaces are specified by implicit equations f (x 1 , x 2 , x 3 ) = 0 with functions f continuous up to the second-order partial derivatives. The functions f describing the smooth surfaces are selected in the form of f (x 1 , x 2 , The Earth surface is flat, x 3 =constant, at the elevation of 0 length units. Surface 2 was interpolated in the form of x 3 = w(x 1 , x 2 ) by the general bicubic splines in the regular rectangular grid of 48×48 gridpoints in IFP. The grid was then extended to 50×50 gridpoints for the interpolation by natural bicubic splines used in the MODEL package. For description of surface 2 by natural splines, the elevations at the original 48 × 48 grid were used without modification, and the elevations at the newly added gridpoints were calculated using the equations presented in Section 2.
Surface 3 is specified in the same way and at the same grid as surface 2. The data for surface 3 has also been transformed from the general bicubic splines used in IFP to the natural bicubic splines. The structural interfaces 2 and 3 of velocity model L7 are displayed in Fig. 3.
For a detailed description of the related input data for the MODEL package refer to Bulant and Klimeš (1996). For an example of two-point ray tracing of reflected rays in model L7 refer to Bulant (1999).

CONCLUSIONS
The transformation from the natural to general (bi-) (tri-) cubic splines is straightforward: all normal curvatures at the outermost gridpoints and the above corresponding derivatives at the edge and corner gridpoints, which are implicitly zeros in the natural splines specification, are explicitly specified as zeros in the general splines description.
In this paper, we have presented the algorithm for the transformation from the general to natural (bi-) (tri-) cubic splines. We thus see that the general splines and the natural splines are equivalent for the corresponding boundary conditions. However, they naturally may differ for other boundary conditions for the general splines, analogously as two general splines differ for different boundary conditions. The boundary conditions for the general splines thus correspond one-to-one to the values of the function specified at the additional outermost gridpoints of the new grid for the natural splines, which is extended by 2 additional gridpoints in the direction of each grid axis with respect to the original 1-D, 2-D or 3-D grid used for the general splines. This correspondence is expressed by relations (6), (13), (14), (24) or (25), respectively. The general splines and the natural splines are thus equivalent if the boundary conditions are carefully treated. In the case of inverse problems, a small formal advantage of the natural splines may consist in the number of spline parameters equal to the number of gridpoints.
The proposed method may be used, e.g., to specify a quadratic function in terms of the natural cubic splines, which has been demonstrated in Section 3. We have then illustrated the method on the example of modelling the structural interfaces in a 3-D velocity model in Section 4.