To accurately simulate the atmospheric state within the planetary boundary layer (PBL) by PBL parameterization scheme in different regions with their dominant weather/climate regimes is important for global/regional atmospheric models. In this study, we introduce the turbulence kinetic energy (TKE) and TKE dissipation rate (ε) based 1.5-order closure PBL parameterization (E–ε, EEPS) in the Weather Research and Forecasting (WRF) Model. The performances of the newly implemented EEPS scheme and the existing Yonsei University (YSU) scheme, the University of Washington (UW) scheme, and Mellor–Yamada–Nakanishi–Niino (MYNN) scheme are evaluated over the stratocumulus dominated southeast Pacific (SEP) and over the Southern Great Plains (SGP) where strong PBL diurnal variation is common. The simulations by these PBL parameterizations are compared with various observations from two field campaigns: the Variability of American Monsoon Systems Project (VAMOS) Ocean–Cloud–Atmosphere–Land Study (VOCALS) in 2008 over the SEP and the Land–Atmosphere Feedback Experiment (LAFE) in 2017 over the SGP. Results show that the EEPS and YSU schemes perform comparably over both regions, while the MYNN scheme performs differently in many aspects, especially over the SEP. The EEPS (MYNN) scheme slightly (significantly) underestimates liquid water path over the SEP. Compared with observations, the UW scheme produces the best PBL height over the SEP. The MYNN produces too high PBL height over the western part of the SEP while both the YSU and EEPS schemes produce too low PBL and cloud-top heights. The differences among the PBL schemes in simulating the PBL features over the SGP are relatively small.
The commonly used planetary boundary layer (PBL) parameterizations in atmospheric models can be divided into two categories: the traditional K-theory (eddy-diffusivity) parameterization and the eddy-diffusivity mass-flux (EDMF) approach. The commonly used K-theory parameterization includes first-order and 1.5-order closures (Orlanski et al. 1974; Pielke and Mahrer 1975; Blackadar 1979; Wyngaard 1982; Yamada and Kao 1986). The EDMF approach includes the K-theory turbulence closure to parameterize the local transport by small turbulent eddies and the mass-flux part to represent the nonlocal boundary layer organized eddy fluxes (Köhler et al. 2011; Hourdin et al. 2013; Han et al. 2016).
The first-order K-theory turbulence closure includes many widely used PBL schemes, for example, the medium-range forecast scheme (MRF; Hong and Pan 1996), the Yonsei University scheme (YSU; Hong et al. 2006), the Atmospheric Convective Model version 2 (ACM2; Pleim 2007), and so on. The first-order closure diagnoses the vertical mixing coefficient K from the local Richardson number or a diagnostic mixing length (l) or both (Blackadar 1962; Bodin 1980; Carlson and Foster 1986). The vertical mixing coefficients for heat/moisture (Kh) and momentum (Km) are often different and are treated with different mixing lengths. The commonly used 1.5-order closure involves a prognostic equation for turbulence kinetic energy TKE (or E) and a mixing length for calculating the vertical mixing coefficient K. The equation for TKE includes buoyancy, local wind shear, vertical transport, and TKE dissipation rate ε. The latter is often determined by a diagnostic equation with E and dissipation mixing length (lε). A common assumption is that the mixing lengths l and lε are equal, although Therry and Lacarre’re (1983) derived a relationship between l and lε to show that l is larger than lε for convective conditions. The mixing length l can be prognostic (Mellor and Yamada 1974, 1982) or diagnostic (Miyakoda and Sirutis 1977; Janjić 1990). The combination of the prognostic E and diagnostic (or prognostic) l is often called an E–l parameterization. The most widely used E–l parameterization is the Mellor–Yamada–Janjić scheme (Janjić 1994) and its improved version, the Mellor–Yamada–Nakanishi–Niino scheme (MYNN; Nakanishi and Niino 2006).
An alternative method to determine K in a 1.5-order closure is to use full prognostic equations for both E and ε (Detering and Etling 1985; Duynkerke and Driedonks 1987; Langland and Liou 1996), often called E–ε parameterization scheme (EEPS). The E–ε scheme has proven to be skillful in predicting atmospheric boundary layer features, including the turbulent structure of the stratocumulus-topped atmospheric boundary layer (Duynkerke and Driedonks 1987), the maritime boundary layer (Wang et al. 2004a,b), the mesoscale features in association with coastal fronts (Holt et al. 1990), maritime low-level jets (Gerber et al. 1989), and the evolution of convective boundary layer in a weather forecast model (Langland and Liou 1996). Beljaars et al. (1987) compared the E–l and E–ε schemes and found that the E–ε scheme can better retain “memory effects” in velocity scales when surface condition changes. Wang (2001, 2002) implemented an E–ε scheme into his tropical cyclone model, which has been successfully used to study various aspects of tropical cyclones (e.g., Wang 2001, 2007, 2008a,b, 2009; Wang and Xu 2010). This tropical cyclone model with the E–ε scheme was later further developed into a regional climate model at the International Pacific Research Center (iRAM), which has been shown to perform well in simulating East Asian summer monsoon (Wang et al. 2003; Sen et al. 2004; Xie et al. 2006; Souma and Wang 2010; Qi and Wang 2012), regional climate over the southeast Pacific (Wang et al. 2004a,b; Xu et al. 2004; Wang et al. 2005; Xie et al. 2007; Wang et al. 2011; Lauer et al. 2012), marine boundary layer clouds in the Okhotsk Sea in summer (Koseki et al. 2012), and precipitation diurnal cycle in tropics (Zhou and Wang 2006; Wang et al. 2007).
Recently, we have implemented the latest version of the E–ε turbulence scheme used in iRAM into the Weather Research and Forecasting (WRF) Model. The purpose of this paper is to document this scheme in the WRF Model and compare this scheme with three other widely used PBL schemes: the latest scale-aware version of the MYNN scheme in the EDMF framework (Olson et al. 2019), the scale-aware version of the YSU scheme (Shin and Hong 2013, 2015), and the University of Washington scheme (UW; Bretherton and Park 2009), in simulating the marine boundary layer over the southeast Pacific (SEP) and the diurnal cycle in the boundary layer over the Southern Great Plains (SGP) of the contiguous United States. The rest of the paper is organized as follows. Section 2 provides a description of the three PBL schemes, particularly the E–ε scheme. Section 3 describes data from two field campaigns and the experimental setup. The simulation results are discussed in section 4. Main conclusions are summarized and discussed in the last section.
2. Description of the four PBL schemes
a. The EEPS scheme
The EEPS scheme we implemented in WRF is similar to that in Langland and Liou (1996) with a few modifications and improvements, such as changing the coefficients for ε equation, selecting maximum of shear production versus the sum of the shear and buoyancy productions in ε to avoid oscillation, enhancing the buoyancy term in E and ε equations when there are clouds, and taking TKE dissipation rate ε as an additional heat source (dissipative heating) of the atmosphere. The equations for TKE E and its dissipation rate ε are (Holt et al. 1990)
where Km and Kh are vertical mixing coefficients for momentum and heat/moisture, respectively; u and υ are zonal and meridional wind speeds, respectively; z is height; t is time; C1, C3, C4, and C5 are constants being taken as 1.35, 1.44, 1.92, and 0.77, as suggested by Duynkerke and Driedonks (1987); Pbuoy is the buoyancy, defined later in Eqs. (3)–(8). Four terms on the rhs of Eq. (1) represents shear production, buoyancy production/destruction, vertical transport, and dissipation of E. Three terms on the rhs of Eq. (2) represent the generation, destruction, and vertical transport of ε. Following Duynkerke and Driedonks (1987), the generation term in the angle brackets is taken as the maximum of shear production versus the sum of the shear and buoyancy productions. The minimum allowed values for E and ε are set to be 1 × 10−4 m2 s−2 and 1 × 10−6 m2 s−3, respectively. The background exchange coefficient is 0.1 for momentum and 0.01 for heat.
To better represent the mixing in clouds, the buoyancy production term in both the TKE and its dissipation rate equations are calculated by the formula suggested by Durran and Klemp [1982, see their Eq. (36)] in clouds where the air is saturated with nonzero cloud water mass or cloud ice mass. This buoyancy production calculation plays a role in enhancing the vertical mixing in clouds (Rotunno and Emanuel 1987; Tripoli 1992). Therefore, in the unsaturated region
where θυ is the virtual potential temperature and g is the gravity acceleration, and in clouds
where qt is the sum of the mixing ratios of cloud water, rainwater, cloud ice, snow, graupel and hail, if any of them is available; θ is potential temperature; Cp is specific heat at constant pressure; R is the gas constant for dry air; Rυ is the gas constant for water vapor; T is air temperature; T0 is air temperature at freezing point; Lυ is latent heat of condensation; and Ls is latent heat of sublimation. Note that the above modification to the original expression given by Durran and Klemp (1982) is important in enhancing the vertical mixing in ice clouds (Tripoli 1992).
The upper boundary conditions at the model top are E = 0 and ε = 0. The lower boundary conditions at the lowest model level are (Mailhot and Benoit 1982)
where zkm is the height of the lowest model level; u* is the friction velocity; w* is the convective velocity scale; L is the Monin–Obukhov length; k (=0.4) is von Kármán constant; Ris is the surface-layer Richardson number, which is defined as
where θυs is the surface virtual potential temperature; θυkm, ukm, and υkm are the virtual potential temperature, zonal wind, and meridional wind at the lowest model level, respectively.
The vertical mixing coefficient for momentum Km in the E–ε scheme is calculated by (Holt and Raman 1988)
where the Richardson number Ri is defined as
The time integration of Eqs. (1) and (2) is similar to that used by Langland and Liou (1996). To avoid large time truncation errors due to short time scales of E and ε, the solutions for E and ε are obtained with multiple time integrations with smaller time steps. Both E and ε are treated as tracers in the model dynamics so that three-dimensional advection and horizontal diffusion can be applied to both. Searching one model level by one model level upward from the level nearest to surface, the PBL height is defined as the lowest model level at which the modified bulk (boundary layer) Richardson number exceeds a critical value of 0.25 (Troen and Mahrt 1986).
A nonlocal term γ is included in vertical eddy mixing of potential temperature (γt) and moisture (γq), respectively,
where h is the PBL height, C = 5.0, Q1 is surface heat flux, and Q2 is surface moisture flux. The TKE dissipation rate ε is an additional heat source of the atmosphere (Han et al. 2016), and it can be expressed as
where T is air temperature and qυ is water vapor mixing ratio.
The prognostic equation for vertical mixing can be represented as follows:
where C represents potential temperature or moisture.
One of the problems we have encountered during the implementation of the EEPS scheme into the WRF Model is the oscillation of solution when the ε production term is very small. The advection of TKE and ε can considerably alleviate the oscillation. We also find that the coefficients in Eq. (2) used in Langland and Liou (1996) lead to strong oscillation in the nocturnal stable boundary layer. Since the coefficients used by Duynkerke and Driedonks (1987) can significantly suppress such an oscillation, we used their coefficients. The implemented version of the EEPS in the WRF Model performs quite well without significant oscillations with the above mentioned modifications. Nevertheless, there is still room to further tune C3 and C4 in Eq. (2) (Freedman and Jacobson 2003).
b. The YSU, MYNN, and UW schemes
where is the flux at the top of the boundary layer, and presents an asymptotic entrainment flux term. Shin and Hong (2015) revised γ and K in Eq. (20) to make the scheme scale aware. The YSU scheme follows Troen and Mahrt (1986) to calculate the PBL height. The critical bulk Richardson number is 0.25, except over water where the critical bulk Richardson number is defined as a function of wind speed at 10-m height and surface roughness length, but not larger than 0.3.
where M is mass flux, Cu represents prognostic variables in plumes, and represents the variables for environment. The default MYNN scheme applies multiple plumes to account for the mass flux [i.e., ] (Sušelj et al. 2012; Olsen et al. 2019). The maximum number of plumes is allowed up to 10 in the MYNN scheme. The mass flux part is active when 1) at least one plume is defined, 2) there must be a positive surface buoyancy flux, and 3) the lowest 50 m must be superadiabatic. A K-theory turbulence closure is used everywhere. In the MYNN scheme, eddy-diffusivity is parameterized with E–l 1.5-order closure. The scale-aware factors are applied to vertical mixing coefficient K and mass flux (Honnert et al. 2011; Shin and Hong 2013; Olson et al. 2019). The MYNN scheme employs a hybrid PBL height diagnostic method. Here, the height where the potential temperature first exceeds its minimum value in the PBL by a fixed value (Nielsen-Gammon et al. 2008), and TKE does not exceed a threshold value (Pichugina et al. 2008), are blended using a hyperbolic tangent.
The UW scheme (Bretherton and Park 2009) is derived from the Grenier–Bretherton scheme (Grenier and Bretherton 2001) that is designed particularly for a PBL influenced by stratocumulus clouds. One distinct feature in the UW scheme is that TKE is diagnosed instead of being prognostic. The UW scheme improves the numerical stability and efficiency for the long time step needed in climate models. The calculation of PBL height in the UW is based on the Richardson number but with a more sophisticated approach (Grenier and Bretherton 2001). The critical Richardson number is 0.19 in the default version of the UW scheme in the WRF Model.
3. Descriptions of observational data and model configuration
a. Observational data
Observational data from two field campaigns are used to evaluate our numerical simulations. The first field campaign is the Variability of American Monsoon Systems project (VAMOS) Ocean–Cloud–Atmosphere–Land Study (VOCALS) Regional Experiment in October and November 2008. It is an international field experiment designed to better understand physics and chemical processes central to the climate system over the SEP (Wood et al. 2011), where the subtropical marine stratocumulus and shallow cumulus clouds are dominant. Some fundamental dynamical and physical processes related to cloud formation and variability are still not well understood and are difficult to simulate properly in most weather and climate models (Wyant et al. 2010). The second field campaign is the Land–Atmosphere Feedback Experiment (LAFE) over the SGP (Wulfmeyer et al. 2018) in August 2017. It deployed several state-of-the-art scanning lidar as well as other in situ and remote sensing instruments to provide three-dimensional measurements of numerous dynamic and thermodynamic quantities for evaluating and improving understanding of land–atmosphere interactions.
1) Data from VOCALS 2008 over the SEP
All data are obtained from https://archive.eol.ucar.edu/projects/vocals/rex.html. A detailed study of the stratocumulus clouds and boundary layer structure along 20°S during the VOCALS 2008 was described in Bretherton et al. (2010). The details of data collected and used in this study are presented below.
Dropsonde data were collected for the cross-section missions along 20°S from the west coast of South America to close to the buoy station at 85°W (Fig. 1a). In total, 47 high vertical resolution soundings were obtained from 31 October to 13 November 2008 on board the Facility for Airborne Atmospheric Measurements (FAAM) Bae-146 aircraft. The highest altitude for each sounding profile was around 7200 m above sea level. Loehrer et al. (1996) described data quality control for the dropsonde.
The radiosonde data over the ocean were provided by the NOAA R/V Ronald H. Brown (operated by NOAA) during the cruise between 18° and 22°S. The raw sounding data were interpolated to vertical levels with 10 m intervals. After quality control, there were 86 soundings during the time period of our model simulation (see next subsection). The 123 soundings at Iquique (marked in Fig. 1a) were also used in this study. All radiosonde soundings extended above 10 000 m above sea level. The dropsonde, radiosonde soundings and model outputs are interpolated into 25-m intervals in the vertical for our analyses.
(iii) Liquid water path (LWP)
The 10-min temporal resolution LWP data were collected on board the NOAA R/V Ronald H. Brown. They covered the time period 22 October to 2 November and 11 to 14 November 2008. The LWP data were calculated using radiances from a two-channel 20–30 GHz radiometer, following a physical retrieval updated from that of Zuidema et al. (2005). The other dataset was the daily collocated measurements of LWP data from the Cloud Profiling Radar (CPR) on board CloudSat, the Advanced Microwave Scanning radiometer-EOS on board the Aqua satellite, and the Special Sensor Microwave Imager (SSM/I) on board the Defense Meteorological Satellite Program satellites. This dataset has a resolution of 0.25° × 0.25° (Wentz 1997; Stephens et al. 2008; Brunke et al. 2010) and it is used to evaluate the spatial distribution of LWP in our model simulation.
(iv) Cloud-top height (CTH)
The CTH was computed using the raw data from the Leosphere ALS300 aerosol lidar on board the NERC Dornier-228 aircraft, which was operated by the Airborne Research and Survey Facility (ARSF) of the Natural Environment Research Council (NERC) in the United Kingdom (Wood et al. 2011).
2) Data from LAFE 2017 over the SGP
The data for LAFE 2017 were collected near Lamont, Oklahoma. (All data can be accessed from https://www.arm.gov/research/campaigns/sgp2017lafe.)
(i) The Atmospheric Emitted Radiance Interferometer (AERI)
AERI was designed by the University of Wisconsin Space Science and Engineering Center for the Department of Energy Atmospheric Radiation Measurement (ARM) program (Knuteson et al. 2004). It is a ground-based infrared spectrometer that measures downwelling infrared radiation at approximately 1 cm−1 resolution from 520 to 3000 cm−1. Using a retrieval algorithm named AERIoe, high temporal resolution (5 min used in this study) profiles of temperature and water vapor can be retrieved from AERI spectra in clear and cloudy scenes (Turner and Lohnert 2014). The information content in the spectra is restricted to the lowest 3 km of the thermodynamic profile or up to cloud base if one exists. The vertical interval gradually increases from 20 m near the surface to ~100 m at 1 km height and to ~230 m at 2.5 km height.
(ii) The flux tower measurements
There are three micrometeorological towers installed near Lamont, Oklahoma, by the Atmospheric Turbulence and Diffusion Division of NOAA. Tower 1 (36.611057°N, 97.482262°W) and Tower 3 (36.620662°N, 97.468153°W) were installed over mature soybean crop. Tower 2 (36.616448°N, 97.474176°W) was installed over a field that encompassed a field of grazed pasture. All three towers are located in the same grid cell of the WRF domain. Therefore, we averaged the measurements from three towers in order to evaluate the model simulation at that grid cell. The land-use type at the grid cell is croplands in the IGBP-modified MODIS 20-category land-use categories. The towers measured temperature, humidity, wind speed and direction, pressure, fluxes of momentum, sensible heat, latent heat, TKE, etc.
(iii) Doppler lidar TKE
The soundings were collected four times (0000, 0600, 1200, and 1800 UTC) a day during LAFE 2017. The sounding data provide the measurements of atmospheric temperature, humidity, wind speed, and wind direction with high vertical resolutions. We, therefore, interpolated the data into 25-m intervals in the vertical.
b. Model configuration and experimental design
The Advanced Research WRF (ARW), version 4.0 (Skamarock and Klemp 2008), is configured with a single domain at 3 km grid spacing. The model atmosphere is discretized with 51 vertical levels on the hybrid vertical coordinate with the model top at 50 hPa and 21 vertical levels below 2.5 km. The lowest model level is at 25 m above the surface. In the SEP experiment, there are 1200 grid points in the zonal direction and 800 grid points in the meridional direction for the domain over the SEP (Fig. 1a). In the SGP experiment, a smaller domain with 360 grids for each direction is used to cover part of the SGP (Fig. 1b).
The WRF single-moment 6-class cloud microphysics scheme (WSM6) with six prognostic cloud variables is employed for gridscale cloud microphysical processes (Hong and Lim 2006). The shortwave and longwave radiation fluxes are calculated by the RRTMG scheme (Iacono et al. 2008). The Noah LSM is used for the land surface processes (Chen and Dudhia 2001). Four PBL schemes are compared for each of the focused regions, namely, the scale-aware YSU scheme, scale-aware MYNN scheme, EEPS scheme, and UW scheme, as outlined in section 2. Because the scale-aware factor is not effective at a 3-km horizontal grid spacing (Shin and Hong 2013), we will simply use the names YSU and MYNN for these two schemes in the following discussion. All four PBL schemes share the same surface layer scheme (Jiménez et al. 2012) so that the simulated differences among the four schemes are not caused by the use of different surface layer calculations.
For the SEP region, the model initial and lateral boundary conditions for the atmosphere and land are obtained from the reanalysis data produced by the European Centre for Medium-Range Weather Forecasts (ECMWF) interim reanalysis (ERA-Interim) project (Dee et al. 2011). The sea surface temperature (SST) is given and updated daily using the 0.25° × 0.25° global analysis provided by National Oceanic and Atmospheric Administration (Reynolds and Chelton 2010). The diurnal SST variation is calculated based on the surface energy budget with the method of Zeng and Beljaars (2005). The model is initialized at 0000 UTC 24 October 2008 and integrated continuously until 0000 UTC 15 November 2008, which covers the main VOCALS field campaign period. The first-day model results are considered the spinup and thus are not included in the analysis discussed below.
For the SGP region, the model initial and lateral boundary conditions are provided by the National Centers for Environmental Prediction (NCEP) operational Global Forecast System (GFS) analysis, which has a horizontal grid spacing around 13 km. Due to the unusual wet season in August 2008, we select five consecutive days from 1200 UTC 27 August to 1200 UTC 1 September 2017. During this period, the dry convective boundary layer (CBL) dominated in day time. Initialized at 0000 UTC 27 August 2017, the model forecasts 36 h, and the first 12-h forecasts are regarded as model spinup and discarded in the analysis. Therefore, there are five forecasts for each PBL scheme. To better initialize Noah LSM, we conducted 3-month spinup runs for each of the three PBL schemes, namely the YSU, EEPS, and MYNN schemes. We then averaged the soil temperature and moisture from these runs to initialize the first forecast (e.g., initialized at 0000 UTC 27 August 2017). The soil temperature and moisture for the next forecast is initialized with previous forecast of the corresponding scheme.
a. The SEP
The subtropical marine boundary layer (MBL) is typically well mixed with relatively weak diurnal variations in terms of the height of the MBL (Wulfmeyer and Janjić 2005; Liu and Liang 2010). The diurnal variations of LWP and cloud cover are relatively large but the differences between the maximum and minimum are no more than a factor of 2 (Wood et al. 2002). Therefore, we will mainly focus on the time averages from 0000 UTC 25 October to 0000 UTC 15 November 2008 in our analyses.
Figure 2 shows the EEPS simulated vertical cross sections of TKE and energy dissipation rate (EDR), which is redefined as EDR = ε1/3 (Muñoz-Esparza et al. 2018). TKE and EDR are meridionally averaged between 18° and 22°N (the red box in Fig. 1a). TKE gradually increases from near the coast to west at 87°W (Fig. 2a) in response to the westward SST increase (not shown). EDR has the same pattern as that of TKE (Fig. 2b). Both TKE and EDR slowly decrease with height. The probability distribution functions (PDFs) of TKE and EDR at different heights are shown in Fig. 3. TKE is rarely smaller than 0.1 m2 s−2 in the lower part of the PBL below 355 m above sea level. Small values of TKE are more common in the upper part of the PBL. At the 678-m height above sea level, 8% of TKE values hit the minimum value (1 × 10−4 m2 s−2) that is set in the code in order to make the scheme less oscillatory. The change of PDFs with height for EDR is similar to that of TKE. The minimum allowed value for ε is 1 × 10−6 m2 s−3. Due to the lack of sufficient observations of TKE and ε over the SEP, we are not able to compare the simulation with observations. The magnitude of TKE is comparable to that found from large-eddy simulation (Bretherton et al. 1999) of the stratocumulus and shallow cumulus during the Atlantic stratocumulus transition experiment (Albrecht et al. 1995).
The vertical cross sections of cloud water mixing ratio and virtual potential temperature are shown in Fig. 4. The MBL is very well mixed as indicated by the contour of virtual potential temperature. The height of the simulated low-cloud layer appears to be largely controlled by the location of the inversion base (Fig. 4). The features simulated with the YSU and EEPS schemes, such as the height of the inversion base and the distribution of cloud water are very similar, in spite of the EEPS scheme simulating higher cloud water content than the YSU scheme. In contrast, the MYNN scheme has a higher inversion base but much less cloud water content (Fig. 4c). The UW scheme also produces a higher inversion base.
The diagnosed PBL height by each PBL scheme (Fig. 5) may not be necessarily consistent with the CTH (Figs. 4 and 6). Here CTH is defined as the height above which the simulated cloud liquid water content falls below 0.025 g m−3 (Wang et al. 2011). Note that CTH observations are mainly east of 75°W. The simulated CTHs by the MYNN and UW schemes are the closest to the observations with only about 50 m higher from the MYNN scheme and about 50 m lower from the UW scheme (Fig. 6). The YSU scheme simulated medium CTH of around 1000 m, while the medium CTH from the EEPS scheme is around 950 m, which is the lowest (Fig. 6). As a result, the diagnosed PBL height by all four PBL schemes are a few hundred meters lower than the simulated CTH. The difference between the simulated CTH and the diagnosed PBL height is smallest from the EEPS scheme.
The spatial distributions of LWP from observations and simulations with the YSU, EEPS, MYNN, and UW schemes are compared in Fig. 7. All schemes underestimate LWP, especially the MYNN scheme. Among the four schemes, the EEPS scheme simulates the highest LWP although it is still 20–30 g m−2 less than that observed. The simulated LWP is much closer to the observations if only counting the observations or simulated grid cells with cloud liquid water (Fig. 8). The simulated LWP by the MYNN scheme is closest to the observed, while the UW scheme consistently simulated smaller LWP. Both the YSU and EEPS schemes lack clouds with large LWP, while the MYNN scheme simulates clouds with much larger LWP, which agree better with the observations. Note that both cloud microphysical parameterization and aerosol concentrations might impact the simulated LWP. The impact, however, is unknown.
Figure 9 shows the mean vertical profiles of potential temperature, specific humidity, and wind speed for zone1, zone2, and Iquique (Fig. 1a). Regardless of the deeper boundary layer depth from the MYNN scheme in comparison to observations over zone1, where SST is warmer (not shown) and chances of decoupling of shallow cumulus and stratocumulus are higher, the simulated boundary layer depth by the MYNN scheme is slightly deeper than that in observations over zone2 and Iquique. The boundary layer depth simulated by the UW scheme is very close to observations in all three evaluated regions. The differences between the YSU and EEPS schemes are surprisingly small although both schemes have the boundary layer depth about 200–300 m shallower than that in observations. The underestimated boundary layer depth is consistent with the underestimated CTH for both schemes. The wind speed does not change significantly below the inversion layer. Despite the employment of nonlocal term in the YSU scheme for momentum, the differences of wind speed between the YSU and EEPS schemes are smaller over zone2 and Iquique.
To understand the dependence of the simulated boundary layer depth on the PBL scheme employed, we now examine the vertical profiles of potential temperature, water vapor mixing ratio, and meridional wind speed tendencies from all four schemes (e.g., tendencies due to the PBL parameterizations). Tendencies are calculated as an accumulated variable within 2-h time window at each location of the ship-based radiosonde or aircraft-based dropsonde. In the simulation with the YSU scheme, there is strong cooling near the top of the inversion layer, which is about 1.2 km over zone1 (Figs. 9a and 10a), 0.9 km over zone2 (Figs. 9d and 10d), and 0.8 km over Iquique (Figs. 9g and 10g). Accordingly, there is a strong warming below the cooling layer. We notice that this strong cooling/warming is mainly above the diagnosed PBL height (600–700 m over zone1 and 500–600 m over zone2). In the YSU scheme, the enhancement of vertical mixing in the boundary layer is achieved by the downward heat transport from the inversion layer. Additional enhancement of vertical mixing is achieved by recalculating and lowering the gradient Richardson number in clouds in free atmosphere. Due to the underestimation of PBL height, this enhancement actually occurs in the boundary layer clouds. In contrast, the cooling near the inversion layer is much weaker from the EEPS scheme. In the simulation with the MYNN scheme, the mass flux part of EDMF might be very active due to the continuous positive surface buoyancy flux over the ocean (not shown). The cooling layer is between 1.5 and 2 km over zone1 (Fig. 10a) and 1.0–1.6 km over zone2 (Fig. 10d), consistent with the vertical profiles of potential temperature shown in Fig. 8. The UW scheme has significantly larger potential temperature and moisture tendencies than any of the other PBL schemes studied. The mean entrainment cooling over Iquique is much weaker due to its diurnal variation over land. The entrainment process is typically associated with cooling (Figs. 10a,d) and moistening (Figs. 10b,e) in the previous inversion layer, and thus the deepening of the boundary layer. We only show the tendencies for meridional wind since it is the dominant wind direction for the SEP during VOCALS 2008 (Figs. 10c,f and i). The biggest differences are in the upper part of the PBL, where the MYNN scheme shows larger tendencies, especially over zone2 and Iquique (Figs. 10f,i).
The vertical mixing coefficient for momentum Km and heat/moisture Kh as well as the Prandtl number are shown in Fig. 11. Note that the vertical mixing only represents the eddy-diffusivity part in EDMF of the MYNN scheme. The height of vertical mixing gradually decreases from Zone1 to Iquique in response to the decreasing SST. The magnitude of the vertical mixing over zone1 and zone2 almost remains the same but is significantly reduced at Iquique for the YSU and EEPS schemes. The vertical mixing increases again above the diagnosed PBL height in the YSU scheme. Based on the three-dimensional distribution of the vertical mixing coefficients, we believe that the parameterized enhancement of vertical mixing in clouds predominantly contributes to the increase of the vertical mixing coefficient in the YSU scheme. The additional contribution of heat flux resulting from the entrainment at the top of the PBL, which is parameterized with vertical mixing coefficient K as well, might be ill-placed. Note that the differences in vertical profiles of potential temperature and water vapor mixing ratio are small between the YSU and EEPS schemes as shown in Fig. 9. Consistent with the very large tendencies, the UW scheme has much larger vertical mixing coefficients than the YSU and EEPS schemes. The Prandtl number from the EEPS scheme is around 0.8 in the PBL while it continuously increases from ~0.5 near the surface to 1.0 at the top of the PBL from the YSU scheme. For the MYNN scheme, the Prandtl number in the lower part of the PBL is around 0.75 but it quickly increases to 5–7 in the upper part of the PBL, indicating that the vertical mixing of momentum is more active corresponding to the larger tendencies there (Figs. 10f,i).
b. The SGP
The diurnal variations are significant over the SGP. The performance of the simulated convective boundary layer (CBL) and nocturnal boundary layer (NBL) over land needs to be evaluated separately. All observations and the results from simulations are illustrated at Lamont, Oklahoma.
The temporal and spatial characteristics of the EEPS simulated TKE and ε are analyzed first. Figure 12 presents the temporal variations of the observed TKE at 10-m height and the simulated TKE at the first model level near the surface from 1900 local time (LT) 27 August to 0700 LT 01 September 2017. The EEPS scheme simulates TKE very well compared with observations in terms of the time consistency and intensity. To better evaluate the vertical variation of TKE with time, we compare the EEPS simulated TKE with observations from Doppler lidar in Fig. 13. The EEPS simulation well captures the features of the diurnal variations of TKE. In particular, the height of large TKE is very close to that in observations (Figs. 13a,b). For example, both the observed and the EEPS simulated TKEs reach 2.2 km height on 28 October. Note that Fig. 12 shows underprediction of TKE within the atmospheric boundary layer (ABL) during most nighttime conditions, while Fig. 13 shows overprediction of TKE within the ABL during most nighttime conditions. This inconsistency might be caused by the uncertainties of observations during nighttime. The spatial and temporal pattern of EDR is similar to that of TKE (Fig. 13c) except that the largest EDRs are always near the surface.
The 5-day air temperature at 2-m height (T2m) is shown in Fig. 14a. The simulated T2m is close to observation during the day but 2–4 K warmer at night, compared to observation. The sensible heat flux (SHF, Fig. 14b) is well simulated by the YSU, EEPS, and UW schemes for both daytime and nighttime while the MYNN scheme slightly underestimates SHF during daytime. Overall, all four PBL schemes capture reasonably well the specific humidity at 2-m height (Q2m, Fig. 14c) and latent heat flux (LHF, Fig. 14d). The YSU, UW, and EEPS schemes slightly overestimate LHF during daytime in most of the days. The wind speed at 10-m height (WSP10m) simulated by all four PBL schemes is higher than that observed at night for most of the days, but the surface friction velocity (U*) simulated by all four schemes is slightly larger than that observed (Figs. 14e,f). Although the simulated WSP10m during daytime is close to that observed, U* is still higher in all simulations. The strong WSP10m and the associated large U* by the three PBL schemes (the YSU, EEPS, and MYNN) in a short time period are related to a storm cell in the vicinity of Lamont, Oklahoma, which is obviously not in observation. The YSU, UW, and EEPS schemes perform very similarly in predicting these near-surface variables while the MYNN scheme performs slightly differently although the same surface layer scheme is used in all simulations.
Many methods have been developed to calculate the PBL height from different observational data. Here, we adopted the method developed by Heffter (1980) to analyze potential temperature profiles for the existence of a critical elevated inversion, which is assumed to indicate the top of the PBL that meets the following two criteria: 1) Δθ/Δz ≥ 0.005 K m−1 and 2) θtop − θbase ≥ 2 K, where Δθ/Δz is the potential temperature lapse rate, θtop and θbase are the potential temperature at the top and bottom of the critical inversion layer. Marsik et al. (1995) reported that this method could overestimate the stable boundary layer (SBL) height in some cases due to the second criterion. To alleviate this issue, the Atmospheric Radiation Measurement Project (ARM, https://www.arm.gov/capabilities/vaps/pblht) slightly modified the Heffter method by iteratively relaxing θtop − θbase from 2 to 0.1 K to provide PBL height for both the SBL and CBL (https://escholarship.org/uc/item/6pp1d93m). We followed the algorithm used by ARM and calculated the PBL height from AERIoe algorithm retrieved potential temperature profiles, radiosonde observed potential temperature profiles, and WRF simulated potential temperature profiles (Fig. 15). Keep in mind that each PBL scheme has its own algorithm to calculate the PBL height. Note also that the difference in the vertical resolution between AERIoe, radiosonde observation and WRF output could contribute to the differences in the diagnosed PBL height even with the same algorithm. For the CBL, the EEPS scheme simulates the highest PBL height among the four schemes with their own diagnostic algorithm in each PBL scheme, closer to that in observations. For the SBL, both the EEPS and YSU schemes simulate higher PBL height than that in observations while the PBL height diagnosed in the MYNN scheme is close to that observed. The different methods used to calculate the PBL height for SBL in different schemes could account for the differences because the vertical profiles of potential temperature, specific humidity, and wind speed are very similar among all simulations (Fig. 16). Note that the CBL in the simulation with the MYNN scheme maintains 1–2 h longer than observed, and compared to simulations with other schemes (Fig. 15). However, this could be just due to the different methods used to diagnose the PBL height because the CBL also last much longer if Heffter’s method (1980) is used. Compared to that diagnosed in each scheme, the PBL height calculated with the method of Heffter (1980) is typically higher for CBL and lower for SBL.
Four radiosonde soundings provide daily observations during LAFE 2017. Figure 16 shows the 5-day mean potential temperature, specific humidity, and wind speed at 0700, 1300, 1900, and 0100 LT, respectively. The vertical profiles in all simulations generally have less discrepancy at 0700 LT. All schemes overestimate wind speed below 0.6 km height. The CBL at 1300 LT is well simulated by all schemes. The height of the CBL simulated by the MYNN scheme is the closest to that in observations and is about 200 m lower than that simulated by the YSU and EEPS schemes. The UW scheme even simulates a higher PBL height. The observed potential temperature profile at 1300 LT is weakly unstable in the lower part and weakly stable in the upper part of the CBL while the simulated potential temperature profiles show different characteristics of stability. Both the MYNN and UW schemes simulate a quite unstable CBL, the YSU scheme produces a nearly neutral profile, and the EEPS scheme simulates an unstable lower CBL and a slightly stable upper CBL, which is consistent with observations. The nonlocal mixing in both the YSU and EEPS schemes is relatively small, but can neutralize the gradient by cooling the lower part and warming the upper part of the PBL. The late afternoon CBL at 1900 LT is well simulated by all four schemes (Figs. 16g–i) although the simulated CBLs are less stable than that observed and there are ~0.3 K cold biases by the MYNN scheme, ~0.3 K warm biases by the YSU and EEPS schemes, and ~0.5 K warm bias by the UW scheme below 0.4 km height. The biases of water vapor and wind speed are small at 1900 LT. The potential temperature profile of the SBL at 0100 LT is again well captured by all four schemes. There are some relatively large biases for specific humidity and wind speed below 0.4 km height (Figs. 16k,l). Overall, both the YSU and EEPS schemes again perform very similarly.
The 2-h window centered at 0700 LT (the sunrise time is 0720 LT) includes part of the early buildup of the CBL as shown in Fig. 17a. The near-surface layer starts to warm up and the entrainment cools the base of the inversion layer. Due to the strong heating of the land surface, the CBL starts to grow quickly. In the early afternoon, the CBL has grown to 1.5 km in response to strong warming in the mixed layer and strong cooling in the lower part of the inversion layer. Note that the intensity of warming accompanied by drying in the mixed layer and cooling accompanied by moistening in the inversion layer is the strongest in the UW scheme and the weakest in the MYNN scheme (Figs. 17d,e). The differences between the YSU and EEPS schemes are relatively small. The vertical mixing is strong for meridional wind as well (Fig. 17f). In the late afternoon, the CBL is well mixed and the CBL depth reaches its maximum (Figs. 17g–i). The largest tendencies appear near the surface due to the rapid decrease of surface heat flux. The large tendencies in the surface layer persist overnight (Figs. 17j–l). The accurate representation of the heat and momentum exchanges at the surface is key to realistically simulate/predict the near-surface atmospheric air temperature, humidity, and wind speed, in particular at night.
Figure 18 shows the vertical profiles of the exchange coefficients Km and Kh and the corresponding Prandtl numbers simulated in all four schemes. The exchange coefficients for the stable boundary layer regime are small but very different among these schemes (Figs. 18a,b,j,k). The exchange coefficients are typically less than 1 m2 s−1 (Figs. 18j,k), with the smallest coefficient in the simulation with the MYNN scheme. The magnitude of K depends on how the parameterization deals with the turbulence above the PBL. The TKE equation basically includes both TKE production and dissipation terms in the whole atmospheric column, therefore, the vertical mixing is parameterized with the same equations below and above the boundary layer in the EEPS scheme. The MYNN, on the other hand, calculates the mixing length (exchange coefficient) differently below and above the PBL height (Olson et al. 2019). The YSU also treats the exchange coefficients differently below and above the PBL height (Hong et al. 2006). Therefore, the accurate calculation of the PBL height is very important for the YSU and MYNN schemes. The vertical profiles of the exchange coefficients in the unstable regime show vertical mixing of momentum above the base of the inversion layer and the vertical mixing of heat and moisture below in the MYNN scheme. This indeed leads to an extremely large Prandtl number above the base of the inversion layer in the MYNN scheme (Figs. 18d–i). The large Km in the MYNN scheme causes large wind tendency above the base of the inversion layer (Fig. 17f), making wind speed there larger (Fig. 16f). In the CBL, the Prandtl number in the YSU scheme almost stays at ~0.8 in the mixed layer for a while and it gradually increases from ~0.5 in the lower part to ~1.0 in the upper part of the PBL (Figs. 18f,i). In the SBL, the Prandtl number is larger in the YSU, EEPS, and MYNN schemes. The Prandtl number is 0.8 in the UW scheme.
5. Conclusions and discussion
In this study we have implemented a TKE and TKE dissipation rate ε based 1.5-order closure PBL parameterization, namely the E–ε scheme (or EEPS) into the WRF Model, and compared its performance with three other commonly used PBL parameterization schemes, the YSU, UW, and MYNN schemes. The simulations with the four PBL schemes are conducted over the SEP where stratocumulus and shallow cumulus are dominant and over the SGP where strong diurnal variation of the PBL is common. The simulations are evaluated/compared against various observations during two field campaigns: the VOCALS 2008 over the SEP and the LAFE 2017 over the SGP. The TKE fields simulated with the EEPS scheme are examined over both the SEP and SGP, especially over the SGP where the simulated TKE is compared with that observed. Due to the lack of TKE observations over the SEP, we have only checked the three-dimensional distribution and magnitude of the simulated TKE. They are in reasonable ranges. Nevertheless, the TKE simulated with the EEPS scheme is very close to that observed at Lamont, Oklahoma.
The performance of the EEPS scheme over the SEP is comparable to that of the YSU scheme, while the EEPS scheme outperforms the YSU scheme in terms of the simulated spatial distribution of LWP. The simulated LWP by the EEPS scheme is slightly better than that by the YSU if only the observations or simulated grid cells with cloud liquid water are considered. Both the MYNN and UW scheme well simulate CTH and PBL height, but the MYNN scheme significantly underestimates LWP over the SEP. The diagnosed PBL height in the EEPS scheme is more consistent with the CTH and PBL height. In contrast, the diagnosed PBL height in the MYNN and UW schemes is too low over the SEP, and that in the YSU scheme is close to the cloud base instead of the cloud top, which parameterizes the vertical mixing enhancement through the cloud depth. However, the strong vertical mixing above the diagnosed PBL height does not deepen the boundary layer, probably because the limited cloud regimes could not effectively modify the state of the whole atmosphere and the associated artificial warming or cooling could be compensated by other physical and/or dynamical processes. The underestimated CTH and PBL height by the YSU and EEPS schemes are caused by the insufficient vertical mixing over the SEP. On the other hand, the EDMF-based MYNN scheme has strong enough vertical mixing, especially in the warm SST region, where the vertical mixing is even too strong and leads to a too high PBL height and too small LWP.
All four PBL schemes perform more similarly over the SGP in terms of the near-surface temperature, specific humidity, wind, heat fluxes, and the diagnosed PBL height. The diagnosed PBL height in the EEPS scheme is the highest among the four schemes, which is the closest to that observed in the unstable regimes. The MYNN scheme simulates the CBL that lasts 1–2 h longer than that observed and those simulated by the YSU and UW schemes. The EEPS scheme simulates the slightly unstable lower part and weakly stable upper part of the CBL, which is close to observations. However, the MYNN and UW schemes simulate a completely unstable CBL and the YSU scheme simulates a nearly neutral BL in the early afternoon.
The performances of the four PBL schemes over the SEP seem to suggest that except for the UW scheme, all other schemes have difficulties in accurately simulating the cloud-topped marine boundary layer, however, the UW scheme also underestimates the LWP while the simulated LWP in the EEPS scheme is the largest LWP among the four schemes and is closer to observations. Although the four schemes perform similarly in many aspects over the SGP, the EEPS scheme captures the slightly unstable lower part and weakly stable upper part of the CBL, which is consistent with observations. Therefore, results presented in this study suggest that the EEPS scheme is an alternative option in mesoscale models. Finally, it is worthy of noting that the TKE dissipation rate (ε) from the EEPS scheme is an independent variable and is important for many applications, such as comparing the prognostic ε with the diagnostic ε in other PBL schemes (Muñoz-Esparza et al. 2018) and evaluating the effect of dissipative heating on tropical cyclones (Ming and Zhang 2018). The directly simulating ε is also an advantage of the use of the EEPS scheme in high-resolution application models. Note that because of the data availability and large uncertainty in observations, the simulated ε has not been compared with any observations in this study. We plan to use observational data to evaluate the simulated ε in different boundary layer regimes in a future study to further improve the performance of the EEPS scheme.
This study has been supported in part by the National Key R&D Program of China under Grant 2017YFC1501602 and the National Basic Research and Development Project (973 program) of China under Contract 2015CB452805 and in part by the Special Fund for Meteorological Scientific Research in the Public Interest (Grant GYHY201406006). The authors are grateful to three anonymous reviewers for their constrictive review comments that helped improve the manuscript. All data from the VOCALS field campaign over the SEP are obtained from https://archive.eol.ucar.edu/projects/vocals/rex.html. The data from LAFE 2017 over the SGP at Lamont, Oklahoma can be accessed from https://www.arm.gov/research/campaigns/sgp2017lafe. The ERA-Interim data and the GFS analysis data are available at https://rda.ucar.edu. The NOAA SST data are provided by https://www.ncdc.noaa.gov/oisst/data-access.