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

. One of biggest uncertainties in Numerical Weather Predictions (NWPs) comes from treating the subgrid-scale physical processes. For the more accurate regional weather/climate prediction by improving physics parameterizations, it is important to optimize a combination of physics schemes as well as unknown parameters in NWP models. We have developed an interface system between 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 15 precipitation forecast for heavy rainfall events in Korea. The μ -GA successfully improved simulated precipitation despite the non-linear 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 has demonstrated the combinatorial optimization of physics schemes in the WRF model is one of possible solutions to enhance the forecast skill of precipitation.


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). It should be also provided more accurate initial conditions. One of the biggest uncertainties in Numerical Weather Predictions (NWPs) comes from treating the subgrid-scale physical processes that have not been sufficiently 25 understood. The subgrid-scale physical processes are parameterized in NWP models through empirical evidences, such as the derived value from observation and/or theoretical backgrounds. Therefore, the accuracy of physics parameterizations strongly depends on the followings: the value of parameters in given uncertainty ranges in parameterization schemes; the choice of parameterization schemes for each corresponding physical process. Note that prior to model simulation the 3 2 Background of combinatorial optimization To select a suitable optimization algorithm, we should consider the characteristics of objective functions, control variables as 65 well as 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 70 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 GA has advantages of using discrete variables for optimization and 75 having an insensitive initial solution (i.e., robustness in the initial solution). Azadivar and Tompkins (1999) applied the GA approach to optimize to 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) shows that the GA provided a lower cost design of water distribution networks (e.g., pipe networks), compared to the non-linear programming technique. Weng and Liaw (2005) established a 80 combinatorial optimization model, called the Sewer System Optimization Model for Layout & Hydraulics, to optimize costeffective designs for an urban sewer system. The 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 85 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 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 90 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 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.

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

Observation data
The observation data have strongly affected the 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 110 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 in 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 RAR system was examined for 10 heavy rainfall cases selected 115 during the summer of 2006 in Suk et al. (2013), obtaining a squared correlation coefficient (R 2 ) of 0.84 between RARestimated 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. 120

μ-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, 125 through the generations. The algorithm applies the crossover and mutation operators to avoid local maximum/minimum solutions. The μ-GA (Krishnakumar, 1989) is an improved version of GA with smaller population sizes (e.g., of 5) and simplifies the a 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 130 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. 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 inner loop does not converge, selection is performed and all populations go through a crossover process, 135 then one of populations is altered by the elite. Here, the elite individual from the previous generation is saved as one of 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 exploring specifically the solution space around the elite.
The μ-GA decides that thean inner loop converges upon an optimum when the different chromosomes between thean elite 140 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 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.

μ-GA-WRF interface system 145
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, one of mesoscale NWP systems, has been developed for atmospheric research and operational forecasting applications from the latter 1990s by a collaborative partnership of the  The μ-GA is implemented as input parameters. We set population size of 5, meaning that each generation has 5 individuals (i.e., model simulations). The maximum value of generations is set to 100, typically used in μ-GA experiments. The number of parameters (groups of bits) of each individual for the μ-GA is 3, which is the number of schemes to be optimized. We 160 used single-point crossover.

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; Hamil, 1999) as the fitness function, also called 165 objective function. The fitness is computed by sum of ETSs within each precipitation threshold: Where i isare specified thresholds of accumulated precipitation in mm and the ETS and chance is defined as: .
( 3)  170 ETS hasve values in the range from -1/3 to 1. The closer the ETS is to a unity, the better the forecast skill. On the other hand, While 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. 175 When focusing on heavy rainfall, 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 the precipitation detection, a precipitation threshold of 3 mm is used. In order to detect precipitation, the threshold of precipitation accumulated over 24 or 12 hours is generally used as a value between 0.1 and 0.3 mm (Rossa et al., 2008).  obtained the threshold value of 3 mm h -1 for the station average precipitation rate when the cumulative percentage of warm-season precipitation events in Korea reached 180 approximately 80 % based on 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 the sensitivity test of precipitation accumulation period.

Experimental design
The WRF model (version 4.0.3) was initialized at 00:00 UTC 5 August 2018 with the 6-hourly initial and boundary 185 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 followings: 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; 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 190 (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 Kain-Fritsch scheme (KF), the WRF Double-Moment (WDM) 6-class scheme, and the YonSei University (YSU) scheme as CU, MP, and PBL scheme, respectively, generally used to simulate precipitation system in Korea. The optimization results from the μ-GA-WRF interface system are referred to as OTPOPT. 195 We selected the CU, MP, and PBL physical processes for the combinatorial optimization. The CU parameterization determines the prediction of sub-grid scale processes associated convective clouds and precipitation at a coarse resolution.
Meanwhile, the MP regulates the gridgird-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,200 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, 1994Janjic, , 1996Janjic, , 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) 205 can be calculated not only by combining the land-surface modelLSM, but also by itself. We have found the best scheme combination by the μ-GA as the mechanical and objective optimization method without a model simulations of 2,688 (= 14×16×12) which is the total possible number of scheme combinations.
During the 12-hr period from 12:00 UTC 5 to 00:00 UTC 6 August 2018, including the first and second periods of intense rainfall (see Section 4.1) precipitation was evaluated by fitness functions. We perform the optimization experiments with 5 210 different fitness functions based on ETS in Sec. 3.2.2. Table 3 shows the summary of experiments -OPT-EXP1 for 12 hourly accumulated precipitation with precipitation thresholds ranged from 10 to 300 mm with the interval of 10 mm; OPT-EXP2 for 12 hourly accumulated precipitation with precipitation threshold of 3 mm; OPT-EXP3 for all 6 hourly accumulated precipitation during the evaluation period with precipitation threshold of 3 mm;. OPT-EXP4 for all 3 hourly accumulated precipitation during the evaluation period with precipitation threshold of 3 mm; OPT-EXP5 for all hourly 215 accumulated precipitation during the evaluation period with 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 hours.

Case description
A coastal flood occurred in Korea due to a quasi-stationary Mmesoscale Cconvective Ssystem (MCS) which produced heavy rainfall on 5-6 August 2018. As the unexpected event, the back-building MCS was located in Yeongdong region about a day, thus the heavy rainfall caused damage to property estimated at 177 million won (KMA, 2018). For the period from 11:00 UTC 5 to 14:00 UTC 6 August, 294.5 mm of precipitation was recorded at Sokcho: the first intense rainfall continued for 4 hrs (13:00 UTC -17:00 UTC 5 August) with the maximum precipitation rate of 35.3 mm/h and total rainfall amount of 83.5 225 mm, whereas the second intense rainfall (17:00 UTC 5 -00:00 UTC 6 August) recorded the 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, weather station observation, etc., at the preceding time. The NWP models can capture the important triggers which can be hardly found out through observations to predict rainfall. This heavy rainfall case occurred due to the 230 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.

Combinatorial optimization of the physics schemes for QPF
The combinatorial optimization of the physics schemes in the WRF model is targeting 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 240 generation converged upon the highest fitness score (here, of 4.292), at the 12thnd generation. Before that, the local optimum (i.e., intermediate optimum; IMD-OPT) appeared at 4thrd generation with a fitness of 2.9896. The optimized schemes of the CU, MP, and PBL for this event are MSKF, NSSL-2 moment, and YSU scheme, respectively. For the IMD-OPT-EXP1, only the MP physics scheme selected as WSM 6-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., 245 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 grey zone scales (i.e., 12 to 1 km) as well as at horizontal grid spacing of 25 km .  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 the 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 255 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 observation and model output, the maximum amount of 12 -hour accumulatedion precipitation in thea grid box did not exceed 190 mm.
The spatial distribution of 12-hr accumulated precipitation for observation (RAR), OPT-EXP1, CTL, and REF are shown in 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 the regressioncorrelation coefficient (R) of 1.01, compared to CTL and REF (Fig. 6). In addition, Table 5 shows OPT-EXP1 has lower spatial mean bias and Root Mean Square Error (RMSE) of precipitation (-7.433 and 21.511) and 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 270 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.

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 275 this section, we conduct the sensitivity test of accumulated precipitation time interval used in fitness function calculation to evaluate the precipitation occurrence with thea precipitation threshold of 3 mm (see Table 3). The ETSs for 12-hr accumulated precipitation, calculated by using 12 hourly (OPT-EXP2), 6 hourly (OPT-EXP3), 3 hourly (OPT-EXP4), 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 = 3 mm), so the maximum value of ETS 280 is equal to 1. When the accumulation time interval become shortened (e.g., an hour), the precipitation prediction must also be more accurate on a temporal scale in order to have the 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 similar as 285 OPT-EXP5, having the fitness of 0.2270. The selected schemes for OPT-EXP2, OPT-EXP3, OPT-EXP4, and OPT-EXP5 areis shown in Table 6. Figure 8 shows the spatial distribution of 12-hr accumulated precipitation for OPT-EXP2 to OPT-EXP5. All experiments underestimate the convective system and overestimated very light precipitation over the inland area of Korean Peninsula (see Fig. 5 and Fig. 8). From the ETS perspective, OPT-EXP2 is the best result, but OPT-EXP3 shows the best simulation results 290 in terms of the spatial distribution. Since no method is existed absolutely superior to others in precipitation evaluation methods such as ETS, CSI, POD, and continuous statics indices, several indices including spatial distribution must be examined together. Rain cells located near both Sokcho and 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 hour intervals results in poor accuracy, possibly because 295 of including the time phase error of the model. 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 lead to better agreement with observed precipitation than OPT-EXP2. OPT-EXP3 has the best fit for the RAR, showing the regression coefficient of a highest correlation (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 300 the lowest spatial mean bias (-8.690). In terms of fitness, OPT-EXP2 is superior, but OPT-EXP3 and OPT-EXP4 shows the better simulations in terms of both the spatial distribution and continuous statistics.
In this section, the sensitivity of the accumulation time interval of precipitation used in 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 theas a fitness 305 function when focusing on precipitation detection. However, compared to OPT-EXP1, both the quantitative precipitation and spatial distribution in OPT-EXP1 was 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 evaluation of the precipitation accuracy at various precipitation thresholds rather than assessing the accuracy of precipitation occurrence.

Discussions 310
All physics schemes including the CU, MP, PBL, radiation, and surface schemes are interrelated, and a non-linear relationship among them is appeared due to the complexity of atmospheric system. Thus, in order to accurately predict precipitation, it is necessary to explore the combination of physics schemes rather than focusing only on individual schemeeach one. 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 315 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 occurred in thea 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 sub-grid scale precipitation. In other words, a realistic parameterization of cloud microphysics is crucial for the precipitation forecast in high resolution models. Typical cumulus convection can be represented by the CU schemes at horizontal grid spacing of about 25 km. 320 However, the selected CU scheme (i.e., MSKF) has been improved for use in the so-called grey 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 number concentrations of 325 hydrometeors. Thus, the double-moment schemes (e.g., NSSL 2-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 , thus being superior to the other schemes for the simulated precipitation. In addition, selected CU scheme (i.e., MSKF scheme) also has been improved 330 for use in the so-called grey zone scales (e.g., 5 km used in this study).
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 non-linear relationship between precipitation and the physics schemes as well as among the physics schemes.
Moreover, it can be note 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 335 are fortunately selected 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 of 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 340 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) is more accurate than REF because of the different domain setting. For this case as a localized heavy rainfall, more accurate precipitation simulations can be achieved when a specific region is set as a 345 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.
The uncertainties related to the subgrid-scale parameterizations significantly increase as NWP models become more complex. 350 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 the physics parameterization schemes in the WRF model. The GA is founded on the natural selection and evolution to search the optimum, and the μ-GA is an efficient version of the GA. In weather and climate studies, most the GA applications have focused on optimizing the empirical parameters of NWP models to represent a real atmospheric 355 system, while the current study attempts to find thean optimal combination of the parameterization schemes, 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. 360 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 non-linear relationship between precipitation and the physics schemes as well as among the physics schemes. The μ-GA reached its maximum evolution in the 12thnd generation and led to the significant improvement in the ETSs, especially at thea threshold range of 20 -180 mm. The optimized CU, MP, and PBL schemes of CU, MP, and PBL for this event are 365 MSKF, NSSL-2 moment, and YSU scheme, respectively. During the evolution process, the MP schemes control gridresolving scale precipitation while the CU and PBL schemes determine sub-grid scale precipitation.
We also conduct the sensitivity test of accumulated precipitation time interval used in fitness function (i.e., ETS) with precipitation threshold of 3 mm. The best result of the optimization experiments has obtained by using the 3 or 6 hourly accumulated precipitation. However, in order to improve the simulation of deep convective systems, it is recommend to 370 evaluate the accuracy of precipitation at various precipitation thresholds (i.e., precipitation thresholds = 10, 20, 30, …, 300 mm) rather than assessing the accuracy of precipitation occurrence (e.g., precipitation threshold = 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 the precipitation prediction (e.g., Duan et al., 2017;Di et al., 2018). We address that the optimized scheme set obtained in this study is specific to the selected rainfall case or at the best to the 375 rainfall systems that occur under 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. 380 Note that prior to model simulation unknown parameters and schemes should be fitted to the regional weather/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 of possible solutions to enhance the forecast skill of the regional or local prediction. We also significantly reduced the number 385 of model simulations for optimization using the GA, one of the artificial intelligence methods. Furthermore, the 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 to conduct comprehensive parameter and scheme estimation to improve the model performance.

Code and data availability 390
The current version of the WRF model is available from the github website: https://github.com/wrf-model/WRF. The GA code, used in this study, was developed by David L. Carroll and last updated on 2 April 2001. The current version of the GA driver is available from the website: https://cuaerospace.com/products-services/genetic-algorithm/ga-drive-free-download.
The exact versions of both the WRF model and the GA driver, used to produce the results in this study, are archived on Zenodo (https://doi.org/10.5281/zenodo.5076930), along with the input data, namelist files, and scripts to run the model and 395 produce the plots of all the simulations presented in this study. The NCEP FNL Operational Model Global Tropospheric Analyses data, used for the initial and boundary conditions of the WRF model, can also be downloaded from the website of Research Data Archive of NCAR: https://rda.ucar.edu/datasets/ds083.2/. The RAR-estimated rainfall data were obtained by the Korea Meteorological Administration (KMA). The KMA does not provide this data set through the public service, called "Open MET Data Portal" (https://data.kma.go.kr/resources/html/en/ncdci.html), but one can obtain the data via separate 400 request to the KMA.