A micro-genetic algorithm (GA v1.7.1a) for combinatorial optimization of physics parameterizations in the Weather Research and Forecasting model (v4.0.3) for quantitative precipitation forecast in Korea

One of the biggest uncertainties in numerical weather predictions (NWPs) comes from treating the subgrid-scale physical processes. For more accurate regional weather and climate prediction by improving physics parameterizations, it is important to optimize a combination of physics schemes and unknown parameters in NWP models. We have developed an interface system between a micro-genetic algorithm (µ-GA) and the WRF model for the combinatorial optimization of cumulus (CU), microphysics (MP), and planetary boundary layer (PBL) schemes in terms of quantitative precipitation forecast for heavy rainfall events in Korea. The µ-GA successfully improved simulated precipitation despite the nonlinear relationship among the physics schemes. During the evolution process, MP schemes control grid-resolving-scale precipitation, while CU and PBL schemes determine subgrid-scale precipitation. This study demonstrates that the combinatorial optimization of physics schemes in the WRF model is one possible solution to enhance the forecast skill of precipitation.

1 Introduction

For numerical weather forecasting to be accurate, a numerical model should be able to represent real atmospheric conditions in terms of dynamics (i.e., governing equations), physics (i.e., parameterizations), and numerics (e.g., resolution and coordinate system). More accurate initial conditions should also be provided. One of the biggest uncertainties in numerical weather predictions (NWPs) comes from treating the subgrid-scale physical processes that have not been sufficiently understood. The subgrid-scale physical processes are parameterized in NWP models through empirical evidence, such as the derived value from observations and/or theoretical backgrounds. Therefore, the accuracy of physics parameterizations strongly depends on the following: the value of parameters in given uncertainty ranges in parameterization schemes and the choice of parameterization schemes for each corresponding physical process. Note that prior to model simulation the unknown parameters and schemes should be fitted to the regional weather and climate to reduce considerable uncertainties in models.

NWP models have several categories of subgrid-scale physical processes – shortwave and longwave radiation transfer, cumulus (CU), microphysics (MP), planetary boundary layer (PBL), and land surface processes, among others. In general, a model consists of a determined physics package, such as the Unified Model (UM) (Cullen, 1993; Brown et al., 2012), the Global Forecast System (GFS), and the Korean Integrated Model (KIM) (Hong et al., 2018), whereas some models, such as the Weather Research and Forecasting (WRF) model and the community Noah land surface model (LSM) with multiparameterization options (Noah-MP) (Niu et al., 2011), have each physical category with several optional parameterization schemes so that users can select the schemes.

Sensitivity experiments of the WRF model to physical parameterization schemes in simulating rainfall systems have been conducted (e.g., Crétat et al., 2012; Cohen et al., 2015; Evans et al., 2012; Song and Sohn, 2018). Cohen et al. (2015) introduced the various PBL schemes employed by the WRF model and examined sensitivity to PBL schemes in simulating cold season severe weather events occurring in the southeastern United States. Evans et al. (2012) explore the performance of various combinations of PBL, CU, MP, and radiation schemes for rainfall events near east coast lows. They suggest that the Mellor–Yamada–Janjic PBL scheme and the Betts–Miller–Janjic CU scheme can be selected in this region with some robustness.

Previous studies on numerical weather and climate prediction models have manually examined one or several suitable schemes for rainfall events as a sensitivity test; hence, not all schemes can be investigated – only selected combinations of schemes. This is because it takes a lot of computer resources and time to conduct sensitivity tests with all physics schemes. However, in addition to model development, all physics schemes need to be explored to simulate more accurate local weather and climate systems. Combinatorial optimization for a system or model configuration has been applied to the water distribution system (Behzadian et al., 2009; Gupta et al., 1999; Simpson et al., 1994; Weng and Liaw, 2005), groundwater monitoring design (Babbar-Sebens and Minsker, 2010), and multi-reservoir operation (Chen et al., 2017). Many algorithms have been developed for combinatorial optimization such as simplex, integer programming, simulated annealing, and genetic algorithm (GA).

In weather and climate prediction models, most applications of GA have focused on the optimization of empirical parameters in the model to fit the modeled properties (e.g., precipitation) to observed counterparts (e.g., Lee et al., 2006; Yu et al., 2013). However, in this study, we created a micro-GA–WRF (µ-GA–WRF) interface to seek an optimal set of CU, MP, and PBL schemes in the WRF model for rainfall events in terms of quantitative precipitation forecast (QPF). The µ-GA–WRF interface system is designed to automatically extract the optimal scheme combination from physics options in the WRF model. The present study attempts to find the optimal combination of parameterization schemes, which is a new and challenging task. This paper is organized as follows: Sect. 2 illustrates the background of combinatorial optimization; Sect. 3 explains the data and method, including the µ-GA–WRF interface system; a case study is presented in Sect. 4; and Sect. 5 contains concluding remarks.

2 Background of combinatorial optimization

To select a suitable optimization algorithm, we should consider the characteristics of objective functions, control variables, and optimization problems. Jamil and Yang (2013) reviewed and compiled benchmark functions found from all the available literature for global optimization problems. They focused on the diverse properties of objective functions such as continuity, linearity, modality, separability, and dimensionality. In terms of combinatorial optimization, control variables can be discretized and indexed values, and the value of the control variable itself can be meaningless. These discrete control variables make the solution space of the cost function discontinuous. Therefore, it is important to choose an algorithm that can handle these properties.

The GA as an evolutionary algorithm is based on the natural selection of genes (i.e., parameters in the algorithm) to search for the optimum. Research has adopted the GA to solve network system design optimization problems with a growing trend from the end of the 20th century (e.g., Simpson et al., 1994; Halhal et al., 1997; Savic and Walters, 1997; Pilar et al., 1999; Dandy et al., 2001). Gupta et al. (1999) emphasized the fact that the GA has advantages of using discrete variables for optimization and having an insensitive initial solution (i.e., robustness in the initial solution). Azadivar and Tompkins (1999) applied the GA approach to optimize qualitative variables (e.g., structural design) in a manufacturing system as simulation optimization. The GA coupled with a simulation model generator searches for the different combinations of design configurations and evaluates the simulations. Gupta et al. (1999) show that the GA provided a lower cost design of water distribution networks (e.g., pipe networks) compared to the nonlinear programming technique. Weng and Liaw (2005) established a combinatorial optimization model, called the Sewer System Optimization Model for Layout & Hydraulics, to optimize cost-effective designs for an urban sewer system. Better alternate network layouts were created more productively by applying the GA. Davis et al. (2019) optimized a malaria model with the GA by clustering locations based on the relationships between malaria and environmental drivers (e.g., temperature, precipitation, and vegetation index). To predict environmentally driven malaria outbreaks across a heterogeneous region, the GA optimized the number of clusters and the environmental predictors for the districts in each cluster in the malaria model.

Furthermore, the GA was applied for combinatorial optimization to the Noah-Multiparameterization (Noah-MP) land surface model (LSM), which can be coupled with the numerical weather prediction model (e.g., WRF model), in Hong et al. (2014). Noah-MP was augmented with multiple physics options for 10 different land surface processes such as phenology, snow, and groundwater (Niu et al., 2011). Hong et al. (2014) performed scheme-based model optimizations in simulating evapotranspiration and runoff (i.e., water balance) in Noah-MP over the Han River basin in South Korea. In addition, they showed a potential advantage of the Noah-MP and GA coupled system to model diagnosis – the evolutionary process provides information on the sensitivity and interrelationship of physics schemes with regard to further model calibrations and improvements. Hong et al. (2015) further evaluated the applicability of the coupling system of micro-GA (i.e., an efficient version of GA; µ-GA) and Noah-MP to larger and multiple regions in East Asia.

The GA does not perform a random search of the extrema, but performs a gradual search toward the extrema. However, it does not directly use gradient information on an objective function, but instead mimics the evolutionary method to quickly reach the global optima. The gradual search is based on the fact that the best individual stands for the nearest point to the optima. In the case of temperature as a physical quantity with a continuous nature (e.g., real number), an increase or decrease in temperature is meaningful. On the other hand, for physics schemes in NWP models with a discontinuous nature (e.g., integer), the option of schemes as an index has no physical meaning. The best option for physics schemes is not related to the nearest options. Thus, the random search is more appropriate to optimize each of the schemes; however, the use of evolutionary algorithms is reasonable when looking for a combination of physics schemes. Note that combinatorial optimization must consider randomness to avoid falling to the local optimum. The µ-GA conducts a global search through a random number generator and a crossover operator; hence, the µ-GA is a very useful tool for combinatorial optimization.

3 Data and methods

3.1 Observation data

Observation data have sometimes strongly affected verification results (Rossa et al., 2008). Merged gauge–radar precipitation has the greatest advantage of the spatially uniform information available. The composite precipitation data produced by the Korea Meteorological Administration (KMA) Radar-Automatic weather station Rain-rate (RAR) system (Suk et al., 2013) using 11 radars was employed to optimize the combination of the physics schemes in the WRF model. The observational domain covers 1241 km × 1761 km on the Korean Peninsula, centered at 38^∘ N and 126^∘ E in the Lambert conformal conic projection. It has enough horizontal (i.e., 1 km) and temporal resolution (i.e., 10 min) to compare with the precipitation fields obtained by the high-resolution model. The performance of the RAR system was examined for 10 heavy rainfall cases selected during the summer of 2006 in Suk et al. (2013), obtaining a squared correlation coefficient (R²) of 0.84 between RAR-estimated rainfall and the observed daily accumulated rainfall from rain gauges.

For comparison with model output, RAR-estimated rainfall data were aggregated to a 5 km resolution grid. A downscaled grid box represents the average of 25 original grid boxes. For reasonable representativeness of samples, we take the average if more than 80 % of original grid boxes have meaningful values.

3.2 µ-GA–WRF interface system

3.2.1 µ-GA

The GA developed by John Holland in the 1970s is a global optimization approach based on the Darwinian principle of natural selection: stronger individuals in a generation are more likely to produce offspring. The aim of the GA is to find the best individual with either a maximum or minimum fitness by means of a stochastic global search of the solution space, through the generations. The algorithm applies crossover and mutation operators to avoid local maximum and minimum solutions. The µ-GA (Krishnakumar, 1989) is an improved version of GA with smaller population sizes (e.g., five) and simplifies the generation-to-generation evolution, hence efficiently reducing the computational resources. To simplify the algorithm, the mutate operator is not used, but the crossover operator is used to increase the diversity at a rate of 100 %. Furthermore, whenever inner loop convergence is achieved, the new population for the next generation consists of all new random individuals, except one elite individual. Thus, the µ-GA can avoid trapping into the local optimum.

The flowchart of the µ-GA interfaced to the WRF model is provided in Fig. 1. The selection operator in the µ-GA is tournament selection with a shuffling technique to choose random pairs for mating. A fitness function to evaluate each individual is of the utmost importance in the GA, and it should be designed taking into account the perspective of the optimization. If the inner loop does not converge, selection is performed and all populations go through a crossover process; then one of the populations is altered by the elite. Here, the elite individual from the previous generation is saved as one of the populations in the current generation (i.e., the elitism). Since the crossover probability of 1.0 is used without the mutation operator in the µ-GA, each individual quickly resembles the elite through generations (i.e., inner loop). In other words, the optimization within the inner loop has a feature of local search by specifically exploring the solution space around the elite. The µ-GA decides that the inner loop converges upon an optimum when the different chromosomes between the elite individual and all the others are less than 5 % as binary bits. After the inner loop convergence, all individuals in the next generation, except for one surviving elite of the parent generation, are regenerated using a random number generator, thus widening the search space (i.e., outer loop; global search). As a criterion of the outer loop convergence to finalize the algorithm, we commonly set the maximum number of generations.

Figure 1 The flowchart of the µ-GA–WRF interface system.

3.2.2 µ-GA–WRF interface system

We created the µ-GA–WRF interface system to seek the optimal physics set of CU, MP, and PBL schemes in the WRF model for rainfall events in terms of QPF. The WRF model, a mesoscale NWP system, has been developed for atmospheric research and operational forecasting applications from the latter 1990s by a collaborative partnership of the National Center for Atmospheric Research (NCAR), the National Oceanic and Atmospheric Administration (represented by the National Centers for Environmental Prediction (NCEP) and the Earth System Research Laboratory), the US Air Force, the Naval Research Laboratory, the University of Oklahoma, and the Federal Aviation Administration (FAA). The details of model configuration are addressed in Sect. 3.3.

Figure 1 shows the flowchart of the µ-GA–WRF interface system. In the µ-GA–WRF interface system, the µ-GA controls WRF simulations in the process of “Compile and Run WRF”. First, the µ-GA randomly initializes individuals (i.e., combinations of physics schemes) in the first generation. Through the inner and outer loop, the µ-GA evaluates the fitness value calculated by WRF results and RAR observations, then reruns the WRF models with new scheme combinations. Finally, if both the inner loop and outer loop converge, the optimization process in the µ-GA–WRF interface system is finished.

The µ-GA is implemented as input parameters. We set a population size of five, meaning that each generation has five individuals (i.e., model simulations). The maximum value of generations is set to 100, which is typically used in µ-GA experiments. The number of parameters (groups of bits) of each individual for the µ-GA is three, which is the number of schemes to be optimized. We used single-point crossover.

3.2.3 Fitness function

Fitness is the basis for evaluating the superiority among individuals consisting of combinations of chromosomes. Designing a fitness function in the GA is critical for optimizing the model as intended. In this study, we are trying to improve the model simulation in terms of QPF. Thus, we used the equitable threat score (ETS; Hamill, 1999) as the fitness function, also called an objective function. The fitness is computed by a sum of ETSs within each precipitation threshold.

$$\begin{array}{}\text{(1)}& \mathrm{Fitness}=\left\{\begin{array}{ll}{\mathrm{ETS}}_i,& \text{for detection},i=\mathrm{3}.\\ {\sum }_i{\mathrm{ETS}}_i,& \text{for heavy rainfall},\\ & i=\mathrm{10},\mathrm{20},\mathrm{30},\mathrm{\dots },\mathrm{300}.\end{array}\right.\end{array}$$

Here, i is a specified threshold of accumulated precipitation in millimeters, and the ETS and chance are defined as

$$\begin{array}{}\text{(2)}& {\displaystyle }\mathrm{ETS}={\displaystyle \frac{\mathrm{hits}-\mathrm{chance}}{\mathrm{hits}+\mathrm{misses}+\mathrm{false}\mathrm{alarms}-\mathrm{chance}}},\text{(3)}& {\displaystyle }\begin{array}{rl}& \mathrm{chance}=\\ & {\displaystyle \frac{(\mathrm{hits}+\mathrm{misses})(\mathrm{hits}+\mathrm{false}\mathrm{alarms})}{(\mathrm{hits}+\mathrm{misses}+\mathrm{correct}\mathrm{negatives}+\mathrm{false}\mathrm{alarms})}}.\end{array}\end{array}$$

ETS has values in the range from $-\mathrm{1}/\mathrm{3}$ to 1. The closer the ETS is to unity, the better the forecast skill. On the other hand, if ETS is equal to or less than 0, the forecast skill is the same as, or even worse than, that of a random forecast. Here, hits, misses, false alarms, and correct negatives from a 2×2 contingency table are estimated by the joint distribution of binary (yes/no) forecasts and observations (Table 1). Rainfall estimations can be evaluated through the table that explicitly provides prediction capability and types of errors in the prediction.

Table 1 A 2×2 contingency table.

When focusing on heavy rainfall, the total fitness value is calculated by the sum of ETS at the threshold ranging from 10 to 300 mm with an interval of 10 mm, whereas when focusing on precipitation detection, a precipitation threshold of 3 mm is used. In order to detect precipitation, the threshold of precipitation accumulated over 24 or 12 h is generally used as a value between 0.1 and 0.3 mm (Rossa et al., 2008). Park et al. (2018) obtained the threshold value of 3 mm h⁻¹ for the station average precipitation rate when the cumulative percentage of warm season precipitation events in Korea reached approximately 80 % based on automatic weather station (AWS) observation data. As we would like to improve the forecast accuracy of precipitation in Korea, we selected the threshold value of 3 mm for hourly precipitation for the calculation of ETS. In this study, we also conducted a sensitivity test of the precipitation accumulation period.

3.3 Experimental design

The WRF model (version 4.0.3) was initialized at 00:00 UTC on 5 August 2018 with the 6-hourly initial and boundary conditions given by the National Center for Environmental Prediction (NCEP) Final (FNL) Operational Model Global Tropospheric Analyses data on 1^∘ × 1^∘ grids. The WRF model configuration is based on the following: horizontal grid spacings of 25 and 5 km for two nested domains (Fig. 2); horizontal grid points of 60×60 and 116×136; the model top of 50 hPa with 33 vertical levels; and the Dudhia shortwave radiation scheme (Dudhia, 1989), Rapid Radiative Transfer Model (RRTM) longwave radiation scheme (Mlawer et al., 1997), revised fifth-generation National Center for Atmospheric Research (NCAR)/Penn State Mesoscale Model (MM5) surface layer scheme (Jimenez et al., 2012), and Unified Noah LSM (Chen et al., 1996; Koren et al., 1999). The control experiment referred to as CTL is simulated with the Kain–Fritsch scheme (KF), the WRF double-moment (WDM) six-class scheme, and the Yonsei University (YSU) scheme as the CU, MP, and PBL scheme, respectively, which are generally used to simulate the precipitation system in Korea. The optimization results from the µ-GA–WRF interface system are referred to as OPT.

Figure 2 WRF nested domains – Domain 1 (d01) and Domain 2 (d02); d01 is centered at 38^∘ N and 126^∘ E in the Lambert conformal conic projection.

We selected the CU, MP, and PBL physical processes for the combinatorial optimization. The CU parameterization determines the prediction of subgrid-scale processes, associated convective clouds, and precipitation at a coarse resolution. Meanwhile, the MP regulates the grid-resolving processes of clouds. The PBL scheme, which could indirectly influence precipitation by interacting with other physics, can affect temperature and moisture profiles in the lower troposphere via exchanges of moisture, heat, and momentum through the mixing associated with turbulent eddies. The options of the CU, MP, and PBL schemes used for the optimization in the µ-GA–WRF interface system are shown in Table 2. If Mellor–Yamada–Janjic (MYJ), quasi-normal-scale elimination (QNSE), Mellor–Yamada Nakanishi and Niino Level (MYNN) 3, or total energy mass flux (TEMF) is selected as the PBL scheme, eta similarity (Monin and Obukhov, 1954; Janjic, 1994, 1996, 2002), QNSE, MYNN, or TEMF should be set as a surface layer scheme, respectively. The surface layer scheme is the lowest part of the atmospheric boundary layer where the surface fluxes (i.e., surface heat, moisture, and momentum fluxes) can be calculated not only by combining the LSM, but also by itself. We have found the best scheme combination by the µ-GA as the mechanical and objective optimization method without model simulations of 2688 (14 × 16 × 12), which is the total possible number of scheme combinations.

Table 2 Summary of CU, MP, and PBL physics schemes used in the µ-GA–WRF interface system for the optimization.

During the 12 h period from 12:00 UTC on 5 August 2018 to 00:00 UTC on 6 August 2018, including the first and second periods of intense rainfall (see Sect. 4.1), precipitation was evaluated by fitness functions. We perform the optimization experiments with five different fitness functions based on ETS in Sect. 3.2.2. Table 3 shows the summary of experiments – OPT-EXP1 for 12-hourly accumulated precipitation with precipitation thresholds ranging from 10 to 300 mm with the interval of 10 mm; OPT-EXP2 for 12-hourly accumulated precipitation with a precipitation threshold of 3 mm; OPT-EXP3 for all 6-hourly accumulated precipitation during the evaluation period with a precipitation threshold of 3 mm; OPT-EXP4 for all 3-hourly accumulated precipitation during the evaluation period with a precipitation threshold of 3 mm; OPT-EXP5 for all hourly accumulated precipitation during the evaluation period with a precipitation threshold of 3 mm. For OPT-EXP1, the total fitness value is calculated as the sum of the ETSs at all thresholds, while for OPT-EXP2–OPT-EXP5 it is calculated as the average of the ETSs for each accumulated time over 12 h.

Table 3 The summary of the experiments. In the fitness function, i and t represent the precipitation threshold and the number of accumulated time instants, respectively. N is the total number of accumulated time instants.

4 Case study

4.1 Case description

A coastal flood occurred in Korea due to a quasi-stationary mesoscale convective system (MCS), which produced heavy rainfall on 5–6 August 2018. As the event was unexpected, the back-building MCS was located in the Yeongdong region for about a day; thus, the heavy rainfall caused damage to property estimated at KRW 177 million (KMA, 2018). For the period from 11:00 UTC on 5 August to 14:00 UTC on 6 August, 294.5 mm of precipitation was recorded at Sokcho: the first intense rainfall continued for 4 h (13:00–17:00 UTC 5 August) with a maximum precipitation rate of 35.3 mm h⁻¹ and total rainfall amount of 83.5 mm, whereas the second intense rainfall (17:00 UTC 5 August–00:00 UTC 6 August) recorded a maximum precipitation rate of 54.9 mm h⁻¹ and total rainfall amount of 192 mm due to the quasi-stationary MCS. To predict more accurately, forecasters essentially need the mesoscale information from NWPs as well as synoptic weather charts, vertical soundings, satellite observations, and weather station observation, among others, at the preceding time. NWP models can capture the important triggers, which can hardly be found through observations, to predict rainfall. This heavy rainfall case occurred due to mesoscale factors: (1) low-level convergence, (2) strong vertical wind shear, (3) coastal fronts and back-building convection bands, (4) mid-level advection of cold air and positive relative vorticity, and (5) vigorous updraft releasing potential instability (Park and Park, 2020). Therefore, it is necessary to improve the NWP model to more accurately identify these mesoscale factors.

4.2 Results

4.2.1 Combinatorial optimization of the physics schemes for QPF

The combinatorial optimization of the physics schemes in the WRF model targets the improved quantitative forecasting of heavy rainfall. OPT-EXP1 shows the simulated results using the optimized combination of the MP, CU, and PBL schemes, focusing on strong precipitation intensity. Figure 3 depicts the evolution of generations of the µ-GA, represented by fitness values, for OPT-EXP1. The µ-GA reached the maximized evolution, which was the point that the best individual in each generation converged upon the highest fitness score (here, of 4.292) at the 12th generation. Before that, the local optimum (i.e., intermediate optimum; IMD-OPT) appeared at the fourth generation with a fitness of 2.9896. The optimized schemes of the CU, MP, and PBL for this event are the MSKF, NSSL two-moment, and YSU scheme, respectively. For the IMD-OPT-EXP1, only the MP physics scheme selected as WSM six-class is different from the global optimum. The optimum of the PBL scheme is the same as CTL, and that of the CU scheme is the updated scheme from the CU scheme for CTL. The KF (i.e., the CU scheme for CTL) is suitable for a horizontal resolution of ∼ 25 km, at which convective clouds can be represented explicitly. However, MSKF (i.e., the optimum of the CU scheme) has been improved for use in the so-called gray zone scales (i.e., 12 to 1 km) and at a horizontal grid spacing of 25 km (Glotfelty et al., 2019).

Figure 3 Evolution of generations in the optimization process, leading to changes in fitness values. OPT and IMD-OPT represent the final optimized phase and one of the intermediate optimized phases, respectively.

Figure 4 compares the ETSs of each precipitation threshold for CTL, OPT-EXP1, IMD-OPT, the best individual at the first generation (GEN1), and REF. REF is simulated with the CU, MP, and PBL physics schemes validated in Park and Park (2020). The Morrison scheme is chosen as the MP scheme, and the CU and PBL schemes are the same as CTL. The selected schemes for CTL, OPT-EXP1, IMD-OPT, GEN1, and REF are summarized in Table 4. OPT-EXP1 performs significantly better than CTL at precipitation thresholds above 20 mm, indicating a remarkable improvement in the ETSs. Although REF also shows more improved forecasting skill than CTL at all precipitation thresholds, OPT-EXP1 performs better at higher precipitation intensity compared to REF. The optimization process from GEN1 to IMD-OPT, corresponding to the evolution of early generations, shows increases in ETS at precipitation thresholds less than 130 mm, whereas IMD-OPT to OPT-EXP1 shows further enhancements at heavy precipitation thresholds (≥40 mm) and even above 130 mm. For both observations and model output, the maximum amount of 12 h accumulated precipitation in the grid box did not exceed 190 mm.

Figure 4 ETS for 12 h accumulated precipitation with each precipitation threshold for CTL (black line), OPT-EXP1 (red line), IMD-OPT (dashed blue line), GEN1 (dashed gray line), and REF (dashed black line).

Table 4 The summary of selected schemes for CTL, OPT-EXP1, IMD-OPT, GEN1, and REF.

The spatial distribution of 12 h accumulated precipitation for observations (RAR), OPT-EXP1, CTL, and REF is shown in Fig. 5. CTL broadly underestimated precipitation, and REF is better than CTL in terms of both quantity and rain cell development. Rain cells located near both Sokcho and north of Gangneung were well captured in REF and OPT-EXP1, although the amount of precipitation was underestimated over the north of Gangneung. OPT-EXP1 shows improved precipitation simulation for the inland area near Gangneung and Daegwallyeong, but it still underestimated. As the MCS stayed near Sokcho and the north of Gangneung during the period of accumulative precipitation, light rainfall occurred over the inland, whereas heavy rainfall occurred over the sea along the coastal line.

Figure 5 The spatial distribution of 12 h accumulated precipitation for RAR, CTL, OPT-EXP1, and REF. Black dots depict locations of weather stations: Sokcho (S; 38.25^∘ N, 128.56^∘ E; 18.06 m), Gangneung (G; 37.75^∘ N, 128.89^∘ E; 26.04 m), and Daegwallyeong (D; 37.68^∘ N, 128.86^∘ E; 772.57 m).

We also verified the effectiveness of the optimization by the continuous statistics for CTL, REF, and OPT-EXP1 (Table 5 and Fig. 6). The scatter plot for OPT-EXP1 exhibits the best performance, with a regression coefficient (R) of 1.01, compared to CTL and REF (Fig. 6). In addition, Table 5 shows that OPT-EXP1 has lower spatial mean bias and root mean square error (RMSE) of precipitation (−7.433 and 21.511) and a greater Pearson's correlation coefficient (PCC) (0.762) than CTL (−8.696, 25.430, and 0.673, respectively). It performs better than REF as well. In conclusion, combinatorial optimization of the physics schemes has enhanced the forecast skill not only in QPF (i.e., ETS) but also in terms of both spatial distribution and continuous statistics.

Figure 6 Scatter plot of CTL, REF, and OPT-EXP1.

Table 5 Continuous statistics for CTL, REF, and OPT-EXP1. The best results are presented in bold.

4.2.2 Sensitivity of fitness functions based on the assessment of precipitation occurrence

For the accuracy of the deep convective precipitation system, we wonder whether it would be effective to increase the accuracy of the precipitation occurrence or to increase the accuracy of precipitation within each precipitation threshold. In this section, we conduct a sensitivity test of the accumulated precipitation time interval used in the fitness function calculation to evaluate precipitation occurrence with the precipitation threshold of 3 mm (see Table 3). The ETSs for 12 h accumulated precipitation, calculated by using 12-hourly (OPT-EXP2), 6-hourly (OPT-EXP3), 3-hourly (OPT-EXP4), and 1-hourly (OPT-EXP5) accumulated precipitation data, were evaluated at each time interval. In contrast to OPT-EXP1, precipitation thresholds for them are set as one criterion (i.e., precipitation threshold of 3 mm), so the maximum value of ETS is equal to 1. When the accumulation time interval becomes shortened (e.g., an hour), precipitation prediction must also be more accurate on a temporal scale in order to have higher fitness because the fitness is computed by the average of ETSs calculated at each time interval. In other words, the shorter the accumulated time intervals, the more ETSs of predicted precipitation are evaluated. Thus, as expected, OPT-EXP2 shows the highest fitness value (i.e., 0.3482), followed by OPT-EXP3 with the fitness of 0.2862 and OPT-EXP5 with the lowest fitness of 0.2249 (Fig. 7). OPT-EXP4 performs similarly as OPT-EXP5, having the fitness of 0.2270. The selected schemes for OPT-EXP2, OPT-EXP3, OPT-EXP4, and OPT-EXP5 are shown in Table 6.

Figure 7 Same as Fig. 3 but for OPT-EXP2 to OPT-EXP5.

Table 6 The selected schemes for OPT-EXP2, OPT-EXP3, OPT-EXP4, and OPT-EXP5.

Figure 8 shows the spatial distribution of 12 h accumulated precipitation for OPT-EXP2 to OPT-EXP5. All experiments underestimate the convective system and overestimate very light precipitation over the inland area of the Korean Peninsula (see Figs. 5 and 8). From the ETS perspective, OPT-EXP2 is the best result, but OPT-EXP3 shows the best simulation result in terms of the spatial distribution. Since no method is absolutely superior to others in precipitation evaluation methods such as ETS, critical success index (CSI), probability of detection (POD), and continuous statistics indices, several indices including spatial distribution must be examined together. Rain cells located near both Sokcho and the north of Gangneung were well captured only in OPT-EXP3 and OPT-EXP4 but were still underestimated. Rainfall over the sea along the coastal line was simulated in OPT-EXP3 and OPT-EXP4 as well. On the other hand, the evaluation of the fitness at 1 h intervals results in poor accuracy, possibly because of including the time phase error of the model.

Figure 8 Same as Fig. 4 but for OPT-EXP2 to OPT-EXP5.

Figure 9 depicts the scatter plot for OPT-EXP2 to OPT-EXP5. In terms of observed precipitation, OPT-EXP3 and OPT-EXP4 have more accuracy than OPT-EXP2. OPT-EXP3 has the best fit for the RAR, showing a regression coefficient of 1.13, followed by OPT-EXP4. Moreover, OPT-EXP4 has the lowest RMSE and the greatest PCC of precipitation (23.952 and 0.731, respectively) (Table 7). OPT-EXP3 has the lowest spatial mean bias (−8.690). In terms of fitness, OPT-EXP2 is superior, but OPT-EXP3 and OPT-EXP4 show better simulations in terms of both the spatial distribution and continuous statistics.

Figure 9 Scatter plot of OPT-EXP2 to OPT-EXP5.

Table 7 Continuous statistics for OPT-EXP2 to OPT-EXP5.

In this section, the sensitivity of the accumulation time interval of precipitation used in the fitness function calculation (i.e., ETS) to the optimization in the µ-GA–WRF interface system was briefly examined. In the current model performance, the best result of the optimization experiments can be obtained by using the 3- or 6-hourly accumulated precipitation in the fitness function when focusing on precipitation detection. However, compared to OPT-EXP1, both the quantitative precipitation and spatial distribution in OPT-EXP1 were much more improved than other experiments (i.e., OPT-EXP2–OPT-EXP5). Therefore, in order to improve the simulations of deep convective systems, it is recommended to evaluate the precipitation accuracy at various precipitation thresholds rather than assessing the accuracy of precipitation occurrence.

4.3 Discussion

All physics schemes including the CU, MP, PBL, radiation, and surface schemes are interrelated, and a nonlinear relationship among them appeared due to the complexity of the atmospheric system. Thus, in order to accurately predict precipitation, it is necessary to explore the combination of physics schemes rather than focusing on an individual scheme. The evolutionary approach to find the optimum combination of the CU, MP, and PBL schemes can provide insightful understanding of the implemented physical schemes and their interrelationships. The accuracy of precipitation of less than 30 mm in large areas has been improved by fitted CU and PBL schemes. On the other hand, the simulation accuracy of high-intensity precipitation occurring in the small area was improved by the MP schemes. It is because the MP schemes control the grid-resolving-scale precipitation, while the CU schemes determine the subgrid-scale precipitation. In other words, a realistic parameterization of cloud microphysics is crucial for precipitation forecast in high-resolution models. Typical cumulus convection can be represented by the CU schemes at a horizontal grid spacing of about 25 km. However, the selected CU scheme (i.e., MSKF) has been improved for use in the so-called gray zone scales (e.g., 5 km used in this study); thus, it can outperform the other CU schemes. On the other hand, the KFCP scheme that is modified to better account for the presence of shallow clouds was selected for OPT-EXP2 and OPT-EXP5, possibly because their fitness functions were focused on the precipitation occurrence. Note that the single-moment MP schemes predict the mixing ratio of hydrometeors by representing the hydrometeor size, while the double-moment schemes also predict the number concentrations of hydrometeors. Thus, the double-moment schemes (e.g., NSSL two-moment, WDM 6, Morrison) can produce a reasonable concentration of large droplets for a heavy precipitation system compared to the single-moment schemes (Lim and Hong, 2010). In addition, the YSU scheme, representing the PBL process, more accurately simulates a deeper vertical mixing in the thermally induced free convection regime covering multiple vertical levels (Hong et al., 2006), thus being superior to the other schemes for the simulated precipitation.

However, it is difficult to insist that the order of fitting scheme is directly related to the importance of the scheme in QPF because of the nonlinear relationship between precipitation and the physics schemes as well as among the physics schemes. Moreover, it can be noted that the combination of the randomly selected schemes in the first generation approaches the optimal solution, allowing the fitness function to converge quickly. For example, in this study, both the PBL and CU scheme are fortunately selected to be the same as the optimum, and this combination has a higher fitness value than the other combinations. Thus, in the µ-GA evolution process, the information on the optimized CU and PBL scheme in the best individual was inherited by the elitism, and the MP scheme of the best individual was changed to be optimized through generations. The GA with a randomly selected initial population is robust in finding solutions as enough generations pass, but still the initial population affects the convergence velocity (i.e., the generation with optimal solution). As a result of sensitivity tests with different initial populations, convergence occurred after the 50th generation in one test, and the MP scheme converged earlier than the CU and PBL scheme in the other. In summary, the optimized results do not depend on the initial population of the first generation, but the initial population may affect which scheme will be optimized first.

The simulation results of Park and Park (2020) are more accurate than REF because of the different domain setting. For this case with localized heavy rainfall, more accurate precipitation simulations can be achieved when a specific region is set as a model domain or when multiple nested domains are used. However, the selected domain in this study may be reasonable to derive a general scheme combination that accurately simulates precipitation over the Korean Peninsula. It is necessary to derive a general set of physics schemes for accurate precipitation simulations through several case studies as a further study.

5 Conclusions

The uncertainties related to the subgrid-scale parameterizations significantly increase as NWP models become more complex. The accuracy of subgrid-scale parameterizations depends on both parameters in the physics schemes and the choice of the parameterization schemes for each corresponding physical process. In this study, we created the µ-GA–WRF interface to seek the optimal set of physics parameterization schemes in the WRF model. The GA is founded on natural selection and evolution to search for the optimum, and the µ-GA is an efficient version of the GA. In weather and climate studies, most GA applications have focused on optimizing the empirical parameters of NWP models to represent a real atmospheric system, while the current study attempts to find the optimal combination of the parameterization schemes, which is a novel and challenging task. Because of the nonlinear relationship among the physics schemes, it is recommended to optimize several interesting schemes concurrently in the WRF model rather than optimizing the schemes in each physics category sequentially. The GA is an appropriate optimization method in that it can handle the nonlinearity of the parameters to be optimized.

The experiments were conducted on the optimal set of the CU, MP, and PBL schemes in terms of QPF for a heavy rainfall event in Korea through the µ-GA–WRF interface system. The µ-GA successfully improved simulated precipitation in spite of the nonlinear relationship between precipitation and the physics schemes as well as among the physics schemes. The µ-GA reached its maximum evolution in the 12th generation and led to significant improvement in the ETSs, especially at the threshold range of 20–180 mm. The optimized CU, MP, and PBL schemes for this event are the MSKF, NSSL two-moment, and YSU scheme, respectively. During the evolution process, the MP schemes control grid-resolving-scale precipitation, while the CU and PBL schemes determine subgrid-scale precipitation.

We also conduct a sensitivity test of the accumulated precipitation time interval used in the fitness function (i.e., ETS) with a precipitation threshold of 3 mm. The best result of the optimization experiments has been obtained by using the 3- or 6-hourly accumulated precipitation. However, in order to improve the simulation of deep convective systems, it is recommended to evaluate the accuracy of precipitation at various precipitation thresholds (i.e., precipitation thresholds of 10, 20, 30, …, 300 mm) rather than assessing the accuracy of precipitation occurrence (e.g., precipitation threshold of 3 mm).

In recent studies, optimization experiments for parameter estimation for multiple heavy rainfall events have been conducted to obtain a set of optimal parameters to improve precipitation prediction (e.g., Duan et al., 2017; Di et al., 2018). We address the fact that the optimized scheme set obtained in this study is specific to the selected rainfall case or to the rainfall systems that occur under a similar synoptic and geographical environment; thus, it is not robust to all the precipitation cases in Korea, which depend on different mechanisms of initiation and development. As a future study, we plan to perform the combinatorial optimization of physical parameterization schemes for several heavy rainfall cases under the same category in terms of location and synoptic environment, expecting to find an optimal scheme set robust to the heavy rainfall systems in that category.

Note that prior to model simulation unknown parameters and schemes should be fitted to the regional weather and climate to reduce considerable uncertainties in NWP models. In addition, in terms of model development, all physics schemes need to be explored to simulate more accurate local weather and climate systems if sufficient computer resources and time are available. This study has demonstrated that the combinatorial optimization of physics schemes in the WRF model is one possible solution to enhance the forecast skill of regional or local prediction. We also significantly reduced the number of model simulations for optimization using the GA, one of the artificial intelligence methods. Furthermore, experiments for combined scheme-based with parameter-based optimization are essentially required to investigate the effect of parameter calibrations on the model sensitivity to scheme selections. As a further study, we strongly suggest conducting comprehensive parameter and scheme estimation to improve the model performance.