Comparison of GEOS-5 AGCM planetary boundary layer

Accurate models of planetary boundary layer (PBL) processes are important for forecasting weather and climate. The present study compares seven methods of calculating PBL depth in the GEOS-5 atmospheric general circulation model (AGCM) over land. These methods depend on the eddy diffusion coefficients, bulk and local Richardson numbers, and the turbulent kinetic energy. The computed PBL depths are aggregated to the Köppen climate classes, and some limited comparisons are made using radiosonde profiles. Most methods produce similar midday PBL depths, although in the warm, moist climate classes, the bulk Richardson number method gives midday results that are lower than those given by the eddy diffusion coefficient methods. Additional analysis revealed that methods sensitive to turbulence driven by radiative cooling produce greater PBL depths, this effect being most significant during the evening transition. Nocturnal PBLs based on Richardson number are generally shallower than eddy diffusion coefficient based estimates. The bulk Richardson number estimate is recommended as the PBL height to inform the choice of the turbulent length scale, based on the similarity to other methods during the day, and the improved nighttime behavior.


Introduction
The planetary boundary layer (PBL) depth is important for surface-atmosphere exchanges of heat, moisture, momentum, carbon, and pollutants. Several studies have attempted to understand the uncertainty associated with the use of different PBL depth definitions and found the estimated PBL depth to depend substantially on the method chosen. Vogelezang and Holtslag (1996) examined the PBL depth by defining it using both bulk and gradient Richardson numbers and found that the choice of Richardson number, the critical number chosen, and the inclusion of surface friction impacted the results. Seidel et al. (2010) tested seven different PBL depth definition methods on radiosonde profiles. Using a single dataset, the estimated PBL depth was found to differ by up to several hundred meters. The use of different methods in their study also produced different seasonal variations. They concluded that it is necessary to compare different PBL depth estimates from different sources using the same method. In a later study, Seidel et al. (2012) recommended a bulk Richardson number based definition.
In the present study, seven different methods to compute the PBL depth were incorporated into the Goddard Earth Observation System (GEOS-5) atmospheric general circulation model (AGCM) (Rienecker et al., 2008;Molod et al., 2012) and intercompared using a single climate simulation. The seven methods are based on vertical profiles of the eddy diffusion coefficient for heat (K h ), the bulk (Ri b ) and local (Ri) Richardson numbers, and the horizontal, shear-based component of the turbulent kinetic energy (TKE). In order to provide insight into implications on the regional and global climate scale, results were aggregated onto the Köppen-Geiger climate classes over land (Peel et al., 2007).
The purpose of this study is two-fold. First, it analyzes differences among the PBL depth definitions evaluated diagnostically within the GEOS-5 AGCM. Results of this comparison will be used to develop a better state-dependent estimate of the turbulent length scale, which must be specified in the current model's turbulence parameterization. A second purpose of this study is to evaluate the influence of different processes, such as turbulence generated by shear and radiative interactions with cloud, on the PBL depth. The following section provides a model description and a description of the PBL depth diagnostics used. The third section presents results of the comparison and the final section contains the conclusions.

2
Model and PBL diagnostics 2.1

GEOS-5 model description
The GEOS-5 AGCM is a comprehensive model with many uses, including atmosphere-only simulations, atmospheric data assimilation operational analyses and reanalyses, and seasonal forecasting when coupled to an ocean model (Rienecker et al., 2008;Molod et al., 2012). An earlier version was used for the Modern-Era Retrospective Analysis for Research and Applications (MERRA) (Rienecker et al., 2011). The latitude-longitude hydrodynamical core of the GEOS-5 AGCM uses the finite volume dynamical core of Lin (2004) and the cubed sphere version is based on Putman and Lin (2007). The GEOS-5 AGCM includes moist physics with prognostic clouds (Bacmeister et al., 2006). The convective scheme is a modified version of the Relaxed Arakawa-Schubert of Moorthi and Suarez (1992), the shortwave radiation scheme is that of Chou and Suarez (1999), and Chou et al. (2001) describe the longwave radiation scheme.
The Catchment Land Surface Model is used to determine fluxes at the land/atmosphere interface (Koster et al., 2000) and the surface layer is determined as in Helfand and Schubert (1995). The model uses 72 vertical layers that transition from terrain following near the surface to pure pressure levels above 180 hPa.
Since details of the turbulence parameterization in the current version of the GEOS-5 AGCM (Rienecker et al., 2008;Molod et al., 2012) are relevant to the analysis of results of the current study, it is described here. The turbulence parameterization is based on the Lock et al. (2000) scheme, acting together with the Richardson number based scheme of Louis et al. (1982). The Lock scheme represents non-local mixing in unstable layers, either coupled to or decoupled from the surface. The parameterization computes the characteristics of rising or descending parcels of air ("plumes"), initiated due to surface heating or to cloud top cooling of boundary layer clouds.
The GEOS-5 AGCM implementation includes moist heating in the calculation of buoyancy and a shear-dependent entrainment in the unstable surface parcel calculations. It is formulated using moist conserved variables, namely the liquid-frozen water potential temperature and the specific total water content, so that it can treat both dry and cloudy layers. The turbulent eddy diffusion coefficients are computed using a prescribed vertical structure, based on the height of the surface and radiative parcels or "plumes".
The Louis scheme is a first order, local scheme, and the eddy diffusion coefficients are computed using Richardson number based stability functions for stable and unstable layers. The Louis scheme unstable layer stability functions require the specification of a turbulent length scale, which is formulated using a Blackadar (1962) Molod et al. (2012) to compare well with a comprehensive set of observations.

PBL depth diagnostics
Seven different methods for determining the PBL depth are evaluated using the GEOS-5 AGCM based on several different output variables (Table 1). All methods diagnostically evaluate the same atmospheric profiles and all differences are related solely to the difference in definition of PBL depth.
The first method (Method 1) is based on the total eddy diffusion coefficient of heat (K h ) and estimates the PBL depth as the model level below that which K h falls below a threshold value of 2 m 2 s -1 . No vertical interpolation is used for this method and the estimated height is the model level edge. This method is the PBL definition used to determine the PBL depth in MERRA, and it is also used in the current GEOS-5 AGCM as part of the state-dependent estimate of the turbulent length scale. The evaluation of this method is one of the goals of the present study because any error in PBL depth shown to be associated with the use of this method may adversely influence the model's simulated climate.
Methods 2 and 3 use a variable K h threshold that depends on the atmospheric profile rather than a constant value. These methods use a threshold of 10% of the column maximum and linearly interpolate between levels to determine the PBL depth. Method 2 uses the total K h and Method 3 uses the surface buoyancy driven eddy diffusion coefficient (neglecting the contribution from the radiative plume). Method 3 therefore neglects the direct influence of clouds, and comparisons between this method and Method 2 isolate the role of the turbulence due to negative buoyancy at cloud top associated with cloud-topped boundary layers.
The PBL depth definition used by Seidel et al. (2012) is used as Method 4. They selected this method because of its applicability to radiosondes and model simulations and its suitability for convectively unstable and stable boundary layers. This method uses a bulk Richardson number (Ri b ) given by: where g is the gravitational acceleration, θ v is the virtual potential temperature, u and v are the horizontal wind components, and z is height above the ground. The virtual potential temperature, by definition, is based on water vapor, but not condensate. The subscript s denotes the surface. The surface winds are assumed to be zero. This bulk Richardson number is evaluated based on differences between the surface and successively higher levels, assuming that the surface layer is unstable, and the PBL top is identified as the level at which Ri b exceeds a critical value of 0.25. The PBL height is found by linearly interpolating between model levels.
Methods 5 and 6 use different versions of the bulk Richardson number, evaluated between two consecutive levels (rather than between the surface and the current height) that we term the "local" Richardson number. This local Richardson number (Ri) is calculated as: Here, z 1 and z 2 represent the heights of the model levels above and below the current level respectively, and θ v without a subscript is the average virtual potential temperature between heights z 1 and z 2 . The PBL top is found by assuming that the surface is unstable and linearly interpolating between the model levels where the critical value is crossed. We test two critical Richardson numbers to determine the sensitivity of the method to the critical value chosen. We use a scaling approximation of TKE to estimate the PBL depth in Method 7. The Lock scheme is not very sensitive to boundary layer shear so we chose a scaling based only on shear sources of TKE to isolate the shear contribution. The top of the PBL is taken to be the height at which the shear-based TKE falls below a threshold value of 10% of the column maximum, vertically interpolating between model levels. The horizontal TKE method should be more sensitive to the wind profile and seasonal changes to it than the other methods, and the daytime PBL heights based on this method should be expected to be lower than PBL height estimates based on static stability.

Climate classes
The computed PBL depths are aggregated by season onto the Köppen-Geiger climate classes ( Fig. 1). The Köppen-Geiger climate classes have been used to group rivers worldwide for comparisons of runoff characteristics (McMahon et al., 1992;Peel et al., 2004). Molod and Salmun (2002) successfully used this aggregation in their study investigating the implications of using different land surface modeling approaches. Their study aggregated results such as canopy temperature, soil moisture, and turbulent fluxes and they were able to use these results to make generalizations that extend to broad climate regions relevant for global models. Aggregation onto these climate classes is a way to characterize similar remote regions and apply findings globally. Peel et al. (2007) recently updated the Köppen-Geiger climate classification, taking advantage of advances in data availability and computing power. They did this by using monthly mean precipitation and temperature data from over 4000 stations (plus additional data from stations reporting only temperature or only precipitation) and interpolating between them using a twodimensional thin-plate spline with tension. The final map is generated on a 0.1°x0.1° grid. The highest station density is in the USA, southern Canada, northeast Brazil, Europe, India, Japan, and eastern Australia while the lowest station data densities are located in desert, polar, and some tropical regions. Peel et al. (2007) used the same classes as the original classification system, but with an updated distinction criterion between the temperate and cold climate classes. The classification consists of five main climate types: tropical (A), arid (B), temperate (C), cold (D), and polar (E) with further divisions based on seasonal variations in temperature and precipitation. Peel et al. (2007) provide a full description of the climate classifications including details on how the classification was determined. The broad climate types, defined over land, are relatively insensitive to temperature trends, including those from global climate change (Triantafyllou and Tsonis, 1994;Peel et al., 2007) and are intended to represent long term mean climate conditions and not year- to-year variability.

Results
This section describes the results of the comparison of the different PBL depth estimates aggregated to the Köppen climate classes. The first subsection (3.1) provides a quantitative description of the variability within climate classes, explains some of the reasons for this variability, and justifies the reliance on the climate class aggregated analysis The following subsections show the general PBL depth response to the different definitions, describe in detail the results from classes that deviate from this behavior, and examine in detail reasons for the difference between the PBL depths estimated using the K h and bulk Richardson number methods. The final subsection reports on the PBL height differences related to the cloudactivated Lock scheme's radiative plume.

Variability within climate classes
The Köppen-Geiger classification does not explicitly take into account some aspects of the climate system relevant to boundary layer processes such as intensity of precipitation, elevation, terrain, and overlying subsidence. The aggregation of PBL height onto climate classes is therefore useful for examining the behavior of the different estimates globally, but differences in behavior within climate classes are neglected by definition.  deviation being about 39% of the mean PBL depth. This climate class will be discussed in greater detail below. Figure 2b shows the summer mean diurnal cycle of PBL depth and standard deviation for the hot, arid, desert. This climate class also produces fairly uniform standard deviations through the diurnal cycle with a mean ratio of standard deviation to PBL depth of about 39%. Figure 2c shows the summer mean diurnal cycle for the hot summer, dry winter temperate climate class. In this class, the variability has a diurnal cycle in which the standard deviation is smallest at night and larger during the day. The mean standard deviation is about 31% of the PBL depth. However, during the dry winter, the variability is more uniform (not shown), similar to the dry climate class represented in Fig. 2b. Figure  Spatial maps in Fig. 3 show the relationship between PBL depth and surface temperature in the Sahara and Arabian deserts. Figure 3a shows the seasonal mean PBL depth estimated using Method 1 for JJA over the Sahara and Arabian desert part of the BWh climate class that was shown in Fig. 2b. In JJA, the PBLs over the coastal regions of the Saharan and Arabian deserts are more than a kilometer shallower than the PBLs found further inland. This behavior reflects the variability of the surface temperature within the BWh climate class. A spatial map of the JJA skin temperature (Fig. 3b) shows the same pattern as the PBL depth. A scatter diagram (not shown) of PBL heights and skin temperature revealed that >60% of PBL height variability is explained by skin temperature.
The second example of intra-class variability is illustrated in Fig. 4, which shows the relationship between PBL depth and 10-meter temperature for the tropical rainforest climate class (Af), shaded according to 10-meter relative humidity. In this climate class, and in the other tropical climate classes, there is a shift in the relationship between PBL depth and 10 m temperature near 302 K. This temperature is near the wilting point for broadleaf evergreen vegetation, the dominant vegetation type in the tropics. At temperatures above the wilting point, the vegetation experiences moisture stress, thus severely limiting transpiration and more of the net radiation at the surface is lost as sensible heat flux. Since sensible heat is much more efficient at growing the PBL than latent heat (Ek and Holtslag, 2004), the PBL depth increases rapidly with temperature in this drier regime. In the regime below the wilting point, transpiration increases with temperature and proceeds with little resistance, wetting the lower atmosphere. In this wetter regime, PBL depth decreases with temperature.
These different regimes and sensitivities of PBL depth to different variables must be kept in mind when examining climatological boundary layer depth. Although the Köppen-Geiger climate classes are useful for organizing land regions in order to make generalizations and simplify the analysis, they do not capture all the conditions relevant to boundary layer processes.
There will therefore be geographical differences within each climate class that will not be captured by this analysis.

General method behavior
When aggregated by climate class, the PBL depth definitions produce similar results for most classes and seasons. In general, both local Richardson number methods (Methods 5 and 6) estimate PBL depths that are lower than the other methods throughout the diurnal cycle. The bulk (Method 4) Richardson number method estimates shallower nocturnal PBLs than the K h methods (Methods 1, 2, and 3) and wintertime PBLs estimated by the TKE method (Method 7) are generally deeper than the other methods.
The focus of the discussion here is on illustrations of the significant differences based on the behavior of PBL depths from representative climate classes. Figure 5 shows the seasonal mean diurnal cycle for the cold climate class with warm summers and no dry season (Dfb; during summer 5a and winter 5c) and for the hot, arid desert class (BWh; during summer 5b and winter 5d m higher than the other methods (Fig. 5c) associated with the greater wintertime wind shear in the winter storm tracks within the Dfb climate class, and are 500 m higher in the winter (Fig. 5d) due to the wind shear aloft in the desert class. estimate PBL depths that are several hundred meters lower at midday than PBL depths using the other methods. This is the case for all the climate classes studied here. This method does not depend greatly on the critical value chosen as the differences between PBL depths estimated using a critical value of zero are only slightly lower than those estimated using a critical value of 0.2. The low PBL depths estimated by the local Richardson number methods make these methods impractical for AGCM-based PBL depth estimates.
Planetary boundary layers based on Richardson number methods (local and bulk) are lower at night than those based on K h or TKE for most classes in summer and winter. This has implications for estimating the shallow nocturnal boundary layer that has been shown to be relevant for constituent transport (e.g. Denning et al., 1995, Jacob et al., 1997, Lin and McElroy, 2010. For instance, over climate class BWh (Fig. 5b), the bulk Richardson number nocturnal PBL is well under 500 meters while the K h methods estimate a PBL depth between 1000 and climates where PBL depths are low for all methods (Fig. 5c).
The BWh climate class (Fig. 5b, 5d) contains radiosonde observations of the nocturnal boundary layer and during the evening transition from a convective to a stable boundary. The observations are from the American Southwest (one coastal station omitted), each represents a single radiosonde station, and do not sample the large desert regions in Africa and Australia, but they provide some insight into how well the model simulates the nocturnal PBL. The observed boundary layers are lower than those simulated by the model by approximately 100 to 300 m.
The radiosonde based estimates sample the PBL depth over the Dfb climate class (Fig. 5a and    5c) well because much of Eastern Europe and the northern United States belong to this climate class. Each observed point represents between 1 and 14 stations. Similar to the model behavior in the desert climate class, the model estimates higher nocturnal boundary layer depths than the radiosonde-based estimates during summer (mean difference of 210 m), and winter (mean difference of 155m). During the day, the mean difference between the model and radiosonde estimates during both seasons is more variable with differences ranging from approximately 10 m up to 150 m, but model estimates are generally lower.

Bulk Richardson vs. K h methods
The bulk Richardson number and K h methods generally give similar midday results, but under warm, wet conditions the estimated daily maximum PBL depth found using the bulk Richardson number method tends to be lower than the K h methods (Fig. 6). An example of this behavior is shown by examining the tropical rainforest climate class, but this occurs in the other tropical climate classes during their rainy seasons and for temperate climate classes when it is both warm and the climatological precipitation is high (not shown). This difference in estimated PBL depth means that the bulk Richardson number exceeds its critical value at a level below that which K h decreases below its threshold value. This implies either a virtual potential temperature inversion or a change in the wind speed within a layer of relatively high K h . Figure 7 shows the annual mean vertical profiles of total K h and K h from the Louis parameterization (7a) and the bulk Richardson number and virtual potential temperature perturbations (mean value of 307.9 K, 7b) from a typical location within the Amazonian rainforest. The horizontal dashed lines indicate the PBL depth found using the total K h (Method 1, Fig. 7a) and bulk Richardson number (Method 4, Fig. 7b). The bulk Richardson number method detects a stable layer below the level at which K h declines. This is due to the presence of a small inversion in the virtual potential temperature profile evident in Fig. 7c.
This behavior could occur under several different meteorological conditions. There could be a turbulent layer aloft that is not fully decoupled from the surface layer that is being detected by the K h methods, but not by the bulk Richardson number method. Since the Louis turbulence parameterization is dependent upon the local Richardson number (Ri), it contains some information about the vertical profile of temperature and shear. While this is a different form of the Richardson number than the one used in the bulk Richardson number method, the Louis scheme can provide information about what to expect from the bulk Richardson number method.
If the K h predicted by the Louis scheme alone (Fig. 7a) has its maximum in a shallow layer low to the ground before decreasing, it can be expected that the PBL depth found using the bulk Richardson number might also be low. If the Lock scheme is strongly active aloft due to entrainment or radiation, the K h methods will detect a deeper PBL.

Impact of the radiative plume
In order to examine the impact of radiative cooling at cloud top, the K h method using a threshold of 10% of the column maximum was compared diagnostically with (Method 2) and without (Method 3) the contribution from the radiative plume. The difference between these two methods is useful for understanding the influence of clouds on PBL depth in the GEOS-5 AGCM. Figure 8 shows the PBL depth difference between the two methods for JJA. At all locations, the PBL depth estimated using the radiative plume is at least as large as that without the radiative plume. The largest differences occur over land in the summer hemisphere and in the Tropics during the evening transition. This result also holds for December, January, and February (DJF) (not shown). The timing of the largest differences (evening) is due to the sensitivity of the radiative plume to cloud top. At night, the total K h decreases due to the lack of incoming solar radiation, but the diffusivity associated with the radiative plume decreases proportionally less since the cloud does not dissipate during the evening transition. The radiative plume eddy diffusion coefficient thus becomes proportionally more important at night and the PBL depth remains greater. The non-radiative method PBL heights are therefore lower at night, consistent with expectations.
Although this study focuses on the sensitivity of simulated PBL depths over land, there are persistent regions of relatively large radiative plume impact over the oceans as well, occurring around 30°N and 45°S. This is due in part to the behavior of the microphysics parameterization in the GEOS-5 AGCM and perhaps to the nature of low level clouds in these regions. The GEOS-5 AGCM uses an empirical estimate of cloud particle radii based on temperature, pressure, and wind. The large differences over oceans are located in regions where the boundary layer clouds contain condensate with small prescribed effective radii and are thus more radiatively active. Since the radiative plume is more active in these locations, PBL depths based on methods sensitive to its impact are greater than depths computed using methods that ignore it.

Conclusions
Although the PBL depth is important for AGCMs and its realism has implications for climate and weather prediction, observations are limited and no consensus on definition exists. The impact of longwave cooling from clouds on PBL depth was found to have its strongest effect over land during the evening transition. This was due to the persistence of cloud cover through the diurnal cycle. Additionally, regions of influence were found in the marine boundary layer related to the larger radiative impact in these regions.
The local Richardson number methods are relatively insensitive to the critical number used and estimate PBL depths several hundred meters lower than the other methods. These local Richardson number methods were therefore found to be inappropriate for use in an AGCM, probably due to the relatively coarse vertical resolution. The PBL depths found using the local and bulk Richardson number methods are generally lower at night than the PBL depth diagnosed using K h and TKE methods. We speculate that this result is due to the choice of K h threshold and that this threshold is more applicable to daytime convective boundary layers than to nocturnal PBLs.
The bulk Richardson number method (Method 4) provides the best match with radiosonde-based estimates using this method, as expected, and also provides the most credible diurnal cycle, due in great part to its capture of low nocturnal boundary layer heights. It is therefore the method recommended for use in estimating the AGCM turbulent length scale. Future work will include incorporating the PBL depth estimated using the various methods into the calculation of the turbulent length scale in the GEOS-5 AGCM. Through this length scale, the PBL depth is allowed to modify vertical mixing and tracer transport and the implications for air quality and carbon inversion studies will be analyzed. Horizontal TKE Uses the diagnosed horizontal turbulent kinetic energy and a threshold of 10% of the column maximum     summers and no dry season, during summer and winter, 5a and 5c) and BWh (hot, arid desert, during summer and winter, 5b and 5d) using 7 different methods for estimating the PBL depth.
The error bars represent two standard deviations for methods 1, 2, and 4. The green triangles indicate the observed PBL depth from the IGRA dataset and the green circles represent the modeled PBL depth (Method 4, green) at the observation locations.   The dashed line is the shortwave radiation zero contour line.