Units and source conventions¶
A frequent source of confusion is the amplitude and units of the source and the recorded data. The conventions differ between the propagators because they solve different equations, so they are collected here. In all cases the source_amplitudes* arguments have shape [shot, source, time] and are added as a source term at each (internal) time step. The recorded receiver data has shape [shot, receiver, time] and contains the value of the corresponding field sampled at each receiver location.
For most inversion and imaging work the relative scaling between sources and receivers is what matters, and any overall constant cancels or is absorbed into a calibration; but if you need physically meaningful amplitudes, the conventions below apply.
Warning
The scalar source channel is named for the quantity it injects, and the two families inject different quantities. In the acoustic family (acoustic(), acoustic_born()) it is volume_rate, a volume-injection rate per unit volume with units of 1/time; in the elastic family (elastic(), elastic_vti(), elastic_born()) and in acoustic_vti() it is pressure_rate, with units of pressure/time. Viscoacoustic also uses pressure_rate, applied directly to its pressure equation. The volume-rate and pressure-rate quantities differ by a modulus: see Converting between propagators. Moving a wavelet from one family to the other unchanged rescales your source by that modulus (of order \(10^9\) in SI units for rock).
The letter p used to name both, which is why it is deprecated as a source name. On the released acoustic and elastic APIs, source_amplitudes_p and sources={"p": ...} still work with a deepwave.DeepwaveFutureWarning until 0.4.0. p is unchanged as a receiver channel, where it has always meant a pressure.
Calibration table¶
The quantity each source argument represents, derived from the injection code. \(\lambda\) and \(\mu\) are the Lamé parameters, \(K\) the bulk modulus, \(\rho\) the density, \(n\) the number of spatial dimensions, and \(\Delta V\) the product of the grid spacings (\(\Delta z\Delta x\) in 2-D, \(\Delta z\Delta y\Delta x\) in 3-D). Here \(h\) is the internal time step, out.info.internal_dt, which can be smaller than the user sample interval dt. SI units are given for concreteness; Deepwave itself is unit-agnostic (see Choosing units for ultrasound).
Propagator | Argument | Enters | Physical quantity | SI units | 2-D versus 3-D |
|---|---|---|---|---|---|
| the scalar field \(u\), pre-scaled by \(-v^2h^2\) | the right-hand side \(f\) of \(\nabla^2u-\ddot u/v^2=f\) | \([u]/\mathrm{m}^2\) (\(\mathrm{Pa}/\mathrm{m}^2\) when \(u\) is pressure) | Units unchanged. \(f\,\Delta V\) is the point-source strength in 3-D; in 2-D it is the strength per metre of a line source along the omitted axis. | |
| pressure \(p\), pre-scaled by \(h\,K\) | volume-injection rate per unit volume \(s\) in \(K^{-1}\dot p=-\nabla\cdot\mathbf{v}+s\) | \(\mathrm{s}^{-1}\) | Units unchanged. \(s\,\Delta V\) is the injected volume rate: \(\mathrm{m}^3/\mathrm{s}\) in 3-D, \(\mathrm{m}^3/\mathrm{s}\) per metre of the omitted axis in 2-D. | |
| pressure \(p\), pre-scaled by \(h\) | additive external pressure rate \(a\) in the pressure equation; no source is added to the relaxation memories | \(\mathrm{Pa/s}\) | Units unchanged. \(a\,\Delta V\) is an integrated pressure-rate source in 3-D and per omitted-axis length in 2-D. | |
| the matching velocity component, pre-scaled by \(h/\rho\) | body-force density \(\mathbf{f}\) in \(\rho\dot{\mathbf{v}}=-\nabla p+\mathbf{f}\) | \(\mathrm{N}/\mathrm{m}^3\) | Units unchanged. \(f\,\Delta V\) is the total force (\(\mathrm{N}\)) in 3-D and the force per metre of the line source (\(\mathrm{N}/\mathrm{m}\)) in 2-D. | |
| both normal stresses, pre-scaled by \(h\) after a sign flip: \(s_H=s_V=-\dot p\) | pressure rate \(\dot p\) (positive is compressional), exactly as for the elastic family | \(\mathrm{Pa}/\mathrm{s}\) | Synthesised by a wrapper, not by the kernel: it is bit-identical to passing | |
| \(\sigma_H\) / \(\sigma_V\), pre-scaled by \(h\) | stress rate \(s_H\), \(s_V\) (tension positive, so the opposite sign convention to | \(\mathrm{Pa}/\mathrm{s}\) | Units unchanged. \(s\,\Delta V\) is a moment rate, \(\mathrm{N\,m}/\mathrm{s}\) in 3-D and per metre in 2-D. | |
| the matching velocity component, pre-scaled by \(h/\rho\) | body-force density, as for | \(\mathrm{N}/\mathrm{m}^3\) | As for | |
| every normal stress, pre-scaled by \(-h\) | pressure rate \(\dot p\) (positive is compressional), i.e. an isotropic moment-rate source | \(\mathrm{Pa}/\mathrm{s}\) | Units unchanged. \(\dot p\,\Delta V\) is the isotropic moment rate, \(\mathrm{N\,m}/\mathrm{s}\) in 3-D and per metre in 2-D. | |
| the matching velocity component, pre-scaled by \(h/\rho\) | body-force density \(\mathbf{f}\) in \(\rho\dot{\mathbf{v}}=\nabla\cdot\mathbf{\sigma}+\mathbf{f}\) | \(\mathrm{N}/\mathrm{m}^3\) | Units unchanged. \(f\,\Delta V\) is the total force (\(\mathrm{N}\)) in 3-D and the force per metre of the line source (\(\mathrm{N}/\mathrm{m}\)) in 2-D. |
Sources are named by quantity; receivers are named by field. The two propagator families would otherwise use the same letter p for two quantities that differ by roughly the bulk modulus – about \(10^9\) in SI – so each source channel is named for the quantity it injects (volume_rate / pressure_rate). The receiver side keeps the plain receiver_locations_p/receiver_amplitudes_p, deliberately: it names a field the receiver records, there is no second quantity for it to be confused with, and renaming it would be churn with no defect behind it.
acoustic_vti() records its p receiver as \(-(\sigma_V+(n-1)\sigma_H)/n\), minus the mean normal stress, which is the same mean-stress definition used by the elastic family; it is not a claim that different constitutive equations produce identical pressures. It is formed from the two recordings after the run rather than inside the time loop, so it costs nothing that recording the two stresses would not.
Receivers record the corresponding field itself: the scalar field \(u\), pressure, a velocity component, or (for acoustic_vti()) a normal stress. The elastic family’s receiver_amplitudes_p is the negative arithmetic mean of the recorded normal stresses, so it is a pressure in the usual compression-positive sense. In 2-D that mean is over the two in-plane normal stresses only, so it is the in-plane mean normal stress rather than the thermodynamic pressure; the two differ whenever the out-of-plane stress does not equal the in-plane average.
Every amplitude is a density¶
Every entry above is a density – a per-unit-volume quantity – because a source amplitude is added to a single grid cell with no cell-volume normalisation (see Source and receiver convention). Writing a continuum point source of strength \(S(t)\) as \(S(t)\,\delta(\mathbf{x}-\mathbf{x}_s)\) and representing \(\delta\) on the grid as \(1/\Delta V\) at the source cell, the number you supply is \(S(t)/\Delta V\).
Two consequences follow. First, the same amplitude is a different physical source on a different grid. Halve the grid spacing in 2-D and the equivalent point-source strength falls by a factor of four, so an unchanged amplitude gives a receiver amplitude four times smaller. Second, the 2-D and 3-D meanings differ in exactly the way the last column of the table records. The units of the number you supply are the same, but a 2-D model is a 3-D model that is invariant along the omitted axis, so the source is a line source, and \((\text{amplitude})\times\Delta V\) is a strength per unit length of that line rather than a total.
Converting between propagators¶
Source conversion requires specifying which observable is to agree. The relations below are continuum relations for a homogeneous isotropic medium, compatible zero initial fields, no body-force sources, and an unbounded domain (or the direct arrival before different boundary responses matter). They are not identities between finite-difference outputs on different grids: spatial stencils, time staggering, and numerical dispersion can differ.
Let \(s(\mathbf{x},t)\) be acoustic volume_rate (units \(\mathrm{s}^{-1}\)), \(a(\mathbf{x},t)\) elastic pressure_rate (\(\mathrm{Pa/s}\)), and \(f(\mathbf{x},t)\) the scalar PDE source (\(\mathrm{Pa/m^2}\) when \(u\) is pressure). Define
Here \(n\) is the number of propagated spatial dimensions; \(M,K_n\) and \(A\) have pressure units. In 2D, \(K_n=\lambda+\mu\) relates the in-plane mean normal stress to in-plane dilatation; it is not the 3D thermodynamic bulk modulus \(K_3=\lambda+2\mu/3\).
Acoustic pressure to scalar. For acoustic wavespeed \(v_P\) and bulk modulus \(M\), differentiating \(p_t=-M\nabla\cdot\mathbf{v}+Ms\) and using \(\rho\mathbf{v}_t=-\nabla p\) gives
The second relation follows because scalar solves \(Lu=-v_P^2 f\). It matches the continuum pressure, including at the source, under the stated assumptions.
Elastic mean stress to scalar. Write \(p=-\operatorname{tr}_n\sigma/n\) and \(\theta=\nabla\cdot\mathbf{v}\). The elastic equations give
and therefore
\(L^{-1}\) denotes the causal, zero-initial-data wave solution operator. The second term is local to the source support. At receivers outside that support, elastic pressure can consequently be reproduced by scalar using
At the source this formula alone omits the local \(A\) term. A scalar source reproducing the full continuum pressure is instead \(f=-\partial_t a/v_P^2+(M-K_n)\nabla^2 A/M\); its spatial derivative makes it a different distributed source construction. These formulas do not extend unchanged across heterogeneous interfaces or elastic mode-converting boundaries.
Same imposed volume rate versus same pressure waveform. An isotropic injected strain rate \((s/n)I\) contributes \(-K_n s I\) to the elastic stress-rate equation, so \(a=K_n s\) represents that same imposed volume rate. For \(\mu>0\) it does not make elastic mean pressure equal to acoustic pressure at speed \(v_P\). Matching their off-source pressure instead requires \(a=M^2s/K_n\), as follows by equating the two coefficients of \(L^{-1}\partial_t s\) above. For a fluid (\(\mu=0\)), \(K_n=M\), the local term vanishes, and these distinctions collapse to \(a=Ms\) and \(f=-\rho\partial_t s\).
All three sources are still grid densities: apply the same cell-volume normalisation if comparing a fixed integrated point or line source.
When each source sample is applied¶
Distinguish the user sample interval \(dt\) from the internal interval \(h=\texttt{out.info.internal_dt}=dt/R\), where \(R=\texttt{out.info.step_ratio}\). The following describes the leapfrog clock before any optional time tapering or resampling filter effects.
On staggered grids, internal receiver index \(k\) records stress/pressure at \(kh\) and velocity at \((k-\tfrac12)h\). Source index \(k\) drives an update:
a force-density sample drives \(\mathbf v^{k-1/2}\) to \(\mathbf v^{k+1/2}\), and is evaluated at \(kh\);
a pressure/stress-rate sample drives stress from \(k\) to \(k+1\), and is evaluated at \((k+\tfrac12)h\).
The staggered families downsample without an additional component-dependent half-step correction. Thus returned pressure sample \(j\) has nominal time \(j\,dt\), while returned velocity sample \(j\) has nominal time \(j\,dt-h/2\), not \(j\,dt-dt/2\) when \(R>1\). To supply a pressure-rate waveform with a specified physical time origin, evaluate its user-rate samples at \(j\,dt+h/2\); the bandlimited internal interpolation then supplies the corresponding internal midpoint samples. Force samples use \(j\,dt\). The remaining finite-record interpolation and filtering assumptions are described in Discretisation: the exact scheme.
On the regular scalar grid, source and receiver index \(k\) both have nominal time \(kh\); there is no staggered half-step offset. Comparing a staggered solution to a scalar or analytic reference requires accounting for this distinction. A half-internal-step timing error has narrowband relative waveform error of order \(\omega h/2\), where \(\omega\) is angular frequency, not \(\omega dt/2\) unless no substepping occurs.
Scalar¶
The scalar propagator solves the constant-density scalar wave equation \(\nabla^2u-\ddot u/v^2=f\), with \(f\) the supplied source amplitude (internally Deepwave scales each source sample by \(-v^2h^2\) at the source cell; you do not apply this yourself). The receiver data is the scalar wavefield sampled at the receivers. Source and receiver therefore share a consistent set of units, set by the units of your source amplitudes.
The Born (scalar_born) propagator follows the same convention for its background source, and its scattered output is the linearised (single scattering) wavefield in the supplied velocity-perturbation direction scatter, which has the same velocity units as v.
Acoustic (variable density)¶
The acoustic propagator uses a velocity–pressure (first-order) formulation on a staggered grid. You may inject a volume-injection source (source_amplitudes_volume_rate, added to the pressure update after scaling by \(h\,K\)) and force-density sources on any velocity component (source_amplitudes_z/_y/_x), and you may record pressure and/or the particle-velocity components, each with its own locations. See the Acoustic Wave Equation API page for the exact ordering of the returned receiver tensors.
Elastic (isotropic)¶
The elastic propagator injects and records several distinct quantities, each with its own source_amplitudes_* / *_locations / receiver arguments:
*_z,*_y,*_x: force-density sources, with units \(\mathrm{N/m^3}\). Internal sample \(k\) acts at \(kh\), while the velocity receiver of that internal index records \((k-1/2)h\).source_amplitudes_pressure_rate: pressure rate, with units Pa/s, added with a negative sign to every normal stress. Internal sample \(k\) acts at \((k+1/2)h\), while the pressure receiver records \(kh\).
See When each source sample is applied for the user-rate clock when \(h<dt\).
Receivers can record the particle-velocity components and/or pressure. The returned receiver tensors are ordered pressure-first; see the Elastic wave equation API page for the precise return order, which differs from the acoustic propagator.
Viscoacoustic, acoustic VTI and elastic VTI¶
deepwave.viscoacoustic uses a direct pressure-rate source \(a\) (source_amplitudes_pressure_rate, Pa/s). With \(\theta=\nabla\cdot\mathbf v\), the sourced constitutive equations are
Here \(r_l\) has units Pa/s and \(\tau_l\) is dimensionless. The source adds \(h a\) only to pressure. It is not a prescribed eigen-volume rate, which would replace \(\theta\) by \(\theta-s\) in both equations. In the lossless limit a pressure-rate source \(a\) matches an acoustic volume-rate source \(s=a/K\). Velocity force-density sources and pressure/velocity receivers follow the acoustic conventions.
The quality factor q is dimensionless. Single-SLS q_freq and GSLS q_band endpoints are cyclic frequencies (Hz for seconds; inverse microseconds for a microsecond time unit). q_freq is required for one mechanism, whereas multiple mechanisms require q_band. The supplied wavespeed determines the relaxed modulus \(K=\rho v^2\), not the phase speed at the reference frequency; see Viscoacoustic (attenuating) wave equation and Material models and physical scope.
deepwave.acoustic_vti takes stress-rate sources on its two normal stresses rather than a pressure source. Because stress is tension-positive, the same wavelet that is compressional as an elastic source_amplitudes_pressure_rate is dilatational as an source_amplitudes_sigma_*.
deepwave.elastic_vti follows the same source/receiver conventions as the elastic propagator. Its stiffness models (c11, c13, c33, c55, and c66 in 3-D) have units of pressure, e.g. Pa if you are working in SI units, exactly like lamb/mu in the isotropic elastic propagator.
Choosing units for ultrasound¶
The equations allow any consistent system of units, provided material, source, and receiver quantities are converted together. Equivalent physical results are subject to floating-point conditioning and representability, and to numerical near-zero guards in vacuum-safe inverses and material averages; arbitrary changes of unit scale need not give identical numerical results. Seismic work uses metres and seconds; ultrasound and NDT work at millimetre/microsecond scales, for which a millimetre–microsecond–megahertz system gives convenient speed, time, and frequency values (a water speed of 1.5 mm/µs, dt = 0.01 µs, pml_freq = 3.0 MHz). When pml_freq is not given, Deepwave infers it from your source amplitudes, so it comes out in your own frequency unit automatically – 3.0 in the MHz system above, 3e6 in SI – rather than the fixed seismic 25 that Deepwave 0.0.27 and earlier used regardless of scale. Pass it explicitly anyway for anything you publish: inference makes the boundary depend on the wavelet (see Choosing pml_freq). The Ultrasound and NDT (coming from k-Wave) page gives a full worked unit example, the excitation wavelets (tone_burst(), gaussian_pulse()), the transducer-array helpers, and a concept mapping for readers coming from k-Wave.
Staggered grids and time staggering¶
For the acoustic and elastic propagators the field components live at different half-cell positions in space and (for velocities) at half steps in time. This is handled internally, but it is the reason the velocity receivers are recorded at \(j\,dt-h/2\) while pressure receivers are at \(j\,dt\), and the reason the same-channel source samples act half an internal step later (see When each source sample is applied). When placing sources or receivers off-grid with floating point coordinates (see Usage), Deepwave’s Hicks interpolation can be made aware of these half-cell component offsets.
The offsets also set the exact form that source–receiver reciprocity takes when the two channels differ. Swapping a scalar source with a pressure receiver, or one velocity channel with another, reproduces the trace exactly; swapping a scalar source with a velocity receiver reverses the sign and shifts the internal trace by one internal sample, of duration \(h\). After resampling this is a fractional \(1/R\) of a returned sample when \(R>1\), not an integer shift of the returned array. The two fields are recorded half an internal step apart and the source samples sit half an internal step later still. On the elastic family, whose scalar source channel injects a pressure rate rather than a volume rate, that relation also carries a factor of the bulk modulus at the scalar channel’s own point. These are exact properties of the discrete scheme, for sources and receivers anywhere on the grid (including inside the absorbing boundary), and tests/test_reciprocity.py checks every channel pairing against them. They hold as stated only where a single real modulus converts between the scalar source rate and pressure, so on an attenuating medium the mixed pairings depart by O(1/Q), and on an anisotropic one by the departure of the stiffness from isotropy.