# Electromagnetic Scattering Analysis of the Sea Surface with Single Breaking Waves

Authors: Jingjing Wang; Lixin Guo; Yiwen Wei; Shuirong Chai; Ke Li; Anqi Wang
Journal: International Journal of Antennas and Propagation
 (2021)

---

## 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.


---

## Body

## 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. Figure1 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.Figure 1
Discrete two-dimensional height-field representation of the water surface, the ground bottom, and the horizontal water velocity.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. Letz=hx be height of the water surface and let z=bx be the height of the ground. Let dx=hx−bx is the water depth and ux is the horizontal velocity of a vertical column of water. Waves where d>L/2 are called deep water and where d<L/2 are called shallow water waves [17]. In this, d is the water depth. L is the wavelength. The shallow water equations following from the above assumptions [18, 19] may be written as follows:(1)∂u∂t+u∂u∂y+g∂h∂y=0,(2)∂d∂t+∂∂yud=0,where g is the gravitational acceleration. Equation (1) expresses Newton’s law F=ma, 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 h. The resulting equations are then(3)∂u∂t+g∂h∂y=0,(4)∂h∂t+d∂u∂y=0.If we differentiate equation (3) with respect to y, then differentiate equation (4) with respect to t, and finally substitute for the cross-derivatives, we end up with(5)∂2h∂t2=gd∂2h∂y2,which is the one-dimensional wave equation with wave velocity gd.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 is(6)∂hi∂t=di−1+di2Δyui−1−di+di+12Δyui,(7)∂ui∂t=−ghi+1−hiΔy,where Δy is the separation of the samples along the y direction. Putting the above two equations together, we get(8)∂2hi∂t2=−gdi−1+di2Δy2hi−hi−1+gdi+di+12Δy2hi+1−hi,which provides the sought discrete approximation to equation (5).For simplicity, we use a first-order implicit method to solve ordinary differential equation (8). Let hn to denote h at the nth iteration and let dots denote differentiation with time. Then, the first-order implicit equations may be written as(9)hn−hn−1Δt=h˙n,h˙n−h˙n−1Δt=h¨n.Notice that the right-hand sides of these equations are evaluated at timen which corresponds to the end of the iteration rather than time n−1 which corresponds to the beginning of the iteration. This is what makes the iteration implicit and stable. Rearranging the above expression, we arrive at(10)hn=2hn−1−hn−2+Δt2h¨n.By substituting equation (8) into equation (10), we obtain(11)hin=2hin−1−hin−2−gΔt2di−1+di2Δy2hin−hi−1n+gΔt2di+di+12Δy2hi+1n−hin.After these manipulations, equation (11) is still a nonlinear equation. In order to find its solution, it needs to be linearized, i.e., d is regarded as a constant, and the wave velocity is fixed as a function of y. In this linear regime, the next value of h can be calculated from previous values with the symmetric tridiagonal linear system(12)Ahin=2hin−1−hin−2,where the matrix A is given by(13)A=e0f0f0e1f1f1e2⋱⋱⋱⋱⋱en−3fn−3fn−3en−2fn−2fn−2en−1,and the elements of A are as follows:(14)e0=1+gΔt2d0+d12Δy2,ei=1+gΔt2di−1+2di+di+12Δy2,0<i<n−1,en−1=1+gΔt2dk−2+dn−12Δy2,fi=−gΔt2di+di+12Δy2.Now, let us change the equation to(15)Ahin=hin−1+1−τ×hin−1−hin−2,where τ is phenomenologically introduced to describe damping phenomena. If τ=0, 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 hi<bi for some index i. 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 h and b vectors in order and testing whether hi<bi. 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 ofh with respect to y is replaced by the Laplacian.(16)∂2h∂t2=gd∂2h∂x2+∂2h∂y2=gd∇2h.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 y, and the other is independent of x. 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 equation(17)∂2h∂t2=gd∂2h∂y2,and in the second subiteration, we solve the equation(18)∂2h∂t2=gd∂2h∂x2.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 matricesRx, Ry, and Rz about x, y, and z axes are defined as(19)Rx=1000cosϕsinϕ0−sinϕcosϕ,Ry=cosθ0−sinθ010sinθ0cosθ,Rz=cosφsinφ0−sinφcosφ0001,where ϕ, θ, and φ are angles with x, y, and z axes, respectively.(20)Rbe=Rx×Ry×Rz=cosθcosφsinφcosθ−sinθsinϕcosφsinθsinφsinϕsinθ+cosϕcosφsinϕcosθsinθcosϕcosφ+sinφsinϕsinθcosϕsinφ−sinϕcosφcosϕcosθ.Scaling matrices are defined as(21)Rs=xs000ys000zs,where, xs, ys, and zs are scaling factors.World space matrix is defined as(22)M=xscosθcosφyssinφcosθ−zssinθ0xssinϕcosφsinθyssinφsinϕsinθ+cosϕcosφzssinϕcosθ0xssinθcosϕcosφ+sinφsinϕyssinθcosϕsinφ−sinϕcosφzscosϕcosθ0xtytzt1,where xt, yt, and zt are coordinates of the breaking waves. Let us also define M′ as the inverse of M and the vector Pt=xiyizi1, where xi, yi, and zi are position coordinates of the rough sea surface. The link between position coordinates and local-coordinate system is a rotation(23)locPt=Pt×M′.Let us now define(24)nPt=nxlocPt2nz,where nx=2−L×locPt1, nz=0.3−0.1×width×locPt3, L is length of the breaking waves, and width is width of the breaking waves.According to [20, 21], the relationship between wave height Hm and wind speed Um/s is(25)H=17.03⋅exp−α218.3269+2.361⋅exp−β28.8646,where α=U2/3−12.6549 and β=U2/3−6.4463.Wave steepnessS [22] is defined as the ratio of H/L, i.e.,(26)S=HL,where S is the wave steepness.We can get the wavelengthLm according to (25) and (26), i.e.,(27)L=17.03⋅exp−α2/18.3269S+2.361⋅exp−β2/8.8646S.The smoothing factor is defined as(28)cr=exp−nPtLenamplitude,where amplitude = 1. If nPtLen<0, then nPtLen=0, whereas if nPtLen>1e9, then nPtLen=1e9.Convert positional coordinates from local to global coordinates, i.e.,(29)nrPt=nPt×curlM=nxcoscr+nzsincrlocPt2nzcoscr−nxsincr,where curlM=coscr0−sincr010sincr0coscr.To transform local coordinates into global coordinates, we needM with 4 rows and 4 columns, and hence let us define(30)nrPt1=nxcoscr+nzsincrlocPt2nzcoscr−nxsincr1.Therefore, the coordinates of the transformed point position are(31)bwc=nrPt1×M,where bwc are the coordinates of the sea surface with single breaking waves.The sea surface with single breaking waves is shown in Figure2. 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 Δt=0.05s. The size of rough sea surface is 5m×5m with wind speed U=5m/s. S=1/3, 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.Figure 2
The sea surface with single breaking waves at different times. (a)t=△t, (b) t=2△t, (c) t=3△t, (d) t=4△t, (e) t=5△t, (f) t=6△t, (g) t=7△t, (h) t=8△t, (i) t=9△t, (j) t=10△t, (k) t=11△t, and (l) t=12△t.
(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 Figure3, 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.Figure 3
A ray from triangle patch 0 shot toward a target.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].(32)E⟶total=∑E⟶PO+E⟶PTD1+E⟶PTD2+E⟶PTD3.The electric field, approximated by PO integral, may be the written as follows:(33)E⟶PO=jk4πexp−jkRR∬η⋅s^×s^×J⟶expjks^⋅r⟶′ds′,where k=ωμ0ε01/2,η=ε0/μ01/2, s^ is the unit vector aligned along the scattering direction, R is the distance from the specular point r⟶' to the observation point, and J⟶ is the current density on the triangle surface. The diffraction field from each edge of a triangle patch is approximated by the PTD method:(34)E⟶PTD=jk4πexp−jkRR∫ls^×s^×t^Ie+s^×t^Imexpjks^⋅r⟶′dt′,where Ie and Im are as follows:(35)Ie=j2E⟶inc⋅t^ksin2βiDePTD+j2ηH⟶inc⋅t^ksinβiDemPTD,Im=j2ηH⟶inc⋅t^ksinβisinβsDmPTD,where E⟶inc and H⟶inc denote the incident electric field strength of patches’ edges; βi/s stands for the angle between the edge and incident or scattered ray; and DePTD, DmPTD, and DemPTD 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.


## 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. Figure1 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.Figure 1
Discrete two-dimensional height-field representation of the water surface, the ground bottom, and the horizontal water velocity.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. Letz=hx be height of the water surface and let z=bx be the height of the ground. Let dx=hx−bx is the water depth and ux is the horizontal velocity of a vertical column of water. Waves where d>L/2 are called deep water and where d<L/2 are called shallow water waves [17]. In this, d is the water depth. L is the wavelength. The shallow water equations following from the above assumptions [18, 19] may be written as follows:(1)∂u∂t+u∂u∂y+g∂h∂y=0,(2)∂d∂t+∂∂yud=0,where g is the gravitational acceleration. Equation (1) expresses Newton’s law F=ma, 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 h. The resulting equations are then(3)∂u∂t+g∂h∂y=0,(4)∂h∂t+d∂u∂y=0.If we differentiate equation (3) with respect to y, then differentiate equation (4) with respect to t, and finally substitute for the cross-derivatives, we end up with(5)∂2h∂t2=gd∂2h∂y2,which is the one-dimensional wave equation with wave velocity gd.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 is(6)∂hi∂t=di−1+di2Δyui−1−di+di+12Δyui,(7)∂ui∂t=−ghi+1−hiΔy,where Δy is the separation of the samples along the y direction. Putting the above two equations together, we get(8)∂2hi∂t2=−gdi−1+di2Δy2hi−hi−1+gdi+di+12Δy2hi+1−hi,which provides the sought discrete approximation to equation (5).For simplicity, we use a first-order implicit method to solve ordinary differential equation (8). Let hn to denote h at the nth iteration and let dots denote differentiation with time. Then, the first-order implicit equations may be written as(9)hn−hn−1Δt=h˙n,h˙n−h˙n−1Δt=h¨n.Notice that the right-hand sides of these equations are evaluated at timen which corresponds to the end of the iteration rather than time n−1 which corresponds to the beginning of the iteration. This is what makes the iteration implicit and stable. Rearranging the above expression, we arrive at(10)hn=2hn−1−hn−2+Δt2h¨n.By substituting equation (8) into equation (10), we obtain(11)hin=2hin−1−hin−2−gΔt2di−1+di2Δy2hin−hi−1n+gΔt2di+di+12Δy2hi+1n−hin.After these manipulations, equation (11) is still a nonlinear equation. In order to find its solution, it needs to be linearized, i.e., d is regarded as a constant, and the wave velocity is fixed as a function of y. In this linear regime, the next value of h can be calculated from previous values with the symmetric tridiagonal linear system(12)Ahin=2hin−1−hin−2,where the matrix A is given by(13)A=e0f0f0e1f1f1e2⋱⋱⋱⋱⋱en−3fn−3fn−3en−2fn−2fn−2en−1,and the elements of A are as follows:(14)e0=1+gΔt2d0+d12Δy2,ei=1+gΔt2di−1+2di+di+12Δy2,0<i<n−1,en−1=1+gΔt2dk−2+dn−12Δy2,fi=−gΔt2di+di+12Δy2.Now, let us change the equation to(15)Ahin=hin−1+1−τ×hin−1−hin−2,where τ is phenomenologically introduced to describe damping phenomena. If τ=0, 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 hi<bi for some index i. 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 h and b vectors in order and testing whether hi<bi. 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 ofh with respect to y is replaced by the Laplacian.(16)∂2h∂t2=gd∂2h∂x2+∂2h∂y2=gd∇2h.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 y, and the other is independent of x. 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 equation(17)∂2h∂t2=gd∂2h∂y2,and in the second subiteration, we solve the equation(18)∂2h∂t2=gd∂2h∂x2.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 matricesRx, Ry, and Rz about x, y, and z axes are defined as(19)Rx=1000cosϕsinϕ0−sinϕcosϕ,Ry=cosθ0−sinθ010sinθ0cosθ,Rz=cosφsinφ0−sinφcosφ0001,where ϕ, θ, and φ are angles with x, y, and z axes, respectively.(20)Rbe=Rx×Ry×Rz=cosθcosφsinφcosθ−sinθsinϕcosφsinθsinφsinϕsinθ+cosϕcosφsinϕcosθsinθcosϕcosφ+sinφsinϕsinθcosϕsinφ−sinϕcosφcosϕcosθ.Scaling matrices are defined as(21)Rs=xs000ys000zs,where, xs, ys, and zs are scaling factors.World space matrix is defined as(22)M=xscosθcosφyssinφcosθ−zssinθ0xssinϕcosφsinθyssinφsinϕsinθ+cosϕcosφzssinϕcosθ0xssinθcosϕcosφ+sinφsinϕyssinθcosϕsinφ−sinϕcosφzscosϕcosθ0xtytzt1,where xt, yt, and zt are coordinates of the breaking waves. Let us also define M′ as the inverse of M and the vector Pt=xiyizi1, where xi, yi, and zi are position coordinates of the rough sea surface. The link between position coordinates and local-coordinate system is a rotation(23)locPt=Pt×M′.Let us now define(24)nPt=nxlocPt2nz,where nx=2−L×locPt1, nz=0.3−0.1×width×locPt3, L is length of the breaking waves, and width is width of the breaking waves.According to [20, 21], the relationship between wave height Hm and wind speed Um/s is(25)H=17.03⋅exp−α218.3269+2.361⋅exp−β28.8646,where α=U2/3−12.6549 and β=U2/3−6.4463.Wave steepnessS [22] is defined as the ratio of H/L, i.e.,(26)S=HL,where S is the wave steepness.We can get the wavelengthLm according to (25) and (26), i.e.,(27)L=17.03⋅exp−α2/18.3269S+2.361⋅exp−β2/8.8646S.The smoothing factor is defined as(28)cr=exp−nPtLenamplitude,where amplitude = 1. If nPtLen<0, then nPtLen=0, whereas if nPtLen>1e9, then nPtLen=1e9.Convert positional coordinates from local to global coordinates, i.e.,(29)nrPt=nPt×curlM=nxcoscr+nzsincrlocPt2nzcoscr−nxsincr,where curlM=coscr0−sincr010sincr0coscr.To transform local coordinates into global coordinates, we needM with 4 rows and 4 columns, and hence let us define(30)nrPt1=nxcoscr+nzsincrlocPt2nzcoscr−nxsincr1.Therefore, the coordinates of the transformed point position are(31)bwc=nrPt1×M,where bwc are the coordinates of the sea surface with single breaking waves.The sea surface with single breaking waves is shown in Figure2. 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 Δt=0.05s. The size of rough sea surface is 5m×5m with wind speed U=5m/s. S=1/3, 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.Figure 2
The sea surface with single breaking waves at different times. (a)t=△t, (b) t=2△t, (c) t=3△t, (d) t=4△t, (e) t=5△t, (f) t=6△t, (g) t=7△t, (h) t=8△t, (i) t=9△t, (j) t=10△t, (k) t=11△t, and (l) t=12△t.
(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 Figure3, 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.Figure 3
A ray from triangle patch 0 shot toward a target.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].(32)E⟶total=∑E⟶PO+E⟶PTD1+E⟶PTD2+E⟶PTD3.The electric field, approximated by PO integral, may be the written as follows:(33)E⟶PO=jk4πexp−jkRR∬η⋅s^×s^×J⟶expjks^⋅r⟶′ds′,where k=ωμ0ε01/2,η=ε0/μ01/2, s^ is the unit vector aligned along the scattering direction, R is the distance from the specular point r⟶' to the observation point, and J⟶ is the current density on the triangle surface. The diffraction field from each edge of a triangle patch is approximated by the PTD method:(34)E⟶PTD=jk4πexp−jkRR∫ls^×s^×t^Ie+s^×t^Imexpjks^⋅r⟶′dt′,where Ie and Im are as follows:(35)Ie=j2E⟶inc⋅t^ksin2βiDePTD+j2ηH⟶inc⋅t^ksinβiDemPTD,Im=j2ηH⟶inc⋅t^ksinβisinβsDmPTD,where E⟶inc and H⟶inc denote the incident electric field strength of patches’ edges; βi/s stands for the angle between the edge and incident or scattered ray; and DePTD, DmPTD, and DemPTD 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 int=△t stage and t=11△t stage is shown in Figure 4. The frequency of the incident wave is f=3GHz. Therefore, the dielectric parameter is ε=70.6059,39.7343 from the Debye model [25]. The incident angles are θi=−90°∼90°, φi=0°. The size of rough sea surface is 1m×1m with wind speed U=5m/s. 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 θi>0° represents downwind and θi<0° represents upwind, and this is true for all of the following examples.Figure 4
Back-scattering RCS of the sea surface with single breaking waves at different times. (a)t=△t HH, (b) t=△t VV, (c) t=11△t HH, (d) t=11△t VV.
(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 Figure5. A comparison of the bistatic scattering coefficient with and without breaking waves is shown in Figure 6. In this, the t=12△t stage of the sea surface with single breaking waves is selected. The frequency of the incident wave is f=10GHz. Therefore, the dielectric parameter is ε=55.8866,37.6777 from the Debye model. The incident angles are θi=30°, φi=0°. The size of rough sea surface is 5m×5m with wind speed U=5m/s. 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.Figure 5
A locally broken wave resembling a dihedral corner reflector.Figure 6
A comparison of the bistatic scattering coefficient with and without breaking waves. (a) HH. (b) VV.
(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,t=△t, t=5△t, t=8△t, and t=11△t, 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 t=△t, t=5△t, t=8△t, and t=11△t stages is shown in Figure 7. The frequency of the incident wave is f=15GHz. Therefore, the dielectric parameter is ε=43.8890,39.1149 from the Debye model. The incident angles are θi=−90°∼90°, φi=0°. The size of rough sea surface is 5m×5m with wind speed U=5m/s.Figure 7
Back-scattering RCS of the sea surface with single breaking waves as obtained with SBR-PTD. (a)t=△t, (b) t=5△t, (c) t=8△t, (d) t=11△t.
(a)(b)(c)(d)In Figure7, for t=△t 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 t=5△t, t=8△t, and t=11△t stages of the sea surface with single breaking waves, the plume forms. HH is larger than VV at the grazing angles of less than 30°, 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 thet=12△t stage of the sea surface with single breaking waves is shown in Figure 8. The frequency of the incident wave is f=10GHz and the corresponding dielectric parameter is ε=55.8866,37.6777. The size of rough sea surface is 5m×5m and the wind speed U=5m/s. The incident angles are 30° (downwind/upwind), 85° (downwind), and 87° (upwind), respectively.Figure 8
Back-scattering field intensity distribution at different incidence angles of each triangle patch on the sea surface with single breaking waves in thet=12△t stage. (a) θi=30° (downwind). (b) θi=−30° (upwind). (c) θi=85° (downwind). (d) θi=−87° (upwind).
(a)(b)(c)(d)Forθi=30° (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 θi=85°/87° (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 thet=12△t stage of the sea surface with single breaking waves is shown in Figure 9. The frequency of incident wave is f=10GHz, and the dielectric parameter ε=55.8866,37.6777 from the Debye model. The size of the rough sea surface is 5m×5m with wind speed U=5m/s. Incidence angles are 25° (downwind/upwind), 45° (downwind/upwind), and 87° (downwind/upwind), respectively.Figure 9
Influence of incidence angle on EM scattering of sea surface with single breaking waves in thet=12△t stage. (a) θi=25° (downwind). (b) θi=−25° (upwind). (c) θi=45° (downwind). (d) θi=−45° (upwind). (e) θi=87° (downwind). (f) θi=−87° (upwind).
(a)(b)(c)(d)(e)(f)In Figure9, for θi=25° (downwind/upwind), HH is larger than VV at some angles. However, overall, HH is almost equal to VV. For θi=45° (downwind/upwind), HH is larger than VV at some angles, and the sea-spike phenomenon occurs. For θi=87° (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 Figure10. 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 f=1.5GHz, f=3GHz, and f=10GHz, respectively. The time interval is taken as Δt=0.05s. The size of the rough sea surface is 5m×5m with wind speed U=5m/s (dielectric constants of sea wave are obtained with the Debye model).Figure 10
Back-scattering RCS varies with time at different incident angles. (a) 1.5 GHz (downwind). (b) 1.5 GHz (upwind). (c) 3 GHz (downwind). (d) 3 GHz (upwind). (e) 10 GHz (downwind). (f) 10 GHz (upwind).
(a)(b)(c)(d)(e)(f)In Figure10, 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.


## 3.1. Verification of the Effectiveness of the SBR-PTD Method

Back-scattering RCS of the sea surface with single breaking waves both int=△t stage and t=11△t stage is shown in Figure 4. The frequency of the incident wave is f=3GHz. Therefore, the dielectric parameter is ε=70.6059,39.7343 from the Debye model [25]. The incident angles are θi=−90°∼90°, φi=0°. The size of rough sea surface is 1m×1m with wind speed U=5m/s. 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 θi>0° represents downwind and θi<0° represents upwind, and this is true for all of the following examples.Figure 4
Back-scattering RCS of the sea surface with single breaking waves at different times. (a)t=△t HH, (b) t=△t VV, (c) t=11△t HH, (d) t=11△t VV.
(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 Figure5. A comparison of the bistatic scattering coefficient with and without breaking waves is shown in Figure 6. In this, the t=12△t stage of the sea surface with single breaking waves is selected. The frequency of the incident wave is f=10GHz. Therefore, the dielectric parameter is ε=55.8866,37.6777 from the Debye model. The incident angles are θi=30°, φi=0°. The size of rough sea surface is 5m×5m with wind speed U=5m/s. 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.Figure 5
A locally broken wave resembling a dihedral corner reflector.Figure 6
A comparison of the bistatic scattering coefficient with and without breaking waves. (a) HH. (b) VV.
(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,t=△t, t=5△t, t=8△t, and t=11△t, 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 t=△t, t=5△t, t=8△t, and t=11△t stages is shown in Figure 7. The frequency of the incident wave is f=15GHz. Therefore, the dielectric parameter is ε=43.8890,39.1149 from the Debye model. The incident angles are θi=−90°∼90°, φi=0°. The size of rough sea surface is 5m×5m with wind speed U=5m/s.Figure 7
Back-scattering RCS of the sea surface with single breaking waves as obtained with SBR-PTD. (a)t=△t, (b) t=5△t, (c) t=8△t, (d) t=11△t.
(a)(b)(c)(d)In Figure7, for t=△t 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 t=5△t, t=8△t, and t=11△t stages of the sea surface with single breaking waves, the plume forms. HH is larger than VV at the grazing angles of less than 30°, 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 thet=12△t stage of the sea surface with single breaking waves is shown in Figure 8. The frequency of the incident wave is f=10GHz and the corresponding dielectric parameter is ε=55.8866,37.6777. The size of rough sea surface is 5m×5m and the wind speed U=5m/s. The incident angles are 30° (downwind/upwind), 85° (downwind), and 87° (upwind), respectively.Figure 8
Back-scattering field intensity distribution at different incidence angles of each triangle patch on the sea surface with single breaking waves in thet=12△t stage. (a) θi=30° (downwind). (b) θi=−30° (upwind). (c) θi=85° (downwind). (d) θi=−87° (upwind).
(a)(b)(c)(d)Forθi=30° (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 θi=85°/87° (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 thet=12△t stage of the sea surface with single breaking waves is shown in Figure 9. The frequency of incident wave is f=10GHz, and the dielectric parameter ε=55.8866,37.6777 from the Debye model. The size of the rough sea surface is 5m×5m with wind speed U=5m/s. Incidence angles are 25° (downwind/upwind), 45° (downwind/upwind), and 87° (downwind/upwind), respectively.Figure 9
Influence of incidence angle on EM scattering of sea surface with single breaking waves in thet=12△t stage. (a) θi=25° (downwind). (b) θi=−25° (upwind). (c) θi=45° (downwind). (d) θi=−45° (upwind). (e) θi=87° (downwind). (f) θi=−87° (upwind).
(a)(b)(c)(d)(e)(f)In Figure9, for θi=25° (downwind/upwind), HH is larger than VV at some angles. However, overall, HH is almost equal to VV. For θi=45° (downwind/upwind), HH is larger than VV at some angles, and the sea-spike phenomenon occurs. For θi=87° (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 Figure10. 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 f=1.5GHz, f=3GHz, and f=10GHz, respectively. The time interval is taken as Δt=0.05s. The size of the rough sea surface is 5m×5m with wind speed U=5m/s (dielectric constants of sea wave are obtained with the Debye model).Figure 10
Back-scattering RCS varies with time at different incident angles. (a) 1.5 GHz (downwind). (b) 1.5 GHz (upwind). (c) 3 GHz (downwind). (d) 3 GHz (upwind). (e) 10 GHz (downwind). (f) 10 GHz (upwind).
(a)(b)(c)(d)(e)(f)In Figure10, 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.



---

*Source: 1545031-2021-11-28.xml*