Numerical simulation of 2-D seismic wave propagation in the presence of a topographic ﬂuid–solid interface at the sea bottom by the curvilinear grid ﬁnite-difference method

This study simulates seismic wave propagation across a 2-D topographic ﬂuid (acoustic) and solid (elastic) interface at the sea bottom by the ﬁnite-difference method (FDM). In this method, seismic waves in sea water are governed by acoustic wave equations, whereas seismic waves in solid earth are governed by elastic wave equations. The ﬂuid–solid interface condition is implemented on the interface. Body-conforming grids are used to ﬁt the topographic ﬂuid– solid interface which naturally avoids spurious diffractions due to staircase approximation. A collocated-grid MacCormack FDM is utilized to update the waveﬁelds in the ﬂuid and solid media. The ﬂuid–solid interface condition is explicitly implemented by decomposing the velocity and stress components to the normal and tangential directions with respect to the interface within a fourth-order Runge–Kutta time-marching scheme. The algorithm solutions for both ﬂat and topographic ﬂuid–solid interface models are compared with analytical solutions and spectral element solutions to validate the proposed method. Results show a suitable agreement with the reference solutions and hence conﬁrms the validity of this method. The proposed FDM enforces the numerical solutions to satisfy the exact interface condition and it is more accurate than the conventional FDM that uses effective media parameters to approximate the interface condition.


I N T RO D U C T I O N
The fluid-solid problem in this study is about seismic waves that travel through the interface between inviscid fluid and solid media. The fluid-solid interface condition shows that the traction and normal component of particle velocity across the interface are continuous. By contrast, the tangential components of the particle velocity may be discontinuous. A chief benefits of correctly implementing the fluid-solid interface is the opportunity to investigate the physical phenomena associated with the interface. These phenomena include the generation and propagation of T-waves (e.g. de Groot-Hedlin & Orcutt 2001;Okal 2008;Jamet et al. 2013), Scholte and leaky Rayleigh waves (e.g. de Hoop & van der Hijden 1983, 1984Padilla et al. 1999;Zhu & Popovics 2004;Carcione & Helle 2004;Zheng et al. 2013) as well as seismic scattering caused by the seafloor (e.g. Robertsson & Levander 1995;Greaves & Stephen 2000). A thorough understanding of the Scholte wave can help seismic exploration studies infer the reliable shear wave velocity below the seafloor (e.g. Bohlen et al. 2004;Kugler et al. 2007;Muyzert 2007a,b). Another benefit is the ability of the approach to help in analysing specific seismic data close to the ocean (e.g. Noguchi et al. 2007) or data from multicomponent seismic ocean acquisi-tion. The acquired information can then be used to obtain enhanced reservoir images for characterization (Farfour & Yoon 2016).
Analytical (or semi-analytical) and numerical methods have been applied to solve the fluid-solid problem. Analytical (or semianalytical) methods, such as the generalized reflection/transmission coefficient method (GRTM; e.g. Luco & Apsel 1983;Chen 1990), can calculate synthetic waveforms for layered or quasi-layered fluid-solid media (Chen 1990(Chen , 1995Bouchon et al. 1995;Chen 1996;Ge & Chen 2007;Qian & Yamanaka 2012) with high accuracy and efficiency. However, handing complex media remains a significant challenge for these methods. By contrast, numerical methods can treat complex media more easily than the above-mentioned methods. However, the numerical methods cannot readily obtain accurate solutions for the fluid-solid interface problem. Hence many numerical methods have been further investigated to solve the problem that arises in scientific research and engineering studies. These numerical methods include the finite-difference method (FDM; e.g. Bayliss et al. 1986;Graves 1996;Moczo et al. 2002;Robbert et al. 2002;Lombard & Piraux 2004;Okamoto & Takenaka 2005), pseudospectral method (e.g. Tessmer et al. 1992;Hung & Forsyth 1998;Feng et al. 2007), finite-element method (FEM; e.g. Bermúdez et al. 1999;Zhang 2004), discontinuous Galerkin C The Authors 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. method (e.g. Käser & Dumbser 2008;Wilcox et al. 2010), finitevolume method (e.g. Voinovich et al. 2003) and spectral element method (SEM; e.g. Komatitsch et al. 2000;Chaljub et al. 2003).
Among these numerical methods, the FDM is a more direct approximation for numerical solutions of wave equations than the methods based on other formulations. The FDM also demands a lower computing cost and fits parallel computing. However, the accurate calculation of wavefield propagation solutions that involve interfaces, such as free surface topography and fluid-solid interface, is a challenge for finite difference. Previous studies on FDMs that involve fluid-solid and free surface interfaces may be catalogued into different classes in terms of grids (i.e. Cartesian grid and curvilinear grid). For Cartesian grid cases, the fluid-solid interfaces are discretized by stair-case approximation which may cause numerical scattering (Muir et al. 1992). The interface conditions are implemented by effective media (e.g. Robertsson 1996;Moczo et al. 2007), which assume that the tangential component of particle velocity is continuous across the interface and conflicts with the fluid-solid interface conditions. Another approach based on the Cartesian grid is the immersed-interface method (also called embedded-boundary method). This method can handle free-surface (e.g. Kreiss & Petersson 2006;Lombard et al. 2008) and fluid-solid interfaces (Lombard & Piraux 2004) in a way that allows high-order accuracy as conjectured by the approach's proponents. However, a recent study on free surface interface based on the aforementioned method shows that instability may arise in a 3-D scenario (Gao et al. 2015). For curvilinear grid cases, the grid lines align with the topographic interface. This occurrence not only avoids spurious diffraction naturally (e.g. Zhang & Chen 2006;Appelö & Petersson 2009;Tarrass et al. 2011;Puente et al. 2014), but renders implementing boundary conditions convenient. Recent studies on curvilinear grids may be distinguished from one another by the algorithms for implementing boundary conditions and grid schemes. These algorithms include the collocated scheme (also called non-staggered or unstaggered schemes) (e.g. Zhang & Chen 2006;Appelö & Petersson 2009;Tarrass et al. 2011;Zhang et al. 2012b), partly staggered scheme (e.g. Ely et al. 2008), and fully staggered scheme (e.g. Puente et al. 2014). Zhang & Chen (2006) utilized the collocated scheme and traction image method to implement the free surface boundary condition in the first-order wave equation. Appelö & Petersson (2009) also proposed a stable method on the collocated scheme to implement the free surface boundary condition in the second-order wave equation. Puente et al. (2014) employed a fully staggered Lebedev scheme and a mimetic approach to solve the free surface condition in the first-order wave equation. The curvilinear grids and collocated scheme are utilized to solve the fluid-solid interface problem for algorithm proposed in the following sections.
The current study extends the curvilinear grid finite difference (Zhang & Chen 2006) to simulate seismic wave propagation in 2-D fluid-solid media with a topographic interface. The curvilinear grids are adopted to conform to the topographic fluid-solid interface. Seismic waves in sea water are governed by first-order acoustic wave equations, whereas those in solid earth are by firstorder velocity-stress elastic wave equations. The fluid-solid interface condition is implemented explicitly by decomposing the velocity and stress components to normal and tangential directions on the interface. The computing results of this algorithm for flat and topographic fluid-solid interface models are compared by analytical solutions and spectral element solutions. The results show a suitable agreement with the reference solutions, which conforms the validity of the proposed method. The proposed FDM enforces the numerical solutions to satisfy the exact condition and it is more accurate than the conventional FDM that uses effective media pa-rameters (e.g. Graves 1996;Robbert et al. 2002;Moczo et al. 2007) to implement the interface condition; it is also demonstrated in the flat fluid-solid interface model section. Finally, a true sea-bottom model with heterogeneities in the media parameters is built to show wavefield complexity. Its synthetic waveforms are also compared with the SEM solutions. When the fluid-solid interface is complex, the proposed curvilinear grid FDM (Curv. FDM) would employ finer grids to discretize the interface to maintain the accuracy; this scenario makes the Curv. FDM require more computing resources than those of the Cartesian grid FDM. Nevertheless, the discontinuous curvilinear grid techniques proposed by Li et al. (2015) can be applied to this problem so that the proposed method can largely save on computational resources by utilizing fine grids for the interface area and coarse grids for other areas.

B A S I C E Q UAT I O N S
For convenience of description, the Einstein summation convention for repeated indices is adopted and the derivatives with respect to the y-axis (into the paper) are assumed to be zero.

Body-conforming grid and coordinate transformation
When the topographic fluid-solid interface is given, the curvilinear grids are adopted to discretize the geological model such that the interface grid line coincides with the topographic fluid-solid interface ( Fig. 1, left panel). This grid is called the boundary-conforming (or body-fitting) grid, which has been widely applied in seismic wave modelling (e.g. Fornberg 1988;Tessmer et al. 1992;Hestholm & Ruud 1994Zhang & Chen 2006;Appelö & Petersson 2009;Zhang et al. 2012b;Sun et al. 2016). The grid qualities, such as orthogonality of grid points, concentration of grid points and smoothness of the variations in the grid size, can be controlled during the generations (Thompson et al. 1985) by algebraic or partial differential equation methods. The body-conforming grid can be generated by software (e.g. ICEM CDF software of ANSYS Inc. and Gridgen) with specified grid qualities. The proposed algorithm can improve its computing efficiency with several special grid qualities (e.g. orthogonality of grid points). However, the algorithm can still work properly without such additions. For simplicity, the bodyconfirming grids in this study are generated by simple shrinking or stretching the grid lines in vertical direction. For finite difference computing, the generated grids in the physical space (x, z) Lifshitz 1959) z}, ρ f is the mass density of fluid media, κ is the bulk modulus of the fluid media, v f i is the particle velocity vector component of the fluid media and p is the hydrostatic pressure disturbance. The solid is assumed to be linear isotropic elastic and the elastic wave equation may be stated in the following form: where the indices i, j, k ∈ {x, z}, ρ s is the mass density of the solid media, λ and μ are the Lamé constants, v s i is the particle velocity vector component of the solid media and σ ij is the stress tensor component.
The 2-D wave equations of the fluid and solid media in the curvilinear coordinate (ξ , ζ ) can be derived from eqs (3) and (4) by the chain rules. Their full forms are presented as follows: and where χ = λ + 2μ.
Eqs (5) and (6) can also be expressed in the following matrix form: where W is the velocity-stress vector and A, B are coefficient matrices (see Appendix A for details). The momentum equations in eqs (5) and (6) can be written in the conservative form of the ∇ operator (Thompson et al. 1985;Zhang & Chen 2006) as follows: and where the tractions in the fluid and solid media are as follows: and |∇ζ | = (∂ x ζ ) 2 + (∂ z ζ ) 2 .

Fluid-solid interface condition
The fluid-solid interface condition that couples the acoustic wave equation (eq. 3) and elastic wave equation (eq. 4) is described as a dynamic boundary condition and a kinematic continuity condition where, indices i, j ∈ {x, z}, n i denotes the component of the normal vector of the fluid-solid interface and I ij is the unit second-order tensor. The dynamic boundary condition ensures that the traction across the fluid-solid interface is continuous and the kinematic continuity condition describes that only the normal component of the particle velocity across the interface is continuous. This condition implies that the fluid contact with the solid can only slide but not be pulled apart (Chapman 2004). Given the 2-D assumption, the dynamic boundary and kinematic continuity conditions are written in the 2-D curvilinear coordinate as follows: where σ nn , σ nτ , v s n and v f n are the wavefield components in the local coordinate (see Appendix B for details).

C U RV I L I N E A R G R I D F I N I T E -D I F F E R E N C E S C H E M E F O R T H E P RO B L E M
The curvilinear grid finite-difference scheme is a type of collocated grid MacCormack-type finite-difference operators with fourthorder Runge-Kutta time-marching scheme. This scheme has been applied in seismic wave modelling (Zhang & Chen 2006;Zhang et al. 2012a,b) and earthquake dynamics (Zhang et al. 2014b). A brief introduction to the curvilinear grid finite-difference scheme is presented in this section (see Zhang & Chen (2006) and Zhang et al. (2012a,b) for details).

MacCormack finite-difference operators
The collocated and staggered grids are two types of grid schemes, which are widely implemented in finite difference. The collocated grid is preferable over the staggered grid for boundary-conforming problem because this grid can avoid the complex interpolations of the staggered grid (Magnier et al. 1994;Igel et al. 1995). As such, the collocated grid can hold the overall accuracy. The MacCormacktype collocated-grid finite-difference scheme (Bayliss et al. 1986;Xie & Yao 1988;Vdfidus et al. 1992;Dai et al. 1995) is widely utilized in seismic modelling. Its wide usage is due to the inherent dissipation (artificial damping) of one-sided differences it introduces to eliminate the instability of collocated-grid central difference caused by odd-even decoupling. The MacCormack-type scheme alternately utilizes forward and backward one-sided differences at the stages of a multistage Runge-Kutta time integration. As a result, the central difference is recovered when the one-sided differences are added together at the last time stage. The DRP/opt MacCormack, 4-2 Mac-Cormack and 2-2 MacCormack operators (Hixon 1997) are used in the algorithm for spatial derivatives in this study. The forward and backward one-sided difference operators of DRP/opt MacCormack finite difference can be expressed as follows (i.e. ξ -axis): where L ξ represents the spatial difference with respect to the ξaxis; the subscript i is the grid index; the superscript F and B denote the forward and backward operators, respectively; and a −1 = −0.30874, a 0 = −0.6326, a 1 = 1.2330, a 2 = −0.3334 and a 3 = 0.04168 (Hixon 1997). The forward and backward finite onesided difference operators of the 4-2 MacCormack scheme can be expressed as follows (i.e. ξ -axis): where the coefficients a 0 = −7/6, a 1 = 8/6 and a 2 = −1/6 (Hixon 1997). The forward and backward one-sided difference operators of the 2-2 MacCormack scheme can be expressed as follows (i.e. ξ -axis): where a 0 = −1.0 and a 1 = 1.0 (Hixon 1997).
The right-hand side of eq. (7) can be discretized by the following equation: where F or B in the first and second superscripts overL indicates the forward or backward difference, respectively, with respect to the ξ -and ζ -axes, respectively.

Runge-Kutta time-marching scheme
The fourth-order Runge-Kutta scheme is utilized to update the wavefields from time step n to n + 1 by the following equation: where W n and W n+1 are the wavefields at time step n and n + 1, respectively, and α 2 = 0.5, α 3 = 0.5, α 4 = 1.0, β 1 = 1/6, β 2 = 1/3, β 3 = 1/3 and β 4 = 1/6. The aforementioned fourthorder Runge-Kutta scheme can be expressed as follows: The forward and backward one-sided differences are alternately applied at the stages of the Runge-Kutta integration, and the central difference is obtained at the last step of the scheme. The following two-step time marching sequence is employed in this study to avoid numerical bias.

Gaussian filter
Velocity wavefields may be discontinuous across the fluid-solid interface. Thus, the spatial velocity derivatives with respect to the interface for grid points that are exactly on the interface are always calculated by a one-sided (forward or backward MacCormack 2-2) difference operator to avoid utilizing values across the discontinuous interface. The velocity derivatives calculated in this manner are reasonable. However, the numerical experiments show that this effect may cause instability in the interface after a long computing time. A Gaussian filter (Kristek et al. 2010;Zhang et al. 2013) is then applied on the interface grid points along the ξ -axis to stabilize the computing for a long time. The filter can be expressed as follows: where α = 4.2, the grid spacing ratios are n = 1 and A is a constant for normalizing w k as follows: 2n k=−2n The filter works on the ξ -axis direction at each time step and for a wavefield variable with the grid index i 0 on the interface grids points. Its filtered value˜ is obtained as follows: The numerical experiments in Section 4 show that the filter can stabilize the computation and negligibly affects the physical solutions.

Implementation of the fluid-solid interface condition
The curvilinear grid line aligns with the fluid-solid interface, which is assumed to be parallel to the ξ -axis and labelled by the grid index k 0 in Fig. 2. The wavefield components and media parameters for the acoustic wave equation are defined on the grid points with k ≥ k 0 . By contrast, those for the elastic wave equation are defined on the grid points with k ≤ k 0 . The grid points on the fluid-solid interface possess the grid index with k = k 0 and contains the wavefield components and media parameters for the acoustic and elastic wave equations. A grid point on the fluid-solid interface is split into fluid and solid points to explicitly implement the fluid-solid interface condition (Fig. 2). The split fluid point contains the wavefield components and media parameters only for the fluid media.
On the contrary, the split solid point harbours those for the solid media. Notably, the split fluid and solid points for a grid point on the fluid-solid interface share the same spatial location. However, those points infinitely approximate the interface from the fluid and   '+', the receiver is on the fluid side of the interface. '−', the receiver is on the solid side of the interface. solid sides, respectively. For convenience, the grid points on the interface are called interface points and labelled by the grid index k 0 . Meanwhile, the split fluid and solid points along the ξ -axis are labelled as k + 0 and k − 0 , respectively, to avoid confusion. The grid points, which are three or more points away the interface points in the ζ -axis direction, are called interior grid points. By contrast, those that are one or two points away from the interface points in the ζ -axis direction are called near interface points. Special treatments should be applied to handle the interface points and their neighbouring points. The wavefield components in the fluid and solid media are updated by eq. (7), except that the velocity compo-nents on the interface and near interface grid points are updated by eqs (8) and (9).
The spatial derivatives with respect to the ξ -axis on all grid points and those with respect to the ζ -axis on the interior grid points are calculated by the DRP/opt MacCormack finite difference. The spatial derivatives of the velocity components with respect to the ζ -axis on the grid points with two points away from the interface points are calculated by the 4-2 MacCormack finite difference. By contrast, those one point away from the interface points are calculated by the 2-2 MacCormack finite difference. The derivatives of the velocity components with respect to the ζ -axis on the interface grid points are calculated by the forward and backward 2-2 MacCormack finite differences and fluid-solid interface conditions.
Given that the fluid-solid interface condition does not ensure that the spatial derivatives of tractions across the interface are continuous, we calculate their one side (fluid or solid side) spatial derivatives with respect to the ζ -axis near the fluid-solid interface from the same side information. Therefore, we utilize the following central symmetric extrapolated formula for tractions in the ghost points as follows: for the fluid side and for the solid side. The extrapolated formula is inspired by the traction image method proposed by Zhang & Chen (2006), but the values of the central points in this formula may change with time. The spatial derivatives of the traction components with respect to the ζ -axis in eqs (8) and (9) on the interface and near interface grid points are calculated by the DRP/opt MacCormack finite difference. The fluid-solid interface condition is implemented on the interface points within the fourth-order Runge-Kutta integration through the following processors. (1) h (1) is calculated in eq. (19) to obtain the components of h (1) in the local coordinate by rotation transformation in Appendix B.
(2) If h (1) is calculated by the L FF operator, then ∂ t v s n and ∂ t σ nn are replaced by ∂ t v f n and −∂ t pj, and ∂ t σ nτ = 0 is set to ensure the interface condition. (3) If h (1) is calculated by the L BB operator, then ∂ t v f n and ∂ t p are replaced by ∂ t v s n and −∂ t σ nn , and ∂ t σ nτ = 0 is set to ensure the interface condition. (4) The modified components of h (1) are rotated into the Cartesian coordinate to obtain h (1) * . (5) The Runge-Kutta processor is implemented to perform the same operations on h (2) , h (3) and h (4) and to obtain h (2) * , h (3) * and h (4) * . (6) W n + 1 * = W n + β 1 h (1) * + β 2 h (2) * + β 3 h (3) * + β 4 h (4) * is calculated, and the components of W n + 1 * are filtered by the Gaussian filter in eq. (24). (7) The filtered components of W n + 1 * are rotated into the local coordinate, and the values of v f n , v s n Downloaded from https://academic.oup.com/gji/article-abstract/210/3/1721/3869249 by guest on 29 July 2018 Figure 5. Relative errors of receivers 3 to 7 in Fig. 3 with different PPWs for the flat fluid-solid interface model. The left part shows the relative errors of the v x components, whereas the right part shows those of the v z components. The black line is a reference on whose slope is valued at 2. The slope of each line is approximately 2, which shows that the spatial accuracy of the algorithm is a second order.     and σ f nn , −p are replaced by (v f n + v s n )/2 and (σ nn − p)/2 to ensure the interface conditions. The modified components of W n + 1 * are rotated to the Cartesian coordinate to obtain the updated W n + 1 .
The lower-order finite-difference operators (eqs 16 and 17) are utilized to calculate the spatial derivatives of the particle velocity components with respect to the ζ -axis near or on the fluid-solid interface. These operators notably lower the spatial accuracy to the second order.

N U M E R I C A L E X P E R I M E N T S
We compare the results of the flat fluid-solid interface model with the GRTM solutions (Chen 1999) and those of the irregu-lar fluid-solid interface model with SEM solutions (Komatitsch et al. 2000) to verify the algorithm. The comparisons are quantitatively assessed by the global measurement of single-valued envelope misfit (EM), phase misfit (PM), envelope goodness-of-fit (EG) and phase goodness-of-fit (PG) of the time-frequency (TF) misfit criteria proposed by Kristeková et al. (2006Kristeková et al. ( , 2009. EM and PM show the disagreements of the two solutions in terms of the amplitude and phase, respectively. By contrast, EG and PG show the agreements of the two solutions in terms of the amplitude and phase, respectively. The TF misfit criteria can properly distinguish between the amplitude and phase errors. The criteria have been widely employed in seismology to evaluate numerical solutions (e.g. Pérez-Ruiz et al. 2007;Fichtner & Igel 2008;;   The v x snapshot is taken at 0.25 s with an explosive source located in the fluid media. By contrast, the v z snapshot is taken at 0.15 s with an explosive source located in the solid media. Figure 11. Snapshots of the v τ and v n components, which are rotated from their v x and v z components in Fig. 10 by formulas in Appendix B. The black star represents the source location, and the v-symbols with specified numbers represent the location of the receivers. The media parameters are written in the snapshots. Receiver 4 is located on the fluid-solid interface. The v τ snapshot is taken at 0.16 s with an explosive source located in the solid media. Meanwhile, the v n snapshot is taken at 0.30 s with an explosive source located in the fluid media.  Table 4.  Zhang et al. 2012b). The maximum value of the 'excellent' level given by Kristeková et al. (2009) for EM is 0.22, whereas that for PM is 0.2. EM or PM should be under their maximum 'excellent' levels. The minimum value of the 'excellent' level given by Kristeková et al. (2009) for EG or PG is 8, and these factors should be larger than 8. The Ricker wavelet is selected in the following experiments as the source time function. It can be expressed as where f c is the centre frequency and t 0 is the delay time. The points per wavelength (PPW) is computed by v s /(2.5f c h s ), where v s is the minimum velocity and h s is the maximum spatial interval value of the adjacent grid nodes in the x-and z-axes directions. The free surface of the fluid media is flat and the free surface boundary conditions are implemented by the traction image method (Zhang & Chen 2006). The outgoing waves can be absorbed by ADE CFS-PML (Berenger 1994;Chew & Liu 1996;Wang & Liang 2006;Zhang & Yang 2010;Zhang et al. 2014a) or the absorbing boundary layers (Cerjan et al. 1985). For simplicity, 10 absorbing boundary layers are built to absorb the outgoing waves. The grids for the curvilinear finite difference is generated by stretching or shrinking the grid lines in the z-axis direction. Given that fluid-solid wavefields across its interface may be discontinuous, the waveforms of the receivers on the interface are labelled with '+' and '−' to represent the fluid and solid sides of the interface, respectively.

Half-space model with the flat fluid-solid interface
The first test is a half-space with the flat fluid-solid interface at z = −200.0 m in the physical domain. An explosive source is located at (x: 450 m, z: −30 m) in the fluid media with f c = 35 Hz and t 0 = 0.05 s. The spatial locations of the receivers are shown in Table 1. The source-receiver geometry and media parameters are displayed in the snapshot in Fig. 3. The grid interval for this model is 2m in x-and z-axes directions. By contrast, the PPW for this computation is approximately 8.57.
The snapshots of the v x and v z components at 0.3 s are shown in Fig. 3. The synthetic waveforms of the receivers compared with the GRTM solutions are qualitatively shown in Fig. 4 and Figure 13. Synthetic v τ and v n waveforms for the receivers in Fig. 11 with PPW ≈5. The waveforms last up to 100 s by a time loop of 10 6 times. their quantitative comparisons are shown in Table 1. The comparison shows a suitable agreement between the GRTM solutions and curvilinear grid finite-difference solutions. This agreement is due to the maximum value of all the receivers' EMs in Table 1 lower than 7.68 per cent and the maximum value of all the receivers' PMs lower than 1.82 per cent. Moreover, the minimum value of all the receivers' EGs is larger than 9.25, whereas the minimum of all the receivers' PGs is larger than 9.81. These factors are all within the 'excellent' levels of the misfit criteria (Kristeková et al. 2009). The relative L 2 errors of the v x and v z components are calculated as follows to show the convergence rate of the spatial accuracy: The aforementioned equation is plotted in Fig. 5, which supposes that the curvilinear grid finite-difference solutions with PPW = 137.12 are the reference solutions. The said figure shows that the relative errors decrease as the PPW increases, and the algorithm has second-order accuracy. Given the lower-order finitedifference operators near and on the fluid-solid interface, the relative error of the receiver near the interface is larger than that far from the interface (Fig. 5). Fig. 3 clearly shows that the v x snapshot is discontinuous across the fluid-solid interface, whereas the v z snapshot is continuous. This phenomenon is a direct representation of the fluid-solid interface condition (eq. 12), which can also be clearly observed at the receivers on the fluid-solid interface (i.e. receivers 6, 8 and 9) in Fig. 4. The v x waveforms of receiver 6 (6+ and 6−), receiver 8 (8+ and 8−) and receiver 9 (9+ and 9−) are discontinuous across the interface, whereas their v z waveforms are continuous.
The conventional FDM in this study implicitly implements the fluid-solid interface condition by effective media parameters (e.g. Graves 1996;Robbert et al. 2002;Moczo et al. 2007). This FDM is also applied to calculate waveforms in this model. The waveforms computed by the conventional FDM are compared with the GRTM solutions in Fig. 6. Their quantitative comparisons are shown in Table 2. The v x waveforms (tangential component of the particle velocity) on the fluid-solid interface do not match with the GRTM solutions. The EMs and PMs of the v x waveforms on the fluid side of the interface (receivers 6+, 8+ and 9+) can be as large as 86.4 per cent. The EMs and PMs do not decrease as the PPW increases (Table 2). This implies that the fluid-solid interface condition (eq. 12) is incorrectly implemented.
The conventional finite-difference solutions of PPW = 32 are compared with the solutions of the curvilinear grid finite difference of PPW = 64 for the receivers on the grid points from receiver 5 to receiver 7 in this model. This comparison was performed to increase the understanding of the errors near the interface of the conventional FDM. Fig. 7 shows the EM and PM of the solutions on the receivers. The EM and PM values of the v z component almost remain the same across the fluid-solid interface. For their v x component, the solutions on the solid side can maintain a lower misfit as low as 2 per cent. By contrast, those on the fluid side can be as much as 84 per cent when the receiver is less than 0.12 wavelength to the fluid-solid interface.
For the flat fluid-solid model, we conclude that the proposed curvilinear FDM can solve the fluid-solid interface problem with second-order accuracy. Because the curvilinear FDM in this paper enforces the numerical solutions to satisfy the exact interface condition, it is also better than the conventional FDM that utilizes effective media parameters (e.g. Graves 1996;Robbert et al. 2002;Moczo et al. 2007) to implement the interface condition accurately.

Half-space model with canyon-shaped fluid-solid interface
The second experiment is a half-space with a canyon fluid-solid interface at approximately z = −200.0 m in the physical domain. An explosive source is located at (x: 300 m, z: −31.13 m) in the fluid media with f c = 10 Hz and t 0 = 0.15 s. The spatial locations of the receivers are shown in Table 3. Receiver positions are selected to coincide with the grid points, which are not evenly spaced in this canyon-shaped fluid-solid interface model. Hence, the receivers are not exactly evenly spaced. The source-receiver geometry and media parameters are shown in the snapshot in Fig. 8. The maximum grid interval for this model is approximately 6.8 m, and the PPW for this computation is approximately 8.82.
The v x and v z snapshots at 0.3 s are shown in Fig. 8. Synthetic waveforms of the receivers compared with the SEM solutions are qualitatively shown in Fig. 9 and their quantitative comparisons are listed in Table 3. The comparisons show a suitable agreement between the SEM solutions and curvilinear grid finite-difference solutions. This agreement is achieved because the maximum value of all the receivers' EMs in Table 3 is lower than 6.6 per cent, whereas the maximum value of all receivers' PMs is lower than 2.25 per cent. Moreover, the maximum value of all the receivers' EGs is larger than 9.36, whereas the maximum values of all receivers' PGs is larger than 9.77. These factors are all within the 'excellent' levels of misfit criteria (Kristeková et al. 2009).
For this model, we conclude that the Curv. FDM can properly handle the fluid-solid interface and obtain appropriate solutions with a density sampling of 8.82 PPW.

Sinusoidal fluid-solid interface model
This section demonstrates a half-space model with the sinusoidal shaped fluid-solid interface to test the algorithm in terms of numerical stability. An explosive Ricker source with f c = 30.0 Hz and t 0 = 0.05 s is located above and below the interface in two experiments. The source is located at (x: 450.0 m, z: −100.0 m) in the fluid media for the first test, but at (x: 450.0 m, z: −240.0 m) in the solid media for the second test. The receivers are located in a vertical line at x = 500 m with a uniform interval of 20 m from receiver 1 (z = −140 m) to receiver 7 (z = −260 m). The geometry of the sources and receivers in each experiment is shown in Fig. 10. Fig. 10 also shows the v x snapshot of the first test and v z snapshot of the second test. The v x and v z snapshots indicate that the topographic fluid-solid interface significantly affects the wavefields. Fig. 11 shows the v τ and v n snapshots, which are rotated from v x and v z components in Fig. 10 by the formulas in Appendix B. The grid line of the curvilinear coordinates on the interface aligns with the topographic fluid-solid interface. Thus, the v τ and v n components on the interface in Fig. 11 represent the tangential and normal velocity components with respect to the interface. The v τ component in Fig. 11 may be discontinuous across the fluid-solid interface, whereas the v n component is continuous. Fig. 12 shows the synthetic v τ and v n waveforms of the receivers in Fig. 11 with different PPWs (i.e. 20, 10 and 5). Their quantitative comparisons are listed in Table 4. The waveforms clearly converge as the PPW increases in Table 4. Receiver 4 is located on the interface, and the v τ waveforms of the fluid side (receiver 4+) and solid side (receiver 4−) are discontinuous across the interface in Fig. 12. By contrast, the v n components are continuous. Fig. 13 shows the synthetic v τ and v n waveforms of the receivers in Fig. 11. The waveforms last up to 100s by a time loop of 10 6 times. This result shows that the method can be stable for a long time.
We conclude that the solutions of the proposed algorithm for the complex fluid-solid interface not only converge as the PPW increases but also maintain stability for a long time in this complex fluid-solid interface model.

A realistic ocean-bottom model
This section builds a realistic ocean-bottom model with heterogeneities in the P-wave speeds and S-wave speeds, as well as densities to validate the proposed curvilinear FDM in a realistic ocean-bottom model by comparing with synthetic waveforms. The realistic 2-D ocean-bottom profile is from (E126.6, N20.5) to (E129.0, N20.5), which are obtained from the shuttle radar topography mission (SRTM) data. The media parameters of the realistic 2-D ocean-bottom model are obtained from Crust 1.0. The 2-D re-alistic models are compressed into 4.0 km in length and 2.0 km in depth (Fig. 14) to save on computing costs. Two heterogeneous media blocks are also added in the realistic model by an increase of 30 per cent of P-wave speed, S-wave speed and density to increase the model complexities. Fig. 14 shows the S-wave speeds of the model and its two heterogeneous media blocks. An explosive source with f c = 6Hz and t 0 = 0.2 s is located at (2.0 km, −0.15 km) in the ocean. The receivers are located in two lines along the x-axis on the free surface and in the ocean. The first line contains 26 receivers on the sea level and starts from receiver 1 (x = 1.0 km) to receiver 26 (x = 3.0 km) with an equal interval of 0.8 km. The second line contains 26 receivers in the fluid ocean at z = −0.4 km. This line also starts from receiver 27 (x = 1.0 km) to receiver 52 (x = 3.0 km) with an equal interval of 0.8 km. Fig. 15 shows the comparisons of the synthetic waveforms computed by Curv. FDM and SEM for receivers in line 1, whereas Fig. 16 shows those in line 2. The synthetic waveforms are complex in Figs 15 and 16, which are caused by the complex 2-D oceanbottom topography and heterogeneities in media parameters. The synthetic waveforms computed by Curv. FDM and SEM in Figs 15 and 16 can suitably match. This result hence validates the Curv. FDM in a realistic ocean-bottom model with heterogeneities in the P-wave speeds, S-wave speeds and densities.

C O N C L U S I O N S
This study extends the Curv. FDM (Zhang & Chen 2006) to simulate seismic wave propagation across a 2-D fluid-solid interface. Seismic waves in fluid media are governed by the acoustic wave equation, whereas those in solid media are governed by the elastic wave equation. The fluid-solid interface condition is implemented on the interface. The curvilinear grid is employed to describe the topographic fluid-solid interface because its grid line can align with the interface. As a result, the numerical scattering caused by staircase approximation can be avoided. The wavefield components on the fluid-solid interface may be discontinuous. Hence, the grid points on the fluid-solid interface are split into the fluid and solid points to explicitly implement the fluid-solid interface condition. This algorithm is validated by comparing with analytical and SEM solutions, and its robustness is also tested by long-time computations. The differences in the synthetic waveforms between this Curv. FDM and the previous FDM are also demonstrated in a flat fluid-solid interface model. This model shows the necessity of this new method for seismic wave propagation in the presence of a fluidsolid interface. The proposed FDM enforces the numerical solutions to satisfy the exact interface condition and it is more accurate than the conventional FDM that uses effective media parameters (e.g. Graves 1996;Robbert et al. 2002;Moczo et al. 2007) to approximate the interface condition. When the geometry of the fluid-solid interface is complex, the interface needs to be discretized by finer grids to maintain accuracy. However, the discretization would cost additional computing resources. Nevertheless, this issue can be solved by the discontinuous curvilinear grid technique (Li et al. 2015) that simultaneously utilizes a finer grid size and a smaller time step for the interface area, and a coarse grid size and a larger time step for areas out of the interface. Future studies can extend this algorithm to 3-D scenario.

A C K N O W L E D G E M E N T S
The 2-D SEM code is downloaded from the CIG website (http://geodynamics.org/ cig/software/specfem2d/). The SRTM data are downloaded from http://srtm.usgs.gov. The Crust1.0 data are downloaded from http://i-gppweb.ucsd.edu/∼gabi/crust1.html. The program package TF_MISFIT_GOF_CRITERIA for quantitative comparison is downloaded from http://www.nuquake.eu/ Computer_Codes/. This research was supported by the National Natural Science Foundation of China (Grant number: 41474037).