This study investigates the sensitivity of daily rainfall rates in regional seasonal simulations over the contiguous United States (CONUS) to different cumulus parameterization schemes. Daily rainfall fields were simulated at 24-km resolution using the NASA-Unified Weather Research and Forecasting (NU-WRF) Model for June–August 2000. Four cumulus parameterization schemes and two options for shallow cumulus components in a specific scheme were tested. The spread in the domain-mean rainfall rates across the parameterization schemes was generally consistent between the entire CONUS and most subregions. The selection of the shallow cumulus component in a specific scheme had more impact than that of the four cumulus parameterization schemes. Regional variability in the performance of each scheme was assessed by calculating optimally weighted ensembles that minimize full root-mean-square errors against reference datasets. The spatial pattern of the seasonally averaged rainfall was insensitive to the selection of cumulus parameterization over mountainous regions because of the topographical pattern constraint, so that the simulation errors were mostly attributed to the overall bias there. In contrast, the spatial patterns over the Great Plains regions as well as the temporal variation over most parts of the CONUS were relatively sensitive to cumulus parameterization selection. Overall, adopting a single simulation result was preferable to generating a better ensemble for the seasonally averaged daily rainfall simulation, as long as their overall biases had the same positive or negative sign. However, an ensemble of multiple simulation results was more effective in reducing errors in the case of also considering temporal variation.
Atmospheric convection associated with cloud generation plays an important role in the global energy balance, hydrologic cycle, and ocean–land surface interactions with the atmosphere (Arakawa 2004). Small-scale convection with cumulus clouds yields large uncertainties in atmospheric model simulations in terms of predicting the magnitude and timing of the convection, because the spatial scales of such convection are less than the grid intervals used in most climate model configurations. The influence of the subgrid convection is therefore parameterized in the calculation of gridpoint prognostic variables. The numerical representation of subgrid cumulus convection is referred to as cumulus parameterization. This type of parameterization is a mutual complement to large-scale condensation or explicit cloud microphysics in atmospheric modeling frameworks.
To date, many cumulus parameterizations have been designed and developed in the framework of global climate models as well as limited-area models. These parameterizations are characterized by different assumptions in terms of, for example, supplementary assumptions for the closure problem, the type of equation describing advection, the parameterized exchange between clouds and environment, etc. These parameterizations significantly influence the evolution of various prognostic variables by modifying atmospheric structure. Therefore, it is important to evaluate the sensitivity of model simulations to the selection of the parameterization through the comparison with observed characteristics at various temporal and spatial scales.
The discrepancy and spread exhibited across simulations as a result of different cumulus parameterizations are recognized as major factors contributing to uncertainty in regional seasonal/climate simulations, along with other physics parameterizations, for example, cloud microphysics, planetary boundary layer, land surface, and radiative transfer schemes (e.g., Liu et al. 2011). The uncertainty due to cumulus parameterizations has been demonstrated by previous studies and coordinated modeling projects. For example, two types of cumulus parameterizations were tested within the framework of the North American Regional Climate Change Assessment Program (NARCCAP; Mearns et al. 2009). Multiple cumulus parameterization schemes were compared to evaluate their ability to simulate the evolution of the North American monsoon system (Gochis et al. 2002). An optimal ensemble approach to incorporate the benefits of multiple simulations with different parameterizations was suggested to improve the skill of precipitation simulation (Liang et al. 2007). Twelve cumulus parameterization schemes were evaluated to explore their skill in predicting the distribution of seasonally averaged rainfall, the frequency of daily rainfall intensity, and the precipitation diurnal variation for a couple of flooding periods over the central United States (Qiao and Liang 2015).
The present study investigates the sensitivity of simulated summer rainfall over the contiguous United States (CONUS) to the choice of cumulus parameterization scheme. The simulations were conducted using the NASA-Unified Weather Research and Forecasting (NU-WRF) Model (Peters-Lidard et al. 2015). The simulation members were prepared by selecting different deep and shallow cumulus parameterization schemes while holding all other model configuration components constant.
2. Description of the NU-WRF CONUS seasonal simulations
NU-WRF is a modeling system integrating the National Center for Atmospheric Research (NCAR) Advanced Research version of Weather Research and Forecasting (WRF) Model (WRF-ARW) (Skamarock et al. 2008) with multiple modeling components developed by the NASA Goddard Space Flight Center (GSFC). The present study employed a special version of NU-WRF, which was based on WRF-ARW, version 3.5.1, and included a bug fix that removed the accumulation of round-off errors in lateral boundary conditions for better long-term simulations (Dudhia 2015). This version of NU-WRF was then used in NASA’s downscaling simulation project (Ferraro et al. 2017), which consisted of regional simulations over CONUS with horizontal grid spacing of 24, 12, and 4 km. The cumulus scheme selection was varied for the 24-km suite of simulations, which was the sole resolution examined in the present study. The simulations extended from 1 November 1999 through 1 September 2000.
The initial and lateral boundary conditions were calculated in 6-hourly input intervals using the Modern-Era Retrospective Analysis for Research and Applications, version 2 (MERRA2; Bosilovich et al. 2015), with a horizontal resolution of 0.5° for latitude and 0.625° for longitude. Spectral nudging (Miguez-Macho et al. 2004) constrained by MERRA2 was applied to the calculation of the horizontal wind velocities, temperature, and geopotential heights. The nudging featured a relaxation time of 10 000 s and wavelengths of approximately 600 km for both east–west and north–south directions for the components above the boundary layer. NASA’s downscaling simulation project (Ferraro et al. 2017) was aimed at decadal simulations with fine horizontal resolution, so that the strong spectral nudging was configured in the control runs to improve the reproducibility of climatic phenomena as a counterpart of simulations without nudging. Note that the strong nudging possibly obscures the original strength of the sensitivity to the choice of cumulus parameterization scheme by inhibiting feedback effects from the physical process to the model dynamics. In turn, this study avoids direct comparison of nudged and nonnudged simulations, the latter of which can include the thermodynamics feedback with its intrinsic (modeled) strength.
The following configurations for model physics parameterizations were employed in the NU-WRF CONUS simulations. The Goddard cumulus ensemble (GCE) single-moment, 3-ice bulk microphysical scheme (e.g., Tao et al. 2003; Lang et al. 2007) was used as the grid-scale cloud microphysics scheme. Goddard shortwave and longwave radiation schemes (Chou and Suarez 1999; Chou et al. 2001) calculated the atmospheric radiation processes. For the planetary boundary layer and subgrid-scale turbulence, the level 2.5 Mellor–Yamada–Janjić turbulence scheme (MYJ; Janjić 1990, 1996, 2002) was chosen, as this scheme had a long reliable history in the National Weather Service operational models (e.g., Eta and North American Mesoscale Forecast System models). The corresponding Monin–Obukhov–Janjić Eta surface scheme (Janjić 1996, 2002) was required when running the MYJ turbulence scheme. The Noah land surface model (Chen and Dudhia 2001; Ek et al. 2003) within the NASA Land Information System (LIS) model (e.g., Kumar et al. 2006) was used to run multiyear offline spinups of the land surface model prior to coupled NU-WRF initialization to improve coarsely resolved initial soil conditions obtained from reanalysis data alone. The Noah land surface model was also employed to calculate the land surface states within NU-WRF after the initialization.
The simulation domain was discretized by horizontal grid points of 332 × 157 with grid intervals of 24 km. The atmosphere extended to a top pressure of 10 hPa and was divided into 41 vertical layers. Time step intervals of 72 s were employed for the integration of processes except for radiation, which was calculated at 24-min intervals.
Table 1 lists the members of the regional seasonal simulations characterized by different cumulus parameterization schemes. Four deep cumulus parameterizations were employed: the Grell 3D ensemble scheme (G3D; Grell and Dévényi 2002; Grell and Freitas 2014), the Betts–Miller–Janjić scheme (BMJ; Janjić 1994, 2000), the new Kain–Fritsch scheme (NKF; Kain 2004), and the new simplified Arakawa–Schubert scheme (NSAS; Han and Pan 2011). Each deep cumulus scheme includes a companion parameterization for shallow cumulus convection, and only G3D has the option of disabling its native shallow convection component. Among the three members using G3D in Table 1, GO used G3D’s native shallow convection component, GW included the University of Washington shallow cumulus parameterization scheme (UWSC; Bretherton et al. 2004; Bretherton and Park 2009) in place of its native component for shallow convection, and G did not employ any shallow convection component. The uses of the native shallow convection schemes and the UWSC are indicated by the capital letters O and W, respectively, in Table 1.
a. Spatial and frequency distributions of rainfall
Figure 1a illustrates the horizontal distribution of the daily surface rainfall averaged for June–August 2000 (JJA2000) in the Parameter-Elevation Regressions on Independent Slopes Model (PRISM) high-resolution spatial climate data (Daly et al. 2008). The original PRISM daily rainfall data gridded at resolution of 2.5 arc min (approximately 0.04°) were upscaled for horizontal grids of 0.25° (approximately 24 km) using bilinear interpolation for comparison to the NU-WRF simulations. The PRISM data in Fig. 1a show that JJA2000-average rainfall rates higher than 1 mm day−1 are observed mostly over the eastern half of CONUS. Rainfall rates exceeding 5 mm day−1 are observed over parts of the Appalachian Mountains and around Indiana, Illinois, Wisconsin, Iowa, and Minnesota. The PRISM rainfall distribution averaged for JJA2000 in Fig. 1a is similar to the climatology for JJA derived from the 30-yr (1981–2010) PRISM data (NCAR 2015), except for the higher rainfall rates over Florida and the southeastern coast in the 30-yr climatology.
Figures 1b and 1c illustrate the daily rainfall field of MERRA2 averaged for JJA2000 and its difference from PRISM (Fig. 1a), respectively. MERRA2 overestimates the PRISM rainfall over the Atlantic and Gulf Coast, Great Lakes regions, and around the southern edge of the Rocky Mountains. In contrast, it exhibits underestimation around Iowa and Missouri and parts of the Appalachian Mountains. The MERRA2 precipitation data over CONUS are derived originally from the background estimates of the Goddard Earth Observing System, version 5 (GEOS-5; Molod et al. 2012), and corrected using the data of National Oceanic and Atmospheric Administration (NOAA) Climate Prediction Center (CPC) Unified Gauge-Based Analysis of Global Daily Precipitation (CPCU; Xie et al. 2007). The details of the correction algorithm can be found in Reichle and Liu (2014).
Figure 2 shows corresponding daily rainfall averaged for JJA2000 in NU-WRF simulations using the different cumulus parameterization options listed in Table 1. Figure S1 in the supplemental material illustrates their differences from that in PRISM (Fig. 1a). The simulations over- or underestimate the rainfall derived from PRISM according to the cumulus parameterization schemes selected. The JJA2000-mean daily rainfall averaged over the entire CONUS is listed in Table 1. The relative magnitudes are GO > KO > G > MERRA2 > SO > BO > PRISM > GW (the definitions of the abbreviations are shown in Table 1). Schemes G, BO, KO, and SO exhibit relatively similar horizontal distributions in Fig. 2, in contrast to the distributions in GO and GW. These four simulations are relatively reasonable at reproducing the PRISM rainfall field over western CONUS. However, they commonly overestimated rainfall compared to the PRISM rainfall over Florida, the Appalachian Mountains and their eastern side, the southern parts of the Rocky Mountains, and the northern parts of the Great Plains around the 105°W meridian. They all underestimated over a central region of CONUS, which was located in the middle of both regions exhibiting an overestimation. The simulation biases in Fig. S1 are coincident with those of the MERRA2 rainfall against the PRISM data in Fig. 1c. This result demonstrates that common biases in the NU-WRF simulations partially originated from those inherent in MERRA2 through the lateral boundary conditions and the strong spectral nudging.
The use of the UWSC scheme in GW uniformly reduced the simulated rainfall, as compared to the result of G. As a result, it alleviated the positive biases over the east Atlantic regions in Fig. S1c. However, it simultaneously caused negative biases over the central CONUS, as compared with the plots of GO and G. This effect of the UWSC scheme on simulating rainfall is larger than those due to the selection of the deep cumulus parameterizations, that is, among BO, KO, and SO as well as G. In contrast, applying the native shallow convection component of G3D in GO increased the simulated rainfall. The GO figures (Fig. 2a and Fig. S1a) exhibit the greatest overestimation among the simulations throughout CONUS except for the Pacific coast and limited areas of the northern Rocky Mountains.
The regional variability of the daily rainfall averaged for JJA2000 is further discussed on the basis of the regionalization suggested by Bukovsky (2011). The Bukovsky regionalization divides the CONUS domain into 17 subregions. Figure 3 illustrates JJA2000-average daily rainfall rates of PRISM, MERRA2, and NU-WRF simulations and their overall biases averaged for each Bukovsky subregion over CONUS. PacificNW, PacificSW, and GreatBasin subregions (these notations follow the Bukovsky regionalization) are characterized by a dry summer, and all simulations exhibit very similar rainfall rates. In contrast, the simulated rainfall rates are much more sensitive to the cumulus parameterization outside of these three subregions. The relative magnitudes of the simulated rainfalls follow roughly those averaged throughout CONUS, that is, GO > KO > G > SO > BO > GW in Table 1. The overall biases compared to PRISM and MERRA2 in Figs. 3b and 3c differ widely according to the subregions. GO exhibits the largest rainfall rates among the simulations in Fig. 3, and its values are substantially larger than those of the other simulations and reference datasets (except as previously noted for PacificNW, PacificSW, and GreatBasin). KO yields generally larger rainfall rates than other simulations, except over some western subregions, that is, GreatBasin, Southwest, SRockies, and Mezquital. It generally overestimates the PRISM rainfall rates over most of the subregions, whereas it is in good agreement with PRISM over SPlains, Prairie, and DeepSouth. The results of G are similar to those of KO, especially for most of the central subregions. BO and SO provide higher rainfall rates than PRISM on the western side from SPlains; lower rates over SPlains, Prairie, and DeepSouth; and higher rates over MidAtlantic and NorthAtlantic.
Figure 4 shows the frequency distribution (divided by the total number of samples) of all pointwise rainfall rates in the 92 days of JJA2000 on the upscaled PRISM grid points at horizontal intervals of 0.25° over the entire CONUS. The NU-WRF simulation and MERRA2 rainfall fields have been regridded to the upscaled PRISM grid points using bilinear interpolation. The NU-WRF simulations and MERRA2 overestimate the frequency of rainfall less than 2 mm day−1, compared to PRISM. The probable reason is that rainfall simulated with relatively coarse resolution tends to be smoothed out and too widespread. The use of bilinear interpolation possibly enhances the increase in the frequency of the light rainfall. The MERRA2 frequency spectrum is in better agreement with that of PRISM than those of most NU-WRF simulations, at most rainfall ranges except for those less than 5 mm day−1. NU-WRF simulation spectra show different rainfall ranges where simulated frequency agrees well with that of PRISM: G simulates the frequency of the PRISM rainfall intensity reasonably at 25–50 mm day−1, while the rainfall range of KO matching PRISM is limited to 35–50 mm day−1. GW’s spectral frequency is the most similar to that of PRISM at a range less than 10 mm day−1. The spectrums of BO and SO intersect the PRISM spectrum around 20 mm day−1. In general, rainfall rates where the NU-WRF-simulated frequency agrees with that of PRISM are in proportion to the overall biases of the simulations compared to PRISM (Table 1).
The simulated rainfall rate is the total derived from the cloud microphysics scheme and from the cumulus parameterization scheme. The horizontal fields of JJA2000-average daily rainfall of the two components are shown in Figs. 5 and 6; Fig. 7 also illustrates the fraction of cumulus parameterization contribution to the total rainfall rates. Higher microphysics-based rainfall is simulated over some of the Northeast and a limited area of the Northwest coast, whereas the rainfall from cumulus parameterization schemes dominates over most of the central and southeastern parts of CONUS. The rainfall rate from the cloud microphysics scheme is also sensitive to the cumulus parameterization scheme. For example, a considerably wide distribution of the rainfall is simulated in GO only in Fig. 6a. The rainfall rate produced by UWSC in GW is exceedingly small (not shown) compared to those from the cloud microphysics and main cumulus parameterization schemes. The use of UWSC strongly suppresses rainfall from the main cumulus parameterization scheme by interfering with its convection trigger through stabilization, but has little influence on rainfall from the cloud microphysics. The contribution of cumulus parameterization to the total rainfall rate in Fig. 7 is dependent on the geographical location, type of parameterization scheme, and use of UWSC. For example, high fractions close to 1.0 are simulated in southern regions of CONUS commonly in all simulations. In contrast, the northern regions exhibit relatively strong dependency on cumulus parameterizations. KO and SO show larger fractions close to 1.0, whereas G and BO exhibit smaller fractions.
b. Correlation and RMSE analysis
Figure 8 contains Taylor diagrams (Taylor 2001) summarizing the correlation coefficients, standard deviations, and centered root-mean-square error (RMSE) of daily rainfall rates over CONUS against PRISM, MERRA2, and NU-WRF simulation G as the references. Simulation G was selected as a reference for the diagrams to show how statistics change due to different shallow convection components as well as different deep cumulus parameterizations. To calculate the statistics, the NU-WRF simulation and MERRA2 rainfall fields were regridded to the upscaled PRISM grid points at horizontal intervals of 0.25° using bilinear interpolation. The statistics in Figs. 8a–c were calculated from the JJA2000-average daily rainfall rates, which were shown in Figs. 1 and 2. Thus, Figs. 8a–c exhibit spatial pattern statistics of the seasonally averaged daily rainfall. In contrast, the statistics in Figs. 8d–f were calculated from nonaveraged pointwise rainfall rates in the 92 days of JJA2000. Thus, Figs. 8d–f show statistics containing both spatial and temporal (daily) variation. Here, the two types of correlation coefficients for the spatial variation in the seasonally averaged daily rainfall fields and for the spatial–temporal variation of rainfall fields in all days of JJA2000 are referred to as CORS and CORST, respectively.
The comparison of Figs. 8a–c shows that the spatial patterns of the seasonally averaged rainfall in the simulations are more correlated with that of the simulation G than that of PRISM. The values of correlation coefficients (CORS) against MERRA2 are between those of the simulation G and PRISM, except for the case of GO. The higher CORS against MERRA2 than against PRISM indicate that the rainfall spatial patterns of the simulations are constrained by the atmospheric fields of MERRA2 through the influence of lateral boundary conditions and strong spectral nudging. In these three diagrams, GO is commonly displaced from the group of the other members because of rainfall significantly increased by the native shallow cumulus components in G3D. The use of UWSC in GW decreases the standard deviation and centered RMSE of the simulated rainfall by suppressing rainfall. GO and GW are not considerably more correlated with G than BO, KO, and SO in Fig. 8c. The shallow cumulus components in GO and GW thus have some influence on changing not only the magnitude of rainfall but also its spatial pattern, similar to the choice of deep cumulus parameterizations.
From Figs. 8a–c to corresponding Figs. 8d–f, the correlations have been similarly degraded, even from Figs. 8c to 8f. This result suggests that the timing of rainfall, even for the same location, is sensitive to the selection of deep and/or shallow cumulus parameterizations, as reported by previous studies. The diurnal variation of simulated rainfall over CONUS is significantly dependent on cumulus parameterization employed (e.g., Liang et al. 2004; Qiao and Liang 2015). Furthermore, some parts of CONUS show diurnal rainfall peaks around 0000 UTC (e.g., Matsui et al. 2010, their Fig. 2), so that small changes in rainfall timing could have large impacts on daily (0000–2400 UTC) rainfall rates. Among Figs. 8d–f, the simulations are more correlated to G (and to each other) than to PRISM or MERRA2. The simulations have a somewhat common rainfall pattern in both spatial and temporal variation, regardless of the selection of cumulus parameterization.
Figure 9 illustrates the regional variability of the two types of correlation coefficients (CORS and CORST) for the Bukovsky subregions against PRISM and MERRA2. Figure S2 shows the same statistics but against the simulation G as the reference dataset. The correlation for the spatial pattern only in Figs. 9a and 9b exhibits significant regional variation. There is no individual simulation that yields the best or worst correlation over most subregions. The series of CORS against both PRISM and MERRA2 are commonly high over some western subregions, that is, PacificNW, Southwest, SRockies, and NPlains. The CORS of the simulations with each other (an example against G is shown in Fig. S2) are very high over most western subregions. However, the values of CORS are very different according to the selection of cumulus parameterization over regions between the Rocky Mountains and the East Coast, including the Appalachian Mountains. Some simulations yield negative correlation coefficients there. These results demonstrate that the spatial pattern of seasonally averaged rainfall is strongly constrained by the topographical pattern over mountainous regions. In contrast, over the plains regions, the spatial pattern is more dependent on the native thermodynamical processes in the model and hence more sensitive to the cumulus parameterization than topographic effects.
The functions of CORST in Figs. 9c and 9d show quite different structure from those in Figs. 9a and 9b for CORS. The values of CORST against PRISM in Fig. 9c are roughly about 0.3 over all subregions, except for PacificNW and PacificSW, and less sensitive to the selection of cumulus parameterization. The regional variability and sensitivity to cumulus parameterization are reduced by including temporal variation. High positive correlations in Fig. 9a have been uniformly degraded in Fig. 9c, whereas negative correlations over some plains subregions have been eliminated. The CORST against MERRA2 in Fig. 9d is generally higher than that against PRISM in Fig. 9c, except for Southwest, Mezquital, and Southeast. Over these southern-edge subregions, almost all rainfall is produced by cumulus parameterization schemes as shown in Fig. 7, so that the resolution difference between the NU-WRF simulations and MERRA2 may have a stronger impact on CORST than CORS.
Figure 10 illustrates the regional variability of the full RMSE against PRISM and MERRA2. The difference between the definitions of centered and full RMSE can be found in Taylor (2001). The full RMSEs are generally much higher in the case of spatial–temporal variation in Figs. 10c and 10d than those of spatial variation only in Figs. 10a and 10b. The increases in the full RMSE from Figs. 10a and 10b to Figs. 10c and 10d are caused by the increases in the centered RMSE, because the contributions of the overall bias to the full RMSE are the same. In Figs. 10a and 10b, the full RMSEs are mostly attributed to the overall bias error because centered RMSEs are not very large. For example, GW yields the smallest full RMSE from NRockies to NPlains. These smallest full RMSEs correspond to the least absolute values of the overall bias of GW over these subregions in Fig. 3b. The link between the full RMSE and the overall bias is stronger over mountainous subregions because of high CORS among the simulations there (e.g., Fig. S2), whereas the link is slightly weaker over the plains subregions. In Figs. 10c and 10d for spatial–temporal variation, the relative magnitudes of the full RMSEs are somewhat less sensitive to the overall bias, although weak relationships between the bias and the full RMSE can be found over mountainous subregions.
c. Estimation of the best mixture of individual simulation results on the basis of least-RMSE optimal ensemble calculation
The sensitivity to the cumulus parameterization scheme is further explored by the calculation of optimal weighted ensembles of multiple individual simulations. The relative magnitudes of weighting coefficients in optimal weighted ensembles represent relative suitability of each cumulus parameterization scheme. Four types of least-RMSE weighted ensemble listed in Table 2 were calculated using the approaches described by Liang et al. (2007). The daily rainfall rates of the weighted ensembles are given by the following equations:
where R is the daily rainfall rate, x is the position vector, and α, β, and γ represent weighting coefficients. Subscripts G, BO, KO, and SO indicate the values of the corresponding individual simulations in Table 1; here G is selected from the three simulations using G3D, because its rainfall distribution is more similar to BO, KO, and SO. Variable R with an overbar represents daily rainfall rates averaged for JJA2000. Equation (1) optimizes the spatial pattern of simulated horizontal fields averaged for JJA2000. In contrast, Eq. (2) is suited to the prediction of pointwise rainfall in all days of JJA2000 with both spatial and temporal variation; the prime symbols in Eq. (2) represent the weighting coefficients and the ensemble obtained for the latter. The weighting coefficients (α, β, γ, α′, β′, and γ′) are calculated for the sets by minimizing the full RMSE of the ensemble against PRISM or MERRA2 through a simple linear regression.
Tables S1–S4 show a series of the weighting coefficients obtained for the four types of ensembles (α, β, γ, and 1 − α − β − γ for WEP_SP and WEM2_SP; α′, β′, γ′, and 1 − α′ − β′ − γ′ for WEP and WEM2 in Table 2) calculated for the entire CONUS or for each Bukovsky subregion. Figures 11 and 12 illustrate the values in Tables S1–S4 on the CONUS map partitioned by the Bukovsky subregions. The relative magnitudes of the weighting coefficients are roughly inversely correlated with the magnitude of full RMSE in Fig. 10.
Figures 11a–d (spatial variation only against PRISM) show that G performs better in western regions around the Rocky Mountains and the northeastern region. BO and KO outperform the others over NPlain and Prairie Bukovsky subregions, respectively. SO generally shows higher weighting coefficients over the eastern half of CONUS. As compared to Figs. 11e–h (spatial–temporal variation against PRISM), Figs. 11a–d suggest that the optimal ensemble over a subregion tends to coincide with a single suitable simulation rather than a mixture of multiple simulations. In the discussion about the spatial pattern statistics in section 3b, the spatial pattern of the JJA2000-average rainfall is similar among the simulations, especially over mountainous regions (e.g., Fig. S2). A mixture of multiple simulations thus does not yield an improvement of the spatial pattern over these regions. As a result, the overall bias shown in Fig. 3b is the key factor to determine the full RMSE. As long as the biases of the four simulations have the same sign (positive or negative), the smallest-bias (absolute value) simulation is identical to the smallest-RMSE simulation. The smallest-bias simulation outperforms the other single simulations and a mixture of multiple simulations over the region. An example is G over NRockies and SRockies. The situation is different over PacificNW, where BO shows a negative bias but G, KO, and SO have positive biases in Fig. 3b. A mixture of multiple simulations is effective to reduce the full RMSE of the optimal ensemble by offsetting the biases in each case.
Over the plains regions in Figs. 11a–d, the overall bias has less impact on determining the optimal ensemble, because the spatial patterns vary among the simulations. For example, KO exhibits the smallest absolute value of the bias in Fig. 3b and the least full RMSE over SPlain in Fig. 10a, but a mixture of all four individual simulations corresponds to the optimal ensemble in Figs. 10a–d. The mixture of multiple simulations has a positive effect on reducing the full RMSE over the subregion by coupling different spatial patterns in the different simulations.
In Figs. 11e–h (spatial–temporal variation against PRISM), the optimal ensembles are composed of mixtures of multiple simulations over almost all subregions. The optimal ensembles are roughly identical to arithmetic mean ensembles, that is, weights of 0.25 for all four simulations. Since the correlation for spatial–temporal variation is lower than that for the spatial variation only (e.g., Fig. S2), a mixture of multiple simulation results is effective in reducing the full RMSE. Although the overall biases have less impact on the full RMSE in this case, the trend in relative magnitudes of the weighting coefficients in Figs. 11a–d for spatial variation only remains in Figs. 11e–h for spatial–temporal variation, especially over mountainous regions.
There are common characteristics between the weighting coefficients of the optimal ensembles against PRISM in Fig. 11 and those against MERRA2 in Fig. 12, in terms of overall distribution of the weighting coefficients and the difference between the cases of the spatial variation only and spatial–temporal variation. BO generally shows slightly higher coefficients against MERRA2 than those against PRISM. KO also exhibits overall better performance, except for Prairie in the spatial variation only case. These results again indicate that the simulations using BO and KO are more constrained by MERRA2.
Figure 13 illustrates scatterplots for a set of two variables sampled over each Bukovsky subregion: the domain-mean fraction of cumulus parameterization contribution to the total rainfall rates and the weighting coefficients of the smallest-RMSE weighted ensembles for spatial–temporal variation against PRISM and MERRA2 (Figs. 11e–h and 12e–h). Positive correlation in the scatterplots indicates that the cumulus parameterization performs better over subregions where subgrid convection processes greatly contribute to the precipitation formation, and vice versa.
Figure 13 shows that KO exhibits negative correlations in both panels. In contrast, SO shows positive correlations, and G and BO have almost no correlation against MERRA2 in Fig. 13b, whereas they yield negative and positive correlations against PRISM in Fig. 13a, respectively. In general, NKF and G3D perform better over regions where the role of subgrid precipitation processes is limited, and the reverse is true for BMJ and NSAS. These characteristics are consistent with the fact that NKF and G3D have been employed and tested more frequently in midlatitude regional climate simulations (e.g., Liang et al. 2007; Mearns et al. 2009), in which the contribution of subgrid parameterization seems to be small compared with the tropics in general circulation model (GCM) simulations. However, BMJ and NSAS are more suitable for the present simulations with a 24-km horizontal resolution for daily rainfall over CONUS for JJA2000, because the subgrid convection processes dominate the rainfall over most regions as shown in Fig. 7.
The large scatter of the sample points in Fig. 13 suggests large variability in the relationship between the two variables. Many other factors, such as the suitability of the cumulus parameterization with the orographic effects over the subregions, likely affect the determination of the weighting coefficients for the optimal ensembles.
4. Summary and conclusions
A series of regional seasonal simulations with a 24-km horizontal resolution over CONUS were conducted using NU-WRF with six different cumulus parameterization schemes. This study explored the sensitivity of simulated daily rainfall rates in JJA2000 to the cumulus parameterization schemes selected. The simulations were evaluated by comparing them with the daily rainfall products in PRISM and MERRA2. The analysis results are summarized as follows:
The relative magnitudes of the JJA2000-average rainfall rates among the simulations are roughly common to the average across CONUS and to the averages over each Bukovsky subregion. The simulated rainfall distribution is generally similar to the rainfall distribution in MERRA2, because the simulated atmospheric fields are constrained by the variable fields in MERRA2 through lateral boundary conditions and strong spectral nudging. The overall bias over a subregion is thus influenced by the bias in boundary conditions regardless of cumulus parameterization. G3D and NKF produce more rainfall than BMJ and NSAS, except for regions around the Rocky Mountains. The different rainfall rates and responses according to the cumulus parameterization are attributed to the difference in the magnitudes of the cloud-base mass flux assumed in the parameterizations. The relatively strong rainfall of NKF is mostly attributable to the design of its closure algorithm parameterizing cloud-base mass fluxes (Qiao and Liang 2016).
The effects of the selectable G3D shallow components on the rainfall are more significant than the effects of switching deep cumulus parameterization schemes. The use of UWSC for G3D uniformly reduces the simulated rainfall, whereas applying the native shallow convection component in G3D increases it. The results of G3D simulation without either UWSC or the native shallow component are the most similar to those of the others, that is, BMJ, NKF, and NSAS with their native shallow convection components.
Because the spatial pattern of rainfall is largely constrained by the topographical pattern over mountainous regions, the spatial pattern of the seasonally averaged rainfall is similar across all simulations and with PRISM and MERRA2. The overall bias against the reference datasets is the key factor in determining the full RMSE over mountainous regions. In contrast, over the Great Plains, the spatial pattern is more sensitive to the selection of cumulus parameterizations. However, even over the mountainous regions, the correlation for spatial–temporal variation in the case of sampling the nonaveraged pointwise rainfall in all days is relatively very low. The lower correlation is likely due to the poor simulation of temporal variation of the daily rainfall over CONUS regardless of the choice of the parameterizations.
The analysis of optimally weighted ensembles highlights regional variability in the performances of the cumulus parameterizations. In constructing optimal ensembles for the spatial pattern of the seasonally averaged rainfall, the adoption of a single simulation with the least bias is generally preferable to a mixture of multiple simulations, as long as all simulations have the same sign (positive or negative) of biases. A mixture of multiple simulation results, even an arithmetic mean ensemble, is more effective in constructing an optimal simulation result when the focus is on nonaveraged rainfall in all days throughout the season.
Overall, G3D and NKF performed better over regions where the roles of subgrid precipitation processes are minor, and the reverse is true for BMJ and NSAS. These differences are likely attributable to the different parameterization designs optimized for certain spatial scales and climate zones.
The present study has investigated the sensitivity to different cumulus parameterization in regional simulations with a 24-km resolution for the summertime of year 2000. The scale awareness of the sensitivity should be further examined in a future study, because the regional variability of the schemes’ performance could be different if it is dependent on their intended model scales and roles of subgrid parameterization. In addition, assessment of the interannual variation and long-term climatology of the sensitivity are relevant to improved assessment of climate projections.
This study was supported by the NASA Downscaling Project led by Tsengdar Lee at NASA headquarters. Computational resources for this project were provided by two NASA High Performance Computing (HPC) facilities: NASA Center for Climate Simulation (NCCS) at Goddard Space Flight Center (GSFC) and NASA Advanced Supercomputing (NAS) Division at Ames Research Center. The MERRA2 data were provided by the NASA Global Modeling and Assimilation Office. The script for plotting Taylor diagrams was provided by Yannick Copin of the University of Lyon. The authors thank the two anonymous reviewers and the journal editor for their helpful comments in improving this paper.
Supplemental information related to this paper is available at the Journals Online website: http://dx.doi.org/10.1175/JHM-D-16-0120.s1.