Improving WRF Typhoon Precipitation and Intensity Simulation Using a Surrogate-Based Automatic Parameter Optimization Method

Abstract

Typhoon precipitation and intensity forecasting plays an important role in disaster prevention and mitigation in the typhoon landfall area. However, the issue of improving forecast accuracy is very challenging. In this study, the Weather Research and Forecasting (WRF) model typhoon simulations on precipitation and central 10-m maximum wind speed (10-m wind) were improved using a systematic parameter optimization framework consisting of parameter screening and adaptive surrogate modeling-based optimization (ASMO) for screening sensitive parameters. Six of the 25 adjustable parameters from seven physics components of the WRF model were screened by the Multivariate Adaptive Regression Spline (MARS) parameter sensitivity analysis tool. Then the six parameters were optimized using the ASMO method, and after 178 runs, the 6-hourly precipitation, and 10-m wind simulations were finally improved by 6.83% and 13.64% respectively. The most significant improvements usually occurred with the maximum precipitation or the highest wind speed. Additional typhoon events from other years were simulated to validate that the WRF optimal parameters were reasonable. The results demonstrated that the improvements in 6-hourly precipitation and 10-m wind were 4.78% and 8.54% respectively. Overall, the ASMO optimization method is an effective and highly efficient way to improve typhoon precipitation and intensity simulation using a numerical weather prediction model.

1. Introduction

A typhoon is a strong cyclonic system stemming from the tropical or subtropical Western Pacific with a low pressure center and strong central wind velocities [1,2]. Typhoons bring rainstorms, strong winds, and storm surges, resulting in huge natural damage to the landfall area. Therefore, accurate forecasting of typhoon precipitation and intensity (referring to central maximum wind speed) is crucial to disaster prevention and mitigation for the typhoon landfall area.

A typhoon track is mainly affected by large-scale weather systems such as mid-latitude westerly circulations and subtropical highs. In addition, typhoon size and underlying surface conditions influence the typhoon track. Many studies have tried to obtain the best track simulation by improving the model physical description (e.g., seeking a set of representative physical schemes). However, compared with the improvement of track simulation, the ability to improve typhoon precipitation and intensity simulations is very limited [3,4]. Therefore, this paper pays more attention to improving typhoon precipitation and intensity simulations.

Advanced observation tools (especially weather satellites and radar) and powerful supercomputing technology have promoted the development of mesoscale numerical weather prediction (NWP) models. As a representative of this model class, the Weather Research and Forecasting (WRF) model [5] has incorporated numerous alternative physical parameterization schemes to describe each physical process, and schemes from different physical processes can be arbitrarily combined to form a new model [6]. Taking the WRF model as an example, proposals for how to improve typhoon simulation in the existing literature are briefly given below.

The uncertainty of the description of the initial generated typhoon state, including location, structure, and intensity, has a significant impact on the development and extinction of the typhoon in later stages. Kurihara et al. [7] replaced the poorly resolved tropical cyclone included in the analysis field using a bogus vortex constructed by the Geophysical Fluid Dynamics Laboratory method. However, incompatibility usually occurred between the bogus vortex and the model as the typhoon developed. To avoid such incompatibility, bogus data assimilation was used [8,9]. The strategy was to match the bogus vortex with the real one in structure and intensity. More recently, a dynamic initialization method has been developed and widely used. It initializes the axisymmetric component of the initial vortex of the model close to observed values by an iterative run of the model within a certain time window (e.g., a 6-h interval) while keeping the large-scale environmental and asymmetric components unchanged [10,11].

The uncertainty of the WRF physics description was also analyzed to quantify the best parameterization scheme for each physics component to improve the typhoon simulation. Li et al. [12] assessed the impact of the different schemes for cumulus physics on WRF simulation of typhoon precipitation. Some studies evaluated the effects of scheme combinations from microphysics and the planetary boundary layer on track forecasts [13,14]. More studies focused on assessing scheme combinations from the three physics components for cumulus, microphysics, and the planetary boundary layer [15,16,17]. In addition, more physics schemes were combined to obtain more accurate typhoon forecasts. For example, Srinivas et al. [18] added three land-surface schemes to the previous combinations of three physics components. More recently, Di et al. [19] proposed a Tukey-based combinatorial optimization method to select the optimal scheme combination from all available schemes for microphysics, cumulus, and the planetary boundary layer to perform WRF typhoon intensity simulation. This approach was different from previous studies that used a limited list of fixed scheme combinations to select the best one.

In recent years, highly efficient parameter optimization methods have been proposed to improve WRF simulation. For instance, the multiple very fast simulated annealing (MVFSA) method proposed by Jackson et al. [20] has been used to adjust the five parameters of the Kain-Fritsch cumulus scheme [21]. The results showed that the WRF model with optimal parameters reduced over-prediction of the precipitation amount by the model with default parameters. Santanello et al. [22] used the NASA Land Information System (LIS-OPT) to calibrate the parameters of the Noah land-surface module in the WRF model to produce the best estimates of land-surface fluxes in land-atmosphere interaction. Duan et al. [23] conducted parameter optimization on seven physics schemes from the WRF model by combining a parameter screening method and an adaptive surrogate model-based optimization (ASMO) method. The results showed that the daily precipitation simulation was improved by 18% with only about 200 WRF runs using the ASMO method. Next, Di et al. [24] assessed the applicability of the WRF optimal parameters under different precipitation simulation conditions and finally confirmed that the optimal parameters are reliable. However, it is not clear whether parameter optimization for the WRF model can also improve simulations of typhoon precipitation and intensity.

This paper tried to address three problems related to WRF parameter optimization to improve simulations of typhoon precipitation and intensity. The first was the question of which parameters were sensitive to simulations of typhoon precipitation (mm day⁻¹) and central 10-m maximum wind speed (abbreviated as 10-m wind (m s⁻¹) in this paper). The second was how to optimize the screened sensitive parameters to improve the precipitation and 10-m wind simulations. The third was whether the optimal parameters obtained were robust when typhoon events in other years were simulated.

This paper was organized as follows. Section 2 presents the data used and the methodology related to systematic parameter optimization for the complex dynamic model. Section 3 presents the typhoon optimization results and validation analyses. Section 4 discusses the difference between single-objective and multi-objective optimization. Conclusions are presented in the last section.

2. Data and Methodology

2.1. Observed Data

This study focuses on parameter optimization of the WRF model to improve the precipitation and 10-m wind simulations. Precipitation observation data came from Chinese hourly merged precipitation analysis products (CHMPA-Hourly, version 1.0) [25], available at 0.1° × 0.1° horizontal resolution, which combines 30,000 gauge stations in China, and the Climate Precipitation Center’s (CPC) morphing satellite technique precipitation product (CMORPH). The 10-m wind observation data came from the China Meteorological Administration (CMA)-Shanghai Typhoon Institute (STI) best-track dataset for tropical cyclones over the western North Pacific. The high-altitude validation data were the 500-hPa wind field and the 500-hPa potential height field, with a horizontal resolution of 0.125°, from the ERA-interim reanalysis dataset of the European Centre for Medium-Range Weather Forecasts (ECMWF) [26].

2.2. WRF Model Configuration for Typhoon Simulations

WRF model version 3.7.1 was used to simulate six 3-day typhoon events over South China. The study area was a two-grid horizontally nested domain (see Figure 1). The outer domain (i.e., area d01 in Figure 1) was composed of 210 × 180 grid cells with a spatial resolution of 18 km. The inner domain was composed of 210 × 210 grid cells with a spatial resolution of 6 km. To save computing time, two inner domains were set up according to the locations of the observed typhoon tracks. The first inner domain was d02A, mainly covering Guangdong, Guangxi, and Hainan Provinces in China, and the second was d02B, mainly covering Jiangxi and Fujian Provinces in China. The uniform time step was 90s. The lateral boundary data for the d01 domain were obtained from the NCEP FNL (Final) Operational Global Analysis dataset, available at 1° × 1° horizontal resolution and 6 h interval.

To obtain reasonable typhoon simulations, a set of fixed WRF physical parameterization schemes was used, following the operational setup of CMA typhoon forecast system (communicated with CMA experts). In addition, the fixed WRF parameterization schemes were used to simulate typhoon track and intensity in some literature [11,27,28]. The specific schemes were the revised Monin-Obukhow surface layer scheme [29], the Kain-Fritsch (new Eta) cumulus scheme [30], the WRF single moment six-class microphysics scheme (WSM6) [31], the Dudhia shortwave scheme [32], the rapid radiative transfer mode (RRTM) longwave scheme [33], the Noah land-surface model [34], and the Yonsei University (YSU) planetary boundary layer scheme [35].

Twenty-five adjustable parameters were identified based on the fixed physics parameterization schemes. These parameters are thought to be related to simulations of precipitation and wind according to previous studies [23,24,36,37]. Table S1 lists their ranges and physical meanings. Note that the parameter list may very well be incomplete, although the best possible effort was made to find the adjustable physics parameters associated with precipitation and wind simulations.

In this study, the parameter optimization focused mainly on the 6-hourly precipitation and 10-m wind simulations of six 3-day typhoon events from 2013 to 2015. Detailed descriptions of the typhoon tracks and durations are shown in Figure 1 and

Table 1. Selected typhoon events from 2013 to 2015.

No.	Typhoon Events	Simulation Period	Typhoon Period
(1)	201306 Rumbia	2013-07-01-08:00–2013-07-04-08:00	2013-07-01-08:00–2013-07-02-20:00
(2)	201409 Rammasun	2014-07-17-14:00–2014-07-20-14:00	2014-07-17-14:00–2014-07-19-20:00
(3)	201522 Mujigae	2015-10-03-08:00–2015-10-06-08:00	2015-10-03-08:00–2015-10-05-08:00
(4)	201307 Soulik	2013-07-12-20:00–2013-07-15-20:00	2013-07-12-20:00–2013-07-14-02:00
(5)	201410 Matmo	2014-07-22-08:00–2014-07-25-08:00	2014-07-22-08:00–2014-07-24-14:00
(6)	201513 Soudelor	2015-08-08-02:00–2015-08-11-02:00	2015-08-08-02:00–2015-08-10-08:00

respectively. In addition, six new typhoon events from 2016 and 2017 were simulated to validate the robustness of the WRF optimal parameters for improving typhoon simulations.

2.3. Systematic Parameter Optimization Framework

The WRF model describes a number of physical processes involving the atmosphere and the land, which contain numerous adjustable parameters and entail time-consuming computation. Therefore, it is not feasible to optimize the WRF parameters using traditional parameter optimization methods. Here, a systematic parameter optimization framework was proposed to conduct parameter optimization for the WRF model. The framework included a sensitive analysis method involving multivariate adaptive regression splines (MARS) and a highly efficient parameter optimization method involving adaptive surrogate modeling-based optimization (ASMO).

2.3.1. MARS Sensitivity Analysis Method

The MARS method [38] is a sensitivity analysis method based on regression. It requires that a reasonable MARS regression model be first built based on the parameter samples. The precise specification is as follows. The MARS model is a polynomial regression function built by a set of basis functions, which exist in three forms: constant, hinge function, and the product of several hinge functions. Constructing a reasonable MARS model requires two steps: the forward step constructs an overfitted model using the paired basis functions on each parameter sample, and the backward step prunes the overfitted model by deleting the least effective term. Once the reasonable MARS model has been built, a generalized cross-validation (GCV) function is used to evaluate the MARS regression model.

$$GCV(M)=\frac{1}{N}\frac{{\displaystyle \sum _{i=1}^N{(Y_i-\widehat{Y_i})}^2}}{{\left[1-\frac{1+c(M)d}{N}\right]}^2}$$

(1)

where N is the number of parameter samples, $Y_i$ is the output of the i-th parameter sample, $\widehat{Y_i}$ is the MARS simulation value for the i-th parameter sample, d is the effective number of degrees of freedom, and $c(M)$ is a positive penalty factor of model M for adding a low-order function. The lower the GCV, the closer the MARS regression model M is to the real model.

Parameter sensitivity was then evaluated based on the reasonable MARS model and its evaluation formula. The precise specification is as follows. For each parameter, if all terms related to that parameter are removed from the reasonable MARS model $\hat{Y}$, a new GCV value results. The parameter sensitivity relies on the difference between the new and original GCV values. The larger the difference, the more sensitive is the parameter. To normalize the parameter sensitivity score, the GCV difference corresponding to each parameter is divided by the maximum GCV difference. The result is that the score is one for the most sensitive parameter and zero for the most insensitive parameter. After sensitivity analysis for adjustable parameters, the sensitive parameters for model output are screened.

2.3.2. ASMO Parameter Optimization Method

To reduce the dimension of the optimized parameters, only the screened sensitive parameters were optimized by the ASMO method. With the help of a statistical regression model, the ASMO method could speed up its search for the optimal parameters (i.e., the optimal parameter values) of a complex physical model [39]. The main optimization procedure of the ASMO method can be briefly described as follows:

• Sample to the adjustable parameter ranges using a uniform sampling method. Then these parameter samples are respectively put into the physical model instead of the default parameter values to obtain the corresponding model outputs (or the output errors compared with observed data). The parameter samples and the corresponding model outputs constitute an initial sample set;
• Build a statistical regression model in the initial sample set. Then search for the optimal value of the statistical model using a traditional parameter optimization method. The corresponding optimal parameters (i.e., the optimal parameter values) are finally found;
• Put the optimal parameters of the statistical regression model into the physical model to update the model output. A new sample point is generated based on the optimal parameters and the corresponding physical model output. Update the initial sample point set by adding the new sample point;
• Repeat Steps 2 and 3 until the parameter optimization convergence condition for the physical model is met. The globally optimal parameters of the physical model have now been found.

In this study, the uniform model evaluation function for parameter sensitivity analysis and parameter optimization was the root mean square error (RMSE); however, the reference of the model evaluation was different. For sensitivity analysis and optimization, the references were WRF simulation results with default parameters and observations respectively. To reduce computational cost, two inner domains of d02A and d02B were set up in experiment design, so the model evaluation focused on the precipitation and 10-m wind simulations in d02A and d02B. In the ASMO optimization procedure, the quasi-Monte Carlo (QMC) method was selected as the uniform sampling method to produce the uniform sample points for building the reasonable regression model. The optimization search on the statistical regression model was conducted by the shuffled complex evolution (SCE-UA) global optimization method [40]. Here, the convergence condition of ASMO optimization was that the local optimal value of the WRF simulation remains unchanged after 30 searches (equal to five times the dimensionality of the parameters). The overall objective function for ASMO optimization is defined as

$$\mathrm{NRMSE}=\left(\frac{1}{N}{\displaystyle \sum _{i=1}^N\frac{RMS{E}_{p,pre}^i}{RMS{E}_{def,pre}^i}}+\frac{1}{N}{\displaystyle \sum _{i=1}^N\frac{RMS{E}_{p,wind}^i}{RMS{E}_{def,wind}^i}}\right)\times 0.5$$

(2)

where NRMSE is the normalized RMSE for the 6-hourly precipitation and 10-m wind simulations, $RMS{E}_{p,pre}^i$ and $RMS{E}_{p,wind}^i$ are the RMSE of the WRF simulation with perturbed parameters in the i-th event for precipitation and 10-m wind, respectively. Similarly, $RMS{E}_{def,pre}^i$ and $RMS{E}_{def,wind}^i$ are the corresponding WRF simulation errors with the default parameters. N is the total number of optimized typhoon events.

3. Results

3.1. Sensitivity Analysis Results

Using the QMC sampling method, 250 parameter sample sets were taken from the 25 parameter ranges listed in Table S1. The corresponding WRF simulation errors (i.e., RMSE) were also obtained by comparison with the default simulation results. As a pair of input variables, the 250 parameter samples and the corresponding WRF simulation errors were entered into the MARS sensitivity analysis method, and the parameter sensitivity scores were finally obtained.

Figure 2 shows the normalized sensitivity scores for the 25 parameters. The x-axis represents the 25 parameters, and the y-axis represents parameter sensitivity scores for the different regions. Three common sensitivity parameters were found for the precipitation simulations: P5, P6, and P7. In addition, P4 and P23 were sensitive for precipitation simulations in the d02B domain. Finally, five parameters (P4, P5, P6, P7, and P23) were ranked as sensitive parameters for precipitation simulations. For 10-m wind, two common sensitive parameters, P3 and P4, were found. In addition, P6 and P23 were sensitive in the d02A and d02B domains respectively. Overall, four parameters (P3, P4, P6, and P23) were sensitive for the 10-m wind simulations.

To validate the reasonability of the screened sensitive parameters, the responses of model performances to six screened sensitive parameters and four cloud microphysics parameters (i.e., ice_stokes_fac, nor, dimax, and peaut) which were thought to be insensitive using MARS method were compared. The results are shown in Figure 3. Each parameter range was equally divided into 10 bins, and the red curve represented an average of the simulation errors at each bin. The first and second lines represented the responses of precipitation simulation errors to the screened precipitation sensitivity parameters (pd, pe, ph, and pfac corresponding to P5, P6, P7, and P23 respectively in Figure 2) and cloud microphysics parameters respectively. The responses of 10-m wind simulation errors to 10-m wind sensitivity parameters (znt_zf, karman, pe, and pfac corresponding to P3, P4, P6, and P23 respectively in Figure 2) and cloud microphysics parameters were represented in the third and fourth lines of Figure 3 respectively. It is noteworthy that the simulation error was RMSE of WRF simulations with between-default and perturbed parameters.

Figure 3 showed that the perturbation for the sensitive parameters of precipitation and 10-m wind brought more significant model outputs changes than that for the insensitive parameters (cloud microphysics parameters). Compared with sensitive parameters, the variances of the precipitation and 10-m wind simulation errors caused by cloud microphysics parameter perturbation were basically a straight line. Therefore, the cloud microphysics scheme parameters should be insensitive and the variances of their values had no significant influence on the precipitation and 10-m wind simulations. It demonstrated that the parametric sensitivity conclusion (shown in Figure 2) drawn by ASMO method were reasonable. Finally, the six sensitive parameters were quantified and are listed in

Table 2. List of sensitivity parameters for typhoon precipitation and 10-m wind simulations of WRF model version 3.7.1.

Index	Scheme	Parameter	Default	Range	Description
P3	Surface layer (module_sf_sfclayrev.F)	znt_zf	1	0.5–2	Scaling related to surface roughness
P4		karman	0.4	0.35–0.42	Von Kármán constant
P5	Cumulus convection (module_cu_kfeta.F)	pd	1	0.5–2	Scaling related to downdraft mass flux rate
P6		pe	1	0.5–2	Scaling related to entrainment mass flux rate
P7	Planetary boundary layer (module_bl_ysu.F)	ph	150	50–350	Starting height of downdraft above updraft source layer (hPa)
P23		pfac	2	1–3	Profile shape exponent of the momentum diffusivity

.

3.2. Parameter Optimization Results

Figure 4 shows the optimization speed for the total NRMSE of the precipitation and 10-m wind simulations using the ASMO method. The minimum error was about 0.92 in 100 simulations of the initial parameter samples obtained by the QMC sampling method. This meant that the precipitation and 10-m wind simulations were improved on average by 8% compared with the WRF default simulations. Next, the total NRMSE was further reduced to 0.898 by 78 adaptive sampling runs. After this, an additional 30 runs (equal to five times the dimensionality of the parameters) did not change the previous minimum NRMSE. Therefore, the convergence criterion of the optimization was met, and the optimal parameters had finally been found. In a word, the ASMO method improved the precipitation and 10-m wind simulations by 10.2% on average using only 178 model runs. This demonstrated that the ASMO optimization method was highly efficient.

To illustrate the optimization results, comparisons of WRF simulations with default and optimal parameters for 6-hourly precipitation and 10-m wind were conducted. The results are shown in Figure 5. Overall, the simulations for the 6-hourly precipitation and 10-m wind were improved by 6.83% and 13.64% respectively. This demonstrated that the improvement for wind was larger than that for precipitation. For six single events, all precipitation simulations were improved to an extent varying from 0.26% for event (2) to 11.88% for event (6), and four of six 10-m wind simulations were improved to an extent varying from 10.43% for event (4) to 35.04% for event (2). Note also that some event simulations for 10-m wind showed negative improvements (e.g., events (3) and (6)). This is related to the optimization objective function, which focuses on the accumulated simulation errors for all six events rather than on each one. However, note that the negative improvements mainly occurred in events with lower simulation errors. In addition, the larger simulation errors for 10-m wind and precipitation occurred in the d02A and d02B regions respectively. The larger the error, the more significant the improvement was.

Besides comparisons for single-typhoon simulations, WRF typhoon simulations with default and optimal parameters were compared for different lead times. The results are shown in Figure 6. Figure 6a shows that the optimal parameters improved the WRF precipitation simulations at 6-hourly lead time by amounts ranging from 0.07% to 25.20%. Similarly, Figure 6b shows that the optimal parameters improved the 10-m wind simulations at 6-hourly lead time except for the 30th hour (−3.55% improvement), and other positive improvements ranged from 1.54% to 41.22%, which were generally greater than the improvement rate of precipitation simulations at each lead time. In addition, it seemed that the improvement rate had an increasing tendency for both precipitation and 10-m wind simulations as lead time increases.

It was noted in Figure 6b that the difference between 10-m wind default simulation and observation was low (less than 8 m s⁻¹), so the smaller reduction for RMSE by the optimal simulation would demonstrate a larger percentage improvement. To avoid misleading results, it was suggested that the error value and its relative improvement percentage should be combined to accurately evaluate the model performance improvement.

Comparisons of the spatial distribution of 3-day precipitation simulations using the WRF model with default and optimal parameters were also conducted, and the results are shown in Figure 7. The first column represents the spatial distribution of observed daily precipitation for six typhoons in the d02A and d02B regions. It showed that strong precipitation occurred mainly in the northeastern parts of Hainan and Fujian Provinces and the coastal areas of Fujian Province. The second and third columns represent the differences between observations and WRF simulations with default and optimal parameters respectively. By comparing Figure 7b with Figure 7c, it is found that the large positive deviation (marked in red and yellow in Figure 7b) for default precipitation simulations in the d02A region was significantly reduced by the optimal precipitation simulations. A similar situation can also be observed in the d02B region by comparing Figure 7e with Figure 7f. In addition, the improvement was significant for simulation of strong precipitation. For instance, large precipitation amounts (over 60 mm day⁻¹) occurred in the northeastern part of Hainan Province in the d02A region, where deviations were lower for the optimization simulation than for the default simulation. For the d02B region, a similar situation occurred north of the border between Fujian and Jiangxi Provinces.

Comparisons of 10-m wind evolution simulations for each single event were also conducted, and the results are shown in Figure 8. Noted that the 10-m wind was an abbreviation of typhoon central 10-m maximum wind speed (referred in introduction section), so it represented a specific grid point value near typhoon center. The location of 10-m wind in Figure 8 was constantly shifting as the typhoon center moved. Overall, compared with the default simulations, the optimal simulations were closer to observations for most of the lead times. The most significant improvements occurred in events (1) and (2), and therefore the average improvement of the optimal simulations in the d02A region including events (1), (2), and (3) was higher than in the d02B region including events (4), (5), and (6), a result that is consistent with the conclusions shown in Figure 5. The optimal simulations decreased the default simulation values to approximate the observations in event (1), while increased the default simulation values close to observations in event (2). Therefore, the parameter optimization method is different from statistical deviation correction, which provides an overall upward or downward shift for the default simulation results to approximate observations. Hence, the parameter optimization method seemed to be more reasonable than the deviation correction method. In addition, there was some indication that the improvement was more significant as the lead time increased.

The results in Figure 8 seemed that default and optimal experiments were almost identical except event (2), which was mainly caused by the definition of the objective function. According to formula (2), the optimization method was more inclined to improve the total simulation error rather than the error of a single event, so more improvement was achieved with events with larger simulation errors. It was apparent that the event (2), with the largest simulation errors, had the most significant improvement. In addition, the optimal simulations did not deviate from the default simulations too far or even improved them for other events with relatively low simulation errors. Therefore, the optimal parameter was reasonable. Future multi-objective parameter optimization could balance the improvement percentage of all individual events.

3.3. Atmospheric Structure Analysis

Two representative events (events (2) and (5)) were selected to analyze the differences between WRF simulations with default and optimal parameters on high-altitude variables such as 500-hPa geopotential height, 500-hPa wind, and precipitable water.

Figure 9 shows the simulated fields of 500-hPa geopotential height and wind at lead times of 24 h, 48 h, and 72 h using the WRF model with default and optimal parameters compared with ERA-interim observations for event (2). The first line represents the observations. The second and third lines represent the simulated field at 500 hPa using the WRF model with default and optimal parameters respectively. The red lines marked with 588 dagpm represented the scope of the subtropical high. By comparing the location of 588 dagpm lines, it was found that when the optimal simulation was used, the southern scope line of the subtropical high moved slightly toward the observation at 24 h; however, the whole scope of the subtropical high significantly moved northeast closer to the observation at 48 h and 72 h. Overall, the simulation of the scope of the subtropical high (marked by the red line) showed improvement, especially for 48 h and 72 h lead times. When the scope of the subtropical high moved east, the flow on the west side of the subtropical high varied from north to northwest, making the typhoon center move slightly south, closer to observations (Figure 9b,e,h). Accordingly, the simulation of typhoon central wind speed was strengthened (see Figure 9d,g, or Figure 9e,h), which is consistent with the results in Figure 8b.

The same comparisons were also conducted for event (5) and the results are shown in Figure 10. As the conclusion of Figure 9, the simulation for the scope of the subtropical high for event (5) showed improvement, especially for lead times of 48 h and 72 h. The difference was that the optimization expanded the scope of the subtropical high in event (5) and shrank the scope in event (2). The improved simulation at 48 h caused the flow on the east and west sides of the subtropical high to shift from northwest to north, which was close to the observed flow. In the same way, the improvement of the subtropical high simulation at 72 h caused the flow on the east and west sides of the subtropical high to shift from northwest to northeast, which was close to the observed flow.

The precipitable water is the total amount of water vapor in a vertical column with unit cross-sectional area. It is an important indicator of precipitation. Figure 11 and Figure 12 show a comparison of the results of precipitable water simulations with default and optimal parameters for event (2) and event (5) respectively. To better demonstrate the improvement in the precipitable water simulation by the WRF model with optimal parameters, the areas with high water-vapor content were shaded. The shaded areas in Figure 11 and Figure 12 represent values higher than 65 kg m⁻² and 70 kg m⁻² respectively. As shown in Figure 11 and Figure 12, the optimal simulations reduced the shaded area of the default simulations to be closer to observations, especially for the area around the typhoon center.

3.4. Validation Analysis of WRF Optimal Parameters

It is necessary to address the question of whether the optimal parameters will still work on WRF simulations of new typhoon events. Therefore, six new typhoon events from 2016 and 2017 were simulated to verify whether the WRF optimal parameters were still reasonable. The new events were different from the previous optimization events from 2013 to 2015. Other than the simulation events, the validation experiment used the same configuration as the previous optimization experiment, including simulation area, WRF model configuration, and forcing data source. In addition, the six new validation events were equally divided into two d02 domains. The uniform simulation period was three days. Detailed descriptions of the observed typhoon tracks and durations are presented in Figure 13 and

Table 3. Typhoon validation events from 2016 and 2017.

No.	Typhoon Events	Simulation Period	Typhoon Period
(I)	201604 Nida	2016-08-01-12:00–2016-08-04-12:00	2016-08-01-12:00–2016-08-03-00:00
(II)	201713 Hato	2017-08-22-18:00–2017-08-25-18:00	2017-08-22-18:00–2017-08-25-00:00
(III)	201714 Pakhar	2017-08-26-18:00–2017-08-29-18:00	2017-08-26-18:00–2017-08-28-00:00
(IV)	201601Nepartak	2016-07-07-06:00–2016-07-10-06:00	2016-07-07-06:00–2016-07-10-00:00
(V)	201617 Megi	2016-09-26-18:00–2016-09-29-18:00	2016-09-26-18:00–2016-09-29-18:00
(VI)	201709 Nesat	2017-07-28-12:00–2017-07-31-12:00	2017-07-28-12:00–2017-07-31-00:00

respectively.

Comparisons of precipitation (10-m wind) simulations for the six new typhoon events (I)–(VI) using the WRF model with default and optimal parameters were conducted, and the results are shown in Figure 14. Compared with the default WRF simulations, the WRF model with the optimal parameters improved the simulations of 6-hourly precipitation and 10-m wind by 4.78% and 8.54% respectively. It is apparent that the optimal parameters not only improved the typhoon simulations of precipitation and 10-m wind for the optimization period from 2013 to 2015, but also for the validation period from 2016 and 2017. Similarly, the improvement rate for wind simulation was higher than for precipitation simulation. These conclusions on the comparison of precipitation simulations with validation events are consistent with those for previous optimization events, which in turn confirms that the optimal parameters obtained by the ASMO method are robust.

Overall, the optimal simulations reduced the positive deviation of the default simulations for the spatial distribution of precipitation, and the improvement was significant for the simulation of strong precipitation (Figures not shown). Comparative analyses were also performed for the 10-m wind simulations in each validation event. Figure 15 shows the corresponding results. Generally, compared with the default simulations, the optimization simulations were closer to observations for most lead times. Among these, the most significant improvements occurred in events (II) and (V). Note also the slightly negative improvements for events (I) and (III), with lower errors between the default simulations and observations. This happened because the objective function for optimization was the total accumulated error for the simulations of all six events, and therefore some event simulations with lower errors might have been slightly sacrificed. However, the optimal parameters overall improved the 10-m wind simulations by the default WRF model, as shown in Figure 14b. It was also found that the improvement of the optimal simulations for 10-m wind was more significant after 24 h, implying that the effect of parameter optimization may be highlighted when the effect of initialization has weakened.

3.5. Parametric Comparison and Physical Interpretation

The default and optimal parameters were normalized based on the parameter ranges listed in Table 2. The normalized results of the default and optimal parameters are shown in Figure 16. The variations for all the optimal parameters are inconsistent. For the parameters znt_zf (scaling related to surface roughness), karman (von Kármán constant), and pe (scaling related to entrainment mass flux rate), their optimal values decreased compared with their default values. However, the opposite situation occurred for the parameters pd (scaling related to downdraft mass flux rate), ph (starting height of downdraft above updraft source layer, and pfac (profile shape exponent of the momentum diffusivity coefficient).

To further demonstrate the relation between model outputs and the sensitive parameters, the analyses on how the variances of the sensitive parameters affect the model performances were also conducted. Based on observations and the previous 178 model simulations from optimization experiment, the response of three variables (i.e., precipitation, 10-m wind, and track) to six sensitive parameters were obtained, and the results are shown in Figure 17. Here, the range of each parameter was equally divided into 10 bins, and the red curves represented an average of the simulation errors at each bin. The simulation error with default parameters was marked as a red cross. The first and second lines demonstrated the variances of the precipitation and 10-m wind simulation errors as parameter value increased respectively. The variances of the track simulation errors were represented in the third lines.

For znt_zf parameter, the default value was close to the optimal value for the precipitation and 10-m wind simulations, but it was still be found that slightly decreasing the default value would improve the precipitation and 10-m wind simulations. The znt_zf was insensitive to track simulation, and therefore it was difficult to improve the track simulation by adjusting znt_zf value. For karman parameter, the default value in the precipitation and 10-m wind simulations was significantly higher than the optimal value, and an apparent upward trend existed between them. Therefore, decreasing the default value of karman could improve the WRF simulations on precipitation and 10-m wind. Like znt_zf, karman was insensitive to track simulation. For pd parameter, the default value was lower than the optimal value for the simulations of precipitation, 10-m wind, and track. Therefore, increasing the default value would improve the respective default simulations.

For pe parameter, the default value was larger than the optimal value for three variables, and the significant upward trend existed between the input parameter and model output errors, so the simulations for three variables could be improved by decreasing the default value. For ph parameter, there was an obvious downward trend for the 10-m wind and track simulation errors as the parameter value increased; however, a slightly upward trend existed for precipitation simulation errors. Therefore, increasing the default value would improve the 10-m wind and track simulations, and the opposite situation occurred in precipitation simulation. For pfac parameter, the default value was lower than the optimal value for three variables, and the downward trend existed between the input parameter and model output errors. Therefore, increasing the default value will improve three variable simulations.

Previous analyses have found (e.g., Figure 7 and Figure 8) that overall, the default simulations overestimated (underestimated) the amount of precipitation (10-m wind), whereas the optimal simulations partially rectified the overestimation (underestimation). This can be interpreted using the physical meanings of the parameter perturbations. Lower values of znt_zf mean that a decreased roughness length exists in the surface layer, and therefore the simulated wind speed will be increased. The 10-m wind simulation is only sensitive to znt_zf parameter (i.e., P3), and therefore the optimal znt_zf values work mainly on wind simulations. A lower karman value can reduce the momentum exchange coefficient (drag coefficient) at the near-surface layer, and therefore the upper 10-m wind speed increases. For precipitation simulation, a smaller karman value meant a smaller enthalpy exchange coefficient at the near-surface level, which caused less water vapor to rise into the atmosphere, leading to less precipitation.

The pd, pe, and ph are the parameters from the cumulus convection scheme, and therefore it is reasonable that precipitation is sensitive to them, as shown in Figure 2. A larger downdraft mass flux rate (larger pd) leads to more condensed water evaporation and further to less precipitation. A higher downdraft-flux starting height (larger ph) will bring more downdraft flow, leading to more condensed water evaporation and less precipitation. A lower entrainment rate (smaller pe) will reduce the mixed amount of the updraft from cold environmental air, leading to less condensed water production and therefore less precipitation. With lower pe, less cold air is involved in the updraft, prompting further development of convection. Correspondingly, the wind speed is increased.

The pfac has a positive effect on the momentum diffusion intensity of turbulent eddies in the planetary boundary layer. When pfac increases, the vertical momentum and heat diffusion intensity are strengthened, inducing higher wind speed. However, the increase in heat diffusion intensity makes more water vapor from the ground transfer upward, leading to more precipitation. It is noteworthy that the increase of pfac parameter improved the low bias of 10-m wind simulation, but had a negative improvement to precipitation simulation. However, the total objective error combined precipitation and 10-m wind evaluations was finally reduced as shown in Figure 4. It can be inferred that the 10-m wind adjustment is more discriminating in terms of pfac parameters than precipitation, and pfac has an antagonistic effect on precipitation simulation compared with pd and ph.

4. Discussion

This study has focused on improving typhoon precipitation and 10-m wind simulations for the WRF model using the ASMO parameter optimization method. Equation (2) indicates that equal weights were assigned to the simulation errors of precipitation and 10-m wind in the overall objective function. Finally, it has been demonstrated that both precipitation and 10-m wind simulations have been improved. However, it remains unresolved whether optimizing only one variable (e.g., precipitation) will reduce the simulation accuracy of the other variable (e.g., 10-m wind).

To address this question, two single-objective optimizations for precipitation and 10-m wind speed simulations were respectively conducted by the ASMO method. Here, the experimental setup was the same as in the previous optimization experiment.

Table 4. Comparisons of the variable improvements for the two single objective optimization and the one dual objective optimization.

Name	Precipitation Improved	10-m Wind Improved
Run 1	8.5%	6.5%
Run 2	6.1%	14.5%
Run 3	6.8%	13.6%

shows the improvement rates of precipitation and 10-m wind simulations for three optimization experiments, including one optimizing only precipitation (defined as Run 1), a second optimizing only 10-m wind speed (defined as Run 2), and a third simultaneously optimizing precipitation and 10-m wind speed (defined as Run 3). Compared with the default simulations, the improvements in precipitation simulation accuracy for Run 1, Run 2, and Run 3 were 8.5%, 6.1%, and 6.8% respectively. Similarly, the improvements in 10-m wind speed simulation accuracy for Run 1, Run 2, and Run 3 were 6.5%, 14.5%, and 13.6% respectively. Based on these results, two obvious conclusions can be drawn. The first was that single-objective optimization improved the simulation not only of the optimized variable, but also of the other variable. The second was that the improvement achieved through single objective-optimization was greater than that achieved by multi-objective optimization for an optimized variable. Note also that the first conclusion may not be correct if two opposite variables are simultaneously optimized.

5. Conclusions

This study applied a systematic parameter optimization framework to a complex mesoscale WRF model to improve typhoon simulations of precipitation and central 10-m maximum wind speed (10-m wind). It first used the MARS sensitivity analysis method to screen out six sensitive parameters for the WRF precipitation and 10-m wind speed simulations from 25 adjustable parameters of seven physical parameterization schemes. Then the ASMO method was used to optimize the six sensitive parameters to obtain optimal WRF simulations of typhoon precipitation and 10-m wind over South China. By comparing these results with the WRF default simulation, the optimization results were evaluated from six aspects, including six single-event simulations, simulations at 6-hourly lead time, precipitation spatial distribution, 10-m wind evolution, 500-hPa geopotential height and wind field, and precipitable water. Finally, the optimal parameters obtained were validated by comparing the simulations of new typhoon events from other years and discussing the physical meanings of the parameter perturbations.

WRF simulation is very time-consuming due to its high spatiotemporal resolution and detailed physics description, and therefore the thousands of model runs required by traditional optimization methods are difficult to achieve. Now, using only 178 runs, WRF optimal parameters for typhoon precipitation and 10-m wind simulations with a resolution of 6 km have been found by the ASMO parameter optimization method. The 178 runs consisted of 100 initial parameter perturbation runs and 78 adaptive parameter search runs. Obviously, the ASMO approach greatly reduces the number of model runs for WRF parameter optimization and is therefore very highly efficient.

The comparisons of WRF typhoon simulations with default and optimal parameters were divided into two phases: one covered the optimization period from 2013 to 2015 and the other the validation period from 2016 and 2017. During the optimization period, the optimal parameters improved overall the simulations of 6-hourly precipitation and 10-m wind by 6.83 % and 13.64 % respectively. In addition, the precipitation and 10-m wind simulations were improved at 6-hourly lead time. By comparing the results with the spatial distribution of observed precipitation, it was found that the positive deviation of the default precipitation simulation was significantly reduced by the optimal simulation and that the improvement was also significant in strong precipitation areas like the northeastern parts of Hainan Province and the coastal areas of Fujian Province. Overall, the average improvement for the 10-m wind simulations in the d02A domain was greater than in the d02B domain. By comparing the simulations of 500-hPa geopotential height and wind field, the simulated scope of the subtropical high showed significant improvement, especially for lead times of 48 h and 72 h, which made the wind speed of the typhoon center stronger and more realistic. The location of the typhoon center was also slightly moved closer to the observed center. For the simulations of precipitable water, the optimal results were closer to observations than the default results, especially for the region with water vapor greater than 65 kg m⁻². For six typhoon simulations during the validation period, the WRF model with optimal parameters demonstrated better simulation capability for precipitation and 10-m wind than the WRF model with default parameters. The improvements in 6-hourly precipitation and 10-m wind simulations were 4.78 % and 8.54 respectively. This demonstrated that the optimal parameters obtained by the ASMO method also worked for improving the simulations of new typhoon events. In addition, the physical meanings of parameter perturbations were analyzed to validate the reasonableness of the optimal parameters. Therefore, the optimal parameters were thought to be robust, and the ASMO method was effective.

During the optimization and validation periods, the overall improvements in the typhoon precipitation and 10-m wind simulations were demonstrated for the WRF model with the optimal parameters obtained by the ASMO method. However, it is worth noting that there were several events with negative improvement, which was related to the design of the objective function. Here, equal weights were not only assigned to the evaluation of the total precipitation and 10-m wind simulations, but also among the different events for one variable, and then summed to make up the overall objective function. In addition, the different weight allocations between precipitation and 10-m wind evaluation in the overall objective function affected the simulation accuracy of individual variables, as discussed in the fourth section. Therefore, future work should focus on real multi-objective optimization, which searches the optimal parameter set for different weights for multi-objective functions using generalized multi-objective optimization methods such as NSGA-II [41] and MO-ASMO [42].