Evaluation of downscaled meteorological information is crucial to identifying model behaviors that may propagate to end applications such as the simulation of local air quality. This study conducted and assessed yearlong simulations of hourly meteorological conditions over the Terrace–Kitimat Valley of northwestern British Columbia, Canada, at 1-km horizontal gridding for six PBL schemes in the Weather and Forecasting (WRF) Model, version 4.0. In terms of key surface meteorological variables that affect air quality, simulations over land demonstrated better skill for specific humidity and wind direction than for air temperature and wind speed. Spatial differences in modeled atmospheric properties and vertical profiles, especially for moisture content, were used to diagnose the relative capacity of each PBL scheme to represent pollutant dispersion and dilution. Stable conditions at night increased suppression of boundary layer mixing by the nonlocal Yonsei University (YSU) scheme when compared with suppression by the local eddy-diffusion component of the Asymmetric Convective Model, version 2 (ACM2), scheme, resulting in decreased wind speed and ambient temperature but moister air with the YSU scheme. The weakening of mixing by the Mellor–Yamada–Nakanishi–Niino (MYNN3) scheme with inland distance suggested that higher-order, nonlocal transport is sensitive to increasing topographic steepness toward the northern part of the valley. Disparities in mixing strengths among PBL schemes were greater in the summer when conditions were generally less stable with moist, warm air blowing inland than in winter when the valley channels cold, stable air from the interior. Increased convection in daytime led to greater entrainment of air from aloft and a thicker PBL with the YSU scheme than with the ACM2 scheme in summer while increasing countergradient transport in the MYNN3 scheme that reduces dilution.
In air quality modeling, correct representation of the wind, temperature, moisture, and mixing heights is needed to predict the fate of pollutants. It is generally recognized that the treatment of planetary boundary layer (PBL) and surface-layer processes in meteorological models has a direct impact on predicting these variables. Mixing by turbulent eddies within the PBL drives the distribution of atmospheric quantities (Cohen et al. 2015; Nielsen-Gammon et al. 2010). Although several formulations describing atmospheric turbulence and interactions with surface fluxes exist, some are not widely tested. Apart from varying spatial and temporal resolutions of atmospheric models, choice of dynamics and physics parameterizations can impact meteorological fields (Mao et al. 2006) that are used in air quality models.
Many investigations of the sensitivity of coupled numerical models to PBL schemes have been conducted in the last decade. Cheng et al. (2012) simulated ozone episodes over Taiwan and found that discrepancies between results were due to vertical mixing strength and entrainment properties of the tested PBL schemes. Banks and Baldasano (2016) evaluated how outputs from some PBL schemes in the Weather Research and Forecasting (WRF) Model differ from one deployed for operational air quality forecasts over Catalonia, Spain. In India, Gunwani and Mohan (2017) researched the sensitivity of PBL schemes to different climatic zones based on atmospheric model simulation of some meteorological parameters. Others (e.g., Borge et al. 2008; Cuchiara et al. 2014; Reboredo et al. 2015; Xie et al. 2012) have quantified meteorological prediction differences to ascertain which PBL scheme is most appropriate to specific geographical domains. However, much of the related literature (e.g., Cuchiara et al. 2014; Mohan and Gupta 2018; Misenis and Zhang 2010) are for short episodes and do not portray cross-temporal facets of statistical measures. In some regions, where pollution episodes may be uncommon and routine contaminant exposures from all other periods in a year are more of a concern, airshed managers could be more interested in simulation of air quality over annual, seasonal, or even diurnal periods. There is a general lack of studies on the effect of different PBL parameterization schemes in the downscaling of meteorological variables for modeling the air quality over the northwest coast of Canada. Since outputs from meteorological models are commonly fed into air quality models, ultimately weighing on airshed planning and policies, it is imperative to assess the reasonableness of downscaled meteorology prior to their use for current assessments and future projections.
The geographically complex Terrace–Kitimat Valley (TKV; Fig. 1) in northwestern British Columbia, Canada, presents an ideal area to test the ability of a mesoscale model to simulate meteorological variables that are of importance from an air quality perspective. Aligned north–south within the Pacific Coast Mountains and flanked by ridges 1500–2000 m high, the TKV lies at the head of a 90-km fjord. Roughly 70 km in length, the valley is broadest around the small city of Terrace but is as narrow as 8 km at some sections farther south and near the port town of Kitimat. The valley’s weather is dictated by air masses originating over both the Pacific Ocean and continental North America. In summer (June–August), when the Pacific high lingers near its northernmost position (Klock and Mullock 2000), clear, fine weather prevails. Surface wind is mostly southerly from the Douglas Channel, and land–sea breezes and slope–valley winds are experienced. In winter (December–February), landfalling midlatitude cyclones bring heavy snowfall, and the intermittent clash of onshore winds and outflowing continental air results in stormy weather (Lange 2003). The mean annual temperature in Kitimat and Terrace is 6.9° and 7.4°C, respectively, varying from a minimum of −4°C mean in January to a maximum of +22°C mean in July; annual precipitation amount is 2210 mm in Kitimat but 1170 mm in Terrace (Environment and Climate Change Canada 2019), in the partial rain shadow of the Coast Mountains. The TKV is sparsely populated; however, considerable amounts of sulfur dioxide emitted from an aluminum smelter in Kitimat are advected inland by onshore winds. Recent proposals for the construction of more industries, such as liquefied natural gas plants (District of Kitimat 2020), in the area has raised concerns about air pollution and deterioration of the natural environment.
This study conducts high-resolution simulations of meteorological variables with the WRF Model, version 4.0, and the fifth major global reanalysis product of the European Centre for Medium-Range Weather Forecasts (ERA5), with a focus on evaluating performance of different PBL formulations. Its objectives are to quantify modeling uncertainties and to qualify mixing and dispersion capabilities of alternative physical schemes. The remainder of this paper is organized as follows. In section 2, the experimental setup, including the tested PBL schemes being tested, nested domains, and WRF Model configuration, is described. Section 3 evaluates the uncertainty in modeled surface meteorological fields through comparisons with observations. Also, differences among simulations resulting from PBL scheme options are analyzed. Section 4 discusses the implications of simulation similarities and differences for air quality modeling and concludes the paper.
a. Description of the WRF Model and PBL schemes
WRF is a state-of-the-art atmospheric modeling system designed for both meteorological research and numerical weather prediction. The dynamical solver, known as the Advanced Research WRF (ARW) core, is based on the fully compressible, nonhydrostatic Euler equations with terrain-following eta coordinates. WRF can run on a variety of computing platforms and suits a broad range of applications across scales ranging from tens of meters to thousands of kilometers (Skamarock et al. 2019).
Turbulence closure schemes, also referred to as PBL schemes are built in to the meteorological models to represent turbulent fluxes that are not explicitly resolved on the model grid. PBL processes are parameterized either through local closure or nonlocal closure models. In the case of local closure, only those vertical levels that are adjacent to a given point directly affect variables at that point. Nonlocal closure schemes use multiple vertical levels to determine variables at a given point (Cohen et al. 2015). For our research with the WRF Model, we tested six such schemes. These schemes are the Mellor–Yamada–Janjić (MYJ); Asymmetric Convective Model, version 2 (ACM2); Mellor–Yamada–Nakanishi–Niino, level 3 (MYNN3); Shin–Hong (SH); University of Washington (UW); and Yonsei University (YSU) schemes.
The YSU scheme is a first-order nonlocal scheme with explicit entrainment at the PBL top and a parabolic K profile in an unstable mixed layer. It calculates PBL height from the surface and uses a threshold Richardson number of 0.25 for stable cases and 0 for unstable flow (Hong et al. 2006; Hong 2010). The SH scheme is a nonlocal “scale aware” formulation for vertical transport in a convective PBL. Vertical mixing in the stable PBL and free atmosphere follows YSU, in addition to diagnosed turbulent kinetic energy (TKE) and mixing length, but the explicit treatment of heat flux entrainment is replaced by grid size-dependent terms for nonlocal, and local transport components (Shin and Hong 2015). The ACM2 scheme is a hybrid scheme that features nonlocal transport from the lowest level to all other model layers alongside a local eddy-diffusion component, and an asymmetrical layer-by-layer downward transport. It calculates PBL height when the critical bulk Richardson number above the level of neutral buoyancy exceeds a value of 0.25. For stable or neutral flows, the scheme shuts off nonlocal transport and uses local closure (Pleim 2007). The MYJ scheme is a one-and-a-half-order local closure that includes a prognostic equation for TKE. The PBL height is the height at which the TKE reaches 0.2 m2 s−2 (Janjić 1994). The MYNN3 scheme is a second-order closure that expresses stability and mixing length based on the results of large-eddy simulations rather than observations. It prescribes PBL height when TKE reaches 1.0 × 10−6 m2 s−2 (Nakanishi and Niino 2009). The UW scheme is a one-and-a-half-order local closure in which TKE is a diagnosed, rather than a prognostic variable. It uses moist-conserved variables with an explicit entrainment closure for convective layers. A threshold of 0.25 for the critical bulk Richardson number is used to determine PBL height in all stability cases (Bretherton and Park 2009).
The above schemes are a mix of legacy and newer physics modules that have accompanied major updates to the WRF code. For instance, whereas the YSU, MYJ, and ACM2 schemes are common in the literature (e.g., Nielsen-Gammon et al. 2010; Zhang et al. 2013), the UW and SH schemes are fairly recent additions. Consequently, selected PBL schemes represent the state-of-the-art in options for parameterizing atmospheric processes over predominantly pristine (nonurbanized) areas.
Each PBL formulation can be mapped to just one or more surface-layer schemes. The surface-layer schemes compute friction velocity and other exchange coefficients for estimating sensible and latent heat fluxes, and momentum flux from land surface models, and surface stress in the PBL schemes. In this study, the MYJ, MYNN3, and UW PBL schemes are coupled to the Eta Model surface-layer scheme whereas the ACM, YSU, and SH PBL schemes are mapped to a modified MM5 surface-layer scheme (Jiménez et al. 2012). Both surface-layer schemes are based on similarity theory (Monin and Obukhov 1954), but the Eta Model similarity includes a parameterization of viscous sublayer according to Janjić (2002).
b. Domain configuration, model initialization, and settings
Computational domains assumed a telescopic nesting arrangement (Fig. 1). A grid ratio of 1:5 applied to successive parent domains on a Lambert conformal projection was used to attain a horizontal resolution of 1 km for the TKV and surrounding areas (Table 1). Because of numerous lakes in the valley, terrain data were interpolated from the Moderate Resolution Imaging Spectroradiometer (MODIS) International Geosphere–Biosphere Programme (IGBP) 21-category land-use and land-cover fields with representation for lake bodies. This was to effectively distinguish between inland water bodies and seas, for appropriate interpolation of time-varying SST fields. Atmospheric forcing data (for initial and boundary conditions) for the entire setup was ERA5 (ECMWF 2019). Simulations up to a model top of 50 hPa were then performed for the year 2017. The 2017 meteorological year was selected because it best represents the climatological average for the valley from 2006 to the present. For all simulations, physics settings were the Noah land surface model (Tewari et al. 2004), the revised Thompson scheme for microphysics (Thompson et al. 2008), the Rapid Radiative Transfer Model (Iacono et al. 2008) for longwave and shortwave radiation, and the horizontal Smagorinsky first-order closure (Talbot et al. 2012) for mixing terms. The Grell–Freitas ensemble cumulus physics scheme (Grell and Freitas 2014) used for the largest and intermediate domains was not required for the innermost (study) domain (Arakawa et al. 2011) because of its sufficiently high resolution. A one-way nesting (http://www2.mmm.ucar.edu/wrf/OnLineTutorial/CASES/NestRuns/ndown4.php) was implemented. For each PBL scheme, the WRF Model was run between 0000 UTC 20 December 2016 and 0000 UTC 2 January 2018. Model runs were initialized monthly, with model spinup of 1 day, and a further 1 day at the end of each month. The monthly overlap days, including for the period 20–31 December 2016, were discarded when merging hourly output files.
c. Observational data
Three surface meteorological variables were selected for the evaluations. These are 2-m temperature (T2), specific humidity at 2 m (Q2), and 10-m wind (speed and direction). These variables were selected because they influence ambient air quality (Banks and Baldasano 2016; Dennis et al. 2010; Miao et al. 2008) and are monitored at various locations in the valley (Table 2). Observation stations are operated and maintained by several organizations, including Environment and Climate Change Canada, British Columbia Ministry of Transportation and Infrastructure, and British Columbia Ministry of Environment and Climate Change Strategy’s Air Quality Network. (Station data are considered independent of ERA5 and can be downloaded from https://data.pacificclimate.org/portal/pcds/map/.) Because each network is designed to meet the specific need of its operating agency, there is variety in the site characteristics and sensor equipment. For each parameter at any station, only records with at least 80% yearlong completeness in hourly values are used for evaluation.
d. Statistical evaluation indices
For evaluations, a cell–point method is used to compare the simulation results from the WRF Model with observations at different locations. In this regard, values for model grid cells corresponding to ground locations of observations were retrieved. The fitness measures used in this work are the mean bias (MB), the root-mean-square error (RMSE), Pearson’s correlation coefficient r, and index of agreement (IOA) according to Willmott et al. (2012). These indices are defined in appendix A and were selected to reflect a mix of accuracy measurements for paired data. They also enable comparison to performance benchmarks (Table 3) for variables within the boundary layer suggested by Emery et al. (2001). Whereas MB would compensate between equally weighted positive error (overprediction) and negative error (underprediction), the RMSE would penalize large errors since pairwise differences are squared before averaging (Draxler and Chai 2014). MB as well as RMSE are in units of the parameter of interest and are negatively oriented indices, which means values closer to zero are better. For r and IOA, however, scores range from −1 to 1 and values approaching 1 indicate good performance. The identification of the most suitable PBL scheme is from a count of instances in which it is top ranked and also meets the performance benchmark. This is checked per statistical measure for all stations.
3. Results and analyses
a. Synoptic and mesoscale representations
PBL schemes drive the vertical distribution of fluxes in the atmosphere and soundings of temperature, humidity, and wind alongside modeled profiles from each PBL scheme are examined qualitatively. All simulations reasonably mimic the above aerological variables for different times in winter and summer seasons at the Annette Island radiosonde station, Alaska, which is the upper-air sounding location nearest the TKV (Fig. B1 in appendix B). This station is outside the 1-km grid but within the 5-km grid. Modeled temperature and dewpoint profiles by the various schemes virtually match observations. Simulated dewpoint profiles for 0400 Pacific standard time (PST) on a January day replicated a drying lower troposphere, with coincident profiles of the YSU and SH schemes being drier than observations at 700-hPa height above the surface. Representation of the veering wind is also reasonable, given that specific days and times could have unusual atmospheric conditions. It needs emphasizing, however, that the atmospheric environment around the Annette Island station, at least for the lowest 2000 m, could differ markedly from that of the TKV, since the station is west of the ranges of the Coast Mountains, completely exposed to Pacific Ocean air masses. Despite physiographical differences between the radiosonde location and stations in the TKV, the fact of general correspondence for upper-air data attests to a reasonable representation of synoptic information.
Spatial plots of near-surface meteorological conditions for typical summer (22 July) and winter (21 January) days over the study (TKV) domain are presented in Figs. 2 and 3, respectively. Specific times are hours of diurnal temperature maximum and minimum over land, which are relevant for assessing mesoscale circulations such as land–sea breezes and mountain-slope winds. One verification of WRF simulation realism is the altitudinal decrease of daytime air temperature in summer, from around 18°C on the valley floor to about 4°C over the mountains, which is consistent with expected summertime conditions (Fig. 2). Wind vectors are also more organized and aligned with the main and tributary valleys than outside them, indicating substantial wind adjustment to the underlying terrain and the efficacy of simulation at high horizontal resolution. Over the main valley itself, modeled winds are seasonally nuanced. In summer (Fig. 2), the wind is southerly, which is representative of the up-valley flow that occurs for about 70% of noncalms (wind speed > 0.5 m s−1) observed at the Terrace station. This pattern is also dominant overnight (0200 PST) in all the PBL schemes, perhaps due to fewer nighttime hours in summer for reversal of wind direction, or from stronger synoptic forcing on this day, but terrain recognition by the WRF Model appears satisfactory. Southerly wind vectors turn anticlockwise on approaching the Skeena River, resulting in a flow down the river channel, toward the coast. Some convergence of vectors toward watercourses that drain the valley, suggesting downslope winds in tributary valleys is also present in the simulations at 0200 PST. The afternoon hour (1400 PST) mainly indicates wind directions out of stream channels, or upslope. Compared to the nighttime period, the southerly flow in the valley is as would be expected, more intense due to enhancement by sea breezes. Wind vectors within the valley are straighter and better organized, with slightly more buildup on valley sidewalls for the MYNN3 and MYJ schemes.
The winter season (Fig. 3) when the valley is snow covered, differs markedly from summer. At night (0200 PST), the warmest parts are around the shoreline, and over major inland water surfaces. This is as expected since water bodies cool more slowly than land. Over the ridges flanking the valley, the MYNN3, MYJ, and UW PBL schemes depict warmer temperatures than the ACM2 and YSU PBL schemes, probably due to differences in surface-layer schemes. By afternoon, the land areas have warmed such that both the valley and shore areas are mostly above 0°C. Indeed, at 1400 PST, the northern half of the valley is warmer than the southern portion, with nearly similar temperatures at mountain peaks (from −9° to −6°C) among PBL schemes. Wind vectors for the winter season are consistent with the dominance of northerly outflow at this time of the year; however, these are less organized than during summer. While the flow within the valley remains fairly distinct from winds over surrounding ridges, wind vectors in the valley are less organized than in summer suggesting, stronger synoptic forcing in the winter season on this day. It is in winter that winds driven by east-moving midlatitude cyclones frequently encounter colder outflow winds from the interior of British Columbia through the valley. This interaction reduces the organized diurnal up- and down-valley circulations that dominate during summer. Once exiting the valley into the Douglas Channel, wind vectors are more aligned with the fjord, consistent with the expectation of lesser drag over water surfaces.
b. Quantitative and performance evaluations for surface variables
1) Ambient temperature
Near-surface temperature (T2) is an important meteorological parameter as it affects the buoyancy of air pollutants and reaction rates of chemical species. Figure 4 shows diurnal sequences of mean hourly bias for summer and winter seasons. In summer, a warm bias occurs after sunrise, peaking at about 0800 h, becoming a cool bias after noontime, and reaches a negative peak at sunset. The adaptation to convective transport by the nonlocal SH and YSU PBL schemes, the bias profiles of which are virtually the same, is evident in larger daytime bias amplitudes than those of the other schemes. At nighttime, YSU, SH, and ACM2 have marginally colder biases than the others, leading to relatively cooler temperatures over the course of a single day with these schemes.
Table C1 in appendix C tabulates model performance across four seasons in a year. Going from the spring to winter period, T2 simulations change from negative biases to predominantly positive biases. Biases are lower during the summer and winter seasons, compared to the spring and autumn seasons. Over land, diurnal MB ranged from −2.5 to 1.3 K in spring, from −0.6 to −0.2 K in summer, from 0.5 to 4.1 K in autumn, and from −0.3 to 0.9 K in winter. Model errors for T2 are low in summer and winter months, considering that RMSE values at all three land stations meet the criteria in these seasons. Temporal correlations are also very good, and the r benchmark is achieved in 100% of evaluations from spring through autumn (>0.9 in summer) and 67% for winter. Few evaluations meet IOA criteria; however, values for IOA are modest, since across seasons the majority are within 0.60–0.79. Results from all statistical measures point that T2 simulations over land are poorest in autumn. Biases and errors for the Nanakwa station in spring (MB range: from −0.5 to 0.1 K and RMSE range: 1.25–1.35 K) and summer (MB range: from −1.0 to −0.4 K and RMSE range: 1.52–1.72 K) suggest reasonable accuracy of interpolated sea surface temperatures in warmer months. In terms of individual PBL schemes, MYJ marginally outperforms UW with regard to accuracy and precision indices (MB, RMSE, and IOA), and MYNN3 performs best for pattern replication (r). The MYJ scheme also has several satisfactory correlation scores, thus demonstrating fitness for various statistical measures.
2) Specific humidity
Another important meteorological variable from an air quality perspective, especially in controlling aerosol formation and transformation processes, is specific humidity (water vapor content). Performance statistics for specific humidity at 2-m height (Q2) for each PBL scheme per season and at two stations are presented in Table C2 of appendix C. Across PBL schemes and relative to observations, the dominant pattern is a dry bias. Biases are larger at the Onion Lake station, which receives more precipitation than the Terrace station (https://data.pacificclimate.org/portal/pcds/map/) and in summer when more evapotranspiration occurs than in winter. These further attest to the dynamical consistency of model outputs. At the Onion Lake station, for instance, the MB range is from −0.50 to −0.37 g kg−1 in winter, as compared with from −0.68 to −0.19 g kg−1 in summer. Overall, model performance for specific humidity is remarkable; 94% of evaluations for MB, 100% for RMSE, 56% for r, and 65% of IOA fall within respective cutoff values of satisfactory performance. With regard to individual PBL schemes, the MYNN3 scheme has the best fit for Q2 simulation based on the number of top-ranked and satisfactory evaluations for it. Indeed, because all schemes underestimate Q2, the MYNN3 scheme best ameliorates negative moisture biases.
Surface wind impacts the mixing and transport of air contaminants from source points/areas to receptors. Biases in modeled wind direction (WDIR10) at the different stations are presented in Table C3 of appendix C. The localized nature of wind over complex terrain is attested by biases at Onion Lake, which are anticlockwise, unlike the other stations. Most remarkable is that at the Nanakwa buoy, almost all evaluations meet the ±10° criteria for wind direction bias. On land, model skill is reduced, arising from increased surface roughness, and wind flow being affected by multiple terrain elements. Only 8% of evaluations for land stations are within the bias threshold. However, 75% of all evaluations for Terrace are ≤ 20° and all for Terrace and Whitesail are ≤ 30°. The largest biases are at the Onion Lake station, in the central portion of the valley (up to ~−60°). Using the same rule for assessing suitability, the MYNN3 and UW PBL schemes are top ranked for wind direction.
Wind speed uncertainty due to PBL scheme choices is also of interest and Fig. 5 shows diurnal sequences of mean hourly biases for 10-m wind speed (WSPD10) during the summer and winter seasons. Differences between schemes rise to as much as 1.5 m s−1 in summer but are smaller (~1 m s−1) in winter. The valley atmosphere is more prone to decoupling the surface layer from layers aloft in winter due to increased stability, thereby lessening the sensitivity of surface winds to PBL parameterizations. Strong synoptic forcing from landfalling Pacific midlatitude cyclones would account for why modeled WSPD10 has greater biases in winter than in summer. It is reasoned that the general overestimation of WSPD10 is due to WRF not sufficiently weakening winter storms as they transition from the smooth ocean to the rugged interior. Indeed, wind speed biases increase from the coast toward the northern part of the valley. Aside from complications of cross-flow of the Skeena River (Fig. 1) and urban obstacles that attenuate real surface winds in Terrace, the general area to the north is dominated by mountainous topography. Wind speed overestimation has been reported in several experiments using WRF at similar horizontal grid resolutions over complex terrain (e.g., Avolio et al. 2017; Hariprasad et al. 2014; Zhang et al. 2013).
The performance statistics for WSPD10 for all four seasons for each of the PBL schemes are in Table C4 of appendix C. Overestimations occur not only in winter and summer months but also through spring and autumn. Across land stations, diurnal MB ranged from 0.4 to 2.7 m s−1 in spring, from −0.2 to −2.5 m s−1 in summer, from 1.3 to 3.0 m s−1 in autumn, and from 1.3 to 3.7 m s−1 in winter. Just 10 of 72 evaluations (14%) are within the MB benchmark for wind speed. Model errors for WSPD10 are also higher in autumn and winter months, although the proportion of evaluations within the RMSE criteria is much improved (68%). Favorable RMSE and r scores are in contrast with negative 1OA values. Whereas 64% of evaluations for r are ≥ 0.6, only 11%—all for Onion Lake—meet the IOA benchmark. Thus, as compared with air temperature and humidity, WSPD10 is poorly represented, especially around Terrace.
c. Spatial-difference examinations of surface meteorological fields
Because diurnal surface heating and cooling impact atmospheric structure, spatial differences in model outputs are examined for the relative effect of individual WRF PBL schemes on meteorological properties at nighttime and daylight periods. The seasonal attributes of the adjustment of PBL schemes to local weather and their capacities for mixing of surface fluxes are also qualified.
Some turbulence closure approaches are conditionally operable. The ACM2 scheme deactivates nonlocal transport for stable conditions, which would prevail at night. Of interest, therefore, is how meteorological quantities for periods when only the local closure component of this scheme operates compare to those from the fully nonlocal YSU scheme, given that both PBL schemes are mapped to the same surface-layer model. Figure 6a shows that over the valley portion of the domain, the skin temperature (SKINT), T2, and WSPD10 from the ACM2 scheme are greater than those of the YSU scheme. At higher elevations outside the valley channel, the topography is less confined, and the YSU scheme is more capable of moving cold air from the surface to layers aloft and downward momentum transport. The relative weakening of nonlocal transport over the valley, therefore, derives from the confinement of the land surface by steep mountains and is also due to nighttime stable conditions that restrain vertical transport through the entire PBL. In the corresponding wintertime, when the atmosphere is more stable, the lower temperature by the YSU scheme (~0.5°C difference), relative to that from the ACM2 scheme, over the valley is sharper but less widespread, since the severest restriction to nonlocal flux exchange would be at the deepest of cold-air pools.
Since Q2 is the moisture content of surface air and is therefore analogous to pollutant concentration at ground level, its comparison across schemes could be useful as a proxy for contaminant accumulation and dispersion. Differences in Q2 may be linked to relative capacities of PBL schemes to facilitate vertical transfers between model layers. It is noteworthy that patterns of ACM2/YSU paired differences for PBL height (PBLH) and Q2 are similar (Fig. 6a). The PBL simulated by the ACM2 scheme is thicker and its Q2 values are moister over water bodies, and is thinner and its Q2 values are drier over land, when compared with those values from the YSU scheme. In the YSU scheme, vertical diffusion is constrained to the diagnosed PBLH. Thus, it is less diffusive of Q2 over the land areas. Over water, Q2 from the YSU scheme has a steeper vertical gradient between the PBL and the overlying model layers, enabling faster vertical transport of water vapor to the layers aloft.
The constraining of nighttime mixing by steep topography is less with local turbulence closure. The stable stratification that develops leads to predominantly small eddies that do not span several atmospheric layers, hence the actual momentum and heat transfers could resemble the transport between adjacent vertical levels in WRF by local PBL schemes. However, various local mixing approaches exist whose expressions may also differ with alternative surface-layer parameterizations. Figure 6b compares the ACM2 scheme nighttime outputs with those of the local MYJ scheme that is mapped to the Eta Model surface-layer scheme. Note that, unlike the MYJ scheme, the local closure implementation in the ACM2 scheme ignores TKE transport. Hence in the winter season when temperature inversions are frequent, along-valley PBLH simulated by the ACM2 scheme is shallower than that of the MYJ scheme. However, the as much as 1.0°C difference for T2 over the TKV suggests surface-layer scheme influence. Over the valley, MYJ outputs warmer T2 than ACM2 but cooler SKINT, hinting at greater heat flux retention with the Eta Model surface-layer scheme than with the MM5 surface-layer scheme. Unlike the MM5 surface-layer scheme, the Eta Model surface-layer scheme considers thermal roughness length (Janjić 2002), apparently contributing to warmer ambient temperatures with the MYJ PBL scheme. Differences in friction velocities (not shown) were similar to patterns for wind speeds, for which the MYJ scheme values over the valley and coastal channels are less. Friction velocity is a fraction of wind speed (Camuffo 2014). The 0.1 m s−1 friction velocity limit in previous versions of the MM5 surface-layer scheme was considered as preventing undesirable effects such as excessive ground cooling (Jiménez et al. 2012). The revised MM5 surface-layer scheme as used in the present study, permits friction velocities as low as 0.01 m s−1 (Jiménez et al. 2012), while in the Eta Model surface-layer scheme, a correction is applied so that only nonzero values are obtained (Janjić 2002). Consequently, lower friction velocity over the valley and coastal channels with the MYJ scheme possibly derives from its coupling with the Eta Model surface-layer scheme, which in turn, could offset part of the warming effect of thermal roughness length. Cooler SKINT of the MYJ output relative to that of the ACM2 simulation over the valley channel may have resulted from such differences in friction velocities.
For other model runs with the Eta Model surface-layer scheme (Figs. 6c,d), differences among PBL schemes for simulated variables are less than in Fig. 6b, where both the surface-layer scheme and PBL schemes are varied. Nonetheless, the differential subtleties of these formulations are worth examining, given that PBL scheme development, testing, and revising have primarily centered on weather prediction rather than on air pollution modeling. Warmer skin temperatures for MYNN3 relative to MYJ translate to a higher T2 but do not result in greater water vapor retention overall. MYNN3’s summertime Q2 over the southern portion of the domain is drier than that of the MYJ scheme. However, this portion roughly coincides with the parts where the height of MYNN3’s PBL surpasses that of the MYJ scheme (also Fig. D1 in appendix D), pointing to greater mixing with the MYNN3 scheme. This seems to affect nighttime WSPD10 as well, with greater or lesser values for the MYNN3 scheme coastward or landward, respectively. In addition to implementing local eddy diffusion, the MYNN3 scheme, unlike the MYJ scheme, has a nonlocal mixing component. The suppression of nonlocal vertical exchange in the MYNN3 scheme relative to the MYJ scheme is therefore less than for the YSU scheme in relation to the ACM2 scheme in Fig. 6a. There is also a seasonal signal with this behavior. In summer, modeled onshore winds generate instabilities over considerable distances upon which nonlocal mixing is enhanced. In wintertime, when valley air is colder and drier, and simulated circulations are northerly, MYNN3’s nonlocal mixing is restrained: the mixing is more diffusive than in the MYJ scheme only over the Douglas channel that is exposed to warmer, maritime air.
The UW scheme, designed for simulation of marine stratocumulus-capped boundary layers, neglects TKE storage and assumes no TKE transport occurs during stable conditions (Bretherton and Park 2009), hence it simulates the lowest PBL heights and cloud bases that are nearest to the ground (Fig. D1 in appendix D). Some modeling experiments in mountainous regions (e.g., Balzarini et al. 2014) have also reported lower PBL heights with the UW scheme. Its output of lower Q2 than that of the MYJ scheme (not shown) suggests a relatively rapid moisture loss mechanism. It is worth stating that the UW PBL scheme’s verification and integration into the WRF Model mentioned the diffusion of humidity and chemical species from the lowest model grid level, rather than their mixing across layers (Bretherton and Park 2009). Because air pollutants are also mass fluxes, less nighttime ambient moisture with the UW scheme relative to the MY schemes indicate the likelihood of lower air pollutant concentrations.
PBL schemes adjust differently to enhanced convection in the daytime, also impacting meteorological fields over complex terrain. Comparisons between the ACM2 and YSU schemes in summer (Fig. 7a) indicate that for the most part, the ACM2 scheme produces a higher PBL over land areas, in reverse of the nighttime pattern. Unlike single estimation of the PBL top in the YSU scheme, the PBLH in the ACM2 scheme is calculated from separate diagnosis of convective and entrainment layers (Pleim 2007), partly explaining PBLH differences. Convective regimes activate the nonlocal transport component of ACM2, thereby providing a hybrid vertical flux transfer mechanism, as opposed to the YSU scheme that is nonlocal. For parts of the valley adjacent to the marine channel, incursions of maritime air may diminish the development of convective regimes more than at positions farther inland; hence, there are small portions where the ACM2 scheme’s Q2 output surpasses that of the YSU scheme. WSPD10 differences along the Douglas Channel show the ACM2 output to be the windier, also contributing to moister inflow. However, for the most part, the ACM2 scheme generates less Q2 in the summer daytime (see also Fig. 8), which suggests stronger updrafts with the ACM2 scheme than those produced by the YSU scheme in this period. The ACM2 scheme is thus more likely to produce less pollutant concentration at the surface than the YSU scheme, especially in the drier, northern part of the valley. In the cold winter season (Fig. 7a), the daytime Q2 difference is negligible (≤0.1 g kg−1) and the change of water vapor with height (Fig. 8) is comparable between the two PBL schemes. This is because the PBL heights simulated by the ACM2 scheme are mainly lower than those of the YSU scheme, and its nonlocal transport component is inactive in stable/neutral conditions that dominate in winter.
Figure 7b compares functionally similar PBL schemes for daytime mixing processes but mapped to different surface-layer schemes. The MYNN3 scheme, with relatively positive incoming radiation over the valley in summer, simulates less SKINT than the ACM2 scheme. However, over the mountains, T2 differences, even in wintertime, are more positive for the MYNN3 scheme than over the corresponding areas in the SKINT plot. This again could be due to greater sensible heat flux from the Eta Model surface-layer scheme. Apart from differences in surface-layer formulation, dissimilarity in the implementation of hybrid mixing methods would also account for differences in simulated quantities. The K-theory-based ACM2 scheme uses first-order statistics to approximate higher moment terms, while the MYNN3 scheme is a TKE-based second-order closure. The PBLH diagnosed by the ACM2 scheme is predominantly thinner (thicker) in summer (winter) than the MYNN3. Adjustment of vertical gradients of momentum and moisture fluxes with PBLH changes in the ACM2 would also affect WSPD10 and Q2, respectively. The MYNN3 scheme mostly outputs greater Q2 and WSPD10 than the ACM2 scheme but the extent to which PBL parameterization alone accounts for this outcome is confounded by the use of alternative surface-layer schemes.
Simulation differences among the TKE-based PBL schemes (Figs. 7c,d) show the influence of local and nonlocal closures with the Eta Model surface-layer scheme. PBLH is greater with MYJ than with MYNN3 except near the coast. As previously stated, the daytime intrusion of marine air counteracts growing instability arising from solar heating. For this reason, higher PBLH with the MYNN3 scheme relative to that of the MYJ scheme does not extend as far inland as it does at night (Fig. D1 in appendix D). But instead of a moisture profile that is more vertically homogeneous than those of the MYJ and UW schemes, the MYNN3 scheme indicates less mixing (Fig. 8) due to how it implements nonlocal transport. The MYNN3 scheme comprises a partial-condensation model for moist convection, in which the buoyancy effects of subgrid-scale clouds are considered using countergradient terms (Nakanishi and Niino 2009). Countergradient diffusion process in the MYNN3 scheme thus explains why it has greater Q2 and lesser WSPD10 values than the outputs from the MYJ and UW schemes.
Recent revisions to mixing-length formulations in the MYNN3 scheme for which the current code further decreases the downward mixing of momentum (Olson et al. 2019), would also contribute to its generation of WSPD10 and Q2 values that are less biased than the UW and MYJ schemes. Because of increased linking of the MYNN3 scheme to cloud processes and to other physics modules, such as radiation schemes beginning from WRF, version 3.8, feedbacks from simulated boundary layer clouds are probably more influential on surface meteorological variables than with other PBL schemes. For instance, the UW scheme outputs greater cloud coverage than the MYNN3 scheme (Fig. D1 in appendix D), which ought to result in less radiation reaching the surface with the UW scheme. Instead, incoming solar radiation is smaller with the MYNN3 scheme, consequently the lower T2. The consideration of moist convection by the MYNN3 scheme is also evident in it having the biggest gradient for daytime potential temperature (Fig. 9). It is reasoned that the facilitation of cloud condensation in the MYNN3 scheme, enabling latent release and warming of the air aloft, alongside less solar radiation to the surface and cooler near-surface temperature, causes a stronger temperature gradient. The least vertically homogenized moisture gradient within the lowest 500 m for the MYNN3 scheme, may partly be due to this stability-enforcing feedback. In winter when absolute water vapor content is lower and convective activity is reduced, MYNN3’s Q2 surplus compared to those of MYJ and UW in the valley is diminished or in deficit. The vertical moisture gradient at this time is slightly steeper for the MYJ scheme than for the MYNN3 scheme. The above analysis points to the seasonal modulation of MYNN3’s nonlocal mixing.
4. Discussion and conclusions
Performance scores of the various PBL schemes for surface meteorological parameters are generally within or close to values that are found in the literature, either as NWP experiments (García-Díez et al. 2013; Horvath et al. 2012; Hu et al. 2010; Zhang et al. 2013) or in the validation of meteorological inputs to air quality models in mountainous regions (Cheng et al. 2012; Hedley and Singleton 1997; Mues et al. 2018). It is not unusual for performance to vary between models, model versions, grid resolutions, and episode durations, and performance is also subject to operating environments and needs. As an example, the evaluation benchmarks in Emery et al. (2001) were based on MM5 model application at 4-km grids to two ozone episodes, each lasting a week, in Texas. Further, there usually is a minimum amount of scatter or representativeness error that is inherent in meteorological downscaling and cannot be removed. Rife et al. (2004) estimated representativeness error in a model grid box of size 1.33 km × 1.33 km to be about 1 m s−1 for near-surface wind speed in a well-mixed boundary layer over complex terrain. Therefore, error and bias estimations previously presented reflect the status of the WRF Model (version 4.0) and ERA5 to simulate surface variables in the TKV without observation nudging.
As a result of atmospheric circulations being constrained by steep topography, the performances of the PBL schemes for wind direction were reasonably comparable. Therefore, the choice of one scheme over another may not be critical to representing the relative contributions of wind sectors to pollutant concentrations at specific locations. Performances for surface air temperature and wind speeds, on the other hand, were more varied. Overpredicted wind speeds in particular could cause an underestimation of ambient pollutant concentrations in areas near major emissions sources. Indeed, the finding that PBL schemes were able to simulate surface wind with sufficient directional accuracy, yet poorly represent the observed speed, summarizes the challenge of providing reliable winds fields over complex terrains for input to coupled numerical models. As was noted previously, dynamical downscaling aligned gridded coarse meteorological output to the valley and fjord channels—hence the general correspondence of modeled wind direction with observations—whereas ambient wind speed biases reflected uncertainties that can arise from smoothing of rugged relief in the simulations.
Model terrain height errors were from a few to several dozen meters (Fig. D1 in appendix D), and terrain data quality and resolutions may have affected individual simulations differently. It was noticed, however, that performances for the YSU and SH schemes were generally less distinct ( appendix E), relative to that between the YSU scheme and the ACM2 scheme. Hence differences in concentrations of air pollutants that may arise from selecting the YSU scheme in place of the SH scheme would probably be negligible. Both schemes are founded on nonlocal K-profile assumptions, but the SH PBL scheme additionally addresses turbulent transport in the “gray zone” (Wyngaard 2004), which is the range of grid spacing below which parameterizations have been developed but above which the relevant processes are entirely resolved. Since simulations at 1-km horizontal grid spacing as used in the present study may be coarser than gray-zone resolution, the grid-size dependency of subgrid-scale, nonlocal heat and momentum transport that is unique to the SH scheme was probably not optimized. In a real-case verification of the SH scheme by Shin and Hong (2015), the simulation of convective rolls was structurally more distinct from those of the YSU scheme at less than 1-km grid spacing. Because of the 0.281° spatial resolution of the ERA5 forcing dataset, implementing a smaller grid size for differential modeling over the TKV was not done since that would have resulted in more computing costs along with sharp elevation gradients causing numerical instabilities in the simulations. Further, convective periods in the TKV are restricted to daytime in the summer and involve moist turbulence that is not represented by the SH scheme.
Distinctions in model performance generally evolved with season changes. In summer when Pacific frontal passages are rare and local thermally driven circulations are enhanced, the diurnal signatures of individual PBL schemes were fairly unique. Particularly for air pollutants that peak in warm season such as ozone, the selection of one PBL scheme over another may be significant. In winter, synoptic forcing is stronger, the daylight period is shorter and diurnal variability is weaker than in summer. More extended decoupling of the atmosphere in wintertime diminished the influence of model configurations on near-surface variables, and all schemes performed poorly in this period. But because atmospheric phenomena such as temperature inversions are not always surface based, and in the absence of radiosonde measurements, ensemble projections of the air quality in the cold season rather than the output from just one PBL scheme may be better. The UW PBL scheme that consistently produced the shallowest boundary layer heights may need further investigation to know how realistic such representations are for air quality parameters in the region.
In essence, PBL schemes are highly integrated with other physics and dynamics options. It needs stating that some of the PBL schemes that were tested in this work can be mapped to other surface-layer models in WRF. Specifically, the ACM2 PBL scheme can also be mapped to the Pleim–Xiu surface-layer scheme, and the MM5 surface-layer scheme can be used for both the MYNN3 and UW PBL schemes. As an example, relative to the Eta Model surface-layer scheme, MYNN3 mapped to the MM5 surface-layer scheme produced larger negative biases for temperature, specific humidity, and wind speed (Fig. 10). Where a turbulence closure scheme can link to more than one surface-layer model option, the practice was to map TKE-based PBL schemes to the Eta Model surface-layer scheme and first-order closure schemes to the MM5 surface-layer scheme. Figure 10 partially strengthens this practice.
This study refrained from comparing modeled values with domain-wide averaged observation data. One shortcoming in several past studies was limited objective verification of target meteorological parameters through time and nonclarity as to what qualifies as good performance of mesoscale models in regions of strong surface heterogeneity. Our approach has been to assess performance at discrete locations and to supplement conventional metrics such as MB, RMSE, and r scores with a modified IOA (bounded between 1 and −1) that distinguishes the strength of model predicted variability relative to observed variability. As an example, summertime RMSE and r values for wind speed at the Onion Lake and Whitesail stations were comparable and within performance benchmarks. However, while IOA values for Onion Lake were good, IOA for Whitesail was poor, which was attributable to errors (overestimation) that are much larger than the variability in observations. This is despite the fact that the Whitesail station is located nearer the shoreline and should ideally match better with strong breezes that were simulated. It is important to consider such error disparities in downscaled meteorological fields over complex coastal areas when simulating long-term ambient chemical concentrations.
This study was supported by the Natural Sciences and Engineering Research Council of Canada Discovery Grant RGPIN-2017-05784 to Peter Jackson. We acknowledge the support of high-performance computing resources from the University of Northern British Columbia. Special thanks are given to the Pacific Climate Impacts Consortium for providing access to meteorological data from their archive. We thank two anonymous reviewers whose comments significantly improved this paper.
Formulas for Conventional Statistical Measures
Below are the formulas for MB, r, RMSE, and IOA:
where Pi and Oi are pairwise modeled and observed values, respectively, n is the number of pairs, and terms with overbars are respective mean values.
Verification of Regional Atmospheric Soundings
Figure B1 shows the similarity of regional atmospheric soundings simulated with the various PBL schemes with the observed soundings for different times in winter and summer at the Annette Island radiosonde station.
Performance Statistics for Surface Meteorological Variables
This appendix gives performance statistics for surface air temperature at diurnal time scale (Table C1), specific humidity at diurnal time scale (Table C2), wind direction bias using all hourly values (Table C3), and wind speed at diurnal time scale (Table C4) for each season.
Valley Centerline Profiles of Model-Diagnosed Variables
Figure D1 shows profiles of WRF Model–diagnosed variables along the valley and coastal inlet centerlines.
Comparison of SH and YSU PBL Schemes for Surface Meteorological Variables
At 1-km horizontal resolution, the nonlocal transport of subgrid-scale turbulence in convective boundary layers dominates total transport (Shin and Hong 2015). Since the YSU scheme is designed for parameterizing the total heat transport regardless of grid size, there is little difference (±0.1°C; Fig. E1) between its T2 and that of the SH scheme that includes grid-size dependency function (cf. Figs. 6a and 7a for the differences with the ACM2 PBL scheme that is mapped to the same surface-layer scheme).
Differences in WSPD10 and Q2 might also not mainly be from scale-dependent functions in the SH scheme, especially considering the similarity of PBL heights (PBLH difference is approximately ±30 m). In this study, simulations with the YSU PBL scheme were run with topographic correction for surface winds (topo_wind: = 1 in WRF; Jiménez and Dudhia 2012; Skamarock et al. 2019) to represent extra drag from subgrid topography and enhanced flow at hilltops. This unique setting for the YSU PBL scheme must have caused wind speeds that are less than that of the SH scheme and hence slightly lower bias with observation data overall (cf. mean biases across stations in Table C4 of appendix C). In the daytime of summer, the stronger winds of the SH scheme can transport more maritime air and consequently more moisture content with this scheme than the YSU scheme. At night when the wind speed differences are smaller and real winds are weaker, the SH scheme outputs less Q2 than the YSU scheme.