Spatial and temporal patterns of planetary boundary layer height during 1979–2018 over the Tibetan Plateau using ERA5

Planetary boundary layer height (PBLH) is a key meteorological parameter that is affected by and influences exchanges of energy, momentum, and mass between the surface and the atmosphere, but there are few comprehensive studies on the PBLH climatology over the Tibetan Plateau (TP). This study presents a multi‐decadal climatology of PBLH over the TP by using high‐resolution reanalysis data (ERA5) for the period 1979–2018. The ERA5 PBLH was first evaluated against radiosonde‐derived PBLH from four sites. For most of the observation periods and sites, ERA5 captured the diurnal cycles reasonably well, although the highest peaks in the radiosonde‐derived PBLH were often not captured. Next, long‐term means and trends, seasonality and diurnal cycles of the PBLHs were examined. The results show that average PBLHs in the central and southwestern TP are higher than those in the west and east. On an interannual time scale, the PBLH is positively correlated with the North Atlantic Oscillation in western and southeastern TP in winter, eastern TP in summer, and central TP in autumn, while a negative correlation is seen for northeastern TP in winter. The summertime PBLH in the central and southwestern TP is negatively correlated to the Indian summer monsoon. As for the seasonal cycle, the PBLH in the western TP reaches its peak in August, while the PBLH in the eastern TP exhibits an earlier summer peak in June. In the central TP, the PBLH is the highest in April and displays a minimum during summer, corresponding to a precipitation peak and a dip in surface sensible heat flux. The climatological diurnal cycles of PBLH are generally affected by the insolation and reach their peaks in the afternoon when a large range in PBLH variation is also found.


| INTRODUCTION
The planetary boundary layer (PBL) and its depth (or height; PBLH) is fundamental for understanding the structure of the lower troposphere (Seibert, 2000).The PBL is commonly defined as the portion of the troposphere which is directly influenced by the underlying surface of the Earth and responds rapidly (within approximately an hour) to surface forcings such as evapotranspiration, heat transfer, and frictional drag (Stull, 1988).Processes within the PBL determine the exchanges of momentum, heat, water, and chemical trace substances between the surface and the free atmosphere (Seidel et al., 2010).PBL characteristics, therefore, influence air quality as well as cloud formation and thereby the atmospheric heat and radiative budgets (Medeiros et al., 2005).
The PBLH varies over time and space, generally ranging between a few 100 m and a few 1,000 m (Stull, 1988).It is often used to determine the vertical extent of PBL turbulent mixing and the level at which exchange with the free atmosphere occurs (Seibert, 2000;Seidel et al., 2010).Accurate estimates of PBLH are therefore crucial in a wide range of atmospheric studies.However, few PBLH climatologies exist, and this is especially true for complex terrains (Collaud Coen et al., 2014).
Since the PBLH is affected by other meteorological variables, such as atmospheric stability and surface heat fluxes (Stull, 1988;Yang et al., 2004;Chen et al., 2016), it can be sensitive to climate change.For example, PBLH trends have been associated with changes in surface temperature and surface relative humidity (Zhang et al., 2013).Regional studies report positive PBLH trends during recent decades in, for example, Iran (Darand and Zandkarimi, 2019), Europe (Zhang et al., 2013) and Japan (Zhang and Li, 2019).A global study of atmospheric mixed layer height (i.e., the height of the convective boundary layer) found positive trends in Europe, Africa, Central US, East Australia, East Asia and Southeast Asia, and negative trends in coastal US, India and West Australia (Li et al., 2020).Guo et al. (2019), who studied PBLH derived from radiosondes in China from 1979 to 2016, found a shift around 2004 from positive to negative trends.In addition to the sign of the trend, the association with atmospheric and land-surface variables controlling the PBLH (soil moisture, temperature, relative humidity and stability) also changed between the periods.PBLH responses to current climate change thus appear to be highly variable.
The Tibetan Plateau (TP) is the highest and largest plateau in the world.The region, characterized by countless glaciers, lakes, valleys and mountains, is noted for its role in enhancing the Asian monsoons (Wu et al., 2012;Ge et al., 2017), supplying fresh water to much of Asia, and is currently undergoing rapid climate change (Yao et al., 2019).The TP and weather systems originating in its PBL exert profound influences on weather and climate in Asia (Gao et al., 1981;Tao and Ding, 1981;Wu et al., 2015).Moreover, due to large sensible heat flux in combination with weak stability of the free atmosphere (Yang et al., 2004;Chen et al., 2016), the TP PBL can grow very high-even up to the vicinity of the tropopause (i.e., fig. 2 in Chen et al. (2013)).For example, a study using radiosonde data from the Gerze station at the central TP found that on fair-weather winter days the PBLH can grow up to 5 km above ground level, and thereby reach stratospheric intrusions in the form of tropopause folds (Chen et al., 2013).Xu et al. (2001) studied the PBLH of Dangxiong County in the northwest of Lhasa based on the data obtained from the Second Tibetan Plateau meteorological experiment (TIPEX-II) and obtained a PBLH of 2,250 m.The PBLH of Naqu area and base cape area of Mount Everest was analysed by using radiosounding observation data (Li et al., 2004(Li et al., , 2006(Li et al., , 2011)).The results show that the PBLH in northern Tibet has different characteristics in dry and wet seasons, and the PBLH in the dry season is higher than that in the wet season.Due to the influence of glacial wind, the daily variation of PBLH in the Mt.Everest area is significant, and the maximum height can reach 3,888 m.Škerlak et al. (2019) used mesoscale model simulations to study the transport pathways of stratospheric air brought into the troposphere by folds.They found that over the TP, significant vertical transport in the free troposphere is not required for high concentrations of stratospheric tracer to reach the surface.Instead, as the PBL grows progressively during the day, stratospheric intrusion air is entrained at the PBL top (Chen et al., 2011) and transported downwards by turbulence, increasing the surface ozone concentration (Chen et al., 2013).At the same time, the elevated surface heating and subsequent deep moist convection allow for large detrainment of water vapour and other tropospheric constituents at the tropopause.Fu et al. (2006) found that compared to the monsoon area from which the water vapour originates, the temperature of the tropopause layer is higher over the TP, allowing for a greater moistening effect of the convection since less water is condensed out.Convection over the TP, therefore, provides the main pathway for water vapour from the monsoon and TP region to get into the lower stratosphere, indicating that global stratospheric water vapour concentrations could be highly sensitive to TP convection changes (Fu et al., 2006).Thus, the TP PBL can play an important role in cross-tropopause exchanges in both directions.Indeed, the TP is considered a hotspot region for deep stratosphere-troposphere exchanges, that is, exchanges in which stratospheric air reaches the PBL within 4 days or vice versa (Škerlak et al., 2014).
Improved understanding of the TP PBL is therefore crucial to gain a better understanding of its influences on Asian weather and climate as well as tropospherestratosphere mixing.Moreover, PBL characteristics such as the PBLH may vary greatly not only over time but also across the complex surface of the vast plateau, which calls for inquiries into its spatiotemporal variability.However, such investigations have been hindered by the fact that observations from which the PBLH can be estimated exist only for a few specific locations and limited periods of time.In addition, the fact that PBLH is not a variable that can be directly measured makes comparisons between datasets precarious, since different definitions and methods to estimate the PBLH may yield very different results.
Due to the difficulty of obtaining observations from complex-terrain regions, several studies applied Global Positioning System (GPS) Radio Occultation (RO) from COSMIC (Constellation Observing System for Meteorology Ionosphere and Climate) to study PBLH in several regions (e.g., Ao et al., 2012;Basha et al., 2019).It is well known that the RO data with high vertical resolution and low bias can retrieve the conditions of the atmospheric upper air, where the traditional observation is insufficient.It has been shown that RO-retrieved PBLH agrees with the heights estimated from radiosonde measurements over the TP (Zhou et al., 2018;Xiang et al., 2021).However, there is a difficulty in retrieving nighttime PBLH due to lower-resolution measuring close to the surface (Zhou et al., 2018).Moreover, due to the relatively recent implementation, RO cannot yet provide a long time series.
Reanalysis, which combines observations and numerical modelling to provide a physically consistent gridded dataset constrained by the observations, can deliver meteorological variables which can be used to estimate PBLH (von Engeln and Teixeira, 2013).As such, it offers a unique opportunity to investigate the multidecadal PBLH variability across the entire plateau.For this purpose, we use data from the ERA5 reanalysis, which is produced by the European Centre for Medium-Range Weather Forecasts (ECMWF) and replaces the widely used ERA-Interim (Hersbach et al., 2020).While the spatial resolution is not as high as regional datasets like, for example, the High Asia Reanalysis (HAR, Maussion et al., 2014), it does cover a longer period, allowing for assessments of multi-decadal variability.Due to ERA5's relatively recent release, there are not yet many studies evaluating its performance specifically for the TP region.Its predecessor ERA-Interim, however, has been used and assessed in a number of studies focusing on the TP.For example, Bao and Zhang (2013) evaluated data from ERA-Interim and other reanalyses against vertical profiles of horizontal wind, temperature and relative humidity from radiosondes over the TP.They found that for mean temperature and horizontal winds, the reanalysis data are consistent with the soundings, and that ERA-Interim is among the two datasets with the smallest average RMSE and bias.Bao and Zhang (2013) also compared profiles from radiosonde campaigns with data from ERA-Interim and three other reanalyses over the TP.Although they note that caution should be taken in trend assessments due to the existence of non-negligible time-varying biases in all four datasets, they found that overall, the reanalyses reproduce the profiles reasonably well and that ERA-Interim has the smallest temperature bias and RMSE at most vertical levels.
The PBLH from ERA-Interim has been compared with several modelling and observation datasets, and used to study PBL conditions (e.g., Patil et al., 2013;Chen et al., 2016).With respect to ERA-Interim, ERA5 benefits from higher spatial and temporal resolution as well as a decade of refinements to model processes and data assimilation methods and therefore shows considerable improvements in several regards (Hersbach et al., 2020).For example, Martens et al. (2020) found improvements in the surface energy partitioning and near-surface meteorology in a global comparison with in-situ data.Sun et al. (2021) analysed the local land-atmosphere coupling at the TP and found that ERA5 represents the investigated surface variables relatively well.They also found that ERA5 PBLH (in June, August and November 2014) generally shows a satisfying agreement with observation-derived PBLH.Recently, ERA5 was used to create a climatology for air-surface temperature differences over the TP (Wang et al., 2020a).
This study aims at developing a climatology of PBLH over the TP using ERA5, which allows us to reveal spatial and temporal changes in the PBLH in the last decades.Given that the existing evaluations of ERA5 PBLH over the TP are scarce, we first use radiosonde-derived PBLH from several Intense Observation Periods (IOPs) to evaluate ERA5 at four different stations.Thereafter, TP-average annual and seasonal PBLHs and the associated long-term trends, spatial distributions of the seasonal and annual mean PBLHs and their associated trends, as well as seasonal cycles of monthly mean PBLH and diurnal cycles of hourly PBLH variations for three selected regions are determined.

| ERA5, MERRA-2 and CFSR reanalyses
For the climatology presented here, we use ERA5 monthly and hourly boundary layer height during the period 1979-2018, which is readily available as a diagnostic variable in ERA5 and was downloaded from the Climate Data Store (https://cds.climate.copernicus.eu/#!/home),where post-processed ERA5 data on a regular longitudelatitude grid with a horizontal resolution of 0.25 × 0.25 are stored.The PBLH was computed in ECMWF's Integrated Forecast System (IFS) using a bulk Richardson number (Ri b ) method where PBLH is defined as the level where Ri b = 0.25.This method is based on an algorithm proposed by Vogelezang and Holtslag (1996) and deemed the most suitable for convective as well as stable boundary layers by Seidel et al. (2012), who reviewed different methods for estimating PBLH (ECMWF, 2016).To study the surface conditions, we used ERA5 monthly averaged data from 1979 to 2018, including total precipitation, surface net solar radiation as well as surface sensible and latent heat fluxes.
The annual trends and seasonal cycles of ERA5 PBLH were compared to two other reanalysis datasets.The monthly PBLHs from the Modern-Era Retrospective analysis for Research and Applications, Version 2 (MERRA-2) are provided by NASA Global Modeling and Assimilation Office (GMAO) (Gelaro et al., 2017) and can be downloaded from the NASA Goddard Earth Sciences (GES) Data and Information Services Center (DISC) (https://disc.gsfc.nasa.gov/datasets?project=MERRA-2).The PBLHs used in the study are at a resolution of 0.5 × 0.625 from 1980 to 2018.Note that the PBLH calculation in MERRA-2 applied the total eddy diffusion coefficient of heat with a threshold value of 2 m 2 s −1 (McGrath-Spangler and Molod, 2014), which is different from the ERA5 PBLH.
The National Centers for Environmental Prediction (NCEP) Climate Forecast System Reanalysis (CFSR) provides forecast PBLHs at a resolution of 0.312 × 0.312 (Saha et al., 2010) and can be downloaded from the National Center for Atmospheric Research (NCAR) Research Data Archive (RDA) hosted webpage (https:// rda.ucar.edu/datasets/ds093.2/).We used PBLHs from 1-hr forecasts in the regular monthly mean (4 per day) products during 1979-2010, which is calculated according to the bulk Richardson number method with Ri b = 0.25 (Hong and Pan, 1996;Saha et al., 2010).

| Radiosonde data for evaluation
PBLH estimated from radiosonde data collected during several Intensive Observation Periods (IOPs) taking place in 2008, 2013, 2014, and 2019(Chen et al., 2013;;Lai et al., 2021b) at the four TP stations Gerze, Ali (Ngari Station for Desert Environment Observation and Research), Qoms (Qomolangma Station for Atmospheric and Environmental Observation and Research), and Shiquanhe (Figure 1) were compared to ERA5 PBLH at the nearest grid point to investigate to what extent the reanalysis data capture diurnal cycle and day-to-day variability at these locations.The radiosonde measurements were generally conducted three times a day, approximately at 08, 14, and 20 (BST).For Gerze and the two IOPs conducted at Shiquanhe in 2019, additional measurements were made at 02 BST, while data was only conducted twice a day for most of the IOP conducted at Shiquanhe in 2013.From the radiosonde data, PBLH was calculated using a parcel method (Holzworth, 1964;Seibert, 2000;Hennemuth and Lammert, 2006) and a bulk Richardson method.With the parcel method, the PBLH is defined as the height at which the actual potential temperature profile intersects the dry-adiabatic profile, starting from the near-surface temperature.The bulk Richardson approach, based on Vogelezang and Holtslag (1996), uses the bulk Richardson number (Ri b ), which is a measure of stability, to estimate the PBLH (Chen et al., 2016).For the geometric height of each vertical level in the sounding, z where g is the globally averaged gravitational acceleration (9.80665 m s −2 ), θ is the potential temperature at each z and the surface (z 0 ), and u and v are the zonal and meridional wind components at each z.The PBLH is defined as the level where Ri b is the closest to 0.25.The same computing code for PBLH calculation from (Chen et al., 2016) is used in this study.The locations of all the radiosonde observation sites are shown in Figure 1.One spurious outlier was found for the second IOP at Shiquanhe, where the PBLH determined by the bulk Richardson approach was over 18 km above ground level.Since this is judged not to be a realistic value, it was removed prior to any calculations.

| Large-scale atmospheric circulation indices
Climate indices for North Atlantic Oscillation (NAO), Indian summer monsoon (ISM) and East Asian summer monsoon (EASM) were used to investigate the interannual correlation between the PBLH and these largescale climate patterns during 1979-2015.These patterns are known to be important for climate in the study region (Yao et al., 2019;Lai et al., 2021a).The NAO describes the variations of sea-level pressure differences between the Iceland low and the Azores high.Positive values of the NAO index are generally correlated to anonymously strong westerlies over the western Europe and the TP (Xin et al., 2010).Moreover, the winter NAO has been positively related to increased snowfall, snow cover and snow duration at the TP (You et al., 2011;Liu et al., 2018), while the summer NAO has been related to negative (positive) precipitation anomalies in the southern (northern) parts of eastern TP (Liu and Yin, 2001), which could influence the PBLH.We applied the NAO index that was calculated based on normalized sea-level pressure and can be found in the webpage of the UCAR/ NCAR Climate Data Guide (https://climatedataguide. ucar.edu/climate-data/hurrell-north-atlantic-oscillationnao-index-station-based).The variability in TP summer precipitation is related to the ISM and EASMI (Lai et al., 2021a), and precipitation, in turn, may alter the PBLH.The ISM index applied in this study was developed by Wang et al. (2001) who used the differences of zonal wind at 850 hPa between the northern India/TP and the Arabian Sea.Although the ISM index is available for June-September, we only use June-August in order to be consistent throughout our seasonal analyses.The index can be accessed from the Asia-Pacific Data-Research Center of the International Pacific Research Center at the University of Hawai'i at M anoa (http://apdrc.soest.hawaii.edu/projects/monsoon/definition.html).The EASM index for summertime (June-August) defines the EASM strength based on wind fields at 850 hPa within the East Asian monsoon domain (Li and Zeng, 2002).The index is available at http://ljp.gcess.cn/dct/page/65577.

| Analysis
We define the TP as the region which lies within 70 -105 E and 27 -40 N and has an altitude of >3,000 m above sea level, following, for example, Kukulies et al. (2019) and limit the study to the 1979-2018 period.For seasonal analyses, the ERA5 data are split into Northern Hemisphere standard seasons (December-February, March-May, June-August, September-November). Trends in PBLH are calculated with linear regression and tested for significance at the 95% confidence level.For any averaging across space, the grid cell areas are used as weights to account for the changes in grid cell size with latitude.
To assess the agreement between the datasets, relative Root Mean Square Error (rRMSE) and linear correlation (Pearsons' R) are calculated between ERA5 and radiosonde-derived PBLH.The rRMSE is obtained by dividing the Root Mean Square Error with the interquartile range (the difference between 75th and 25th percentiles) of the radiosonde-derived PBLH for a particular station and observation period.Since the diurnal cycle of PBLH most likely dominates the statistical agreement between the datasets, correlations and rRMSEs are also computed for the datasets after the diurnal cycle has been removed.The diurnal cycle is removed by subtracting the mean PBLH for a particular hour h from each observation conducted at h.In similarity to the trends, the correlations are considered statistically significant if p < 0.05, that is, the confidence level is 95%.

| Evaluation of planetary boundary layer height in ERA5
Radiosonde-derived PBLHs at the TP stations Gerze, Ali, Qoms, and Shiquanhe (Figure 1) are plotted along with ERA5 for the corresponding time (nearest hour) and place (nearest grid cell) in Figure 2. The observation dates for each location as well as the rRMSE and correlation between ERA5 and radiosonde PBLH are listed in Table 1.
In cases where radiosonde PBLH from both the parcel method (PBLH p ) and the bulk Richardson method (PBLH r ) are available, they generally agree well with each other, although there are cases of large differences between them, for example at Ali in November (Figure 2e).The ERA5 PBLH tends to agree slightly better with PBLH p than PBLH r .
ERA5 seems to capture the general pattern of the daily cycles relatively well.However, the highest PBLHs are often missed, both in the sense that the daily maximum tends to be lower in ERA5 and that the timing, as inferred from the available observations, appears to differ.For example, at the western TP station Ali (Figure 2d), the ERA5 PBLH often starts its daily decline a few hours before the peak is reached in the radiosonde data.
A recent study of diurnal cycles of summer precipitation at the TP also found that ERA5 tends to peak too early in the day (Ou et al., 2020), which seems to indicate that some processes related to diurnal variability are still not well represented in ERA5.Partly the mismatch can be due to the different scales associated with the in situ observations and the ERA5 grid spacing.Downscaling simulations may be helpful to understand the details and interconnections between relevant processes and improve their representation (Ou et al., 2020;Zhou et al., 2021).
The correlations between radiosonde-derived PBLH and ERA5 PBLH are mainly strong and statistically significant (Table 1).Some of this correlation arises in response to the diurnal cycle, but the significance remains when the correlations are performed after removing the diurnal cycle (correlation coefficients in brackets in Table 1).The highest correlations between PBLH and ERA5 are seen for Shiquanhe in IOP 2 when the correlation for PBLH p is 0.95 and the correlation for PBLH r is 0.94.The lowest correlation coefficients are seen for the two latter IOPS (August and November) at the southern TP location Qoms, when neither PBLH p nor PBLH r is significantly correlated with the ERA5 PBLH.Since topographic variations influence the PBL flows (Lai et al., 2021b), it is possible that the poor agreement with ERA5 at Qoms is related to the complex terrain around this station.Qoms is situated in a relatively narrow (about 1 km) valley approximately 40 km northwest of Mt.Everest (Han et al., 2015), where there is large variability in elevation over distances much too small to be resolved by ERA5.
Moreover, the area encompassed by the ERA5 grid cell which is centred on the coordinates closest to the Qoms station likely includes part of a glaciated region south of the station, which may affect the grid cell average.For such complex topography with highly heterogeneous surface conditions, a higher-resolution dataset would arguably be better suited to represent the PBLH at specific locations.However, the observations from the first IOP at Qoms (June) agree very well with the ERA5 PBLH, so the agreement is hardly determined by topography alone, but perhaps rather the interplay of topography and the variations in weather and surface conditions.As for the temporal variability in agreement at the other stations, the PBLH at Shiquanhe displays greater agreement during the first two IOPs (summer and spring) than during the third (late summer).In contrast, Gerze displays similar correlation coefficients (0.82-0.84) as well as rRMSE (0.35-0.37) in all three IOPs (which represents winter, monsoon onset and monsoon season, respectively (Chen et al., 2013)).However, its ability to capture the Note: Observations were conducted 2-4 times a day at four different radiosonde stations at the Tibetan Plateau.The columns display the name of each radiosonde station, the dates during which the observations were conducted, the relative root mean square error (rRMSE; the root mean square error between the datasets divided by the interquartile range of the radiosonde-derived PBLH) and coefficients from a linear correlation (Pearsons' R) between the datasets.The rRMSE and correlations were computed between ERA5 PBLH and radiosonde-derived PBLH calculated with a parcel method (rRMSE p and R p ) as well as a bulk Richardson method (rRMSE rib and R rib ).The rRMSE and correlation coefficients calculated after removing the diurnal cycle by subtracting the mean of all observations for a particular hour from each observation conducted at that hour are shown in brackets (not included for the first observation period at Ali and the third observation period at Qoms, where some observations were conducted at irregular times).Stars denote correlations that were significant at the 95% confidence level.One outlier, a bulk Richardson-derived PBLH of over 18 km above ground level at the Shiquanhe station, was removed prior to any calculations.
daily peak, as inferred from Figure 2a-c, seem to be considerably better during the first and last IOPs than during IOP 2. Despite often missing the highest peaks at other stations, it thus appears that ERA5 is able to represent the extremely high peaks that are sometimes reached at Gerze in winter, as well as the more moderate peaks in IOP 3. At Ali, PBLH r agrees relatively well with the ERA5 PBLH in both IOPs (June and November), while PBLH p shows a very poor agreement in the second IOP, with a non-significant correlation of 0.37 and RMSE p exceeding 13 m.As seen in Figure 2e, the PBLH p seems unable to capture the daily cycle as derived from PBLH r and ERA5 for this particular IOP and station.Thus, the bad agreement for PBLH p in this particular case may indicate the general difficulties in deriving and defining PBLH rather than the inability of ERA5 to reproduce it.
As noted above, deriving PBLH is not straightforward and differences between the datasets are likely to arise simply because of the differences in calculation methods and spatial coverage (the radiosondes are released at a specific point while ERA5 is averaged over the entire grid cell).Keeping this in mind, ERA5 reproduces the daily cycles reasonably well, except at Qoms, where the bad agreement is likely related to the complex topography and surface cover surrounding this station.However, ERA5 generally seems to be unable to capture the full height of the daily PBLH peak.

| Climatology and trends of planetary boundary layer height
When spatially averaged across the entire TP, the seasonal PBLH trends from ERA5 over the 1979-2018 period are positive in all seasons except summer (Figure 3).The strongest trends are seen for the winter (+5.2 m decade −1 ) and summer (−2.6 m decade −1 ), but statistical significance is lacking for all seasons as well as for the annual trend of 1.1 m decade −1 .The annual trends from ERA5, MERRA-2 and CFSR over the period covered by all three datasets (1980-2010) are −2.35,−6.87 and −13.1 m decade −1 , respectively, however, none of the trends are statistically significant at the 95% level.The correlations of annually and spatially averaged detrended PBLHs from ERA5 and MERRA-2/CFSR are 0.46 (statistically significant) and 0.08 (not significant), respectively.The results indicate that the spatially averaged 1980-2010 PBLH trends over the TP are not pronounced and that the interannual variability is inconsistent among reanalyses.
All three datasets show their seasonal peak in spring and early summer, with the highest median in May and the lowest median in December (Figure 4).In MERRA-2, there is a distinct local maximum in October which is not present in the other datasets, although the large spread in ERA5 PBLH in this month indicates that some of the years exhibit this small October peak in ERA5 too.The variation among the individual years is largest during spring and winter, and overall considerably greater in CFSR than in ERA5 and MERRA-2.The difference in the absolute heights is very large, with the median MERRA-2 PBLH being at least twice as high as the median ERA5 PBLH every month.The differences of heights between MERRA-2 and ERA5 could be possibly due to the differing definitions of the PBLH and the overestimation of snow cover in ERA5 (Orsolini et al., 2019).
In summary, while the high correlations between ERA5 and radiosonde-derived PBLH indicate a F I G U R E 3 ERA5 planetary boundary layer height (PBLH) time series and trends at the Tibetan Plateau (TP) during 1979-2018.Shown are the TP-average annual mean (thick line with round markers) as well as seasonal means for December-February (dashed line), March-May (solid line without markers), June-August (dotted line) and September-November (dash-dot line).The legend gives the size of the trends (m decade −1 ), of which none are significant at the 95% confidence level [Colour figure can be viewed at wileyonlinelibrary.com] reasonable performance of ERA5 for most of the investigated locations and periods, the comparisons among ERA5, MERRA-2 and CFSR point to large disagreement in the interannual variability and absolute heights of the reanalysis datasets.The disagreement suggests that variations on an interannual time scale from either of the three datasets should be interpreted with care since longer-term observation data for evaluation is lacking.For the long-term climatology, we choose to stick to one of them (ERA5) in this study.More studies are needed to determine the causes of the differences.
To investigate the spatial characteristics of the ERA5 PBLH, the 1979-2018 medians of the seasonal and annual means as well as the linear trend over the period are shown in Figure 5.In all seasons, the PBLH is greater in the interior of the central plateau and lower towards its rims as well as the western and eastern parts.At least for winter and spring, these spatial characteristics of the PBLH are similar to the pattern of surface-air temperature differences which are in turn dominated by downward radiation (Wang et al., 2020a).While observations from the westernmost TP are lacking, the tendency for decreasing PBLHs from the central-western plateau to the eastern has been seen in studies of PBLH from RO data (Zhou et al., 2018) and soundings (Che and Zhao, 2021).The spatial differences across the plateau are most pronounced during winter (Figure 5a) and spring (Figure 5b), which is also when the highest PBLHs (well over 1,500 m) are found.Previous studies have found that the TP's PBLH is higher during the dry season and lower during the wet season (e.g., Li et al., 2011).This may explain why the westernmost TP, which receives most of its precipitation during winter, displays its lowest PBLH during this season, while the central TP, where precipitation peaks in late summer (Lai et al., 2021a), exhibits its lowest PBLH in summer.The eastern TP, in spite of receiving most of its precipitation in early summer (Lai et al., 2021a), displays the lowest PBLH during winter, and may thus be more influenced by the seasonal variations in insolation and temperature.
The trends differ in size and direction between regions and seasons, with positive winter trends of up to 100 m decade −1 in the central, southern, and southeastern TP (Figure 5f) and negative spring and summer  5g,h).However, statistical significance is lacking, except in a few small parts of the regions with positive trends, for example, in southeastern TP in winter and autumn (Figure 5f,i).Although determining the causes for the trends are beyond the scope of this study, it can be noted that the southern parts of the region having increasing winter PBLH correspond to areas of increasing surfaceair temperature differences in ERA5 (Wang et al., 2020a), which likely drives increasing heat flux and PBL growth.As for the decreasing PBLH in central TP in summer, the wet season (July-September) precipitation has a positive trend in this region (Lai et al., 2021a).Given that the PBLH in this region, at least on a seasonal scale, is higher when dry conditions prevail (Figure 5a,b), it is plausible that the negative PBLH trend is related to the increased precipitation.
The seasonal PBLHs were correlated to the NAO and Monsoon indices to explore their potential linkages.The wintertime (December to February) PBLH is significantly (95% confidence level) positively correlated to the variations of the NAO in the western, southwestern and southeastern TP, and significantly negatively correlated in the northeastern TP (Figure 6a).For the summertime (June-August), a positive correlation is seen in the eastern TP (Figure 6c), and during the autumn (September-November), a positive correlation is present over much of the central plateau (Figure 6d).While investigating the mechanisms linking the NAO and Monsoon indices to the PBLH is beyond the scope of this study, previous studies have shown that the NAO is correlated to changes in summer precipitation (Liu and Yin, 2001;Wang et al., 2018) and snow cover (Immerzeel and Bierkens, 2010) over the TP.During the summer, weather, especially precipitation in the eastern TP is also affected by the Asian monsoons.The PBLH in the central and southwestern TP is significantly negatively associated with the variations in ISM (Figure 6f), in line with the seasonal variability discussed above with high (low) PBLH in this region coinciding with dry (wet) conditions (Figure 5a-d).Note that the summertime PBLH is not significantly correlated to the EASM index in any region (Figure 6e).
Based on the spatial distribution of the median PBLH over the entire period, three regions within which the PBLH is relatively homogenous were chosen to represent western, central, and eastern TP PBLH (red boxes in Figure 5e).The regions were selected to contain spatially coherent median PBLH, aiming to capture the cores of distinct PBLH regimes, and to be of similar size.The selected regions also represent three different precipitation regimes (Lai et al., 2021a) and are reasonably homogenous in terms of climatic conditions (Lu and Liu, 2010;Gao and Liu, 2013).In the following section, the seasonal cycle of the PBLH is discussed and compared to the seasonal cycles of the ERA5 variables total precipitation (i.e., the sum of large-scale and convective precipitation), surface net solar radiation, surface sensible heat flux, and surface latent heat flux.The surface net solar radiation represents the amount of solar radiation that reaches a horizontal plane at the surface of the Earth, minus the amount reflected by the surface and is given here as positive downwards.The turbulent heat fluxes are given as positive upwards.
In the western region (WTP), the PBLH exhibits a clear summertime peak with the median (of all the monthly mean values in the time span) reaching around 550 m in August (Figure 7a).The summertime peak suggests that PBLH at this location is strongly affected by the seasonal variation in surface turbulent heat fluxes, which are also highest during summer in this region, with the sensible heat flux peaking in August (Figure 8b).During winter and spring, the PBLH remains low, with values ranging from around 35 m in January and December to just over 100 m in May.During the part of the year when the PBLH is low, the statistical dispersion is also very small, indicating a small interannual variability.The persistently low PBLH during winter and spring could be affected by the large extent of snow cover (Wang et al., 2020b), which is known to be significantly overestimated in ERA5 (Orsolini et al., 2019).
In the central region (CTP), the yearly peak is reached in spring with median values of around 1,500-1,700 m (Figure 7b).The springtime maximum is followed by a summertime minimum with median values under 700 m in August before a second, smaller peak is reached in autumn.The seasonal pattern of PBLH resembles that of the surface sensible heat flux, which also displays high spring values and a minimum in August.In contrast, the precipitation as well as the surface latent heat flux show relatively low values during spring, followed by pronounced summertime maxima which coincide with a small depression in the surface solar radiation in July and August (Figure 8c,d).It thus appears that during spring, dry conditions and relatively large surface solar radiation give rise to large sensible heat flux promoting PBL growth.In summer, increased cloudiness likely blocks some of the radiation, acting to reduce surface heating, while the precipitation leads to larger surface moisture availability, further reducing the sensible heat flux since a larger fraction of the energy that does reach the surface is used for evaporation.The resulting decrease in surface sensible heat flux then leads to a lower PBLH.As the precipitation decreases again in September, the surface sensible heat flux increases somewhat and may thereby permit the small autumn peak in PBLH.The tendency for PBLH to peak during the dry season and drop during the summer monsoon season is common in tropical to subtropical latitudes (<35 ), especially near monsoon action centres (Chan and Wood, 2013).Chan and Wood (2013) also attribute this seasonal pattern to the effects of water availability on the partitioning of surface energy into latent and sensible heat fluxes.In addition, they note that while latent heat flux can also drive mixing through the production of thermals, these tend to vent mass out of the PBL rather than add to its build up during the wet season (Chan and Wood, 2013).Moreover, the cycle may be affected by seasonal variations in atmospheric stability and jet stream position as argued by Chen et al. (2016), who found that the local surface fluxes were not enough to explain the high winter and spring PBLH at the Gerze station (which is located within the CTP region as the regions are defined here).In their study, weak atmospheric stability and great PBLH (in ERA-Interim) were associated with high upper-level potential vorticity, corresponding to a strong and southerly positioned jet stream.During the summer when the jet stream is weaker and positioned farther to the north, the atmospheric stability is not as weak due to weaker potential vorticity and the PBL does not grow so deep.
In the eastern region (ETP), PBLH starts increasing in early spring and reaches its greatest heights in May and June, with medians just under 600 m (Figure 7c).The PBLH remains relatively stable throughout summer and then falls to reach its minimum in December, when the median is just over 100 m.The surface latent heat flux and the precipitation increase gradually during spring and reach their highest values in summer.The surface sensible heat flux also remains relatively high during summer, although its highest values are found in May and June (Figure 8e,f).The early summer PBLH thus seems to be controlled by the surface sensible heat flux.In July and August, not only insolation but latent heat release from summer rainfall could contribute to nearsurface heating (Yanai et al., 1992) and possibly lead to the growth of PBLH.As for the seasonal variation in statistical dispersion, the summertime PBLH exhibits the smallest interannual variability, while the winter and early spring have a considerably larger spread.Outliers are values that are found more than 1.5 times the interquartile range away from the top or the bottom of the box.The areal extent of the regions are shown in Figure 5e daily cycle of insolation.In the WTP region in winter and spring (Figure 9a,b), although many outliers reach above 1,000 m, the median PBLH displays a small diurnal variation with the daily peak centred around 15:00 local time (UTC +6).The intra-daily differences become more pronounced during summer when the values climb from just a few tens of metres in the early morning to around 1,250 m in the afternoon (Figure 9c).In autumn the median values reach considerably lower and peak around 700 m, while the extreme values are found over 2,000 m, which is only slightly lower than in the summertime (Figure 9d).

| Diurnal cycle of planetary boundary layer height
In the CTP region, the median value of the daily peak, which is generally reached at 15:00 (UTC + 6) in all seasons, ranges from just over 2,000 m in summer (Figure 9g) to over 3,700 m in spring (Figure 9f).In the winter (Figure 9e) the variability in daytime PBLH is large, with afternoon values ranging from about 10-5,000 m.For those high PBLHs, previous studies have shown that large-scale flows such as westerlies could enhance the growth of PBLH up to 5,000 m above ground (Chen et al., 2013;Lai et al., 2021b).The highest individual values are found in spring, when there are occurrences of PBLHs reaching around 6 km above ground level.As for the summer, Che and Zhao (2021) studied the summer PBLH at the TP and also found afternoon PBLH (at 12:00 UTC + 6) of 2,000 m in western TP (roughly corresponding to CTP as the regions are defined here).
Finally, in the ETP region, the diurnal cycles in all seasons are characterized by a comparatively rapid decline after the daily peak, which is reached by 13:00-14:00 (UTC + 6).The intra-daily differences are smallest in winter (Figure 9i) when the median PBLH stays below 530 m throughout the day.During summer in contrast (Figure 9k), the median PBLH ranges from below 100 m in the nighttime to almost 1,500 m in the afternoon.The afternoon values are again comparable to Che and Zhao (2021), who found that the eastern TP PBLH is around 1,500 m at 12:00 (UTC + 6).In spring (Figure 9j), although the median PBLH is lower, there are occurrences of afternoon PBLHs well over 3,000 m.
As for the difference in summer afternoon PBLH between the CTP and ETP, Che and Zhao (2021) points to west-east differences in sensible heat flux, downward solar radiation, soil moisture and cloud cover as important factors.In the western parts (CTP), low cloud cover and high downward solar radiation at the dry, bare soil surface promotes high sensible heat flux and thereby PBL growth.In contrast, eastern TP is characterized by higher cloud cover, less downward solar radiation and different F I G U R E 9 Diurnal cycles of hourly 1979-2018 ERA5 planetary boundary layer height (PBLH) for a western (left column), a central (middle column) and an eastern (right column) TP region in winter (top row), spring (second row), summer (third row) and autumn (bottom row).Each box gives the median (central line in box), the 25th percentile (bottom of box), the 75th percentile (top of box), the most extreme value that is not considered an outlier (the whiskers) and possible outliers (circles).Outliers are values that are found more than 1.5 times the interquartile range away from the top or the bottom of the box.The areal extent of each region is shown in Figure 5e surface properties (alpine meadow or steppe with higher soil moisture than the dry bare soil further west), leading to smaller PBLH.However, Che and Zhao (2021) only considered the summer PBLH, and the same mechanisms may not determine PBL growth in other seasons.For example, as noted above, Chen et al. (2016) found that the great wintertime PBLH (at Gerze) is related to weak atmospheric stability influenced by large-scale circulations from above.Previous studies also show that the growth of PBLH was generally later than the sensible heat flux due to the accumulation of energy (Zhang et al., 2011).Besides, the development of PBLH could also be eased by deep convection related to interaction between mountain-valley circulation and evaporative cooling from precipitation during afternoon to evening (Yang et al., 2004).The peak summer rainfall amount is overall earlier in the ETP than that in the CTP (Xu and Zipser, 2011;Kukulies et al., 2020), which is similar to the pattern of PBLH.Overall, this section of our study identified that the PBLH peak of the climatological diurnal cycles is earlier in the ETP than that in the CTP as well as WTP, and this pattern does not vary much in seasons.

| Uncertainties
When interpreting trends from reanalysis data, it should be noted that the amount, quality and accuracy of the assimilated observations is not constant over time (Bao and Zhang, 2019).The longer-term evolution inferred from reanalysis data for an observation-scarce region like the TP is therefore highly uncertain, and at the same time, the lack of observations makes it difficult to confirm any trends in the reanalysis data.
Differences in the methods used to identify the PBLH also introduce uncertainty and a potential source of disagreement between datasets.To estimate the PBLH from the radiosonde data, we applied a bulk Richardson approach, which is similar to the method used in ECMWFs IFS, as well as a parcel method.However, we should also note that there are still different ways to determine the PBLH from both observation and reanalysis data based on temperatures, humidity and refractivity, resulting in variability in spatiotemporal patterns of PBLH (von Engeln and Teixeira, 2013).As indicated by our comparisons in Section 3.2, the methods applied in different reanalyses can lead to large differences in the absolute values of the PBLH (Figure 4).Even using the same method, the differences in the simulated fields (e.g., temperature) that are used to derive PBLH are essential to the PBLH differences among reanalyses (Wang et al., 2020a).The PBLH climatological patterns revealed by ERA5 explored here thus need to be taken as a first step and further studies are needed.

| CONCLUSIONS
This study evaluated ERA5 PBLH with radiosonde measurements and investigated the spatiotemporal variations of 40-year seasonal and diurnal PBLH over the TP.ERA5 performs reasonably well in reproducing the diurnal cycles of radiosonde-derived PBLH at three out of four investigated locations in the western-central TP.At the fourth station, Qoms, which is located in the Himalayas at the southern TP, the agreement between the estimates from ERA5 and in situ observations is lower, possibly due to the complex terrain with large differences in altitude and surface properties including the complex topography at sub-grid cell horizontal scale of the model used to produce the reanalysis data.In addition, a comparison between the TP-average seasonality and trends in the PBLH from ERA5, CFSR and MERRA-2 showed that there are large differences among the three reanalyses in terms of the absolute heights and interannual variability.
The following results are obtained for the climatology based on ERA5: • The 1979-2018 PBLH varies greatly with time as well as location within the TP.On average, the PBLH is higher over the interior of the plateau with a median height of about 1,100 m and lower towards its rims as well as its westernmost (~280 m) and easternmost (~400 m) regions, with more pronounced spatial differences during winter and spring.• Overall, the seasonal and annual trends in PBLH are small and insignificant over the study period.TPaverage PBLH trends range from −2.6 m decade −1 in summer to +5.2 m decade −1 in winter.As for the spatial distribution, the magnitude and direction of the trends vary between regions and seasons, but statistical significance is lacking except in a few small regions with increasing PBLH.• The interannual PBLH variations are positively correlated with the strength of the NAO in western and southeastern TP in winter, eastern TP in summer, and central TP in autumn.A negative correlation is seen for northeastern TP in winter.Besides, the interannual changes in the ISM are negatively associated with the summertime PBLH in the central and southwestern TP. • The seasonal cycles show very different characteristics in different parts of the plateau, with a late-summer peak in WTP, a spring peak in CTP, and an earlysummer peak in ETP.In the CTP, the summertime minima in PBLH appears to be related to precipitation and cloudiness acting to suppress surface sensible heat fluxes and PBL growth.
• The diurnal cycles show that the PBLH can reach extremely high for individual hours.This is especially true for the CTP region, which features occurrences of high afternoon peaks in all seasons, with the maximum PBLH (around 6 km) arising in spring.The PBLH reaches its peak at about 15:00 (UTC + 6) in the CTP and WTP region while the peak is earlier at 14:00 (UTC + 6) in the ETP region.The largest range of PBLHs tend to occur at the PBLH peak time.
The Tibetan Plateau (TP) and the radiosonde stations.The TP is defined here as the region that spans 70 -105 E, 27 -40 N and has an altitude of more than 3,000 m above sea level.The colormap shows elevation (surface geopotential height, m) from the ERA5 reanalysis.The circles mark four radiosonde stations: Ali (79.70 E, 33.39 N), Shiquahne (80.10 E, 32.50 N), Gerze (84.03 E, 32.17 N) and Qoms (86.95 E, 28.36 N).At these stations, data were collected during a few different intensive observation periods between 2008 and 2019 [Colour figure can be viewed at wileyonlinelibrary.com]

F
I G U R E 4 Seasonal cycles of 1980-2010 monthly mean planetary boundary layer height (PBLH) from ERA5 (no hatching), CFSR (single line hatching) and MERRA-2 (crossed line hatching).Each box gives the median (central line in box), the 25th percentile (bottom of box), the 75th percentile (top of box), the most extreme value that is not considered an outlier (the whiskers) and possible outliers (circles).Outliers are values that are found more than 1.5 times the interquartile range away from the top or the bottom of the box.The dotted black lines show the medians and are added to make the seasonal pattern of variability clearer F I G U R E 5 ERA5 planetary boundary layer height (PBLH) at the Tibetan Plateau.In the left column, the median PBLH for the entire 1979-2018 period is displayed for winter (December-February; a), spring (March-May; b), summer (June-August; c), autumn (September-November; d) and the entire year (e).In the right column, linear trends in PBLH are shown for the same periods.Black lines are drawn over regions where the trend is not statistically significant at the 95% level.The red boxes on subplot e show the extent of the western (WTP), central (CTP) and eastern (ETP) TP regions for which seasonal and diurnal data are further explored trends of up to 60-70 m decade −1 in central TP (Figure

F
I G U R E 6 Correlation coefficient between the detrended ERA5 planetary boundary layer height (PBLH) and the North Atlantic Oscillation (NAO) indices in (a) DJF, (b) MAM, (c) JJA, and (d) SON as well as correlation coefficient between the detrended (e) PBLH (JJA) and East Asian Summer Monsoon (EASM), (JJA) indices, and (f) PBLH (JJA) and Indian Summer Monsoon (ISM, JJA) indices during 1979-2015.Black dots represent correlations significant at the 95% confidence level according to a Student's t test The climatological diurnal PBLH cycles for each of the three regions are presented in Figure9.The maximum values are reached in the afternoon, as expected given the Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Seasonal cycles of 1979-2018 monthly mean ERA5 planetary boundary layer height (PBLH) for a western (a), a central (b) and an eastern (c) TP region.Each box gives the median (central line in box), the 25th percentile (bottom of box), the 75th percentile (top of box), the most extreme value that is not considered an outlier (the whiskers) and possible outliers (circles).

F
I G U R E 8 Seasonal cycles of 1979-2018 monthly mean ERA5 total precipitation, surface net solar radiation and surface sensible and latent heat fluxes.Precipitation (bars, left axis) and surface net solar radiation (lines, right axis) are shown in the left panel, surface sensible heat flux (solid lines, left axis) and surface latent heat flux (dashedlines, right axis) are shown in the right panel.The cycles are for the western (top panel), central (middle panel) and eastern (bottom panel) TP regions.The areal extent of the three regions are shown in Figure 4e [Colour figure can be viewed at wileyonlinelibrary.com] Qoms Jun 2014 r F I G U R E Planetary boundary layer height (PBLH) at the Tibetan Plateau for different observation periods.The PBLH is shown for the radiosonde stations Gerze (a, b and c; 32.17 N, 84.03 E), Ali (d and e; 79.70 E, 33.39 N), Qoms (f, g and h; 86.95 E, 28.36 N) and Shiquanhe (i, j and k; 80.10 E, 32.50 N).For each station and observation period, hourly PBLH from the nearest grid cell in the ERA5 reanalysis is shown assolid lines.Radiosonde-derived PBLH is shown as rings (parcel method) and stars (bulk Richardson method).The locations of the radiosonde station are shown in Figure 1 [Colour figure can be viewed at wileyonlinelibrary.com] Summary statistics for the agreement between planetary boundary layer height (PBLH) from radiosonde data and the ERA5 reanalysis.
T A B L E 1