# Development and Evaluation of a Hydrometeorological Forecasting System Using the Coupled Ocean-Atmosphere-Wave-Sediment Transport (COAWST) Model

Authors: Jing Zou; Chesheng Zhan; Haiqing Song; Tong Hu; Zhijin Qiu; Bo Wang; Zhiqian Li
Journal: Advances in Meteorology
 (2021-1-04)

---

## Abstract

In this study, an experimental hydrometeorological forecasting system was developed based on the Coupled Ocean-Atmosphere-Wave-Sediment Transport (COAWST) model. The system downloads global real-time ocean, atmosphere, and wave forcing data, producing regional forecasts every day. A coastal area in South China, encompassing Hainan Island, Leizhou Peninsula, and surrounding sea areas, was chosen as the study domain. A series of 72-hour forecasting simulations were conducted in the area, lasting from July 27 to August 31, 2019. The forecasts throughout August were chosen for evaluation with station observations, along with two sets of reanalysis data, ERA5 and CLDAS. The evaluation results revealed that the COAWST model had high potential for routine forecasting operations. The 24 h forecasts, with a lead time of 24 hours, had high accuracy, while the 48 h and 72 h forecasts did not differ greatly in terms of performance. The distributions of bias between forecast and reanalysis data showed obvious differences between land and sea, with more forecasted precipitation and lower temperatures in land grids than in sea grids. In most cases, the forecasts were closer to ERA5 in terms of means and other statistical measures. The forecasts enlarged the land-sea differences of temperature when compared with ERA5 and strengthened summer monsoon with more moisture transported to land areas. Resulting from that, a forecasted bias of lower surface pressure, higher air humidity, stronger south wind, and so forth was also detected over the domain but at low values.


---

## Body

## 1. Introduction

The forecast skill of numerical weather prediction (NWP) has consistently increased over the past 40 years [1, 2]. Bauer et al. found that the forecast skill, which is defined as the correlation between forecasts and verified analyses at a height of 500 hPa, has been increasing by about one day per decade [3]. However, NWP still has several limitations in terms of representing physical processes, ensemble modeling, and model initialization. Further, improvements are required in observation and computing [4–6].To improve the representation of physical processes in hydrometeorological simulations, many studies have shown that the sea-surface temperature (SST) and wave processes at the air-sea interface influence air-sea surface roughness, momentum, and heat fluxes, improving the accuracy of NWP in coastal regions [7–9]. Traditional uncoupled atmosphere models usually adopt observations or modeling data as one-way inputs of SST and use idealized parameters to describe changes in air-sea surface variables (i.e., surface roughness) during the forecast process. In contrast to these uncoupled models, coupled forecasting models—which combine mesoscale atmospheric models with ocean and wave models—can consider the dynamic interactions between atmospheric boundary layers and marine upper layers. It has been proven that these coupled models improve air-sea interactions and facilitate more accurate predictions than uncoupled models, especially in extreme events, that is, tropical cyclones (TCs) [10–13].A number of fully coupled ocean-atmosphere-wave models have been developed to simulate air-sea interactions [14–16]. Several of these fully coupled models have been implemented for the simulation of TCs or for operational forecasts. For example, the 3-dimensional Coupled Ocean/Atmosphere Mesoscale Prediction System (COAMPS) was developed by the United States Naval Research Laboratory (NRL) in the mid-1990s [17]. The COAMPS model included a data assimilation system, wave models, and other improvements to previous models and was used for operational forecasting in the U.S. navy [18–20]. The University of Miami Fully Coupled Atmosphere-Wave-Ocean Modeling (UMCM) system was developed with two modes, including different atmosphere-wave-ocean model options [21]. UMCM-MWP, one mode of UMCM, included the fifth-generation Pennsylvania State University-National Center for Atmospheric Research Mesoscale Model (MM5), the third-generation wave model (WAVEWATCH III), and the three-dimensional Price-Weller-Pinkel ocean model (3DPWP). It used an improved wind-wave coupling scheme that was verified using Hurricane observations [22]. Liu et al. established a coupled atmosphere-wave-ocean modeling system (CAWOMS) that included a sea-spray model and studied the effects of atmosphere-wave-ocean coupling on TC intensity.Among these fully coupled ocean-atmosphere-wave models, the Coupled Ocean-Atmosphere-Wave-Sediment Transport (COAWST) model, developed by Warner et al., has a high potential for application in forecasting [23, 24]. The model is open-access on the website (https://code.usgs.gov/coawstmodel/COAWST) and is being updated continuously. The COAWST model has been widely applied in a number of fields, such as TC hindcasts [25, 26], air-sea interaction mechanisms during TCs or other extreme events [27–29], coastal hydrodynamic processes [30, 31], and sediment-dynamics simulations [32]. According to studies on TC simulation, the fully Coupled Ocean-Atmosphere-Wave-Sediment Transport (COAWST) model provided lower SSTs and better general performance regarding TC intensity in comparison to uncoupled atmosphere models.Owing to the good performance of the COAWST model, some agencies or research groups implemented this model in operational forecasting products. For example, the Ocean Observing and Modeling Group (OOMG) in North Carolina State University developed two forecasting systems, the Coupled Northwest Atlantic Prediction System (CNAPS) and the South Atlantic Bight-Gulf of Mexico (SABGOM) model, and provided the three-dimensional forecasting products [33]. The Madeira Oceanic Forecasting System, which was developed by the Oceanic Observatory of Madeira, Portugal, utilized the COAWST model with an offline river routing model for operational predictions of hydrometeorological environment and suspended particulate matter (SPM) concentration near the Madeira Island [34, 35]. In addition, the Woods Hole Coastal and Marine Science Center, U.S. Geological Survey, where the COAWST model was developed, and the Hydro and Agro Informatics Institute, Thailand, also released regional operational products based on the COAWST model.This paper attempts to establish an experimental forecasting system over coastal areas of China based on the COAWST model. One real-time forecast test, lasting for a month, was conducted in a coastal area of South China to evaluate the forecasting potential of the COAWST model during routine operation. The rest of this paper introduces the descriptions of the COAWST model, forecasting system flow, and experimental design and subsequently provides results, conclusions, and discussions.


## 2. Methodology

### 2.1. Model Description

The COAWST model includes one atmosphere model, the Weather Research and Forecasting (WRF) model; one ocean model, the Regional Ocean Modeling System (ROMS); one wave model, the Simulating Waves Nearshore (SWAN) model; and one sediment-transport model, the Community Sediment Transport Model (CSTM).The WRF model is a nonhydrostatic, quasi-compressible, mesoscale atmosphere model with a number of physical schemes. In the COAWST model, the dynamic core of the Advanced Research WRF (ARW) model is used [36]. The ROMS model is a free-surface, topography-following-coordinate model that solves the three-dimensional Reynolds-averaged Navier-Stokes equations using hydrostatic and Boussinesq approximations [37, 38]. The SWAN model is a spectral-wave model that simulates the generation and propagation of wind-induced waves in coastal waters by solving the spectral-density evolution equation [39]. The CSTM is a sediment-hydrodynamic model that simulates the erosion, deposition, and transport processes across different sizes and types of sediments [24]. In this study, the CSTM was not utilized.The COAWST model utilizes the model-coupling toolkit (MCT) to achieve communication among the submodels via the message passing interface (MPI) [40]. In addition, the Spherical Coordinate Remapping Interpolation Package (SCRIP) is used to support variable interpolations in different coordinate systems [41]. In the coupling process, WRF receives sea-surface roughness from SWAN and the SST from ROMS, while it provides wind at a 10 m level (W10 m) to SWAN and surface stress and net heat fluxes to ROMS. To facilitate this variable exchange between ROMS and SWAN, ROMS receives the surface- and bottom-wave direction, height, length, period, percentage breaking, energy dissipation, and bottom orbital velocity from SWAN, while it provides bathymetry, bottom elevation, sea-surface height, and depth-averaged currents to SWAN.


### 2.2. Forecasting System Flow Description

Based on the COAWST model, an experimental forecasting system was established in this study, the system flowchart of which is shown in Figure1. The system is controlled by executable shell scripts with timers and consists of one data server, one compute server, and the user client.Figure 1
Flowchart of the developed real-time forecast system.The data server downloads and preprocesses the global forecast data every day from the public data server of the National Centers for Environmental Prediction (NCEP). The downloaded data includes the Global Forecast System (GFS) data and the Global Data Assimilation System (GDAS) data for WRF (available online:https://www.nco.ncep.noaa.gov/pmb/products/gfs/), the Real-Time Ocean Forecast System (RTOFS) data for ROMS (available online: https://www.nco.ncep.noaa.gov/pmb/products/rtofs), and NCEP operational wave-product data from the global WAVEWATCH III model for SWAN (available online: https://www.nco.ncep.noaa.gov/pmb/products/wave/). The resolution information of the downloaded data is listed in Table 1.Table 1
Resolution information of the downloaded forecast data.
GFSRTOFSWaveSpatial resolution0.25° × 0.25°0.083° × 0.083°0.5° × 0.5°Time resolution3 hours3 hours3 hoursStart timeUTC 0 : 00UTC 0 : 00UTC 0 : 00Time rangeSubsequent 72 hoursSubsequent 72 hoursSubsequent 72 hoursIn the compute server, the GFS and GDAS data are processed by the WRF Preprocessing System (WPS) and 3-dimensional Variational data assimilation system (3DVar), respectively, to produce WRF’s initial and boundary conditions. The GFS data, RTOFS data, and NCEP operational wave-product data are processed in the ROMS and SWAN preprocessing module to produce the boundary files for SWAN, along with the initial and boundary files of for ROMS.In this study, the initial file for SWAN was provided by the previous restart file from 72 hours earlier—instead of the default idealized initialization—as many previous studies on the COAWST model have done [26, 42, 43]. The bathymetry dataset for ROMS used the DTU10 dataset, with a resolution of 0.125° × 0.125°, from the Technical University of Denmark (available online: https://www.space.dtu.dk/english/Research/Scientific_data_and_models/Global_Bathymetry_Model). The tide dataset for ROMS was derived from the TPXO7.0 dataset, with a resolution of 0.25° × 0.25°, which was obtained using the Oregon State University Tidal Prediction Software (OTPS). The TPXO global tide model provides eight primary (M2, S2, N2, K2, K1, O1, P1, and Q1) and two long (Mf and Mm) period harmonic constituents as amplitudes of the earth-relative sea-surface elevation [44, 45]. The TPXO7.0 dataset can be found online at http://www−po.coas.oregonstate.edu/∼poa/www−po/research/po/research/tide/global.html.The initial and boundary files for WRF, ROMS, and SWAN are transferred to the COAWST model for forecast simulations; subsequently, the simulation results are processed using a postprocessing module and are visualized for the end user.


### 2.3. Experimental Design

The configuration parameters for model gridding and operation are listed in Table2. The COAWST model in this study used two nested grids to obtain results with a higher resolution via dynamic downscaling. The model was compiled with fully coupling options and three submodels exchanged variables every 1800 seconds, controlled by the main program via the MCT.Table 2
Configuration parameters for each submodel.
WRFROMSSWANTime step30 s60 s180 sGrid nestingYesYesYesOuter grid number100 × 10090 × 9090×90Inner grid number100 × 10090 × 9090×90Horizontal grid resolution18 km for outer grids, 6 km for inner grids18 km for outer grids, 6 km for inner grids18 km for outer grids, 6 km for inner gridsVertical layer number3916NoneInitial data sourceGFS data assimilated by GDASGFS and RTOFS dataRestart file of simulation 3 days earlierVariable exchange frequency1800 s−11800 s−11800 s−1The default scheme options in the COAWST model were used without optimization in this study. The schemes used in WRF were the WRF Single-Moment 3-class (WSM3) microphysics scheme [46], the Rapid Radiative Transfer Model (RRTM) longwave scheme [47], the Dudhia shortwave scheme [48], the Noah Land Surface Model (LSM) [49], the modified Kain-Fritsch cumulus scheme [50], the Eta surface layer scheme [51], and the Mellor-Yamada-Janjic (MYJ) planetary boundary-layer scheme [52, 53]. For ROMS setup, the Mellor-Yamada scheme was used to compute vertical turbulent mixing [53], and the Flather boundary-condition method was used for barotropic currents to allow the free propagation of wind-generated currents and tides [54]. For the SWAN model, the Madsen scheme was used for wave bottom dissipation parameterization [55], and the Komen formulation was used to simulate wind-induced wave growth [56].The forecast simulation tests began at UTC 0 : 00 every day with a forecast lead time of 72 hours and lasted from July 27 to August 31, 2019. The first three tests (from July 27 to 29) commenced using the default initialization of SWAN, which regarded the initial wave height as zero. In the following 33 daily forecasts, the initial files of SWAN were derived from the restart files of simulations three days earlier.The analysis period in this study ran from UTC 0 : 00 on August 1 to UTC 23 : 00 on August 31. The 72 h forecast simulations during the analysis period were divided into three groups, 24-hour, 48-hour, and 72-hour forecast results, based on their forecast lead times. The three groups of forecast simulations are evaluated separately in the following sections.The model domains, as well as their topography and bathymetry, are shown in Figure2. The outer and inner domains are named D01 and D02, respectively. The land area of the inner domain with a spatial resolution of 6 km × 6 km covers Hainan Island, Leizhou Peninsula, and the southern coastal areas of Guangxi and Guangdong Province in South China. The elevation in most land areas is below 500 m, except for the southern mountains of Hainan Island, with an elevation of 500–1500 m. The sea area of D02 is located in the northwestern corner of the South China Sea. The water depth increases from northwest to southeast, with a maximum depth of 3000–3500 m.Figure 2
(a) Topography of the WRF submodel; (b) bathymetry for the ROMS and SWAN submodel. The entire colored domain is D01 and the area enclosed by the red frame is D02.
(a)(b)Owing to space- and observation-related constraints, this study only focuses on the evaluation of surface meteorological and hydrological elements in the D02 domain.


## 2.1. Model Description

The COAWST model includes one atmosphere model, the Weather Research and Forecasting (WRF) model; one ocean model, the Regional Ocean Modeling System (ROMS); one wave model, the Simulating Waves Nearshore (SWAN) model; and one sediment-transport model, the Community Sediment Transport Model (CSTM).The WRF model is a nonhydrostatic, quasi-compressible, mesoscale atmosphere model with a number of physical schemes. In the COAWST model, the dynamic core of the Advanced Research WRF (ARW) model is used [36]. The ROMS model is a free-surface, topography-following-coordinate model that solves the three-dimensional Reynolds-averaged Navier-Stokes equations using hydrostatic and Boussinesq approximations [37, 38]. The SWAN model is a spectral-wave model that simulates the generation and propagation of wind-induced waves in coastal waters by solving the spectral-density evolution equation [39]. The CSTM is a sediment-hydrodynamic model that simulates the erosion, deposition, and transport processes across different sizes and types of sediments [24]. In this study, the CSTM was not utilized.The COAWST model utilizes the model-coupling toolkit (MCT) to achieve communication among the submodels via the message passing interface (MPI) [40]. In addition, the Spherical Coordinate Remapping Interpolation Package (SCRIP) is used to support variable interpolations in different coordinate systems [41]. In the coupling process, WRF receives sea-surface roughness from SWAN and the SST from ROMS, while it provides wind at a 10 m level (W10 m) to SWAN and surface stress and net heat fluxes to ROMS. To facilitate this variable exchange between ROMS and SWAN, ROMS receives the surface- and bottom-wave direction, height, length, period, percentage breaking, energy dissipation, and bottom orbital velocity from SWAN, while it provides bathymetry, bottom elevation, sea-surface height, and depth-averaged currents to SWAN.


## 2.2. Forecasting System Flow Description

Based on the COAWST model, an experimental forecasting system was established in this study, the system flowchart of which is shown in Figure1. The system is controlled by executable shell scripts with timers and consists of one data server, one compute server, and the user client.Figure 1
Flowchart of the developed real-time forecast system.The data server downloads and preprocesses the global forecast data every day from the public data server of the National Centers for Environmental Prediction (NCEP). The downloaded data includes the Global Forecast System (GFS) data and the Global Data Assimilation System (GDAS) data for WRF (available online:https://www.nco.ncep.noaa.gov/pmb/products/gfs/), the Real-Time Ocean Forecast System (RTOFS) data for ROMS (available online: https://www.nco.ncep.noaa.gov/pmb/products/rtofs), and NCEP operational wave-product data from the global WAVEWATCH III model for SWAN (available online: https://www.nco.ncep.noaa.gov/pmb/products/wave/). The resolution information of the downloaded data is listed in Table 1.Table 1
Resolution information of the downloaded forecast data.
GFSRTOFSWaveSpatial resolution0.25° × 0.25°0.083° × 0.083°0.5° × 0.5°Time resolution3 hours3 hours3 hoursStart timeUTC 0 : 00UTC 0 : 00UTC 0 : 00Time rangeSubsequent 72 hoursSubsequent 72 hoursSubsequent 72 hoursIn the compute server, the GFS and GDAS data are processed by the WRF Preprocessing System (WPS) and 3-dimensional Variational data assimilation system (3DVar), respectively, to produce WRF’s initial and boundary conditions. The GFS data, RTOFS data, and NCEP operational wave-product data are processed in the ROMS and SWAN preprocessing module to produce the boundary files for SWAN, along with the initial and boundary files of for ROMS.In this study, the initial file for SWAN was provided by the previous restart file from 72 hours earlier—instead of the default idealized initialization—as many previous studies on the COAWST model have done [26, 42, 43]. The bathymetry dataset for ROMS used the DTU10 dataset, with a resolution of 0.125° × 0.125°, from the Technical University of Denmark (available online: https://www.space.dtu.dk/english/Research/Scientific_data_and_models/Global_Bathymetry_Model). The tide dataset for ROMS was derived from the TPXO7.0 dataset, with a resolution of 0.25° × 0.25°, which was obtained using the Oregon State University Tidal Prediction Software (OTPS). The TPXO global tide model provides eight primary (M2, S2, N2, K2, K1, O1, P1, and Q1) and two long (Mf and Mm) period harmonic constituents as amplitudes of the earth-relative sea-surface elevation [44, 45]. The TPXO7.0 dataset can be found online at http://www−po.coas.oregonstate.edu/∼poa/www−po/research/po/research/tide/global.html.The initial and boundary files for WRF, ROMS, and SWAN are transferred to the COAWST model for forecast simulations; subsequently, the simulation results are processed using a postprocessing module and are visualized for the end user.


## 2.3. Experimental Design

The configuration parameters for model gridding and operation are listed in Table2. The COAWST model in this study used two nested grids to obtain results with a higher resolution via dynamic downscaling. The model was compiled with fully coupling options and three submodels exchanged variables every 1800 seconds, controlled by the main program via the MCT.Table 2
Configuration parameters for each submodel.
WRFROMSSWANTime step30 s60 s180 sGrid nestingYesYesYesOuter grid number100 × 10090 × 9090×90Inner grid number100 × 10090 × 9090×90Horizontal grid resolution18 km for outer grids, 6 km for inner grids18 km for outer grids, 6 km for inner grids18 km for outer grids, 6 km for inner gridsVertical layer number3916NoneInitial data sourceGFS data assimilated by GDASGFS and RTOFS dataRestart file of simulation 3 days earlierVariable exchange frequency1800 s−11800 s−11800 s−1The default scheme options in the COAWST model were used without optimization in this study. The schemes used in WRF were the WRF Single-Moment 3-class (WSM3) microphysics scheme [46], the Rapid Radiative Transfer Model (RRTM) longwave scheme [47], the Dudhia shortwave scheme [48], the Noah Land Surface Model (LSM) [49], the modified Kain-Fritsch cumulus scheme [50], the Eta surface layer scheme [51], and the Mellor-Yamada-Janjic (MYJ) planetary boundary-layer scheme [52, 53]. For ROMS setup, the Mellor-Yamada scheme was used to compute vertical turbulent mixing [53], and the Flather boundary-condition method was used for barotropic currents to allow the free propagation of wind-generated currents and tides [54]. For the SWAN model, the Madsen scheme was used for wave bottom dissipation parameterization [55], and the Komen formulation was used to simulate wind-induced wave growth [56].The forecast simulation tests began at UTC 0 : 00 every day with a forecast lead time of 72 hours and lasted from July 27 to August 31, 2019. The first three tests (from July 27 to 29) commenced using the default initialization of SWAN, which regarded the initial wave height as zero. In the following 33 daily forecasts, the initial files of SWAN were derived from the restart files of simulations three days earlier.The analysis period in this study ran from UTC 0 : 00 on August 1 to UTC 23 : 00 on August 31. The 72 h forecast simulations during the analysis period were divided into three groups, 24-hour, 48-hour, and 72-hour forecast results, based on their forecast lead times. The three groups of forecast simulations are evaluated separately in the following sections.The model domains, as well as their topography and bathymetry, are shown in Figure2. The outer and inner domains are named D01 and D02, respectively. The land area of the inner domain with a spatial resolution of 6 km × 6 km covers Hainan Island, Leizhou Peninsula, and the southern coastal areas of Guangxi and Guangdong Province in South China. The elevation in most land areas is below 500 m, except for the southern mountains of Hainan Island, with an elevation of 500–1500 m. The sea area of D02 is located in the northwestern corner of the South China Sea. The water depth increases from northwest to southeast, with a maximum depth of 3000–3500 m.Figure 2
(a) Topography of the WRF submodel; (b) bathymetry for the ROMS and SWAN submodel. The entire colored domain is D01 and the area enclosed by the red frame is D02.
(a)(b)Owing to space- and observation-related constraints, this study only focuses on the evaluation of surface meteorological and hydrological elements in the D02 domain.


## 3. Results

### 3.1. Evaluation Based on Station Observations

The hourly observations of air temperature at a 2 m level (T2 m), the air relative humidity at a 2 m level (RH2 m), and sea level pressure (SLP) from 29 stations over the study domain were derived from China Meteorological Data Service Center (CMDC) (http://data.cma.cn/en). The locations of the 29 stations are shown in Figure 3(a), and the Taylor diagrams of forecasted T2 m, RH2 m, and SLP are shown in Figures 3(b)–3(d).Figure 3
(a) Locations of the meteorological stations. The stations colored in red dots are the stations where the NSDs of T2 m and RH2 m for the 48 h and 72 h forecasts are higher than 0.5; those in blue are the others (NSDs lower than 0.5). Taylor diagram of (b) T2 m, (c) RH2 m, and (d) SLP for 24 (h), 48 (h), and 72 h forecasts.
(a)(b)(c)(d)The scatter distributions of T2 m and RH2 m in the Taylor diagram (Figures3(b) and 3(c), respectively) are similar. The 24 h forecasts (in blue), which are the first 24 h results of the 72 h simulation conducted every day, have the highest accuracy. In comparison to hourly station observation, most of the 24 h forecast series at the corresponding grids have a temporal correlation coefficient (CC) of 0.8–0.9 and a normalized standard deviation (NSD) of around 0.75–1.0. However, the 48 h and 72 h forecasts (in red and green, respectively), which are produced one and two days before the analysis, respectively, have lower accuracies than the 24 h forecasts. In addition, there are no obvious differences in correlation and NSD for the 48 h and 72 h forecasts.As shown in Figures3(b) and 3(c), the scatter distributions of the 48 h and 72 h forecasts can be divided into two groups using the circle line of an NSD of 0.5. To facilitate differentiation, the stations with NSDs greater than 0.5 were colored in red in Figure 3(a), while those with NSDs lower than 0.5 were colored in blue. Most of the stations with NSDs lower than 0.5 were located in the southern Hainan Island, indicating that the 48 h and 72 h forecasts of T2 m and RH2 m vary in a lower range than those in observations.During the analysis period, the Hainan Island and Leizhou Peninsula were affected by Typhoon Wipha (20190731–20190802), which passed from the southeast to the northwest. Further, the inner domain was affected by the periphery of Typhoon Bailu (20190824–20190826) and Typhoon Podul (20190828–20190829). To confirm whether the poorer performances of 48 h and 72 h forecasts were related to the typhoon weather, the Taylor diagrams in Figure3 were regraphed using the data from August 4 to 16, when no strong convection events occurred (figure not shown). The new diagrams showed no obvious differences compared to the diagrams derived from the entire monthly data. Thus, irrespective of strong convection events, the 48 h and 72 h forecasts of T2 m and RH2 m perform more poorly than the 24 h forecast in NSD and CC over the southern part of Hainan Island.Unlike the T2 m and RH2 m, the SLP forecasts in Figure3(d) have a high concentration with an NSD of 1.0. The correlation coefficients drop from 0.95 to about 0.85 for the 24 h, 48 h, and 72 h forecasts.The diagram of precipitation is unavailable in this section, because significant data on precipitation is missing in the station observations. Instead, the reanalysis dataset is derived in the following section to evaluate the performance of T2 m and precipitation, which are two important elements of NWP.


### 3.2. Evaluation of T2 m and Precipitation Based on Reanalysis Datasets

In this section, two sets of reanalysis data, from ERA5 and CLDAS (Chinese Land Data Assimilation System), were derived for the evaluation of T2 m and precipitation. ERA5 is the fifth generation of reanalysis datasets produced by the European Center for Medium-Range Weather Forecasts (ECMWF), with a spatial resolution of 0.25° × 0.25° and a temporal resolution of 1 hour. The dataset is available online athttps://cds.climate.copernicus.eu/#!/search?text=ERA5&type=dataset. The CLDAS data is produced by the National Meteorological Information Center of China and has a spatial resolution of 0.0625° × 0.0625° in Asia and a temporal resolution of 1 hour. The data is available online at http://data.cma.cn/data/cdcdetail/dataCode/NAFP_CLDAS2.0_NRT.html. According to related evaluation studies, the CLDAS data had better performance when compared with the other reanalysis products in the region of China [57–59].


#### 3.2.1. Spatial Distributions of Mean Bias

Figure4 shows the spatial distributions of the mean precipitation bias between the forecasts and reanalysis data. All of the data was interpolated into the same resolution of 0.05° × 0.05° through bilinear interpolation. As shown in Figure 4, all of the forecasts indicate more precipitation in the land grids than the reanalysis data. For the bias of forecasts minus ERA5 in the sea grids (Figures 4(a)–4(c)), lower values are detected, especially in offshore areas near the northern continental region. As for the southern open-sea areas, the forecast bias is weak for both sets of reanalysis data. The mean biases averaged over the entire domain are –0.01 mm/h, 0.08 mm/h, and 0.03 mm/h.Figure 4
Spatial distributions of mean precipitation biases of (a) 24 (h), (b) 48 (h), and (c) 72 h forecasts minus ERA5 data and biases of (d) 24 (h), (e) 48 (h), and (f) 72 h forecasts minus CLDAS data.
(a)(b)(c)(d)(e)(f)In Figures4(d)–4(f), the bias distributions of forecasts minus CLDAS are similar to those in Figures 4(a)–4(c). However, the land-sea differences in Figures 4(d)–4(f) are not as strong as those in Figures 4(a)–4(c), and a higher bias is detected over the northern and eastern sea areas. The mean biases in Figures 4(d)–4(f) are 0.11 mm/h, 0.19 mm/h, and 0.14 mm/h, which are slightly higher than those in Figures 4(a)–4(c). The bias distributions of precipitation indicate that the performances of ERA5 and CLDAS reanalysis data have greater differences in the sea grids, especially in offshore regions.Further, the bias of 24 h forecasts is weaker than that of the others, while the bias distributions of 48 h and 72 h forecasts show no obvious differences.The spatial distributions of mean T2 m bias are shown in Figure5. Similar bias distributions—but with more obvious land-sea differences than those in Figure 4— are detected; that is, a lower bias is present in the land grids and a higher bias is present in the sea grids, especially in the northern offshore regions. Over the domain, the mean T2 m of seas is about 2–3 K higher than that over land in August (figure not shown). Mean lower surface pressure center locates in the northwestern area of the domain, while the higher-pressure center locates in the southern sea areas. The distributions of cooler land and warmer seas for the forecasts will enlarge the land-sea thermal differences and the surface pressure gradient from northwest to southeast. As a result, the prevailing South China Sea summer monsoon, which flows from southern seas, will be strengthened. More moisture will also be transported to the northern continental, inducing more precipitation there.Figure 5
Spatial distributions of T2 m precipitation biases of (a) 24 (h), (b) 48 (h), and (c) 72 h forecasts minus ERA5 data and biases of (d) 24 (h), (e) 48 (h), and (f) 72 h forecasts minus CLDAS data.
(a)(b)(c)(d)(e)(f)In Figure5, it is important to note that the two sets of reanalysis data do not coincide with each other at all points. In southern areas of Hainan Island, for example, the mean T2 m differences between the ERA5 and CLDAS data exceed 1 K, leading to a higher spatial heterogeneity. In general, the mean forecasts are closer to the CLDAS data. The mean biases of forecasts minus CLDAS are 0.09 K, 0.05 K, and 0.06 K, while the biases of forecasts minus ERA5 are 0.22 K, 0.18 K, and 0.19 K.In addition, the mean bias distributions of the three forecasts do not differ greatly owing to their differing forecast lead times, and the slight differences among the three forecasts can be accounted for by the greater land-sea differences during the averaging process.


#### 3.2.2. Time Variations of Bias

Figures6 and 7 show the hourly bias series of precipitation and T2 m, respectively. Owing to the land-sea differences in bias, Figures 6 and 7 adopted the series averaged in the land and sea grids separately, instead of the ones averaged over the entire domain.Figure 6
Precipitation bias series of 24 h forecasts averaged over (a) land grids and (b) sea grids; bias series of 48 h forecasts averaged over (c) land grids and (d) sea grids; and bias series of 72 h forecasts averaged over (e) land grids and (f) sea grids.
(a)(b)(c)(d)(e)(f)Figure 7
T2 m bias series of 24 h forecasts averaged over (a) land grids and (b) sea grids; bias series of 48 h forecasts averaged over (c) land grids and (d) sea grids; and bias series of 72 h forecasts averaged over (e) land grids and (f) sea grids.
(a)(b)(c)(d)(e)(f)As shown in Figure6, the precipitation bias series of forecasts minus ERA5 are in accordance with the series of forecasts minus CLDAS. During the analysis period, the bias series of forecasts minus ERA5 were generally lower in amplitude than the series of forecasts minus CLDAS. Taking the series in the land grids (Figures 6(a), 6(c), and 6(e)) as an example, the mean biases of forecasts minus ERA5 are 0.24 mm/h, 0.30 mm/h, and 0.21 mm/h for the 24 h, 48 h, and 72 h forecast series, respectively, while the mean biases of forecasts minus CLDAS are 0.28 mm/h, 0.34 mm/h, and 0.25 mm/h.The amplitudes of bias series in the sea grids are smaller than those in the land grids. For Figures6(b), 6(d), and 6(f), the mean biases of forecasts minus ERA5 are –0.09 mm/h, 0.01 mm/h, and –0.05 mm/h for the series of forecasts minus ERA5, while the biases of forecasts minus CLDAS are 0.05 mm/h, 0.14 mm/h, and 0.09 mm/h. The precipitation forecast series are closer to the ERA5 data than the CLDAS data, irrespective of whether they were over land or sea.Figure7 shows the bias series of T2 m in land and sea grids. For the land grids in Figures 7(a), 7(c), and 7(e), the bias series are majorly negative. The series of forecasts minus ERA5 have higher amplitudes and greater mean biases in comparison to the series of forecasts minus CLDAS, which is in accordance with the spatial distributions in Figure 5. The mean biases of forecasts minus ERA5 in land grids are −0.38 K, −0.33 K, and −0.24 K, while the biases of forecasts minus CLDAS are −0.27 K, −0.21 K, and −0.12 K.For the sea grids in Figures7(b), 7(d), and 7(f), the forecasted T2 m is overall higher than that in the ERA5 and CLDAS data. The mean biases of the two series are 0.44 K, 0.37 K, and 0.37 K for forecasts minus ERA5 and 0.24 K, 0.17 K, and 0.17 K for forecasts minus CLDAS. The largest T2 m bias appears during the period of 2019-08-20–2019-08-25, when a lower precipitation bias is detected in Figures 6(b), 6(d), and 6(f). In fact, the time when the extreme values of T2 m bias appear in Figures 7(b), 7(d), and 7(f) roughly coincides with the time when the extreme values of precipitation bias appear in Figures 6(b), 6(d), and 6(f). The bias series of T2 m in the sea grids indicate that the COAWST forecast system still needs further improvement for simulating the quick processes of local cooling in strong convective weather over the sea. This drawback of COAWST can also be found in other studies. For example, according to the hindcast results comparing observations from a buoy made by Liu et al. [26], the simulated T2 m by the COAWST model was around 1.5 K–3.5 K higher when a typhoon passed by the sea area near this buoy, while the bias reduced to less than 1 K after the typhoon.


#### 3.2.3. Statistical Indices over the Domain

The statistical indices averaged over the entire domain are presented in Table3, and the indices in the land and sea grids are also presented in Tables 4 and 5, respectively. The temporal CC, root mean square error (RMSE), and standard deviation (STD) are employed for this evaluation. As shown in Table 3, the indices of precipitation and T2 m reduce with increasing forecast lead times. The indices derived from forecasts and ERA5 data are slightly better than those from forecasts and CLDAS data in general; however, the STDs of forecast tests are closer to those of CLDAS data. In addition, the STDs of forecasted T2 m are generally lower than those for the reanalysis data over the entire domain, indicating that the forecasted T2 m has lower deviations and slower amplitude changes in the diurnal cycle.Table 3
Statistical indices of forecasts averaged over the entire domain.
ElementForecast testObservationCCRMSESTDForecastObservationPrecipitation24 h forecastERA50.281.24 mm/h1.15 mm/h0.82 mm/hCLDAS0.231.53 mm/h−1.18 mm/h48 h forecastERA50.161.35 mm/h1.25 mm/h−CLDAS0.181.65 mm/h−−72 h forecastERA50.121.34 mm/h1.15 mm/h−CLDAS0.071.64 mm/h−−2 m air temperature24 h forecastERA50.621.02 K1.03 K1.20 KCLDAS0.591.00 K−1.28 K48 h forecastERA50.561.06 K0.93 K−CLDAS0.541.06 K−−72 h forecastERA50.491.11 K0.90 K−CLDAS0.481.12 K−−CC: correlation coefficient;CC=∑i=1nxmi−xm¯xoi−xo¯/∑i=1nxmi−xm¯2⋅∑i=1nxoi−xo¯2; xmi and xoi are the simulated and observed hourly precipitation, respectively, or 2 m air temperature in the time i; and n is the number of hours. RMSE: Root Mean Square Error; RMSE=1/n∑i=1nxmi−xoi2. STD: standard deviation; STD=1/n∑i=1nxi−x¯2.Table 4
Statistical indices of forecasts averaged over land grids.
ElementForecast testObservationCCRMSESTDForecastObservationPrecipitation24 h forecastERA50.321.67 mm/h1.67 mm/h0.96 mm/hCLDAS0.261.95 mm/h−1.41 mm/h48 h forecastERA50.301.69 mm/h1.65 mm/h−CLDAS0.231.97 mm/h−−72 h forecastERA50.221.66 mm/h1.50 mm/h−CLDAS0.171.93 mm/h−−2 m air temperature24 h forecastERA50.821.33 K2.03 K1.90 KCLDAS0.811.46 K−2.38 K48 h forecastERA50.741.44 K1.88 K−CLDAS0.731.65 K−−72 h forecastERA50.701.52 K1.85 K−CLDAS0.691.74 K−−Table 5
Statistical indices of forecasts averaged over sea grids.
ElementForecast testObservationCCRMSESTDForecastObservationPrecipitation24 h forecastERA50.271.11 mm/h0.99 mm/h0.79 mm/hCLDAS0.161.44 mm/h−1.15 mm/h48 h forecastERA50.191.25 mm/h1.12 mm/h−CLDAS0.091.59 mm/h−−72 h forecastERA50.141.22 mm/h1.01 mm/h−CLDAS0.031.56 mm/h−−2 m air temperature24 h forecastERA50.560.89 K0.62 K0.91 KCLDAS0.520.82 K−0.86 K48 h forecastERA50.500.91 K0.55 K−CLDAS0.470.83 K−−72 h forecastERA50.420.96 K0.51 K−CLDAS0.410.89 K−−For the indices averaged over land grids (Table4), the CCs, RMSEs, and STDs are overall higher than the domain-averaged indices in Table 3. The averaged CCs of T2 m reach 0.7–0.8, which is in accordance with the station observations in Figure 3(b). The STDs of T2 m are obviously higher than those in Table 3, because of the more intense diurnal cycle on land. The indices in sea grids (Table 5) are generally lower than those in Table 3, in contrast to the values for land grids found in Table 4. This indicates that the forecasts in sea grids have lower mean biases, flatter fluctuations, and weaker linear correlations with reanalysis data, in comparison to the indices in Tables 3 and 4.


### 3.3. Evaluation of Other Surface Elements

In addition to the precipitation and T2 m, the spatial distributions of mean bias for other surface elements are also presented in Figures8 and 9; these elements include RH2 m, SLP, wind at a 10 m level (W10 m), significant wave height (Hwave), SST, and surface seawater salinity. All observation data for comparison was obtained from the ERA5 reanalysis dataset, other than the salinity data, which was derived from the daily ocean reanalysis data using the GLO-CPL weakly coupled ocean-atmosphere data assimilation system of the Met Office, UK [60]. The salinity reanalysis data, with a spatial resolution of 0.25° × 0.25°, is available online at https://resources.marine.copernicus.eu/?option=com_csw&task=results. The spatial resolution of the Hwave reanalysis data is 0.5° × 0.5=°, while the resolutions of other elements are 0.25° × 0.25°.Figure 8
Spatial distributions of SST biases of (a) 24 (h), (b) 48 (h), and (c) 72 h forecasts minus ERA5 data; SLP biases of (d) 24 (h), (e) 48 (h), and (f) 72 h forecasts minus ERA5 data; and W10 m biases of (g) 24 (h), (h) 48 (h), and (i) 72 h forecasts minus ERA5 data.
(a)(b)(c)(d)(e)(f)(g)(h)(i)Figure 9
Spatial distributions of RH2 m biases of (a) 24 (h), (b) 48 (h), and (c) 72 h forecasts minus ERA5 data; Hwave biases of (d) 24 (h), (e) 48 (h), and (f) 72 h forecasts minus ERA5 data; and surface seawater salinity biases of (g) 24 (h), (h) 48 (h), and (i) 72 h forecasts minus ERA5 data.
(a)(b)(c)(d)(e)(f)(g)(h)(i)As shown in Figures8(a)–8(c), the distributions of SST are similar to the distributions of T2 m in Figure 5. Higher bias is detected in the northern offshore regions and lower bias is detected in the surrounding offshore regions of Hainan Island. The mean bias in the southern open seas, however, is closer to zero. The mean bias for the three forecast tests is 0.03 K, 0.08 K, and 0.08 K over the domain, which is weaker than the T2 m bias.The simulated thermal differences have induced the changes of surface pressure. For the SLP bias in Figures8(d)–8(f), the bias is overall negative and decreases from the northwest to the southeast direction of the domain. The mean bias is −0.33 hPa, −0.71 hPa, and −1.15 hPa, respectively, for the 24 h, 48 h, and 72 h forecasts. The lowest bias is located in the northwestern areas of the domain. In the month of August, as referred above, a cooler and lower-pressure center locates in the northwestern areas of the domain on average and the warmer and higher-pressure center locates in the southeast. The distribution of SLP bias induces more pressure dropping over the low-pressure area, enlarging the surface horizontal pressure gradient over the domain. The changes of pressure gradient will further enhance the prevailing south wind and transport more moisture from the sea to the land.Figures8(g)–8(i) show the W10 m bias between forecasts and ERA5 data. A bias of the south wind is detected over the domain, and the mean bias of wind speed is 0.58 m/s, 1.30 m/s, and 1.64 m/s. The south wind bias enhances the moisture transport from the southern seas and induces higher local humidity. As shown in Figures 9(a)–9(c), the forecasted RH2 m is generally higher than that in the ERA5 data. The mean bias is 3.10%, 4.14%, and 4.00%. Higher bias is detected in the southern sea grids than in the northern land grids.The bias of Hwave (Figures9(d)–9(f)) is slightly lower than that in the ERA5 data and decreases roughly from the southwest to the northeast. The mean bias over the domain is −0.12 m, −0.06 m, and −0.01 m. The Hwave here is denoted as the combined height of wind-induced waves with a high frequency and swells with a low frequency. The wind-induced wave height is known to be positively correlated with the wind speed at the sea surface, and the mean bias of wind-induced wave height is 0.34 m, 0.38 m, and 0.42 m for the 24 h, 48 h, and 72 h forecasts, respectively (figure not shown). The overall higher bias of wind-induced wave height indicates that the low Hwave bias results from the offsetting of a higher bias of wind-induced wave height and a lower one of swell height.The ERA5 data has no ocean products; therefore, Figures9(g)–9(i) show the bias distributions of surface seawater salinity between the forecasts and the ocean reanalysis data from the Met Office. The mean bias is 0.31, 0.33, and 0.35. The distributions of salinity bias do not differ in the three forecast tests, and the bias in the southern open seas is obviously weaker than that in the northern offshore regions. The highest bias of simulated salinity appears in the northeastern corner of the domain, which is close to the mouth area of the Pearl River, a river with the second largest discharge in China. In the COAWST model, the freshwater horizontal input is inactive by default, which is probably the primary reason for the bias distributions and requires further improvements in the future.


## 3.1. Evaluation Based on Station Observations

The hourly observations of air temperature at a 2 m level (T2 m), the air relative humidity at a 2 m level (RH2 m), and sea level pressure (SLP) from 29 stations over the study domain were derived from China Meteorological Data Service Center (CMDC) (http://data.cma.cn/en). The locations of the 29 stations are shown in Figure 3(a), and the Taylor diagrams of forecasted T2 m, RH2 m, and SLP are shown in Figures 3(b)–3(d).Figure 3
(a) Locations of the meteorological stations. The stations colored in red dots are the stations where the NSDs of T2 m and RH2 m for the 48 h and 72 h forecasts are higher than 0.5; those in blue are the others (NSDs lower than 0.5). Taylor diagram of (b) T2 m, (c) RH2 m, and (d) SLP for 24 (h), 48 (h), and 72 h forecasts.
(a)(b)(c)(d)The scatter distributions of T2 m and RH2 m in the Taylor diagram (Figures3(b) and 3(c), respectively) are similar. The 24 h forecasts (in blue), which are the first 24 h results of the 72 h simulation conducted every day, have the highest accuracy. In comparison to hourly station observation, most of the 24 h forecast series at the corresponding grids have a temporal correlation coefficient (CC) of 0.8–0.9 and a normalized standard deviation (NSD) of around 0.75–1.0. However, the 48 h and 72 h forecasts (in red and green, respectively), which are produced one and two days before the analysis, respectively, have lower accuracies than the 24 h forecasts. In addition, there are no obvious differences in correlation and NSD for the 48 h and 72 h forecasts.As shown in Figures3(b) and 3(c), the scatter distributions of the 48 h and 72 h forecasts can be divided into two groups using the circle line of an NSD of 0.5. To facilitate differentiation, the stations with NSDs greater than 0.5 were colored in red in Figure 3(a), while those with NSDs lower than 0.5 were colored in blue. Most of the stations with NSDs lower than 0.5 were located in the southern Hainan Island, indicating that the 48 h and 72 h forecasts of T2 m and RH2 m vary in a lower range than those in observations.During the analysis period, the Hainan Island and Leizhou Peninsula were affected by Typhoon Wipha (20190731–20190802), which passed from the southeast to the northwest. Further, the inner domain was affected by the periphery of Typhoon Bailu (20190824–20190826) and Typhoon Podul (20190828–20190829). To confirm whether the poorer performances of 48 h and 72 h forecasts were related to the typhoon weather, the Taylor diagrams in Figure3 were regraphed using the data from August 4 to 16, when no strong convection events occurred (figure not shown). The new diagrams showed no obvious differences compared to the diagrams derived from the entire monthly data. Thus, irrespective of strong convection events, the 48 h and 72 h forecasts of T2 m and RH2 m perform more poorly than the 24 h forecast in NSD and CC over the southern part of Hainan Island.Unlike the T2 m and RH2 m, the SLP forecasts in Figure3(d) have a high concentration with an NSD of 1.0. The correlation coefficients drop from 0.95 to about 0.85 for the 24 h, 48 h, and 72 h forecasts.The diagram of precipitation is unavailable in this section, because significant data on precipitation is missing in the station observations. Instead, the reanalysis dataset is derived in the following section to evaluate the performance of T2 m and precipitation, which are two important elements of NWP.


## 3.2. Evaluation of T2 m and Precipitation Based on Reanalysis Datasets

In this section, two sets of reanalysis data, from ERA5 and CLDAS (Chinese Land Data Assimilation System), were derived for the evaluation of T2 m and precipitation. ERA5 is the fifth generation of reanalysis datasets produced by the European Center for Medium-Range Weather Forecasts (ECMWF), with a spatial resolution of 0.25° × 0.25° and a temporal resolution of 1 hour. The dataset is available online athttps://cds.climate.copernicus.eu/#!/search?text=ERA5&type=dataset. The CLDAS data is produced by the National Meteorological Information Center of China and has a spatial resolution of 0.0625° × 0.0625° in Asia and a temporal resolution of 1 hour. The data is available online at http://data.cma.cn/data/cdcdetail/dataCode/NAFP_CLDAS2.0_NRT.html. According to related evaluation studies, the CLDAS data had better performance when compared with the other reanalysis products in the region of China [57–59].


### 3.2.1. Spatial Distributions of Mean Bias

Figure4 shows the spatial distributions of the mean precipitation bias between the forecasts and reanalysis data. All of the data was interpolated into the same resolution of 0.05° × 0.05° through bilinear interpolation. As shown in Figure 4, all of the forecasts indicate more precipitation in the land grids than the reanalysis data. For the bias of forecasts minus ERA5 in the sea grids (Figures 4(a)–4(c)), lower values are detected, especially in offshore areas near the northern continental region. As for the southern open-sea areas, the forecast bias is weak for both sets of reanalysis data. The mean biases averaged over the entire domain are –0.01 mm/h, 0.08 mm/h, and 0.03 mm/h.Figure 4
Spatial distributions of mean precipitation biases of (a) 24 (h), (b) 48 (h), and (c) 72 h forecasts minus ERA5 data and biases of (d) 24 (h), (e) 48 (h), and (f) 72 h forecasts minus CLDAS data.
(a)(b)(c)(d)(e)(f)In Figures4(d)–4(f), the bias distributions of forecasts minus CLDAS are similar to those in Figures 4(a)–4(c). However, the land-sea differences in Figures 4(d)–4(f) are not as strong as those in Figures 4(a)–4(c), and a higher bias is detected over the northern and eastern sea areas. The mean biases in Figures 4(d)–4(f) are 0.11 mm/h, 0.19 mm/h, and 0.14 mm/h, which are slightly higher than those in Figures 4(a)–4(c). The bias distributions of precipitation indicate that the performances of ERA5 and CLDAS reanalysis data have greater differences in the sea grids, especially in offshore regions.Further, the bias of 24 h forecasts is weaker than that of the others, while the bias distributions of 48 h and 72 h forecasts show no obvious differences.The spatial distributions of mean T2 m bias are shown in Figure5. Similar bias distributions—but with more obvious land-sea differences than those in Figure 4— are detected; that is, a lower bias is present in the land grids and a higher bias is present in the sea grids, especially in the northern offshore regions. Over the domain, the mean T2 m of seas is about 2–3 K higher than that over land in August (figure not shown). Mean lower surface pressure center locates in the northwestern area of the domain, while the higher-pressure center locates in the southern sea areas. The distributions of cooler land and warmer seas for the forecasts will enlarge the land-sea thermal differences and the surface pressure gradient from northwest to southeast. As a result, the prevailing South China Sea summer monsoon, which flows from southern seas, will be strengthened. More moisture will also be transported to the northern continental, inducing more precipitation there.Figure 5
Spatial distributions of T2 m precipitation biases of (a) 24 (h), (b) 48 (h), and (c) 72 h forecasts minus ERA5 data and biases of (d) 24 (h), (e) 48 (h), and (f) 72 h forecasts minus CLDAS data.
(a)(b)(c)(d)(e)(f)In Figure5, it is important to note that the two sets of reanalysis data do not coincide with each other at all points. In southern areas of Hainan Island, for example, the mean T2 m differences between the ERA5 and CLDAS data exceed 1 K, leading to a higher spatial heterogeneity. In general, the mean forecasts are closer to the CLDAS data. The mean biases of forecasts minus CLDAS are 0.09 K, 0.05 K, and 0.06 K, while the biases of forecasts minus ERA5 are 0.22 K, 0.18 K, and 0.19 K.In addition, the mean bias distributions of the three forecasts do not differ greatly owing to their differing forecast lead times, and the slight differences among the three forecasts can be accounted for by the greater land-sea differences during the averaging process.


### 3.2.2. Time Variations of Bias

Figures6 and 7 show the hourly bias series of precipitation and T2 m, respectively. Owing to the land-sea differences in bias, Figures 6 and 7 adopted the series averaged in the land and sea grids separately, instead of the ones averaged over the entire domain.Figure 6
Precipitation bias series of 24 h forecasts averaged over (a) land grids and (b) sea grids; bias series of 48 h forecasts averaged over (c) land grids and (d) sea grids; and bias series of 72 h forecasts averaged over (e) land grids and (f) sea grids.
(a)(b)(c)(d)(e)(f)Figure 7
T2 m bias series of 24 h forecasts averaged over (a) land grids and (b) sea grids; bias series of 48 h forecasts averaged over (c) land grids and (d) sea grids; and bias series of 72 h forecasts averaged over (e) land grids and (f) sea grids.
(a)(b)(c)(d)(e)(f)As shown in Figure6, the precipitation bias series of forecasts minus ERA5 are in accordance with the series of forecasts minus CLDAS. During the analysis period, the bias series of forecasts minus ERA5 were generally lower in amplitude than the series of forecasts minus CLDAS. Taking the series in the land grids (Figures 6(a), 6(c), and 6(e)) as an example, the mean biases of forecasts minus ERA5 are 0.24 mm/h, 0.30 mm/h, and 0.21 mm/h for the 24 h, 48 h, and 72 h forecast series, respectively, while the mean biases of forecasts minus CLDAS are 0.28 mm/h, 0.34 mm/h, and 0.25 mm/h.The amplitudes of bias series in the sea grids are smaller than those in the land grids. For Figures6(b), 6(d), and 6(f), the mean biases of forecasts minus ERA5 are –0.09 mm/h, 0.01 mm/h, and –0.05 mm/h for the series of forecasts minus ERA5, while the biases of forecasts minus CLDAS are 0.05 mm/h, 0.14 mm/h, and 0.09 mm/h. The precipitation forecast series are closer to the ERA5 data than the CLDAS data, irrespective of whether they were over land or sea.Figure7 shows the bias series of T2 m in land and sea grids. For the land grids in Figures 7(a), 7(c), and 7(e), the bias series are majorly negative. The series of forecasts minus ERA5 have higher amplitudes and greater mean biases in comparison to the series of forecasts minus CLDAS, which is in accordance with the spatial distributions in Figure 5. The mean biases of forecasts minus ERA5 in land grids are −0.38 K, −0.33 K, and −0.24 K, while the biases of forecasts minus CLDAS are −0.27 K, −0.21 K, and −0.12 K.For the sea grids in Figures7(b), 7(d), and 7(f), the forecasted T2 m is overall higher than that in the ERA5 and CLDAS data. The mean biases of the two series are 0.44 K, 0.37 K, and 0.37 K for forecasts minus ERA5 and 0.24 K, 0.17 K, and 0.17 K for forecasts minus CLDAS. The largest T2 m bias appears during the period of 2019-08-20–2019-08-25, when a lower precipitation bias is detected in Figures 6(b), 6(d), and 6(f). In fact, the time when the extreme values of T2 m bias appear in Figures 7(b), 7(d), and 7(f) roughly coincides with the time when the extreme values of precipitation bias appear in Figures 6(b), 6(d), and 6(f). The bias series of T2 m in the sea grids indicate that the COAWST forecast system still needs further improvement for simulating the quick processes of local cooling in strong convective weather over the sea. This drawback of COAWST can also be found in other studies. For example, according to the hindcast results comparing observations from a buoy made by Liu et al. [26], the simulated T2 m by the COAWST model was around 1.5 K–3.5 K higher when a typhoon passed by the sea area near this buoy, while the bias reduced to less than 1 K after the typhoon.


### 3.2.3. Statistical Indices over the Domain

The statistical indices averaged over the entire domain are presented in Table3, and the indices in the land and sea grids are also presented in Tables 4 and 5, respectively. The temporal CC, root mean square error (RMSE), and standard deviation (STD) are employed for this evaluation. As shown in Table 3, the indices of precipitation and T2 m reduce with increasing forecast lead times. The indices derived from forecasts and ERA5 data are slightly better than those from forecasts and CLDAS data in general; however, the STDs of forecast tests are closer to those of CLDAS data. In addition, the STDs of forecasted T2 m are generally lower than those for the reanalysis data over the entire domain, indicating that the forecasted T2 m has lower deviations and slower amplitude changes in the diurnal cycle.Table 3
Statistical indices of forecasts averaged over the entire domain.
ElementForecast testObservationCCRMSESTDForecastObservationPrecipitation24 h forecastERA50.281.24 mm/h1.15 mm/h0.82 mm/hCLDAS0.231.53 mm/h−1.18 mm/h48 h forecastERA50.161.35 mm/h1.25 mm/h−CLDAS0.181.65 mm/h−−72 h forecastERA50.121.34 mm/h1.15 mm/h−CLDAS0.071.64 mm/h−−2 m air temperature24 h forecastERA50.621.02 K1.03 K1.20 KCLDAS0.591.00 K−1.28 K48 h forecastERA50.561.06 K0.93 K−CLDAS0.541.06 K−−72 h forecastERA50.491.11 K0.90 K−CLDAS0.481.12 K−−CC: correlation coefficient;CC=∑i=1nxmi−xm¯xoi−xo¯/∑i=1nxmi−xm¯2⋅∑i=1nxoi−xo¯2; xmi and xoi are the simulated and observed hourly precipitation, respectively, or 2 m air temperature in the time i; and n is the number of hours. RMSE: Root Mean Square Error; RMSE=1/n∑i=1nxmi−xoi2. STD: standard deviation; STD=1/n∑i=1nxi−x¯2.Table 4
Statistical indices of forecasts averaged over land grids.
ElementForecast testObservationCCRMSESTDForecastObservationPrecipitation24 h forecastERA50.321.67 mm/h1.67 mm/h0.96 mm/hCLDAS0.261.95 mm/h−1.41 mm/h48 h forecastERA50.301.69 mm/h1.65 mm/h−CLDAS0.231.97 mm/h−−72 h forecastERA50.221.66 mm/h1.50 mm/h−CLDAS0.171.93 mm/h−−2 m air temperature24 h forecastERA50.821.33 K2.03 K1.90 KCLDAS0.811.46 K−2.38 K48 h forecastERA50.741.44 K1.88 K−CLDAS0.731.65 K−−72 h forecastERA50.701.52 K1.85 K−CLDAS0.691.74 K−−Table 5
Statistical indices of forecasts averaged over sea grids.
ElementForecast testObservationCCRMSESTDForecastObservationPrecipitation24 h forecastERA50.271.11 mm/h0.99 mm/h0.79 mm/hCLDAS0.161.44 mm/h−1.15 mm/h48 h forecastERA50.191.25 mm/h1.12 mm/h−CLDAS0.091.59 mm/h−−72 h forecastERA50.141.22 mm/h1.01 mm/h−CLDAS0.031.56 mm/h−−2 m air temperature24 h forecastERA50.560.89 K0.62 K0.91 KCLDAS0.520.82 K−0.86 K48 h forecastERA50.500.91 K0.55 K−CLDAS0.470.83 K−−72 h forecastERA50.420.96 K0.51 K−CLDAS0.410.89 K−−For the indices averaged over land grids (Table4), the CCs, RMSEs, and STDs are overall higher than the domain-averaged indices in Table 3. The averaged CCs of T2 m reach 0.7–0.8, which is in accordance with the station observations in Figure 3(b). The STDs of T2 m are obviously higher than those in Table 3, because of the more intense diurnal cycle on land. The indices in sea grids (Table 5) are generally lower than those in Table 3, in contrast to the values for land grids found in Table 4. This indicates that the forecasts in sea grids have lower mean biases, flatter fluctuations, and weaker linear correlations with reanalysis data, in comparison to the indices in Tables 3 and 4.


## 3.2.1. Spatial Distributions of Mean Bias

Figure4 shows the spatial distributions of the mean precipitation bias between the forecasts and reanalysis data. All of the data was interpolated into the same resolution of 0.05° × 0.05° through bilinear interpolation. As shown in Figure 4, all of the forecasts indicate more precipitation in the land grids than the reanalysis data. For the bias of forecasts minus ERA5 in the sea grids (Figures 4(a)–4(c)), lower values are detected, especially in offshore areas near the northern continental region. As for the southern open-sea areas, the forecast bias is weak for both sets of reanalysis data. The mean biases averaged over the entire domain are –0.01 mm/h, 0.08 mm/h, and 0.03 mm/h.Figure 4
Spatial distributions of mean precipitation biases of (a) 24 (h), (b) 48 (h), and (c) 72 h forecasts minus ERA5 data and biases of (d) 24 (h), (e) 48 (h), and (f) 72 h forecasts minus CLDAS data.
(a)(b)(c)(d)(e)(f)In Figures4(d)–4(f), the bias distributions of forecasts minus CLDAS are similar to those in Figures 4(a)–4(c). However, the land-sea differences in Figures 4(d)–4(f) are not as strong as those in Figures 4(a)–4(c), and a higher bias is detected over the northern and eastern sea areas. The mean biases in Figures 4(d)–4(f) are 0.11 mm/h, 0.19 mm/h, and 0.14 mm/h, which are slightly higher than those in Figures 4(a)–4(c). The bias distributions of precipitation indicate that the performances of ERA5 and CLDAS reanalysis data have greater differences in the sea grids, especially in offshore regions.Further, the bias of 24 h forecasts is weaker than that of the others, while the bias distributions of 48 h and 72 h forecasts show no obvious differences.The spatial distributions of mean T2 m bias are shown in Figure5. Similar bias distributions—but with more obvious land-sea differences than those in Figure 4— are detected; that is, a lower bias is present in the land grids and a higher bias is present in the sea grids, especially in the northern offshore regions. Over the domain, the mean T2 m of seas is about 2–3 K higher than that over land in August (figure not shown). Mean lower surface pressure center locates in the northwestern area of the domain, while the higher-pressure center locates in the southern sea areas. The distributions of cooler land and warmer seas for the forecasts will enlarge the land-sea thermal differences and the surface pressure gradient from northwest to southeast. As a result, the prevailing South China Sea summer monsoon, which flows from southern seas, will be strengthened. More moisture will also be transported to the northern continental, inducing more precipitation there.Figure 5
Spatial distributions of T2 m precipitation biases of (a) 24 (h), (b) 48 (h), and (c) 72 h forecasts minus ERA5 data and biases of (d) 24 (h), (e) 48 (h), and (f) 72 h forecasts minus CLDAS data.
(a)(b)(c)(d)(e)(f)In Figure5, it is important to note that the two sets of reanalysis data do not coincide with each other at all points. In southern areas of Hainan Island, for example, the mean T2 m differences between the ERA5 and CLDAS data exceed 1 K, leading to a higher spatial heterogeneity. In general, the mean forecasts are closer to the CLDAS data. The mean biases of forecasts minus CLDAS are 0.09 K, 0.05 K, and 0.06 K, while the biases of forecasts minus ERA5 are 0.22 K, 0.18 K, and 0.19 K.In addition, the mean bias distributions of the three forecasts do not differ greatly owing to their differing forecast lead times, and the slight differences among the three forecasts can be accounted for by the greater land-sea differences during the averaging process.


## 3.2.2. Time Variations of Bias

Figures6 and 7 show the hourly bias series of precipitation and T2 m, respectively. Owing to the land-sea differences in bias, Figures 6 and 7 adopted the series averaged in the land and sea grids separately, instead of the ones averaged over the entire domain.Figure 6
Precipitation bias series of 24 h forecasts averaged over (a) land grids and (b) sea grids; bias series of 48 h forecasts averaged over (c) land grids and (d) sea grids; and bias series of 72 h forecasts averaged over (e) land grids and (f) sea grids.
(a)(b)(c)(d)(e)(f)Figure 7
T2 m bias series of 24 h forecasts averaged over (a) land grids and (b) sea grids; bias series of 48 h forecasts averaged over (c) land grids and (d) sea grids; and bias series of 72 h forecasts averaged over (e) land grids and (f) sea grids.
(a)(b)(c)(d)(e)(f)As shown in Figure6, the precipitation bias series of forecasts minus ERA5 are in accordance with the series of forecasts minus CLDAS. During the analysis period, the bias series of forecasts minus ERA5 were generally lower in amplitude than the series of forecasts minus CLDAS. Taking the series in the land grids (Figures 6(a), 6(c), and 6(e)) as an example, the mean biases of forecasts minus ERA5 are 0.24 mm/h, 0.30 mm/h, and 0.21 mm/h for the 24 h, 48 h, and 72 h forecast series, respectively, while the mean biases of forecasts minus CLDAS are 0.28 mm/h, 0.34 mm/h, and 0.25 mm/h.The amplitudes of bias series in the sea grids are smaller than those in the land grids. For Figures6(b), 6(d), and 6(f), the mean biases of forecasts minus ERA5 are –0.09 mm/h, 0.01 mm/h, and –0.05 mm/h for the series of forecasts minus ERA5, while the biases of forecasts minus CLDAS are 0.05 mm/h, 0.14 mm/h, and 0.09 mm/h. The precipitation forecast series are closer to the ERA5 data than the CLDAS data, irrespective of whether they were over land or sea.Figure7 shows the bias series of T2 m in land and sea grids. For the land grids in Figures 7(a), 7(c), and 7(e), the bias series are majorly negative. The series of forecasts minus ERA5 have higher amplitudes and greater mean biases in comparison to the series of forecasts minus CLDAS, which is in accordance with the spatial distributions in Figure 5. The mean biases of forecasts minus ERA5 in land grids are −0.38 K, −0.33 K, and −0.24 K, while the biases of forecasts minus CLDAS are −0.27 K, −0.21 K, and −0.12 K.For the sea grids in Figures7(b), 7(d), and 7(f), the forecasted T2 m is overall higher than that in the ERA5 and CLDAS data. The mean biases of the two series are 0.44 K, 0.37 K, and 0.37 K for forecasts minus ERA5 and 0.24 K, 0.17 K, and 0.17 K for forecasts minus CLDAS. The largest T2 m bias appears during the period of 2019-08-20–2019-08-25, when a lower precipitation bias is detected in Figures 6(b), 6(d), and 6(f). In fact, the time when the extreme values of T2 m bias appear in Figures 7(b), 7(d), and 7(f) roughly coincides with the time when the extreme values of precipitation bias appear in Figures 6(b), 6(d), and 6(f). The bias series of T2 m in the sea grids indicate that the COAWST forecast system still needs further improvement for simulating the quick processes of local cooling in strong convective weather over the sea. This drawback of COAWST can also be found in other studies. For example, according to the hindcast results comparing observations from a buoy made by Liu et al. [26], the simulated T2 m by the COAWST model was around 1.5 K–3.5 K higher when a typhoon passed by the sea area near this buoy, while the bias reduced to less than 1 K after the typhoon.


## 3.2.3. Statistical Indices over the Domain

The statistical indices averaged over the entire domain are presented in Table3, and the indices in the land and sea grids are also presented in Tables 4 and 5, respectively. The temporal CC, root mean square error (RMSE), and standard deviation (STD) are employed for this evaluation. As shown in Table 3, the indices of precipitation and T2 m reduce with increasing forecast lead times. The indices derived from forecasts and ERA5 data are slightly better than those from forecasts and CLDAS data in general; however, the STDs of forecast tests are closer to those of CLDAS data. In addition, the STDs of forecasted T2 m are generally lower than those for the reanalysis data over the entire domain, indicating that the forecasted T2 m has lower deviations and slower amplitude changes in the diurnal cycle.Table 3
Statistical indices of forecasts averaged over the entire domain.
ElementForecast testObservationCCRMSESTDForecastObservationPrecipitation24 h forecastERA50.281.24 mm/h1.15 mm/h0.82 mm/hCLDAS0.231.53 mm/h−1.18 mm/h48 h forecastERA50.161.35 mm/h1.25 mm/h−CLDAS0.181.65 mm/h−−72 h forecastERA50.121.34 mm/h1.15 mm/h−CLDAS0.071.64 mm/h−−2 m air temperature24 h forecastERA50.621.02 K1.03 K1.20 KCLDAS0.591.00 K−1.28 K48 h forecastERA50.561.06 K0.93 K−CLDAS0.541.06 K−−72 h forecastERA50.491.11 K0.90 K−CLDAS0.481.12 K−−CC: correlation coefficient;CC=∑i=1nxmi−xm¯xoi−xo¯/∑i=1nxmi−xm¯2⋅∑i=1nxoi−xo¯2; xmi and xoi are the simulated and observed hourly precipitation, respectively, or 2 m air temperature in the time i; and n is the number of hours. RMSE: Root Mean Square Error; RMSE=1/n∑i=1nxmi−xoi2. STD: standard deviation; STD=1/n∑i=1nxi−x¯2.Table 4
Statistical indices of forecasts averaged over land grids.
ElementForecast testObservationCCRMSESTDForecastObservationPrecipitation24 h forecastERA50.321.67 mm/h1.67 mm/h0.96 mm/hCLDAS0.261.95 mm/h−1.41 mm/h48 h forecastERA50.301.69 mm/h1.65 mm/h−CLDAS0.231.97 mm/h−−72 h forecastERA50.221.66 mm/h1.50 mm/h−CLDAS0.171.93 mm/h−−2 m air temperature24 h forecastERA50.821.33 K2.03 K1.90 KCLDAS0.811.46 K−2.38 K48 h forecastERA50.741.44 K1.88 K−CLDAS0.731.65 K−−72 h forecastERA50.701.52 K1.85 K−CLDAS0.691.74 K−−Table 5
Statistical indices of forecasts averaged over sea grids.
ElementForecast testObservationCCRMSESTDForecastObservationPrecipitation24 h forecastERA50.271.11 mm/h0.99 mm/h0.79 mm/hCLDAS0.161.44 mm/h−1.15 mm/h48 h forecastERA50.191.25 mm/h1.12 mm/h−CLDAS0.091.59 mm/h−−72 h forecastERA50.141.22 mm/h1.01 mm/h−CLDAS0.031.56 mm/h−−2 m air temperature24 h forecastERA50.560.89 K0.62 K0.91 KCLDAS0.520.82 K−0.86 K48 h forecastERA50.500.91 K0.55 K−CLDAS0.470.83 K−−72 h forecastERA50.420.96 K0.51 K−CLDAS0.410.89 K−−For the indices averaged over land grids (Table4), the CCs, RMSEs, and STDs are overall higher than the domain-averaged indices in Table 3. The averaged CCs of T2 m reach 0.7–0.8, which is in accordance with the station observations in Figure 3(b). The STDs of T2 m are obviously higher than those in Table 3, because of the more intense diurnal cycle on land. The indices in sea grids (Table 5) are generally lower than those in Table 3, in contrast to the values for land grids found in Table 4. This indicates that the forecasts in sea grids have lower mean biases, flatter fluctuations, and weaker linear correlations with reanalysis data, in comparison to the indices in Tables 3 and 4.


## 3.3. Evaluation of Other Surface Elements

In addition to the precipitation and T2 m, the spatial distributions of mean bias for other surface elements are also presented in Figures8 and 9; these elements include RH2 m, SLP, wind at a 10 m level (W10 m), significant wave height (Hwave), SST, and surface seawater salinity. All observation data for comparison was obtained from the ERA5 reanalysis dataset, other than the salinity data, which was derived from the daily ocean reanalysis data using the GLO-CPL weakly coupled ocean-atmosphere data assimilation system of the Met Office, UK [60]. The salinity reanalysis data, with a spatial resolution of 0.25° × 0.25°, is available online at https://resources.marine.copernicus.eu/?option=com_csw&task=results. The spatial resolution of the Hwave reanalysis data is 0.5° × 0.5=°, while the resolutions of other elements are 0.25° × 0.25°.Figure 8
Spatial distributions of SST biases of (a) 24 (h), (b) 48 (h), and (c) 72 h forecasts minus ERA5 data; SLP biases of (d) 24 (h), (e) 48 (h), and (f) 72 h forecasts minus ERA5 data; and W10 m biases of (g) 24 (h), (h) 48 (h), and (i) 72 h forecasts minus ERA5 data.
(a)(b)(c)(d)(e)(f)(g)(h)(i)Figure 9
Spatial distributions of RH2 m biases of (a) 24 (h), (b) 48 (h), and (c) 72 h forecasts minus ERA5 data; Hwave biases of (d) 24 (h), (e) 48 (h), and (f) 72 h forecasts minus ERA5 data; and surface seawater salinity biases of (g) 24 (h), (h) 48 (h), and (i) 72 h forecasts minus ERA5 data.
(a)(b)(c)(d)(e)(f)(g)(h)(i)As shown in Figures8(a)–8(c), the distributions of SST are similar to the distributions of T2 m in Figure 5. Higher bias is detected in the northern offshore regions and lower bias is detected in the surrounding offshore regions of Hainan Island. The mean bias in the southern open seas, however, is closer to zero. The mean bias for the three forecast tests is 0.03 K, 0.08 K, and 0.08 K over the domain, which is weaker than the T2 m bias.The simulated thermal differences have induced the changes of surface pressure. For the SLP bias in Figures8(d)–8(f), the bias is overall negative and decreases from the northwest to the southeast direction of the domain. The mean bias is −0.33 hPa, −0.71 hPa, and −1.15 hPa, respectively, for the 24 h, 48 h, and 72 h forecasts. The lowest bias is located in the northwestern areas of the domain. In the month of August, as referred above, a cooler and lower-pressure center locates in the northwestern areas of the domain on average and the warmer and higher-pressure center locates in the southeast. The distribution of SLP bias induces more pressure dropping over the low-pressure area, enlarging the surface horizontal pressure gradient over the domain. The changes of pressure gradient will further enhance the prevailing south wind and transport more moisture from the sea to the land.Figures8(g)–8(i) show the W10 m bias between forecasts and ERA5 data. A bias of the south wind is detected over the domain, and the mean bias of wind speed is 0.58 m/s, 1.30 m/s, and 1.64 m/s. The south wind bias enhances the moisture transport from the southern seas and induces higher local humidity. As shown in Figures 9(a)–9(c), the forecasted RH2 m is generally higher than that in the ERA5 data. The mean bias is 3.10%, 4.14%, and 4.00%. Higher bias is detected in the southern sea grids than in the northern land grids.The bias of Hwave (Figures9(d)–9(f)) is slightly lower than that in the ERA5 data and decreases roughly from the southwest to the northeast. The mean bias over the domain is −0.12 m, −0.06 m, and −0.01 m. The Hwave here is denoted as the combined height of wind-induced waves with a high frequency and swells with a low frequency. The wind-induced wave height is known to be positively correlated with the wind speed at the sea surface, and the mean bias of wind-induced wave height is 0.34 m, 0.38 m, and 0.42 m for the 24 h, 48 h, and 72 h forecasts, respectively (figure not shown). The overall higher bias of wind-induced wave height indicates that the low Hwave bias results from the offsetting of a higher bias of wind-induced wave height and a lower one of swell height.The ERA5 data has no ocean products; therefore, Figures9(g)–9(i) show the bias distributions of surface seawater salinity between the forecasts and the ocean reanalysis data from the Met Office. The mean bias is 0.31, 0.33, and 0.35. The distributions of salinity bias do not differ in the three forecast tests, and the bias in the southern open seas is obviously weaker than that in the northern offshore regions. The highest bias of simulated salinity appears in the northeastern corner of the domain, which is close to the mouth area of the Pearl River, a river with the second largest discharge in China. In the COAWST model, the freshwater horizontal input is inactive by default, which is probably the primary reason for the bias distributions and requires further improvements in the future.


## 4. Discussion and Conclusion

In this study, an experimental hydrometeorological forecasting system was developed based on the ocean-atmosphere-wave coupled model COAWST, which includes the atmosphere model WRF, the ocean model ROMS, and the wave model SWAN. The system downloads the global forecast data every day from the NCEP website and produces regional forecasts for the next 72 hours, controlled by shell scripts with timers. The system adopts two sets of nesting grids for dynamic downscaling. The inner grids in southern coastal areas of South China, covering Hainan Island, Leizhou Peninsula, and the northwestern corner of the South China Sea, were chosen as the study domain. The 72 h forecast tests were conducted from July 27 to August 31, 2019, and the forecasts for the entire month of August were chosen for model evaluation. These forecasts, produced once a day, were then divided into three groups with different forecast lead times, that is, 24 h, 48 h, and 72 h, for further analysis.To facilitate evaluation, the hourly observations of T2 m, RH2 m, and SLP within the domain were derived from 29 stations. The Taylor diagrams indicated that the forecast system performed well in terms of the CCs and NSDs of the three surface elements. The 24 h forecast series had the highest accuracy with CCs between 0.8 and 0.9 and NSDs between 0.75 and 1.0 in the forecast tests. The 48 h and 72 h forecasts, however, had similar accuracy, with the 48 h forecasts performing slightly better. In the diagrams of T2 m and RH2 m, it can be seen that the 48 h and 72 h forecasts series in the southern stations of Hainan Island performed more poorly in terms of time variability and could be easily discerned from those of other stations owing to their low CCs and NSDs.In addition to the station observations, two sets of reanalysis data, ERA5 and CLDAS, were derived for evaluating forecasted precipitation and T2 m. The spatial distributions of bias had obvious land-sea differences for both elements; that is, more precipitation and lower temperatures were present on land than on sea. For the sea grids, the bias in the northern offshore regions was generally stronger than that in the southern open seas. As for the land grids, the bias in the southern areas of Hainan Island had high spatial heterogeneity, which was obvious from the T2 m bias distributions. In general, the forecasted precipitation was closer to the ERA5 data and the T2 m was closer to the CLDAS data in terms of their mean.To evaluate the time variability of forecasted precipitation and T2 m, the bias series of forecasts minus reanalysis data averaged over the land and sea grids have also been presented. The bias series revealed that the precipitation bias of forecasts minus ERA5 had lower amplitudes than others, but, for the T2 m bias, the series of forecasts minus CLDAS had weaker amplitudes. The extreme values of the bias series roughly corresponded with the strong convective weather over the domain; that is, the bias of precipitation and T2 m was relatively greater during strong convection processes. The forecasted T2 m did not vary as fast as that in the reanalysis data. The forecasting system did not perform well in severe fluctuations of T2 m over the domain, especially in the sea areas.Statistical indices, including CC, RMSE, and STD, were employed to evaluate the statistics of forecasted precipitation and T2 m. The domain-averaged indices, as well as the ones averaged over land and sea, were listed in tables. The indices derived from forecasts and the two sets of reanalysis data did not differ greatly. In general, the indices from forecasts and ERA5 were slightly better than those from forecasts and CLDAS. The indices also exhibited obvious land-sea differences. The CCs and STDs in land grids were, in general, better than those in sea grids, while the RMSEs in land grids were usually poorer. This phenomenon indicated that the forecasted precipitation and T2 m in sea grids performed slightly better in simulating the mean values but performed worse in simulating their fluctuation characteristics. Additionally, the performance level of forecasted statistical indices decreased with the forecast lead time, which was in accordance with the mean bias.In addition to the evaluations of precipitation and T2 m, the bias distributions of other surface elements have also been presented, including SST, SLP, W10 m, RH2 m, Hwave, and surface seawater salinity. The reanalysis data for comparison was still ERA5 data, except for the salinity, for which the ocean product released by Met Office was used. Higher SST bias was detected, especially in northern offshore areas, which was similar to the distribution pattern of T2 m bias in sea grids. The land-sea thermal differences also induced changes of surface pressure. Lower SLP bias was detected over the domain and the lower bias decreased most in the northwestern areas of the domain, which were a center of low pressure. The distribution pattern of SLP bias enhanced the horizontal pressure gradient and induced a higher bias in wind speed at the surface. The W10 m flowing from southwest to northeast was strengthened with more moisture transported to land areas, increasing RH2 m over the domain.The mean Hwave bias of forecasts minus ERA5 was rather weak, varying from −0.12 m to −0.01 m for different lead times. This weak bias resulted from the offsetting of the higher bias of wind-induced waves and lower bias of swells. For the salinity bias at 0 m depth between forecasts and reanalysis data from the Met Office, a systematic difference seemed to exist, and there were no obvious differences among the bias distributions for 24 h, 48 h, and 72 h forecasts. Stronger bias was detected in the offshore regions near coastal lines because the freshwater input was inactive in this study.Overall, the proposed forecast system provided a high accuracy for simulating the surface meteorological and hydrological elements. The COAWST model had a high potential for operational NWP.It is important to note that there were still some uncertainties or limitations in the forecasting results of this study, including the reliability of reference data, model structural deficiencies, and issues with the experimental design.With regard to the reference data, the hourly precipitation data of station observations was unavailable. This lack of station precipitation data led to the absence of a basis of real data in subsequent analyses. In the future, more observational data should be obtained for sea areas, such as buoy observations, to make up for the lack of marine data. Although the ERA5 and CLDAS reanalysis data had high accuracy, there was still some bias between the two datasets and station observations. For example, in Section3.1, the mean T2 m biases of forecasts in the corresponding grids of 29 stations were –0.65 K, 0.37 K, and 0.39 K for 24 h, 48 h, and 72 h forecasts, respectively. In comparison, the mean bias of ERA5 minus station observations was 0.47 K, while the bias of CLDAS was 0.56 K. Thus, the bias in Sections 3.2 and 3.3, which was the difference between forecasts and reanalysis data, differed in accuracy from the bias in Section 3.1. According to the assessment report, the salinity reanalysis data used for comparison in Figures 9(g)–9(i) demonstrated some issues regarding accuracy, because their system adopted the surface and subsurface runoff generated from a land-surface model as the freshwater input of the ocean model, which introduces some uncertainties [60]. In the future, more precise data should be obtained for further evaluation.The deficiencies of the model proposed herein include the absence of some functional modules and the defects of characterization schemes in COAWST. The experimental forecast system developed in this study lacked certain data assimilation modules for ROMS and SWAN. The default initialization of ROMS may lead to lower accuracy in the beginning periods of ocean simulation, while the procedure of reading SWAN’s restart files continuously could also lead to an accumulation of errors during wave simulation. Numerous studies have informed the ocean and wave-assimilation methods of ROMS and SWAN [61, 62]. In the future, the assimilation modules should be introduced, as done in these previous studies, to constrain the simulation of the ocean and waves. In addition, as shown in Figures 9(g)–9(i), the freshwater input has not been factored in and the salinity simulation needs further improvements, too.There were also some limitations regarding the characterization schemes in the COAWST model. Firstly, the atmosphere-wave coupling in the COAWST model only considered momentum-related factors, without considering thermal processes. However, the wave-dissipative heating and sea-spray processes could increase the sensible and latent heat fluxes at the sea surface, which should be included in future studies. Secondly, as mentioned in Section3.2.2, the COAWST model still requires improvements when simulating in extreme events, such as TCs. The Hurricane WRF model (code available online at http://www.dtcenter.org/community−code/hurricane−wrf−hwrf/download), for instance, provides an excellent example of cyclone simulation by introducing vortex-related schemes, forcing data correction, and utilizing characterization optimization. The COAWST model should optimize its dynamic core in WRF as the Hurricane WRF model did, in order to promote the simulation level of TCs.For the experimental design in this study, the COAWST model was configured with its default parameters, without optimization. Some new scheme choices and parameter adjustments were required for further optimization in the simulations over South China. We noted that the forecast simulations in this study are not sufficiently generalizable, because they only lasted for one month within a small area. More simulations, conducted over a longer period and in a larger domain, should be conducted for more detailed evaluations. Further, in addition to surface elements, atmospheric and ocean elements—at different heights and depths—should be addressed in future studies.



---

*Source: 6658722-2021-01-19.xml*