The interaction between electromagnetic/acoustic waves and natural surfaces in radar/sonar remote sensing is strongly influenced by interface roughness, which has to be considered relative to the wavelength λ (Engman and Wang, 1987), as expressed by the dimensionless roughness parameter at incidence θ 1 g = 2 k   ζ   cos ⁡ θ 1 (4) where k = 2 π / λ is the incident signal wavenumber and ζ is the root mean square (RMS) height of interface elevations around the average depth value. This roughness parameter (aka “Rayleigh parameter”) actually scales the surface physical relief to the signal wavelength.

Large-scale seafloor relief features (i.e., at topographic scale) are to be treated deterministically and are not considered here; for acoustic signal backscatter, the concept of roughness concerns relief scales that cannot be resolved by bathymetry sonar measurements. According to interface roughness magnitude, different reflection/scattering regimes occur, corresponding to various physical models. For rough interfaces ( ζ >> λ ), the wavelength is much smaller than the dominant relief features, so that the interface can be considered as a continuum of apparent “facets” acting like locally plane reflectors (see, e.g., Brekhovskikh and Lysanov, 1982). This is the fundamental concept of the “Facet Method” based on Kirchhoff’s approximation, a geometric approach considering the local interface curvature and assuming that specular reflection occurs then on the locally tangent plane; it is valid for large curvature radii of the interface compared to wavelength and works most effectively for steep incidence angles (Beckmann and Spizzichino, 1963). At the other end of the roughness scale ( ζ << λ ), the interface roughness just infringes small statistical phase perturbations to the incident wavefront; this situation is addressed by the Small Perturbation Method, involving small phase shifts of the incident wave caused by surface roughness and modeled as a stationary function; the main contribution to scattering is then Bragg’s phenomenon. This case corresponds to small values of g and oblique incidences (Rice, 1951; Valenzuela, 1978).

In the model proposed here, the interface scattering component is based on a two-scale approach involving both the facets (or Kirchhoff’s) theory at “small scale” and the small perturbation method (SPM) at “micro-scale”, combined through a physically motivated transition function. The detailed formulation of this transition is given in Section 2.4.2. Following radar literature as well as previous acoustical models (e.g., Jackson et al., 1986), this approach acknowledges that different scattering mechanisms dominate at different roughness scales and angles of incidence, a long-established concept both in radar and sonar remote sensing.

As a preliminary to the presentation of the two models, it is useful to establish the relationship between the components of interface roughness, linking its roughness spatial spectrum with the variances of interface slopes and elevations.

The seafloor roughness is characterized by its power spectrum Ω κ usually written in the form: Ω κ = a 2 κ − γ (5) where κ is the interface roughness spatial frequency and γ takes values typically between 3.0 and 3.5 (de Moustier and Alexandrou, 1991). The exponent γ directly and strongly influences acoustic backscatter modeling and sonar data interpretation based on the physical characteristics of the seabed. The value γ = 3.5 is often used for natural sedimentary seabed and areas with gradual variations in roughness, while γ = 3 may be more appropriate for rockier seabeds and areas with more uniform roughness scales. The scaling factor a 2   gives the magnitude of the elevation distribution corresponding to the roughness spectrum.

The spatial spectrum Ω κ has to be split in order to address the different roughness regimes. Disregarding its lower-frequency part associated with deterministic large-scale topography, one can consider a medium-frequency part (“small-scale”) adapted to the Kirchhoff’s approach and a high-frequency part (“micro-scale”) associated with the SPM. These two parts split around a transition “cut-off” value κ _C that is not intrinsic to the interface relief alone but also depends on the acoustical frequency. Following Novarini and Caruthers (1994), the cut-off spatial frequency corresponds to the condition that the Rayleigh parameter (i.e., Formula 4 considered at normal incidence) for microroughness equals unity: g 2 = 4 k 2   ζ μ 2 = 16 π 2 f 2 c 2   ζ μ 2 = 1 (6)

Formula 6 makes it possible to obtain the effective ζ μ 2 at a given frequency f. Moreover, the microroughness elevation variance ζ μ 2 is defined as ζ μ 2 = 2 π ∫ κ c κ H Ω κ κ d κ (7) and the small-scale facet-slope variance as: δ f 2 = 2 π ∫ κ L κ C Ω κ κ 3 d κ (8) where κ _L and κ _H are respectively the low- and high-limit values of the spatial frequency spectrum.

Formulas 5-8 will now be used in order to detail the roughness backscatter models, at both scales in a consistent way, as described in the next paragraphs.

For the small-scale facet (or Kirchhoff) regime, the backscattering cross-section σ f is primarily given (Brekhovskikh and Lysanov, 1982) by the classical formula: σ f θ 1 = V 2 0   exp − tan 2 ⁡ θ 1 / 2 δ 2 / 8 π δ 2 ⁡ cos 4 ⁡ θ 1 (9) where V 2 0 is the intensity reflection coefficient at normal incidence; the exponential term exp − tan 2 ⁡ θ 1 / 2 δ 2 / 2 π δ 2 is a Gaussian distribution depicting the assumed probability distribution of surface slopes, in which the facet-slope variance δ 2 controls the distribution width; and the c o s 4 θ 1 term accounts for geometric spreading and projection effects.

In this “facets model”, the key parameter is the standard deviation of the facet-slope distribution, which directly controls the angular width of the specular lobe around normal incidence θ 1 = 0. Importantly, this model usually assumes a Gaussian distribution of slopes (Brekhovskikh and Lysanov, 1982), which simplifies the mathematical treatment while still providing a reasonable approximation for many natural surfaces ¹ . The key point of this facets theory is that the seafloor backscatter angular response near normal incidence (typically up to 35°) is proportional to the slope distribution, as the backscatter strength in this region is dominated by facets oriented normally to the transducer location. This relationship between slope statistics and angular response provides a strong physical meaning to the model parameters, allowing for an intuitive interpretation of seafloor characteristics from acoustic data. However, although this is frequently a good approximation, the roughness-slope distribution is not systematically Gaussian, as non-Gaussian angular responses are frequently observed in field data.

It must be emphasized that the relationship between interface-scattering and acoustic-wavelength is fundamental to the model’s physical relevance. In this respect, the variance of the facet-slope distribution has to be scaled by the normalized frequency, reflecting a key insight about rough-interface scattering: for a given sonar configuration, the backscatter effectively corresponds to only fractions of the full distributions of slopes and elevations, limited both by the sonar footprint extent on the interface and by the frequency-scaled roughness. This wavelength-dependent filtering effect has been well documented in previous research works (Novarini and Caruthers, 1998; Shaw and Smith, 1990) demonstrating that acoustic backscatter is mostly sensitive to roughness components with scales comparable to the acoustic wavelength. Using both the existing theories and the experimental data from our Concarneau dataset (Fezzani et al., 2025), we tried to determine frequency-scaling factors in order to minimize the frequency dependence of the derived roughness parameters. This approach (also applied here below to the volume backscatter contribution) makes it possible for our ESAB model to extract physically meaningful seabed parameters that remain consistent across the operational frequency range of typical seafloor-mapping sonars.

In order to build the ESAB model, it is now proposed to improve the classical and widely accepted expression of the facets model on three important points, often disregarded:

At the micro-roughness scale, the backscatter component generated by the interaction between the acoustical wave and the interface roughness is dominated by the “Bragg effect” resulting from the resonance between the incident acoustical wavelength projected on the interface and the corresponding spatial frequency spectrum of the interface relief (Thorsos and Jackson, 1989). For one signal frequency and incidence angle, the obtained backscatter cross-section is classically expressed (e.g., Novarini and Caruthers, 1998) as: σ B   θ 1 = V 2 θ 1 π 2 k 4 co s 4   θ 1   Ω k π sin   θ 1   (13) using the same notations as above.

When assuming (as in Equation 5) a negative-exponential roughness spectrum Ω κ = a 2 κ − γ , a direct analytical relationship can be found between the roughness-slope standard deviation and the RMS roughness elevation. This property is particularly valuable in our case since it allows the expression of the Bragg scattering component without the need to introduce an additional independent roughness parameter for this particular regime. This simplification respects the theoretical integrity of the model while reducing its parametric complexity, making it more practical for inversion problems and preserving its physical meaning.

Using the relationships established above between δ 2 f and ζ μ 2 , the micro-scale backscatter cross-section given by Equation 13 can be developed into a function of the facet-slope variance   δ 1 2 alone: σ B   θ 1 = D 3   V 2 θ 1   co s 4   θ 1 si n γ   θ 1   δ 1 2 D 4 (14) with: D 3   = 2 γ − 5 π γ − 3 4 − γ D 4 γ − 2 1 − D 4     D 4 = γ − 2 2

This microroughness component σ μ shows a characteristic angular dependence proportional to co s 4   θ 1 /  si n γ   θ 1 , leading it to strongly diverge at vertical incidence (in 1 / si n γ   θ 1 ), and to rapidly decrease at grazing angles (in co s 4   θ 1 ); hence its practical influence is restricted to intermediate oblique angles. It also depends, as expected, on the interface reflection coefficient, on frequency, and on the roughness spectrum reference power. Practically, for numerical implementation the divergence around vertical incidence has to be controlled by a dedicated weighting term (see §2.4.2 below).

This microroughness backscatter component σ μ can interestingly be expressed for particular values of the spectrum exponent γ resulting in some simplification:

• γ = 3.0   σ B   θ 1   ∼   V 2   θ 1   co s 4   θ 1 si n 3   θ 1   δ 1 2 0.5

Finally, the ESAB backscattering strength (Equation 16) combines both interface ( σ I ) and volume ( σ V ) contributions before converting their intensity summation to decibels: B S θ 1 , f = 10  lo g 10   σ I θ 1 , f + σ V   θ 1 , f (16) where each term comes from Equations 12‐15 and accounts for the various improvements, simplifications and approximations presented and discussed in the sections above. All these terms are now explicitly dependent on frequency.

As a reminder from the previous paragraphs, it must be emphasized that this final formulation (Equation 16) only requires four input parameters related to the seafloor configuration: the water-sediment impedance contrast z; the sediment volume scattering strength parameter μ; and the two interface roughness facet-slope standard deviations δ 1 and δ 2 .

The transition between the two roughness scattering regimes (facets and Bragg) presented above is managed through a simple sigmoid angle-dependent function (Equation 17) featuring a frequency-dependent crossing angle. This ensures physical continuity between the two regimes and hence better represents the overall seafloor scattering behavior. As proposed in (Mourad and Jackson, 1989), a crossing angle θ x can be defined by σ f 1 (Equation 9) dropping by −15 dB from its normal incidence value σ f 1 0 . The weighting function writes: Σ θ 1 = 1   / 1 + exp − 180 / π   θ 1 − θ x (17) where θ 1 and θ x are expressed in radians.

The resulting form of the interface roughness scattering σ I s (featuring the sigmoid weighting) is finally expressed in Equation 18. Note that, before being combined with σ f 1 through the sigmoid (Equation 17), σ f 2 is summed with σ B , as both are triggered beyond the crossing angle θ x . σ I s   θ 1 =   σ B   θ 1 + σ f 2   θ 1   ]   Σ   θ 1 + [ 1 − Σ   θ 1 ] σ f 1   θ 1 (18)

The ESAB model uses a second sigmoid function in order to reduce the volume scattering contribution near normal incidence, complementing the constrained model inversion approach (Section §3) where all parameters (z, δ 1 , δ 2 , m 0 ) are bounded within physically meaningful ranges typical of seafloor sediments. At normal incidence, both decreasing roughness and increasing impedance contrast can produce higher backscatter strength, creating an inherent ambiguity (Figures 5A,C). Without sigmoid weighting, adding volume scattering in this angular region could lead the inversion process to converge on parameter combinations that, while mathematically valid, would violate physical expectations for seafloor properties. By effectively separating the angular domains where different scattering mechanisms predominate (Fonseca and Mayer, 2007), the sigmoid weighting allows the interface parameters (z, δ 1 , δ 2 ) to be primarily determined from near-nadir returns ( θ < θ v , while the volume parameter m 0 is mainly controlled by oblique angle responses ( θ > θ v , improving the numerical stability of the simulated annealing inversion process (see §3.1.1). A limit angle value of θ v = 0.1745 radians (or 10^o) was found appropriate for the inversion process. The division by 2 forces a smoother transition in the sigmoid. The volume backscatter cross-section featuring the sigmoid term finally writes: σ V s   θ 1 = σ V   θ 1   /   1 + exp − 180 / π   θ 1 − θ v / 2   (19)

The final expression for the ESAB backscatter strength is given in Equation 20; it includes the interface and volume components corrected by the sigmoid functions presented above. B S θ 1   = 10  lo g 10   σ I s θ 1   + σ V s θ 1   (20)

Although these transition terms are, strictly speaking, out of the physical modelling, they are of paramount importance for the numerical stability of the model, especially when applying it to inversion purposes.

The simplified adaptation of both sonar- and radar-inspired modeling techniques presented above is expected to provide a comprehensive description of seafloor backscattering while ensuring computational efficiency, minimizing the formal complexity (including the input physical parameters) and flexibly, and preserving physical consistency across different angular and frequency regimes. The model versatility makes it able to represent a wide range of seafloor types through appropriate parameter selection, while its theoretical foundation in scattering theories justifies its relevance and ensures robust physical behavior.

A major interest of direct backscattering computations is the demonstration (see, e.g., Figure 5) of how each parameter distinctly affects the angular response, making possible a clear phenomenological interpretation of intricate processes. The impedance contrast z primarily influences backscatter level near normal incidence; this is illustrated in Figure 5A, where increasing z from 1.3 to 4.1 produces significantly stronger returns (+13 dB) in the specular lobe while maintaining comparable responses (within 3 dB) at oblique angles. Notably, this enhancement of normal incidence backscatter by increasing impedance is similar to the effect of decreasing surface roughness, through different physical mechanisms (Figure 5C). The volume parameter μ shows opposite behavior (Figure 5B), mainly affecting oblique angles (>30°) with minimal impact on near-nadir returns (where it has been, by the way, minimized by the sigmoid weighting). The facet-slope standard deviations ( δ 1 and δ 2 ) shape the overall angular response: increasing δ 1 reduces the specular peak while broadening the angular distribution, a characteristic property of diffuse scattering from small-scale roughness. Meanwhile and as expected from its introduction and justification, δ 2 more strongly influences the response beyond 20° incidence angles (Figures 5E,F) and significantly departs the overall angular response from the canonical Gaussian shape. Interestingly, the red curve (z = 1.7) in Figure 5A shows higher backscatter at oblique angles due to the impedance-dependent attenuation model. Low impedance ratios result in lower attenuation coefficients, allowing deeper acoustic penetration and enhanced volume scattering contributions, which dominate the oblique-angle response and produce the elevated off-nadir backscatter seen in this curve.

The frequency dependence shows behavior similar to roughness variation since the model uses the classical Rayleigh parameter, which relates wavelength and surface roughness (Equation 4). At the lowest frequencies, the angular response shows a pronounced specular lobe with rapid decay at oblique angles (Figure 5D). As frequency increases, the angular response broadens and the specular component diminishes. Unlike pure interface roughness dependence, the frequency increase also enhances volume scattering, which scales with the k 0 0.65 term discussed in §2.3, demonstrating how the model captures the distinct frequency dependence of surface and volume scattering mechanisms.

A predecessor to the ESAB presented here, the Generic Seafloor Acoustic Backscatter (GSAB) model (Lamarche et al., 2011) already proposed a more practical approach aiming at a functional description of the main backscatter regimes (near-specular lobe, oblique-angle plateau and grazing-angle fall-off) using simple mathematical expressions and a limited set of six input parameters (possibly only four significant ones), and effectively fitting a majority of backscatter ARCs. The GSAB model comprises three components: a normal law A   exp − θ 1 2 / 2 B 2   describing near-specular reflection; a powered-cosine component C  co s D   θ 1 for the backscatter plateau and fall-off; and a second normal law E   exp − θ 1 2 / 2 F 2   effective at intermediate angle regimes. While really practical and versatile for empirical data fitting, the GSAB parameters lacked direct physical interpretation leading to quantified seafloor properties, and missed an explicit frequency dependence (in the sense that each frequency needs one specific parameter set).

ESAB builds upon the GSAB’s functional approach while providing clear physical interpretation through a small number of parameters representing three fundamental seafloor properties: acoustic impedance for water-sediment contrast; slope values for interface roughness; and a volume scattering parameter for sediment inhomogeneity. This limitation to a remarkably low number of input parameters was obtained by using ad hoc relationships linking the various factors involved in the physical processes.

Moreover, unlike GSAB, the ESAB model incorporates explicit frequency dependence of interface and volume scattering phenomena, defining sediment intrinsic physical properties valid across frequencies. Applied to the Concarneau dataset (Figure 6), ESAB demonstrates physical consistency at 70 and 290 kHz with stable impedance values for sandy mud (z = 1.72) and medium sand (z∼1.90), while GSAB parameters expectedly show greater variability between frequencies. (See Supplementary Material ‐ Comparison of GSAB and ESAB).

Despite its intrinsic limitations compared to ESAB, the GSAB model still features interesting properties. Its excellent versatility in fitting very different ARC shapes makes it very useful for a parametric description of backscatter; the resulting parameters, although not directly physically interpretable as in ESAB, are however usable for some seafloor type classification or backscatter map segmentation. Moreover, GSAB is still usable for data obtained from non-calibrated echosounders, provided that the result interpretation bears only on parameters not depending on absolute levels: for instance, the specular lobe apertures (B or F); or the grazing angle decrement (D); or the contrast between levels at specular (A) and at oblique (C); all these are expected to be efficient descriptors of the seafloor type. Finally, GSAB can be used as a preliminary analysis tool to be applied before a full ESAB inversion for (1) preconditioning the ARCs under an averaged smoothed shape, and (2) possibly pre-selecting data from seafloor configurations found to be ill-adapted to ESAB.

Note that the final inversion results presented in Section 3 (Figure 8) used GSAB-fitted curves rather than the sparse experimental data points shown in Figure 6. The GSAB fitting provides complete angular response curves covering all angles, which better represent typical multibeam echosounder data that continuously samples the entire angular range. This approach allows the same simulated annealing procedure to be applied to both single-beam and multibeam datasets.

The ESAB model implementation follows a structured approach that computes each scattering component separately before combining them into the total backscatter response. The algorithm operates on input angles and includes safeguards against numerical singularities. The implementation calculates reflection and transmission coefficients using Snell’s law with the sound speed ratio relationship given in Equation 3. For the interface scattering components, facet contributions at both scales ( σ f 1 and σ f 2 ) are computed using frequency-dependent slope variances following (Equation 11). The roughness effect on reflection is incorporated through the coherence loss factor 1/e discussed in Section 2.4.2.

The Bragg scattering component is considered at oblique angles, with special handling of near-normal incidences to avoid mathematical singularities. The transition between facet- and Bragg-dominated regimes is managed through a crossing angle ( θ x ) determined by the specular lobe width at −15 dB. This transition is implemented using the sigmoid function described in Equation 17. Volume scattering follows the formulation in Equation 15, with frequency-dependent attenuation, again calculated from impedance-based relationships. An optional sigmoid function (Equation 19) reduces volume contribution near normal incidence, ensuring stable inversion convergence by separating the angular domains where different mechanisms dominate.

The optimization cost function uses a weighted RMS approach, with additional emphasis on near-nadir angles (0°–20°) to enhance the accuracy of interface parameter determination. Parameter constraints are enforced throughout the optimization process, ensuring physical realism by maintaining impedance within realistic bounds (1.0–14.0), volume backscatter parameter between −15 dB and +15 dB, and enforcing the condition δ 2 > δ 1 . The roughness parameters are bounded according to frequency-dependent limits that reflect physical constraints on surface slopes.

The ESAB model inversion uses simulated annealing optimization to determine the three physical parameters (z, μ, s ₁) and two distribution parameters (δ ₁, δ ₂) by fitting theoretical angular response curves to measured backscatter data. This optimization technique, inspired by metallurgical annealing processes (Kirkpatrick et al., 1983), can escape local minima through a so-called “temperature” parameter that gradually decreases, making it suitable for, e.g., the multimodal nature of seafloor parameter inversion. The present implementation uses an initial temperature of T = 0.1, N _c = 15 cooling cycles, and 10,000 iterations per cycle, with temperature decreasing according to T/log₂(1+n _c ), where n _c is the current cycle number from 1 to N _c . Parameter perturbations follow normal distributions with controlled scaling to ensure efficient exploration of the solution space while respecting physical constraints. The acceptance probability of the fitting process follows the Metropolis criterion (Metropolis et al., 1953), allowing occasional acceptance of suboptimal configurations to escape local minima and explore the broader solution space.

One of the objectives in building ESAB was to keep the number of model’s parameters as low as possible, in order to facilitate the inversion from experimental data, and to classify the seafloor types. In the latter respect, the objective is a classification along only three parameters: impedance ratio, roughness parameter and volume parameter.

The parameter inversion process uses a two-step approach to determine the optimal roughness spectrum exponent γ. First, a model inversion is performed for all seven sites across 18 frequencies using γ = 10/3 ≈ 3.333, which ensures a linear dependence between slope variance and frequency (see Section 2.2.3). After obtaining the inversion results, the frequency dependence of δ 1 is analyzed and a residual slope in the regression of log ⁡ δ 1 against log ⁡ f is determined for each sediment site (Sites 1–5). An average residual exponent e ₁ is determined, and since δ 1 2 is proportional to f 2 e 1 , this value is used to calculate the corrected γ through the relationship γ = 8 + 2 e 1 / e 1 + 2 , derived from Equation 10. This analysis yielded here γ = 3.508.

A. Sediment bulk properties

Regarding the sediment properties, the model hypothesizes a viscous fluid medium that can be described using three parameters (density, velocity, absorption). However, the model inversion only uses one parameter (impedance, i.e., the density-velocity product); the three original parameters have to be retrieved from geoacoustic considerations (see Supplementary Material ‐ Geoacoustical modelling applied in ESAB), and can be obtained as intermediate results of the inversion process.

B. Interface roughness

Regarding roughness, the classical theories (“facets” and “small perturbations”) explicitly need several input parameters, namely, the reflection coefficient and the variances of roughness slopes and elevations. Following the works by Novarini and Caruthers (1994), Novarini and Caruthers (1998) we have expressed these two roughness variances as frequency-dependent quantities, both defined by the roughness power spectrum that is mainly controlled by the exponent γ . The cut-off frequency separating the “facets” and “Bragg” regimes is defined by the Rayleigh parameter equalling one. Hence, for a given γ value, all the other quantities can be obtained. A practical difficulty is to invert the γ value together with the other parameters. The strategy proposed here is to process the inversion in several steps:

1. Fix the γ value at one arbitrary value, e.g., 10/3 (see §2.2.4). Run the inversion until obtaining the s ₁ parameter as a function of frequency

Of course, γ can be extracted as another intermediate parameter; it is possibly an interesting descriptor of the local roughness.

C. Sediment volume

The volume backscatter model used in ESAB (Equation 15) features the water-sediment transmission coefficient, the absorption coefficient, and the volume scattering parameter. Since the latter one is expected to be frequency dependent while a constant quantity is desired, this dependence is a priori described by a power of frequency that is here taken equal to an average value extracted from the analysis by Fezzani et al. (2021). The other model parameters are available from the geoacoustical modelling. Finally, the quantity extracted from the model inversion is the volume backscatter parameter μ normalized at frequency f ₀.