Development of a Numerical Prediction Model for Marine Lower Atmospheric Ducts and Its Evaluation across the South China Sea

Abstract

The Coupled Ocean–Atmosphere–Wave–Sediment Transport (COAWST) model serves as the foundation for creating a forecast model to detect lower atmospheric ducts in this study. A set of prediction tests with different forecasting times focusing on the South China Sea domain was conducted to evaluate the short-term forecasting effectiveness of lower atmospheric ducts. The assessment of sounding observation data revealed that the prediction model performed well in predicting the characteristics of all types of ducts. The mean values of the forecasting errors were slightly lower than the reanalysis data but had lower levels of correlation coefficients. At an altitude of about 2000 m, the forecasted error of modified atmospheric refractivity reached peak values and then decreased gradually with increasing altitude. The accuracy of forecasted surface ducts was higher than that of elevated ducts. Noticeable land–sea differences were identified for the spatial distributions of duct characteristics, and the occurrence rates of both the surface and elevated ducts were high at sea. As for the differences among the forecasts of 24, 48, and 72 h ahead, the differences primarily occurred at altitude levels below 20 m and 500 m~1500 m, which are consistent with the differences in the duct height.

1. Introduction

The atmospheric duct phenomenon is a unique weather event that can substantially impact the propagation process of electromagnetic waves. Atmospheric ducting leads to the confinement of electromagnetic waves in a specific area because of the rapid vertical changes in the atmospheric refractivity, as if these waves are propagating within a duct, giving rise to its name [1]. According to the characteristics of vertical profiles, atmospheric ducts are generally divided into evaporation, surface, and elevated ducts. Evaporation ducts mostly occur in the offshore atmosphere below a few tens of meters, and the humidity decreases sharply with altitude due to the evaporation of water vapor. Surface ducts occur in the boundary layer atmosphere, with the lower boundary connected to the surface. Elevated ducts mostly occur in the lower troposphere atmosphere, with the lower boundary hanging in the air.

Atmospheric duct events exhibit significant effects on radio communications and detections. In communication, atmospheric ducts may generate large-scale interferences from signals with the same frequencies, thereby posing challenges for local terminals to access base stations and triggering communication issues [2]. In radar detection, atmospheric ducts may also lead to aberrant detections beyond the range and occurrences of blind spots for radars [1]. In the marine environment, atmospheric ducting events at sea are more frequent, and the conditions requisite for different types of ducts differ in a non-trivial way [3,4]. Specifically, evaporation ducts only occur over water and may be affected by the horizontal homogeneity of the sea surface. The surface and elevated ducts, however, are defined more by the surface layer characteristics and overlying atmospheric columns, yielding the vertical gradient of refractivity. Thus, the efficient prediction of maritime atmospheric duct events will be crucial in ship navigation, communication security, and radar detection.

A previous investigation [5] suggested that the atmospheric duct phenomenon is essentially the aberrant distributions of atmospheric temperature and humidity vertical profiles, leading to drastic alterations in the atmospheric refractivity with height. The generation of atmospheric ducts can be inferred based on vertical changes in atmospheric refractivity gradients. The following equation can be used to calculate the atmospheric refractivity as the curvature of the propagation path in air for electromagnetic waves [5]:

$$N=\frac{77.6}{T}\left(P+\frac{4810e}{T}\right),$$

(1)

where N represents the refractivity (N-units); T represents the air temperature (K); P represents the air pressure (hPa); and e represents the water vapor pressure, which can be calculated by the following equation [6]:

$$e=\frac{qP}{\epsilon +\left(1-\epsilon \right)q},$$

(2)

where q represents the specific humidity (g/kg) and $\epsilon$ is a constant (0.622).

For ground devices of the electromagnetic wave transceiver, the modified atmospheric refractivity, which considers the effects of Earth’s curvature, has been applied more widely. Its formula is expressed as follows:

$$M=N+\frac{h}{R}\times {10}^6=N+0.157h,$$

(3)

where M represents the modified atmospheric refractivity (M-units); R represents the radius of the Earth (6,371,000 m); and h represents the height (m).

The occurrence of atmospheric duct events can be further diagnosed based on the outputs of numerical prediction models. At present, several atmospheric duct prediction models based on numerical weather models have been established in previous studies, and many simulation studies focusing on model performance evaluation and spatiotemporal variations of atmospheric ducts have been conducted [6,7,8,9,10].

For example, Atkinson et al. [11] created a duct prediction model based on the fifth-generation Pennsylvania State University National Center for Atmospheric Research Mesoscale Model (MM5) and performed simulation evaluations over the Persian Gulf area. Their evaluation test demonstrated that the simulated atmospheric ducts under strong wind conditions are more consistent with the observations after assimilating the aircraft observations in the boundary layer. Haack et al. [12] assessed the duct simulations of comparative experiments conducted in Wallops Island using four different models: the Coupled Ocean–Atmosphere Mesoscale Prediction System (COAMPS), Met Office’s Unified Model (MetUM), MM5, and Global Environmental Multiscale Model (GEM). They discovered that regardless of any model, there was a certain degree of systematic error. Moreover, the frequency and strength of the simulated atmospheric ducts were less and weaker than those observed in the offshore atmospheric boundary layer. Chai et al. [13] applied the LightGBM evaporation duct height prediction model (LGB-PHY) to perform evaporation duct simulation tests over the South China Sea. Their findings revealed that the root mean square error (RMSE) value of the evaporation duct height decreased by 68%, and the R-square value increased by 6.5% compared with the Extreme Gradient Boosting (XGB) model.

With the progress in duct research, researchers have found that the phenomenon of atmospheric ducts at sea is not only related to changes in atmospheric stratifications but also closely related to alterations in the underlying sea surface [14,15,16]. The dynamic changes at the sea surface, including the generation and propagation of ocean surface waves and eddies, as well as the distribution of ocean surface temperature, influence the distributions of winds, turbulence, temperature, and humidity in the lower atmosphere via energy and water vapor exchanges between the atmosphere and the sea, thereby affecting the occurrence and characteristics of atmospheric ducts. In this way, it is crucial to comprehensively consider the processes of air–sea interactions for simulating and predicting the atmospheric duct events at sea. A previous study developed a prediction system based on the Coupled Ocean–Atmosphere–Wave–Sediment Transport (COAWST) model while comprehensively considering the coupling processes of atmosphere, ocean, and waves [17]. Compared with an atmosphere-only model, COAWST can provide a more accurate and comprehensive simulation of marine meteorological and hydrological elements [15], enhancing the reliability of atmospheric duct simulations at sea.

Due to the low height of the evaporation ducts, which ranges from a few meters to several tens of meters [18], the evaporation ducts will not be discussed in this study. The changes only observed in the surface and elevated ducts (commonly referred to as lower atmospheric ducts) will be the research objects in this study. Here, the authors will establish a diagnostic prediction model for low atmospheric ducts at sea based on the COAWST model and conduct forecasting tests over the South China Sea area. The South China Sea is a marginal sea in the Western Pacific Ocean, with frequent shipping and maritime operational activities. Atmospheric ducting over the South China Sea is believed to play a crucial role in radar detection and radio communication, among others. The researchers attempt to evaluate the performance of the established model and its potential applications in the South China Sea by comparing the forecasts with the sounding and reanalysis data.

Section 2 will introduce the model, duct diagnosis process, and validation data. Section 3 will discuss the design and model configuration of this experiment. Next, Section 4 will present a specific evaluation of the simulation results. Finally, Section 5 will summarize this work and briefly discuss the uncertainties in the simulation results.

2. Model and Data Description

2.1. The COAWST Model

Warner et al. [19] developed the COAWST model, which coupled ocean, atmosphere, surface waves, and sediment transport components through the Model Coupling Toolkit (MCT). The ocean sub-model is the Regional Ocean Modeling System (ROMS), the atmosphere sub-model is the Advanced Research Weather Research and Forecasting (WRF) model, and the wave model is the Simulating Waves Nearshore (SWAN). The sediment model is the Community Sediment Transport Modeling System (CSTMS), which exchanges air–sea surface variables regularly through the MCT.

2.2. Diagnosis Scheme of Lower Atmospheric Ducts

Zhu et al. [8] classified the lower atmospheric ducts into four categories based on their profile characteristics (Figure 1): surface duct without a base layer (surface duct close to the sea/land surface), surface-based duct with a base height (surface-based duct), elevated duct, and complex duct (where both surface and elevated ducts exist simultaneously or have more than two layers of an elevated duct). Among them, h or h_t represents the top height of the duct layer, h_b represents the bottom height of the duct layer, ΔM represents the duct strength, and ΔT represents the thickness of the duct layer.

Figure 2 depicts the flowchart of the lower atmospheric duct diagnosis based on the outputs of the COAWST model, which primarily includes the following steps:

• The temperature, humidity, and other related variables were read from the model, and the modified atmospheric refractivity was calculated according to Formulas (1), (2), and (3). The height h₀ = 0 was initialized, and whether a duct occurred starting from the lowest layer was checked;
• Starting from height h₀, the first occurrence of a maximum peak of the modified refractivity was searched upwards, and the current height h₁ (shown in blue in Figure 1) and the current modified refractivity M₁ were then recorded. If there is no h₁, it is determined that there is no duct, and the diagnosis is completed;
• If h₁ exists, the search is continued upwards to determine the first minimum value of the modified refractivity. Next, the current height h₂ and modified refractivity M₂ were recorded;
• If there is no increased modified refractivity within the height of h₀~h₁, then the current duct is determined as a surface duct. At that time, the bottom height of the duct is 0 m. The top height of the duct is h₂. The duct strength is M₁–M₂. The duct thickness is h₂–h₀;
• If there is an increased modified refractivity below the height of h₁, the first point where the modified refractivity is less than M₂ starting from h₁ downwards is searched, and the height h₃ at that point is recorded;
• If no points of less than M₂ are found, and h₀ is still the initial height of 0 m, the current duct is then identified as a surface-based duct. At that time, the bottom height of the duct is 0 m. The top height of the duct is h₂. The duct strength is M₁–M₂. The duct thickness is h₂–h₀;
• If a point of less than M₂ is found, the height h₃ is recorded there, and it is determined that the current duct is an elevated duct. The bottom height of the duct is h₃. The top height of the duct is h₂. The duct strength is M₁–M₂, and the duct thickness is h₂–h₃;
• After the diagnosis process of the duct in this layer is completed, the initial height h₀ is replaced with the height h₂. Step b is now repeated, and the search is continued upwards to diagnose ducts in the next layer until the maximum height is reached.

Due to space limitations, this paper will only evaluate the characteristics of surface and elevated ducts to facilitate simulation results. In the future, the authors will investigate the four sub-types of ducts in detail. The following analysis includes surface ducts with a base height and those without a base layer. Elevated ducts remain unchanged. Complex ducts are divided and included in the statistics of surface and elevated ducts according to their duct types at different levels.

2.3. Observation Data for Evaluation

This study used the sounding data from the radiosonde stations around the South China Sea. The sounding data were derived from the data website of the University of Wyoming (http://weather.uwyo.edu/, accessed on 5 January 2024). This dataset was compiled from daily measurements performed twice (at 0000 h UTC and 1200 h UTC), primarily including air pressure, height, temperature, dew point temperature, relative humidity, and water vapor mixing ratio. The locations of sounding stations around the South China Sea over the study domain are marked in Figure 3, with the yellow dots representing the main stations in the following evaluations.

During the evaluation process, the sounding data were used as the references. The data from the European Centre for Medium-Range Weather Forecasts Reanalysis v5 (ERA5) (https://www.ecmwf.int/en/forecasts/dataset/ecmwf-reanalysis-v5, accessed on 5 January 2024) were also used as comparisons. The ERA5 data are produced by a General Circulation Model (GCM) for global climate and weather research. Its horizontal resolution is 0.25° × 0.25°, with a vertical resolution of 37 pressure layers ranging from 1000 hPa to 1 hPa and a time resolution of one hour. Many studies have proved the ERA5 reanalysis data to be one of the most reliable regional datasets with high accuracy in Southern China [20,21].

3. Experimental Design

This study selected the South China Sea and its surrounding areas as the research domain (Figure 3). The WRF sub-model was used for the entire domain. In contrast, the area framed by a red box is the domain used by the ROMS and SWAN sub-models, which are slightly smaller than the WRF domain. The domain design for the three sub-models aimed to maintain computational stability during the grid interaction process. The water depth around the edge of the South China Sea is about 200 m. The average water depth in the central area can reach about 4000 m. The surrounding terrain primarily comprises mountains and hills with a height of about 4000 m or less.

Table 1. Model configurations for the WRF, ROMS, and SWAN sub-models.

	WRF	ROMS	SWAN
Number of grid points	180 × 260	170 × 240	170 × 240
Horizontal resolution	12 km × 12 km	12 km × 12 km	12 km × 12 km
Vertical levels	45	16	1
Driving data	GFS	RTOFS	GFS_wave
Time step	30 s	30 s	600 s

presents the model configurations for WRF, ROMS, and SWAN. Here, the number of grid points in WRF is 180 × 260 with a horizontal resolution of 12 km × 12 km and 45 vertical layers. The eta values of the lowest 10 vertical levels in WRF are 1, 0.9998, 0.9995, 0.9992, 0.9985, 0.9981, 0.9979, 0.9975, 0.997, and 0.9965. The number of grid points in ROMS and SWAN is 170 × 240, with the same horizontal resolution of WRF and 16 vertical layers in ROMS.

The parameterization schemes used by WRF included the WSM6 microphysics scheme [22], the RRTM longwave radiation scheme [23], the Dudhia shortwave radiation scheme [24], the MM5 surface layer scheme [25], the Noah Land Surface Model [26], the YSU planetary boundary layer scheme [27], and the GF cumulus scheme [28]. For the configuration of ROMS, the Mellor–Yamada scheme was used to calculate the turbulent mixing in the vertical direction [29]. The Flather boundary condition method was applied to calculate anisotropic ocean currents to allow for the free propagation of wind-generated currents and tides [30].

The forecasting simulations started at 00:00 UTC on 29 July 2021 and ended at 00:00 UTC on 1 September 2021. The forecasts were conducted every day for the next 72 h. The forcing data of WRF were obtained from the Global Forecast System (GFS, available at https://www.ftp.ncep.noaa.gov/data/nccf/com/gfs/prod/, accessed on 5 January 2024). The forcing data of ROMS came from the Global Real-Time Ocean Forecast System (RTOFS, available at https://www.nco.ncep.noaa.gov/pmb/products/rtofs/, accessed on 5 January 2024), and the forcing data of SWAN were generated from the GFS_wave dataset. For the initialization of the models, WRF and ROMS used the default way for a cold start, whereas SWAN used the restart files from three days ago for a hot start.

The forecasting simulations lasted for one month, and they were divided into three groups of tests with different forecasting times ahead. Those are the forecasts for the next 24, 48, and 72 h. First, the sounding observations from all the radiosonde stations in the domain were collected to evaluate the bias of forecasted modified refractivity. Next, based on the sequence of bias in mean values and the situation of missing data in each station, observations from three stations were selected to represent the high, medium, and low levels of error for the following duct validation. In such a way, the researchers attempted to cover the different error levels to enhance the reliability of the findings.

4. Results

4.1. Evaluation of the Mean Bias at All Radiosonde Stations

Figure 4 depicts the mean RMSEs of the 24 h forecasted modified atmospheric refractivity at all observation stations below the height of 10,000 m over the study domain. The RMSE can be calculated by [31]

$$RMSE=\sqrt{\frac{1}{n}{\displaystyle {\sum }_{i=1}^n{\left(x_{mi}-x_{oi}\right)}^2}},$$

(4)

where x_mi and x_oi are the simulated and observed modified atmospheric refractivity in the time i, and n is the number of samples. Both observations and simulations were interpolated to the unified height using linear interpolation. The figure shows that the mean RMSE at each station in August was between 7 M and 11.5 M. This error level is basically consistent with the simulated error of Haack et al. [12], indicating the reliability of the predictions. Considering the station locations and the amount of valid data there, the sounding data at Haikou, Xisha, and Tanay stations were selected for the subsequent duct evaluations (the locations of the three stations are marked in yellow dots in Figure 3). Each of the three stations has data for 62 time points at 00:00 UTC and 12:00 UTC per day during the forecasting simulations. The mean RMSEs of the modified refractivity at these three stations were 9.6 M, 8.5 M, and 7.1 M, respectively. In Figure 4, they were ranked 3rd, 8th, and 13th, respectively, which could cover different error levels.

4.2. Evaluation of the Forecasted Duct Characteristics at Typical Stations

Figure 5 depicts the vertical profiles of RMSEs and the Pearson Correlation Coefficients (CCs) for the forecasted modified atmospheric refractivity at three stations. The vertical coordinate’s minimum value corresponded to each station’s altitude.

The CC can also be calculated by [32]

$$COR=\frac{{\displaystyle {\sum }_{i=1}^n(x_{mi}-\overline{x_m})(x_{oi}-\overline{x_o})}}{\sqrt{{{\displaystyle {\sum }_{i=1}^n\left(x_{mi}-\overline{x_m}\right)}}^2\cdot {\displaystyle {\sum }_{i=1}^n{\left(x_{oi}-\overline{x_o}\right)}^2}}},$$

(5)

Figure 5A–C shows that among the three groups of prediction tests, the mean RMSEs of the forecasts for 24 h ahead were the lowest, whereas the RMSEs of the forecasts for 72 h ahead were the highest. However, the errors of the 72 h forecasts were only 0.31 M~2.21 M larger than those of the 24 h forecasts, indicating that with the prolonging of the forecasting time, the increased prediction error could be maintained at a lower level within the forecasting period of 72 h. In addition, except for the differences in values, the profile patterns of the three groups of results were basically the same, suggesting that the vertical distributions of error did not change considerably with the prolonging of forecasting time.

As for the profile shapes of RMSEs, the simulation error usually reached peak values below 2000 m. With the further increase in height, the error gradually decreased, indicating that the simulations of COAWST within the boundary layer still needed further improvements. This phenomenon also appeared in the CC profiles (Figure 5D–F). A low correlation was detected below the heights of 2000 m, but the CC profiles did not exhibit an apparent changing trend above 2000 m. On average, the CCs for the 24 h forecasts at the three stations were 0.5, 0.68, and 0.4, respectively, whereas the CCs of the 72 h forecast results dropped to 0.28, 0.45, and 0.22, respectively.

Compared with ERA5 reanalysis data, the prediction results had lower RMSEs and correlation coefficients in most cases. In order to further investigate the statistical differences between ERA5 data and forecast results, Figure 6 also provides vertical profiles of mean bias, covariance, and standard deviation.

The standard deviation can be calculated by [33]

$$STD=\sqrt{\frac{1}{n}{\displaystyle {\sum }_{i=1}^n{\left(x_i-\overline{x}\right)}^2}}.$$

(6)

Figure 6A–C shows the vertical profiles of mean bias in modified atmospheric refractivity during the simulation period. The mean bias of forecasts was generally lower than that of ERA5 data, and the bias differences between them were larger than the RMSE differences in Figure 5A–C because some forecasting biases canceled each other out, while ERA5 data showed a constant negative bias.

According to Formula (5), the correlation coefficient is essentially the standardized covariance. Therefore, Figure 6D–I actually represent the numerator and denominator of Formula (5). It could be seen that both the covariances and standard deviations of the forecasting results were higher than those of ERA5 data, and the standard deviations of the forecasting results were basically more consistent with the observations. When calculating the correlation coefficient, the increase in the denominator (i.e., covariance) was obviously lower than the increase in the numerator (product of simulated and observed standard deviations). This difference ultimately led to the lower correlation coefficients of the forecasting results compared to ERA5 data.

Compared with ERA5 data, the RMSEs of COAWST prediction results were generally lower at the three stations. Especially at Xisha and Tanay stations, the RMSEs of 72 h forecasts with the highest error level were still lower by 2.54 M and 4.84 M compared to the RMSEs of the ERA5 data. Interestingly, there were also high-value zones in the RMSE profiles of the ERA5 data; however, the height zones of high values were higher than those of prediction results. The CCs of the ERA5 data were generally higher than the prediction results, which were 0.24, 0.11, and 0.19 higher than the 24 h forecasts for the Haikou, Xisha, and Tanay stations, respectively. It indicated that the reanalysis data are closer to the observation than the model prediction in terms of time variability. In addition, the vertical profiles of CC of ERA5 data were basically consistent with the prediction results but with a smaller variation range.

The characteristic parameters of predicted ducts were also validated. Figure 7 displays the histograms of occurrence rate, mean height, thickness, and strength for the predicted surface ducts by COAWST at the three stations. The surface ducts here included the surface-based duct with a base layer and the surface duct without base layers mentioned in Section 2.2. They were compared with the duct parameters diagnosed from sounding and the ERA5 data. Figure 7A–C shows the mean occurrence rates of surface ducts at the three stations during the entire simulation period. As the test with the highest accuracy, there were still some differences between the 24 h forecasts and the observations in the occurrence rates of the surface duct. The occurrence rate differences between them were –32.8%, 41.9%, and –6.5% for Haikou, Xisha, and Tanay stations, respectively. The differences between the three stations might be partly attributed to the size limits of sounding data. During the validation period of one month, there were only 62 data groups at each station without missing values. The data size might introduce some randomness in the statistics. For the mean height (Figure 7D–F), thickness (Figure 7G–I), and strength (Figure 7J–L) of surface ducts, the predicted ducts were lower in height, thinner in thickness, and weaker in strength. This result was consistent with the general error of the numerical models mentioned by Haack et al. [12]. In this study, the temporal output frequency of the WRF sub-model in COAWST was once per hour. The profiles of modified atmospheric refractivity were hourly-mean profiles, and the smoother profiles of modified atmospheric refractivity partly contributed to these errors. For the sounding data, some random events often occurred with the rise in sounding balloons, resulting in oscillations in the temperature and humidity measurement, which were likely to be diagnosed as ducts. In addition, probably due to the insufficient vertical resolution of the air temperature and pressure in the ERA5 data to capture the surface ducts, the occurrence of surface ducts was not diagnosed at the three stations using the ERA5 data and could not provide references on duct parameters.

Similar to Figure 7, Figure 8 also shows the occurrence rate, mean height, thickness, and strength of the predicted elevated ducts. Unlike the surface ducts, which highly depend on the thermal variations of the underlying surface, the elevated ducts originate entirely from the abnormal distributions of temperature and humidity stratification and, therefore, have more substantial randomness. The occurrence rates of elevated ducts were far lower than the diagnostic results based on observations. On average, the duct occurrence rates from the 24 h forecasts with the highest accuracy (Figure 8A–C) were lower by 68.8%, 45.2%, and 39.4% compared with the soundings for the Haikou, Xisha, and Tanay stations, respectively. These results are not as accurate as the surface ducts in Figure 7A–C. In addition to the more extensive altitude range for diagnosis and the weaker correlation between atmospheric stratification and the underlying surface, the prediction error of modified atmospheric refractivity gradually increased within the altitude below 2000 m (Figure 5A–C), which might also be one of the sources of diagnostic error of ducts. Compared with the sounding data, the predicted mean height (Figure 8D–F), thickness (Figure 8G–I), and strength (Figure 8J–L) of the elevated ducts also presented the characteristics of low height, thin thickness, and weak strength as the surface ducts. Furthermore, the characteristic parameters of the elevated ducts are more random. For example, at Xisha Station, the thickness and strength of the predicted elevated ducts were generally higher than the sounding data. In addition, with the prolonging of the forecasting time, the prediction accuracy of elevated ducts declined, and the declined degrees were more noticeable than those of the predicted surface ducts. This outcome is probably related to the quality of the driving field data and developing levels of current numerical models. The essence of an elevated duct is an extreme weather phenomenon such as temperature inversion and wind shear. With the prolonging of forecasting time, the increasing rates of prediction error for extreme weather features are much higher than the changing rates of error for general atmospheric variables. This might be why the prediction error of elevated ducts changed rapidly within 72 h in this study. In addition, a few elevated ducts were diagnosed using the ERA5 data at the Xisha Station (but not found at other stations). The mean height and thickness were basically similar to the prediction results but with weaker duct strength.

The sounding data would likely miss some ducts because of its severe limitations in the vertical resolution of the data and measurement error. Referring to previous studies,

Table 2. Hit rate, FAR and CSI of ducts from the forecasts and ERA5 data at Haikou, Xisha, and Tanay stations.

		Surface Duct			Elevated Duct
		Haikou	Xisha	Tanay	Haikou	Xisha	Tanay
Hit rate	24 h	26.32%	100%	20%	2.33%	31.11%	7.14%
	48 h	18.42%	96.67%	10%	4.65%	24.44%	3.57%
	72 h	18.42%	93.33%	0%	2.33%	22.22%	0%
	ERA5	0%	0%	0%	0%	4.44%	0%
FAR	24 h	44.44%	46.43%	66.67%	0%	17.65%	50%
	48 h	46.15%	50.85%	66.67%	0%	8.33%	50%
	72 h	58.82%	52.54%	100%	0%	23.08%	0%
	ERA5	0%	0%	0%	0%	0%	0%
CSI	24 h	0.22	0.54	0.14	0.02	0.29	0.07
	48 h	0.16	0.48	0.08	0.05	0.24	0.03
	72 h	0.15	0.46	0	0.02	0.21	0
	ERA5	0	0	0	0	0.04	0

lists the hit rates, false alarm ratio (FAR), and critical success index (CSI) of predicted ducts when the sounding data for observed ducts were analyzed [34]. Table 2 also compares the diagnostic results based on the ERA5 data. The hit rates of predicted surface ducts from 24 h forecasts at three stations were 26.32%, 100%, and 20%, respectively. The hit rates of predicted elevated ducts were 2.33%, 31.11%, and 7.14%, respectively. Haack et al. [12] demonstrated that the highest hit rate for the lower atmospheric ducts from the COAMPS was 43%, and the lowest rate from MM5 was 6%. The hit rate of the prediction model in this study was roughly equivalent to levels of MetUM, 25%. As a comparison, there were still a few hit times while using the ERA5 data for duct diagnosis.

Overall, the performance of forecasts at the Xisha station was better than those at the Haikou and Tanay stations in the coastal area. As island environments are less affected by the interaction between sea and land processes, they are closer to the marine environment. This indicated that the duct forecasts based on the ocean–atmosphere coupling model established in this study performed better in the marine environment. To enhance the credibility of the results, we also selected two island stations, Bach Longvi and Ranai (see Figure 3 for locations), and calculated the hit rate, FAR, and CSI there. For the 24 h surface duct forecasts at the Bach Longvi station, these three indices were 75.0%, 74.43%, and 0.26, respectively, while the indices for elevated duct forecasts were 28.57%, 33.33%, and 0.25, respectively. The indices at the Ranai station were not significantly different from those at the Bach Longvi station, with 57.14%, 60.0%, and 0.31 for surface duct forecasts, and 75%, 62.5%, and 0.33 for elevated duct forecasts, respectively. The evaluation indices at both island sites were better than those at the coastal area Haikou and Tanay stations.

4.3. Spatial Distributions of the Forecasted Duct Characteristics over the Domain

The evaluations above at the station demonstrated that the duct prediction model established in this study had high reliability. The following analysis will focus on the spatial distributions of the characteristics of the ducts and the differences among the forecasts with different forecasting times. Figure 9 depicts the spatial distributions of the occurrence rate of surface and elevated ducts from the forecasts 24 h ahead, accompanied by the differences between 48 h forecasts minus 24 h forecasts and between 72 h forecasts minus 24 h forecasts. Generally, no matter whether they are the surface or elevated ducts, the occurrence rates of ducts at sea were higher than that on land, which was consistent with previous studies [3,35].

Figure 9A shows the spatial distributions of the occurrence rates of surface ducts from the 24 h forecasts. The occurrence rate of surface ducts on land was less than 70%, and the occurrence rate of surface ducts at sea was more than 80%. In the whole domain of the South China Sea, the mean occurrence rate of surface ducts was 68.78%, with a maximum rate of about 99%. Except for the Taiwan Strait, there was a high occurrence rate over the entire sea area with no apparent spatial differences. Figure 9B,C shows the differences in the occurrence of surface ducts from the 48 h and 72 h forecasts minus the 24 h forecasts. The figures demonstrate that the distributions of the occurrence difference presented polarized distributions. In the northern coastal areas of the domain, the occurrence rate of surface ducts from forecasts of 24 h ahead was low, and the 48 h and 72 h forecasts further reduced the occurrence of surface ducts in this area. In contrast, in the eastern part of the South China Sea, the 48 h and 72 h forecasts further increased the occurrence rate of surface ducts there, with a maximum increase of about 15%.

Figure 9D demonstrates the spatial distributions of the occurrence rate of the elevated ducts from the 24 h forecasts, with a mean occurrence rate of 11.78%. The land–sea difference here was also rather noticeable. The occurrence rates on lands were generally less than 15% and gradually decreased from coastal areas to inland. The occurrence rates of offshore ducts were generally more than 15%, and there were two high-value centers, including the seas between Hainan and Xisha Islands and the seas between Nansha and Kalimantan Islands, with a maximum rate of more than 30%. On average, the mean occurrence rate of elevated ducts over the domain was 11.78%, with the highest rate of 39.7% and the lowest rate of 0%. Figure 9E,F depicts the different distributions of elevated ducts from the 48 h and 72 h forecasts minus the 24 h forecasts. Unlike surface ducts, the differences were not apparent on land but in a polarized distribution pattern at sea. In the area with high occurrence between Hainan and Xisha Islands, the occurrence rates of 48 h and 72 h forecasts significantly decreased by 5~15%. However, in another area with high duct occurrence between the Nansha and Kalimantan Islands, the occurrence rate of elevated ducts rose in an extensive range, with an increase of more than 15% in some areas.

Based on the differences between the 48 h, 72 h, and 24 h forecasts in Figure 9, the Pearl River Estuary (21°~23.8° N, 112°~116° E), the Xisha Islands (14°~17° N, 111°~114° E), and the eastern Nansha Islands (8°~11.5° N, 115.8°~118° E) with noticeable changes were selected for further analysis to explore the reasons for the differences among the three groups of forecasting tests. The three selected areas were framed by black lines in Figure 9. Figure 10 demonstrates the vertical profiles of the occurrence time counts of ducts from the 24 h, 48 h, and 72 h forecasts in the three selected areas during one month. To show the differences between surface and elevated ducts, the vertical axis was divided into two sections: 0 m~100 m and 100 m~2500 m. Overall, within the 100 m height near the sea/land surface (Figure 10A–C), the number of ducts below the height of 20 m reached the most, most of which were surface ducts near the sea surface. With the increase in height (Figure 10D,E), the number of duct occurrences reached another peak near the height of 500 m~1000 m, which was usually an elevated duct. Considering the distributions in Figure 9, the numbers of surface and elevated ducts from the 24 h forecasts in the Pearl River Estuary area (Figure 10A,D) were the most; therefore, the occurrence rates of ducts from the 24 h forecasts in this area were higher than those from the other two groups of forecasting tests. In addition, the height where the maximum number of occurrences from the 24 h forecasts was lower than that from the other two groups of tests, resulting in a relatively lower height of ducts. The duct occurrence times from the 24 h forecasts in the Xisha Islands (Figure 10B,E) were slightly higher than that in the other tests, and there were no significant differences in the peak height among the three forecasting tests. Furthermore, the duct occurrence times from the 24 h forecasts in the eastern Nansha Islands (Figure 10C,F) were relatively less, leading to a lower occurrence rate than that from the other tests.

Given the close relationship between the conditions of the underlying surface and atmospheric ducts, hydrometeorological variables at the ocean/land surface, including the forecasted sea surface temperature (SST), 2 m air temperature (T2m), 2 m relative humidity (RH2m), sea level pressure (SLP), and 10 m wind speed, were collected to calculate their correlations with the occurrence rates of atmospheric ducts.

Table 3. Spatial correlation coefficients between surface variables and duct occurrence rates over the domain.

Surface Variable	Surface Duct			Elevated Duct
	24 h Forecasts	48 h Forecasts	72 h Forecasts	24 h Forecasts	48 h Forecasts	72 h Forecasts
SST	0.86	0.84	0.84	0.68	0.63	0.62
T2m	0.56	0.53	0.52	0.41	0.33	0.30
RH2m	−0.15	−0.16	−0.14	−0.04	0.04	0.07
SLP	0.58	0.55	0.55	0.52	0.48	0.48
Wind speed	0.41	0.44	0.43	0.47	0.36	0.27

presents the spatial correlation coefficients between each variable and the duct occurrence rates during the simulation period. As shown in the table, SST presented a higher spatial correlation with the duct occurrence rates compared to other variables. The spatial correlation coefficient was denoted here as the correlation coefficient of the mean spatial distributions between two variables. The mean spatial distribution of each variable (two-dimensional array) was reshaped into a one-dimensional array, and then, the spatial correlation coefficient was calculated using Formula (5). Specifically, the spatial correlation between surface duct occurrence rates and SST exceeded 0.8. Although elevated ducts did not directly interact with the underlying surface, their correlation with SST could also reach 0.6 or above. This high correlation indicated a strong relationship between the thermal conditions of the underlying surface and the duct occurrence. In regions with higher SST, the ducts occurred more often [36]. On the other hand, variables such as air temperature, humidity, and atmospheric pressure played important roles in the formation of ducts but exhibited spatial correlation coefficients mostly below 0.6, lacking the strong correlations observed in SST.

Furthermore, the spatial distributions of duct characteristic parameters were also analyzed in this study. We attempted to investigate the changes in the forecasted duct characteristic parameters by comparing the differences in results for different forecasting times in order to confirm whether the forecasting model established in this study was suitable for conducting longer-term forecasts. Figure 11 demonstrates the spatial distributions of the mean height, thickness, and strength of the surface ducts in 24 h forecasts and the differences in the 48 h and 72 h forecasts minus the 24 h forecasts. Generally, although the occurrence rates of surface ducts on land were far lower than at sea, the mean height, thickness, and strength of surface ducts on land were higher, thicker, and stronger than those at sea. It probably indicated that the thermal and dynamical conditions in the land surface seemed unable to maintain a weaker duct phenomenon.

Figure 11A presents the distributions of mean height of surface ducts from the 24 h forecasts. Within the domain, the duct height ranged from 0 m to 1200 m. The distributions of duct height at sea were relatively uniform at the levels of 100 m below. The duct height on lands ranged from 200 m to 1200 m. The topographic elevation of the terrain affected these land–sea differences in the duct height. If the elevation factor on land was excluded, the height distributions of surface ducts still demonstrated the distributions of low at sea and high on land. The heights of surface ducts on land were about 5 m~10 m higher than those at sea. The mean height differences between the 48 h, 72 h, and 24 h forecasts (Figure 11B,C) ranged from –3 m to 3 m. The main differences were concentrated in the coastal areas of the Chinese mainland, where the duct heights from the 48 h and 72 h forecasts were larger.

The spatial distributions of the mean thickness (Figure 11D) and strength (Figure 11G) of the surface ducts from the 24 h forecasts were similar, indicating a higher consistency of the two parameters. The mean thickness and strength of the surface ducts in the northeast areas of the South China Sea were slightly higher than those in other sea areas. Compared with the 24 h forecasts, the difference distributions of duct thickness and strength were similar to the mean height differences (Figure 11E,F,H,I). With the prolonging of the forecasting time, the characteristics of the surface ducts were obviously enhanced in the coastal and northeast areas of the South China Sea. In the other areas, however, the changes in duct parameters were relatively weak.

Similar to the surface ducts in Figure 11, Figure 12 shows the spatial distributions of the mean height, thickness, and strength of the elevated ducts from the 24 h forecasts and the differences among the 48 h, 72 h, and 24 h forecasts. Generally, the duct heights at sea were about 600 m~1100 m, with a mean thickness of more than 100 m and a strength of 3 M~7 M. The heights of the elevated ducts at sea were higher than those on land, and the mean thickness and strength were also significantly stronger at sea. The areas with high values coincided with the areas with high occurrence rates in the northwest of the South China Sea and the coastal area of Kalimantan Island (Figure 9D). This outcome indicated that elevated ducts were found to easily occur in those regions and had higher strength. Ships passing by these areas needed to pay more attention to the interference of elevated ducts on radio equipment within a height of hundreds of meters. In most areas, the parameters of the elevated ducts from the 48 h and 72 h forecasts are enhanced. Only in a few areas, such as the areas south of Hainan Island and east of the Nansha Islands, were the thickness and strength of ducts weakened.

5. Discussion

Although this study has conducted a preliminary evaluation of the established duct prediction model, there are still some uncertainties in model physics, grid resolution, experimental design, and data selection. These limitations may increase the errors in the simulations. In the future, the model should be improved to increase the representativeness and credibility of the simulation results, including the following aspects.

• In terms of model improvements, the physical parameterization schemes used in this study were obtained from sensitivity tests conducted in the previous studies for the China Seas. When applied in the South China Sea, the error of this scheme configuration may not be minimized. Therefore, a local optimization test specific to the South China Sea is needed to develop a scheme configuration with more minor errors. In addition, the assimilation process was not included in the forecasting simulation at this time. In the next step, the WRF data assimilation module will assimilate observation and satellite data to improve the duct forecasting accuracy further;
• In the simulation process, the model resolution also plays a crucial role in simulation accuracy. Previous studies have shown that high-resolution numerical models are capable of simulating finer-scale processes for structures or phenomena generated within a few hours, even without assimilation [37]. In this study, the horizontal resolution of the model was set to 12 km without using grid nesting to enhance the resolution. Such a resolution setting was sufficient for simulating a homogeneous marine environment but fell short in describing finer-scale processes such as coastal sea-land breezes. To validate this point, an additional nested forecasting test was conducted. The test ran from 00:00 UTC on 15 August 2021 to 00:00 UTC on 20 August 2021. The external grid remained at a 12 km resolution, while the internal grid had a resolution of 4 km. By comparing the mean bias of the modified atmospheric refractivity between the nested and original simulations, it was found that the higher-resolution nested simulation resulted in reduced errors by 12.2% and 3.2% at the Haikou and Tanay stations in coastal areas. However, for the Xisha Island station, which is closer to the marine environment, the nested test only reduced the error by 1.3%. This variation among stations suggested that marine weather forecasting was not highly sensitive to changes in model resolution, but coastal areas affected by finer-scale processes such as sea-land breezes were more sensitive to resolution changes. To enhance credibility, the analysis was further conducted at two island stations (Bach Longvi and Ranai) and two coastal stations (Kowloon and Kota Kinabalu) (see Figure 3 for station locations). The results demonstrated that at a resolution of 4 km, the modified atmospheric refractivity bias at the two island stations changed by only 1.1% and −2.1% compared to the original test at 12 km resolution, while the two coastal stations changed by 10.4% and −6.9%. It is important to note that this study primarily focuses on the South China Sea region, and a resolution of 12 km is sufficient for the analysis. However, future research will involve conducting more high-resolution simulations to analyze the effects of small-scale processes in coastal areas on the modified atmospheric refractivity and atmospheric ducting processes.
• In addition to horizontal resolution, the vertical resolution of the model also plays a significant role in the atmospheric ducting diagnosis process. To investigate the influence of model layering on duct diagnosis, three additional simulations were conducted with different vertical layering. One test reduced the model layers from 45 layers in the original test to 37 layers, and another test further reduced the vertical layers to 29 layers, which corresponds to the layers below 50 hPa in ERA5 data. The last one increased the layers to 53 layers. In the test with 37 layers, the eta values of the lowest 10 vertical levels in WRF are 1, 0.999, 0.998, 0.997, 0.995, 0.994, 0.993, 0.992, 0.990, and 0.987. In the test with 29 layers, the eta values of the lowest 10 vertical levels in WRF are 1, 0.998, 0.996, 0.993, 0.990, 0.974, 0.946, 0.905, 0.880, and 0.835. In the test with 53 layers, the eta values of the lowest 10 vertical levels in WRF are 1, 0.9998, 0.9995, 0.9992, 0.9985, 0.9981, 0.9979, 0.9975, 0.997, and 0.9965. The three tests ran from 00:00 UTC on 15 August 2021 to 00:00 UTC on 21 August 2021 for a total of 6 days. The numbers of ducts diagnosed from the tests and the original test are shown in

  Table 4. Numbers of ducts diagnosed from the tests with different vertical layers and ERA5 data at Haikou, Xisha, and Tanay stations from August 15 to 21, 2021.

  		Surface Duct			Elevated Duct
  		Haikou	Xisha	Tanay	Haikou	Xisha	Tanay
  Test with 53 layers		78	149	14	6	44	3
  Original test (45 layers)		76	144	12	2	40	2
  Test with 37 layers		32	93	0	7	0	0
  Test with 29 layers		6	0	0	2	0	0
  ERA5 (29 layers)		0	0	0	0	2	0

  . The results indicated that the number of ducts obtained from the numerical model was indeed related to the vertical layering of the model. The test with 29 layers diagnosed the lowest number of ducting cases, similar to the ERA5 data. The test with 37 layers diagnosed more ducting cases, while the test with 53 layers did not show significant differences from the original test. This suggested that a smaller number of vertical layers was not conducive to ducting diagnosis, and as the number of model layers increased, the vertical layering configuration was no longer a barrier to ducting diagnosis. In addition to the impact of changes to vertical layer density, it is worth noting that there may be some additional non-negligible impact from changes to the vertical location of those layers. It is still an outstanding unknown if more added layers are concentrated in the lowest few kilometers of the atmospheric column, and more detailed investigations should be made in the future.
• In terms of experimental design, the forecasting test only lasted for one month, which was too short to determine the long-term application effects of this duct prediction model. During the forecasting period, there were a few intense cyclones in the South China Sea. Therefore, the changes in ducts’ characteristics under extreme weather conditions were not discussed. In the future, it is necessary to conduct longer-term evaluations covering various types of weather to enhance the representativeness of the forecasting results;
• The observation data for evaluations in this study also need improvements in representativeness. Due to the inability of ERA5 reanalysis data to serve as an evaluation reference, this study only selected the sounding data from three stations to verify the ducts’ characteristics. The selection way of stations still had some randomness. Due to the scarcity of publicly available data, it is difficult for us to collect atmospheric-sounding data over the ocean. Therefore, this study could only be verified using publicly available sounding data from coastal areas or islands. Radiosonde stations located on several islands are less influenced by the interaction between sea and land processes, and their data are believed to be the closest to the marine environment. It was a pity that the sounding data were only available at two specific moments per day, and even so, there were some missing values in the sounding data. Therefore, more and denser marine observation data need to be collected in the future for more detailed evaluations;
• Due to space limitations, this study only considered the overall characteristics of surface and elevated ducts. It did not provide a more detailed classification and evaluation for each duct type, nor did it involve an analysis of the formation mechanism of ducts. The authors will conduct more in-depth research in the future.

6. Conclusions

In conclusion, this study established a regional prediction model for lower atmospheric ducts based on the COAWST model, which coupled the atmospheric, ocean, and wave processes. Using the GFS, RTOFS, and GFS_wave global forecasting fields as forcing, one prediction simulation test was conducted over the South China Sea for the next 72 h daily. The test lasted for the whole month of August 2021. According to the different forecasting times, the test was divided into three groups: 24 h, 48 h, and 72 h forecasts ahead. The prediction model’s performance was evaluated by comparing it with the sounding data and ERA5 reanalysis data, and the differences in duct characteristics among tests with different forecasting times were also analyzed. The evaluation results demonstrated that the prediction model established in this study had an excellent forecasting ability for lower atmospheric ducts and had a potential for operational applications. More details of the conclusion are presented below.

• During the simulation period in August 2021, the RMSEs of the predicted modified atmospheric refractivity were maintained between 7 M and 11.5 M at the 13 radiosonde stations around the South China Sea, with a high accuracy. The sounding data from three radiosonde stations with different error levels were selected for further evaluation. The results demonstrated that the predicted modified atmospheric refractivity could maintain a high accuracy in the forecasting time of 72 h, with minor changes in error. Compared with the observations, the maximum error of the predicted modified atmospheric refractivity appeared near the height of 2000 m, and the error gradually decreased with the increase in height. At the three stations, the RMSEs of the predicted modified atmospheric refractivity were generally lower than those of the ERA5 reanalysis data. However, the CCs of the forecasts with observations were slightly lower than the ERA5 data;
• The atmospheric duct characteristic parameters evaluations demonstrated that the duct forecasts 24 h ahead had the highest accuracy compared to the other two groups of tests. In the locations of Haikou, Xisha, and Tanay radiosonde stations, the predicted occurrence rates were generally lower than the diagnostic results from observations. The mean duct height, thickness, and strength were generally low, thin, and weak. This was because the predicted hourly modified refractivity gradient at each layer was smoother than the instantaneous observations with more randomness. The prediction accuracy was equivalent to the levels of duct forecasts using MetUM and other models in previous studies. The duct forecasting model established in this study, based on the ocean–atmosphere coupling model COAWST, was more suitable for marine duct forecasting. The forecasting accuracy at island radiosonde stations closer to the marine environment was higher than that in coastal areas. In addition, the prediction accuracy of surface ducts was generally higher than that of elevated ducts. This study also used the ERA5 reanalysis data for diagnosing ducts to enhance the reliability of the results. However, unfortunately, the ERA5 data could hardly be diagnosed with ducts and could not provide references for duct evaluation. This is due to the limited vertical resolution of ERA5 data below 2000 m, where atmospheric duct activity was more frequent. Subsequent sensitivity tests on vertical layers indicated that the diagnosed number of ducts increased with an increase in the number of vertical layers. Despite some uncertainties in these tests, the diagnosed duct number was proved not to continue to grow rapidly as the vertical layer number further increased. After surpassing 50 model layers, the growth of diagnosed duct numbers slowed down, and the number of vertical layers was not always the primary limiting factor for duct diagnosis;
• The spatial distributions of predicted surface and elevated ducts indicated that the duct occurrence rates at sea were significantly higher than on land. The occurrence rates of surface ducts were higher than those of elevated ducts in the same area. Throughout August, the occurrence rates of surface ducts accounted for about 85.5% of all duct events. The subsequent statistics over typical areas demonstrated that compared to the 24 h forecasts, fewer ducts occurred in the Pearl River Estuary and the Xisha Islands, with more ducts in the eastern Nansha Islands for the 48 h and 72 h forecasts. The most apparent differences in duct occurrence occurred at heights below 20 m and 500~1000 m. Furthermore, even after deducting terrain height, surface ducts occurring on lands were still higher, thicker, and stronger than at sea. The situation was reversed for elevated ducts, where the duct height, thickness, and strength at sea were generally higher than those on land. From the differences among the three groups of tests with different forecasting times, the differences in the duct height, thickness, and strength were consistent. Compared to the 24 h forecasts, the characteristics of surface ducts from the 48 h and 72 h forecasts were enhanced in the coastal areas of the northern South China Sea. The characteristics changes in elevated ducts were complicated. The duct parameters were enhanced in the northern and central regions of the South China Sea and weakened in the southern and western regions.