Abstract
A new electromagnetic (EM) scattering model of the sea surface with single breaking waves is proposed based on the high-frequency method in this paper. At first, realistic breaking wave sequences are obtained by solving the fluid equations which are simplified. Then, the rough sea surface is established using the linear filtering method. A new wave model is obtained by combining breaking waves with rough sea surface using a 3D coordinate transformation. Finally, the EM scattering features of the sea surface with breaking waves are studied by using shooting and bouncing rays and the physical theory of diffraction (SBR-PTD). It is found that the structure that is similar to a dihedral corner reflector between the breaking wave and rough sea surface exhibits multiple scattering, which leads to the sea-spike phenomenon that the scattering result of horizontal (HH) polarization is larger than that of vertical (VV) polarization, especially at low-grazing-angle (LGA) incidents with upwind. The sea-spike phenomenon is also closely related to the location of strong scattering.
1. Introduction
Target recognition in a complex ocean background is a focus problem that received large attention. Sea clutter influences the detection of targets in extreme sea conditions, which in turn increases the probability of false alarms. Therefore, it is critical to construct a realistic breaking wave model and study its electromagnetic (EM) scattering characteristics.
LONGTANK profiles [1] laid the foundation for research on breaking waves. Waves at different times are spliced together in two-dimensional LONGTANK profiles to generate three-dimensional breaking waves [2, 3]. In particular, three-dimensional breaking waves have been generated by the finite element method combined with the ocean dynamics [4]. One-dimensional time-varying breaking waves have been discussed [5]. A limited size dihedral corner reflector split structure has been used to approximate the overflow breakage [6]. Breaking waves have been determined by solving the governing equation of fluid mechanics [7, 8]. However, most of the previous methods are not realistic enough and often inefficient. Moreover, they may be only employed to study small areas of sea surface containing breaking waves.
Concerning EM scattering, scholars have studied EM scattering characteristics of breaking waves based on the LONGTANK model. Holliday et al. used the method of moments (MoM) with backward and forward iteration to calculate the scattering of breaking waves with different geometrical shapes [9, 10]. West and Zhiqin Zhao studied EM scattering of breaking water waves with rough faces [11]. Zhiqin Zhao and West [12] used the multilayer fast multipole method to investigate scattering characteristics of three-dimensional breaking waves. Superevents where scattering due to horizontal (HH) is larger than that due to vertical (VV) have been found in the above numerical simulations.
In this work, realistic breaking wave sequences are obtained by solving the fluid equations which are simplified. The model includes the detailed characteristics of the wave at different times and is therefore inherently efficient. Moreover, a new sea surface with the single breaking wave model is obtained by combining breaking waves with rough sea surface using a 3D coordinate transformation. It is necessary to take into account the multiple scattering produced by breaking wave itself, and the coupling effect between breaking wave and sea surface. There are two methods for doing this. The first is the numerical method where it is suitable for calculating small objects. The second one is a high-frequrency approximation method which suites large objects. The given sea surface with a breaking wave is a large target; therefore, the shooting and bouncing rays and physical theory of diffraction (SBR-PTD) method is used here. Indeed, the SBR method effectively combines the advantages of geometrical (GO) and physical optics (PO). In particular, it is necessary to take into account the multiple scattering produced by breaking wave itself, and the coupling effect between breaking wave and sea surface.
2. Mathematical Model and Theoretical Formulations
2.1. Three-Dimensional Simulation of Breaking Waves
Many popular methods for modeling water surfaces work well in producing static images, but these methods only attempt to depict a surface with roughness, and they are not suitable for the realistic dynamics of the surface over time [13, 14]. Moreover, traditional modeling methods based on the shallow water approximation are inefficient. Therefore, an improved breaking wave model is established here building on reference [15].
In practice, breakers are normally associated with shorelines but can occur anywhere in the ocean. The major disturbing force in the open ocean is wind. In this paper, the plunging breakers are studied. They occur when the wave encounters an abrupt transition from deep to shallow water. The base of the wave decelerates rapidly, while the top of the wave continues moving at a higher speed. With this large-speed differential, the top of the wave pitches out in front, forming a curl or tube. Moreover, we just make an assumption and establish a geometric model of the breaking waves, which is combined with the rough sea surface. The model is applied to the study of EM scattering.
The process is divided into three steps:(1)The realistic breaking wave sequences are obtained by solving the fluid equations which are simplified.(2)The sea surface is modeled by the Monte Carlo method based on the Elfouhaily [16] sea spectrum.(3)A new sea surface with the single breaking wave model is obtained by combining breaking waves with the rough sea surface using a 3D coordinate transformation.
Step 1 and Step 3 will be described in detail in the following.
We begin with a vastly simplified set of equations that has been widely used for shallow water. The simplification arises from three approximations:(a)We assume that the water surface is a height field. Figure 1 shows the discrete representation of the height field in two dimensions. This, of course, has some obvious limitations. The water cannot splash and waves cannot break. However, so long as the forces on the water are sufficiently gently, the height-field assumption will not introduce error.(b)The vertical component of the velocity of the water particles can be ignored. Once again, the limitations of this assumption are fairly clear. If a disturbance creates very steep waves on the water surface, the model will cease to be accurate.(c)The horizontal component of the velocity of the water in a vertical column is assumed constant.

If there is turbulent flow or unusually high friction on the bottom, this assumption will break down.
Those assumptions lead to some obvious limitations, but the experience of hydrodynamicists suggests that they provide a rather useful model to describe breaking sea wave [15].
For simplicity, we begin with a height-field curve in two dimensions. Later, the same techniques will be extended to a height-field surface in three dimensions. Let be height of the water surface and let be the height of the ground. Let is the water depth and is the horizontal velocity of a vertical column of water. Waves where are called deep water and where are called shallow water waves [17]. In this, is the water depth. is the wavelength. The shallow water equations following from the above assumptions [18, 19] may be written as follows:where is the gravitational acceleration. Equation (1) expresses Newton’s law , while equation (2) accounts for volume conservation. Notice that even with the above three simplifying assumptions, the resulting differential equations are still nonlinear. A further simplification which is often used is to ignore the second term in equation (1) and linearize the expression around a constant value of . The resulting equations are then
If we differentiate equation (3) with respect to , then differentiate equation (4) with respect to , and finally substitute for the cross-derivatives, we end up withwhich is the one-dimensional wave equation with wave velocity .
In order to solve equation (5), we need to construct a discrete representation of the continuous partial-differential equation. Here, the finite-difference technique works particularly well because of the simple height-field representation. After experimenting with a number of finite-difference approximations to equations (3) and (4), the most stable version we have found iswhere is the separation of the samples along the direction. Putting the above two equations together, we getwhich provides the sought discrete approximation to equation (5).
For simplicity, we use a first-order implicit method to solve ordinary differential equation (8). Let to denote at the iteration and let dots denote differentiation with time. Then, the first-order implicit equations may be written as
Notice that the right-hand sides of these equations are evaluated at time which corresponds to the end of the iteration rather than time which corresponds to the beginning of the iteration. This is what makes the iteration implicit and stable. Rearranging the above expression, we arrive at
By substituting equation (8) into equation (10), we obtain
After these manipulations, equation (11) is still a nonlinear equation. In order to find its solution, it needs to be linearized, i.e., is regarded as a constant, and the wave velocity is fixed as a function of . In this linear regime, the next value of can be calculated from previous values with the symmetric tridiagonal linear systemwhere the matrix is given byand the elements of are as follows:
Now, let us change the equation towhere is phenomenologically introduced to describe damping phenomena. If , then equation (15) reduces to equation (12), whereas if is between zero and one, it will make the waves damp out over time. The effect is that observed in a viscous fluid. There is one further subtlety of importance in the two-dimensional case. Even though equation (15) was derived from equation (6) which specifies conservation of volume, there is no guarantee that the results of the iteration will precisely conserve volume. The primary cause of departures from volume-conserving behavior is that the iteration may leave for some index . To compensate for this negative volume, the iteration will create excess positive volume elsewhere. While the effect is small, it can accumulate over time and create substantial drift. If the entire surface acquires a small net upwards velocity, it will very quickly create noticeable amounts of water. To combat this effect, the following simple projection appears to be adequate. After each iteration, we find the connected pieces of the fluid, and this can be done by scanning the and vectors in order and testing whether . For each connected piece of the fluid, we calculate the old volume and the new volume. If the new volume is different, we distribute the difference uniformly over the samples in the connected region.
The above treatment is valid to a height field in two dimensions. The three-dimensional case can be approximated by a series of two-dimensional equations. The basic wave equation for water in three dimensions is the same as the two-dimensional case except that the second derivative of with respect to is replaced by the Laplacian.
In order to solve the equations in three dimensions, we rely on the alternating-direction method [15]. The basic idea of the method is to take equation (16) and split the right-hand side into the sum of two terms, one of which is independent of , and the other is independent of . We then divide the iteration into two subiterations. In the first, we replace the right-hand side of equation (17) with the first term, and in the second one, we replace the right-hand side of equation (17) with the second term. More specifically, in the first subiteration, we solve the equationand in the second subiteration, we solve the equation
For the first subiteration, we compute the update as before on each row of the height field. For the second subiteration, we do the same for each column in the height field.
Then, the breaking waves generated are combined with the rough sea surface by 3D coordinate transformation. This method will be described in detail in the following.
Rotation matrices , , and about x, y, and z axes are defined aswhere , , and are angles with , , and axes, respectively.
Scaling matrices are defined aswhere, , , and are scaling factors.
World space matrix is defined aswhere , , and are coordinates of the breaking waves. Let us also define as the inverse of and the vector , where , , and are position coordinates of the rough sea surface. The link between position coordinates and local-coordinate system is a rotation
Let us now definewhere , , L is length of the breaking waves, and width is width of the breaking waves.
According to [20, 21], the relationship between wave height and wind speed iswhere and .
Wave steepness S [22] is defined as the ratio of , i.e.,where S is the wave steepness.
We can get the wavelength according to (25) and (26), i.e.,
The smoothing factor is defined aswhere amplitude = 1. If , then , whereas if , then .
Convert positional coordinates from local to global coordinates, i.e.,where .
To transform local coordinates into global coordinates, we need with 4 rows and 4 columns, and hence let us define
Therefore, the coordinates of the transformed point position arewhere are the coordinates of the sea surface with single breaking waves.
The sea surface with single breaking waves is shown in Figure 2. The evolution process describes the breaking waves from the initial moment to gradual curl. It is worth mentioning that the wave is a model that changes all the time, and Figure 2 is only a realistic description of the state of the wave probably. The time interval is taken as . The size of rough sea surface is with wind speed . , and the width of the breaking waves is 3 m. The height and length of breaking waves can be calculated by equations (25) and (27), respectively. Therefore, the location of the breaking wave can be determined. The aim of this paper is to study EM scattering characteristics in the presence of the sea surface with single breaking waves, so this paper only studies the breaking waves before involution.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)
2.2. SBR-PTD Method
The SBR method is a hybrid method of GO combined with PO, and the PTD method can calculate edge diffraction. Calculation accuracy can be improved using PTD together with SBR. The SBR-PTD method uses GO to determine the triangles illuminated by the incident wave and each order of reflected wave. Moreover, recording the IDs of each triangle’s three neighbors is also useful in the calculation of PTD fields.
The tracing process may be understood by looking at the case depicted in Figure 3, where a ray from triangle patch 0 is shot toward a target, which consists of 16 triangle patches. The following steps should be carried out to determine which triangles are illuminated by the ray from triangle patch 0. At first, forward ray tracing is performed to detect one of these illuminated triangles (i.e., triangle 7). Then, after finding illuminated triangle 7, it is necessary to determine whether the adjacent triangles (i.e., triangles 6, 8, and 13) are illuminated as well (e.g., triangle 13). Backward ray tracing is used to determine whether triangle 13 intersects triangle patch 0 and is not occluded. If so, triangle 13 is illuminated; otherwise, it is discarded. Then, backward ray tracing is used to determine whether triangle 6 and triangle 8 are illuminated or not, and the procedure continues until some triangles outside the illuminated area occur. Moreover, a binary tree is built to accelerate the process of ray tracing.

For each triangle patch, the scattered far field is obtained by vector-summing each order of the PO fields from the triangle surface and the PTD fields from the three edges of this triangle [23].
The electric field, approximated by PO integral, may be the written as follows:where , is the unit vector aligned along the scattering direction, is the distance from the specular point to the observation point, and is the current density on the triangle surface. The diffraction field from each edge of a triangle patch is approximated by the PTD method:where and are as follows:where and denote the incident electric field strength of patches’ edges; stands for the angle between the edge and incident or scattered ray; and , , and represent the PTD diffraction coefficients [24], which are closely related to the angle between two adjacent triangles whose IDs have already been stored in the computer memory.
3. Numerical Simulations and Discussion
3.1. Verification of the Effectiveness of the SBR-PTD Method
Back-scattering RCS of the sea surface with single breaking waves both in stage and stage is shown in Figure 4. The frequency of the incident wave is . Therefore, the dielectric parameter is from the Debye model [25]. The incident angles are , . The size of rough sea surface is with wind speed . It may be readily seen that the result of SBR-PTD agrees well with the result of MoM, which verified the effectiveness of the SBR-PTD method. Note that represents downwind and represents upwind, and this is true for all of the following examples.

(a)

(b)

(c)

(d)
3.2. A Comparison of the Bistatic Scattering Coefficient with and without Breaking Waves
A locally breaking wave resembling a dihedral corner reflector is shown in Figure 5. A comparison of the bistatic scattering coefficient with and without breaking waves is shown in Figure 6. In this, the stage of the sea surface with single breaking waves is selected. The frequency of the incident wave is . Therefore, the dielectric parameter is from the Debye model. The incident angles are , . The size of rough sea surface is with wind speed . In Figure 6, for the sea surface with single breaking waves, there is an obvious enhancement in the backward direction. This is caused by the fact that the structure between the breaking waves and the rough sea surface (i.e. Figure 5) is similar to a dihedral corner reflector, and thus it exhibits multiple scattering, resulting in the change of scattering field near the backward nonspecular direction which is small, while the sea surface scattering has a dominant contribution at the specular direction.


(a)

(b)
3.3. Back-Scattering RCS of Sea Surface with Single Breaking Waves for Different Stages
To numerically demonstrate the sea-spike phenomenon of breaking waves, , , , and , which denote four different time evolution histories, are chosen to analyze the EM scattering mechanism. Back-scattering RCS of the sea surface with single breaking waves in the , , , and stages is shown in Figure 7. The frequency of the incident wave is . Therefore, the dielectric parameter is from the Debye model. The incident angles are , . The size of rough sea surface is with wind speed .

(a)

(b)

(c)

(d)
In Figure 7, for stage of the sea surface with single breaking waves, it is at the very beginning of the plume forming. The multipath effects are not serious. HH is almost equal to VV. For the , , and stages of the sea surface with single breaking waves, the plume forms. HH is larger than VV at the grazing angles of less than , and the sea-spike phenomenon has occurred. This is caused by the fact that the structure is similar to a dihedral corner reflector between the breaking waves and the rough sea surface (i.e., Figure 5), and thus it exhibits multiple scattering which, in turn, enhances the Brewster effect, which will cause the reflection coefficient to decrease at VV polarization, and the VV multipath is greatly attenuated.
3.4. Back-Scattering Field Intensity Distribution of the Sea Surface with Breaking Waves
The back-scattering field intensity distribution for different incidence angles of each triangle patch in the stage of the sea surface with single breaking waves is shown in Figure 8. The frequency of the incident wave is and the corresponding dielectric parameter is . The size of rough sea surface is and the wind speed . The incident angles are (downwind/upwind), (downwind), and (upwind), respectively.

(a)

(b)

(c)

(d)
For (downwind/upwind), most of triangle patches on the model are illuminated. Only the inside of the breaking waves is obscured. At this time, the scattering of the rough sea surface is dominant. For (downwind/upwind), a few triangle patches on the model are illuminated by EM waves. Most of the surface is partially obscured by the breaking waves of the model at low grazing angle (LGA). At this time, the scattering of the breaking waves is dominant. In Figure 8, for moderate incident angles, the main scattering comes from the rough sea surface of the model. For LGA, the main scattering comes from the breaking waves of the model. Moreover, a few triangles patches with warm colors have a stronger scattering effect compared to other triangle patches with cool colors, which reflects the strong scattering location of the model, i.e., the structure is similar to a dihedral corner reflector formed between the breaking waves and the rough sea surface (i.e., Figure 5).
3.5. Influence of Incidence Angle on EM Scattering of Sea Surface with Single Breaking Waves
Back-scattering RCS for different incidence angles in the stage of the sea surface with single breaking waves is shown in Figure 9. The frequency of incident wave is , and the dielectric parameter from the Debye model. The size of the rough sea surface is with wind speed . Incidence angles are (downwind/upwind), (downwind/upwind), and (downwind/upwind), respectively.

(a)

(b)

(c)

(d)

(e)

(f)
In Figure 9, for (downwind/upwind), HH is larger than VV at some angles. However, overall, HH is almost equal to VV. For (downwind/upwind), HH is larger than VV at some angles, and the sea-spike phenomenon occurs. For (downwind/upwind), HH is larger than VV, even at some moderate incident angles, and this phenomenon is even more evident in the upwind. Therefore, in Figures 8 and 9, it can be inferred that the main reason of sea-spike occurrence is due to the multiple scattering of dihedral corner reflector between the breaking waves and the rough sea surface (i.e. Figure 5), which will enhance the Brewster effect, and the VV multipath is greatly attenuated.
3.6. Scattering Echo Time Series
The time behavior of wave back-scattering RCS for different incident angles and HH polarization is shown in Figure 10. In Figure 10, the incidence angles of downwind and upwind are set as 30°, 45°, 60°, 85°, 86°, and 87°, respectively, and 64 time points of the model are selected, which describe the breaking waves gradually curling. The incidence wave frequencies are set to , , and , respectively. The time interval is taken as . The size of the rough sea surface is with wind speed (dielectric constants of sea wave are obtained with the Debye model).

(a)

(b)

(c)

(d)

(e)

(f)
In Figure 10, one can see that for increasing frequency, the back-scattering RCS changes with the angle (downwind/upwind) are more and more pronounced. For a certain time point, the back-scattering RCS of the model decreases gradually as angle increases in the downwind case. However, the change of the back-scattering RCS is not regular when the angle increases at the upwind incident. Moreover, for a certain angle, back-scattering RCS at each angle is gradually increasing as time goes by at the upwind incident. The main reason for this behavior is again the fact that the structure is similar to a dihedral corner reflector between the wave and the sea surface. The structure is gradually formed, and this enhances the Brewster effect.
4. Conclusions
In this work, a sea surface with single breaking waves which contains realistic details has been established. Moreover, the SBR-PTD method is firstly applied to study the EM scattering of breaking waves. The existence of superevents where the scattering of HH polarization is larger than that of VV polarization has been confirmed by studying the scattering field of the sea surface with the single breaking wave model. Those superevents are more likely to occur at LGA with upwind. Moreover, the strong scattering localization of the model is also analyzed. The scattering of the sea with multiple breaking waves will be studied in the future.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This study was supported by the National Natural Science Foundation of China (grant nos. 61871457, 61971002, and 61801362) and the Foundation for Innovative Research Groups of the National Natural Science Foundation of China (grant no. 61621005).