**2. First-Order Bistatic Electric Field Scattered from an Ocean Surface with Arbitrary Heights**

For the purposes of this work, a conductive rough surface defined as *z* = *f*(*x*, *y*; *t*) is considered. *f*(*x*, *y*; *t*) is a zero-mean, time-varying, two-dimensional random variable representing the ocean surface displacement from the sea level *z* = 0 [27]. In general, the equation for the electromagnetic propagation over a rough surface is defined as [16]

$$\mathcal{NL}^{-1}\left[\frac{\mathcal{LN}\mathcal{L}^{-1}\left[2u\mathcal{F}\_{\text{xy}}(\mathbf{E}\_{s}^{\Xi})\right]}{u+jk\Delta}\right] = \mathbf{E}\_{n}^{+} + \mathcal{NL}^{-1}\left[\frac{\mathcal{LN}\left[\frac{\nabla\_{\text{xy}}\left(|\mathbf{n}|E\_{n}^{+}\right)}{|\mathbf{n}|^{2}}\right]}{u+jk\Delta}\right],\tag{1}$$

where **n** is a vector normal to the surface, understood here to be the ocean surface, and is defined as

$$\mathbf{n} = \dot{\mathbf{z}} - \nabla f(\mathbf{x}, \mathbf{y}; \mathbf{t}). \tag{2}$$

<sup>F</sup>*xy*(·) is the spatial Fourier transform in the *xy*-plane, **<sup>E</sup>***z*<sup>−</sup> *<sup>s</sup>* is the source electric field at the point *<sup>z</sup>* <sup>=</sup> *<sup>z</sup>*<sup>−</sup> <sup>&</sup>lt; *<sup>f</sup>*(*x*, *<sup>y</sup>*; *<sup>t</sup>*), <sup>∀</sup>(*x*, *<sup>y</sup>*; *<sup>t</sup>*), **<sup>E</sup>**<sup>+</sup> *<sup>n</sup>* is the normal electric field immediately above the ocean surface, *k* is the radar wavenumber, *u* is defined as

$$
\mu = \sqrt{K^2 - k^2},
\tag{3}
$$

where *K* is the surface wavenumber, Δ is the surface impedance, N is the normalizing operator defined as

$$\mathcal{N}\{\mathbf{A}\} = \mathfrak{fin}\mathfrak{i}\cdot\mathbf{A}, \,\forall \mathbf{A}\_{\prime}$$

and L is an invertible operator, defined as [16]

$$\mathcal{L}\{\mathbf{A}\} = \mathcal{F}\_{\mathbf{xy}} \left\{ |\mathbf{n}|^{2} \mathbf{A} e^{(z^{-} - f(\mathbf{x}, \mathbf{y}; t)) \mathbf{n}} \right\} , \,\forall \mathbf{A}. \tag{4}$$

Now, defining the inverse operator L−<sup>1</sup> such that it supports an ocean surface with arbitrary heights [28] and proceeding with the derivations for a vertical dipole source in the far-field region defined as

$$\mathbf{E}\_s = \frac{I(\omega)\Delta\ell k^2}{j\omega\epsilon\_0} \mathbf{G}\_0 \mathbf{\hat{z}} \equiv \mathbf{C}\_0 \mathbf{G}\_0 \mathbf{\hat{z}},\tag{5}$$

where *G*<sup>0</sup> is the Green's function solution for the Helmholtz equation in free space [16], Δ is the dipole length, and *I*(*ω*) is the current flowing on the antenna, the first-order electric field for a rough surface with arbitrary heights can be written as [28]

$$\left(E\_n^+\right)\_1 \sim -jkC\_0 \frac{\varepsilon^{\mathbb{Z}(x\_1, y\_1; t)}}{(2\pi)^2} \int\_{x\_1} \int\_{y\_1} \mathfrak{d}\_1 \cdot [\nabla\_{x\_1 y\_1} f(x\_1, y\_1; t)] F(\rho\_1) F(\rho\_2) \frac{\varepsilon^{-jk(\rho\_1 + \rho\_2)}}{\rho\_1 \rho\_2} d\mathbf{x}\_1 d\mathbf{y}\_1. \tag{6}$$

where *F*(*ρ*) is the Sommerfeld attenuation function as presented in [29], and *ζ*(*x*1, *y*1; *t*) is the arbitrary height factor, defined as [28]

$$\mathcal{Z}(\mathbf{x}, y; t) \equiv \mathcal{F}\_{\mathbf{x}\mathbf{y}}^{-1} \left\{ \mathcal{F}\_{\mathbf{x}\mathbf{y}} \left\{ f(\mathbf{x}, y; t) \right\} u \right\} = f(\mathbf{x}, y; t) \mathop{\rm sk}\_{\mathbf{x}\mathbf{y}} \mathcal{F}\_{\mathbf{x}\mathbf{y}}^{-1} \left\{ u \right\},\tag{7}$$

where ∗ *xy* indicates a two-dimensional spatial convolution in the *xy*-plane. The scattering geometry for the first-order electric field is depicted in Figure 1.

**Figure 1.** Scattering geometry for the bistatic first-order electric field.

*Remote Sens.* **2020**, *12*, 667

If a Fourier series representation of the ocean surface is considered, the surface displacement *f*(*x*, *y*; *t*) can be written as [30]

$$f(\boldsymbol{\varrho};t) = \sum\_{\mathbf{K},\omega\_{\mathbf{K}}} f(\mathbf{K},\omega\_{\mathbf{K}})e^{j(\mathbf{K}\cdot\boldsymbol{\rho}-\omega\_{\mathbf{K}}t)} = \sum\_{\mathbf{K},\omega\_{\mathbf{K}}} f(\mathbf{K},\omega\_{\mathbf{K}})e^{-j\omega\_{\mathbf{K}}t}e^{j\mathbf{K}\cdot\boldsymbol{\rho}\cos(\theta\_{\mathbf{K}}-\theta)},\tag{8}$$

with **K** = (*Kx*, *Ky*)=(*K*, *θ***K**) being the wave vector for the ocean surface and *ω***<sup>K</sup>** being the angular frequency obtained from the dispersion relation of ocean surface gravity waves [27]. By substituting (8) into (6), expanding the exponential that contains the arbitrary height factor into a power series and applying an asymptotic perturbation expansion to both the Fourier components of the ocean surface and the electric fields using the surface slope as the perturbation parameter [31], it can be shown that the received first-order electromagnetic and second-order hydrodynamic electric fields contain terms equivalent to those for the small-height case, but an extra term appears in the second-order derivation which relates the first-order cross-section and the arbitrary height factor. The new term can be written as [28]

$$\langle \mathbb{P}(\mathbb{E}\_{n}^{+})\_{\mathrm{II}2} \rangle \sim \frac{k\mathbb{C}\_{0}}{(2\pi)^{2}}\mathbb{f}\_{1}(\boldsymbol{\rho}\_{1};t)\sum\_{\mathbf{K},\omega\_{\mathbf{K}}}f\_{1}(\mathbf{K},\omega\_{\mathbf{K}})e^{-j\omega\_{\mathbf{K}}t}\mathbb{K}\iint\limits\_{\mathbf{x}\_{1}\mathbf{y}\_{1}}\cos(\theta\_{\mathbf{K}}-\theta\_{1})\frac{\mathrm{F}(\boldsymbol{\rho}\_{1})\mathrm{F}(\boldsymbol{\rho}\_{2})}{\rho\_{1}\rho\_{2}}e^{-j\theta\rho\_{2}}e^{j\rho\_{1}[\mathbf{K}\cos(\theta\_{\mathbf{K}}-\theta\_{1})-k]} \cdot d\boldsymbol{x}\_{1}d\boldsymbol{y}\_{1}.\tag{9}$$

Comparing the double integrals in (9) with those in the first-order electric field expression in [16,17], it is evident that they are identical. Therefore, following the same procedure detailed in [17] for the first-order bistatic electric field, the following expression is obtained as

$$\delta\_{\mathbf{k}}(\mathbf{E}\_{\mathbf{n}}^{+})\_{\mathrm{l}\mathbf{2}c} \sim \frac{k\mathbb{C}\_{0}}{(2\pi)^{3/2}}\mathbb{I}\_{\mathbf{1}}(\rho\_{\mathbf{1}};t)\sum\_{\mathbf{K},\omega\_{\mathbf{K}}}f\_{\mathbf{1}}(\mathbf{K},\omega\_{\mathbf{K}})e^{-j\omega\_{\mathbf{K}}t}\sqrt{\mathbf{K}}e^{j\frac{\mathbf{K}}{2}\mathbf{K}}\int\_{\rho\_{\mathbf{s}}}\frac{\mathrm{F}(\rho\_{\mathbf{1}})\mathrm{F}(\rho\_{\mathbf{2}})}{\sqrt{\rho\_{\mathbf{s}}\left[\rho\_{\mathbf{s}}^{2}-\left(\frac{\mathbf{s}}{2}\right)^{2}\right]}}e^{\frac{-j\mathrm{s}\_{\mathbf{s}}\cdot\left[\pm\mathcal{K}\cos\boldsymbol{\mathfrak{q}}-2\mathrm{i}\right]}{2}d\rho\_{\mathbf{s}}}d\rho\_{\mathbf{s}},\tag{10}$$

where *ρ<sup>s</sup>* is defined as

$$
\rho\_5 = \frac{\rho\_1 + \rho\_2}{2},
$$

*φ* is the bistatic angle, defined as the bisection of the angle between *ρ*<sup>1</sup> and *ρ*2, and *ρ* is the vector between the transmitter and receiver, shown in Figure 1.

If, similarly to [16,18], an inverse Fourier transform with respect to the radar frequency is applied to (10) while using the associative property of the convolution, the following expression can be obtained:

$$\begin{split} \mathbb{P}(\boldsymbol{E}\_{n}^{+})\_{12}(t) \sim \frac{1}{(2\pi)^{3/2}} \Bigg\{ \int\_{I} \boldsymbol{\mathcal{F}}\_{l}^{-1} \left\{ \mathbf{k} \mathbb{C}\_{0} \right\} \ast\_{l} \boldsymbol{\mathcal{F}}\_{l}^{-1} \Bigg\{ \sum\_{\mathbf{K},\omega\_{\mathbf{K}}} f\_{1}(\mathbf{K},\omega\mathbf{e}) e^{-j\omega\_{\mathbf{K}}t} \sqrt{\mathbf{K}} e^{j\frac{\theta}{2} \cdot \mathbf{K}} \int\_{\rho\_{\mathbf{s}}} \frac{\overline{\mathbf{r}}(\rho\_{1}) \overline{\mathbf{r}}(\rho\_{2})}{\sqrt{\rho\_{\mathbf{s}} \left[\rho\_{\mathbf{s}}^{2} - \left(\frac{t}{2}\right)^{2}}} e^{\overline{\boldsymbol{x}}/\boldsymbol{\pi}/4} \\ \qquad \cdot \left(\pm\sqrt{\cos\boldsymbol{\phi}}\right) e^{\left[\rho\_{\mathbf{s}}\left[\pm\boldsymbol{K}\cos\boldsymbol{\phi} - 2\boldsymbol{k}\right] \right]} d\rho\_{\mathbf{s}} \right) \Bigg\} \times {}\_{l} \boldsymbol{\mathcal{F}}\_{l}^{-1} \left\{ \zeta\_{1}(\rho\_{1};t) \right\} \,. \end{split} \tag{11}$$

Again, it can be easily observed that the time convolution (∗*t*) inside the braces is similar to the one presented in Equation (5) of [18], with the exception of the time-dependent exponential term for the Fourier series expansion of the ocean surface; although this term will not affect the inverse Fourier transform, it might have an effect on the final convolution. It is easy to show that the additional time-dependent exponential term does not affect the resulting expression, since the added terms in the final expression are significantly smaller than the rest of the terms in the expression. Therefore, substituting the resulting expression for the first-order time-varying electric field in [18] for a pulsed radar source, the expression in (11) becomes

$$
\begin{split} \mathcal{E}(E\_{n}^{+})\_{12}(t) \sim \left\{ \frac{-\langle \eta\_{0}\Delta t\rho\_{0} \rangle^{2}}{(2\pi)^{3/2}} \sum\_{\mathbf{K},\omega\_{\mathbf{K}}} f\_{1}(\mathbf{K},\omega\_{\mathbf{K}}) \sqrt{\mathcal{K}\cos\phi\_{0}} e^{j\omega\_{\mathbf{K}}t} e^{j\xi\_{\mathbf{K}}\Delta\rho\_{0}} e^{j\frac{\theta}{2}\cdot\mathbf{K}} e^{j\rho\_{0}\left(\mathcal{K}\cos\phi\_{0}\right)} \frac{F(\rho\_{0},\omega\_{0})F(\rho\_{0},\omega\_{0})}{\sqrt{\rho\_{0}\left[\rho\_{0}^{2}-\left(\frac{\theta}{2}\right)^{2}\right]}} \\ \qquad \cdot e^{j\pi/4} \Delta\rho\_{0} \operatorname{Sa} \left[ \frac{\Delta\rho\_{\varepsilon}}{2} \left(\frac{\K}{\cos\phi\_{0}} - 2k\_{0}\right) \right] \end{split} \tag{12}
$$

where Sa(·) is the sampling function, defined as

$$\text{Sa}(\mathfrak{x}) = \frac{\sin \mathfrak{x}}{\mathfrak{x}}, \,\,\forall \mathfrak{x},$$

and Δ*ρ<sup>s</sup>* is the patch width on the ocean surface, defined for a pulse radar as

$$
\Delta \rho\_s = \frac{c \tau\_0}{2},
$$

with *c* being the speed of light and *τ*<sup>0</sup> being the radar pulse width. Here, it should be noted that the zero-subscripts in *φ*0, *ρ*01, and *ρ*<sup>02</sup> indicate that the scattering patch is considered to be sufficiently small, allowing variable values at the center of the scattering patch to be taken as representative of their values on the whole patch [32]. Consequently, *ρ*0*<sup>s</sup>* is defined as

$$
\rho\_{0\mathfrak{s}} = \frac{c(t - \frac{\tau\_0}{2})}{2} = \frac{\rho\_{01} + \rho\_{02}}{2}.
$$

On the other hand, the zero-subscripts in *k*0, *ω*0, and *ζ*01(*ρ*1; *t*) indicate that the radar transmitting frequency is considered constant during the radar operation. As explained in [16,18], the Sommerfeld attenuation function does not present significant variations with respect to radar frequency for typical bandwidths in a pulsed HF radar operation, allowing the frequency in *F*(*ρ*) to be considered constant and equal to the center-transmitting angular frequency *ω*0. Similarly, it is easy to verify through dimensional analysis that variations in the arbitrary height factor *ζ*1(*ρ*1; *t*) with respect to radar frequency are very small, allowing it to be redefined as *<sup>ζ</sup>*01(*ρ*1; *<sup>t</sup>*), where *<sup>k</sup>* <sup>=</sup> *<sup>k</sup>*<sup>0</sup> <sup>=</sup> *<sup>ω</sup>*<sup>0</sup> *c* .

In order to obtain the expression for the time-varying bistatic electric field over the ocean surface, the arbitrary height factor must be further addressed. From the definition in (7), and taking the first-order expression for the ocean surface expansion appearing in (8), it can be shown that

$$\zeta\_{01}(\boldsymbol{\rho}\_{1};t) = \sum\_{\mathbf{K'},\omega\_{\mathbf{K'}}} f\_{1}(\mathbf{K'},\omega\_{\mathbf{K'}}) \sqrt{\mathbf{K'}^{2} - k\_{0}^{2}} e^{-j\omega\_{\mathbf{K'}}t} e^{j\boldsymbol{\rho}\_{1}\cdot\mathbf{K'}}.\tag{13}$$

Substituting (13) into (12) and performing the convolution, the time-varying bistatic electric field for arbitrary heights is obtained:

$$\begin{split} \left( (\mathbf{E}\_{n}^{+})\_{12} \langle t, t\_{0} \rangle \sim \frac{-|\eta\_{0} \Lambda \ell \rangle \mathbf{s}\_{\mathbf{k}}^{2}}{(2\pi)^{1/2}} \sum\_{\mathbf{K}', \mathbf{k} \mathbf{k}'} \sum\_{\mathbf{K} \mathbf{k} \mathbf{k}} f\_{1}(\mathbf{K}, \omega\_{\mathbf{K}}) f\_{1}(\mathbf{K}', \omega\_{\mathbf{K}'}) \sqrt{\mathbf{K}'^{2} - \mathbf{k}\_{0}^{2}} \sqrt{\mathbf{K} \cos \phi\_{0}} e^{-|\omega\_{\mathbf{k}} \cdot \mathbf{t}'|} \delta(\omega\_{\mathbf{k}} - \omega\_{\mathbf{k}'}) \\ \cdot e^{|\rho\_{1} \cdot \mathbf{K}'} e^{|\omega\_{\mathbf{k}} \mathbf{t} \cdot \mathbf{t}} e^{|\mathbf{k} \mathbf{k} \mathbf{k} \rho\_{1}} e^{|\frac{\mathbf{K}}{2} \cdot \mathbf{K}} e^{|\rho\_{0} \cdot (\mathbf{K} \cos \phi\_{0})} \frac{F(\rho\_{0} \mathbf{u}, \omega\_{0}) F(\rho\_{0} \omega\_{0})}{\sqrt{\rho\_{0} \left[ e^{2}\_{0} - \left( \frac{\mathbf{s}}{2} \right)^{2} \right]}} e^{|\pi/4} \Delta \rho\_{\mathbf{s}} \operatorname{Sa} \left[ \frac{\Lambda \rho\_{1}}{2} \left( \frac{\mathbf{K}}{\cos \phi\_{0}} - 2k\_{0} \right) \right], \end{split} \tag{14}$$

where *ρ*0*s*, *ρ*01, and *ρ*<sup>02</sup> are functions of *t*0.

Here, a differentiation must be made between the two different time-related variables *t* and *t*0. In [16], Walsh and Gill differentiated between the "time of observation" *t*<sup>0</sup> and the "experiment time" *t* for successive pulses, while here, *t*<sup>0</sup> will be treated as the "radar time" and *t* as the "ocean surface time". Since radar-dependent events occur on a different time scale from that of the events of the ocean surface, they can be treated as two independent variables, even though both variables refer to time.

Now that (14) has been derived, it is possible to obtain the radar cross-section from the Fourier transform of the autocorrelation of the electric field, as shown in [18].

## **3. First-Order Radar Cross-Section of the Ocean Surface with Arbitrary Heights**

From [18], the autocorrelation for a general electric field, denoted as *E*<sup>+</sup> *<sup>n</sup>* , can be defined as

$$\mathcal{R}(\tau) = \frac{A\_r}{2\eta\_0} \mathbb{E}\left\{ E\_n^+(t) \overline{E\_n^+(t-\tau)} \right\},\tag{15}$$

such that R(0) coincides with the average power at the receiver. In (15), *Ar* is the effective free-space aperture of the receiver, *<sup>η</sup>*<sup>0</sup> is the intrinsic impedance of the free space, <sup>E</sup>{·} is the expected value operator, and the bar over the electric field indicates its conjugate. Expanding (15) into its different perturbation orders and knowing that, for the ocean surface [27]

$$\mathbb{E}\left\{f\_{\mathfrak{m}}(\mathbf{K},\omega\_{\mathbb{K}})\overline{f\_{\mathfrak{n}}(\mathbf{K}',\omega\_{\mathbb{K}'})}\right\} = \begin{cases} \mathbb{S}\_{\mathfrak{m}}(\mathbf{K},\omega\_{\mathbb{K}})d\mathbf{K}d\omega\_{\mathbb{K}'} \text{ if } \mathfrak{m} = \mathfrak{n}, \ \mathbf{K} = \mathbf{K}', \ \omega\chi = \omega\_{\mathbb{K}'}\\ 0, \text{ otherwise,} \end{cases}$$

where *fm*,*n*(·) are the *m*, *n*-th-order terms of the asymptotic expansion of the ocean surface and *Sm*,*n*(·) are their corresponding ocean wave spectra, it is easy to show that the only surviving terms of the autocorrelation are the ones multiplying fields of the same order:

$$\begin{split} R(\tau) &= \frac{A\_r}{2\eta\_0} \mathbb{E}\left\{ (E\_n^+)\_{11}(t) \overline{(E\_n^+)\_{11}(t-\tau)} \right\} + \frac{A\_r}{2\eta\_0} \mathbb{E}\left\{ (E\_n^+)\_{12\varepsilon}(t) \overline{(E\_n^+)\_{12\varepsilon}(t-\tau)} \right\} \\ &+ \frac{A\_r}{2\eta\_0} \mathbb{E}\left\{ (E\_n^+)\_{21}(t) \overline{(E\_n^+)\_{21}(t-\tau)} \right\} + \dotsb \\ &\equiv R\_{11}(\tau) + R\_{12\varepsilon}(\tau) + R\_{21}(\tau) + \dotsb, \end{split} \tag{16}$$

where

$$R\_{11}(\tau) = \frac{A\_r}{2\eta\_0} \mathbb{E}\left\{ (E\_n^+)\_{11}(t) \overline{(E\_n^+)\_{11}(t-\tau)} \right\},$$

$$R\_{12c}(\tau) = \frac{A\_r}{2\eta\_0} \mathbb{E}\left\{ (E\_n^+)\_{12c}(t) \overline{(E\_n^+)\_{12c}(t-\tau)} \right\},$$

and

$$R\_{21}(\tau) = \frac{A\_r}{2\eta\_0} \mathbb{E}\left\{ (E\_n^+)\_{21}(t) \overline{(E\_n^+)\_{21}(t-\tau)} \right\},$$

are respectively the autocorrelations of the first-order electric field, second-order correction to the first-order electric field, and second-order hydrodynamic electric field, with the first-order and second-order hydrodynamic electric fields defined as in [17]. Therefore, the correction term for the bistatic electric field for the ocean surface with arbitrary heights presented here does not affect the form of the radar cross-section expressions previously derived in [18]. Thus, taking the autocorrelation of (14) according to the expression presented in (16) and proceeding with the derivations, the following expression is obtained:

$$R\_{12\varepsilon}(\tau) = \frac{A\_{\tau}\pi^{2}\eta\_{0}|\Delta\ell\_{0}|^{2}k\_{0}^{4}|F(\rho\_{01}\omega\_{0})F(\rho\_{02}\omega\_{0})|^{2}}{2(2\pi)^{3}\rho\_{0\ast}\left[\rho\_{0\ast}^{2}-\left(\frac{\nu}{2}\right)^{2}\right]}\sum\_{m=-\pm 1}\int\int\mathcal{S}(m\mathbf{K})\mathcal{S}(m\mathbf{K})|K^{2}-k\_{0}^{2}|\varepsilon^{-j\omega\_{0}\tau}\frac{K^{\frac{\delta}{2}}}{\sqrt{\varepsilon}}\\ \tag{17}$$
 
$$\therefore \Delta\rho\_{s}^{2}\operatorname{Sa}\left[\frac{\Delta\rho\_{s}}{2}\left(\frac{\mathbf{K}}{\cos\phi\_{0}}-2k\_{0}\right)\right]dKd\theta\_{\mathbf{K}}.$$

From [18], it is known that

$$\frac{d\mathcal{P}(\omega\_d)}{dA} = \frac{A\_r \eta\_0 |\Delta \ell I\_0|^2 k\_0^2 |F(\rho\_{01}, \omega\_0) F(\rho\_{02}, \omega\_0)|^2}{16 (2\pi)^3 (\rho\_{01} \rho\_{02})^2} \sigma(\omega\_d),\tag{18}$$

where P(*ωd*) is the power spectral density of the electric field, defined as the Fourier transform of the autocorrelation with respect to *τ*, and *σ*(*ωd*) is the radar cross-section of the scattering object. After obtaining the power spectral density from the autocorrelation in (17), knowing from the bistatic scattering geometry that *θ***<sup>K</sup>** = *θN*, where *θ<sup>N</sup>* is the direction normal to the scattering ellipse at the scattering patch, and that [18]

$$\frac{\Delta\rho\_{\rm s}d\theta\_{N}}{\rho\_{0\rm s}\left[\rho\_{0\rm s}^{2}-\left(\frac{\rho}{2}\right)^{2}\right]}=\frac{dA}{(\rho\_{01}\rho\_{02})^{2}}\gamma$$

the second-order correction to the first-order bistatic radar cross-section for an ocean surface with arbitrary heights can be obtained by comparison with (18) as

$$\sigma\_{12c}(\omega\_d) \quad = \left. 2^5 \pi^3 k\_0^2 \Delta \rho \sum\_{m=\pm 1} \mathcal{S}(m\mathbf{K}) \mathcal{S}(m\mathbf{K}) \left| \boldsymbol{K}^2 - k\_0^2 \right| \frac{\boldsymbol{K}^4}{\mathcal{S}} \cos \phi\_0 \text{Sa} \left[ \frac{\Delta \rho\_5}{2} \left( \frac{\boldsymbol{K}}{\cos \phi\_0} - 2k\_0 \right) \right] . \tag{19}$$

## **4. Simulation Results**

In order to assess the impact of the newly-derived term in the total radar cross-section of the ocean surface, simulations were conducted in both high and low sea states. In both cases, two different bistatic configurations were chosen, with two different dominant wave directions, such that both the effects of wave directions and bistatic configurations on the results could be analyzed.

Before proceeding with the radar cross-section simulations, the validity of the second-order correction to the first-order term must be investigated, since the small-slope approximation still applies to (19). For this purpose, a number of total mean-square slope models proposed in the literature were used to compute the root-mean-square slope for the different ocean conditions used in the simulations. These models were developed empirically, using ocean surface measurements obtained with different instruments such as aerial photographs [33] and GPS-R [34]. The resulting total root-mean-square slopes for each of the simulated meteorological conditions are presented in Table 1, where *U*19.5 is the wind speed measured at 19.5 m above the ocean surface.

**Table 1.** Total root-mean-square slopes for simulated meteorological conditions using different slope models.


For the radar cross-section simulations, the Pierson–Moskowitz (PM) spectrum was chosen as the nondirectional spectral model of the ocean surface [38], with a cosine-power model for the directional factor [39] using the frequency-dependent wave-spreading factor proposed in [40]. Figures 2 and 3 present the radar cross-section simulation results for low and high sea states, respectively. In the presented results, *σ*<sup>11</sup> and *σ*2*<sup>P</sup>* respectively indicate the first- and second-order radar cross-sections of the ocean surface.

*θ*<sup>01</sup> = 15 deg.

**Figure 2.** Total bistatic radar cross-section and its components for three different bistatic geometries and wave directions at low sea states. Wind speed *U*19.5 = 10 m/s, radar frequency *f*<sup>0</sup> = 13.385 MHz, and roughness scale *k*0*Hs* = 0.60.

*θ*<sup>01</sup> = 60 deg.

*θ*<sup>01</sup> = 15 deg.

**(a)** *θ<sup>W</sup>* = 90 deg., *φ*<sup>0</sup> = 12.5 deg., *θ*<sup>01</sup> = 15 deg.

**(b)** *θ<sup>W</sup>* = 120 deg., *φ*<sup>0</sup> = 12.5 deg., *θ*<sup>01</sup> = 15 deg.

**(c)** *θ<sup>W</sup>* = 120 deg., *φ*<sup>0</sup> = 20 deg., *θ*<sup>01</sup> = 60 deg.

**Figure 3.** Total bistatic radar cross-section and its components for three different bistatic geometries and wave directions at high sea states. Wind speed *U*19.5 = 20.7 m/s, radar frequency *f*<sup>0</sup> = 13.385 MHz, and roughness scale *k*0*Hs* = 2.56.

#### **5. Discussion**

In [26], the transition zone for the validity of the asymptotic perturbation method for roughness scales was defined for perturbation parameters between 0.4 and 0.7. Using the same intervals for the perturbation parameter chosen in the present work, it can be observed that all of the empirical models shown in Table 1 yielded total root-mean-square slopes below the lower bound of the transition zone, meaning that the validity condition for the perturbation theory has not been violated under either of the meteorological conditions used in the simulations.

In observing the results in Figures 2 and 3, it is clear that the proposed cross-section has little to no impact on the total cross-section at low sea states, with a maximum difference of less than 0.1 dB, as shown in Figure 2b. This is an expected result, since the traditional second-order scattering theory is still valid for *k*0*Hs* < 0.7, even though the roughness scale is within the transition zone [26]. However, when observing the total radar cross-section at high sea states in Figure 3, where the roughness scale is above the upper limit of the transition zone, the correction term has an evident impact on the total radar cross-section, with a maximum difference of 8 dB in Figure 3c.

The enhanced part in the total cross-section depends on the dominant wave direction of the ocean surface *θW*, as observed when comparing Figure 3a with Figure 3b, and on the bistatic radar configuration, as evidenced in analyzing Figure 3b,c. The maximum difference between the total cross-sections with and without the additional term on the presented results is 8 dB.

In addition, the effects of the additional term are mostly evident in the central part of the total cross-section, at Doppler frequencies close to 0 rad/s; this is due to natural limitations on the steepness of large ocean waves, as well as to restrictions on the wave slope that are still imposed in the present analysis.

#### **6. Conclusions and Future Work**

In the present work, the scattered electric field and radar cross-section for an ocean surface with arbitrary heights using a narrow-beam bistatic HF radar have been derived. Previously derived electric field expressions presented in [17] for the small-height condition still appear in the final result, but due to the removal of the small-height approximation, a new term appears. The new term is interpreted as a correction to the first-order cross-section for arbitrary heights. The radar cross-section due to the correction term is also derived following the procedure presented in [18], and is then simulated for low and high sea states, showing an impact on the total radar cross-section for high sea states. A similar analysis is currently under review for the monostatic case.

The present work shows that changes in the dominant wave direction of the ocean surface and the bistatic configuration impact the contribution of the correction term to the radar cross-section. Since the other terms of the radar cross-section are saturated at sufficiently high sea states, the correction term may provide input for determining the dominant wave direction at high sea states. More extensive work needs to be dedicated to the understanding of interactions between the bistatic configuration and dominant wave direction, as well as the correction term, determining thresholds above which the correction term should be considered, and how different bistatic configurations can be used in observing the ocean surface at large roughness scales.

It can also be noted that the simulations in the current work were carried out using a Pierson–Moskowitz spectral model, which implies the assumption of a fully-developed sea and does not depend on fetch [38]. Since changes in fetch are known to affect the mean squared slope of the ocean surface [27], the work of determining how changes in fetch affect the correction term, allowing its application to developing sea states, is ongoing.

Since the present work represents the first attempt to overcome the small-height constraint imposed on previously derived bistatic high-frequency radar cross-sections of the ocean surface, no practical applications of the proposed theory have been suggested in the present work, as more work must be done in the future for these applications to be devised. In addition, due to the lack of narrow-beam bistatic HF radar data available to the authors, which were measured under conditions

that violate the limiting roughness scales of the theories proposed in the literature, or the lack of any other derived model for radar cross-section of the ocean surface at large roughness scales, a validation of the presented theory is not available in the current work. To validate the current theory, collaboration from the radio oceanography community is necessary, as multiple observations need to be conducted by users of narrow-beam bistatic HF radars at ocean conditions that violate the small-height assumption. In addition, efforts to establish bistatic operation of existing monostatic HF radars are ongoing on the coastline of Placentia Bay, Newfoundland, Canada in order to collect data for the validation of the presented theory.

**Author Contributions:** Conceptualization: M.T.S., W.H., and E.W.G.; methodology: M.T.S.; formal analysis: M.T.S.; investigation: M.T.S.; resources: W.H. and E.W.G.; data curation: M.T.S.; writing—original draft preparation: M.T.S.; writing—review and technical editing: M.T.S., W.H., and E.W.G.; visualization: M.T.S.; project administration: W.H.; supervision: W.H. and E.W.G.; funding acquisition: W.H. and E.W.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the Natural Sciences and Engineering Research Council of Canada (NSERC) under Discovery Grants to W. Huang (NSERC RGPIN-2017-04508 and RGPAS-2017-507962) and E. W. Gill (NSERC RGPIN-2015-05289).

**Conflicts of Interest:** The authors declare no conflict of interest.
