# Sea Surface Radar Scattering at L-Band Based on Numerical Solution of Maxwell's Equations in 3-D (NMM3D)

Tai Qiao, Leung Tsang, Douglas Vandemark, Simon H. Yueh, Tien-Hao Liao, Frederic Nouguier, Bertrand Chapron

► To cite this version:

Tai Qiao, Leung Tsang, Douglas Vandemark, Simon H. Yueh, Tien-Hao Liao, et al. Sea Surface Radar Scattering at L-Band Based on Numerical Solution of Maxwell's Equations in 3-D (NMM3D). IEEE Transactions on Geoscience and Remote Sensing, 2018, 56 (6), pp.3137-3147. 10.1109/TGRS.2018.2792432. hal-04202139

HAL Id: hal-04202139

https://hal.science/hal-04202139v1

Submitted on 26 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Sea Surface Radar Scattering at L-Band Based on Numerical Solution of Maxwell's Equations in 3-D (NMM3D)

Tai Qiao<sup>D</sup>, Student Member, IEEE, Leung Tsang, Fellow, IEEE, Douglas Vandemark, Simon H. Yueh, Fellow, IEEE, Tien-Hao Liao, Frédéric Nouguier<sup>D</sup>, and Bertrand Chapron

**Abstract**— Radar scattering from ocean surfaces is investigated by 3-D numerical solution of Maxwell's equations [numerical Maxwell's model in 3-D (NMM3D)] using the ocean surface profiles stochastically generated from a 3-D Durden–Vesecky ocean spectrum. The surface integral equations (SIEs) are formulated for dielectric surfaces using Green's functions of the air and the ocean permittivities with the surface tangential electric and magnetic fields as the unknowns. In solving the SIEs using the method of moment, a fast matrix solver of the sparse matrix canonical grid is used in conjunction with Rao–Wilton–Glisson basis functions. The computation has been implemented on a high-performance parallel computing cluster for problems with up to six million surface unknowns. Unlike the two-scale model (TSM) approximation, NMM3D does not require division of the surface spectrum into large- and small-scale ocean waves. The results of backscattering simulations are compared to Aquarius satellite radar measurements for wind speeds of 5, 8, and 10 m/s and for incidence angles of 29°, 39°, and 46°. The results show that NMM3D ocean backscattering solutions at L-band are in good agreement with Aquarius satellite radar data for co-polarized VV, HH, and cross-polarized VH returns as well as for the VV/HH ratio. The azimuthal dependence of L-band backscatter is also assessed. Finally, NMM3D results are compared to TSM solutions and are shown to lie close to Aquarius data in observed VV/HH ratio, and their azimuthal dependencies.

**Index Terms**—Aquarius, Durden–Vesecky (DV) spectrum, numerical simulations, ocean remote sensing, radar scattering.

## I. INTRODUCTION

PHYSICAL models of electromagnetic scattering from random rough surfaces have broad application in microwave remote sensing of land surfaces and ocean surfaces.

Manuscript received March 4, 2017; revised August 7, 2017 and November 14, 2017; accepted December 11, 2017. This work was supported by the NASA Remote Sensing Theory Program under Project NNX15AQ92G. (Corresponding author: Tai Qiao.)

T. Qiao and L. Tsang are with the Radiation Laboratory, Department of Electrical Engineering and Computer Science, University of Michigan, Ann Arbor, MI 48109 USA (e-mail: qiaot@umich.edu).

D. Vandemark is with the Institute for the Study of Earth, Oceans and Space, University of New Hampshire, Durham, NH 03824 USA.

S. H. Yueh and T.-H. Liao are with the Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA 91125 USA.

F. Nouguier and B. Chapron are with the Laboratoire d'Océanographie Spatiale, Institut Français de Recherche pour l'Exploitation de la Mer, 29280 Plouzané, France.

Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org.

Digital Object Identifier 10.1109/TGRS.2018.2792432

Analytical approximate methods include the small perturbation method (SPM), geometrical optics (GO), physical optics (PO), advanced integral equation method, the small-slope approximation, and the composite surface model [1]–[6].

Ocean surfaces have a range of scales, with gravity waves of order of meters and short gravity and gravity-capillary waves with scales in the range of centimeters to millimeters. A common scattering model is the composite surface model, also known as the two-scale model (TSM) [7], [8]. In this approach, the surface height spectrum is divided into two in order to separate large- and small-scale roughness elements. PO/GO and SPM are applied, respectively, to the two scales, and the total radar cross section is obtained by summing the two contributions. The dividing line is usually chosen to be in the range of 1/5 to a half of the free-space wavenumber of the incident microwave [9]–[11]. The TSM is a heuristic model as the decomposition of an object into parts followed by coherently or incoherently combining the scattering of each does not follow from Maxwell's equations; otherwise, Maxwell's equations can be readily solved for any large object. Thus, TSM requires tuning of dividing line between large scales and small scales. The optimal choice depends on incidence angles [11].

In recent years, numerical simulations of 3-D rough surface scattering have become feasible by combining fast numerical methods and high-performance computing [12]–[16]. The simulations have been applied to land and ocean surfaces. Both land and sea surfaces have fine-scale (millimeter to meter) features that can cause significant microwave radar backscattering. For example, gravity-capillary and capillary waves have small radii of curvature compared with microwave wavelength that can lead to strong backscattering and cross-polarized returns. Because of the fine-scale surface features, the surface integral equation (SIE) approach, with discretizations of surfaces, is more accurate than using the finite-difference time-domain or finite-element method approaches. The latter two methods require volumetric discretization. The SIE approach combines fine-scale surface discretization with the method of moments (MoM) to capture the impact of all relevant surface structure on the scattered electromagnetic radiation. Examples of fast dense matrix solvers include the sparse matrix canonical grid method (SMCG) [15] and the multilevel fast multiple algorithm [16]. In order to reduce
the number of iterations, a near-field preconditioner has been used [13]. Some past work in numerical simulation treated the ocean surface as a perfect electric conductor [14]. Several past studies also implemented impedance boundary conditions (IBCs) [16], [17]. The IBC states that the tangential electric field is equal to the product of the wave impedance of the ocean permittivity and the tangential magnetic field. The approximation is valid if the radius of curvature of the surface is large compared with the microwave wavelength. Compared to previous works which are making approximations in the modeling process, in this paper, Maxwell's equations are solved with a finite permittivity. Some earlier works are simulating 1-D surfaces or 2-D surfaces with pulse basis functions, while in this paper, 2-D surfaces with Rao-Wilton-Glisson (RWG) basis functions are being studied [18]–[23].

In this paper, we apply numerical Maxwell's model in 3-D (NMM3D) [12] to simulate L-band radar scattering from a wind-roughened ocean. The SIEs of Poggio-Miller-Chang-Harrington-Wu-Tsai (PMCHWT) [13] are used. The ocean surface profiles are stochastically generated using 3-D ocean spectra. In PMCHWT formulation, two separate Green's functions are used, one for air permittivity and another for ocean permittivity. The SMCG [15] is used in conjunction with RWG basis functions. For modeling of the rough surface scattering, the physical problem is 3-D vector wave scattering. However, the rough surface boundary has 2 dimensions with 2 degrees of freedom. The SMCG method uses a Taylor expansion of Green's function about the flat surface so that 2-D fast Fourier transform (FFT) can be applied to the 3-D electromagnetic scattering problem. The SMCG method has been implemented in a parallel computation framework [24]. The scattering problems in this paper include cases with up to six million surface unknowns in the tangential electric and magnetic fields. In NMM3D simulations, a finite surface size is used and the surfaces are discretized. This means that the ocean spectrum is truncated between $k_l$ and $k_u$, where $k_l$ is inversely proportional to the surface size while $k_u$ is inversely proportional to the discretizing sampling. Note that in NMM3D simulations, there is no need to impose a surface scale separation within this domain of $k_l$ to $k_u$. By changing $k_l$ or surface size, one can investigate the influence of adding roughness due to longer ocean surfaces. By varying $k_u$, one can assess the impact of small-scale capillary waves. This paper will focus on L-band ocean backscattering simulations and their comparison to results obtained from the Aquarius satellite radar system for wind speeds of 5, 8, and 10 m/s. Aquarius measurements provide estimates for VV, HH, VH, and the VV/HH ratio, and NMM3D backscattering results are produced for each case. Comparisons will also be made with TSM backscatter estimates to evaluate the efficacy of this widely used approximation method.

This paper is organized as follows. In Section II, we describe the ocean spectrum used in this paper. In Section III, we describe the formulations of SIE and the numerical implementations. In Section IV, the influences of $k_l$ and $k_u$, the lower and upper limits of truncated ocean spectrum on backscatter, are illustrated. In Section V, we illustrate numerical results and make comparisons to data and the TSM model.

## II. OCEAN SPECTRUM AND SURFACE GENERATION

Ocean roughness spectra have been extensively used in electromagnetic modeling of ocean surface scattering. Common ocean spectra include the Pierson-Moskowitz, Durden-Vesecky (DV) [16], [26], Donelan-Banner-Jahne, and Elfouhaily spectra [25], [27]–[29]. These wind-dependent ocean wave height models are written in the form of the product of an isotropic part and an azimuthal anisotropic (wave direction) part

$$W(k, \phi, u) = \frac{1}{k} S(k, u) \Phi(k, \phi, u) \quad (1)$$

$$\Phi(k, \phi, u) = 1 + \Delta(k, u) \cos(2\phi) \quad (2)$$

where $S(k, u)$ is the omnidirectional spectrum with $k$ being the spectral components and $u$ being the surface wind speed. $\phi$ is the azimuthal angle with respect to the wind direction. $\Phi(k, \phi, u)$ is the directional spreading function or angular dependence function which denotes the nonuniform distribution of wave roughness in azimuth. In this paper, we used the DV spectrum for the NMM3D simulations. Appendix B gives the expression of the DV spectrum. Two versions are used. One version is the isotropic, or omnidirectional, spectrum defined by setting $\Delta(k, u) = 0$, and is referred to as isotropic DV in this paper. The other version used is the full directional DV spectrum, referred to as the anisotropic DV spectrum. Although other spectrum can be used such as the Elfouhaily spectrum, the purpose of this paper is to use NMM3D simulations to compare the results of TSM and SPM. Thus, we limit ourselves to the isotropic and anisotropic DV spectra [26].

In this paper, the microwave frequency is 1.26 GHz, which is the frequency of the Aquarius L-band scatterometer. The wavenumber of the microwave $k_0$ is $26 \text{ m}^{-1}$. According to the Bragg scattering theory [29], [30], the spectral component of the ocean spectrum that contributes to the backward scattering is

$$k_B = 2k_0 \sin \theta \quad (3)$$

where for the incidence angle ($\theta$) of 39°, the Bragg scattering wavenumber is $33 \text{ m}^{-1}$.

In the rough surface simulations computed using DV as an input for the SIE approach [15], we specify the surface area to be of an extent $L_x$ and $L_y$ in x- and y-directions, respectively. We also set surface lengths $L_x = L_y = L$. In the discretization of the surfaces, sampling intervals of $\Delta x = \Delta y$ are used. These limits set surfaces that contain spectral components $W(k, \phi, u)$ in the range of $k_l \le k \le k_u$, where

$$k_u = \frac{\pi}{\Delta x} \quad (4)$$

$$k_l = \frac{2\pi}{L} \quad (5)$$

In Fig. 1, we show a log-log scale of the isotropic DV spectrum for wind speeds of 5, 8, and 10 m/s. The wavenumber of the incident microwave, the Bragg wavenumber, $k_l, k_u$ for several surface lengths, and sampling are also shown. In Fig. 1, the lower limit $k_l$ corresponds to a surface by $64\lambda \times 64\lambda$ and the upper limit $k_u$ corresponds to $\Delta x = \Delta y = 1/32\lambda$.
<figure><img src="image_2.png" /><figcaption>Fig. 1. Isotropic DV surface elevation spectrum (solid lines) for wind speeds of 5, 8, and 10 m/s. Magenta and blue lines denote lower and upper limits. The lower limit $k_l$ corresponds to a surface size of $64\lambda \times 64\lambda$. The upper limit $k_u$ corresponds to $\Delta x = \Delta y = 1/32\lambda$. Red and green lines denote the incident microwave wavenumber and Bragg wavenumber, respectively.</figcaption></figure>

<table><caption>TABLE I
UPPER AND LOWER SPECTRAL LIMITS FOR VARYING SAMPLING POINT DENSITIES AND SURFACE SIZES</caption><thead><tr><th>Surface</th><th>Δx, Δy (λ)</th><th>k<sub>l</sub> (m<sup>-1</sup>)</th><th>k<sub>u</sub> (m<sup>-1</sup>)</th></tr></thead><tbody><tr><td>Length (λ)</td><td></td><td></td><td></td></tr><tr><td>16</td><td>1/16</td><td>1.65</td><td>211</td></tr><tr><td>32</td><td>1/16</td><td>0.82</td><td>211</td></tr><tr><td>16</td><td>1/32</td><td>1.65</td><td>422</td></tr><tr><td>64</td><td>1/16</td><td>0.41</td><td>211</td></tr></tbody></table>

The incident L-band and Bragg wavenumbers lie between the upper limit and lower spectral limits.

In Table I, we tabulate $k_l$ and $k_u$ values corresponding to the surface sizes and the discretization intervals used in this paper. With the increase of surface length, large-scale ocean waves are incrementally increased in the simulations. The peak wavenumbers for 5, 8, and 10 m/s are around 0.27, 0.1, and 0.07 m⁻¹. Because of the finite length in these surface simulations, longer ocean wave scales (especially at higher wind speeds) and their tilt impacts are not be included in the results.

### III. PMCHWT SURFACE INTEGRAL EQUATION AND NUMERICAL IMPLEMENTATION USING RWG AND SMCG

In this section, we summarize the methodology of NMM3D simulations [13]–[15].

For scattering by dielectric surfaces, the PMCHWT formulations of SIEs are given in the following form:

$$ \begin{bmatrix} -T - \frac{\eta_1}{\eta} T_1 & K + K_1 \\ -K - K_1 & -T - \frac{\eta_1}{\eta} T_1 \end{bmatrix} \begin{bmatrix} \tilde{J} \\ M \end{bmatrix} = \begin{bmatrix} \hat{n} \times \bar{E}^{\text{inc}} \\ \hat{n} \times \bar{H}^{\text{inc}} \end{bmatrix} \quad (6) $$

where

$$ T_i(\mathbf{X}) = ik_i \hat{n} \times \int_S \left( \bar{I} + \frac{\nabla\nabla}{k_i^2} g_i(\bar{r}, \bar{r}') \right) \cdot \mathbf{X}(\bar{r}') d\bar{r}' \quad (7) $$

$$ K_i(\mathbf{X}) = \hat{n} \times \text{P.V.} \int_S \nabla g_i(\bar{r}, \bar{r}') \times \mathbf{X}(\bar{r}') d\bar{r}' \quad (8) $$

and

$$ g_i(\bar{r}, \bar{r}') = \frac{e^{ik_i |\bar{r} - \bar{r}'|}}{4\pi |\bar{r} - \bar{r}'|} \quad (9) $$

$$ \nabla g_i(\bar{r}, \bar{r}') = \frac{(-1 + ik_i |\bar{r} - \bar{r}'|) \exp(ik_i |\bar{r} - \bar{r}'|)}{4\pi |\bar{r} - \bar{r}'|} (\bar{r} - \bar{r}'). \quad (10) $$

In (9) and (10), $k$ and $g(\bar{r}\bar{r}')$ are the wavenumber and Green's function in air, respectively, while $k_1$ and $g_1(\bar{r}\bar{r}')$ are the wavenumber and Green's function for ocean permittivity. Next, we apply the MoM to the above SIEs. RWG basis functions are used. Curl-conforming RWG functions are used as testing functions. A property of the discretized equations is that the impedance matrix elements are symmetric, e.g.

$$ Z_{mn} = Z_{nm} \quad (11) $$

We have used this property to reduce CPU memory and computation time.

A finite surface area of $L$ by $L$ is used in the SIEs. To eliminate the edge effects for a surface with finite size, a tapered plane wave [31], in the form of 2-D Fourier transform, is used

$$ \bar{E}_i(\bar{r}) = \int_{-\infty}^{\infty} d\bar{k}_{\rho} e^{i(\bar{k}_{\rho} \cdot \bar{r} - k_z z)} \psi(\bar{k}_{\rho}) \bar{e}(\bar{k}_{\rho}) \quad (12) $$

$$ \bar{H}_i(\bar{r}) = \int_{-\infty}^{\infty} d\bar{k}_{\rho} e^{i(\bar{k}_{\rho} \cdot \bar{r} - k_z z)} \frac{\psi(\bar{k}_{\rho})}{\eta} \bar{h}(\bar{k}_{\rho}). \quad (13) $$

For fast matrix vector multiplication in the iterative solution of the matrix equation, we used the SMCG method. In the SMCG, a Taylor series expansion of Green's function is applied [16] about a flat surface defined by the xy plane. Subsequently, a 2-D FFT is applied. The Taylor series expansions are performed for Green's functions and the gradient of Green's functions

$$ g_i(\bar{r}, \bar{r}') = \frac{e^{ik_i |\bar{r} - \bar{r}'|}}{4\pi |\bar{r} - \bar{r}'|} = \sum_{m=0}^{M} a_m^{(1,2)}(\rho_R) \left(\frac{z_d^2}{\rho_R^2}\right)^2 \quad (14) $$

$$ \nabla g_i(\bar{r}, \bar{r}') = \frac{(-1 + ik_i |\bar{r} - \bar{r}'|) \exp(ik_i |\bar{r} - \bar{r}'|)}{4\pi |\bar{r} - \bar{r}'|} (\bar{r} - \bar{r}') = \sum_{m=0}^{M} b_m^{(1,2)}(\rho_R) \left(\frac{z_d^2}{\rho_R^2}\right)^2 \quad (15) $$

where $z_d = f(x, y) - f(x', y')$. The coefficients $a_m^{(1,2)}(\rho_R)$ and $b_m^{(1,2)}(\rho_R)$ are translationally invariant in the horizontal
directions. We apply 2-D FFT to the coefficients, and these are performed for the product of the Toeplitz matrix and column vector [16]. The first four coefficients of (14) and (15), $a_m$ and $b_m$, are listed in Appendix A. In Appendix A, we also describe improvements to use SMCG on NMM3D compared to that used in previous implementations.

To implement the SMCG in conjunction with RWG basis function into the surface integral equations, the discretized surface equation must be written in a form that the field points and sources points are separated apart. The details are illustrated in Appendix A.

In this paper, the Taylor series are expanded to fifth order. We have checked that five terms in the Taylor series expansion could give accurate results compared to direct Green's function multiplication such that the error norm is less than 0.2% for a wind speed of 5 m/s. The details of numerical implementation of SMCG on RWG basis functions are also described in Appendix A.

## IV. OCEAN SURFACE PROFILE GENERATIONS AND L-BAND BACKSCATTERING VARIATIONS WITH $k_l$ AND $k_u$

For L-band NMM3D simulations, we account for roughness due to both short gravity waves (5 m and less) [32] and shorter gravity-capillary waves (down to 2 cm or so) [33]. Since the wavelength of wind-generated gravity waves will increase with wind forcing, a larger sea surface for the NMM3D analysis is necessary for greater wind speeds. To account for the influence of fine capillary waves on backscatter, a small discretization of the surface is necessary. By varying surface length and discretizations as shown in Table I, we examine the effects of longer gravity waves and shorter gravity-capillary waves on backscattering at L-band. At a wind speed of 5 m/s, Aquarius data show the negative upwind-crosswind asymmetry [33], [37]. This means that upwind backscatter is smaller than crosswind backscatter. At other wind speeds (>8 m/s), upwind backscatter is stronger than crosswind backscatter. This asymmetry cannot be reproduced based on the DV anisotropic spectrum or using other existing spectrum. We will investigate the spreading function in the future. In the meantime, we use the isotropic DV spectrum for NMM3D for 5 m/s, but upwind data for Aquarius. For a wind speed of 8 m/s, the Aquarius data are weakly dependent on wind direction, especially for beam 1 and beam 2. Then, we used isotropic surfaces for 8 m/s. For 10 m/s, we used anisotropic surfaces. We have also calculated surface sample statistics from randomly generated surface profiles [34]–[36]. The rms heights are provided in Table II for various surface lengths and wind speeds. We have made comparisons between numerical computation and analytical integration and found that the ocean statistics within the same spectrum range are the same between the two methods.

For a wind speed of 8 m/s, the increase of surface length from 16 wavelengths to 64 wavelengths, corresponding to the decrease of $k_l$, results in an increase of rms height from 5.1 to 13.8 cm. This indicates inclusion of longer surface waves with increased surface size.

In order to study the influence of the upper and lower limits of the spectral range used in the simulation, we conducted

<table><caption>TABLE II
RMS FOR ISOTROPIC DV SPECTRUM AT DIFFERENT WIND SPEEDS IN TERMS OF LONGEST SURFACE WAVELENGTH</caption><thead><tr><th>Ocean wind speed (m/s)</th><th>Surface Length L(λ)</th><th>k<sub>l</sub> (m<sup>-1</sup>)</th><th>rms in centimeters</th></tr></thead><tbody><tr><td>5</td><td>16</td><td>1.65</td><td>4.8</td></tr><tr><td>5</td><td>32</td><td>0.82</td><td>5.5</td></tr><tr><td>5</td><td>64</td><td>0.41</td><td>5.5</td></tr><tr><td>8</td><td>16</td><td>1.65</td><td>5.1</td></tr><tr><td>8</td><td>32</td><td>0.82</td><td>10</td></tr><tr><td>8</td><td>64</td><td>0.41</td><td>13.8</td></tr></tbody></table>

<table><caption>TABLE III
VV AND HH FOR VARYING DISCRETIZING DENSITIES AT WIND = 5 m/s WITH 39° INCIDENCE ANGLE AND ISOTROPIC DV SPECTRUM IN COMPARISON WITH UPWIND AQUARIUS DATA</caption><thead><tr><th>Δx, Δy (λ)</th><th>k<sub>u</sub> (m<sup>-1</sup>)</th><th>VV (dB)</th><th>HH (dB)</th><th># of realizations</th></tr></thead><tbody><tr><td>1/16</td><td>211</td><td>-15.42</td><td>-20.45</td><td>30</td></tr><tr><td>1/32</td><td>422</td><td>-15.49</td><td>-20.27</td><td>30</td></tr><tr><td>Aquarius</td><td>N/A</td><td>-15.50</td><td>-20.00</td><td>N/A</td></tr></tbody></table>

NMM3D simulations using different surface lengths and discretizations.

Table III shows the backscatter results for wind speed = 5 m/s at several sampling densities. The results show that 16 points per wavelength and 32 points per wavelength yield very similar backscattering results. The difference is within 0.3 dB. The Aquarius upwind data are also listed in Table III. Based on extensive simulations, sampling at 16 points per wavelength yields sufficient accuracy to account for the small-scale roughness. Including smaller capillary wave roughness does not alter the backscattering results for this case. For emissivity calculations, denser sampling is required. And results for cross-wind and upwind are compared for the anisotropic DV spectrum and are shown in Fig. 5.

Next, we conducted numerical simulations for several surface sizes at wind speeds of 5 and 8 m/s. These correspond to changing $k_l$. The backscatter results are shown in Table IV. We also include in Tables III and IV the number of realizations used in the simulations as shown in the last column of Tables III and IV. We have conducted convergence test with the number of realizations. The results show that 30 realizations can give stable results for backscatter. For the $16λ$ case, the computation time is about 3–4 min for each realization. For the $32λ$ case, it is about 20 min for each realization. For the
<figure><img src="image_3.png" /><figcaption>Fig. 2. Comparison of L-band NMM3D scattering amplitude distributions for 300 realizations for co-polarized cases (Left) VV and (Right) HH and for a wind speed of 5 m/s using the isotropic DV spectrum and an incidence angle 29°. The red curve indicates that for a Rayleigh-distributed process.</figcaption></figure>

<table><caption>TABLE IV
VV AND HH FOR L-BAND 39° INCIDENCE ISOTROPIC DV SPECTRUM
COMPARED WITH UPWIND AQUARIUS DATA. (a) VV AND HH FOR
VARYING SURFACE SIZES AT WIND = 5 m/s. (b) VV AND HH
FOR VARYING SURFACE SIZES AT WIND = 8 m/s
(a)</caption><thead><tr><th>Surface Size</th><th>k<sub>l</sub> (m<sup>-1</sup>)</th><th>VV (dB)</th><th>HH (dB)</th><th># of realizations</th></tr></thead><tbody><tr><td>16 × 16λ<sup>2</sup></td><td>1.65</td><td>-15.42</td><td>-20.45</td><td>30</td></tr><tr><td>32 × 32λ<sup>2</sup></td><td>0.82</td><td>-15.54</td><td>-20.71</td><td>30</td></tr><tr><td>64 × 64λ<sup>2</sup></td><td>0.41</td><td>-15.50</td><td>-20.50</td><td>30</td></tr><tr><td>Aquarius</td><td>N/A</td><td>-15.50</td><td>-20.00</td><td>N/A</td></tr></tbody></table>

<table><thead><tr><th colspan="5" style="text-align: center;">(b)</th></tr><tr><th>Surface Size</th><th>k<sub>l</sub> (m<sup>-1</sup>)</th><th>VV (dB)</th><th>HH (dB)</th><th># of realizations</th></tr></thead><tbody><tr><td>16 × 16λ<sup>2</sup></td><td>1.65</td><td>-14.72</td><td>-19.75</td><td>30</td></tr><tr><td>32 × 32λ<sup>2</sup></td><td>0.82</td><td>-14.65</td><td>-19.51</td><td>30</td></tr><tr><td>64 × 64λ<sup>2</sup></td><td>0.41</td><td>-14.06</td><td>-19.34</td><td>29</td></tr><tr><td>Aquarius</td><td>N/A</td><td>-14.80</td><td>-18.65</td><td>N/A</td></tr></tbody></table>

64λ case, it is about 1 h for each realization. The convergence rate of iterations is influenced by the roughness of the ocean surface.

The NMM3D backscatters of various simulated surface sizes are shown in Table IV(a) and (b) for wind speeds of 5 and 8 m/s, respectively. It is noted that there is a less than 1-dB difference for VV and a 1.4-dB difference for HH at 5 and 8 m/s. One observes that for the two wind speeds, the maximum difference for backscatter for varying surface sizes is 0.7 dB. There is also an increase of backscatter with increasing wind speed from Table IV(a) to (b). It is to be noted from Table IV(b) that the three surface sizes at a wind

speed of 8 m/s have approximately the same backscatter. But in Table II, the three surface sizes have very different rms height for the same wind speed of 8 m/s. With the increase of simulated surface size, longer waves (waves that are meters long and have a large rms height) are included in the simulations. However, these longer waves do not cause significant differences on backscatter at a wind speed of 8 m/s. The backscatter arises from the spectral components in the ocean spectrum that are comparable to the wavelength of the incoming microwaves. For L-band, the decimeter and centimeters ocean waves cause backscatter. The above convergence analysis may also be influenced by the random errors inherent to the Monte Carlo simulation approach. In future studies, the use of larger surface sizes will be studied.

## V. NMM3D RESULTS AND DISCUSSION

In this section, we present NMM3D numerical ocean scattering results at L-band and their comparison with Aquarius data. The assumed surface relative permittivity is 75+61i, corresponding to a complex index of refraction of 9.26 + 3.29i. This is the case of seawater with a salinity of 35 parts per thousand at a water temperature of 10 °C at 1.26 GHz. For each realization, a single sea surface is generated, and Maxwell's equations are solved numerically. The scattering coefficients represent results averaged over 30–300 realizations.

Isotropic and anisotropic DV ocean spectra are used for simulations of wind-dependent L-band backscatter. The isotropic DV spectrum is used for 5 and 8 m/s. The anisotropic DV spectrum is used for 10 m/s [38], [40].

Fig. 2 shows the statistical distributions of the simulated amplitude of HH and VV computed using 300 realizations with a wind speed of 5 m/s and an incidence angle of 29°. The histogram of computed backscatter values is consistent with that shown for a Rayleigh distribution.

Fig. 3 shows the in-plane bistatic L-band scattering coefficients at 5-m/s wind speed and 39° incidence. In these simulations, we decompose the bistatic scattering coefficients into coherent and incoherent waves [13], [15]. For the surface sizes used, the coherent waves are only significant in the specular direction while scattering in other directions arises from the incoherent waves [13]–[15], [39], [40]. The backscattering direction is highlighted with a dashed line in Fig. 3.
<figure><img src="image_4.png" /><figcaption>Fig. 3. In-plane bistatic VV/HH for wind = 5 m/s isotropic DV spectrum.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. L-band backscatter at a wind speed of 5 m/s for different incidence angles: 29°, 39°, and 46° and isotropic DV spectrum compared with upwind Aquarius data.</figcaption></figure>

Fig. 4 shows the backscatter dependence at the three incidence angles of the Aquarius scatterometer 29°, 39°, and 46° at a wind speed of 5 m/s. Backscattering decreases with incidence angle. The upwind Aquarius data are also shown. For co-polarization, the agreement between simulations and upwind Aquarius data are within 1 dB at the three incidence angles.

In Fig. 5, we illustrate the backscatter for a wind speed of 10 m/s using an anisotropic DV spectrum. The incidence angle is 46°. There are four curves in Fig. 5: NMM3D with error bars derived from the Monte Carlo simulations, first-order SPM, Aquarius data, and composite surface model (TSM). There are several TSM models because TSM is a heuristic approximation, subdividing the objects into parts and then performing incoherent combining with tilts. We have referred to several TSM models in our computation. The uncertainty of numerical simulation is defined as

$$ \text{error (dB)} = 10 \log 10 \left( 1 + \frac{1}{\sqrt{N}} \right) \quad (16) $$

for VV and HH. In (16), N is the number of realizations. The four panels in Fig. 5 correspond to VV, HH, VV/HH (polarization ratio), VH, respectively. We can see that HH of SPM is lower than Aquarius data, while NMM3D, TSM, and SPM show VV results that are all comparable with Aquarius observations. The polarization ratio, VV/HH, derived using SPM is 1.5 dB larger than the data. NMM3D and

<figure><img src="image_6.png" /><figcaption>Fig. 5. L-band backscatter for wind speed = 10 m/s and incidence angle = 46°. (Top) VV and HH. (Bottom) VV/HH. TSM, NMM3D, and SPM results use the anisotropic DV spectrum.</figcaption></figure>

TSM for $k_d = (k_0/2)$ can both provide polarization ratio in reasonable agreement with data, while the polarization ratio given by SPM is larger than data. For the range of wind speeds (wind speed <10 m/s), we are studying, TSM and NMM3D both agree well with cross-polarization data. Fig. 5 also indicates that Aquarius data show a directional dependence in the VV/HH polarization ratio. The wind direction dependence of this ratio is also exhibited in NMM3D, while the ratio is direction invariant for SPM. With a dividing wavenumber equals to $(k_0/2)$, we see that TSM can predict scattering coefficients that are as consistently close the data as for NMM3D.

Fig. 6 provides the backscattering polarization ratio for wind speeds of 8 and 10 m/s and for the three incidence Aquarius angles [40]. NMM3D simulation uncertainties are also shown. For a wind speed of 8 m/s, the surfaces used in the computations are derived from the isotropic DV spectrum. The results of TSM and data shown in Fig. 6 are averaged over all wind directions. For wind = 10 m/s, the anisotropic DV spectrum is used in NMM3D. The results of NMM3D, TSM, and Aquarius data are also averaged over all wind directions. The results indicate that the VV/HH ratio increases with increasing incidence angle. From Fig. 6, we see that the difference for both TSM and NMM3D compared to Aquarius data is both within 0.6 dB [42].

In Fig. 7, we plotted the polarization ratio for a wind speed of 10 m/s and at an incidence angle of 46° for a range of cutoff wavenumbers, $(k_0/2)$, $(k_0/3)$, and $(k_0/5)$. Vertical axis is VV/HH ratio in dB. Each data point is obtained by averaging over all wind directions. It is apparent that the polarization ratio computed using the TSM varies with cutoff wavenumbers. For $k_d$ in the range of $(k_0/5)$ to $(k_0/2)$, the VV/HH ratio can vary by about 0.4 dB, indicating the sensitivity of TSM predictions within the nominal range of cutoff wavenumber selection.
<figure><img src="image_7.png" /><figcaption>Fig. 6. Comparison of NMM3D and TSM of VV/HH with Aquarius at 29°, 39°, and 46°. NMM3D is using isotropic DV spectrum for 8 m/s and anisotropic DV spectrum for 10 m/s. TSM with kd = k₀/2 is derived using the anisotropic DV spectrum for wind speeds of 8 and 10 m/s. Each data point is obtained by averaging over all wind directions.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7. Polarization ratio for a wind speed of 10 m/s and an incidence angle of 46°. The horizontal axis is different cutoff wavenumbers. Vertical axis is VV/HH ratio in dB. Each data point is obtained by averaging over all wind directions.</figcaption></figure>

## VI. CONCLUSION

With an efficient SMCG method implemented to solve the PMCHWT SIEs, we have performed NMM3D simulations for ocean radar scattering at L-band frequencies. Surface areas used are up to 15.2 m by 15.2 m (64 wavelengths). For the largest surfaces, there are more than six million surface unknowns in the tangential electric and magnetic fields. This use of a large surface for simulations allows us to examine the effects of gravity and capillary waves. The numerical simulations, unlike the two-scale scattering model, do not require a wavenumber dividing line between large- and small-scale roughness. Simulations are shown for 5, 8, and 10 m/s using isotropic DV spectrum and anisotropic DV spectrum. The backscatter results are in good agreement with Aquarius observation data for VV and HH co-polarization, cross-polarization, and VV/HH ratio. A rigorous approach to the problem based on NMM3D is to incrementally increase the surface size until convergence is reached. Recent computational electromagnetic efforts about scattering by large objects are on using rigorous domain decomposition which should be considered in remote sensing simulations. Without a large domain, it is possible that tilt effects are not included. To obtain tilt effects require a definition of decomposition and tilts can both decrease or increase local incident angles, requiring a heuristic use of tuning parameters. For inclusion of tilt effects, the readers can refer to [7].

For NMM3D to move forward, the increase of domain size via domain decomposition and the increase of sampling points per wavelength are essential steps. Also, the authors have used (k₀/2) as the dividing TSM wavenumber for separation in application of SPM and GO approximations. The results produced by TSM can agree reasonably well with Aquarius. But if another dividing wavenumber (k₉) is chosen, the scattering coefficients will vary. This uncertainty in the TSM approach remains.

## APPENDIX A

In order to integrate SMCG with the PMCHWT equation, it is necessary to write the equation in a form where the field and source points are separated. We will use the subscript *m* to denote the field point and subscript *n* to denote the source point. The superscripts *x*, *y*, and *z* represent the three axes. The following discretized impedance matrix expressions are obtained after performing moment method operations to the kernel of the electric field integral equation:

$$ \int_{S_m} \int_{S_n} \left( \bar{f}_m \cdot \bar{f}_n - \frac{1}{k^2} \nabla \cdot \bar{f}_m \nabla' \cdot \bar{f}_n \right) g(\bar{r}, \bar{r}') d\bar{r} d\bar{r}' \quad (17) $$

where $\bar{f}_m$ and $\bar{f}_n$ are the RWG basis functions for the field and source edges. Also

$$ \bar{f}_m(\bar{r}) g(\bar{r}, \bar{r}') \cdot \bar{f}_n(\bar{r}') = f_m^x(\bar{r}) g(\bar{r}, \bar{r}') f_n^x(\bar{r}') \\ + f_m^y(\bar{r}) g(\bar{r}, \bar{r}') f_n^y(\bar{r}') \\ + f_m^z(\bar{r}) g(\bar{r}, \bar{r}') f_n^z(\bar{r}') \quad (18) $$

$$ \nabla \cdot \bar{f}_m g(\bar{r}, \bar{r}') \nabla' \cdot \bar{f}_n = \left( \pm \frac{l_m}{S_m^{\pm}} \right) g(\bar{r}, \bar{r}') \left( \pm \frac{l_n}{S_n^{\pm}} \right). \quad (19) $$

The kernel of the magnetic field integral equation is

$$ \bar{f}_m(\bar{r}) \times \nabla g(\bar{r}, \bar{r}') \cdot \bar{f}_n(\bar{r}') = \bar{f}_m(\bar{r}) \times \frac{(-1 + ik|\bar{r} - \bar{r}'|) \exp(ik|\bar{r} - \bar{r}'|)}{4\pi|\bar{r} - \bar{r}'|} (\bar{r} - \bar{r}') \cdot \bar{f}_n(\bar{r}'). \quad (20) $$

We use

$$ G(\bar{r}, \bar{r}') = \frac{(-1 + ik|\bar{r} - \bar{r}'|) \exp(ik|\bar{r} - \bar{r}'|)}{4\pi|\bar{r} - \bar{r}'|} \quad (21) $$
Then

$$
\begin{align}
\tilde{f}_m(\bar{r}) \times \nabla g(\bar{r}, \bar{r}') \cdot \tilde{f}_n(\bar{r}') & \\
&= \bar{f}_m(\bar{r}) \times (\bar{r} - \bar{r}') \cdot \bar{f}_n(\bar{r}') G(\bar{r}, \bar{r}') \nonumber \\
&= f_m^y r_m^x G(\bar{r}, \bar{r}') f_n^z - f_m^z r_m^x G(\bar{r}, \bar{r}') f_n^y - f_m^y G(\bar{r}, \bar{r}') f_n^z r_n^x \nonumber \\
&\quad + f_m^z G(\bar{r}, \bar{r}') f_n^y r_n^x + f_m^z r_m^y G(\bar{r}, \bar{r}') f_n^x - f_m^x r_m^y G(\bar{r}, \bar{r}') f_n^z \nonumber \\
&\quad - f_m^z G(\bar{r}, \bar{r}') f_n^x r_n^y + f_m^x G(\bar{r}, \bar{r}') f_n^z r_n^y + f_m^x r_m^z G(\bar{r}, \bar{r}') f_n^y \nonumber \\
&\quad - f_m^y r_m^z G(\bar{r}, \bar{r}') f_n^x - f_m^x G(\bar{r}, \bar{r}') f_n^y r_n^z + f_m^y G(\bar{r}, \bar{r}') f_n^x r_n^z. \tag{22}
\end{align}
$$

Equations (18), (19), and (22) are in a form such that the source point is on the right side of Green's function while the field point resides on the left side of Green's function. Next, Green's function is expanded in a Taylor series with respect to the flat surface. The first four coefficients are listed in the following:

$$
a_0^{(1,2)}(\rho_R) = (1 - ik_{1,2}\rho_R) \frac{\exp(ik_{1,2}\rho_R)}{4\pi\rho_R} \quad (23)
$$

$$
a_1^{(1,2)}(\rho_R) = \frac{\exp(ik_{1,2}\rho_R)}{4\pi\rho_R} \left\{ \frac{k_{1,2}^2}{2\rho_R} + \frac{3ik_{1,2}}{2\rho_R^2} - \frac{3}{2\rho_R^3} \right\} \quad (24)
$$

$$
a_2^{(1,2)}(\rho_R) = \frac{\exp(ik_{1,2}\rho_R)}{4\pi\rho_R} \left\{ \frac{ik_{1,2}^3}{8} - \frac{6k_{1,2}^2}{8\rho_R} - \frac{15ik_{1,2}}{8\rho_R^2} + \frac{15}{8\rho_R^3} \right\} \quad (25)
$$

$$
a_3^{(1,2)}(\rho_R) = \frac{\exp(ik_{1,2}\rho_R)}{4\pi\rho_R} \left\{ -\frac{k_{1,2}^4 \rho_R}{48} - \frac{10ik_{1,2}^3}{48} + \frac{105k_{1,2}^2}{112\rho_R} \right. \\
\left. \qquad + \frac{35ik_{1,2}}{16\rho_R^2} - \frac{35}{16\rho_R^3} \right\} \quad (26)
$$

$$
b_0^{(1,2)}(\rho_R) = \frac{\exp(ik_{1,2}\rho_R)}{4\pi\rho_R} \qquad (27)
$$

$$
b_1^{(1,2)}(\rho_R) = \exp(ik_{1,2}\rho_R) \left\{ \frac{ik_{1,2}}{8\pi} - \frac{1}{8\pi\rho_R} \right\} \quad (28)
$$

$$
b_2^{(1,2)}(\rho_R) = \exp(ik_{1,2}\rho_R) \left\{ -\frac{k_{1,2}^2 \rho_R}{32\pi} - \frac{3ik_{1,2}}{32\pi} + \frac{3}{32\pi\rho_R} \right\} \quad (29)
$$

$$
b_3^{(1,2)}(\rho_R) = \exp(ik_{1,2}\rho_R) \left\{ -i \frac{k_{1,2}^3 \rho_R^2}{192\pi} + \frac{k_{1,2}^2 \rho_R}{32\pi} + \frac{13ik_{1,2}}{192\pi} - \frac{13}{192\pi\rho_R} \right\}. \quad (30)
$$

We have used the centroid as the integration point within each triangular patch [26].

As is shown in Fig. 8(a), the local coordinates of the two integration points A and B are (1/4, 3/4) and (3/4, 1/4) which are not centroids of the two triangles. In Fig. 8(b), the local coordinates of C and D are (1/3, 2/3) and (2/3, 1/3), respectively. The points C and D fall onto the center of each triangle. After making the modification from Fig. 8(a) to (b), the accuracy of SMCG is greatly improved.

We define the error of SMCG by

$$
\text{error} = \frac{\|z_{mn} - z_{mn}^{\text{SMCG}}\|}{\|z_{mn}\|} \qquad (31)
$$

<figure><img src="image_9.png" /><figcaption>Fig. 8. SMCG with centroid as integration points. (a) Local coordinates of A, B: (1/4, 3/4) and (3/4, 1/4). (b) Local coordinates of C, D: (1/3, 2/3) and (2/3, 1/3).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9. SMCG in conjunction with RWG basis functions.</figcaption></figure>

where $z_{mn}$ are the impedance matrix elements calculated directly using seventh-order quadrature and $z_{mn}^{\text{SMCG}}$ are the matrix elements calculated indirectly by the SMCG method. For a $4\lambda$ by $4\lambda$ surface with 16 points per wavelength, the error is about 0.2%. In addition to SMCG, we have also used Green's function interpolation in the numerical simulations [9].

In Fig. 9, we illustrate how SMCG based on RWG functions works. First, the currents inside one source triangle patch aggregate to the centroid. Next, the accumulated currents transfer to the field triangle through Green's function. Then, the transferred field distributes to each inner edge inside the source triangle.

APPENDIX B

The ocean surface height spectrum for a fully developed wind sea, as proposed by Johnson et al. [21], has the following form:

$$
W_{\text{DV}}(k \cdot \phi) = \frac{a_0}{2\pi k^4} \Phi_{\text{DV}}(k, \phi) \begin{cases} \exp\left(-\frac{\beta g^2}{k^2 U_{19,5}^4}\right), & k < k_j \\ \left(\frac{b k u_*^2}{g_*}\right)^{a \log_{10}(k/2)}, & k \ge k_j \end{cases} \tag{32}
$$

where $k_j = 2m^{-1}$. The subscript DV denotes the DV spectrum. The constant $a_0$ is 0.004 in [21], but was adjusted upward to 0.008 in [20]. In order to best fit the L-band measurement data in this paper, a value of $a_0 = 0.008 \times 1.3$ is applied. Other parameters are $\beta = 0.74$, $b = 1.25$, $a = 0.225$, $g_* = g + \gamma k^2$, $g = 9.81$ m/sec is acceleration due to gravity, and $\gamma = 7.25 \times 10^{-5}$ m$^3$/sec$^2$ is the ratio of surface tension to
water density. The wind friction velocity at the ocean surface $u_*$ is related to the wind speed at a height ($h$) of 19.5 m above sea level $U_{19.5}$ using the following form:

$$U_h = \frac{u_*}{0.4} \log \frac{h}{0.0000684/u_* + 0.00428u_*^2 - 0.000443}. \quad (33)$$

The wave directional spreading factor $\Phi(k, \phi)$ is

$$\Phi(k, \phi) = 1 + \Delta(k)(1 - \exp(-sk^2)) \cos(2(\phi - \phi_{\text{wind}})). \quad (34)$$

In (34), the upwind direction is defined as $\phi_{\text{wind}} = 0^{\circ}$. In this paper, $\Delta(k)$ is defined as

$$\Delta(k) = c \frac{(1 - R)}{(1 + R)} \frac{2}{1 - D} \quad (35)$$

with

$$R = \frac{0.003 + 1.92 \times 10^{-3} U_{12.5}}{3.16 \times 10^{-3} U_{12.5}} \quad (36)$$

$$D = \frac{\int_{0}^{\infty} dkk^{2}S(k) \exp(-sk^{2})}{\int_{0}^{\infty} dkk^{2}S(k)}. \quad (37)$$

## ACKNOWLEDGMENT

The high-performance computing facilities of XSEDE Stampede at NSF and the FLUX at the University of Michigan, Ann Arbor, MI, USA, were used. The part of work carried out by the Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA, USA, is under a contract with the National Aeronautics and Space Administration, Washington, DC, USA.

## REFERENCES

[1] A. K. Fung and K. S. Chen, Microwave Scattering and Emission Models for Users. Norwood, MA, USA: Artech House, 2010.

[2] L. Tsang, J. A. Kong, and R. T. Shin, Theory of Microwave Remote Sensing. Hoboken, NJ, USA: Wiley, 1985.

[3] K. S. Chen, T.-D. Wu, L. Tsang, Q. Li, J. Shi, and A. K. Fung, "Emission of rough surfaces calculated by the integral equation method with comparison to three-dimensional moment method simulations," IEEE Trans. Geosci. Remote Sens., vol. 41, no. 1, pp. 90-101, Jan. 2003.

[4] F. L. Teixeira, R. J. Burkholder, J. T. Johnson, and C. Yardim, "An efficient implementation of the nonlocal small-slope approximation for bistatic scattering from two-dimensional rough surfaces," in Proc. USNC-URSI Radio Sci. Meeting, Jun./Jul. 2016, pp. 95-96.

[5] C.-A. Guérin and J. T. Johnson, "A simplified formulation for rough surface cross-polarized backscattering under the second-order small-slope approximation," IEEE Trans. Geosci. Remote Sens., vol. 53, no. 11, pp. 6308-6314, Nov. 2015.

[6] A. Voronovich, "Small-slope approximation for electromagnetic wave scattering at a rough interface of two dielectric half-spaces," Waves Random Media, vol. 4, no. 3, pp. 337-367, 1994.

[7] W. J. Plant, "A two-scale model of short wind-generated waves and scatterometry," J. Geophys. Res., Oceans, vol. 91, no. C9, pp. 10735-10749, Sep. 1986.

[8] S. H. Yueh, "Modeling of wind direction signals in polarimetric sea surface brightness temperatures," IEEE Trans. Geosci. Remote Sens., vol. 35, no. 6, pp. 1400-1418, Nov. 1997.

[9] S. Huang, L. Tsang, E. G. Njoku, and K. S. Chan, "Backscattering coefficients, coherent reflectivities, and emissivities of randomly rough soil surfaces at L-band for SMAP applications based on numerical solutions of Maxwell equations in three-dimensional simulations," IEEE Trans. Geosci. Remote Sens., vol. 48, no. 6, pp. 2557-2568, Jun. 2010.

[10] W. Yang, Z. Zhao, and Z. Nie, "Fast Fourier transform multilevel fast multipole algorithm in rough ocean surface scattering," Electromagnetics, vol. 29, no. 7, pp. 541-552, Sep. 2009.

[11] S. L. Durden and J. F. Vesecky, "A numerical study of the separation wavenumber in the two-scale scattering approximation (ocean surface radar backscatter)," IEEE Trans. Geosci. Remote Sens., vol. 28, no. 2, pp. 271-272, Mar. 1990.

[12] S. Huang and L. Tsang, "Electromagnetic scattering of randomly rough soil surfaces based on numerical solutions of Maxwell equations in three-dimensional simulations using a hybrid UV/PBTG/SMCG method," IEEE Trans. Geosci. Remote Sens., vol. 50, no. 10, pp. 4025-4035, Oct. 2012.

[13] T.-H. Liao et al., "Copolarized and cross-polarized backscattering from random rough soil surfaces from L-band to Ku-band using numerical solutions of Maxwell's equations with near-field precondition," IEEE Trans. Geosci. Remote Sens., vol. 54, no. 2, pp. 651-662, Feb. 2016.

[14] J. Yang, Y. Du, and J. Shi, "Polarimetric simulations of bistatic scattering from perfectly conducting ocean surfaces with 3 m/s wind speed at L-band," IEEE. Sel. Topics Appl. Earth Observ. Remote Sens., vol. 9, no. 3, pp. 1176-1186, Mar. 2016.

[15] L. Tsang, J. A. Kong, and R. T. Shin, Scattering of Electromagnetic Waves: Numerical Simulations, vol. 2. Hoboken, NJ, USA: Wiley, 2000.

[16] J.-W. Hao, W. Song, and X.-Q. Sheng, "Efficient algorithm for calculating backscattering from two-dimensional rough sea surface under low grazing angle," in Proc. IEEE Int. Geosci. Remote Sens. Symp. (IGARSS), Jul. 2016, pp. 3695-3698.

[17] J. T. Johnson and H. T. Chou, "Numerical studies of low grazing angle backscatter from 1D and 2D impedance surfaces," in Proc. Geosci. Remote Sens. Symp., vol. 4. 1998, pp. 2295-2297.

[18] D. Miret, G. Soriano, F. Nouguier, P. Forget, M. Saillard, and C.-A. Guerin, "Sea surface microwave scattering at extreme grazing angle: Numerical investigation of the Doppler shift," IEEE Trans. Geosci. Remote Sens., vol. 52, no. 11, pp. 7120-7129, Nov. 2014.

[19] G. Soriano, C.-A. Guérin, and M. Saillard, "Microwave ocean scattering at low-grazing angles with the GMoM," in Proc. Eur. IEEE Radar Conf. (EuRAD), Sep./Oct. 2010, pp. 5-8.

[20] P. Spiga, G. Soriano, and M. Saillard, "Scattering of electromagnetic waves from rough surfaces: A boundary integral method for low-grazing angles," IEEE Trans. Antennas Propag., vol. 56, no. 7, pp. 2043-2050, Jul. 2008.

[21] J. T. Johnson, R. T. Shin, J. A. Kong, L. Tsang, and K. Pak, "A numerical study of the composite surface model for ocean backscattering," IEEE Trans. Geosci. Remote Sens., vol. 36, no. 1, pp. 72-83, Jan. 1998.

[22] Q. Li, C. H. Chan, L. Tsang, and G. Zhang, "Wave scattering from a lossy dielectric rough surface using PBTG-BMIA/CAG and applications to passive remote sensing," in Proc. IEEE Int. Geosci. Remote Sens. Symp. (IGARSS), vol. 3. Jul. 1998, pp. 1469-1471.

[23] D. Miret, G. Soriano, and M. Saillard, "Rigorous simulations of microwave scattering from finite conductivity two-dimensional sea surfaces at low grazing angles," IEEE Trans. Geosci. Remote Sens., vol. 52, no. 6, pp. 3150-3158, Jun. 2014.

[24] S.-Q. Li, C.-H. Chan, L. Tsang, Q. Li, and L. Zhou, "Parallel implementation of the sparse-matrix/canonical grid method for the analysis of two-dimensional random rough surfaces (three-dimensional scattering problem) on a Beowulf system," IEEE Trans. Geosci. Remote Sens., vol. 38, no. 4, pp. 1600-1608, Jul. 2000.

[25] S. H. Yueh, W. J. Wilson, and S. Dinardo, "Polarimetric radar remote sensing of ocean surface wind," IEEE Trans. Geosci. Remote Sens., vol. 40, no. 4, pp. 793-800, Apr. 2002.

[26] S. Durden and J. Vesecky, "A physical radar cross-section model for a wind-driven sea with swell," IEEE J. Ocean. Eng., vol. OE-10, no. 4, pp. 445-451, Oct. 1985.

[27] N. Pinel, B. Chapron, C. Bourlier, N. de Beaucoudrey, R. Garello, and A. Ghaleb, "Statistical analysis of real aperture radar field backscattered from sea surfaces under moderate winds by Monte Carlo simulations," IEEE Trans. Geosci. Remote Sens., vol. 52, no. 10, pp. 6459-6470, Oct. 2014.

[28] T. Elfouhaily, B. Chapron, K. Katsaros, and D. Vandemark, "A unified directional spectrum for long and short wind-driven waves," J. Geophys. Res., vol. 102, no. C7, pp. 15781-15796, Jul. 1997.

[29] F. Nouguier, C.-A. Guérin, and B. Chapron, "'Choppy wave' model for nonlinear gravity waves," J. Geophys. Res., vol. 114, no. C9, p. C09012, Sep. 2009.

[30] A. K. Fung, Backscattering from Multiscale Rough Surfaces With Application to Wind Scatterometry. Norwood, MA, USA: Artech House, 2015.

[31] H. Braunisch et al., "Tapered wave with dominant polarization state for all angles of incidence," IEEE Trans. Antennas Propag., vol. 48, no. 7, pp. 1086-1096, Jul. 2000.
[32] A. G. Voronovich and V. U. Zavorotny, "Sensitivity of modeled polari-
metric radar ocean scattering to wind direction," in *Proc. IEEE Geosci.
Remote Sens. Symp. (IGARSS)*, Jul. 2014, pp. 5021–5024.

[33] P. A. Hwang and D. W. Wang, "An empirical investigation of source term balance of small scale surface waves," *Geophys. Res. Lett.*, vol. 31, no. 15, p. L15301, Aug. 2004.

[34] Y. Liu, M.-Y. Su, X.-H. Yan, and W. T. Liu, "The mean-square slope of ocean surface waves and its effects on radar backscatter," *J. Atmos. Ocean. Technol.*, vol. 17, no. 8, pp. 1092–1105, Aug. 2000.

[35] C. Cox and W. Munk, "Statistics of the sea surface derived from sun glitter," *J. Marine Res.*, vol. 13, no. 2, pp. 198-227, 1954.

[36] J. Wright, "A new model for sea clutter," *IEEE Trans. Antennas Propag.*, vol. AP-16, no. 2, pp. 217-223, Mar. 1968.

[37] O. Isoguchi and M. Shimada, "An L-band ocean geophysical model function derived from PALSAR," *IEEE Trans. Geosci. Remote Sens.*, vol. 47, no. 7, pp. 1925-1936, Jul. 2009.

[38] M. V. Yurovskaya, V. A. Dulov, B. Chapron, and V. N. Kudryavtsev, "Directional short wind wave spectra derived from the sea surface photography," *J. Geophys. Res.*, vol. 118, no. 9, pp. 4380-4394, Sep. 2013.

[39] J. T. Johnson et al., "Backscattering enhancement of electromagnetic waves from two-dimensional perfectly conducting random rough surfaces: A comparison of Monte Carlo simulations with experimental data," *IEEE Trans. Antennas Propag.*, vol. 44, no. 5, pp. 748-756, May 1996.

[40] S. H. Yueh et al., "L-band passive and active microwave geophysical model functions of ocean surface winds and applications to aquar-ius retrieval," *IEEE Trans. Geosci. Remote Sens.*, vol. 51, no. 9, pp. 4619-4632, Sep. 2013.

[41] Z. Zhao and J. C. West, "Low-grazing-angle microwave scattering from a three-dimensional spilling breaker crest: A numerical investigation," *IEEE Trans. Geosci. Remote Sens.*, vol. 43, no. 2, pp. 286-294, Feb. 2005.

[42] D. R. Thompson, T. M. Elfouhaily, and J. L. Garrison, "An improved geometrical optics model for bistatic GPS scattering from the ocean surface," *IEEE Trans. Geosci. Remote Sens.*, vol. 43, no. 12, pp. 2810-2821, Dec. 2005.

<figure><img src="image_11.png" /></figure>

**Douglas Vandemark** received the B.S. degree in physics from Hope College, Holland, MI, USA, in 1986, the M.S. degree in electrical engineering from the University of Massachusetts, Amherst, MA, USA, in 1988, and the Ph.D. degree in oceanography from the University of New Hampshire (UNH), Durham, NH, USA, in 2005.

From 1990 to 2005, he was with the NASA's Goddard Space Flight Center, Greenbelt, MD, USA. Since 2005, he has been with the Institute for the Study of Earth, Oceans, and Space, UNH, where he is currently a Research Professor. He currently serves as a Principal Investigator with the Ocean Surface Topography, Ocean Vector Wind, and Ocean Surface Salinity Science Teams, NASA. He is also a Federal Advisory Committee Member of the U.S. Integrated Ocean Observing System. His research interests include air-sea interaction with a focus on ocean remote sensing of wind waves and the exchange of momentum, heat, and mass including greenhouse gases.

Dr. Vandemark is a member of the American Geophysical Union, the American Meteorological Society, and the IEEE Geoscience and Remote Sensing and Oceanic Engineering Societies.

<figure><img src="image_12.png" /></figure>

**Tai Qiao** (S'15) received the B.S. degree in electrical engineering from the University of Electronic Science and Technology of China, Chengdu, China, in 2011, and the M.S. degree in electrical engineering from Southeast University, Nanjing, China, in 2014. He is currently pursuing the Ph.D. degree with the University of Michigan, Ann Arbor, MI, USA.

His research interests include electromagnetic the-
ory, wave scattering, and remote sensing of ocean.

<figure><img src="image_13.png" /></figure>

**Simon H. Yueh** (F'09) received the Ph.D. degree from the Massachusetts Institute of Technology, Cambridge, MA, USA, in 1991.

He joined the Jet Propulsion Laboratory (JPL), Pasadena, CA, USA, in 1991. He served as the Project Scientist of the National Aeronautics and Space Administration (NASA) Aquarius Mission from 2012 to 2013, the Deputy Project Scientist of the NASA Soil Moisture Active Passive Mission from 2013 to 2013, and has been the SMAP Project Scientist since 2013.

Dr. Yueh received the 2014 IEEE GRSS TRANSACTIONS Prize Paper Award, the 2010 IEEE GRSS TRANSACTIONS Prize Paper Award, the 2002 IEEE GRSS TRANSACTIONS Prize Paper Award, the 2000 Best Paper Award in the IEEE International Geoscience and Remote Symposium, the 1995 IEEE GRSS TRANSACTIONS Prize Paper Award for a paper on polarimetric radiometry, the JPL Lew Allen Award in 1998, and the Ed Stone Award in 2003. He was an Associate Editor of the *Radio Science* from 2002 to 2006. He is the Editor of the IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING.

<figure><img src="image_14.png" /></figure>

**Leung Tsang** (F'90) was born in Hong Kong. He received the S.B., S.M., and Ph.D. degrees from the Massachusetts Institute of Technology, Cambridge, MA, USA, in 1971, 1973, and 1976, respectively.

From 1983 to 2014, he was a Professor at the University of Washington, Seattle, WA, USA, where he was the Chair of the Department of Electrical Engineering from 2006 to 2011. From 2001 to 2004, he was on leave and was a Professor at the City University of Hong Kong, Hong Kong. He is currently a Professor with the Department of Electrical Engineering and Computer Science, University of Michigan, Ann Arbor, MI, USA. He has co-authored the *Theory of Microwave Remote Sensing* and the *Scattering of Electromagnetic Waves* Volumes 1-3. His research interests include waves in random media, rough surfaces, remote sensing, computational electromagnetics, signal integrity, electromagnetic compatibility, and plasmonics.

Dr. Tsang was the President of the IEEE Geoscience and Remote Sensing Society from 2006 to 2007. He was a recipient of the William Pecora Award cosponsored by USGS and NASA in 2012, and the IEEE Electromagnetics Award in 2013. He is on the Editorial Board of the IEEE ACCESS and is the Chair of *PIERS, Progress in Electromagnetics Research*.

<figure><img src="image_15.png" /></figure>

**Tien-Hao Liao** was born in Taichung, Taiwan. He received the B.S. degree in mathematics from National Taiwan University, Taipei, Taiwan, in 2003, the M.S. degree in communication engineering from National Chiao Tung University, Hsinchu, Taiwan, in 2005, and the Ph.D. degree in electrical engineering from the University of Washington, Seattle, WA, USA, in 2015.

He was a Research Fellow with the Radiation Laboratory, University of Michigan, Ann Arbor, MI, USA. In 2016, he received the Fellowship from NASA Post-Doctoral Program and was affiliated with the Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA, USA, where he is currently a Post-Doctoral Researcher. His research interests include electromagnetic scattering model for vegetated surface, random rough surface for microwave remote sensing, radar polarimetry and radar interferometry using airborne and satellite remote sensing data for earth science, and microwave waveguide modeling for broadband Green's function.
<figure><img src="image_16.png" /></figure>

Frédéric Nouguier received the “Agrégation” and M.S. degrees in applied physics from the École Nor- male Supeèrieure de Cachan, Cachan, France, in 2005 and 2006, respectively, the M.S. degree in physical methods for remote sensing from the University of Paris-Diderot, Paris, France, in 2006, and the Ph.D. degree in physics from the University of Marseille, Marseille, France, in 2009.

He is currently a Research Scientist with the Laboratoire d'Oceianographie Physique et Spatiale, Institut Francais de Recherche pour là'Exploitation de la Mer, Plouzanei, France. His research interests include applied math- ematics, physical oceanography, and electromagnetic wave theory and its application to ocean remote sensing.

<figure><img src="image_17.png" /></figure>

Bertrand Chapron received the Ph.D. degree in fluid mechanics from Aix-Marseille University, Marseille, France, in 1988.

He is currently a Senior Research Scientist with the Laboratoire d'Océanogr aphie Physique Spatiale, Institut Francais de Recherche pour là'Exploitation de la Mer, Plouzaneï, France. He has been the Co-Investigator or a Principal Investigator in several ESA (ENVISAT RA2, ASAR, and SMOS), NASA, and CNES (TOPEX and JASON) projects. He con- tributed to the development of several algorithms for geophysical parameter retrieval from altimeters, radiometers, or SAR. He is co-responsible for the ENVISAT ASAR-Wave Mode algorithms and scientific preparation for the ENVISA T and S1 wind, wave, and current. Recently, he collaborated in studies that demonstrated the high potential of radiometers for wind retrieval in extreme conditions and the possibility of doing directional ocean waves spectrum from Sentinel-2. He has authored over 80 publications in refereed journals in applied mathematics, physical oceanography (upper ocean dynamics), and electromagnetic wave theory and its application to ocean surface remote sensing.