Acoustic and elastic wave equations are routinely used in geophysical and engineering studies to simulate the propagation of waves, with a broad range of applications, including seismology, near surface characterization, non-destructive structural evaluation, etc.

Finite difference methods remain popular choices for these simulations due to their simplicity and efficiency.

In particular, the family of finite difference methods based on the summation-by-parts operators and the simultaneous-approximation-terms technique have been proposed for these simulations, which offers great flexibility in addressing boundary and interface conditions.

For the applications mentioned above, surface of the earth

is usually associated with the free surface boundary condition.

In this study, we demonstrate that the weakly imposed free surface boundary condition through the simultaneous-approximation-terms technique can have issue when the source terms, which introduces abrupt disturbances to the wave field, are placed too close to the surface.

In response, we propose to build the free surface boundary condition into the summation-by-parts finite difference operators and hence strongly and automatically impose the free surface boundary condition to address this issue.

The procedure is very simple for acoustic wave equation, requiring resetting a few rows and columns in the existing

difference operators only.

For the elastic wave equation, the procedure is more involved and requires special design of the grid layout and summation-by-parts operators that satisfy additional requirements, as revealed by the discrete energy analysis.

In both cases, the energy conserving property is preserved.

Numerical examples

are presented to demonstrate the effectiveness of the proposed approach.

65M06, 35Q86, 86-08

§ INTRODUCTION

The phenomena of wave propagation has broad applications in geoscience and civil engineering such as earthquake studies, site characterization, non-destructive structural testing, etc.; see [1, 2, 3, 4] for a few examples.

Among the many discretization techniques that have been proposed for wave simulations, finite difference methods are popular choices for these applications due to their simplicity and efficiency.

In particular, a family of finite difference methods based on the summation-by-parts (SBP) operators, which dates back to [5], has found its success in wave simulations; see [6, 7, 8, 9, 10] for a few examples.

In recent years, several studies (see, e.g., [11, 12, 13, 14, 15]) have applied these methods to wave simulations on staggered grids, which is a popular choice of grid layout in seismic studies.

With the SBP operators, a common approach to impose boundary conditions is by appending penalty terms (i.e., weakly), sometimes referred to as the simultaneous-approximation-terms (SAT) technique; see [16].

One of the advantages of the SAT technique is its flexibility, since it enables the separation of concern between designing the discrete operators that approximate differential operators and addressing the the boundary conditions, which allows a modular procedure when approaching a discretization task and the reusability of designed operators.

In this study, we consider a particular simulation setting pertinent to the practical applications mentioned above, namely, a point source is placed near or on the boundary associated with the free surface boundary condition.

In seismic studies for example, this mimics the survey setting where controlled
seismic sources

In fact, this setting is also relevant to seismic studies involving natural (i.e., earthquake) sources, which happens deep beneath the surface.

This is because in commonly used imaging techniques such as reverse time migration or full waveform inversion, the signals or data misfits at the sensor locations, which are placed on or near the earth surface, are used as the source terms in a back propagation or adjoint simulation that follows the forward propagation driven by the earthquake source to eventually form the image or update of the subterranean media.

In this second simulation, the source location (i.e., the sensor location) is on or near the surface.

(e.g., those generated by airgun in marine environment or explosive in land environment) are used to inject energy to the subterranean media that drive the wave propagation.

Signals generated by these sources are collected at sensor locations to infer the subterranean media.

Similar application settings appear in other applications such as site characterization and non-destructive testing.

In this setting, we demonstrate that the weakly imposed boundary condition can have severe violation of the free surface boundary conditions and lead to inaccuracies in the simulation results.

To address this issue, we propose to build the free surface boundary condition into the SBP operators and hence strongly impose the boundary condition.

Both acoustic and elastic waves are considered.

For the acoustic case, the proposed procedure is very simple, involving only resetting a few rows and columns in the SBP operators to impose the free surface boundary condition strongly.

In the elastic case, the procedure is more involved and requires special design of the grid layout and SBP operators that satisfy additional requirement, as revealed by discrete energy analysis.

Numerical examples are presented to demonstrate the advantages of the strong approach over the weak approach for both acoustic and elastic waves, including improved accuracy and relaxed time step restriction.

In the remainder, we start with the 1D case to explain the procedures for both approaches and compare their outcomes in section <ref>.

We then extend the study to the 2D acoustic case in section <ref> and to the 2D elastic case in section <ref>.

Finally, we conclude in section <ref>.

§ THE 1D CASE

§.§ Background

For clarity, we start our discussion with the 1D case, where both the acoustic and elastic wave equations reduce to the following system:

[left=]alignat = 2
ρ∂v/∂t     =   ∂σ/∂x  +  s^v   ;

1/ρc^2 ∂σ/∂t     =   ∂v/∂x   +  s^σ  ,

where $\rho$ and $c$ are given physical parameters, standing for density and wave-speed, respectively; $s^v$ and $s^\sigma$ are source terms that drive the wave propagation; $\sigma$ and $v$ are the sought solution variables.

When interpreted as the acoustic wave system, $\sigma$ represents the negative of pressure;
when interpreted as the elastic wave system, $\sigma$ represents the stress.

In both cases, $v$ represents the velocity.

In seismic survey involving artificially generated sources,
$s^v$ and $s^\sigma$ are often abstracted as point sources.

In this study, we consider the simple case where the point source location coincides with a grid point.

The term $s^v$ that applies on velocity is often referred to as the directional source, which can be generated by, e.g., striking a sledgehammer on a metal plate on the ground;

the term $s^\sigma$ that applies on pressure or stress is often referred to as the compressional source, which can be generated using, e.g., airgun or explosive.

Supposing that wave system (<ref>) is defined over an interval $(x_L,x_R)$, its associated physical energy can be expressed as:

\begin{equation}
\label{1D_continuous_energy}
\mathscr E \ = \ \frac{1}{2} \int_{x_L}^{x_R} \rho v^2 dx \ + \ \frac{1}{2} \int_{x_L}^{x_R} \frac{1}{\rho c^2} \sigma^2 dx \, ,
\end{equation}

where the two terms represent the kinetic and potential energy in the system, respectively.

Omitting the source terms in (<ref>) for now,

taking the time derivative on both sides of (<ref>) and substituting the equations from (<ref>), we arrive at:

\begin{equation}
\label{1D_continuous_energy_derivative}
\frac{d \mathscr E}{d t} \ = \ - \ \sigma(x_L) \cdot v(x_L) \ + \ \sigma(x_R) \cdot v(x_R) \, .
\end{equation}

In other words, in the absence of source terms, time derivative of the physical energy depends on boundary data only. If free surface boundary condition is associated with both boundaries, i.e., $\sigma(x_L) = \sigma(x_R) = 0$, we have $\frac{d \mathscr E}{d t} = 0$, i.e., the physical energy is conserved.

The continuous wave system (<ref>) can be discretized in space using the SBP-SAT approach.

As an example, one such discretization on staggered grids has been presented in [12].

The underlying grid layout is illustrated in Figure <ref>, where the two sub-grids are referred to as the $N$-grid
and $M$-grid,

For brevity, we only sketch the outline of the SBP-SAT discretization below. For a more detailed discussion, the readers are referred to [12].

width=1, font=small,labelfont=small

Illustration of the grid layout underlying the discretization presented below.

Omitting the source terms and boundary condition, (<ref>) can be discretized in space using the SBP operators as follows:

[left=]alignat = 2
𝒜^M ρ^M d V/d t
=  𝒜^M 𝒟^N Σ  ;

𝒜^N β^N d Σ/d t
=  𝒜^N 𝒟^M V   ,

where $V$ and $\Sigma$ are the discrete solution vectors, occupying the $M$-grid and $N$-grid, respectively,

$\boldsymbol \rho^M$ and $\boldsymbol \beta^N$ are the discretization of physical parameters $\rho$ and $\nicefrac{1}{\rho c^2}$ on the $M$-grid and $N$-grid, respectively.

Matrices $\mathcal D^N$ and $\mathcal D^M$ are the difference operators; $\mathcal A^M$ and $\mathcal A^N$ are the norm matrices; together, they constitute the SBP operators.

The discrete energy associated with (<ref>) is defined as

\begin{equation}
\label{1D_discrete_energy}
E \ = \ \frac{1}{2} V^T \left(\mathcal A^M \boldsymbol \rho^M \right) V
\ + \ \frac{1}{2} \Sigma^T \left(\mathcal A^N \boldsymbol \beta^N \right) \Sigma \, ,
\end{equation}

where the norm matrices $\mathcal A^M$ and $\mathcal A^N$ play the roles of the integrals from (<ref>).

Taking the time derivative on both sides of (<ref>) and substituting the equations from (<ref>), we arrive at:

\begin{equation}
\label{1D_discrete_energy_derivative}
\frac{d E}{d t} \ = \ \Sigma^T \left[ \mathcal A^N \mathcal D^M \ + \ \left( \mathcal A^M \mathcal D^N \right)^T \right] V \, .
\end{equation}

To simplify the notation, we define

\begin{equation}
\label{1D_Q_definition}
Q \ = \ \mathcal A^N \mathcal D^M \ + \ \left( \mathcal A^M \mathcal D^N \right)^T .
\end{equation}

When designing the SBP operators, $Q$ is asked to satisfy the following property:

\begin{equation}
\label{1D_Q_property}
Q \ = \ - \ \mathcal E^L \left(\mathcal P^L\right)^T + \ \mathcal E^R \left(\mathcal P^R\right)^T ,
\end{equation}

where $\mathcal E^L$ and $\mathcal E^R$ are the canonical basis vectors whose first and last entries are 1, respectively, and therefore, select the first and last entries of $\Sigma$ when applied on it;

$\mathcal P^L$ and $\mathcal P^R$ are projection operators that can be applied on vector $V$ to provide approximations to $v(x_L)$ and $v(x_R)$, respectively.

We now can rewrite (<ref>) as

\begin{equation}
\label{1D_discrete_energy_derivative_boundary}
\frac{d E}{d t} \ = \ - \ \left( \mathcal E^L \right)^T \! \Sigma \cdot \left( \mathcal P^L \right)^T \! V
\ + \ \left( \mathcal E^R \right)^T \! \Sigma \cdot \left( \mathcal P^R \right)^T \! V \, .
\end{equation}

that $\left( \mathcal E^L \right)^T \! \Sigma$ and $\left( \mathcal P^L \right)^T \! V$ are approximations to $\sigma(x_L)$ and $v(x_L)$ and
that $\left( \mathcal E^R \right)^T \! \Sigma$ and $\left( \mathcal P^R \right)^T \! V$ are approximations to $\sigma(x_R)$ and $v(x_R)$, (<ref>) is the discrete equivalent to (<ref>).

Moreover, the free surface boundary condition

can be imposed weakly by appending two trailing terms i.e, SATs, to (<ref>), leading to the following system:

[left=]alignat = 2
𝒜^M ρ^M d V/d t
=  𝒜^M 𝒟^N Σ +  𝒫^L [ ( ℰ^L )^T  Σ- 0 ]
-  𝒫^R [ ( ℰ^R )^T  Σ- 0 ]   ;

𝒜^N β^N d Σ/d t
=  𝒜^N 𝒟^M V   .

One can easily verify that $\tfrac{d E}{d t} = 0$ for this system, i.e., the discrete energy is conserved, just as the continuous energy $\mathscr E$ is conserved under the free surface boundary condition; see (<ref>).

For the SBP-SAT approach described above, the design of SBP operators and the design of SATs can be detached from each other.

When designing the SBP operators, one only needs to be careful with leaving some suitable residual terms on the boundary; see (<ref>).

The boundary conditions are addressed later via the SATs and are not part of the consideration when designing SBP operators.

For a different type of boundary condition, one only needs to change the SATs accordingly and can reuse the same SBP operators.

This flexibility offers great advantage in practice.

Below we give a concrete example of the SBP operators (originally presented in [12]) in (<ref>) for the grid points near the left boundary (see Figure <ref>).

These operators correspond to the case of unit grid spacing, i.e., $\Delta x=1$.

For general cases, $\mathcal D^N$ and $\mathcal D^M$ need to be scaled by $\tfrac{1}{\Delta x}$ whereas $\mathcal A^N$ and $\mathcal A^M$ need to be scaled by $\Delta x$.

The bottom right corners of $\mathcal A^M$ and $\mathcal A^N$ are the mirror opposite of their top left counterparts; the bottom right corners of $\mathcal D^N$ and $\mathcal D^M$ are the negative mirror opposite of their top left counterparts.

They are omitted here to conserve space.

We note here that the set of operators presented in (<ref>) is not the unique choice.

Other choices exist; see [11] for another example.

However, numerical experiments reveal that the operators presented in [11] can lead to severe restriction on time step size.

Specifically, using the same 1D test with homogeneous medium that leads to (<ref>) and (<ref>), it is revealed that the maximally allowed time step size allowed by the operators presented in [11] is less than half ($\approx 43.3\%$) of that associated with the interior stencil; cf. (<ref>) and (<ref>).

For this reason, we omit this alternative choice in this study.

\begin{align}
\label{SBP_matrices_1D_extrapolating_A}
\mathcal A^M = \left[
\arraycolsep=3.2pt%\def\arraystretch{0.75}
\begin{array}{r r r r r r r r r r r r r r r r r}
\nicefrac{13}{12} &                 & \multicolumn{1}{r|}{} \\
& \nicefrac{7}{8} & \multicolumn{1}{r|}{} \\
&                 & \multicolumn{1}{r|}{ \nicefrac{25}{24} }  \\ \cline{1-3} \\[-1.75ex]
&  			    &    			& 		   1 &			  & \\
& 			    &    			& 			 & 			1 & \\[-0.75ex]
&  			    &    			& 			 &   		  & \ddots
\end{array}
\right],
\mathcal A^N = \left[
\arraycolsep=3.2pt%\def\arraystretch{0.75}
\begin{array}{r r r r r r r r r r r r r r r r r}
\nicefrac{7}{18} & 				   & 		   		& \multicolumn{1}{r|}{} \\
& \nicefrac{9}{8} &     			& \multicolumn{1}{r|}{} \\
&	               & 1  			& \multicolumn{1}{r|}{} \\
&    			   &    			& \multicolumn{1}{r|}{ \nicefrac{71}{72} }    \\ \cline{1-4} \\[-1.75ex]
&   			   &    			&    			&  1  			&    		  \\
&    			   &  				&    			&    			&  1  		& \\[-0.75ex]
&  			   &    			&   		    &    			&      		& \ddots
\end{array}
\right];
\allowdisplaybreaks
\\
\label{SBP_matrices_1D_extrapolating_D}
\mathcal D^N = \left[
\arraycolsep=3.2pt%\def\arraystretch{0.75}
\begin{array}{r r r r r r r r r r r r r r r r r}
-\nicefrac{79}{78} &  \nicefrac{27}{26} & -\nicefrac{1}{26}  &  \nicefrac{1}{78}  & \multicolumn{1}{r|}{0}      \\
\nicefrac{2}{21}  & -\nicefrac{9}{7}   &  \nicefrac{9}{7}   & -\nicefrac{2}{21}  & \multicolumn{1}{r|}{0}      \\
\nicefrac{1}{75}  &  0			        & -\nicefrac{27}{25} &  \nicefrac{83}{75} & \multicolumn{1}{r|}{ -\nicefrac{1}{25} } & \\ \cline{1-5}
&   					&  \nicefrac{1}{24}  & -\nicefrac{9}{8}   & \nicefrac{9}{8}   						 & -\nicefrac{1}{24} 		  & \\
&     				&    				 &  \nicefrac{1}{24}  & -\nicefrac{9}{8} 	  					 &  \nicefrac{9}{8}  		  & -\nicefrac{1}{24} \\[-0.5ex]
& 					& 					 & 					  & \multicolumn{1}{l}{\ddots} 				 & \multicolumn{1}{l}{\ddots} & \multicolumn{1}{l}{\ddots} & \multicolumn{1}{l}{\ddots}
\end{array}
\right],
\mathcal D^M = \left[
\arraycolsep=3.2pt%\def\arraystretch{0.75}
\begin{array}{r r r r r r r r r r r r r r r r r}
-2       		  &  3    			  & -1    			   &  0    			    & \multicolumn{1}{r|}{0}   \\
-1       		  &  1    			  &  0    			   &  0    			    & \multicolumn{1}{r|}{0}   \\
\nicefrac{1}{24} & -\nicefrac{9}{8}  &  \nicefrac{9}{8}   & -\nicefrac{1}{24}  & \multicolumn{1}{r|}{0} & \\
-\nicefrac{1}{71} &  \nicefrac{6}{71} & -\nicefrac{83}{71} &  \nicefrac{81}{71} & \multicolumn{1}{r|}{ -\nicefrac{3}{71} } & \\ \cline{1-5}
&    				  &  \nicefrac{1}{24}  & -\nicefrac{9}{8}   &  \nicefrac{9}{8}   		 			   & -\nicefrac{1}{24} 		    & \\
&     			  &    				   &  \nicefrac{1}{24}  & -\nicefrac{9}{8}   		 			   &  \nicefrac{9}{8}  		    & -\nicefrac{1}{24} \\[-0.5ex]
& 				  & 				   & 					& \multicolumn{1}{l}{\ddots} 			   & \multicolumn{1}{l}{\ddots} & \multicolumn{1}{l}{\ddots} & \multicolumn{1}{l}{\ddots}
\end{array}
\right].
\end{align}

With the above matrices, $Q$ takes the following form:

\begin{equation}
\label{Q_1D_matrix_weak}
Q =
\footnotesize
\left[
\arraycolsep=3.2pt%\def\arraystretch{0.75}
\begin{array}{r r r r r r r r r r r r r r r r}
-\nicefrac{15}{8} & \nicefrac{5}{4} & -\nicefrac{3}{8} \\
&    &    &    &    &    &    &    &    &    &    &    \\
&    &    &    &    &    &    &    &    &    &    &    \\
&    &    &    &    &    &    &    &    &    &    &    \\
&    &    &    &    &    &    &    &    &    &    & \nicefrac{3}{8} & -\nicefrac{5}{4} & \nicefrac{15}{8}
\end{array}
\right],
\end{equation}

i.e., only the first and last rows are nonzero and correspond to $-\left(\mathcal P^L\right)^T$ and $\left(\mathcal P^R\right)^T$, respectively.

§.§ Numerical experiments

In the following, we conduct a series of numerical experiments to illustrate the behavior of the above discretization in the presence of source terms near the boundary.

For clarity, homogeneous medium with unit density and wave-speed (i.e., $\rho=1~\text{kg}/\text{m}^3$ and $c=1~\text{m}/\text{s}$) is considered.

The free surface boundary condition is associated with both boundaries (i.e., $\sigma=0$ at $x_L$ and $x_R$).

A compressional point source is used to drive the wave propagation.

Temporal profile of the point source is specified as the Ricker wavelet with central frequency 5 Hz and time delay 0.25 s; see <cit.> for more detail about the source specification.

The maximal frequency in the source content is considered as 12.5 Hz, which leads to a minimal wavelength of 0.08 m.

The length of the interval $(x_L,x_R)$ is specified as 20 of such minimal wavelengths (i.e., $1.6$ m).

The stencil used for these experiments is the standard fourth-order staggered grid stencil in the interior (i.e.,
\nicefrac{ [\nicefrac{1}{24}, \enskip -\nicefrac{9}{8}, \enskip \nicefrac{9}{8}, \enskip -\nicefrac{1}{24}] }{ \Delta x }
see, e.g., [17, 18])
with adaptations near the boundaries, as illustrated in (<ref>), to satisfy the SBP property.

The staggered leapfrog scheme is used for time integration; see, e.g., [19] for more detail.

Denoting the Courant number as $C = \nicefrac{\Delta t}{\Delta x}$ (the unit wave-speed $c$ has been omitted), the CFL restriction for the aforementioned interior stencil and time integration scheme on unbounded domain stipulates

\begin{equation}
\label{Interior_CFL_limit}
C \leq C_{\max} = \nicefrac{6}{7}.
\end{equation}

Interested readers may consult [17] or <ref> for more detail.

Four experiments have been conducted, which use 10, 20, 40, and 80 grid points per minimal wavelength (ppw),

The grid layouts and source receiver locations in these experiments are illustrated in Figure <ref>.

Specifically, the point sources in these experiments are placed at a fixed physical location, which are 1, 2, 4, and 8 grid point(s) away from the left boundary $x_L$, respectively.

The time step length is specified as 2.5e-4 s for all four simulations, and the simulation time is specified as 6 s, which amounts to 24000 time steps.

width=1, font=small,labelfont=small

Grid layouts and source receiver locations for the 1D experiments.

The four grid lines, from top to the bottom, correspond to the left ends of the staggered grids with 10, 20, 40, and 80 ppw, respectively.

Source and receiver locations are placed at fixed physical locations for all four simulations, and indicated by the arrows and letters ($S$ for source; $R$ for receiver).

Time histories of the solution variable $\Sigma$ are recorded at three fixed physical locations ($R_0$, $R_1$, and $R_2$ in Figure <ref>) and displayed in Figure <ref>.

We note here that if the free surface boundary condition is strictly satisfied, $\Sigma$ at $R_0$ (i.e., the left boundary) should remain zero for the entire simulation duration for all four experiments.

From Figure <ref>, we observe that when the point source is placed too close to the free surface (in terms of the number of grid points; see Figure <ref>), which corresponds to the smaller ppw cases, the weakly imposed free surface boundary condition is severely violated, leading to significant inaccuracy in the recorded time histories.

Two additional experiments are presented in the Supplementary Material (<ref>supp_section_finer_grid and <ref>supp_section_interior_source) to demonstrate that such violation is not due to the lack of grid resolution and that such violation subsides when the source is placed far away from the free surface.

width=1, font=small,labelfont=small

width=1, font=small,labelfont=small

Time histories of $\Sigma$ at $R_0$.

width=1, font=small,labelfont=small

Time histories of $\Sigma$ at $R_1$.

width=1, font=small,labelfont=small

Time histories of $\Sigma$ at $R_2$.

Time histories of the solution variable $\Sigma$ at the three receiver locations depicted in Figure <ref> for the four experiments where the distances, in terms of number of grid points, between the source and the free surface boundary are 1, 2, 4, and 8, respectively.

Moreover, numerical experiments reveal that the CFL restriction for the above discretization setting (encompassing the difference operators, weakly imposed boundary terms, and time integration scheme) is

\begin{equation}
\label{1D_extrapolating_CFL_limit_weak}
C \leq C_{\max} \approx 0.6355,
\end{equation}

which is approximately 74.1% of that associated with the interior stencil on an unbounded domain; see (<ref>).

§.§ Strong imposition of the free surface boundary condition

To address the severe violation of the free surface boundary condition illustrated above, we propose to impose the free surface boundary condition strongly by incorporating it into the difference operators (cf. Remark <ref>).

While this can be achieved by starting anew the process of designing the SBP operators, it can also be achieved by adapting the existing operators originally designed with weakly imposed boundary conditions in mind.

We find this later approach easier to comprehend and illuminate the connections and contrasts between the weak and strong approaches better and will present it in the following.

Taking the SBP operators illustrated in (<ref>) and its associated grid layout (see Figure <ref>) as an example,
the first entry of the solution vector $\Sigma$, denoted as $\Sigma(1)$ and corresponding to the left boundary $x_L$, is updated by the first row of $\mathcal D^M$; see (<ref>).

From the free surface boundary condition at $x_L$, i.e., $\sigma(x_L) = 0$, we know that $\Sigma(1)$ is supposed to remain zero during the simulation.

We can simply set the first row of $\mathcal D^M$ to zero to achieve this.

On the other hand, since $\Sigma(1)$ will be multiplied by the first column of $\mathcal D^N$ to update $V$, we can also set the first column of $\mathcal D^N$ to zero to achieve the same effect.

In doing so, information of the free surface boundary condition is included in the simulation, and therefore, the SATs in (<ref>) are no longer needed.

For all other grid points, their stencils remain (effectively) unchanged and the accuracy of their approximations inherits from the original design of $\mathcal D^M$ and $\mathcal D^N$.

Finally, if we set both the first row of $\mathcal D^M$ and the first column of $\mathcal D^N$ to zero, it is easy to verify that the first row of $Q$ is also zero; see (<ref>) and (<ref>).

Performing the corresponding operations for the right boundary $x_R$ (i.e., setting the last row of $\mathcal D^M$ and the last column of $\mathcal D^N$ to zero), we have that

\begin{equation}
\label{Q_1D_property_strong}
Q = \boldsymbol 0.
\end{equation}

Recalling the discrete energy analysis result from (<ref>), we have that $\frac{d E}{d t} = 0$ because of (<ref>), i.e., the semi-discretized system (<ref>),

with operators modified as described above,
is also energy-conserving.

Conducting the same experiments as those in section <ref> with the free surface boundary condition imposed strongly, the recorded signals are displayed in Figure <ref>, from where we observe that the violation of the free surface boundary condition disappears and that results from the different simulations agree well with each other (cf. Figure <ref>).

width=1, font=small,labelfont=small

width=1, font=small,labelfont=small

Time history of $\Sigma$ at $R_0$.

width=1, font=small,labelfont=small

Time history of $\Sigma$ at $R_1$.

width=1, font=small,labelfont=small

Time history of $\Sigma$ at $R_2$.

Time histories of the solution variable $\Sigma$ at the three receiver locations depicted in Figure <ref>.

The experimental setup is identical to that for Figure <ref>.

The free surface boundary condition is imposed strongly.

Moreover, numerical experiments reveal that the CFL restriction for the above discretization with strongly imposed boundary conditions is

\begin{equation}
\label{1D_extrapolating_CFL_limit_strong}
C \leq C_{\max} = \nicefrac{6}{7},
\end{equation}

i.e., the same as that allowed by the interior stencil; see (<ref>) and compare to (<ref>).

For a cross validation of the simulation results from different discretizations, including the two discretizations described above (i.e., the weak and strong cases, respectively) and a finite element discretization, the readers are referred to the Supplementary Material <ref>supp_cross_comparison.

For the spectral radii information associated with the semi-discretized systems, which provide confirmation to the above claims on CFL restrictions, the readers are referred to the Supplementary Material <ref>supp_spectral_radii.

§ THE 2D ACOUSTIC CASE

Extension of the above procedure to the 2D acoustic wave equation is straightforward.

We can simply set the rows of the 2D equivalence of $\mathcal D^M$ that correspond to grid points on the free surface to zero and set the columns of the 2D equivalence of $\mathcal D^N$ that correspond to grid points on the free surface to zero to strongly impose the free surface boundary condition, provided that the pressure grid points are placed on the free surface.

The acoustic equation and semi-discretized systems are included in the Supplementary Material <ref>supp_acoustic_background and omitted here to conserve space. A numerical example is presented below to demonstrate the effectiveness of the strongly imposed free surface boundary condition, in comparison to the weak approach.

Specifically, homogeneous medium with unit density and wave-speed is considered, as in the 1D case. The source specification is also the same as in the 1D case, which leads to a minimal wavelength of 0.08 m.

The simulation domain is 6 by 6 of such minimal wavelengths (i.e., 0.48 m on both directions).

The free surface boundary condition is associated with all boundaries.

Three experiments are conducted for the strong and weak case each, which use 10, 30, and 50 grid points per minimal wavelength (ppw). The grid layout is illustrated in Figure <ref>.

The point source and receivers are placed at fixed physical locations at one third and two thirds of the simulation domain horizontally (i.e., 0.16 m and 0.32 m).

The source is placed at one grid point below the surface for the case where ppw is 10.

Three receivers are placed at the surface, at 3 and 6 grid points below the surface for the case where ppw is 10.

The time step length is specified as 5e-4 for all simulations, and the simulation time is specified as 3 s.

width=1, font=small,labelfont=small

Grid layout for the 2D acoustic wave equation.

Time histories of the solution variable $\Sigma$ at the three receiver locations are presented in Figures <ref> and <ref>, for the strong and weak case, respectively.

As a reminder, variable $\sigma$ and its discretization $\Sigma$ represents the negative of pressure for the acoustic case.

From these figures, we observe that the strongly imposed boundary condition delivers satisfactory results in all three simulations.

We note here that because of the free surface boundary condition, $\Sigma$ is supposed to be zero at the surface; see Figure <ref>.

On the other hand, when the point source is placed too close to the surface (in terms of the number of grid points in between), the weakly imposed free surface boundary condition is severely violated, leading to inaccuracies in the simulated results.

Such inaccuracy diminishes as the receiver depth increases.

Time histories of the solution variable $V_x$ and $V_y$ are included in the Supplementary Material <ref>supp_additional_figures_2D_acoustic.

Moreover, numerical experiments reveal that the CFL restriction for the strong and weak cases are $\nicefrac{6}{7}$ and 0.6355, respectively, the same as in the 1D case.

In other words, for the strong case, there is no additional penalty on the time step size compared to that from the interior stencil; see (<ref>), while for the weak case, there is a penalty of around 26% on the time step size allowed.

width=1, font=small,labelfont=small

width=1, font=small,labelfont=small

Time history of $\Sigma$ at the surface.

width=1, font=small,labelfont=small

Time history of $\Sigma$ at 3 grid points below the surface for the case ppw = 10.

width=1, font=small,labelfont=small

Time history of $\Sigma$ at 6 grid points below the surface for the case ppw = 10.

Time histories of the 2D acoustic experiments. The free surface boundary condition is imposed strongly.

width=1, font=small,labelfont=small

width=1, font=small,labelfont=small

Time history of $\Sigma$ at the surface.

width=1, font=small,labelfont=small

Time history of $\Sigma$ at 3 grid points below the surface for the case ppw = 10.

width=1, font=small,labelfont=small

Time history of $\Sigma$ at 6 grid points below the surface for the case ppw = 10.

Time histories of the 2D acoustic experiments. The free surface boundary condition is imposed weakly.

§ THE 2D ELASTIC CASE

Unlike the acoustic case, extension of the above procedure on strongly imposing the free surface boundary condition to the 2D elastic case is less straightforward and requires intricate design in the grid layout and the operators, which will be the focus of this section.

§.§ Background

The 2D isotropic elastic wave system is presented in the following:

[left=]alignat = 2

ρ∂v_x/∂t     =  ∂σ_xx/∂x + ∂σ_xy/∂y    ;

ρ∂v_y/∂t     =  ∂σ_xy/∂x + ∂σ_yy/∂y    ;

s_xxkl ∂σ_kl/∂t     =  ∂v_x/∂x;

s_xykl ∂σ_kl/∂t     =  1/2 ( ∂v_y/∂x + ∂v_x/∂y )   ;

s_yykl ∂σ_kl/∂t     =  ∂v_y/∂y,

where velocity components ($v_x$ and $v_y$) and stress components ($\sigma_{xx}$, $\sigma_{xy}$, and $\sigma_{yy}$) are the sought solution variables; density ($\rho$) and compliance tensor ($s_{ijkl}$) are given physical parameters.

In (<ref>), the Einstein summation convention applies to the subscript indices $k$ and $l$, which go through $x$ and $y$.

We note here that in (<ref>), the constitution relation is expressed via the compliance tensor, which is convenient for analysis and derivation.

Its equivalent form, where the constitution relation is expressed via the stiffness tensor, is convenient for implementation.

The analysis and derivation based on (<ref>) can be easily translated to its equivalent system.

To conserve space, we omit their link here. Interest readers may referred to <cit.> for more detail.

width=1, font=small,labelfont=small

Standard grid layout for the 2D elastic wave equation without considering boundaries.

As illustrated in Figure <ref>, when discretized on staggered grids using the standard grid layout, the normal stress components $\sigma_{xx}$ and $\sigma_{yy}$ are placed on the same location,

the horizontal velocity $v_x$ and vertical velocity $v_y$ are shifted from the normal stress components on $x$- and $y$- directions, respectively,

the shear stress component $\sigma_{xy}$ is shifted from the normal stress components on both directions.

Now considering the top boundary for example, the free surface boundary condition stipulates $\sigma_{xy} = \sigma_{yy} = 0$.

Strongly imposing these boundary conditions is non-trivial in the standard staggered grid setting for two reasons.

First, due to the grid staggering, the two stress components involved in the free surface boundary condition (i.e., $\sigma_{xy}$ and $\sigma_{yy}$) are not on the surface simultaneously.

Second, among the two normal stress components that share the same grid (i.e., $\sigma_{xx}$ and $\sigma_{yy}$) and hence share the same operators, one of them, namely, $\sigma_{xx}$, is not constrained by the free surface boundary condition. Strongly imposing the condition $\sigma_{yy} = 0$ by modifying the operators as described in section <ref> will lead to erroneous update to $\sigma_{xx}$.

§.§ Grid layout

In order to strongly impose the free surface boundary condition, we have to modify the grid layout slightly at the boundary.

We find the grid layout illustrated in Figure <ref> a suitable choice.

In Figure <ref>, the shear stress component $\sigma_{xy}$ is aligned with the $N$-grid points on both directions.

Again taking the top boundary for example, we have inserted an extra grid point occupied by the normal stress components $\sigma_{xx}$ and $\sigma_{yy}$ on the boundary.

Now both $\sigma_{xy}$ and $\sigma_{yy}$ are present on the top surface. This addresses the first hindrance mentioned above.

For the purpose of imposing the free surface boundary condition strongly, modifying the operators as described in section <ref> will still lead to erroneous updates to $\sigma_{xx}$ on the top surface.

However, noticing that $\sigma_{xx}$ is only used to update $v_x$; see (<ref>), which is not present on the top surface.

This means that the erroneous updates to $\sigma_{xx}$ will not have an effect on the rest of the simulation; $\sigma_{xx}$ on the top surface is essentially isolated from the solution variables at the other grid points.

This addresses the second hindrance mentioned above.

In fact, these extra grid points can be considered as only conceptually there, i.e., they can be omitted in practical implementations.

This is because, on one hand, $\sigma_{xx}$ on the top surface is never used.

On the other hand, $\sigma_{yy}$ is zero on the top surface and hence always has zero contribution to the variables it is used to update. This is equivalent to not taking its contribution by truncating the difference operators (instead of setting the rows and columns to zero as described in section <ref>).

The omission mentioned above lead to some inconsistency between how the grid columns for $\sigma_{xy}$ and $v_x$ and the grid columns for $\sigma_{xx}$, $\sigma_{yy}$, and $v_y$ are treated.

For the convenience of discussion, we refer to the former as $\mathbb N$-grid columns and the later as $\mathbb M$-grid columns henceforward.

Namely, on $\mathbb N$-grid columns, the operators are reset to zero but not truncated; on $\mathbb M$-grid columns, the operators are truncated and grid points on the surface are omitted.

To reflect this choice, these extra grid points are shaded in Figure <ref>.

This may cause slight confusion but brings the benefit that the grid points and hence the sizes of the operators on these two grid columns are the same, which greatly simplifies the analysis and implementation.

As we will see later, by design, the operators associated with these two grid columns can be made to be effectively the same, which is convenient for practical implementation.

width=1, font=small,labelfont=small

Grid layout for the 2D elastic wave equation suitable for strongly imposing the free surface boundary condition.

§.§ SBP operators

As alluded to in Remark <ref>, we still need to design the suitable SBP operators for this grid layout with the strongly imposed free surface boundary condition in mind.

We note here that although $\sigma_{yy}$ on the top surface is omitted in numerical implementation, it still needs to be formally included in the stencil. Otherwise, the effect of the boundary condition $\sigma_{yy} = 0$ will not be imposed.

To start, we conduct discrete energy analysis for the elastic wave system (<ref>) to gain insights on what are the requirements for these SBP operators.

The semi-discretized system is presented below:

[left=]alignat = 2

𝒜^V_x ρ^V_x d V_x/d t     =  𝒜^V_x 𝒟^Σ_xx_x Σ_xx + 𝒜^V_x 𝒟^Σ_xy_y Σ_xy   ;

𝒜^V_y ρ^V_y d V_y/d t     =  𝒜^V_y 𝒟^Σ_xy_x Σ_xy + 𝒜^V_y 𝒟^Σ_yy_y Σ_yy   ;

𝒜^Σ_xx S^Σ_kl_xxkl d Σ_kl/d t     =  𝒜^Σ_xx 𝒟^V_x_x V_x   ;

𝒜^Σ_xy S^Σ_kl_xykl d Σ_kl/d t     =  1/2 𝒜^Σ_xy ( 𝒟^V_y_x V_y + 𝒟^V_x_y V_x )   ;

𝒜^Σ_yy S^Σ_kl_yykl d Σ_kl/d t     =  𝒜^Σ_yy 𝒟^V_y_y V_y   ,

In (<ref>), superscripts such as $^{V_x}$ indicate the grids with which the underlying quantities or operators are associated; the Einstein summation convention still applies, but not to those appearing in these superscripts.

The 2D SBP operators appearing in (<ref>) are constructed from their 1D counterparts via tensor product as

𝒜^Σ_xy                            = ℕ𝒜^N_x ⊗ℕ𝒜^N_y ,
𝒜^V_x                                    = 𝕄𝒜^N_x ⊗ℕ𝒜^M_y ,
𝒜^V_y                                    = ℕ𝒜^M_x ⊗𝕄𝒜^N_y ,
𝒜^Σ_xx = 𝒜^Σ_yy = 𝕄𝒜^M_x ⊗ℕ𝒜^M_y ,

𝒟_x^Σ_xy = ℕ𝒟^N_x ⊗ℐ^N_y ,
𝒟_x^V_x         = 𝕄𝒟^N_x ⊗ℐ^M_y ,
𝒟_x^V_y         = ℕ𝒟^M_x ⊗ℐ^N_y ,
𝒟_x^Σ_xx = 𝕄𝒟^M_x ⊗ℐ^M_y ,

𝒟_y^Σ_xy = ℐ^N_x ⊗ℕ𝒟^N_y ,
𝒟_y^V_x         = ℐ^N_x ⊗ℕ𝒟^M_y ,
𝒟_y^V_y         = ℐ^M_x ⊗𝕄𝒟^N_y ,
𝒟_y^Σ_yy = ℐ^M_x ⊗𝕄𝒟^M_y .

Since the $\mathbb N$ and $\mathbb M$ grid rows and columns have different grid layouts, they are associated with different sets of 1D SBP operators.

In (<ref>) and (<ref>), the overhead symbols $\mathbb N$ and $\mathbb M$ are attached to distinguish.

The discrete energy associated with (<ref>) is

\begin{equation}
\label{Discrete_energy_2D_elastic_wave_system}
E \ = \ \tfrac{1}{2} V_i^T \left( \mathcal A^{V_i}  \boldsymbol{\rho}^{V_i} \right) V_i^{\phantom{T}}
\! + \,
\tfrac{1}{2} \Sigma_{ij}^T \left( \mathcal A^{\Sigma_{ij}}  S^{\Sigma_{kl}}_{ijkl} \right) \Sigma_{kl}^{\phantom{T}}\ ,
\end{equation}

Taking the time derivative on both sides of (<ref>) and substituting the equations from (<ref>), we arrive at

d E/d t
V_x^T [
𝒜^V_x 𝒟^Σ_xx_x
(𝒜^Σ_xx 𝒟^V_x_x)^T
] Σ_xx

V_y^T [
𝒜^V_y 𝒟^Σ_xy_x
( 𝒜^Σ_xy 𝒟^V_y_x )^T
] Σ_xy

V_y^T [
𝒜^V_y 𝒟^Σ_yy_y
( 𝒜^Σ_yy 𝒟^V_y_y )^T
] Σ_yy

V_x^T [
𝒜^V_x 𝒟^Σ_xy_y
( 𝒜^Σ_xy 𝒟^V_x_y )^T
] Σ_xy   .

Substituting the definitions from (<ref>) and (<ref>) and invoking the tensor product properties, we arrive at

d E/d t
V_x^T [
( 𝕄𝒜^N_x 𝕄𝒟^M_x )
( 𝕄𝒜^M_x 𝕄𝒟^N_x )^T
] Σ_xx

Σ_xy^T [
( ℕ𝒜^N_x ℕ𝒟^M_x )
( ℕ𝒜^M_x ℕ𝒟^N_x )^T
] V_y

V_y^T [
ℕ𝒜^M_x ⊗( 𝕄𝒜^N_y 𝕄𝒟^M_y )
𝕄𝒜^M_x ⊗( 𝕄𝒜^M_y 𝕄𝒟^N_y )^T
] Σ_yy

Σ_xy^T [
ℕ𝒜^N_x ⊗( ℕ𝒜^N_y ℕ𝒟^M_y )
𝕄𝒜^N_x ⊗( ℕ𝒜^M_y ℕ𝒟^N_y )^T
] V_x

If we have the following properties

ℕ𝒜^N_x    = 𝕄𝒜^N_x    = 𝒜^N_x

ℕ𝒜^N_y    = 𝕄𝒜^N_y    = 𝒜^N_y

( ℕ𝒜^N_x ℕ𝒟^M_x )
( ℕ𝒜^M_x ℕ𝒟^N_x )^T    = 0

( 𝕄𝒜^N_x 𝕄𝒟^M_x )
( 𝕄𝒜^M_x 𝕄𝒟^N_x )^T    = 0

( ℕ𝒜^N_y ℕ𝒟^M_y )
( ℕ𝒜^M_y ℕ𝒟^N_y )^T    = 0

( 𝕄𝒜^N_y 𝕄𝒟^M_y )
( 𝕄𝒜^M_y 𝕄𝒟^N_y )^T    = 0

in the 1D operators, we have that $\frac{d E}{d t} = 0$ following (<ref>) and the property of tensor product, i.e., the discrete energy is conserved, which mimics the behavior of the continuous elastic wave system.

We note her that while (<ref>) is expected from common practice, (<ref>) is the additional requirement related to the special grid layout, namely, the two sets ($\mathbb N$ and $\mathbb M$) of 1D SBP operators need to share the same norm matrices, which implies that the design of these two sets of operators need to be interwined.

Such operators do exist.

Moreover, it turns out that we can ask for an additional requirement on these operators, namely,

ℕ𝒟^N_x    = 𝕄𝒟^N_x    = 𝒟^N_x

ℕ𝒟^N_y    = 𝕄𝒟^N_y    = 𝒟^N_y   ,

so that the implementation can be greatly simplified.

These operators are presented below

Convergence studies of these operators are included in the Supplementary Material <ref>supp_convergence_tests.

\begin{align}
\label{SBP_matrices_1D_intertwined_A}
\mathcal A^M = \left[
\arraycolsep=3.2pt%\def\arraystretch{0.75}
\begin{array}{r r r r r r r r r r r r r r r r r}
\nicefrac{13}{12} &                 & \multicolumn{1}{r|}{} \\
& \nicefrac{7}{8} & \multicolumn{1}{r|}{} \\
&                 & \multicolumn{1}{r|}{ \nicefrac{25}{24} }  \\ \cline{1-3} \\[-1.75ex]
&  			    &    			& 		   1 &			  & \\
& 			    &    			& 			 & 			1 & \\[-0.75ex]
&  			    &    			& 			 &   		  & \ddots
\end{array}
\right],
\mathcal A^N = \left[
\arraycolsep=3.2pt%\def\arraystretch{0.75}
\begin{array}{r r r r r r r r r r r r r r r r r}
\nicefrac{7}{18} & 				   & 		   		& \multicolumn{1}{r|}{} \\
& \nicefrac{9}{8} &     			& \multicolumn{1}{r|}{} \\
&	               & 1  			& \multicolumn{1}{r|}{} \\
&    			   &    			& \multicolumn{1}{r|}{ \nicefrac{71}{72} }    \\ \cline{1-4} \\[-1.75ex]
&   			   &    			&    			&  1  			&    		  \\
&    			   &  				&    			&    			&  1  		& \\[-0.75ex]
&  			   &    			&   		    &    			&      		& \ddots
\end{array}
\right];
\allowdisplaybreaks
\\
\label{SBP_matrices_1D_intertwined_D}
\mathcal D^N = \left[
\arraycolsep=3.2pt%\def\arraystretch{0.75}
\begin{array}{r r r r r r r r r r r r r r r r r}
-\nicefrac{79}{78} &  \nicefrac{27}{26} & -\nicefrac{1}{26}  &  \nicefrac{1}{78}  & \multicolumn{1}{r|}{0}      \\
\nicefrac{2}{21}  & -\nicefrac{9}{7}   &  \nicefrac{9}{7}   & -\nicefrac{2}{21}  & \multicolumn{1}{r|}{0}      \\
\nicefrac{1}{75}  &  0			        & -\nicefrac{27}{25} &  \nicefrac{83}{75} & \multicolumn{1}{r|}{ -\nicefrac{1}{25} } & \\ \cline{1-5}
&   					&  \nicefrac{1}{24}  & -\nicefrac{9}{8}   & \nicefrac{9}{8}   						 & -\nicefrac{1}{24} 		  & \\
&     				&    				 &  \nicefrac{1}{24}  & -\nicefrac{9}{8} 	  					 &  \nicefrac{9}{8}  		  & -\nicefrac{1}{24} \\[-0.5ex]
& 					& 					 & 					  & \multicolumn{1}{l}{\ddots} 				 & \multicolumn{1}{l}{\ddots} & \multicolumn{1}{l}{\ddots} & \multicolumn{1}{l}{\ddots}
\end{array}
\right],
\mathcal D^M = \left[
\arraycolsep=3.2pt%\def\arraystretch{0.75}
\begin{array}{r r r r r r r r r r r r r r r r r}
\nicefrac{79}{28} & -\nicefrac{3}{14} & -\nicefrac{1}{28}  &  0    			    & \multicolumn{1}{r|}{0}   \\
-1       		  &  1    			  &  0    			   &  0    			    & \multicolumn{1}{r|}{0}   \\
\nicefrac{1}{24} & -\nicefrac{9}{8}  &  \nicefrac{9}{8}   & -\nicefrac{1}{24}  & \multicolumn{1}{r|}{0} & \\
-\nicefrac{1}{71} &  \nicefrac{6}{71} & -\nicefrac{83}{71} &  \nicefrac{81}{71} & \multicolumn{1}{r|}{ -\nicefrac{3}{71} } & \\ \cline{1-5}
&    				  &  \nicefrac{1}{24}  & -\nicefrac{9}{8}   &  \nicefrac{9}{8}   		 			   & -\nicefrac{1}{24} 		    & \\
&     			  &    				   &  \nicefrac{1}{24}  & -\nicefrac{9}{8}   		 			   &  \nicefrac{9}{8}  		    & -\nicefrac{1}{24} \\[-0.5ex]
& 				  & 				   & 					& \multicolumn{1}{l}{\ddots} 			   & \multicolumn{1}{l}{\ddots} & \multicolumn{1}{l}{\ddots} & \multicolumn{1}{l}{\ddots}
\end{array}
\right].
\end{align}

Comparing with (<ref>), the only difference is in the first row of $\mathcal D^M$, which requires some explanation.

When interpreted as the $\mathbb M$ operators, these operators can be used directly since they are already truncated; see Remark <ref>.

Taking the top surface for example, the first row of $\mathcal D^M$ applies on $\sigma_{yy}$ to update $v_y$ on the surface.

The full stencil before truncation is actually $\small[-\nicefrac{18}{7} \,\, \nicefrac{79}{28} \,\, -\nicefrac{3}{14} \,\, -\nicefrac{1}{28} \,\, 0 \,\, 0 \,\, \cdots ]$ with an extra $-\nicefrac{18}{7}$ on the left, which corresponds to $\sigma_{yy}$ on the top surface. It can be omitted because of the free surface boundary condition, i.e., $\sigma_{yy} = 0$, on the top surface. This stencil is second order accurate, same as those for the other grid points within the boundary region.

When interpreted as $\mathbb N$ operators, we need to reset the rows in $D^M$ and columns in $D^N$ that corresponding to the free surface, as described in section <ref>.

With this resetting in mind, the stencil mentioned above (i.e., the first row of $\mathcal D^M$) is not asked to provide accurate update for $\sigma_{xy}$ on the surface.

We note here that the above operators are not the unique choice.

Specifically, the additional requirement (<ref>) is optional and included here for the convenience of implementation.

One can potentially use the extra degrees of freedom for other properties that one may desire.

§.§ Numerical example

Numerical experiments similar to that presented for the 2D acoustic case in section <ref> is conducted for the 2D elastic case.

The specification is mostly the same with a few exceptions.

First, the grid points per minimal wavelength used here are increased to 20, 60, and 100.

This follows the wisdom that elastic wave equation admits a surface wave mode that requires finer grid spacing to resolve than for the body waves; see [21].

(However, experiments using 10, 30, 50 ppw reveal the same pattern.)

Moreover, the source is applied on $\sigma_{xx}$ and $\sigma_{yy}$ at one third of the simulation domain horizontally and half grid spacing below the surface vertically.

Time history of the solution variable $\Sigma_{xy}$ is recorded at two thirds of the simulation domain horizontally and at 0, 1 and 2 grid spacings below the surface vertically, which are shown in Figures <ref> and <ref>, for the strong and weak case, respectively.

The time step length is specified as 2e-4 s for all simulations, and the simulation time is specified as 2 s.

To conserve space, time histories of the other solution components are omitted here and included in the Supplementary Material <ref>supp_additional_figures_2D_elastic.

width=1, font=small,labelfont=small

width=1, font=small,labelfont=small

Time history of $\Sigma_{xy}$ at the surface.

width=1, font=small,labelfont=small

Time history of $\Sigma_{xy}$ at 1 grid points below the surface for the case ppw = 20.

width=1, font=small,labelfont=small

Time history of $\Sigma_{xy}$ at 2 grid points below the surface for the case ppw = 20.

Time histories of the 2D elastic experiments. The free surface boundary condition is imposed strongly.

width=1, font=small,labelfont=small

width=1, font=small,labelfont=small

Time history of $\Sigma_{xy}$ at the surface.

width=1, font=small,labelfont=small

Time history of $\Sigma_{xy}$ at 1 grid points below the surface for the case ppw = 20.

width=1, font=small,labelfont=small

Time history of $\Sigma_{xy}$ at 2 grid points below the surface for the case ppw = 20.

Time histories of the 2D elastic experiments. The free surface boundary condition is imposed weakly.

From the above figures, we observe that both the strong and weak approaches deliver consistent and accurate results, with the strong approach offers only minor benefit at or near the free surface for the $\Sigma_{xy}$ component.

Their performance for the other solution components are even less distinguishable.

A “stress test” of both approaches is conducted, where the source is applied on $V_y$ (i.e., a directional source) on the top surface and half grid spacing away from the left surface while the receivers are placed at 5 grid points (for the case ppw = 20) away from the source.

In this case, the strong approach demonstrates noticeable improvement in the $\Sigma_{xy}$ component, as can be observed by comparing Figures <ref> and <ref>.

There are little meaningful differences in the other recorded solution components, which are omitted here to conserve space but included in the Supplementary Material <ref>supp_additional_figures_2D_elastic_stress_test_Vy.

width=1, font=small,labelfont=small

width=1, font=small,labelfont=small

Time history of $\Sigma_{xy}$ at the surface.

width=1, font=small,labelfont=small

Time history of $\Sigma_{xy}$ at 1 grid points below the surface for the case ppw = 20.

width=1, font=small,labelfont=small

Time history of $\Sigma_{xy}$ at 2 grid points below the surface for the case ppw = 20.

Time histories of the 2D elastic experiments. The free surface boundary condition is imposed strongly.

width=1, font=small,labelfont=small

width=1, font=small,labelfont=small

Time history of $\Sigma_{xy}$ at the surface.

width=1, font=small,labelfont=small

Time history of $\Sigma_{xy}$ at 1 grid points below the surface for the case ppw = 20.

width=1, font=small,labelfont=small

Time history of $\Sigma_{xy}$ at 2 grid points below the surface for the case ppw = 20.

Time histories of the 2D elastic experiments. The free surface boundary condition is imposed weakly.

Another “stress test” is included in the Supplementary Material <ref>supp_additional_figures_2D_elastic_stress_test_Sxy, where the source is applied on $\Sigma_{xy}$ at the corner, more specifically, 1 grid spacing away from the top and left surfaces.

In this case, the strong approach performs markedly better than the weak approach for all five solution components.

Moreover, to accompany the above numerical experiments, convergence test results using a smooth manufactured solution for the elastic wave system (<ref>) is included in Supplementary Material <ref>, where the strong approach exhibits smaller error and higher convergence rates.

Finally, numerical experiments reveal that the CFL restriction for the strong and weak cases are $\nicefrac{6}{7}$ and 0.849, respectively. In other words, the weak case has a slight penalty on the time step size allowed while the strong case has none.

We note here that the results presented in this work is oriented towards a particular simulation setting (although appearing broadly in practice for various important wave-related applications), and should not be interpreted beyond its scope.

Indeed, some authors have argued for the benefit of weakly imposed boundary conditions over the strongly imposed ones; see [22] for an example in the fluid setting.

§ CONCLUSIONS

We have considered the imposition of the free surface boundary condition for both the acoustic and elastic wave simulations in the framework of summation-by-parts finite difference discretizations.

Within this framework, the standard approach of imposing boundary conditions, i.e., weakly through appending simultaneous approximation terms, despite its flexibility, can have issues in certain practical situations, such as when a point source is placed near the surface, inducing abrupt changes in the wave field, which have been demonstrated in this work through numerical experiments.

As an alternative, we have demonstrated the possibility of imposing the free surface boundary condition strongly, i.e., by incorporating it into the summation-by-parts finite difference operators.

While easy to achieve for the acoustic case, this approach requires specially designed grid layout for the elastic case and operators that need to satisfy additional requirements, as revealed by the discrete energy analysis.

One possible set of such operators is presented in this work, which also takes into account the ease of implementation in its design.

Numerical experiments reveal that the strong approach can have significant advantages over the weak approach, in terms of both accuracy and time step size restrictions.

Such advantages is particularly pronounced in the acoustic case.

The strong approach has its obvious disadvantage, namely, the boundary condition needs to be built into the operators and, therefore, renders the approach less flexible and less generally applicable.

However, in the particular case considered, i.e., free surface boundary condition for wave propagation, which has broad applications in geoscience and civil engineering, the additional effort is relatively small for the elastic case and trivial for the acoustic case.

For this reason, we recommend practitioners consider this option

in their practical settings.

§ CFL RESTRICTION FOR THE INTERIOR STENCIL AND STAGGERED LEAPFROG TIME INTEGRATION SCHEME ON UNBOUNDED DOMAIN

In this appendix, we consider the prototype 1D wave system (<ref>) with unit density and wave-speed on unbounded domain and derive the maximal Courant number associated with the standard fourth-order staggered grid spatial stencil
\nicefrac{ [\nicefrac{1}{24}, \enskip -\nicefrac{9}{8}, \enskip \nicefrac{9}{8}, \enskip -\nicefrac{1}{24}] }{ \Delta x }
and the staggered leapfrog time integration scheme via von Neumann analysis.

Figure <ref> illustrates the staggered spatial grids and the action of the stencil.

width=1, font=small,labelfont=small

Illustration of the staggered spatial grids and the action of the stencil.

In the following, we assume that variable $\sigma$ is discretized on the $N$-grid (which is associated with integer indices) and updated at integer time steps, while variable $v$ is discretized on the $M$-grid (which is associated with half-integer indices) and updated at half-integer time steps.

Considering the approximation of $\sigma$ at the $i$th $N$-grid point, i.e., $\Sigma_i$, when updated from time step $n-1$ to $n$, the following formula is used:

\begin{equation}
\label{update_nm1_to_n}
\Sigma_{i}^{(n)} - \Sigma_{i}^{(n-1)} = \left( \tfrac{1}{24} V_{i-\frac{3}{2}}^{(n-\frac{1}{2})} - \tfrac{9}{8} V_{i-\frac{1}{2}}^{(n-\frac{1}{2})} + \tfrac{9}{8} V_{i+\frac{1}{2}}^{(n-\frac{1}{2})} - \tfrac{1}{24} V_{i+\frac{3}{2}}^{(n-\frac{1}{2})} \right) \tfrac{ \Delta t }{ \Delta x }.
\end{equation}

Similarly, when $\Sigma_i$ is updated from time step $n$ to $n+1$, the following formula is used:

\begin{equation}
\label{update_n_to_np1}
\Sigma_{i}^{(n+1)} - \Sigma_{i}^{(n)} = \left( \tfrac{1}{24} V_{i-\frac{3}{2}}^{(n+\frac{1}{2})} - \tfrac{9}{8} V_{i-\frac{1}{2}}^{(n+\frac{1}{2})} + \tfrac{9}{8} V_{i+\frac{1}{2}}^{(n+\frac{1}{2})} - \tfrac{1}{24} V_{i+\frac{3}{2}}^{(n+\frac{1}{2})} \right) \tfrac{ \Delta t }{ \Delta x }.
\end{equation}

Subtracting (<ref>) from (<ref>) and substituting the update formulas for $V_{i-\frac{3}{2}}$, $V_{i-\frac{1}{2}}$, $V_{i+\frac{1}{2}}$, and $V_{i+\frac{3}{2}}$ from time step $n-\frac{1}{2}$ to $n+\frac{1}{2}$, we arrive at the following relation:

\begin{equation}
\label{update_relation_P}
\Sigma_i^{(n+1)} - 2 \Sigma_i^{(n)} + \Sigma_i^{(n-1)} = \left(
\tfrac{1}{576} \Sigma_{i-3}^{(n)}
- \tfrac{3}{32} \Sigma_{i-2}^{(n)}
+ \tfrac{87}{64} \Sigma_{i-1}^{(n)}
- \tfrac{365}{144} \Sigma_{i}^{(n)}
+ \tfrac{87}{64} \Sigma_{i+1}^{(n)}
- \tfrac{3}{32} \Sigma_{i+2}^{(n)}
+ \tfrac{1}{576} \Sigma_{i+3}^{(n)} \right)
\left( \tfrac{ \Delta t }{ \Delta x } \right)^2.
\end{equation}

Considering a particular wave mode $p = e^{\alpha t} e^{\imath k x}$, where $\imath$ denotes the imaginary unit, we have, after simplification, the following relation:

\begin{equation}
\label{error_relation}
e^{\alpha \Delta t} - 2 + e^{- \alpha \Delta t} =
\left(
\tfrac{1}{72} \cos(k \Delta x)^3
-\tfrac{3}{8} \cos(k \Delta x)^2
+\tfrac{65}{24} \cos(k \Delta x)
-\tfrac{169}{72} \right) \left( \tfrac{ \Delta t }{ \Delta x } \right)^2.
\end{equation}

Using $\beta$ to denote $\left( \tfrac{1}{72} \cos(k \Delta x)^3 - \tfrac{3}{8} \cos(k \Delta x)^2 + \tfrac{65}{24} \cos(k \Delta x) - \tfrac{169}{72} \right)$ for brevity, we have that $\beta \in \left[-\tfrac{49}{9},0\right]$.

Further, denoting $\tfrac{ \Delta t }{ \Delta x }$ as $C$, we have $C>0$ for sensible discretizations.

Finally, we denote $e^{\alpha \Delta t}$ as $G$, which is the growth factor between two neighboring time steps for
the wave mode

under consideration.

With these notations, we now have the following quadratic equation concerning $G$:

\begin{equation}
\label{quadratic_equation_G}
G^2 - \left(2 + \beta C^2 \right) G + 1 = 0\, ,
\end{equation}

which has two roots:

\begin{equation}
\label{roots_quadratic_equation_G}
G_1 = 1 + \tfrac{\beta C^2}{2} - \sqrt{ \left( 1 + \tfrac{\beta C^2}{2} \right)^2 - 1}
\text{\quad\enskip and \quad\enskip}
G_2 = 1 + \tfrac{\beta C^2}{2} + \sqrt{ \left( 1 + \tfrac{\beta C^2}{2} \right)^2 - 1}\, .
\end{equation}

Supposing that $\left( 1 + \tfrac{\beta C^2}{2} \right)^2 - 1 > 0$, we have that either $G_1\! < \!-1$, or $G_2\! > \!1$, i.e., this particular mode can grow out of bounds.

In the opposite case $\left( 1 + \tfrac{\beta C^2}{2} \right)^2 - 1 \leq 0$, we have that $\left| G_1 \right| = \left| G_2 \right| = 1$,

i.e., this particular mode is stable.

For the overall simulation to be stable, $C$ needs to be chosen such that $\left( 1 + \tfrac{\beta C^2}{2} \right)^2 - 1 \leq 0$ holds for all $\beta \in \left[-\tfrac{49}{9},0\right]$.

This constraint is at its strictest when $\beta = -\tfrac{49}{9}$, which leads to $C \leq \nicefrac{6}{7}$, as stated in (<ref>).

[1]
K. Aki and P. Richards.
Quantitative Seismology.
University Science Books, 2002.

[2]
J. Tromp.
Seismic wavefield imaging of Earth’s interior across scales.
Nature Reviews Earth & Environment, 1(1):40–53, 2020.

[3]
L. Kallivokas, A. Fathi, S. Kucukcoban, K. Stokoe II, J. Bielak, and
O. Ghattas.
Site characterization using full waveform inversion.
Soil Dynamics and Earthquake Engineering, 47:62–82, 2013.

[4]
B. W. Drinkwater and P. D. Wilcox.
Ultrasonic arrays for non-destructive evaluation: A review.
NDT & E International, 39(7):525–541, 2006.

[5]
H.-O. Kreiss and G. Scherer.
Finite element and finite difference methods for hyperbolic
partial differential equations. In: Mathematical aspects of finite elements
in partial differential equations.
Academic Press, 1974.

[6]
S. Nilsson, N. A. Petersson, B. Sjögreen, and H.-O. Kreiss.
Stable difference approximations for the elastic wave equation in
second order formulation.
SIAM Journal on Numerical Analysis, 45(5):1902–1936, 2007.

[7]
D. Appelö and N. A. Petersson.
A stable finite difference method for the elastic wave equation on
complex geometries with free surfaces.
Communications in Computational Physics, 5(1):84–107, 2009.

[8]
B. Sjögreen and N. A. Petersson.
A fourth order accurate finite difference scheme for the elastic wave
equation in second order formulation.
Journal of Scientific Computing, 52(1):17–48, 2012.

[9]
S. Wang and G. Kreiss.
Convergence of summation-by-parts finite difference methods for the
wave equation.
Journal of Scientific Computing, 71(1):219–245, 2017.

[10]
S. Wang and N. A. Petersson.
Fourth order finite difference methods for the wave equation with
mesh refinement interfaces.
SIAM Journal on Scientific Computing, 41(5):A3246–A3275, 2019.

[11]
O. O'Reilly, T. Lundquist, E. M. Dunham, and J. Nordström.
Energy stable and high-order-accurate finite difference methods on
staggered grids.
Journal of Computational Physics, 346:572–589, 2017.

[12]
L. Gao, D. C. D. R. Fernández, M. Carpenter, and D. Keyes.
SBP-SAT finite difference discretization of acoustic wave
equations on staggered block-wise uniform grids.
Journal of Computational and Applied Mathematics, 348:421–444,

[13]
L. Gao and D. Keyes.
Combining finite element and finite difference methods for isotropic
elastic wave simulations in an energy-conserving manner.
Journal of Computational Physics, 378:665–685, 2019.

[14]
L. Gao, O. Ghattas, and D. Keyes.
Nonuniform 3D finite difference elastic wave simulation on
staggered grids.
Geophysics, 87(4):1–79, 2022.

[15]
O. O’Reilly, T.-Y. Yeh, K. B. Olsen, Z. Hu, A. Breuer, D. Roten, and C. A.
A high-order finite-difference method on staggered curvilinear grids
for seismic wave propagation applications with topography.
Bulletin of the Seismological Society of America, 112(1):3–22,

[16]
M. H. Carpenter, D. Gottlieb, and S. Abarbanel.
Time-stable boundary conditions for finite-difference schemes solving
hyperbolic systems: methodology and application to high-order compact
Journal of Computational Physics, 111(2):220–236, 1994.

[17]
A. R. Levander.
Fourth-order finite-difference P-SV seismograms.
Geophysics, 53(11):1425–1436, 1988.

[18]
B. Fornberg and M. Ghrist.
Spatial finite difference approximations for wave-type equations.
SIAM Journal on Numerical Analysis, 37(1):105–130, 1999.

[19]
M. Ghrist, B. Fornberg, and T. Driscoll.
Staggered time integrators for wave equations.
SIAM Journal on Numerical Analysis, 38(3):718–741, 2000.

[20]
L. Gao and D. Keyes.
Explicit coupling of acoustic and elastic wave propagation in
finite-difference simulations.
Geophysics, 85(5):T293–T308, 2020.

[21]
H.-O. Kreiss and N. A. Petersson.
Boundary estimates for the elastic wave equation in almost
incompressible materials.
SIAM Journal on Numerical Analysis, 50(3):1556–1580, 2012.

[22]
Y. Bazilevs, C. Michler, V. Calo, and T. Hughes.
Weak Dirichlet boundary conditions for wall-bounded turbulent
Computer Methods in Applied Mechanics and Engineering,
196(49-52):4853–4862, 2007.

Supplementary Materials

§ ADDITIONAL 1D EXPERIMENT WITH FINER GRID RESOLUTION

An additional experiment is presented here to supplement the discussion in section <ref>, Figure <ref> in particular.

The experimental setup is identical to that has been illustrated in Figure <ref>, except that the ppw associated with the four grid lines, from top to bottom, are 40, 80, 160, and 320, respectively.

In other words, the grid spacing is a quarter of that used in the experiment from section <ref>.

The time step length is specified as 5e-5 s for all four simulations, and the simulation time is specified as 6 s, which amounts to 120000 time steps.

Time histories of the solution variable $\Sigma$ at the three receiver locations are displayed in Figure <ref>, from which we observe that the violation of the weakly imposed free surface boundary condition persists when the source is placed too close to the free surface, despite the refined grid resolution.

width=1, font=small,labelfont=small

width=1, font=small,labelfont=small

Time histories of $\Sigma$ at $R_0$.

width=1, font=small,labelfont=small

Time histories of $\Sigma$ at $R_1$.

width=1, font=small,labelfont=small

Time histories of $\Sigma$ at $R_2$.

Time histories of the solution variable $\Sigma$ at the three receiver locations depicted in Figure <ref> with the four grid lines now representing the cases of 40, 80, 160, 320 ppw, respectively.

§ ADDITIONAL 1D EXPERIMENT WITH THE SOURCE PLACED FAR AWAY FROM THE FREE SURFACE

An additional experiment is presented here to supplement the discussion in section <ref>, Figure <ref> in particular.

The experimental setup is identical to that has been illustrated in Figure <ref>, except that the source is placed at the middle of the simulation interval (10 minimal wavelength away from both boundaries, which correspond to 100, 200, 400, and 800 grid points for the four simulations, respectively).

Time histories of the solution variable $\Sigma$ at the three receiver locations are displayed in Figure <ref>, from which we observe that there is no severe violation of the free surface boundary condition

and that the four simulation results agree much better (cf. Figure <ref>).

width=1, font=small,labelfont=small

width=1, font=small,labelfont=small

Time histories of $\Sigma$ at $R_0$.

width=1, font=small,labelfont=small

Time histories of $\Sigma$ at $R_1$.

width=1, font=small,labelfont=small

Time histories of $\Sigma$ at $R_2$.

Time histories of the solution variable $\Sigma$ at the three receiver locations depicted in Figure <ref> with the source now placed far ($10$ times ppw grid points) away from the free surface.

§ CROSS COMPARISON

Simulation results from three different discretizations are compared against each other to validate the implementation, which includes the finite difference discretization described in section <ref> with weakly imposed free surface boundary condition, the finite difference discretization described in section <ref> with strongly imposed free surface boundary condition, and a finite element discretization.

The finite element discretization is based on a different formulation of the wave system from (<ref>).

Essentially, we reduce (<ref>) to a single equation by substituting (<ref>) into the time derivative of (<ref>), leading to an equation that involves only second order temporal and spatial derivatives of $\sigma$.

Starting with a different formulation gives better confidence in the validation.

In space, the finite element simulation uses second-order Lagrange polynomials; in time, it uses second-order central difference.

width=1, font=small,labelfont=small

width=1, font=small,labelfont=small

Time histories of $\Sigma$ at $R_0$.

width=1, font=small,labelfont=small

Time histories of $\Sigma$ at $R_1$.

width=1, font=small,labelfont=small

Time histories of $\Sigma$ at $R_2$.

Time histories of the solution variable $\Sigma$ at the three receiver locations depicted in Figure <ref> from three different discretizations, including the finite difference discretization described in section <ref> (denoted as weak), the finite difference discretization described in section <ref> (denoted as strong), and a finite element discretization (denoted as FEM). All three simulations use $\text{ppw}=80$.

For clarity, we only compare the simulation results with $\text{ppw}=80$ for all three discretizations. The grid spacing and time step length are the same for all three simulations.

Time histories of the solution variable $\Sigma$ at the three receiver locations depicted in Figure <ref> are displayed in Figure <ref>, from where we observe that results from all three simulations agree well with each other, hence providing confidence in our implementations used for the experiments.

§ SPECTRAL RADII

In the following, we consider the 1D wave system (<ref>) defined on interval $[0,1]$ with unit density and wave-speed; spectral radii of the semi-discretized systems described in sections <ref> and <ref> are evaluated numerically to confirm the claims on CFL restrictions made therein.

Both semi-discretized systems are represented using the following generic notations:

\begin{equation}
\label{generalized_eigenvalue_system_large}
\frac{d }{d t}
\left[
\begin{array}{c c}
\mathcal A^M & \boldsymbol 0 \\
\boldsymbol 0 & \mathcal A^N
\end{array}
\right]
\left[
\begin{array}{c}
V \\
\Sigma
\end{array}
\right]
\enskip = \enskip
\left[
\begin{array}{c c}
\boldsymbol 0 & \mathcal A^M \widetilde{\mathcal D}^N \\
\mathcal A^N \widetilde{\mathcal D}^M & \boldsymbol 0
\end{array}
\right]
\left[
\begin{array}{c}
V \\
\Sigma
\end{array}
\right],
\end{equation}

where $\widetilde{\mathcal D}^M$ and $\widetilde{\mathcal D}^N$ are the modified finite difference operators, either by appending SATs to (<ref>) as described in section <ref> (the weak case), or by setting the specific rows and columns to zero as described in section <ref> (the strong case).

In both cases, $\widetilde{\mathcal D}^M$ and $\widetilde{\mathcal D}^N$ satisfy the following relation:

\begin{equation}
\label{relation_modified_operators}
\mathcal A^N \widetilde{\mathcal D}^M \ + \ \left( \mathcal A^M \widetilde{\mathcal D}^N \right)^T \ = \ \boldsymbol 0
\end{equation}

by construction.

The generalized eigenvalues associated with the above system, together with the chosen time integration scheme, determine the stability requirement on time step length.

We note here that since the two matrices from (<ref>) are diagonal and skew-symmetric, all generalized eigenvalues are on the imaginary axis.

Moreover, because of the skew-symmetric nature, the generalized eigenvalues calculated based on (<ref>) tend to be more reliable than those from the system where the norm matrices $\mathcal A^M$ and $\mathcal A^N$ are omitted on both sides.

A related note is that system (<ref>) can be further reduced to

\begin{equation}
\label{generalized_eigenvalue_system_small}
\mathcal A^M \frac{d^2 V}{d t^2}
\enskip = \enskip
- \left( \widetilde{\mathcal D}^M \right)^T \!\! \mathcal A^N \widetilde{\mathcal D}^M V \, ,
\end{equation}

by differentiating with respect to $t$, eliminating $\Sigma$ by substitution, and invoking (<ref>).

It can be shown that the generalized eigenvalues from (<ref>) are the square roots of the generalized eigenvalues from (<ref>).

Since (<ref>) is smaller in size and involves only symmetric matrices, the associated generalized eigenvalue problem is easier to solve numerically with more reliable outcome.

The spectral radii (i.e., maximal magnitude of the generalized eigenvalues) are shown in Table <ref> for $\Delta x = \nicefrac{1}{40}$,  $\Delta x = \nicefrac{1}{160}$, and $\Delta x = \nicefrac{1}{640}$.

All calculated generalized eigenvalues associated with (<ref>) are indeed on the imaginary axis.

For staggered leapfrog scheme, its associated stability region is $[-2\imath, 2\imath]$; see [19].

In this case, the spectral radii dictates the restriction on time step length.

The spectral radii scaled by $C_\text{max} \cdot \Delta x$ (i.e., the maximally allowed $\Delta t$ according to $C_\text{max}$) are presented in Table <ref>, which need to be less than or equal to 2 for stable simulations with the staggered leapfrog scheme.

The values in Table <ref> are right beneath 2, which confirms the claims made in sections <ref> and <ref> on time step restrictions.

width=1, font=small,labelfont=small

1|c|spectral radii      				    $\Delta x = \nicefrac{1}{40}$     $\Delta x = \nicefrac{1}{160}$     $\Delta x = \nicefrac{1}{640}$

Extrapolating (weak)       $125.871385897805$     $503.485543591221$     $2013.942174364883$

Extrapolating (strong)     $93.256016056512$      $373.311280516844$     $1493.327619757113$

Periodic 			       $93.333333333333$      $373.333333333333$     $1493.333333333333$

Spectral radii corresponding to the semi-discretized systems described in sections <ref> (weak) and <ref> (strong).

The spectral radii associated with periodic boundary condition are also presented as references, which amounts to applying the interior stencil on unbounded domain.

width=1, font=small,labelfont=small

1|c|spectral radii (scaled)           				    $\Delta x = \nicefrac{1}{40}$     $\Delta x = \nicefrac{1}{160}$     $\Delta x = \nicefrac{1}{640}$

Extrapolating (weak)       $1.999781643451$     $1.999781643451$     $1.999781643451$

Extrapolating (strong)     $1.998343201211$     $1.999881859912$     $1.999992347889$

Periodic 			       $2.000000000000$     $2.000000000000$     $2.000000000000$

Spectral radii presented in Table <ref> scaled by
$C_\text{max} \cdot \Delta x$, where $C_\text{max}$ is the numerically measured maximal Courant number. Such $C_\text{max}$ are $0.6355$, $\nicefrac{6}{7}$, and $\nicefrac{6}{7}$, as mentioned in (<ref>), (<ref>), and (<ref>) of section <ref>.

§ ACOUSTIC WAVE EQUATION AND THE SEMI-DISCRETIZED SYSTEM

Below, we present the acoustic wave equation and its semi-discretization as the background information omitted in the main text, section <ref>.

The acoustic wave equation can be expressed as

[left=]alignat = 2
ρ∂v_x/∂t     =   ∂σ/∂x;

ρ∂v_y/∂t     =   ∂σ/∂y;

β∂σ/∂t     =   ∂v_x/∂x + ∂v_y/∂y,

where $\rho$ and $\beta$ are the physical parameters, representing density and compressibility, respectively, $v_x$ and $v_y$ are the particle velocities, $\sigma$ is the negative of pressure.

We note here that for consistency within this study, we use $\sigma$, i.e., the negative of pressure, instead of pressure itself, so that all wave equations under discussion will have the same sign in front of the spatial derivatives.

Moreover, $\beta = \nicefrac{1}{\rho c^2}$, where $c$ is wave-speed.

The semi-discretization of (<ref>) without concerning the boundary conditions can be expressed as

[left=]alignat = 2

𝒜^V_x ρ^V_x d V_x/d t     =  𝒜^V_x 𝒟^Σ_x   Σ  ;

𝒜^V_y ρ^V_y d V_y/d t     =  𝒜^V_y 𝒟^Σ_y   Σ  ;

𝒜^Σ β^Σ d Σ/d t     =  𝒜^Σ 𝒟^V_x_x V_x + 𝒜^Σ 𝒟^V_y_y V_y   .

Following the standard procedure, the 2D SBP operators appearing in (<ref>) are constructed from their 1D counterparts via tensor product as

\mathcal A^{P}   \! = \mathcal A^N_x \otimes \mathcal A^N_y, \enskip
\mathcal A^{V_x} \! = \mathcal A^M_x \otimes \mathcal A^N_y, \enskip
\mathcal A^{V_y} \! = \mathcal A^N_x \otimes \mathcal A^M_y;

\mathcal D_{x}^{P}   \! = \mathcal D^N_x \otimes \mathcal I^N_y, \enskip
\mathcal D_{x}^{V_x} \! = \mathcal D^M_x \otimes \mathcal I^N_y, \enskip
\mathcal D_{y}^{P}   \! = \mathcal I^N_x \otimes \mathcal D^N_y, \enskip
\mathcal D_{y}^{V_y} \! = \mathcal I^N_x \otimes \mathcal D^M_y.

For the numerical experiment considered in section <ref> of the main text, where the free surface boundary condition is associated with all boundaries, the above difference operators need to be modified to account for these boundary conditions.

For the strong case, these modifications involve setting the corresponding rows and columns in the 1D operators to zero as follows:

setting the first and last rows of $D_x^M$ to zero;

setting the first and last rows of $D_y^M$ to zero;

setting the first and last columns of $D_x^N$ to zero;

setting the first and last columns of $D_y^N$ to zero;

setting the first and last entries of $I_x^N$ to zero;

setting the first and last entries of $I_y^N$ to zero.

For the weak case, these modifications involve adjusting for the penalty terms as follows:

𝒟_x^P ⟵𝒟_x^P
[ ( ( 𝒜^M_x )^-1 ℰ^L_x ) ⊗ℐ^N_y ]
[ ( 𝒫^L_x )^T ⊗ℐ^N_y - 0 ]

[ ( ( 𝒜^M_x )^-1 ℰ^R_x ) ⊗ℐ^N_y ]
[ ( 𝒫^R_x )^T ⊗ℐ^N_y  - 0 ];

𝒟_y^P ⟵𝒟_y^P
[ ℐ^N_x ⊗( ( 𝒜^M_y )^-1 ℰ^L_y ) ]
[ ℐ^N_x ⊗( 𝒫^L_y )^T  - 0 ]

[ ℐ^N_x ⊗( ( 𝒜^M_y )^-1 ℰ^R_y ) ]
[ ℐ^N_x ⊗( 𝒫^R_y )^T  - 0 ],

where (<ref>) - (<ref>) accounts for the left, right, bottom, and top boundaries, respectively.

§ ADDITIONAL FIGURES FOR THE 2D ACOUSTIC EXPERIMENTS

Below, we present the time histories of the solution variables $V_x$ and $V_y$ for the 2D acoustic experiment from section <ref> of the main text.

For $V_x$, the strong and weak cases are presented in Figures <ref> and <ref>, respectively.

For $V_y$, the strong and weak cases are presented in Figures <ref> and <ref>, respectively.

Because of grid staggering, the receiver locations for $V_x$ is shifted rightwards for half grid spacing from those for $\Sigma$ in the main text; the receiver locations for $V_y$ is shifted downwards for half grid spacing from those for $\Sigma$ in the main text.

We note here that the first receiver of $V_x$ is on the surface, which is supposed to remain zero.

§ 1D CONVERGENCE TESTS

In the following, 1D convergence test results using a smooth manufactured solution are shown for the experiment settings presented in section <ref> (weak), section <ref> (strong) of the main text, and for the operators presented in (<ref>) of the main text applied on an $\mathbb M$ grid column of Figure <ref> of the main text.

These operators applied on an $\mathbb N$ grid column delivers the same results as the strong case corresponding to section <ref> since after resetting, the operators become the same, and are omitted below.

Specifically, the manufactured solution reads:

\begin{equation}
\label{smooth_manufactured_solution}
\begin{array}{rcl}
p & \!\!\! = \!\!\! & \displaystyle \phantom{-}  \sin(8 \pi x) \sin(8 \pi t); \\
v & \!\!\! = \!\!\! & \displaystyle - \cos(8 \pi x) \cos(8 \pi t).
\end{array}
\end{equation}

where $x \in [0,1]$.

For all simulations involved in the following tests, the time step length is fixed at 1e-6 s, while the total number of time steps is chosen as 666667, leading to a final time at approximately 2/3 s.

The error between the simulated solution and the manufactured solution at the final time is measured using the energy norm (i.e., the weighted $L_2$ norm with the norm matrices serving as the weights) and displayed below.

width=1, font=small,labelfont=small

10 ppw         20 ppw         40 ppw         80 ppw         160 ppw

error  	    2.5174e-02     2.1005e-03     1.6514e-04     1.4079e-05     1.2281e-06

rate        —            3.5832         3.6689         3.5521         3.5191

Error and rate for the weak case, corresponding to section <ref> of the main text.

width=1, font=small,labelfont=small

10 ppw         20 ppw         40 ppw         80 ppw         160 ppw

error       2.4014e-02     1.2913e-03     6.2780e-05     3.4777e-06     2.2958e-07

rate        —            4.2169         4.3624         4.1741         3.9210

Error and rate for the strong case, corresponding to section <ref> of the main text.

width=1, font=small,labelfont=small

10 ppw         20 ppw         40 ppw         80 ppw         160 ppw

error       1.0505e-02     5.5674e-04     3.2667e-05     2.5932e-06     2.1248e-07

rate        —            4.2379         4.0911         3.6550         3.6093

Error and rate for the strong case, corresponding to applying the operators presented in (<ref>) of the main text on an $\mathbb M$ grid column of Figure <ref> of the main text.

§ ADDITIONAL FIGURES FOR THE 2D ELASTIC EXPERIMENTS

Below, we present the time histories of the solution variables $\Sigma_{xx}$, $\Sigma_{yy}$, $V_x$, and $V_y$ for the 2D elastic experiment from section <ref> of the main text, corresponding to Figures <ref> and <ref> therein, where a compressional source (i.e., applied on $\Sigma_{xx}$ and $\Sigma_{yy}$) is considered.

For $S_{xx}$, the strong and weak cases are presented in Figures <ref> and <ref>, respectively.

For $S_{yy}$, the strong and weak cases are presented in Figures <ref> and <ref>, respectively.

For $V_x$, the strong and weak cases are presented in Figures <ref> and <ref>, respectively.

For $V_y$, the strong and weak cases are presented in Figures <ref> and <ref>, respectively.

§ ADDITIONAL FIGURES FOR THE 2D ELASTIC “STRESS TEST” WITH DIRECTIONAL SOURCE

Below, we present the time histories of the solution variables $\Sigma_{xx}$, $\Sigma_{yy}$, $V_x$, and $V_y$ for the 2D elastic experiment from section <ref> of the main text, corresponding to Figures <ref> and <ref> therein, where a directional source (i.e., applied on $V_y$) is considered and placed at the corner.

For $S_{xx}$, the strong and weak cases are presented in Figures <ref> and <ref>, respectively.

For $S_{yy}$, the strong and weak cases are presented in Figures <ref> and <ref>, respectively.

For $V_x$, the strong and weak cases are presented in Figures <ref> and <ref>, respectively.

For $V_y$, the strong and weak cases are presented in Figures <ref> and <ref>, respectively.

§ ANOTHER 2D ELASTIC “STRESS TEST” WITH SOURCE APPLIED ON $\SIGMA_{XY}$

Below, we present the time histories of another “stress test” for the strong and weak approaches.

Specifically, the source is applied on $\Sigma_{xy}$ (although not common, this can still happen in practice, e.g., when complex source pattern is being inverted) and placed at 1 grid spacing away from the top and left surfaces, i.e., the first interior $\Sigma_{xy}$ point at the top left corner.

The receiver is placed at 5 grid spacing away from the left surface and at 0, 1, and 2 grid spacing at depth for $\Sigma_{xy}$ and $V_y$, with an additional half grid downward shift for the other components due to grid staggering.

The recorded time histories are presented in Figures <ref> - <ref> , in the order of $\Sigma_{xx}$, $\Sigma_{yy}$, $\Sigma_{xy}$, $V_x$, and $V_y$, from which we observe that the strong approach performs markedly better than the weak case for all five solution components.

§ 2D CONVERGENCE TESTS

In the following,

2D convergence test results using a smooth manufactured solution for the elastic wave system (see (<ref>) of the main text) are shown for the strong (Table <ref>) and weak (Table <ref>) approaches to accompany the discussion and comparison made in section <ref> of the main text.

Specifically, the physical domain is a square with unit length on both $x$- and $y$-directions.

The manufactured solution reads:

[left=]alignat = 4

σ_xx =                2 k sin(w t) sin(k x) sin(k y);

σ_xy =                0;

σ_yy =      -         2 k sin(w t) sin(k y) sin(k x);

v_x 		=      -         w cos(w t) sin(k y) cos(k x);

v_y 		=                w cos(w t) cos(k y) sin(k x),

where $k$ is set to $2\pi$ and $w$ is set to $\sqrt{2}k$.

For all simulations involved in the following tests, the time step length is fixed at 1e-6 s, while the total number of time steps is chosen as 666667, leading to a final time at approximately 2/3 s.

The error between the simulated solution and the manufactured solution at the final time is measured using the energy norm, which is derived from (<ref>) of the main text.

width=1, font=small,labelfont=small

10 ppw         20 ppw         40 ppw         80 ppw         160 ppw

error  	    1.7389e-01     1.1894e-02     7.9771e-04     6.3477e-05     5.4416e-06

rate        —            3.8699         3.8982         3.6515         3.5441

Error and rate for the strong approach

width=1, font=small,labelfont=small

10 ppw         20 ppw         40 ppw         80 ppw         160 ppw

error  	    2.8617e-01     3.3790e-02     5.1094e-03     8.6026e-04     1.5148e-04

rate        —            3.0822         2.7253         2.5703         2.5057

Error and rate for the weak approach.