2.1. Study zone

The state of Puebla is located in the central part of Mexico. It borders to the north with the states of Hidalgo and Veracruz; to the east also with Veracruz and Oaxaca; to the south with the latter and Guerrero, and to the west with this state, Morelos, Mexico, Tlaxcala and Hidalgo (^{Tamayo, 1996}). It has an area of 34 290 km², which represents 17% of the total national space. It is characterized by a wide topographical heterogeneity because it houses four major biogeographical provinces: The Sierra Madre Oriental, the coastal plain of the North Gulf, the Neovolcanic Axis, and the Sierra Madre del Sur. This geomorphological diversity causes marked changes in altitude, which give rise to a wide variety of climates, dominating the temperate climates that cover most of the territory, followed by warm and semi-arid climates. The climatic heterogeneity is due, in part, to the fact that as altitude increases, temperature decreases and 65% of the territory of Puebla is composed of mountainous topography and hills.

2.2. Description of the WRF model

The performance of numerical models of weather prediction has increased in the last 40 years due mainly to four factors (^{Kalnay, 2003}): (1) The increase in computing power of supercomputers, allowing a finer numerical resolution and fewer approximations in operational atmospheric models; (2) the improvement of the representation of small-scale physical processes within the models (clouds, precipitation, turbulent heat transfer, humidity, momentum and radiation); (3) the use of more accurate data assimilation methods, which results in an improvement of the initial conditions for the model, and (4) the increase in data availability, especially satellite and aircraft data on the oceans and the southern hemisphere.

The main features of the WRF model revolve around its non-hydrostatic dynamics and its ability to process spatial resolutions of a few kilometers (^{Moya-Álvarez and Ortega-León, 2015}). About its structure, the WRF model has two dynamic cores (Advanced Research WRF [ARW] and Nonhydrostatic Mesoscale Mode [NMM]), a data assimilation system, and a software architecture that allows the application of parallel computing to perform simulations (^{Skamarock et al., 2008}). In this study, the dynamic core WRF-ARW version 3.9.1.1 was used. The WRF modeling system requires external data sources and consists of three main modules: (1) the WRF preprocessing system (WPS); (2) the ARW solver, and (3) third-party postprocessing and visualization tools.

2.3. Dataset

The forecast accuracy of the models lies primarily in a good description of the initial state of the atmosphere. This initial state, analysis or first approximation is created by an optimal combination between observed data and a short-term forecast derived from a previous analysis through a process known as data assimilation. According to ^{Zepka (2011)}, for satisfactory results regarding the predictability of a storm or any adverse weather phenomenon characterized by very small spatial and time dimensions, high-quality input data with high temporal and spatial resolutions are necessary, as well as a high-resolution model.

The WRF model requires knowledge of the initial and boundary conditions for the simulation period at constant time intervals, information which is usually incorporated from global data. For this study, we chose to use the data from the Global Forecast System (GFS) model, which is a reference for operational forecasting and research studies. GFS data was downloaded from ftp://ftp.ncep.noaa.gov/pub/data/nccf/com/gfs/prod.

Daily rainfall records were obtained from the Hydrological Information System (HIS) developed in the Surface Water and River Engineering Management of CONAGUA, through the site http://148.204.8.145. The HIS is a system sponsored by the World Meteorological Organization (WMO) through the Water Management Modernization Program.

Fifty-four weather stations with at least 80% data in the study period were selected as sampling units. Table I lists the 54 selected weather stations with their corresponding keys, and Figure 1 their geographical location within the state of Puebla.

Fig. 1 Geographical location of the selected weather stations in the state of Puebla (right). The station number corresponds to the list of stations in Table I. Experimentation domains with the WRF model (D1 [mother domain, red] and D2 [nested domain, orange]) are marked on the map of Mexico (left). 

2.4. Experiment design

Two domains were defined for the execution of the WRF model, the mother domain with a spatial resolution of 16 km and a nested domain with a resolution of 8 km. Only simulated precipitation in the nested domain was used to evaluate the performance of the WRF model in the state of Puebla. Table II shows the fixed parameters selected for experimentation with the WRF model.

As mentioned above, 768 experiments with the WRF model were executed, each one corresponding to one of the possible combinations of the different physical settings selected. Each parameterization was varied in two options, except for the convection parameterization, from which six different ones were selected since precipitation is directly associated with the formation of clouds. However, all the different physical parameterizations affect the simulation of precipitation. The combination was selected on a discretionary basis of the compatibility at the run time between different combinations. Table III lists the physical parameters considered in the 768 experiments.

The results of the experiments were sorted by date and weather station, generating a text file consisting of 81 lines which corresponds to the analysis days of the period from June 1 to August 20, 2017. Each line contains 770 columns, the first one containing the date; the second column contains the observed precipitation value, and columns 3 through 770 contain the precipitation value simulated by the WRF model in each of the 768 experiments.

2.5. Metrics to measure the performance of the WRF model

In systems modeling, an essential stage that presents both conceptual and practical difficulties is the validation of the models. An important part of this process is empirical validation, which according to ^{Reynolds (1984)} and ^{Mitchell (1997)} is done to compare the predictions of the model with observations from the real world. According to ^{Aguilar (1997)} and ^{Rauscher et al. (2000)}, these comparisons should ideally be carried out using appropriate statistical methods, with an acceptable level of confidence, so that the inferences are correct (^{Barrales et al., 2004}).

There are different methods to quantitatively validate numerical forecast models, highlighting the use of simple statistics of bias, root mean square error (RMSE) and Pearson correlation. These parameters were selected because they are not exclusive and can be used together. Among the works that propose and use these statistical parameters for the evaluation of forecast models are ^{Willmott (1982)}, ^{Pielke (1984)}, ^{Willmott et al. (1985)}, ^{Carbonell et al. (2003)} and ^{Das et al. (2015)}.

The selected statistics were applied to each observation site for both the set of observed data (O) and model predictions (P). To the extent that the statistical indicators are favorable and show that the simulated data are approximate to observed data, and that the behavior over time of simulated variables is similar to that observed, it can be concluded that the simulation provides representative data and that the WRF model can simulate the precipitation in the state of Puebla (^{Gavidia, 2012}).

RMSE is a measure of quantitative performance commonly used to evaluate forecasting methods. In this context, RMSE consists in the square root of the sum of the quadratic errors, which captures positive errors as negative; therefore, it expresses both systematic and random errors. RMSE amplifies and penalizes with greater force those errors of greater magnitude (Eq. 1).

where P _i is the prediction of the model in position i, O _i is the value observed in position i and n is the sample size.

The bias is a consistent, constant, and one-way error. It can be positive (overestimation) or negative (underestimation), which does not depend on the number of elements but on many other factors that can be controlled in general if their existence is suspected and the research is carefully planned and executed (^{Moya de Madrigal, 2005}). This statistic measures the reliability of the model, revealing the systematic error. Not having a bias or having a minimal one is a desirable property which indicates that the forecast is close to reality (Eq. 2).

where P _i is the prediction of the model in position i, O _i is the value observed in position i, and n is the sample size.

Pearson correlation, denoted by the letter r, is a normalized measure of the linear relationship between two continuous quantitative variables, that is, it measures the dependence of one variable with respect to another independent variable. The Pearson correlation coefficient allows establishing similarities or dissimilarities between variables, to make evident the joint variability and therefore typify what happens with the data (^{Mondragón-Barrera, 2014} ). The coefficient can score values ranging from -1.0 to 1.0. Values close to 1.0 indicate that there is a strong association between the variables, that is, when one increases the other as well. On the other hand, values close to -1.0 indicate that there is a strong negative association between the variables, implying that when one variable increases, the other decreases. When the value is 0.0, it indicates that there is no correlation, or the correlation is null (^{Anderson et al., 2008}). In the correlation, the dependent variable is not distinguished from the independent one, so the correlation of O with respect to P is the same as the correlation of P with respect to O (Eq. 3).

where O _i is the value observed in position i, P _i is the prediction of the model in position i, P̅ is the mean observed value, and P̅ is the mean model prediction.

The efficiency multiparameter index (EMI) was used to identify exposed cases with better performance at the state level, which involves four steps:

where EMI _e is the efficiency multiparameter index for experiment e, ne the number of experiments, nl the number of locations, nm the number of metrics, and E _eml the weight of the experiment e for the metric m and the location l.

3.1. Precipitation variability during the study

Given that the year 2017 presented a neutral ENSO condition, rainfall in the state of Puebla during the study period was normal with the exception of August 10, 2017 when in the northern and eastern regions of the state, precipitation from very strong to intense occurred with values exceeding 75 mm accumulated in 24 h. This extraordinary precipitation was caused by hurricane Franklin (category 1), which hit land on August 10 at 00:00 LT in the state of Veracruz and subsequently weakened at 4:00 LT to a tropical storm in the northern highlands of Puebla.

3.2. Analysis and interpretation of results

A first analysis was performed with the results of the statistical evaluation parameters in the 54 locations. Table IV shows the minimum and maximum values in the 768 experiments by weather station and statistical metrics.

Overall, the 768 experiments tend to underestimate precipitation with negative biases in 33 observation sites, positive biases with a tendency to overestimate precipitation in only six sites, and a mixed trend in 15 locations. Extreme biases in the 54 locations were presented in the interval [-41.96, 18.06].

On the other hand, the different experiments with the WRF model gave RMSEs in the interval [3.14, 69.73]. The lowest RMSE was obtained at ACHPB, AHUPB, CNGPB, CDSPB, TEHPB and TPYPB stations with intervals [3.18, 13.52], [4.71, 20.51], [5.0, 26.36], [4.65, 19.66], [4.02, 18.83], and [3.14, 15.92], respectively. At the opposite end, the WRF model presented the highest RMSE at the HJTPB and TLPPB stations with intervals [15.53, 55.27] and [8.25, 69.73], respectively.

Regarding the Pearson correlation parameter, the results are very varied in the same location. In the 768 experiments both negative and positive correlations were observed, the latter predominating. Under this metric, the WRF model achieved extremely high positive correlations with values greater than 0.9 at the AHUPB, CDSPB, OYAPB, TECPB, TEOPB, TPYPB, ZCPPB, and ZRGPB observation sites. Meanwhile, through the different 768 configurations, the lowest correlations were obtained at stations CTZPB, TEZPB and ZOQPB stations with intervals of [-0.04, 0.27], [-0.08, 0.06] and [-0.04, 0.27], respectively.

A second analysis of results was performed by grouping the observation sites into six groups with ranges of 500 masl, as shown in Table V. It was observed that at altitudes above 2000 masl, corresponding to groups five and six (which together concentrate 50% of the weather stations), the WRF model has better performance, especially in the central and northern region of the state of Puebla, with the exception of the ZOQPB station, where moderate performance was observed. The region in which the WRF model generally presented the lowest performance was the southwest of the state. The contrast between the 768 experimental groups in the same observation site is evident, starting with performances from low to very high. These results emphasize the importance of conducting experiments to determine the appropriate configuration according to the study area. The exception is the CNGPB station, located to the north of the state of Puebla at an altitude of 2300 masl, where the WRF model in all experiments presented medium to very high performance, with bias in the interval [-11.85, 1.61], an RMSE maximum of 26.36 and r greater than 0.29; however, a bad configuration of the WRF model physical parameters in this station would not present, in statistical terms, serious errors in the simulation of precipitation.

3.3. Optimal configuration of the WRF model

Three optimal configurations were determined to execute the WRF model by location, each one corresponding to the experiment which obtained the better performance in the three evaluation metrics. Tables VI, VII and VIII describe the optimal configurations according to the statistical parameters of bias, RMSE, and Pearson correlation, respectively.

The results show that no single experiment has the best performance for the WRF model in all the observation sites, nor in the three statistical metrics. Although some experiments with high performance in different metrics and observation sites were identified, they were scarce.

The experiment with control number 571 obtained the best performance at the state level according to the EMI, with a very good value of 0.76. Table IX shows the 10 experiments with highest EMI values. After analyzing the 10 best configurations according to the EMI (Table X), it was identified that the settings of the microphysics WRF Single-Moment (WSM) 5-class scheme (option 4), the convection parameterization Grell-Devenyi ensemble scheme (option 3), the planetary boundary layer parameterization Mellor-Yamada-Janjic (Eta) TKE scheme (option 2), and the long wave radiation Goddard Shortwave Scheme (option 2) are kept constant to achieve a comparable performance at the state level.

Figures 2 3, 4, 5, 6, 7 and 8 present the time series corresponding to observed precipitation in comparison to simulated precipitation by the WRF model, which were analyzed in four experiments corresponding to the following parameters: (1) less BIAS, (2) lower RMSE, (3) greater Pearson correlation, and (4) higher statewide EMI. The following four observations were derived from this time series analysis:

Fig. 2 Time series for the ACHPB, ACJPB, ADOPB, AFRPB, AHAPB, AHUPB, APAPB, and ATXPB stations, showing the comparison between observed precipitation and precipitation simulated by the WRF model in the four experiments with better performance in each locality. 

Fig. 3 Time series for the AVCPB, CDSPB, CEMPB, CHIPPB, CHLPB, CHSPB, CNAPB, and CNGPB stations, showing the comparison between observed precipitation and precipitation simulated by the WRF model in the four experiments with better performance in each locality. 

Fig. 4 Time series for the COAPB, CPLPB, CTZPB, ECHPB, ELCPB, HJTPB, HQCPB, and HUAPB stations, showing the comparison between observed precipitation and precipitation simulated by the WRF model in the four experiments with better performance in each locality. 

Fig. 5 Time series for the IZMPB, LIBPB, MAYPB, NNEPB, OYAPB, PCRPB, PIXPB, and PTLPB stations, showing the comparison between observed precipitation and precipitation simulated by the WRF model in the four experiments with better performance in each locality. 

Fig. 6 Time series for the PUOPB, QUIPB, RNAPB, SCRPB, SMCPB, SOLPB, TDRPB, and TECPB stations, showing the comparison between observed precipitation and precipitation simulated by the WRF model in the four experiments with better performance in each locality. 

Fig. 7 Time series for the TEHPB, TEOPB, TEPPB, TEZPB, TLPPB, TLXPB, TPYPB, and TYAPB stations, showing the comparison between observed precipitation and precipitation simulated by the WRF model in the four experiments with better performance in each locality. 

Fig. 8 Time series for the UDSPB, VENPB, XDJPB, ZCPPB, ZOQPB, and ZRGPB stations, showing the comparison between observed precipitation and precipitation simulated by the WRF model in the four experiments with better performance in each locality. 

During the study period, the state of Puebla presented a strong rainfall concentration (greater than 25 mm per day) from June 20 to July 10. Intense to extraordinary rainfall (greater than 75 mm per day) also took place at 23 observation sites in the northern part of the state due to Hurricane Franklin on August 10, exceeding the highest historical rainfall at SOLPB, ZCPPB, and ZRGPB stations with records of 225, 281 and 198 mm, respectively.

It was observed that in the CTZPB station, the observed data of extraordinary rainfall occurred on August 10 show a two-days forward displacement, while at the TEZPB station data show a one-day backward change . Therefore, there might be an error in the capture dates at these two observation sites, which would explain the poor performance of the WRF model at these locations.

Although the three best experiments identified by location (corresponding to the lower bias and RMSE, and the higher Pearson correlation) show an outstanding performance of the WRF model in the precipitation simulation, they do not display the same performance in a homogeneous way in the 54 observation sites.

The performance of experiment 571, with the highest EMI at the state level in the 54 observation sites, is comparable to that obtained in the experiments that presented better performance locally for the bias, RMSE and r parameters.

Finally, Figure 9 shows the spatial distribution of the observed and simulated precipitation in experiment 571, together with its corresponding anomaly for August 9, 10 and 11, 2017. On day 10, the highest rainfall of the study period occurred in the northern region of the state of Puebla, associated with Hurricane Franklin. It was observed that the WRF model had the ability to simulate the intensity of precipitation during the extreme event; however, rainfall patterns show variations in detail. These variations are mainly associated with the fact that the 8-km resolution of the model in the nested domain was insufficient to characterize the terrain in the complex topography of the state.

Fig. 9 Spatial distribution of rainfall on August 9-11, 2017. Left: observed rainfall; center: rainfall simulated with the WRF model in experiment 571; right: anomaly between observed rainfall and rainfall simulated in experiment 571.