A merged continental planetary boundary layer height dataset based on high-resolution radiosonde measurements, ERA5 reanalysis, and GLDAS

The planetary boundary layer (PBL) is the lowermost part of the troposphere that governs the exchange of momentum, mass and heat between surface and atmosphere. To date, the radiosonde measurements have been extensively used to estimate PBL height (PBLH); suffering from low spatial coverage and temporal resolution, the radiosonde data are incapable of providing a diurnal description of PBLH across the globe. To fill this data gap, this paper aims to produce a temporally continuous PBLH dataset during the course of a day over the global land by applying machine learning algorithms to integrate high-resolution radiosonde measurements, ERA5 reanalysis, and the Global Land Data Assimilation System (GLDAS) product. This dataset covers the period from 2011 to 2021 with a temporal resolution of 3 h and a horizontal resolution of $\mathrm{0.25}{}^{\circ }\times \mathrm{0.25}{}^{\circ }$. The radiosonde dataset contains around 180 million profiles over 370 stations across the globe. The machine learning model was established by taking 18 parameters derived from ERA5 reanalysis and GLDAS as input variables, while the PBLH biases between radiosonde observations and ERA5 reanalysis were used as the learning targets. The input variables were presumably representative regarding the land properties, near-surface meteorological conditions, terrain elevations, lower tropospheric stabilities, and solar cycles. Once a state-of-the-art model had been trained, the model was then used to predict the PBLH bias at other grids across the globe with parameters acquired or derived from ERA5 and GLDAS. Eventually, the merged PBLH can be taken as the sum of the predicted PBLH bias and the PBLH retrieved from ERA5 reanalysis. Overall, this merged high-resolution PBLH dataset was globally consistent with the PBLH retrieved from radiosonde observations in terms of both magnitude and spatiotemporal variation, with a mean bias of as low as −0.9 m. The dataset and related codes are publicly available at https://doi.org/10.5281/zenodo.6498004 (Guo et al., 2022), and are of significance for a multitude of scientific research endeavors and applications, including air quality, convection initiation, climate, and climate change, to name but a few.

1 Introduction

The planetary boundary layer (PBL), the lowermost part of the troposphere where turbulence and convection mainly occur, is of significance in modulating the exchange of momentum, heat, moisture, and mass between the surface and the free atmosphere over a range of scales (Stull, 1988; Cooper and Eichinger, 1994; Edson et al., 2013). Turbulence in the PBL is largely generated mechanically owing to both wind shear and friction, whereas convectively, it is generated by buoyancy and surface heating (Degrazia et al., 2020). Within the PBL, vertical turbulent mixing of air masses is rapid and constant, in the order of 30 min or less (Wallace and Hobbs, 2006). Therefore, reliable parameterization of the PBL is crucial for accurate representations of vertical diffusion, cloud formation and development, and pollutant deposition in numerical weather prediction (NWP), climate, air quality, and coupled atmosphere–hydrosphere–biosphere models (Seibert et al., 2000; Hu et al., 2010; Baklanov et al., 2011). It is well recognized that the variation of PBL height (PBLH) significantly impacts the near-surface air quality (Petäjä et al., 2016; Wang and Wang, 2016; Lou et al., 2019; Li et al., 2021) and climate system as well (Esau and Zilitinkevich, 2010; Davy and Esau, 2016).

Development of the PBL is subject to changes in the energy balance near the ground surface, largely through the linkages between soil moisture and sensible heat flux, latent heat flux, and net radiation (Dirmeyer et al., 2014; Xu et al., 2021). In particular, the sensible heat flux is closely associated with the variation in evapotranspiration, land type, and cloud cover. Also, the daytime convective PBL is modulated by cloud radiative effects, particularly in the early afternoon (Guo et al., 2016; Zhang et al., 2018; Davis et al., 2020). Furthermore, the aerosol radiative effect (due to both aerosol scattering and absorption) indirectly affects evolution of the PBL by changing the atmospheric heating rate and the solar radiation reaching the surface (Wang et al., 2013; Li et al., 2017; Yang et al., 2016). Besides, the entrainment of air from above the PBL can also significantly drive evolution of the PBL (Hu et al., 2010).

To date, a variety of methods have been applied to vertical profiles of aerosol properties, water vapor, temperature, refractivity, and wind to estimate PBLH (e.g., Holzworth, 1964; Seibert et al., 2000; Lammert and Bösenberg, 2006; McGrath-Spangler and Denning, 2012; Chan and Wood, 2013; Su et al., 2018; Liu et al., 2019; Ding et al., 2021). The estimates vary considerably with data sources, algorithms, and vertical data resolution (Seibert et al., 2000; Seidel et al., 2010). For instance, PBLH determined by the minimum vertical gradient relative humidity is about 1 km larger than that from the parcel method, even though the latter algorithm is generally thought to be one of the most reliable methods for estimation of convective boundary layer (CBL) height (Hennemuth and Lammert, 2006; Seidel et al., 2010). In addition, different data sources, such as ceilometer lidar, COSMIC GPS RO satellite, radiosonde, and the fifth-generation ECMWF (European Centre for Medium-Range Weather Forecasts) atmospheric reanalysis (ERA5) dataset can reach quite different estimates of PBLH (Saha et al., 2022). Recently, as suggested by Teixeira et al. (2022), the PBLH should ideally be estimated using direct observations of vertical profiles of turbulent quantities, which is due in a large part to the turbulent nature of the PBL. However, only a few places have such observations. A wide range of complex physical and chemical processes involved in the PBL contribute to make PBLH estimates quite elusive and tricky (Seidel et al., 2010; Teixeira et al., 2022).

Table 1 Basic information on the data used in the present study, including data source, number of stations, vertical resolution, and the years with data curation.

Figure 1 Spatial distribution of sample number (color circles) for each radiosonde station at 00:00 (a), 06:00 (b), 12:00 (c), and 18:00 UTC from the years 2011 to 2021. Stations with less than 10 samples are not indicated.

Among the instruments, radiosonde is the most accepted instrument for deriving PBLH for both CBL and the stable boundary layer (SBL), due to its unprecedented capability to provide in situ observations of the thermodynamic and dynamic states of the PBL (Seidel et al., 2010; de Arruda Moreira et al., 2018, Guo et al., 2019). In addition, the bulk Richardson number method has been proven to be the most suitable PBLH algorithm for application to a large radiosonde dataset (Seidel et al., 2012). The dataset with full vertical resolution (5–8 m) has previously been used to study PBLHs over China and the near globe (Guo et al., 2016, 2021). The limitation of this dataset is its poor coverage over the ocean and some continental areas without high-resolution radiosonde observations.

By contrast, reanalysis datasets, such as the ERA5 reanalysis and the Modern-Era Retrospective Analysis for Research and Applications version 2 product (MERRA-2), have a unique advantage in terms of spatiotemporal coverage. Our recent study (Guo et al., 2021) suggests that ERA5 is the most promising reanalysis data source in terms of characterizing the evolution of PBLH, with an underestimation of daytime PBLH of around 130 m when compared to high-resolution radiosonde. Nevertheless, the underestimation of PBLH in the ERA5 reanalysis can be as high as 500 m in the afternoon, when the PBL is fully developed. This underestimation could be attributed to, but not limited to, the gradient of terrain elevation and the lower tropospheric stability. Particularly a higher terrain gradient or a more unstable troposphere generally lead to a lower PBLH in ERA5 reanalysis.

By exploiting both the advantages of in situ atmospheric measurements from radiosonde and the high-resolution model products from ERA5 reanalysis, it is desirable to generate a new PBLH dataset by seamlessly blending these versatile products. The biases between PBLHs retrieved from the ERA5 and radiosonde could be represented by the land properties and near-surface meteorological conditions, among other things, and be minimized or optimized via a machine learning model. The Global Land Data Assimilation System (GLDAS) incorporates satellite- and ground-based observations and produces a global, high-resolution product regarding land states and fluxes (Rodell et al., 2004). To this end, the present analyses used the radiosonde dataset that contains around 180 million profiles over 370 stations across the world in combination with the ERA5 reanalysis and GLDAS data. A long-term merged PBLH dataset covering the period 2011 to 2021 was generated, which could have crucial implications for the development and evaluation of weather and climate, environmental meteorology, and boundary layer parameterization. The rest of the paper is organized as follows: Sect. 2 describes the fundamental datasets and the PBLH methodology we used in this study, Sects. 3 and 4 report on the machine learning algorithm used to generate the merged PBLH dataset and also reveal the data quality, Sect. 5 represents the climatological merged continental PBLH, and Sect. 6 ends with a brief summary and conclusion.

Figure 2 Evolution of the difference between PBLH_ERA5 and PBLH_RS at various time scales: different years (a), different seasons (b), and at different local times (c). MAM, March–April–May; JJA, June–July–August; SON, September–October–November; DJF, December–January–February. The mean bias is labeled in the upper right corner of panel (a). Note that the Southern Hemisphere DJF (JJA) is combined with Northern Hemisphere JJA (DJF).

2 Data sources and conventional PBLH determination method

2.1 High-resolution radiosonde measurements

As described in Guo et al. (2021) and Zhang et al. (2022), a high-resolution radiosonde dataset gained from several organizations was adopted, spanning the years from 2011 to 2021. The organizations include the China Meteorological Administration (CMA), the National Oceanic and Atmospheric Administration (NOAA), the Global Climate Observing System (GCOS) Reference Upper-Air Network (GRUAN), the Centre for Environmental Data Analysis of the United Kingdom (CEDA), the University of Wyoming, and the Deutscher Wetterdienst (German Weather Service). Detailed information on the provided data is listed in Table 1. In total, over 185 million radiosonde profiles were collected to determine PBLH, 95 % of which were released at regular synoptic times of 00:00 and 12:00 UTC, and the rest of which were irregularly launched at other times during the intensive observational periods. Note that those soundings with the lowest burst height, lower than 10 km above ground level (a.g.l.), were eliminated. In addition, all the original soundings were evenly interpolated into the profiles with a vertical resolution of 10 m by cubic spline interpolation.

Table 2 Summary of input parameters of machine learning algorithms and the corresponding statistical metrics for their correlation analyses with PBLH bias between radiosonde and ERA5 reanalysis, including correlation coefficient and confidence level.

Figure 3 Joint distribution of the difference in PBLH_RS and PBLH_ERA5 and the surface sensible heat flux (a), the lower tropospheric stability (b), transpiration (c), and the near-surface temperature (d). The box-and-whisker plots in 10 even intervals are overlaid in each panel, and the correlation coefficients are marked in the upper right corner of each panel, wherein the asterisks indicate that the values are statistically significant (p<0.05).

The spatial distribution of sample numbers over each radiosonde station at four different synoptic times (00:00, 06:00, 12:00, 18:00 UTC) is presented in Fig. 1. It is noticeable that the radiosonde stations over Europe, the US, China, and Australia have an unprecedented, rich geographic coverage. Furthermore, the radiosonde measurements over China and the US have a fair temporal continuity at 00:00 and 12:00 UTC, with a total sample size reaching up to 3000 for each station. In comparison, the stations are poorly distributed over regions or countries such as South America, the Pacific islands, Russia, the Middle East, India, and Africa.

2.2 ERA5 and GLDAS

ERA5 is the latest version of the ECMWF reanalysis, benefiting from a decade of developments in model physics, core dynamics, and data assimilation (Hersbach et al., 2020). The PBLH product is resolved by the ERA5 reanalysis on a 1440×721 longitude and latitude grid, with a spatial resolution of $\mathrm{0.25}{}^{\circ }\times \mathrm{0.25}{}^{\circ }$ and a temporal resolution of 1 h, which is realistically simulated by the bulk Richardson number method. In addition, the parameters, such as the lower tropospheric stability (LTS), the standard deviation of digital elevation model (SDDEM), 10 m surface wind speed, 2 m air temperature, and 2 m pressure, are either computed or directly extracted from the ERA5 reanalysis. LTS is defined as the difference in potential temperature between 700 and 1000 hPa (Guo et al., 2016). As a result, a total of six parameters were obtained from the ERA5 reanalysis.

The land property parameters were taken from GLDAS, which include downward shortwave radiation (DSWR), downward longwave radiation (DLWR), surface heat net flux (SHF), surface latent heat net flux (LHF), evapotranspiration, transpiration, soil moisture in 0–10, 10–40, 40–100, and 100–200 cm, and total precipitation. In total, 11 parameters were extracted from the GLDAS product. GLDAS has a temporal resolution of 3 h and the same spatial resolution as the ERA5 reanalysis. However, GLDAS has no data over Antarctica. It should be noted that there exists a 0.125^∘ lag between the start latitude and longitude of GLDAS and those of ERA5 and, therefore, the latitude and longitude of GLDAS minus 0.125^∘ have to be used to match with ERA5 reanalysis.

According to the methods proposed by Guo et al. (2021), the collocation procedures between the grid products from ERA5 and GLDAS and station-based radiosonde observations were mainly implemented as follows: (1) the grid should contain the radiosonde station and (2) the UTC time (hour) of the grid product and radiosonde stay the same.

Table 3 Basic information on evaluation indices. MSE, mean squared error; RMSE, root mean square error; ABSmean, mean of the absolute bias; STD, standard derivation; RMS, root mean square.

Figure 4 Similar to Fig. 3, but for the difference between PBLH_RS and PBLH_merged.

2.3 PBLH determination using the bulk Richardson number method

The bulk Richardson number (Ri) is widely used for climatological study of PBLH from radiosonde measurements thanks to its applicability and reliability for all atmospheric conditions (Anderson, 2009; Seidel et al., 2012). Ri, a good indicator of turbulence and thermodynamic stability, is calculated as the ratio of turbulence due to buoyancy to that due to mechanical shear, which is formulated as

$$\begin{array}{}\text{(1)}& Ri\left(z\right)={\displaystyle \frac{\left(\frac{g}{{\mathit{\theta }}_{\mathrm{v}s}}\right)\left({\mathit{\theta }}_{\mathrm{v}z}-{\mathit{\theta }}_{\mathrm{v}s}\right)z_{\mathrm{AG}}}{{\left(u_z-u_s\right)}^{\mathrm{2}}+{\left(v_z-v_s\right)}^{\mathrm{2}}+\left(b{u}_{\ast }^{\mathrm{2}}\right)}},\end{array}$$

where g is the gravitational acceleration, z_AG the AGL, θ_v the virtual potential temperature, u_∗ the surface friction velocity, u and v the horizontal wind component, and b the constant which is usually set to zero since friction velocity is much weaker compared to the horizontal wind (Seidel et al., 2012). The subscripts of z and s denote the parameters at z height above ground and the ground level, respectively.

The critical value of Ri(z) can be used to identify a statically stable layer atop the PBL (Seibert et al., 2000), and it is commonly taken as 0.25. Meanwhile, PBLH estimates were found to vary little by differing the input of critical values ($Ri=\mathrm{0.2};\mathrm{0.25};\mathrm{0.3}$) (Guo et al., 2016). Therefore, the PBLH here is identified as the interpolated height at which the Ri(z) profile crosses the critical value of 0.25. The determined PBLH was set invalid in the following two scenarios: (1) Ri(z) in Eq. (1) exceeds 0.25, where z is the second level of radiosonde measurement; (2) the estimated PBLH is extremely high (for instance, 10 km) and apparently overestimates the well recognized climatological averaged PBLH.

3 Methodology

As shown in Fig. 2, there exist discernable biases between PBLH retrieved from radiosonde (hereinafter referred to as PBLH_RS) and PBLH determined from ERA5 reanalysis (hereinafter referred to as PBLH_ERA5). The match procedures between PBLH_RS and PBLH_ERA5 follow Guo et al. (2021). Noticeably, the PBLH bias (PBLH_RS minus PBLH_ERA5) is less dependent on years, with a mean bias of 95.7 m, indicative of a possible systematic PBLH underestimation by the ERA5 reanalysis. By contrast, the underestimation is around 137 m during the daytime (Guo et al., 2021), which is systematically larger than that during all days obtained in the present study. However, the bias is found to vary with season and local solar time (LST). More precisely, the mean bias varies from 150 m in March–April–May (MAM) to 64 m in September–October–November (SON), and from 309 m at 17:00 LST to 1.8 m at 00:00 LST. Moreover, the standard deviation of bias greatly changes from 64 m at 01:00 LST to 807 m at 17:00 LST. The large uncertainty raised by PBLH_ERA5 during the daytime motivated this study to establish a new PBLH dataset that would be more consistent with observations.

Figure 5 Spatial variations of PBLH differences between PBLH_RS and PBLH_merged. Panel (d) indicates the overall spatial distribution, and panels (c) and (e) illustrate its longitudinal and latitudinal variations, respectively. Panels (a), (b), (f), and (g) represent the seasonal variations over the four regions of interest, including North America, Europe, East Asia, and Australia, respectively. MAM, March–April–May; JJA, June–July–August; SON, September–October–November; DJF, December–January–February.

Previous studies have indicated that the bias could be physically attributed to variables such as SDDEM and LTS (Guo et al., 2021). However, the potential correlations with other variables, including DLWR, DSWR, SHF, LHF, evapotranspiration, transpiration, total precipitation rate (TPR), and soil moisture (SM), as well as wind speed, pressure, and air temperature at the near surface, have yet to be systematically investigated. Figure 3 shows that the bias is positively correlated with SHF, transpiration, LTS, and 2 m near-surface temperature, with a correlation coefficient ranging from 0.39 to 0.9 based on 10 evenly split bins. However, these parameters could be independent. For instance, evapotranspiration is determined by surface features which include plant physiology, land cover, and soil moisture, and it is the most important non-radiative process transmitting latent heat from the surface to the atmosphere (Cuxart and Boone, 2020). In addition, soil moisture probably contributes to decreases in the surface sensible flux locally (Basha and Ratnam, 2009). We further performed correlation analyses between the aforementioned variables and PBLH biases between radiosonde and ERA5 reanalysis, and the statistical results are shown in Table 2.

It was found that the PBLH bias is highly associated with variations in land properties, near-surface meteorological conditions, terrain elevations, LTS, and solar cycles. Consequently, it is possible to predict the PBLH bias based on these potential influential variables. Once the spatially resolved bias is available, a bias-corrected PBLH dataset, namely a merged PBLH product (denoted as PBLH_merged hereafter), can be acquired by perturbating PBLH_ERA5 with the addition of predicted bias. This process can be formulated as

$$\begin{array}{}\text{(2)}& {\mathrm{PBLH}}_{\mathrm{merged}}={\mathrm{PBLH}}_{\mathrm{bias}}+{\mathrm{PBLH}}_{\mathrm{ERA}\mathrm{5}},\end{array}$$

where PBLH_bias denotes the PBLH bias to be predicted. Under this philosophy, we here established a data-driven PBLH_bias prediction model with the abovementioned factors used as the potential input variables, while the PBLH bias over radiosonde sites was the learning target. Considering the possible dependence on the magnitude of PBLH_ERA5 and its corresponding LST, these two factors were also used as covariates in predicting PBLH bias.

After testing with several machine learning models, such as the ridge regression, the decision tree regressor, the support vector regressor, the multilayer perceptron regression, and random forest (RF), we found the latter method to give the most proper and robust prediction. Therefore, an RF regressor was established to give a prediction of PBLH_bias, and can be described as

$$\begin{array}{}\text{(3)}& \begin{array}{rl}{\mathrm{PBLH}}_{\mathrm{bias}}& =\mathrm{RF}(\text{DSWR, DLWR, LHF, SHF, EP, TP,}\\ & \text{SM10, SM40, SM100, SM200, TPR, PBLHE,}\\ & \text{LTS, SDDEM, NSP, NST, NSWS, LST})\end{array},\end{array}$$

where the abbreviation RF represents the random forest regressor, and the other acronyms and abbreviations are listed in Table 2. In the RF model, the hyper-parameters of the maximum depth of the tree and the random state of the bootstrapping of the samples are compiled to 20 and 5 in this analysis, respectively. The dataset that contains the input array and the learning target is randomly divided into two parts, with 70 % for training and 30 % for validation. All the data from 2011–2021 were included in the model training stage. Statistical metrics, including the mean squared error (MSE), root mean square error (RMSE), arithmetic mean, and arithmetic mean of the absolute difference, were applied to evaluate the performance of the prediction model.

4 Validation

Table 3a presents the prediction accuracy on the training and testing sets. Overall, the RMSE and arithmetic mean on the training subset are 243 and −0.2, respectively. In comparison, these two metrics are 370 and −2.8 on the testing subset, implying the presence of slight overfitting. To demonstrate the merit of PBLH_merged, we further compared the PBLH bias before and after merging. As illustrated in Fig. 4a, the mean bias between PBLH_RS and PBLH_merged is −0.9 m, which is smaller than the bias between PBLH_RS and PBLH_ERA5. In addition, the mean of absolute bias decreases from 260 m (PBLH_RS minus PBLH_ERA5) to 168 m (PBLH_RS minus PBLH_merged), and the standard derivation declines from 472 to 241 m, as listed in Table 3b. Moreover, the correlation coefficient between PBLH_RS and PBLH_ERA5 is 0.59, and it increases to 0.92 between PBLH_RS and PBLH_merged. More importantly, the bias between PBLH_RS and PBLH_merged during the daytime is dramatically decreased, to 20 m, compared to the bias between PBLH_RS and PBLH_ERA5 (300 m). These metrics clearly demonstrate a better accuracy of PBLH_merged than PBLH_ERA5, indicative of the merit of correcting modeling biases in PBLH_ERA5.

Figure 6 (a) Spatial distributions of PBLH_merged at 06:00 UTC across China for the years 2011 to 2021. The scatter plot in the left bottom of the panel illustrates the statistical correlation between PBLH_merged and PBLH_RS, where the asterisks indicate that the values are statistically significant (p<0.05). Also shown is the temporal evolution of annual average PBLH_merged, PBLH_RS, and PBLH_ERA5 during the period 2011 to 2021 (b).

Figure 7 Temporal variations of PBLH_merged (red) and PBLH_RS (blue) at Beijing (39.8^∘ N, 116.47^∘ E) (a), the Urad Zhongqi station (41.3^∘ N, 108.3^∘ E) (b) in the Nei Monggol Autonomous Region, and (c) the Heyuan station (23.7^∘ N, 114.7^∘ E) in the Guangdong province. Panels (b), (d), and (f) demonstrate the joint distributions of PBLH_RS and PBLH_merged, and correlation coefficients (R) and the fitted linear functions are given in the bottom right corner, where the asterisks indicate that the values are statistically significant (p<0.05).

Furthermore, the overview of PBLH bias (PBLH_RS minus PBLH_merged) in terms of spatial variation, and the seasonal variations over the four regions of interest are presented in Fig. 5. As compared to the finding in Guo et al. (2021), the bias dramatically decreases to dozens of meters for all stations (Fig. 5d), many of which slightly overestimate PBLH. More specifically, the PBLH over East Asia is overestimated by around 6 m (Fig. 5f), whereas it is underestimated by around 1 m over North America (Fig. 5a). Based on the bias with near-global coverage, we could infer that the merged model gives a more realistic PBLH estimate.

Figure 8 Spatial distribution of PBLH at 00:00 (a–d) and 12:00 UTC (e–h) in four seasons over land produced by the merged algorithms proposed here (i–l). The solid colored circles indicate the PBLH retrieved from high-resolution radiosondes. The shadow zones show nighttime regions, depending on the solar zenith angle on 15 April 2019 (MAM), 15 July 2019 (JJA), 15 October 2019 (SON), and 15 January 2019 (DJF). MAM, March–April–May; JJA, June–July–August; SON, September–October–November; DJF, December–January–February.

Intensive radiosonde observation is conducted across China in the boreal summer season at 06:00 UTC (14:00 Beijing Time), when the PBL is fully developed (Zhang et al., 2018). In addition to the overall near-global spatial distribution, a deeper investigation of PBLH_merged across China at 06:00 UTC is presented in Fig. 6. The spatial distribution of PBLH_merged exhibits a pronounced “northwest high southeast low” spatial pattern (Fig. 6a), which generally agrees with Zhang et al. (2018). The correlation coefficient between PBLH_merged and PBLH_RS is as high as 0.99, indicating their extreme consistencies in terms of spatial variations. The annual variations in PBLH_merged, PBLH_RS, and PBLH_ERA5 followed a similar trend, achieving a maximum in 2013 and a minimum in 2019 (Fig. 6b). The variations in PBLH_merged and PBLH_RS are rather close to each other. However, PBLH_ERA5 creates a different temporal variation, and it is systematically underestimated compared to PBLH_RS.

Figure 9 Spatial distributions of PBLH differences between the merged dataset and the ERA5 reanalysis from the years 2011 to 2021 at 00:00 (a), 06:00 (b), 12:00 (c), and 18:00 UTC (d).

As a good case in point for the comparison of fine structures, we show the diurnal variation of PBLH_merged and PBLH_RS at 06:00 UTC over three stations in Fig. 7. Three sites, including one in northwestern China where the highest PBLH is usually obtained, one in northern China where the most intensive observations can be found, and one in southern China where the lowest PBLH can be detected. The diurnal variations of PBLH_merged and PBLH_RS are strongly correlated with the lowest correlation of 0.88 (Fig. 7d). From Figs. 5–7, we can observe that the spatiotemporal variations of PBLH_merged and PBLH_RS are in good agreement.

5 Merged continental planetary boundary layer height

The climatological mean of PBLH_merged in four seasons at 00:00 and 12:00 UTC during the years 2011 to 2021 is illustrated in Fig. 8, and PBLH_RS at the same UTC and in the same season is overlaid as filled circles. At all UTCs and in all seasons, the PBLH_merged is high during the daytime as compared to the nighttime, and reaches a maximum of around 2 km, especially in the afternoon. In addition, PBLH_merged experiences noticeable seasonal variation. For instance, over Australia, the PBLH_ERA5 in SON and December–January–February (DJF) seasons is about 400 m larger than in the other two seasons (Fig. 8a–d), and vice versa in the Northern Hemisphere. Moreover, we observe that PBLH_merged has a clear latitude and elevation dependence: it decreases from approximately 2 km at low and middle latitudes to around 0.8 km at high latitudes during the daytime. At similar latitudes, PBLH_merged over terrain with a high elevation could be substantially larger than for that with a low elevation. For example, in the DJF season and at 00:00 UTC, the PBLH_ERA5 over the Andes Mountain is about 0.4 km higher than that over the surrounding flat region (Fig. 8d). In a short conclusion, the spatiotemporal variability of the PBLH_merged is inevitably associated with local times, seasons, latitudes, terrain elevations, and hemispheres.

In general, PBLH_merged is remarkably consistent with PBLH_RS in terms of seasonal variation and diurnal cycle, especially at 00:00 and 12:00 UTC, when the radiosonde measurement is comparatively sufficient. These findings suggest that PBLH_merged could adequately resolve the climatological variation in PBLH.

The difference in PBLH_merged and PBLH_ERA5 during the years 2011–2021 at four typical times is further illustrated in Fig. 9. Compared to PBLH_ERA5, PBLH_merged is overestimated overall, with a mean overestimation of approximately 90 m. The overestimation appears very close to the difference in PBLH_RS and PBLH_ERA5. The overestimation over North America at 00:00 UTC, over East Asia and South Asia at 12:00 UTC, and over Africa at 18:00 UTC can be as high as 500 m. However, PBLH over some areas, such as the Middle East at 06:00 UTC and the Western United States at 18:00 UTC, is slightly underestimated by around 200 m.

6 Data availability

The merged PBLH dataset and the related codes can be accessed at https://doi.org/10.5281/zenodo.6498004 (Guo et al., 2022). ERA5 data is publicly accessible at https://cds.climate.copernicus.eu/#/search?text=ERA5&type=dataset (ECMWF, 2019). NASA GLDAS can be accessed at: https://disc.gsfc.nasa.gov/datasets/GLDAS_NOAH025_3H_2.1/summary?keywords=GLDAS (NASA, 2021).

7 Conclusions and summary

The general underestimation of PBLH by the reanalysis dataset, especially during the daytime, motivated the present analysis to generate a merged long-term high-resolution seamless continental PBLH dataset, i.e., PBLH_merged, by integrating multi-modal data products including 185 million high-resolution radiosondes from the years 2011 to 2021, ERA5 reanalysis, and the GLDAS product. The PBLH_merged dataset generated in this study has a horizontal resolution of $\mathrm{0.25}{}^{\circ }\times \mathrm{0.25}{}^{\circ }$ and a temporal resolution of 3 h, identical to PBLH_ERA5 but with much higher data accuracy.

Compared to PBLH_RS, PBLH_merged is overestimated by around −0.9 m, which is considerably smaller than the bias between PBLH_RS and PBLH_ERA5 (95.7 m). During the daytime, the mean and the standard derivation of bias are remarkably decreased from 300 and 600 m (PBLH_RS minus PBLH_ERA5) to 20 and 300 m (PBLH_RS minus PBLH_merged), respectively. In addition, the climatological variation of the merged PBLH dataset is highly correlated with PBLH_RS, both in terms of magnitude and in terms of spatiotemporal variation. Moreover, the climatological mean of continental PBLH_merged is around 90 m higher than that of PBLH_ERA5, which is quantitatively consistent with the comparison result of PBLH_RS and PBLH_ERA5. Overall, the merged dataset closely agrees with the radiosonde-derived PBLH in terms of magnitude and spatiotemporal variation.

In conclusion, the PBLH_merged dataset is outstanding in terms of both spatiotemporal coverage and good accuracy. This dataset could be of importance for advancing our understanding of the PBL processes involved in air quality prediction, weather forecasting, and climate projection under global warming. In the future, with more data available over the ocean, a global seamless PBLH dataset is warranted, and this needs more field campaigns to be deployed over the open ocean or islands in the ocean in which more intensive radiosonde balloons are launched. Besides, it is imperative to improve the observational capability of satellite-based instruments in characterizing the temperature and humidity profiles in the PBL, which will no doubt help to fill the gaps in atmospheric sounding over the ocean.