Modeling a water layer (Vs=0) above a solid with `deepwave.elastic`

[z981201]

Hi Deepwave authors,

I’m using deepwave.elastic and trying to model a fluid–solid configuration where the upper layer is fluid (i.e., Vs = 0) and the lower layer is an elastic solid (Vs > 0), similar to a Marmousi-II–type setup with a water layer over an elastic medium.

In practice, setting Vs = 0 in the top layer seems to cause numerical issues. Specifically, when I place receivers near the bottom of the fluid layer (close to the fluid–solid interface), I observe noticeable oscillations in the recorded data.

Could you please let me know whether deepwave.elastic is expected to support models with Vs = 0 regions (μ = 0) adjacent to elastic solids, and whether there are recommended settings or workarounds for stable fluid–solid simulations?

Thanks a lot for your help!

[ar4]

Hi @z981201

Thank you for your question. I am sorry to hear that you experienced poor results while using Deepwave.

May I suggest that you try using a higher number of grid cells per wavelength (smaller grid spacing) to see whether that helps? If you have not already tried it, using 8th order accuracy (accuracy=8) might also be worth testing. Please let me know whether either of these suggestions help.

[z981201]

Hi Dear Deepwave authors,

Thanks for the suggestions, and for taking the time to look into this.

I tried increasing the grid resolution, and I also looked for an 8th-order option. However, as far as I can tell, deepwave.elastic does not accept accuracy=8 (it seems to only support accuracy=2 or 4). Please correct me if I’m missing something.

To make the issue easier to reproduce, I put together a small minimal example with a constant model (Vs = 0 everywhere). Even with a relatively fine spacing (dx = 10 m) and freq = 10 Hz, I still observe clear oscillations/ringing in the recorded trace, especially in a later time window (shown in the final zoomed plot).

Minimal reproducible example:

import torch
import matplotlib.pyplot as plt

import deepwave
from deepwave import elastic


def my_ricker(freq_hz: float, nt: int, dt: float, peak_time: float):
    t = np.arange(nt, dtype=np.float32) * dt - peak_time
    r = np.pi * freq_hz * t
    w = (1.0 - 2.0 * r**2) * np.exp(-r**2)
    return w.astype(np.float32)

device = torch.device("cuda:0" if torch.cuda.is_available() else "cpu")
torch.manual_seed(0)

ny, nx = 350, 1000
dx = 10.0
dt = 0.5e-3
nt = 9000

pml_width = 100
pml_freq = 14

vp = torch.full((ny, nx), 1500.0, device=device, dtype=torch.float32)
vs = torch.zeros((ny, nx), device=device, dtype=torch.float32)  # Vs=0
rho = torch.full((ny, nx), 1000.0, device=device, dtype=torch.float32)

n_shots = 6
n_sources_per_shot = 1
d_source = 160
first_source = 7
source_depth = 2

n_receivers_per_shot = 200
d_receiver = 5
first_receiver = 2
receiver_depth = 40

freq = 10.0
peak_time = 1.0 / freq

source_locations = torch.zeros(n_shots, n_sources_per_shot, 2, dtype=torch.long, device=device)
source_locations[..., 0] = source_depth
source_locations[:, 0, 1] = torch.arange(n_shots, device=device) * d_source + first_source

receiver_locations = torch.zeros(n_shots, n_receivers_per_shot, 2, dtype=torch.long, device=device)
receiver_locations[..., 0] = receiver_depth
receiver_locations[:, :, 1] = (
    torch.arange(n_receivers_per_shot, device=device) * d_receiver + first_receiver
).repeat(n_shots, 1)

w = torch.from_numpy(my_ricker(freq, nt, dt, peak_time)).to(device)
source_amplitudes = w[None, None, :].repeat(n_shots, n_sources_per_shot, 1).contiguous()

out = elastic(
    *deepwave.common.vpvsrho_to_lambmubuoyancy(vp, vs, rho),
    dx, dt,
    source_amplitudes_x=source_amplitudes,
    source_amplitudes_y=source_amplitudes,
    source_locations_x=source_locations,
    source_locations_y=source_locations,
    receiver_locations_x=receiver_locations,
    receiver_locations_y=receiver_locations,
    accuracy=4,
    pml_width=pml_width,
    pml_freq=pml_freq,
    model_gradient_sampling_interval=15,
)
data0, data1 = out[-2], out[-1]

src_idx = 3
rec_idx = 100
trace = data1[src_idx, rec_idx].detach().cpu().numpy()
t = np.arange(nt) * dt

plt.figure(figsize=(10, 3))
plt.plot(t, trace)
plt.xlabel("Time (s)")
plt.ylabel("Amplitude")
plt.title(f"data1 trace (shot={src_idx}, rec={rec_idx})")

i0, i1 = 2000, 4000
plt.figure(figsize=(10, 3))
plt.plot(t[i0:i1], trace[i0:i1])
plt.xlabel("Time (s)")
plt.ylabel("Amplitude")
plt.title(f"Zoomed window [{i0}:{i1}] (shot={src_idx}, rec={rec_idx}) - oscillation/ringing")

plt.show()

Thank you again for your help.

[ar4]

Thank you for your code, and for producing a minimal example.

Your statement about the elastic propagator only supporting 2nd and 4th order accuracy might reveal the cause of the problem, as it indicates that you are using an older version of Deepwave. Would it be possible for you to update it to a new version? Better handling of vs=0 in the elastic propagator was one of the changes made over the past few months, so I hope that it might help.

When I run your code using the current version of Deepwave, I get the following:

I see some ringing (the following is the zoomed-in plot produced by your code):

The amplitude is quite small, however. Is the ringing that you experienced more severe, and the result with the current version would be acceptable, or is even this level of ringing a problem for your application?

[z981201]

Thanks again for testing my minimal example and for the suggestion to update Deepwave.

I’ve updated to the latest version. While the minimal example looks better, I’m still seeing an issue in my real application: I’m running EFWI on the Marmousi model, and I consistently observe small “micro-artifacts” / ringing-like noise in the inverted result, as shown below:

The corresponding 1D profile at x = 6 km is:

Do you have any recommended settings or best practices for EFWI with Deepwave? I’m happy to share a minimal reproducible EFWI script and additional plots if that would help.

Thanks again for your time and help.