Front. Mar. Sci. Frontiers in Marine Science Front. Mar. Sci. 2296-7745 Frontiers Media S.A. 10.3389/fmars.2025.1611884 Marine Science Original Research Numerical modeling of electromagnetic wave propagation in spatially-varying evaporation duct conditions via 3D parabolic equation method Ji Hanjie 1 2 Guo Lixin 1 Zhang Yan 1 Nie Tianhang 1 Wei Yiwen 1 * Zhang Jinpeng 2 * Li Qingliang 2 Guo Xiangming 2 Zhang Yusheng 2 1 School of Physics, Xidian University, Xi’an, China 2 National Key Laboratory of Electromagnetic Environment, China Research Institute of Radiowave Propagation, Qingdao, China

Edited by: Jian Wang, Tianjin University, China

Reviewed by: Zi He, Nanjing University of Science and Technology, China

Cheng Yang, Tianjin University, China

*Correspondence: Yiwen Wei, ywwei@xidian.edu.cn; Jinpeng Zhang, zhangjp@crirp.ac.cn

06 06 2025 2025 12 1611884 15 04 2025 21 05 2025 Copyright © 2025 Ji, Guo, Zhang, Nie, Wei, Zhang, Li, Guo and Zhang 2025 Ji, Guo, Zhang, Nie, Wei, Zhang, Li, Guo and Zhang

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

This study numerically investigates electromagnetic (EM) wave propagation in spatially-varying evaporation ducts over rough sea surfaces. Conventional two-dimensional (2D) models assume homogeneous refractive index distribution along the cross-range dimension in a single propagation plane, limiting their ability to capture the 3D spatial heterogeneities present in real-world scenarios. Under significant horizontal gradient variations in evaporation ducts, EM wave propagation effects across the cross-range dimension become significant. We investigate an advanced 3D parabolic equation (3DPE) framework that synergistically integrates anisotropic refractive profiles with sea-surface roughness characterization. An even-odd splitting Fourier transform algorithm enables efficient computational analysis of EM wave propagation across azimuthal planes. Quantitative analysis reveals that the 3DPE framework delivers over 40% performance improvement compared to the 2D model. This approach significantly enhances predictive accuracy for over-the-horizon radar assessments in maritime environments, providing crucial support for optimizing next-generation communication systems.

 electromagnetic wave propagation spatially-varying evaporation ducts 3D parabolic equation rough sea surface numerical modeling section-in-acceptance Ocean Observation Introduction

Evaporation ducts (Zhang et al., 2011) are meteorological structures formed by marine surface evaporation combined with vertical temperature/humidity gradients. Typically confined to the lower atmospheric boundary layer near the ocean surface, these ducts exhibit characteristic heights between 0-40m (Zhao et al., 2021). The modified refractive index profile (M-profile) governs their electromagnetic (EM) trapping capability, enabling over-the-horizon (OTH) propagation through vertical gradient constraints (Wang et al., 2023a). When the M-profile gradient satisfies trapping criterion, EM waves in specific bands (e.g., S/C/X) undergo multipath reflections between sea surface and duct upper boundary (Yang et al., 2022b), creating divergent propagation trajectories that overcome Earth’s curvature limitations (Zhao et al., 2009). This energy confinement simultaneously degrades radar detection in supra-duct region, generating surveillance blind zones (Anderson, 1995). Duct-induced propagation paths exhibit marked deviations from standard atmospheric refraction patterns, inducing systematic errors in radar range, altitude, and velocity measurements (Douchin et al., 1994).

Understanding wave propagation characteristics in evaporation duct environments is crucial for optimizing maritime radio system designs (Yang et al., 2024). Since the 1980s, the parabolic equation (PE) methodology has become the dominant approach due to its concise mathematical formulation and superior capability in characterizing EM propagation in ducting environments (Hardin and Tappert, 1973). established the standard PE framework through Taylor series approximation, extending angular validity to 15° elevation while introducing the split-step Fourier transform (SSFT) algorithm—still the most computationally efficient technique for large-scale EM propagation modeling. Subsequent advancements extended PE applications to rough sea surface environments (Kuttler and Janaswamy, 2002). developed numerical stabilization techniques to address computational instabilities in sea surface modeling (Thomson and Quach, 2005). systematically investigated wave propagation in Arctic conditions, while (Fabbro et al., 2006) expanded PE methodology for rough sea surface (Guo et al., 2023). further advanced this field through comparative analysis of sea surface roughness effects.

Traditional two-dimensional (2D) PE (2DPE) models exhibit critical limitations in resolving EM wave propagation in complex marine environments due to oversimplified spatial dimensionality and boundary conditions. Restricted to a single vertical propagation plane (e.g., range-altitude plane), these models unrealistically assume uniform cross-range refractive index distributions—an assumption directly contradicted by the 3D spatial heterogeneity of M-profile fields induced by air-sea interactions. Furthermore, 2D models inherently neglect lateral reflection and scattering effects caused by realistic sea surface boundaries. These limitations highlight the necessity for advanced 3D models capable of resolving interactions between spatially-varying refractivity gradients and anisotropic sea surfaces. The 3D PE (3DPE) framework overcomes these constraints by establishing a full Cartesian coordinate system that incorporates lateral dimensions. This advancement enables comprehensive characterization of 3D propagation dynamics, especially under scenarios with significant lateral refractive index gradients. When simulating wave propagation over rough sea surfaces with lateral heterogeneity, the model effectively addresses boundary truncation issues through transverse boundary condition implementation.

This study presents a 3DPE-based numerical framework for simulating EM wave propagation in spatially-varying evaporation ducts over rough sea surfaces. By integrating meteorological reanalysis data with the Naval Postgraduate School (NPS) model, we establish a spatially heterogeneous M-profile field and generate rough sea surface boundary conditions through sea wave spectrum theory. The SSFT algorithm is applied iteratively during numerical solving, while an even-odd decomposition method transforms governing equations into strict 2D Fourier forms, enabling efficient computation via fast Fourier transform (FFT) techniques. Systematic validation using shipborne radar-measured propagation loss (PL) demonstrates the 3DPE framework’s computational accuracy improvement of more than 40% over conventional 2D models, establishing a multiphysics-coupled simulation tool for OTH radar performance evaluation.

 Methodology 3D vector PE: solutions, initial field, and boundary conditions

Assuming the EM wave propagates along the positive x-axis and considering only forward wave propagation, with the time-harmonic factor e − i ω t , the electric field vector E → and the magnetic field vector H → in the Cartesian coordinate system (Levy, 2000). For an arbitrary scalar field component Π , the ansatz Π = Π e + Π m is postulated where Π e and Π m correspond to the scalar potentials of transverse magnetic (TM) and transverse electric (TE) modes respectively. This scalar field inherently satisfies the homogeneous scalar wave equation in the source-free region, expressed as ∇ 2 Π + k 0 2 Π = 0 (Harrington, 2001). By substituting Π ( x , y , z ) = e − i k 0 x Π ( x , y , z ) into ∇ 2 Π + k 0 2 Π = 0 , we can derive (Equation 1):

∂ 2 Π ∂ x 2 + 2 i k 0 ∂ Π ∂ x + ∇ t 2 Π = 0

where ∇ t 2 denotes the transverse Laplace operator and k 0 is the free space propagation constant. Under uniform atmospheric distribution conditions ∂ n / ∂ x ≈ 0 , the governing wave equation decomposes into forward and backward propagation operators. For forward-propagating EM waves, the Feit-Fleck approximation method (Feit and Fleck, 1978) constructs a pseudo-differential operator Q = 1 + ∇ t 2 / k 0 2 to reformulate 3D wave propagation, expressed as Equation 2:

∂ Π ( x , y , z ) ∂ x = i k 0 [ 1 + 1 k 0 2 ∇ t 2 − 1 ] Π ( x , y , z )

The 3DPE numerical solution is an iterative advancement process: wavefronts in the ( y , z ) plane propagate stepwise along the x-axis, with each subsequent wavefront configuration determined algorithmically from its predecessor via operator splitting. In scenarios considering only reflected waves above the surface, during each iteration of the 3DPE implementation, there exists (Equation 3):

Π ( x + Δ x , y , z ) = e − i k 0 Δ x · ℑ 2 − 1 [ e i k x Δ x Π ˜ ˜ ( x 0 , k y , k z ) ]

where the x-component of the wave vector k 0 is defined by k x = k 0 2 − k y 2 − k z 2 , and ℑ 2 − 1 is the 2D inverse Fourier transform operator. The analytical solution of (Equation 6) is impeded by intrinsic complexity of differential operator 1 + ∇ t 2 / k 02 , requiring numerical approximation methods for implementation. This study implements the SSFT algorithm for 3DPE solutions, mirroring the computational framework used in 2DPE approaches. To balance computational efficiency and algorithmic simplicity, 2D FFT techniques are routinely applied (Janaswamy, 2001). When ground creeping waves are neglected and only reflected wave effects are considered, the wave function Π ˜ ˜ ( x , k y , k z ) at each propagation step reduces to a superposition of direct and reflected wave components (Janaswamy, 2003):

Π ˜ ˜ ( x , k y , k z ) = ∫ − ∞ + ∞ ∫ − ∞ + ∞ Π ( x , y ′ , z ′ ) e − i k y y ′ · [ e − i k z z ′ + Γ ( k z ) e i k z z ′ ] d y ′ d z ′

where Γ ( k z ) denotes the surface reflection coefficient. To accelerate computations via FFT in (Equation 4), this study employs a field parity decomposition method where Π is decomposed into odd and even components. This transformation converts (Equation 4) into a rigorous 2D Fourier transform expression across y-z coordinates, with the decomposed wave function expressed as (Hu et al., 2005):

{ Π o d d ( x , y , z ) = { Π ( x , y , z ) , z > 0 − Π ( x , y , − z ) , z < 0 Π e v e n ( x , y , z ) = { Π ( x , y , z ) , z > 0 Π ( x , y , − z ) , z < 0

By substituting (Equation 5) and the functional relationships between Π and its even/odd components into (Equation 4), simplification of the resulting integral yields the final expression:

Π ˜ ˜ e v e n , o d d ( x , k y , k z ) = ∫ − ∞ + ∞ ∫ − ∞ + ∞ Π e v e n , o d d ( x , y ′ , z ′ ) e − i ( k y y ′ + k z z ′ ) d y ′ d z ′

where Π ˜ ˜ e v e n and Π ˜ ˜ o d d denote the 2D Fourier transforms even/odd components, respectively.

Establishing initial field conditions is crucial for 3D EM propagation solutions using the 3DPE method. This approach utilizes the Fourier transform relationship between antenna radiation patterns and the model’s initial configuration. By applying inverse Fourier transforms to the antenna’s directional pattern A ( p ) , the initial field distribution u ( 0 , z ) can be systematically derived, as shown in Equation 7:

u ( 0 , z ) = ℑ − 1 [ U ( 0 , p ) ] = k 0 2 π e i π / 4 ∫ − ∞ + ∞ A ( p ) ( k 0 2 − p 2 ) 1 / 4 e i p z d p

where U ( 0 , p ) = 2 π k 0 e i π / 4 A ( p ) / ( k 0 2 − p 2 ) 1 / 4 . The Gaussian antenna pattern, commonly used in EM propagation modeling. This model accurately models practical radiation characteristics in the paraxial region, mathematically expressed as A ( p ) = e − p 2 ω 2 / 4 , where ω = 2 ln 2 / sin ( θ b w / 2 ) , θ b w denotes the half-power beamwidth (Barrios, 1994). This equation defines the initial free-space field distribution. For practical EM propagation analysis, the image theory must transform this distribution into upper half-space. The derived relationship follows (Barrios, 1994):

u g ( 0 , z ) = u ( 0 , z − z 0 ) + Γ H / V u ( 0 , z + z 0 )

Here, z 0 represents the transmitting antenna’s height, and Γ H / Γ V denote the Fresnel reflection coefficient for horizontal/vertical polarization (Ray et al., 2019). For Equation 8, the transformation in Equations 9 and 10 can be derived by applying the FFT’s shift property:

U ( 0 , p ) = U ( 0 , p ) e − i p z 0 + Γ H / V U ( 0 , − p ) e i p z 0 u g ( 0 , z ) = 1 2 π ∫ − ∞ + ∞ [ U ( 0 , p ) e − i p z 0 + Γ H / V U ( 0 , − p ) e i p z 0 ] e i p z d p

In the SSFT algorithm, to truncate the computational domain, the 3DPE model must implement absorbing boundary conditions. Physically, these boundary conditions aim to fully absorb EM waves reaching the ± y and ± z boundaries of the computational domain. This study implements a Cosine-taper window function along both ± y and ± z boundaries. The window configuration in Equation 11 is derived from the even-odd decomposition methodology and the sequential arrangement of FFT sampling points N:

w ( i ) = { 0.5 + 0.5 cos [ ( 8 π ( i − 1 − 3 N / 8 ) / N ) ] , 1 + 3 N / 8 ≤ i ≤ 1 + N / 2 0.5 + 0.5 cos [ ( 8 π ( i − 2 − 5 N / 8 ) / N ) ] , 2 + N / 2 ≤ i ≤ 2 + 5 N / 8 1 , others

This study investigates numerical modeling of 3DPE under lossy boundary conditions, focusing on EM wave propagation over rough sea surfaces where seawater’s intrinsic material properties induce significant EM dissipation. Let n → ( n x , n y , n z ) be the unit vector of the outer normal of the boundary surface, and Δ g be the normalized impedance. Then we have: n → × E → = Δ g η 0 · n → × ( n → × H → ) . For the horizontal plane boundary, n → = z ^ , that is, n x = n y =0 and n z = 1 . From this, we can obtain:

E y = Δ g η 0 H x ,   E x = − Δ g η 0 H y

In the vectorial 3DPE model, any EM field component arises from the coupled contributions of both Π e and Π m , as described by the unified formulation:

E x , y , z = E x , y , z e + E x , y , z m ,   H x , y , z = H x , y , z e + H x , y , z m

For horizontal planar boundaries, the general boundary condition formulated in Equation 12 can be systematically transformed into a set of simplified interface constraints in Equation 14 by invoking the constitutive relationships defined in Equation 13:

∂ Π e ∂ z + i k 0 Δ g Π e = 0 ,   ∂ Π m ∂ z + i k 0 Δ g Π m = 0 3D EM-PL modeling

Tropospheric EM wave PL consists of two components: free-space loss L f from natural spherical wave spreading and medium-induced loss L A from wave absorption, scattering, and reflection in ducting environments (Zhang, 2012). Total PL is expressed as their summation, as shown in Equation 15:

L = L f + L A = 32.45 + 20 lg f + 20 lg r − 20 lg F

where f is frequency, r represents the distance between transmitting and receiving antennas, and F corresponds to the propagation factor (PF). In the 3DPE model, the PF is expressed as (Harrington, 2001) (Equation 16):

F = | Π ( x , y , z ) Π y − f r e e s p a c e ( x , y , z ) |

where Π ( x , y , z ) represents the recursive field solution from 3DPE marching computations, and Π y − f r e e s p a c e ( x , y , z ) denotes the free-space distribution of the y-directed magnetic field component H y , specifically formulated as Equation 17:

Π y − f r e e s p a c e ( x , y , z ) = i k 0 e − i k 0 r 4 π r sin θ cos φ e − ( k 0 cos θ ) 2 σ z 2 / 2

In the formula, R = ( z − h t ) 2 + x 2 + y 2 denotes the transmitter-receiver distance, where h t is the transmit altitude, sin θ = x 2 + y 2 / R , cos φ = x / x 2 + y 2 , and x, y represent the coordinates of the receiver’s projection onto the ( x , y ) plane.

 NPS evaporation duct model

The NPS model (Babin and Dockery, 2002), the U.S. Navy’s operational evaporation duct prediction framework, integrates the Coupled Ocean-Atmosphere Response Experiment (COARE) bulk flux algorithm with empirical marine experiment data. This synthesis represents a substantial advancement over the foundational Paulus–Jeske model (Paulus, 1985) (Babin and Dockery, 2002). rigorously validated the NPS model against in situ measurements, benchmarking it against three alternatives: the Naval Warfare Assessment model, the Naval Research Laboratory model, and the Babin–Young–Carton (BYC) model (Babin et al., 1997). Results showed that both NPS and BYC models achieve optimal refractivity profile estimation, with NPS exhibiting superior algorithmic stability. Given its proven reliability across diverse marine environments, this study employs the NPS model for evaporation duct computations.

The NPS model derives M-profiles through vertical modified refractivity profiling using five near-surface parameters–air temperature (AT), air pressure (AP), wind speed (WS), relative humidity (RH), sea surface temperature (SST). These parameters, obtainable synchronously or asynchronously above sea level, are processed through the model’s Monin-Obukhov Similarity (MOS) theory-based boundary layer solver (Babin et al., 1997). Simultaneously, AT and specific humidity (SH) profiles are calculated using MOS theory-derived stability functions, as shown in Equations 18 and 19: (Wang et al., 2023b):

T ( z ) = T 0 + θ * κ [ ln ( z z 0 θ ) − Ψ t ( ζ ) ] − Γ d z q ( z ) = q 0 + q * κ [ ln ( z z 0 θ ) − Ψ t ( ζ ) ]

where T 0 and q 0 are the AT and SH of the sea surface, respectively; ζ = z / L , L is the Monin-Obukhov length; θ ∗ and q ∗ are the related MOS scaling parameters of potential temperature (PT) and SH, respectively; z 0 θ and z 0 q are PT and SH roughness lengths, respectively; Γ d is dry adiabatic lapse rate; κ is the von Karman constant; and Ψ t is the temperature correction function (Yang and Wang, 2022).

The original NPS model employed COARE2.6 bulk flux algorithm to compute MOS scaling parameters (Fairall et al., 1996). This study adopts the improved COARE3.0 algorithm (Fairall et al., 2003) for more accurate computations in complex air-sea conditions. Unlike conventional evaporation duct models that rely on land-based correction functions-which perform poorly over marine surfaces-the COARE3.0 model incorporates maritime observation-derived temperature correction functions for neutral atmospheric conditions, as shown in Equation 20:

Ψ t = − ξ − 2 3 [ ξ − 5 0.35 ] exp ( − c ) − 2 3 × 5 0.35

where c = min ( 50 , 0.35 ξ ) . For unstable conditions, the correction function accounts for the convective-limit scenario, as shown in Equations 21-25:

Ψ t = Ψ t k + ξ 2 Ψ k t 1 + ξ 2 Ψ t k = 2 ln [ 1 + z p t 2 ] z p t = ( 1 − 15 ξ ) 0.5 Ψ k t = 1.5 ln [ z p g 2 + z p g + 1 3 ] − 3 arctan [ 2 z p g + 1 3 ] + π 3 z p g = ( 1 − 34.15 ξ ) 0.333

The COARE3.0 algorithm implements Grachev’s method (Grachev et al., 1997) for stability parameter determination, achieving significant computational efficiency improvements, as shown in Equation 26:

ξ = { C R i b 1 − α R i b , R i b ≥ 0 C R i b ( 1 + R i b R i b c ) − 1 , R i b < 0

where α ≈ 5 is an empirical coefficient, and C is obtained from the overall flux conversion coefficient. R i b and R i b c are the bulk Richardson number and the saturation Richardson number respectively.

The COARE3.0 algorithm quantifies wind effects on sea surface roughness through its aerodynamic roughness formulation, as shown in Equations 27 and 28:

z 0 = a u ∗ 2 g + 0.11 ν u ∗ 2 a = { 0.11 U ≤ 10 0.11 + 0.007 ( U 10 − 10 ) 8 10 < U ≤ 18 0.018 U > 18

where v denotes the sea surface viscosity coefficient, and U 10 represents the WS at 10m above the sea surface. The roughness lengths for temperature and humidity are set equal, taken as Equations 29 and 30:

z 0 t = z 0 q = min ( 1.1 × 10 − 4 , 5.5 × 10 − 5 R r − 0.6 ) R r = u ∗ 2 z 0 v

where u ∗ is friction velocity.

The original NPS model employs the temperature correction function proposed by Beljaars and Holtslag (B&H) for stable atmospheric stratification. However, empirical analyses reveal that while the B&H function performs well under weakly stable conditions, it significantly overestimates evaporation duct heights (EDHs) in strongly stable regimes. To address this limitation, the correction function developed by Gorbachev et al.—a widely validated approach derived from observational data in the Surface Heat Budget of the Arctic Ocean Experiment (SHEBA). The revised functional form is expressed as follows based on SHEBA-derived parameters (Grachev et al., 2007), as shown in Equation 31:

Ψ h   S H E B A ( ζ ) = − b h 2 ln ( 1 + c h ζ + ζ 2 ) + ( − a h B h + b h c h 2 B h ) × ( ln 2 ζ + c h B h 2 ζ + c h + B h − ln c h − B h c h + B h )

where a h ≡ b h = 5 , c h = 3 , and B h = c h 2 − 4 = 5 .

Quantifying the near-surface atmospheric M-profile requires vertical profiles of AP and water vapor pressure (WVP). The AP profile utilizes combined solutions from the ideal gas law and hydrostatic equations (Babin and Dockery, 2002), while the WVP profile involves thermodynamic integration of precalculated SH profile, as shown in Equations 32 and 33.

P ( z 2 ) = P ( z 1 ) exp ( g ( z 1 − z 2 ) R T ν ¯ ) e ( z ) = q ( z ) P ( z ) ϵ + ( 1 − ϵ ) q ( z )

where T ν ¯ is the mean value of the virtual temperature at height z 1 and z 2 , and ϵ is a constant with a value of 6.22. Atmospheric refraction arises from vertical gradient variations in atmospheric properties, governed by meteorological parameters including AT, AP, and WVP. To address the geometric distortions caused by Earth’s curvature in EM wave propagation, a simplified flat Earth model is commonly used, which approximates Earth’s spherical surface as a flat plane. Consequently, the modified atmospheric refractivity can be formulated as follows (Hitney et al., 1985):

M ( z ) = 77.6 T × ( P + 4810 e T ) + 0.157 z

The marine atmospheric M-profile is computed using Equation 34 based on vertical profiles of AT, AP, and WVP obtained from Equations 18-33. Using the NPS model with five predefined meteorological parameters, we calculated two sets of evaporation duct refractive index profiles ( Figure 1 ). As can be seen from the figure, the evaporation duct profiles are different for different environmental parameters, which indicates that the evaporation duct variations are closely related to the environment. The EDH corresponding to the minimum modified refractive indices of 341.3 M-units and 389.3 M-units were determined as 13.5 m and 27.3 m, respectively (Zhang et al., 2011).

The M-profiles for EDHs of (a) 13.5m and (b) 27.3m, respectively.

 Pierson-Moskowitz wave spectrum

The sea wave spectrum (Guerin and Johnson, 2015) mathematically characterizes the energy distribution of sea surface undulations across frequency and direction domains. This spectral representation quantifies the variance of sea surface displacement as a spectral density function in wavenumber space (spatial frequency domain). Through Fourier transformation, it decomposes stochastic sea surfaces into constituent wave components with distinct frequencies and propagation directions. In practical applications, the 2D ocean wave spectrum is typically expressed as the product of a directional spreading function D ( k ) and a wavenumber spectrum S ( k ) , as shown in Equation 35:

S ( k ) = 1 k S ( k ) D ( k )

where k = ( k x , k y ) constitutes a wavenumber vector. The wavenumber k = ∥ k ∥ = k x 2 + k y 2 quantifies spatial wave frequency, defined as the number of complete wave cycles per unit length. The PM spectrum (Pierson and Moskowitz, 1964), empirically derived from 1964 North Atlantic wave measurements, serves as the canonical model for energy distribution in fully-developed seas under constant wind forcing. This spectral formulation establishes quantitative relationships between wave energy density and spectral components in equilibrium wind-wave systems. For modeling surface waves, the equation of the PM spectrum is defined as Equation 36:

S P M ( k ) = α 2 k 3 · exp ( − β g 2 k 2 U 19.5 4 )

where α =0.0081 and β =0.74 are dimensionless constants, g denotes gravitational acceleration, and U 19.5 denotes the WS value at an elevation of 19.5 m above the sea surface.

After obtaining the PM spectral density, a Monte Carlo methodology can be used to simulate 2D rough sea surfaces (Meng et al., 2024). This method models the sea surface as a superposition of stochastic-amplitude harmonic components with stochastic amplitudes following independent Gaussian distributions. In the frequency domain, Gaussian white noise is filtered through a PM spectral model incorporating wind-wave dynamics, which imposes wave number-dependent energy constraints. An inverse Fourier transform then converts these filtered spectral components into spatially correlated sea surface height fluctuations that statistically match the target PM spectrum. For a 2D random rough surface with spatial dimensions L x = M · Δ x and L y = N · Δ y , where M and N denote the number of grid points along the x- and y-axes, respectively, Δ x and Δ y represent the grid resolutions, the height f ( x m , y n ) at each surface coordinate ( x m = m Δ x , y n = n Δ y ) can be formulated as Equation 37:

f ( x m , y n ) = 1 L x L y ∑ m k = − M / 2 + 1 M / 2 ∑ n k = − N / 2 + 1 N / 2 F ( k m k , k n k ) · exp [ i ( k m k x m + k n k y n ) ]

The Fourier transform of the discrete function is denoted as F ( k m k , k n k ) , where k m and k n denote the spatial frequency components along the x- and y-axes. Its explicit mathematical form is defined by Equation 38:

F ( k m k , k n k ) = 2 π [ L x L y S ( k m k , k n k ) ] 1 / 2 · { [ N ( 0 , 1 ) + i N ( 0 , 1 ) ] 2 , m k ≠ 0 , M / 2 ,   a n d n k ≠ 0 , N / 2 N ( 0 , 1 ) , m k ≠ 0 , M / 2 ,   a n d n k ≠ 0 , N / 2

Figure 2 presents PM spectral model-generated 2D sea surface profiles under varying WSs (3–10 m/s). Analysis shows a strong correlation between wind forcing and sea surface dynamics: lower wind velocities produce mild surface undulations, while increased WSs trigger nonlinear height fluctuation amplification, reflecting enhanced wave-turbulence interactions.

Rough sea surfaces generated using the PM spectrum under varying WS conditions.

 Construction of spatially-varying evaporation duct refractivity field using meteorological reanalysis data

High-fidelity modeling of spatially-varying marine evaporation duct environments requires high-resolution meteorological parameter fields as fundamental inputs. Meteorological reanalysis data, particularly through their multiscale coupling and physically consistent boundary layer parameterization schemes, offer an effective technical approach for deriving accurate spatial distributions of meteorological elements (Ji et al., 2024a). ERA5 (Ji et al., 2024b), the fifth-generation global atmospheric reanalysis dataset published by the European Centre for Medium-Range Weather Forecasts, employs an advanced data assimilation framework that synergizes multi-source observations (satellite retrievals, surface stations, radiosondes) with numerical weather models. This integration generates spatiotemporally continuous meteorological fields with enhanced fidelity compared to its predecessor ERA-Interim (Yang et al., 2022a).

The South China Sea (SCS, 4°N–23°N, 105°E–118°E), a semi-enclosed tropical marginal sea influenced by complex atmospheric dynamics, maintains year-round high temperatures and humidity due to the interplay of meteorological and topographical factors. Using ERA5 reanalysis data, we extract six key meteorological parameters for evaporation duct modeling: AT, AP, zonal/meridional wind velocities, dewpoint temperature (DT), and SST. Table 1 details these parameters with their corresponding vertical levels. Figure 3 illustrates their spatial distribution across the SCS on March 1, 2025 at 00:00 UTC, revealing significant spatial heterogeneity in all meteorological parameters between different geographic locations.

Retrieved parameters along with corresponding reanalysis height levels.

Meteorological element Reanalysis height Unit
AT 2 m °C
DT 2 m °C
zonal wind velocity 10 m m/s
meridional wind velocity 10 m m/s
AP Surface hPa
SST Surface °C

Spatial distribution of meteorological parameters in the SCS.

The WS is calculated by vectorially summing the zonal and meridional wind velocities, as shown in Equation 39:

W S = u 2 + v 2

where u and v denote the zonal and meridional wind velocity components, respectively. RH is calculated using the Magnus-Tetens approximation (Alduchov and Eskridge, 1996), which employs AT and DT to determine saturation vapor pressure (Lawrence, 2005), as shown in Equations 40-42:

e s ( T ) = 6.112 · exp ( 17.62 · T T + 243.12 ) e ( T d ) = 6.112 · exp ( 17.62 · T d T d + 243.12 ) R H = e ( T d ) e s ( T ) × 100 %

where T and T d denote AT and DT in degrees Celsius, respectively. This study incorporates derived WS and RH parameters along with AT, AP, and SST into the NPS model to simulate evaporation ducts.

Figure 4 illustrates the spatial distribution of evaporation duct refractivity calculated using the NPS model with reanalysis meteorological data from Figure 3 . The projection of the study area onto longitude-latitude coordinates is simultaneously displayed. Significant spatial heterogeneity in refractivity is observed due to variations in meteorological parameters across geographic locations. For this investigation, we select an EM emission source at 111°E, 19.75°N. This location was explicitly configured as the radiation origin with EM wave propagation azimuthally fixed at 90° (true east). We focus on the evaporation duct along this eastern bearing where spatial refractivity heterogeneity is identified. Subsequent phases will implement numerical modeling of ducting effects along this bearing, with specific focus on quantifying PL anomalies under non-standard atmospheric conditions.

3D evaporation duct refractivity field and its horizontal projection in longitude-latitude plane.

 Computational modeling of EM propagation in evaporation duct environments

We perform EM wave PL simulations using a 3DPE framework, incorporating both the computed 3D refractivity field (constructed along the EM emission direction) and a numerically generated rough sea surface under 10 m/s WS condition. Beyond spatially-varying ducts, we simulate 3D PL characteristics under both standard atmospheric conditions and uniform evaporation duct environments, where the latter employs the refractive index profile shown in Figure 1 . The simulation framework incorporates a radar system operating at 10 GHz with the following configuration: antenna altitude of 5 m, horizontal polarization, elevation angle of 0°, and azimuthal half-power beamwidth of 15°.

Numerical simulations model EM wave propagation over smooth sea surfaces ( Figure 5 ), revealing distinct characteristics between standard atmospheric and evaporation duct environments. Under conventional atmospheric refraction ( Figure 5A ), EM attenuation stems primarily from free-space spreading and medium absorption, showing exponential signal decay with distance. In contrast, evaporation ducts create vertically stratified refractive gradients that enable anomalous propagation. This ducting effect confines EM energy within the atmospheric layer, reducing dielectric losses and permitting OTH signal transmission with non-exponential attenuation patterns. While standard environments cause rapid signal divergence, duct propagation preserves signal integrity beyond standard refractive limits, though partial energy leakage still occurs at duct boundaries. As EDH increases ( Figures 5B, C ), the duct structure’s EM wave trapping capability intensifies. Elevated EDH configurations demonstrate enhanced confinement of EM waves, enabling lower-loss propagation with extended transmission distances. This improved confinement mechanism allows sustained wave guidance while significantly reducing energy dissipation over extended propagation paths.

3D PL characterization over smooth sea surfaces in: (A) standard atmosphere, (B) uniform evaporation duct with 13.5m EDH, (C) uniform evaporation duct with 27.3m EDH, and (D) spatially-varying evaporation duct environments.

The propagation characteristics in spatially-varying duct media ( Figure 5D ) demonstrate fundamentally distinct behaviors compared to those within homogeneous evaporation ducts. This divergence primarily stems from the inherent structural inhomogeneity of the medium, which induces multiscale interactions including localized reflection, refraction, and scattering phenomena. The presence of spatially varying refractive index profiles across different regions causes alterations in phase velocity and ray path divergence during propagation. Such refractive index gradients generate complex interference patterns through multipath propagation mechanisms, thereby introducing heightened complexity in predicting signal propagation trajectories and energy distribution characteristics. These anisotropic propagation effects significantly deviate from the deterministic path confinement observed in uniform ducting structures.

Experimental observations demonstrate that EM wave propagation over rough sea surfaces produces striped shadow patterns in PL field distributions ( Figure 6 ). Compared to smooth sea conditions, rough sea surfaces exhibit accelerated attenuation rates and expanded blind zones due to surface roughness effects. This phenomenon primarily results from wave crest-induced shadow effects, which amplify multipath interference through enhanced reflection, scattering, and attenuation processes. Significant PL intensification occurs under conditions of pronounced sea surface height variations and active wave motion, substantially compromising transmission reliability.

3D PL characterization over rough sea surfaces in: (A) standard atmosphere, (B) uniform evaporation duct with 13.5m EDH, (C) uniform evaporation duct with 27.3m EDH, and (D) spatially-varying evaporation duct environments.

 Empirical PL observation-based performance validation of 3DPE framework

The China Research Institute of Radiowave Propagation (CRIRP) executed an EM wave OTH propagation observation in the SCS ( Figure 7 ), using the research vessel Qiongsha 3, equipped with an S-band radar system (Ji et al., 2024b). Key system parameters are summarized in Table 2 . During the 15-hour maritime campaign between Wenchang City (19°33’N, 110°49’E) and Yongxing Island (16°84’N, 112°33’E) as shown in Figure 7 , the mobile platform transmitted S-band radar signals while shore-based stations along Hainan coast continuously received OTH signals. The coastal reception system incorporated a low-noise amplifier (LNA) in its front-end circuitry to optimize signal-to-noise ratio during prolonged monitoring. The EM PL is quantified by rigorously applying radar equations to OTH signals received during experimental configurations (Guo et al., 2019), as shown in Equation 43:

EM wave OTH propagation observation campaign: (A) vessel navigation trajectory, (B) Qiongsha 3 research vessel, (C) S-band radar system and auxiliary experimental instruments.

Primary operational parameters.

Configuration Value
Frequency 3.1 GHz
Transmitting Power 41.8 dBm
Transmitting Antenna Gain 28.0 dB
Bandwidth 5.0 MHz
Pulse Width 3.0 μs
Polarization HH
L o b s = P t + G t + G r + G L N A − P r − L t − L r

where P t and P r denote the transmit power and received power, respectively; G t and G r represent the antenna gains at the transmitter and receiver ends, respectively; while L t and L r indicate the feeder line losses at the transmitter and receiver sides, respectively.

During the Qiongsha 3 research vessel’s transect from Wenchang City to Yongxing Island, we deployed a multi-sensor evaporation duct monitoring system (EDMS) to collect high-resolution meteorological and oceanographic data (Zhang et al., 2020). This integrated system incorporated meteorological sensors and infrared SST radiometers, continuously measuring five key variables—including AT, RH, WS, AP, and SST. Following data acquisition, parameters underwent preprocessing for outlier removal and temporal alignment before being fed into the NPS model to generate modified refractivity profiles (Ji et al., 2024c). This methodology enabled quantitative characterization of spatially heterogeneous evaporation duct structures.

To rigorously validate the 3DE framework’s efficacy and precision, we conduct systematic verification experiments using the 2DPE propagation model. Shipborne meteorological parameters are processed to generate non-uniform refractive index profiles and integrated into both computational frameworks, ensuring atmospheric condition consistency. Validation strictly followed Table 2 parameters, with results compared against the 3DPE model implementation. The 2DPE model configuration utilizes a Gaussian antenna profile for initial field distribution calculation to accurately simulate practical EM emission scenarios. This implementation employs a 2D window function for upper boundary processing, while enhancing lower boundary conditions through a roughness attenuation factor that modifies Fresnel reflection coefficients at impedance surfaces. The Miller-Brown (Guo et al., 2023) formulation is adopted to quantify surface roughness characteristics.

The comparative analysis reveals distinct morphological differences between measured and simulated PL characteristics ( Figure 8 ). Field observations exhibit stochastic nonlinear fluctuations attributable to complex marine environmental factors, whereas parabolic equation method simulations under idealized assumptions generate smoothed attenuation curves. This constitutes the primary source of discrepancy between measured and simulated PLs. The established 3DPE framework demonstrates superior alignment with experimental data, accurately reproducing both the attenuation trend and numerical magnitudes. In contrast, 2DPE calculations show significant deviations in magnitude and fail to capture observed attenuation patterns, with discrepancies amplifying with transmission range. These findings empirically validate the 3DPE framework’s enhanced capability in modeling EM wave propagation within complex environments featuring anisotropic refractivity gradients and rough sea surface boundaries. In contrast, the 2DPE methodology employs simplified environmental assumptions (Zhang et al., 2021).

Comparison of PL values calculated by 3DPE framework and 2DPE with actual observation.

We quantitatively compare PL observations with two simulation results using three statistical metrics: root mean square error (RMSE), mean absolute error (MAE), and mean absolute percentage error (MAPE). These metrics provide multidimensional accuracy assessment: RMSE highlights larger deviations by squaring errors, MAE measures average error magnitude linearly, while MAPE calculates relative error percentages unaffected by absolute values. The corresponding formulas are defined as Equations 44–46. Quantitative evaluation demonstrates notable differences in prediction accuracy between the models. The 2DPE model shows significant PL prediction errors with metrics of 6.173 dB (RMSE), 5.424 dB (MAE), and 2.617% (MAPE). In comparison, the 3DPE model achieves substantially improved performance under identical test conditions, yielding 3.450 dB RMSE, 2.821 dB MAE, and 1.355% MAPE. Through comprehensive metric analysis, the 3DPE framework demonstrates significant advantages over traditional 2DPE method. Specifically, it achieves performance improvements of 54.301%, 47.990%, and 48.223% in RMSE, MAE, and MAPE respectively, with an overall error reduction exceeding 40%–quantitatively confirming the superiority of 3D EM modeling in complex propulsion system analysis.

R M S E = 1 n ∑ i = 1 n ( L o b s ( i ) − L s ( i ) ) 2 M A E = 1 n ∑ i = 1 n | L o b s ( i ) − L s ( i ) | M A P E = 100 % n ∑ i = 1 n | L s ( i ) − L o b s ( i ) L o b s ( i ) |

where L s denotes simulated PL, n specifies the number of observational samples.

 Conclusion

This study examines EM wave propagation through spatially-varying evaporation ducts above rough sea surfaces, focusing on addressing limitations of conventional 2D modeling approaches. Traditional 2D models impose restrictive homogeneity assumptions within propagation planes, failing to account for critical refractive index variations characteristic of real maritime environments. Our analysis reveals that the 2DPE model shows significant prediction errors in PL estimation, with RMSE, MAE, and MAPE attaining 6.173dB, 5.424dB, and 2.617% respectively. In contrast, the 3DPE model demonstrates substantially improved accuracy under identical conditions, achieving 3.450dB (RMSE), 2.821dB (MAE), and 1.355% (MAPE). This represents an overall error reduction exceeding 40% compared to the 2D baseline, conclusively validating 3D modeling’s superior capability in handling complex wave propagation phenomena. These advancements establish new capabilities for precise EM environmental characterization, paving the way for enhanced maritime communication systems.

 Data availability statement

The original contributions presented in the study are included in the article/supplementary material. Further inquiries can be directed to the corresponding authors.

 Author contributions

HJ: Funding acquisition, Writing – original draft. LG: Funding acquisition, Writing – original draft, Project administration, Conceptualization, Supervision. YZ: Software, Writing – original draft, Resources. TN: Writing – review & editing, Software. YW: Writing – review & editing, Funding acquisition. JZ: Writing – review & editing, Funding acquisition, Validation, Project administration. QL: Writing – review & editing, Formal Analysis. XG: Writing – review & editing, Investigation. YSZ: Writing – review & editing, Data curation.

 Funding

The author(s) declare that financial support was received for the research and/or publication of this article. This work was supported in part by the National Natural Science Foundation of China under Grant U21A20457, Grant 62231021, Grant 62371380, and Grant 62271457; and in part by the Stabled-Support Scientific Project of China Research Institute of Radiowave Propagation under Grant A132301214; and in part by the Fundamental Research Funds for the Central Universities under Grant YJSJ25020; and in part by the Innovation Fund of Xidian University.

 Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

 Generative AI statement

The author(s) declare that no Generative AI was used in the creation of this manuscript.

 Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

 References Alduchov O. A. Eskridge R. E. (1996). Improved Magnus form approximation of saturation vapor pressure. J. Appl. Meteorol. 35, 601–609. doi: 10.1175/1520-0450(1996)035<0601:IMFAOS>2.0.CO;2 Anderson K. D. (1995). Radar detection of low-altitude targets in a maritime environment. IEEE Trans. Antennas Propag. 43, 609–613. doi: 10.1109/8.387177 Babin S. M. Dockery G. D. (2002). LKB-Based evaporation duct model comparison with buoy data. J. Appl. Meteorol. 41, 434–446. doi: 10.1175/1520-0450(2002)041<0434:LBEDMC>2.0.CO;2 Babin S. M. Young G. S. Carton J. A. (1997). A new model of the oceanic evaporation duct. J. Appl. Meteorol. 36, 193–204. doi: 10.1175/1520-0450(1997)036<0193:ANMOTO>2.0.CO;2 Barrios A. (1994). A terrain parabolic equation model for propagation in the troposphere. IEEE Trans. Antennas Propag. 42, 90–98. doi: 10.1109/8.272306 Douchin N. Bolioli S. Christophe F. Combes P. (1994). Theoretical study of the evaporation duct effects on satellite-to-ship radio links near the horizon. IEE Proc. - Microwaves Antennas Propagation 141, 272–278. doi: 10.1049/ip-map:19941168 Fabbro V. Bourlier C. Combes P. F. (2006). Forward propagation modeling above gaussian rough surfaces by the parabolic shadowing effect. Prog. Electromagnetics Res. 58, 243–269. doi: 10.2528/PIER05090101 Fairall C. W. Bradley E. F. Hare J. E. Grachev A. A. Edson J. B. (2003). Bulk parameterization of air-sea fluxes: Updates and verification for the COARE algorithm. J. Clim. 16, 571–591. doi: 10.1175/1520-0442(2003)016(0571:BPOASF)2.0.CO;2 Fairall C. W. Bradley E. F. Rogers D. P. Edson J. B. Young G. S. (1996). Bulk parameterization of air-sea fluxes for Tropical Ocean-Global Atmosphere Coupled-Ocean Atmosphere Response Experiment. J. Geophys. Res. 101, 3747–3764. doi: 10.1029/95JC03205 Feit M. D. Fleck J. A. (1978). Light propagation in graded-index optical fibers. Appl. Optics. 17, 3990–3998. doi: 10.1364/AO.17.003990 Grachev A. A. Andreas E. L. Fairall C. W. Guest P. S. Persson P. O. G. (2007). SHEBA flux-profile relationships in the stable atmospheric boundary layer. Bound.-Layer Meteor. 124, 315–333. doi: 10.1007/s10546-007-9177-6 Grachev A. A. Fairall C. W. (1997). Dependence of the Monin-Obukhov stability parameter on the bulk Richardson number over the ocean. J. Appl. Meteor. 36, 406–414. doi: 10.1175/1520-0450(1997)036<0406:DOTMOS>2.0.CO;2 Guerin C. Johnson J. T. (2015). A simplified formulation for rough surface cross-polarized backscattering under the second-order small-slope approximation. IEEE Trans. Geosci. Remote Sensing. 53, 6308–6314. doi: 10.1109/TGRS.2015.2440443 Guo X. Li Q. Zhao Q. Kang S. Wei Y. Yang L. (2023). A comparative study of rough sea surface and evaporation duct models on radio wave propagation. IEEE Trans. Antennas Propag. 71, 6060–6071. doi: 10.1109/TAP.2023.3276571 Guo X. Zhao D. Zhang L. Wang H. Kang S. Lin L. (2019). C band transhorizon signal characterisations in evaporation duct propagation environment over Bohai Sea of China. IET Microwaves Antennas Propagation 13, 407–413. doi: 10.1049/iet-map.2018.5040 Hardin R. H. Tappert F. D. (1973). Applications of the split-step Fourier method to the numerical solution of nonlinear and variable coefficient wave equations. SIAM Rev. 15, 423. Harrington R. F. (2001). Time-harmonic electromagnetic fields (John Wiley & Sons). Hitney H. V. Richter J. H. Pappert R. A. Anderson K. D. Baumgartner G. B. (1985). Tropospheric radio propagation assessment. Proc. IEEE. Wiley-IEEE Press 73, 265–283. doi: 10.1109/PROC.1985.13138 Hu H. Chai S. Mao J. (2005). “ Field strength prediction behind a lossy dielectric block by using the 3D parabolic equation model ,” in 2005 Asia-Pacific Microwave Conference Proceedings. (Suzhou, China: IEEE) 1–4. doi: 10.1109/APMC.2005.1606416 Janaswamy R. (2001). Radio wave propagation over a nonconstant immittance plane. Radio Sci. 36, 387–405. doi: 10.1029/2000RS002338 Janaswamy R. (2003). Path loss predictions in the presence of buildings on flat terrain: a 3-D vector parabolic equation approach. IEEE Trans. Antennas Propag. 51, 1716–1728. doi: 10.1109/TAP.2003.815415 Ji H. Guo L. Zhang J. Wei Y. Guo X. Zhang Y. (2024a). EDH-STNet: An evaporation duct height spatiotemporal prediction model based on Swin-Unet integrating multiple environmental information sources. Remote Sens. 16, 4227. doi: 10.3390/rs16224227 Ji H. Guo L. Zhang J. Wei Y. Guo X. Zhang Y. . (2024b). IPILT–OHPL: An over-the-horizon propagation loss prediction model established by incorporating prior information into the LSTM–Transformer structure. IEEE J. Sel. Top. Appl. Earth Observ. Remote Sens. 17, 10067–10082. doi: 10.1109/JSTARS.2024.3395630 Ji H. Wei Y. Zhao Q. Guo L. Zhang J. (2024c). EKDInformer-LTEDH: an Informer-based environmental knowledge-driven prediction model for long-term evaporation duct height. IEEE Geosci. Remote Sens. Lett. 21, 1–5. doi: 10.1109/LGRS.2024.3468290 Kuttler J. R. Janaswamy R. (2002). Improved Fourier transform methods for solving the parabolic wave equation. Radio Sci. 37, 1–5. doi: 10.1029/2001RS002488 Lawrence M. G. (2005). The relationship between relative humidity and the dewpoint temperature in moist air: A simple conversion and applications. Bull. Amer. Meteorol. Soc 86, 225–233. doi: 10.1175/BAMS-86-2-225 Levy M. (2000). Parabolic equation methods for electromagnetic wave propagation (IET Digital Library: IEE Electromagnetic Wave Series). Meng X. Su J. Li K. Dong C. Hou G. Liu Y. . (2024). An effective composite scattering model for the sea surface with a target based on the TSBR-TSM algorithm. IEEE Antennas Wirel. Propag. Lett. 23, 2411–2415. doi: 10.1109/LAWP.2024.3393773 Paulus R. A. (1985). Practical application of an evaporation duct model. Radio Sci. 20, 887–896. doi: 10.1029/RS020i004p00887 Pierson W. J. Moskowitz L. (1964). A proposed spectral form for fully developed wind seas based on the similarity theory of S. A. Kitaigorodskii. J. Geophys. Res. 69, 5181–5190. doi: 10.1029/JZ069i024p05181 Ray A. Ray M. Basuray A. (2019). Plane wave spectrum analysis of transverse discontinuities in dielectric waveguide I. Reflection coefficients. Opt. Eng. 58, 117103. doi: 10.1117/1.OE.58.11.117103 Thomson A. D. Quach T. D. (2005). Application of parabolic equation methods to HF propagation in an Arctic environment. IEEE Trans. Antennas Propag. 53, 412–419. doi: 10.1109/TAP.2004.838751 Wang S. Yang K. Shi Y. Yang F. Zhang H. Ma Y. (2023a). Prediction of over-the-horizon electromagnetic wave propagation in evaporation ducts based on the gated recurrent unit network model. IEEE Trans. Antennas Propag. 71, 3485–3496. doi: 10.1109/TAP.2023.3240998 Wang S. Yang K. Shi Y. Zhang H. Yang F. Hu D. . (2023b). Long-term over-the-horizon microwave channel measurements and statistical analysis in evaporation ducts over the Yellow Sea. Front. Mar. Sci. 10. doi: 10.3389/fmars.2023.1077470 Yang C. Shi Y. Wang J. Feng F. (2022a). Regional spatiotemporal statistical database of evaporation ducts over the South China Sea for future long-range radio application. IEEE J. Sel. Top. Appl. Earth Observ. Remote Sens. 15, 6432–6444. doi: 10.1109/JSTARS.2022.3197406 Yang C. Wang J. (2022). The investigation of cooperation diversity for communication exploiting evaporation ducts in the South China Sea. IEEE Trans. Antennas Propag. 70, 8337–8347. doi: 10.1109/TAP.2022.3177509 Yang C. Wang J. Ma J. (2022b). Exploration of X-band communication for maritime applications in the South China Sea. IEEE Antennas Wirel. Propag. Lett. 21, 481–485. doi: 10.1109/LAWP.2021.3136044 Yang F. Wang S. Yang K. Shu Y. Shi Y. (2024). Marine high-speed over-the-horizon communications and channel sensing in evaporation ducts over the South China Sea. IEEE Antennas Wirel. Propag. Lett. 23, 2870–2874. doi: 10.1109/LAWP.2024.3411057 Zhang H. Li Q. Guo X. Hou C. (2021). “Research on radio wave propagation prediction model using 3-D parabolic equation over rough sea surface,” in 2021 International Applied Computational Electromagnetics Society (ACES-China) Symposium. 1–4. doi: 10.23919/ACES-China52398.2021.9581721 Zhang J. Wu Z Zhang Y Wang B . (2012). Evaporation duct retrieval using changes in radar sea clutter power versus receiving height. Prog. In Electromagnetics Res. 126, 555–571. doi: 10.2528/PIER11121307 Zhang J. Wu Z. Zhu Q. Wang B. (2011). A four-parameter m-profile model for the evaporation duct estimation from radar clutter. Prog. In Electromagnetics Res. 114, 353–368. doi: 10.2528/PIER11012204 Zhang J. Zhang Y. Xu X. Li Q. Wu J. (2020). Estimation of sea clutter inherent Doppler spectrum from shipborne S-band radar sea echo. Chin. Phys. B. 29, 68402. doi: 10.1088/1674-1056/ab888a Zhao X. Huang J. Gong S. (2009). Modeling on multi-eigenpath channel in marine atmospheric duct. Radio Sci. 44, 106–110. doi: 10.1029/2008RS003847 Zhao W. Zhao J. Li J. Zhao D. Huang L. Zhu J. . (2021). An evaporation duct height prediction model based on a long short-term memory neural network. IEEE Trans. Antennas Propag. 69, 7795–7804. doi: 10.1109/TAP.2021.3076478