Knee Point‐Based Multiobjective Optimization for the Numerical Weather Prediction Model in the Greater Beijing Area

Determination of the optimal parameter values in numerical weather prediction (NWP) models has a significant impact on predictions. Here, we propose a knee point‐based multiobjective optimization (KMO) method to find an optimal solution of the NWP model parameters. We apply it to optimize the Weather Research and Forecasting (WRF) model's summer precipitation and temperature simulations for the Greater Beijing area. The results showed that it required fewer than 125 samples (i.e., 25 times the number of dimensions of the parameter space) to obtain the WRF model's optimal parameter values. The optimal parameters determined by KMO outperform the default parameters in WRF simulations for summer precipitation and temperature prediction in the Greater Beijing area, across all periods (calibration, validation, and testing). Additionally, clear physical interpretations are provided to explain why the optimal parameters lead to improved precipitation and temperature forecasting. Overall, the proposed method is effective and efficient to improve NWP.

• A knee point-based multiobjective optimization (KMO) method was proposed for numerical weather prediction (NWP) • Optimal parameter set can be found by running the NWP model 25 times the number of dimensions of the parameter space • The optimal parameters determined by KMO outperform the default parameters, across all periods (calibration, validation, and testing) Supporting Information: Supporting Information may be found in the online version of this article.
parameter uncertainty is critical for assessing the ability of the NWP model to emulate actual weather system behaviors (Aksoy, 2015;Duan et al., 2017).
The NWP model parameters are physical constants that describe the constitutive laws of a weather system.In general, the parameters cannot be accurately determined from observations or physical principles but need to be optimized.Parameter optimization (i.e., parameter estimation, parameter calibration, or parameter tuning) refers to a process in which model parameters are tuned by minimizing the objective function to achieve the best match between the model predictions and the actual observations.Initial studies typically determined the parameter values based on the experience of the modeler or specialist (i.e., trial and error).But the manual trial and error method based on experience is subjective, in order to overcome the defect, the utility of automatic parameter optimization method of the NWP model that tunes model parameters automatically to achieve accurate model predictions has been suggested (Afshar et al., 2020;Chinta & Balaji, 2020;Chinta et al., 2021;Di et al., 2015Di et al., , 2018Di et al., , 2019;;Duan et al., 2017;Ollinaho et al., 2014;Yang et al., 2012Yang et al., , 2015)).Duan et al. (2017) used the adaptive surrogate model-based optimization (C.Wang et al., 2014) to optimize nine WRF model's parameters to improve precipitation and temperature simulations.This is a method of converting multiobjective optimization (MO) to single objective optimization that only optimizes the transformed objective function.However, during the optimization process, the change state of each objective function is unknown, and only one objective function may improve, whereas all other objective functions worsen.In a study conducted by Chinta and Balaji (2020), they employed multiobjective adaptive surrogate model-based optimization (MO-ASMO), as proposed by Gong et al. (2016), to calibrate nine sensitive parameters of the WRF model.The objective was to enhance the accuracy of the WRF model's predictions for the Indian summer monsoon.However, it is important to note that the utilization of MO-ASMO came with the trade-off of a significantly larger sample size, which exceeded the number of dimensions in the parameter space by more than 40 times.Additionally, despite the optimization efforts, not all output objectives demonstrated improvement.
Research on parameter optimization for the NWP model has made some progress.However, there are several challenges that make MO not widely practiced in the NWP field.There are two primary challenges that need to be addressed: (a) MO: The NWP model involves numerous meteorological variables, and each parameter value needs to be optimized to achieve a Pareto optimal set that balances different objectives.This complexity increases the difficulty of parameter optimization and impacts the effectiveness of the optimization process.(b) Computational cost: Running an NWP model is computationally expensive, requiring significant time to complete long-term simulations over large geographical areas.Multiobjective optimization approaches often require a substantial number of model runs to obtain the Pareto-optimal set.This computational burden adds to the challenge of parameter optimization in the NWP model and affects the efficiency of the optimization procedure.
Therefore, the goal of this study is to develop a multiobjective method that can find the optimal solution for NWP model parameters while ensuring a reasonable optimization time.The remainder of this paper is organized as follows.In Section 2, the detailed methodology specified for the NWP model's MO is described.Section 3 provides a brief introduction to the experimental design.Section 4 presents the results and some discussions of MO.Finally, a brief conclusion is presented in Section 5.

Methodology
An NWP model can be represented by mathematical functions of numerous variables, including model input forcing, physical functions, model parameters, and model outputs.The formula is as follows: where y is the model output,   = (1, 2, . . ., )  is the model parameter vector, D is the number of model parameters, y 0 represents the initial states of the model, I is the model input variable, and f is the physical function.Model parameters are often estimated through MO.Model parameters can typically define subject to lower and upper bounds based on their physical meanings p (Lower) ≤ p ≤ p (Upper) , the MO process seeks to find the Pareto optimal parameter vector p (opt) with the parameter ranges that satisfied: where Y(p) is the objective function, Y(p) = y − y (obs) = f(y 0 ,I,p) − y (obs) + ϵ, y (obs) is the observation data corresponding to y, M is number of model outputs and ϵ is the model estimation error.For MO of the NWP model, a flowchart of the MO specified for the NWP model is shown in Figure S1 in Supporting Information S1.
First, the lower and upper bounds of the parameters, based on physical meanings and their probability distributions, need to be clarified.Here, a uniform distribution was assumed.A design of experiment (DoE) should be established before MO (Wu & Hamada, 2009), and model simulations should be performed based on this design.
The final result of this step is a sequence of design points p i = [p i1 ,p i2 ,⋯,p iD ],i = 1,2,⋯,m where m is the number of points in the design.The model under consideration was then evaluated for each design point, as previously described.This step creates a sequence of results y i = f(p i ) = f(p i1 ,p i2 ,⋯,p iD ),i = 1,2,⋯,m.Here, the Quasi Monte Carlo sampling, specifically the Sobol' sequence (proposed by Sobol' (1967)), was utilized.QMC sampling is applied by generating a deterministic sequence of points based on the Sobol' sequence.The Sobol' sequence is a low-discrepancy sequence that possesses uniform distribution properties and exhibits low discrepancy.It is generated using prime numbers as bases to generate sample points for each dimension.
Because the NWP model requires high computational resources, surrogate modeling should be used to construct a cheaper-to-run surrogate model of the WRF model.The Gaussian process (GP) model is a surrogate modeling method that uses prediction with noise-free responses (interpolation) or noisy responses (regression) to represent the original model (Gong et al., 2016;Jones et al., 1998;Knowles, 2006).A GP model is described by the following equation (Rasmussen & Williams, 2006;Williams & Rasmussen, 1996): the first term β T f(x) represents the mean value of the GP, which is also referred to as the trend.It is composed of regression coefficients β and arbitrary functions f.The second term consists of the variance of the GP, denoted as σ 2 , and Z(x,ω), which is a stationary GP with zero mean and unit variance.The underlying probability space is represented by ω and is defined in terms of a correlation function R (i.e., correlation family) and its hyperparameters θ.The correlation function R = R(x,x′,θ), in turn, describes the correlation between two sample points in the output space, that depends on x, x′, and the hyperparameters θ.
In this study, the surrogate model of the WRF was constructed using samples of input parameters and the mean values of the GP model output with noise-free responses.For GP, the Matérn-5/2 covariance kernel was considered with a constant yet unknown trend.The hyperparameters of the GP model, estimated through maximum likelihood estimation and optimized using a hybrid genetic algorithm, significantly impact the behavior and performance of the GP model.
For population generation of MO, a reference point adaptation-based evolutionary algorithm (AR-MOEA) was applied (Tian et al., 2018).The main idea of the AR-MOEA lies in the reference point adaptation method, which can adaptively adjust the distribution of the reference points according to the shape of the current population in the objective space, thus making the reference points capable of capturing different shapes of Pareto fronts.The AR-MOEA has promising versatility for problems with regular or irregular Pareto fronts.
The performance of a MO algorithm can be evaluated from the perspectives of convergence and diversity.Convergence refers to the minimization of the average distance between the realized Pareto optimal points and the real Pareto front, while diversity refers to the uniformity of the distribution matched with the real Pareto front (Gong et al., 2016).Diversity can be measured by crowding distance, which is defined as the sum of the edge lengths of the cuboid (Deb et al., 2002), as shown in Figure S2a in Supporting Information S1.The crowding distance directs the population toward unexplored regions, and this mechanism can effectively improve the diversity of the popula tion (Gong et al., 2016).A larger distance indicates better diversity.If the next iteration is selected in P 1 and P 2 in Figure S2a in Supporting Information S1 based on the crowding distance, P 1 will be selected and P 2 will be ignored.
For population selection (i.e., adaptive sampling in Figure S1 in Supporting Information S1), a knee point was suggested, as shown in Figure S2b in Supporting Information S1.We assume that four solutions are selected from the five nondominated solutions for the next population in Figure S2b in Supporting Information S1.If a crowding-distance-based criterion is chosen, then solutions P 0 , P 1 , P 3 , and P 4 are selected, and P 2 is ignored.However, P 2 can be considered as a knee point among the five nondominated solutions.When selecting solution P 2 , the knee point can be more beneficial than others in terms of the hypervolume.Therefore, in this study, we first select the knee point to prevent excellent solutions from being ignored, and then select solutions with larger crowding distances for adaptive sampling.
The knee point means the "concave point" of the nondominated fronts in the current population during the search process (Bechikh et al., 2011;Zhang et al., 2015).The basic idea is that knee points are naturally the most preferred among nondominated solutions if no explicit user preferences are given (Das, 1999).First, find M boundary point in the population, that is, the solution with the maximum objective value in each dimension, and then determine the hyperplane formed by these M points (line in 2 dimension, plane in 3 dimension).Then, calculate the vertical distance from each solution in the population to L. In this study, there are two objectives: precipitation and temperature, so the equation of L is ax + by + c = 0, so it is easy to know that the distance between solution P(x P ,y P ) and hyperplane L is For the minimization problem, the solution on the side closer to the origin is considered priority as the knee point, while the solution on the side away from the origin is not considered priority.Then the distance formula from P to L is changed to Finally, the solution P with the maximum distance d is chosen as the knee point.
Overall, we mainly use the knee point to promote the convergence of Pareto optimal solutions and the crowding distance to maintain the diversity of Pareto optimal solutions.Finally, the optimal solution set was updated using the NWP model simulation results obtained by adaptive sampling.The main loop is repeated until the termination conditions are satisfied such as performance indicators (e.g., hypervolume), iteration limit, sum of objective functions, and total original model evaluation limit.Here, the hypervolume was selected as the termination criterion.
The Pareto solution corresponding to the final knee point was an optimal solution for the NWP model parameters.This is the entire process of the KMO.

WRF Model Configuration of the Study Area and Weather Events Election
The Advanced Research WRF model (Skamarock et al., 2019) was used in this study as a representative of the NWP model.The Greater Beijing Area in North China was selected as our research area and was shown in Figure 1, which is located at 38.35°-42.25°Nand 113.35°-119.55°E.Here, a two-grid horizontally nested simulation area was designed.In the outer layer (d01 area), it was a horizontal resolution of 9 km.In the inner layer (d02 area), it was a horizontal resolution of 3 km.The 6-hr interval meteorological data from the Climate Forecast System Reanalysis, including surface and radiative flux data with a Gaussian T382 spatial resolution and three-dimensional pressure level data with 0.5° spatial resolution in latitude and longitude, were used to drive the WRF model as the initial and lateral boundary fields.Here, each simulated event spanned 2 days with additional 6-hr spin-up period.
Simulations of precipitation and temperature in the WRF model were considered in this study.The observed data were obtained from the China Meteorological Administration Land Data Assimilation System (CLDAS-V2.0),which includes a series of 1-hourly, 0.0625° meteorological element data (Shi et al., 2019).Based on the observed data, Figure S3 in Supporting Information S1 shows the daily accumulated precipitation for the Greater Beijing Area during the summer season from 2008 to 2010.Nine events were selected to obtain the WRF parameter optimal value for the precipitation and temperature simulations, and another six events were selected to validate the optimal parameter values.In order to increase the reliability of the results obtained by our method, a testing period of nine events from 2015 to 2017 was also considered (Figure S4 in Supporting Information S1).

Tunable Parameters Selection and Experimental Setup for KMO
Based on the operational setup of the Beijing Meteorological Bureau, the following physics parameterization schemes are applied: the Dudhia shortwave radiation scheme (Dudhia, 1989), the RRTM longwave radiation scheme (Mlawer et al., 1997), the Monin-Obukhov surface layer scheme (Dudhia et al., 1999), the unified Noah land-surface model (Chen & Dudhia, 2001), the Kain-Fritsch Eta cumulus scheme (Kain, 2004), the Yonsei University planetary boundary layer scheme (Hong et al., 2006), and the WSM six-class Graupel microphysics scheme (Hong & Lim, 2006).Among them, the Kain-Fritsch Eta cumulus scheme was only used for the outer-layer simulation with a spatial resolution of 9 km.Di et al. (2017) conducted a sensitivity analysis study of chosen 18 tunable parameters for precipitation and temperature in the WRF model.It identified the sensitive parameters (P 1-5 in Table S1 in Supporting Information S1).As a follow-up study of Di et al. (2017), we optimize the most sensitive parameters to improve the precipitation and temperature to test the effectiveness and efficiency of the newly developed method.For the MO, the objective function for precipitation is the normalized root-mean-square error (NRMSE), while the objective function for temperature is also NRMSE.
where sim i,t , obs i,t , and def i,t are the simulated, observed and default values of a variable at the ith grid cell and at time t, respectively, that is, numerator is the RMSE of the particular sample and denominator is the RMSE obtained from the default parameter run.N is the number of grid cells in domain d02, and T is the total number of forecasted time intervals.M is the number of simulation events.The number of initial samples used to construct the GP model was set to 10 times the number of sensitive parameters; that is, 5 × 10 = 50.For the AR-MOEA setup, the population size was set to 100, and the number of generations was set to 1,000; that is, the maximum number of evaluations for the AR-MOEA was 100 × 1,000 = 100,000.For adaptive sampling, the sample size in each iteration was 5 (i.e., one knee point and four Pareto solutions with the largest crowding distances), and the maximum number of iterations for MO was set to 20; that is, the maximum number of adaptive WRF evaluations was 5 × 20 = 100.Overall, the total WRF model evaluation is the number of initial samples plus the number of adaptive sampling in iterations, that is, 50 + 100 = 150.

Multiobjective Optimization for Precipitation and Temperature
We compared the KMO algorithm with another popular method (Gong et al., 2016), MO-ASMO, and discovered that our algorithm can uncover superior solutions.The settings for the MO-ASMO algorithm are referenced from Gong et al. (2016).The convergence trend for the MO algorithms is presented in Figure S5 in Supporting Information S1.We found that the number of iterations required for the KMO's convergence was 15, that is, a set of optimal solutions was obtained after 15 iterations.As the initial sample set was 50 samples, 125 (i.e., 50 + 15 × 5 = 125) WRF model evaluations were performed to find the final optimal solution.
The RMSE values of the different calibration events with the default and optimal parameter set were shown in Figure S6 in Supporting Information S1.It could be seen that there was a trend of improvement in both precipitation and temperature.In particular, for event e, precipitation improved by about 32%.Although there were a few events in which the RMSE values increased slightly, the overall trend showed that the optimal parameter set for precipitation and temperature performed well.Seventy percent of the precipitation events and 90% of the temperature events improved.It could be seen that the average percentage improvements of RMSE values for precipitation and temperature were 4.81% and 18.15%, respectively.Additionally, we have incorporated additional evaluation metric (i.e., frequency distribution of bias differences between the optimal and default simulations of precipitation/temperature) to provide a more comprehensive assessment.The results are shown in Figure S7 in Supporting Information S1.When the optimal parameter values are substituted for the original parameter values, approximately 52.5% of the total grids experience a decrease in precipitation error, while around 61.8% of the grids observe a decrease in temperature error.
Figures 2 and 3 shown the comparison of spatial patterns of the daily average precipitation and temperature for all calibration events with the default and optimal parameter set, respectively.Figure 2a represents the observed precipitation, and Figures 2b and 2c represent the default and optimal spatial patterns, respectively.Figures 2d and 2e show the spatial patterns of the bias between observation and simulation.The observed precipitation indicated it had heavy rainfall over the eastern coastline of the Greater Beijing Area, whereas the default parameter set failed to capture some parts of the coastline.The optimal parameter set predicted the heavy rainfall of the east coastline slightly better than the default parameter set, that is, the optimized length of the rainfall coastline was closer to the observed data, as also seen in the difference's spatial patterns.For the temperature simulations, it could be seen from Figures 3a-3c that WRF simulations using the optimal parameter set have significant improvements compared to the default parameter set, for instance, the southernmost central high-temperature region predicted by the optimized parameters is closer to the observation data, as also seen in the difference's spatial patterns in Figures 3d and 3e.
During calibration, the superiority of the optimal parameter set obtained by KMO for improving the WRF precipitation and temperature simulation was demonstrated compared with the default parameter set.However, for the other periods, the question of whether the optimal parameters are still effective needs to be further studied.During validation, Figure S8 in Supporting Information S1 shows the RMSE values and Figure S9 in Supporting Information S1 displays the bias differences between the optimal and default simulations.The average percentage improvements in RMSE for precipitation and temperature are 18.05% and 24.5%, respectively.Approximately 58.3% of grid points exhibit a decrease in precipitation error, while around 81.9% show a decrease in temperature error.Figures S10 and S11 in Supporting Information S1 compare spatial patterns of daily average precipitation and temperature for all validation events using the default and optimal parameter sets.
During testing, Figure S12 in Supporting Information S1 presents the RMSE values and Figure S13 in Supporting Information S1 shows the bias differences between the optimal and default simulations.The average percentage improvements in RMSE for precipitation and temperature are 3.76% and 6.67%, respectively.Approximately 61.7% of grid points exhibit a decrease in precipitation error, while around 66.3% show a decrease in temperature error.Figures S14 and S15 in Supporting Information S1 compare spatial patterns of daily average precipitation and temperature for all testing events using the default and optimal parameter sets.
Overall, the KMO method improved the simulation accuracy of precipitation and temperature, particularly for temperature.

Physical Interpretation of the Parameters' Optimal Value
The default and optimal parameter values are listed in Table S2 in Supporting Information S1.Compared to the observations, WRF simulations with default parameters showed an overestimation of both precipitation and temperature (Figures 2 and 3), particularly in regions characterized by heavy rainfall (e.g., maritime region) and high temperatures (e.g., the southernmost central high-temperature region).However, after optimization, the WRF simulations using the optimized parameters exhibited improved agreement with the observations, especially in the regions of maritime heavy rainfall and the southernmost central high-temperature region.By examining the tail behavior of the relative precipitation error frequency histograms (e.g., Figure S7 in Supporting Information S1), improvements in extreme values due to optimization can also be observed.In summary, compared to simulations using default parameters, WRF simulations with optimized parameters resulted in reduced predictions of both precipitation and temperature.In the microphysics scheme, scaling factor applied to ice fall velocity P1 was optimized from the original values of 14,900 to 8,000 s −1 .Limiting maximum value for the cloud-ice diameter P2 was optimized from the original values of 5•10 −4 to 7•10 −4 m.A reduction in the falling speed of ice crystals with larger diameters weakens the conversion of cloud ice into rainwater, potentially leading to a decrease in precipitation.
In the short-wave scheme, scattering tuning parameter P3 was optimized from the original values of 1•10 −5 to 2•10 −5 m 2 kg −1 .The higher the value of P3, the more scattered solar radiation in the atmospheric layer, which reduces the amount of shortwave radiation that reaches the surface and evaporates from it, potentially leading to a decrease in precipitation.Simultaneously, the amount of shortwave radiation that reaches the surface is relatively low, which may also lead to a decrease in temperature.
In the land-surface scheme, multiplier for the saturated soil water content P4 was optimized from the original values of 1 to 0.5080.The lower the value of P4, the smaller the soil porosity between the land surface and the groundwater surface, which leads to weaker upward transport of soil moisture, thereby reducing surface evapotranspiration.This affects the development of rainfall.At the same time, a lower porosity also leads to a higher soil water content and wetter soil, which are not helpful for the development of temperature.This may have led to a decrease in precipitation and temperature.
In the planetary boundary-layer scheme, profile shape exponent used to calculate the momentum diffusivity coefficient P5 was optimized from its original value of 2 to 1.4107.Because P5 is the profile shape exponent used to calculate the diffusivity coefficient for momentum (Hong & Pan, 1996;Noh et al., 2003).
where k is the von Kármán constant, ω s is the mixed-layer velocity scale, z is the height from the surface, and h is the height of the planetary boundary-layer.Given the same vertical shear and surface heating (i.e., ω s ), as the parameter is smaller, turbulent mixing is weaker (Equation 7) and is confined to lower levels, suppressing the development of the planetary boundary layer (thus leading to a lower planetary boundary layer height).This may suppress convection, as indicated by decreased precipitation.

Conclusions
The study proposed the KMO method to improve summer NWP in the Greater Beijing Area, China.The results showed that it required fewer than 125 samples (i.e., 25 times the number of dimensions of the parameter space) to obtain the WRF model's optimal parameter values.The calibrated WRF model, using these optimal parameters from simulations conducted during calibration period, showed enhanced precipitation and temperature results compared to default simulations, across all periods (calibration, validation, and testing).This improvement was validated through metrics such as RMSE and bias differences.Therefore, the optimized parameters were considered suitable for precipitation and temperature simulations in the WRF model.
The knee point on the Pareto front represents a locally maximal marginal utility.It signifies a "bulge" on the Pareto front toward the direction of the ideal point, indicating its proximity to the origin (representing a solution where simulations perfectly match observations) compared to neighboring solutions.The KMO method holds significant implications for NWP, enhancing accuracy and reliability by optimizing multiple objectives in weather modeling.These findings have practical applications in sectors such as agriculture, transportation, and disaster management, where precise and timely weather information is crucial.
However, the KMO method encounters challenges in effectively modeling dynamic input-output maps.When applied to NWP models across extensive spatial and temporal domains, it is unlikely to generate consistent spatial patterns and dynamic relationships between model inputs and outputs.As a result, the model's performance may become unbalanced across different spatiotemporal domains.This issue requires further attention and investigation in future research.
Furthermore, it is worth noting that the presence of structural uncertainties in the WRF model and variations in model structure, particularly in microphysical processes, have the potential to affect the accuracy of precipitation and temperature simulations in the WRF model.It is recommended that future research explores the evaluation of physics parameterization schemes in different regions before conducting parameter optimization.This integrated approach of model structure selection along with parameter optimization methods holds promise for achieving improved results.

Figure 1 .
Figure 1.The two-level nested study domain for the Weather Research and Forecasting model, with the outer layer (d01) being North China and the inner layer (d02) being the Greater Beijing Area.

Figure 2 .
Figure 2. Spatial distributions of daily average precipitation simulation: (a) observation of precipitation; (b) Weather Research and Forecasting (WRF) simulation of precipitation using default parameter set; (c) WRF simulation of precipitation using optimal parameter set; (d) difference between WRF simulations using default parameter set and observation; and (e) difference between WRF simulations using optimal parameter set and observation.The average field values are displayed after the subfigure titles.

Figure 3 .
Figure 3. Same as in Figure 2, but for temperature.