Long-term trends of global maximum atmospheric mixed layer heights derived from radiosonde measurements

The height of the atmospheric mixed layer is a critical parameter controlling the vertical dispersion of pollutants. Here we calculate daily maximum mixed layer height (MMLH) using operational radiosonde and surface meteorological measurements made at 219 carefully selected WMO weather stations and analyze their long-term trends from 1973 to 2018. We found that 74 stations showed significant increases in MMLH, whereas 48 sites showed negative trends. Positive trends are mainly found in Central US, Europe, Africa, East and Southeast Asia and East Australia. Stations over the coastal US, India and West Australia generally exhibit negative trends. The trends can be attributed to changes in vertical temperature gradient ∇ θ v between 950 and 700 hPa and diurnal temperature range (DTR). ∇ θ v in general decreased and caused positive MMLH trends over the majority of the regions. DTR decreases globally, causing negative MMLH trends (corresponding to decreased DTR), and plays a more important role over India and Central Asia.


Introduction
The atmospheric mixed layer, also known as the convective boundary layer, is the lowest part of the atmosphere characterized by vigorous turbulence. It is usually formed in the early afternoon when the surface air is heated and becomes unstable. Because the vertical distribution of many atmospheric components such as water vapor, trace gases or aerosols, can be nearly uniform within the mixed layer, the mixed layer height (MLH) is a critical parameter controlling the surface concentration of these components. For example, the MLH has been found to be closely related to surface PM 2.5 concentrations (Tang et al 2016, Su et al 2018. Therefore, continuous monitoring of the MLH is very important in understanding the formation of air pollution, as well as improving parameterization of climate, weather and air quality models. The MLH can be obtained from several observational techniques including radiosondes, groundbased or space born lidars, climate towers and airplane surveys. There are also various methods developed to calculate the MLH from these measurements, such as the parcel method (Holzworth 1964), the maximum vertical gradient of potential temperature or aerosol backscatter (Hennemuth and Lammert 2006), the bulk Richardson number threshold (Richardson et al 2013), etc. Till now, most of previous studies of the MLH are regional and/or of short duration. Holzworth (1964) derived a climatology of MLH in the US using radiosondes. Seidel et al (2012) reported climatological boundary layer heights (BLHs) in Europe and the US. Guo et al (2016) constructed BLH climatology in China with both radiosondes and reanalysis data. Although Seidel et al (2010) gave a global estimate of BLH, the study period was only 10 years while no longer-term analysis is available. Zhang et al (2013) estimated the near 30 year BLH trends but only for Europe. Guo et al (2019) investigated the temporal trends of BLH in China and found a trend shift since ∼2004. Therefore, there is still lack of a study of longterm MLH trends and variability globally.
A major difficulty in estimating global long-term trends of the MLH is the lack of spatially and temporal resolved measurements of the atmospheric vertical profile. Radiosondes that are operationally carried out at worldwide WMO weather stations provide the most extensively covered vertical measurements with relatively long-time series. However, the standard launching times of radiosondes are 0000 UTC and 1200 UTC each day, which do not necessarily represent the fully developed mixed layer locally. For example, in East China these two times correspond to 0800 and 2000 local time respectively, when convection in the lower atmosphere is usually prohibited by strong temperature inversion due to nighttime cooling. In fact, because the maximum development of the daily MLH usually occurs in the early afternoon, the operational launch times are sub optimal for MLH related studies over most places. Holzworth (1964) proposed the parcel method to calculate daily maximum mixed layer height (MMLH) using morning and evening radiosondes and diurnal potential temperature observations to characterize the convective mixing of the lower troposphere. This method has the advantage of not requiring continuous sounding profiles during the entire day, and has been adopted by various studies to estimate the MLH thereafter (Seidel et al 2010, Du et al 2013, Yang et al 2013, Karimian et al 2016, Su et al 2017.
In this paper, we attempt to generate long-term estimates of daily MMLH using the parcel method globally and analyze its trends and variability. We also try to explain the trends by examining the changes of vertical temperature gradient and diurnal temperature variability. Although uncertainties still remain, our study is the first to document and to explain long-term global MLH changes, and we hope that it will facilitate the evaluation of climate change impacts on atmospheric structure and surface pollution, as well as provide references for model parameterization.

Data and method
2.1. Radiosonde data The radiosonde measurements are downloaded from the Integrated Global Radiosonde Archive (IGRA, https://ncdc.noaa.gov/data-access/weatherballoon/integrated-global-radiosonde-archive). We use the updated IGRA Version 2 dataset , Durre and Yin 2008, 2011 which contains more stations and parameters compared to the previous version. The soundings of pressure, temperature and dewpoint depression are used to calculate profiles of potential temperature (q) virtual potential temperature (q v )used for MLH estimation. Although radiosonde measurement at some stations can date back as early as 1905, the majority of the stations came into operation after 1973. We thus use data from 1973 to 2017 in our study.

Surface meteorological data
The surface meteorological measurements are obtained from the National Centers for Environmental Information (NCDC, http://1.ncdc.noaa.gov/ pub/data/noaa/) of the National Oceanic and Atmospheric Administration. The measurements are generally made in 3 h intervals. We use station pressure, temperature and dewpoint to calculate surface q and q .
v For the few occasions when no station pressure is reported, we use sea level pressure, temperature and station elevation to approximate it as follows (Jacob 1999): where P is station pressure, P 0 is sea level pressure, h is station elevation and T is temperature in K.

Data selection and MMLH estimation
Although there are more than 2000 radiosonde stations and even more surface stations, many of them do not have sufficiently long data records for trend estimation. Moreover, calculating MMLH using the parcel method requires matching of radiosonde and surface measurements and relatively frequent diurnal sampling. Therefore, we develop a set of careful data selection and quality control criteria to reduce analysis uncertainty as much as possible: 1. There should be at least four surface q v measurements during the day time; 2. Because local maximum q v typically occurs in the early afternoon, we require each valid day have at least two surface observations of pressure, temperature and dew point between 1100 and 1800 local time. Note that in the presence of clouds and precipitation, convection is no longer primarily driven by solar heating. We therefore remove all measurements when precipitation is observed or when the sky cover indicator shows 'overcast', 'obscure' or 'partial obscuration'; 3. There should be at least one q v sounding within 6 h of the daily maximum surface q . v When there are more than one sounding profiles, the one closest to the surface maximum q v in time is selected; 4. For the sounding(s) identified by step 3, there must be at least five levels below 5000 m above ground level (agl). Because the standard pressure levels consistent of approximately five measurements below 5000 m agl (made at 1000 hPa, 950 hPa, 850 hPa, 700 hPa and 500 hPa respectively), we set this value as the lower limit. As pointed out by Seidel et al (2010), MLH estimate is sensitive to the vertical resolution of the sounding data, which actually changes over time. We therefore only use the standard pressure levels to calculate MLH, even though for a particular sounding the vertical resolution can be much higher; 5. Using available data selected from the above four steps, we calculate daily MMLH by identifying the level at which q v is equal to surface daily maximum q . v We restrict our search within 5000 m, because MLH can rarely reach above this level. This only removes less than 0.1% from the total results. 6. For long-term trend analysis at a certain station, we require there be at least 10 days of data each month, at least 9 m of data each year and at least 30 years of data within the 45 year study period.
After processing the data according to these six steps, 219 stations remain for further analysis. The data coverage is reasonable for North America, Europe, Asia, South America, Australia and the islands, as can be seen in figure 1.
According to the parcel method (Holzworth 1964), nocturnal radiative cooling of the ground results in a stable lapse rate at night. In the daytime, as the Sun heats the surface, air mass is assumed to be lifted adiabatically. The level at which potential temperature matches surface maximum potential temperature is considered as the MMLH. Here we also consider moist mixing and uses virtual potential temperature (q v ) instead. The MMLH thus depends upon the vertical temperature structure and daily maximum q . v A spline interpolation is applied to the discrete surface temperature measurements to estimate daily maximum q .
v Note that only the effect of vertical convection is considered in this method, while several other factors, such as wind shear, horizontal mixing and mechanical turbulence is neglected. The parcel method has been shown to perform reasonably well (Seibert et al 2000, Seidel et al 2010. The seasonal climatology of derived global daily MMLH is shown in figure S1, available online at stacks.iop.org/ERL/15/ 034054/mmedia. We can see that in general, MLH is higher in low latitudes than in high latitudes, which well follows the pattern of solar heating. Southern hemisphere also has comparatively higher MLH. The island sites have lower MLH than inland sites to some extent reflecting the difference between land and marine boundary layers. For mid to high latitudes, higher MLH is observed during local spring (MAM for northern hemisphere and SON for southern hemisphere) and summer (JJA for northern hemisphere and DJF for southern hemisphere), whereas tropical sites exhibit little seasonal variation. Qualitatively, these patterns agree well with regional BLH climatology reported by Seidel et al (2012) for Europe and Guo et al (2016) for China. The seasonal patterns of MMLH also qualitatively agree with those derived using ceilometers over Shanghai (Peng et al 2017), London (Kotthaus and Grimmond 2018), Vienna (Lotteraner and Piringer 2016), and Vancouver (van der Kamp and McKendry 2010), and with that derived using micro-pulse lidar derived MMLH over Hong Kong (Yang et al 2013) and Beijing (Chu et al 2019). These studies suggest higher MMLH in the local summer and fall and lowest MMLH is found in winter. Our results in general appear higher than ceilometer or MMLH results, which should be attributed to differences in the observation and MMLH calculation method.
Another uncertainty is that at most sites, the sounding data are not obtained when the MLH reaches daily maximum. For example, for East Asia, the routine radiosondes are performed at ∼0800 and 2000 local time, whereas the daily MMLH typically occurs in the early afternoon. Regarding this issue, we examined the validity of using 0800 radiosonde measurements to estimate 1400 MLH by comparing the results with that calculated using soundings at exactly 1400 at selected sites in China during the summer time. All MLHs are calculated using the parcel method (Holzworth 1964). Figure S2 shows the comparison between the 1400 MLH estimated using 0800 and 1400 sounds respectively, which agree very well with a high correlation coefficient of 0.74. This supports the assumption that upper atmosphere q v profile changes little during the day and indicates that using soundings at 0800 can estimate daily MMLH reasonably well.

Results
To estimate long-term MLH trends, we use the Sen's slope method (Sen 1968) combined with Mann-Kendall test (MK test, Mann 1945, Kendall 1975) for significance test. These two techniques have been successfully used in trend estimation of atmospheric variables (Li et al 2014 and are less prone to outliers compared to linear regression. The multiyear averaged seasonal cycle is first removed from the MMLH time series for each station to reduce serial correlation. In the MK test, serial correlation is further removed using the pre-whitening process suggested by Yue et al (2002). The MMLH trends reported in m/ decade are shown in figure 1. Out of the 219 stations, 122 show statistically significant (p value<0.05) trends, in which 74 stations show positive trends. Increases of MMLH are mainly observed over Europe, East and South Asia, Central US, East Australia and many island sites. The trends at most stations are on the order of 100 m/decade. Stronger trends reaching 150 m/decade are observed in Northern hemisphere low latitudes, i.e. India, Southeast Asia, and a few sites in Africa and Central America. Given that the magnitudes of the MMLH over these regions are also larger, their relative changes are comparable to the other parts of the world. Negative trends are mainly found in Central Asia, East and West coasts of US and West Australia. Negative trends appear generally weaker than positive trends.
We also examined the seasonal MMLH trends in figure 2. Overall, the seasonal trend patterns are consistent with global trends. Both hemispheres exhibit the largest trend in local spring (MAM for the NH and SON for the SH). Also, more sites show significant trends in these two seasons. The seasonality of the trends is less obvious for the sites close to the equator ( i.e. South Asia). This is reasonable since at these sites, there is less seasonal variation of insolation and surface temperature.
Next, we will investigate the possible mechanisms causing these observed MMLH trends. As mentioned above, the MMLH depends on vertical q v profile and daily maximum q .
v Indeed, we find that de-seasonalized, de-trended MMLH is negatively correlated with 700/950 hPa q v gradient ( q  , v calculated as q q  - v v ,700 hPa ,950 hPa ) and positively correlated with diurnal temperature range (DTR, figure 3). The first relationship is obvious, as larger q  v indicates less buoyancy and thus lower convection height. The second one can be understood by the fact that stronger DTR is usually associated with stronger diurnal solar heating of the surface which increases the instability of surface air. DTR also reflects the surface q v difference between daily maximum value and the value at the radiosonde launch time. The impact of q  v on the MLH or BLH was also noticed by Guo et al (2016), whereas the effect of DTR is less documented. For clearer illustration, a schematic figure is shown in figure S3. It is seen that the increase of q  v will decrease MMLH whereas the increase of DTR serves to increase it. We also included the q  v profiles at three example sites showing positive, negative and neutral changes in q  v during the 1973-2018 period ( figure S4). According to the parcel method demonstrated by figure S3, the site with positive q  v change will show decreased MMLH and that with negative q  v will show increased MMLH. Nonetheless, there are still a few outliers in figure 3 that show positive correlations with q  v or negative correlations with DTR.
One reason is that the 950-700 hPa gradient may not fully represent the q v structure within the mixed layer. For example, the mixed layer top might lie below 700 Figure 2. Distribution of seasonal MMLH trends. Larger dots with black circle mean that the trend is significant at 95% level according the MK test.
hPa or well above it, or there are temperature inversions. Uncertainties in the measurement and associated with interpolation may also lead to the reversed correlations.
Compared with MMLH, more sites show statistically significant trends for q  v and DTR (figures S5 (a) and (b) respectively). The trends for the latter two parameters also appear more spatially uniform, especially for the DTR. Globally, DTR exhibit a prevailing decline over the past four decades, with stronger trends typically observed over India, East Asia and East Australia ( figure S5). q  v also shows more negative trends than positive trends. Decreases are found over most North America, South America, Europe, Central Asia, Central China and Australia. Several sites in India, northeast Asia and West Australia show increased q  v ( figure S5(b)). Because negative q  v trends cause increased MMLH, whereas a weaker DTR tends to cause decreases in the MMLH, the final MMLH trends thus reflect the competing result of these two effects.
To draw a more quantitative conclusion, we need to decompose the observed MMLH trends into those associated with these two effects respectively. Therefore, we perform two MMLH trend decomposition experiments, one using climatological q  v gradient and the other using climatological DTR. Specifically, for the first experiment, we first calculate a q  v profile climatology by averaging all q  v profiles for each month. We then scale this climatological profile with the observed q  v at the lowest level of each raw profile. This results in a set of q  v profiles with the same vertical gradient but still captures the temperature change with time (e.g. global warming signal). We then use this newly generated set of q  v profiles to recalculate the daily MMLH (hereafter referred as MMLH θv ) and its trends. Because q  v is fixed in this experiment, the MMLH trends represent those caused by the DTR. In the second experiment, the DTR is fixed at its climatology level. This is done by first calculating the multi-year averaged difference between daily maximum temperature and daily mean temperature for each month, then adding this value to the observed daily mean temperature for each day to obtain a set of new daily maximum temperature, and finally calculating the MMLH using this new daily maximum temperature (hereafter referred as MMLH DTR ) and its trends. Therefore, these new trends reflect the effect of temperature gradient changes on MMLH trend. Figure 4 shows the scatter plot between MMLH trend and the sum of q MMLH v and MMLH DTR trends. We can see that these two quantities agree reasonably well, with the majority of the dots distributed close to the 1:1 line and a R 2 value of 0.77 (correlation coefficient R=0.88). This indicates that the effects of q  v and DTR on the MMLH is linear, and thus justifies our trend decomposition analysis. Consistent with figure S5(a), q MMLH v also exhibit extensive negative trends throughout the globe ( figure 5(a)). The largest trends can reach −100 m/ decade and mostly observed over Indian subcontinent and Australia. These are also regions with large DTR trends. In contrast, there are more sites with MMLH DTR trends than negative trends ( figure 5(b)). MMLH DTR typically increased over Europe, East Asia, South America and East Australia, well corresponding to places with negative q  v trends ( figure S5(b)).
Compared with figures 5(a), (b) shows better resemblance with MMLH trends pattern shown in figure 2, indicating that changes in q  v might be primarily responsible for the MMLH changes at most stations. We further analyzed the relative contribution of these two effects and summarized the results into figure 6. The relationship between the sign and magnitude of the trends for MMLH, q MMLH v and MMLH DTR can be divided into four categories: (1) MMLH DTR and MMLH trend agree in sign and both are opposite to q MMLH v trend; (2) q MMLH v trend and MMLH trend agree in sign and both are opposite to MMLH DTR trend; (3) All trends agree in sign but MMLH DTR trend is greater than q MMLH v trend (4) All trends agree in sign but q MMLH v trend is greater than MMLH DTR trend. These four categories are respectively marked Figure 3. Correlation between MMLH and 950/700 hPa θ v gradient (calculated as 700 hPa θ v minus 950 hPa θ v , (a)) and between MMLH and diurnal range of θ v (calculated as the daily maximum θ v minus daily minimum θ v , (b)). Larger dots with black circles mean the correlation coefficient is significant at 95% level using Student's t test. Note the correlations are estimated using de-seasonalized and de-trended data, and serial correlation is taken into account when calculating the t statistic. Note that MMLH DTR trends reflect those caused by the change of θ v gradient whereas MMLH θv trends reflect those caused by the change of diurnal temperature range. by red, blue, orange, green dots on figure 6. Note that categories 1 and 3 both mean that changes in q  v contribute more to the observed MMLH trend, whereas for categories 2 and 4, changes in DTR play a more important role. It is clearly seen from figure 6 that the numbers of red and orange dots significantly outweigh those of blue and green dots (67% versus 33%). For the majority of the sites over Europe, the Americas, Africa, Central China and Australia, the MMLH trends can be primarily attributed to changes in q  , v and mostly to the decreased q  v ( figure S5(b)). Only for some East US, Central America, Central Asia and Indian sites, DTR decrease is mainly responsible for the decreased MMLH trends there. Also, for 62% of the sites, q MMLH v and MMLH DTR trends are opposite in sign (red and blue dots), indicating that temporal changes of q  v and DTR are causing different changes in convection and thus MMLH.
In addition, since we consider moist convection here, apart from q v gradient and DTR, changes of   water vapor can also be a potential factor. We thus examine its effect by calculating MMLH using q instead of q , v which means only dry convection is concerned. The trends of q calculated MMLH are shown in figure S6. It is clearly seen that these trends are almost the same as those shown in figure 2. Indeed, the scatter plot in figure S7 suggest that the two trends agree very well. This indicates that change of water vapor is a negligible factor.
The final word to mention is that apart from q  v and DTR considered here, daily MMLH can also be affected by other factors. In particular, since we assume the upper air virtual potential temperature profile does not change during the day, there is uncertainty in the daily MMLH estimation if there are sudden changes of q  v due to synoptic events (e.g. cold fronts, old/warm advection, etc). However, these events are occasional and should not seriously impact the long trend of daily MMLH.

Conclusions and discussion
In this paper, we use radiosonde and surface meteorological measurements from 219 carefully selected WMO weather stations to derive global long-term daily MMLH and analyze its trend and variability over the 46 year period from 1973 to 2018. To calculate daily MMLH, the parcel method is applied to the morning/afternoon virtual potential temperature profiles and surface maximum virtual potential temperature. Results show that the spatial distribution of MMLH trends displays certain regional variability, with more than half of the sites with significant trends exhibiting an upward MMLH trend. the spatially varying MMLH trends can generally be considered as the competing effect of decreased q  v (more convection) and decreased DTR (less convection).
Our current work still suffers from several sources of uncertainty. First, the parcel method for daily MMLH calculation assumes environmental atmospheric profile unchanged for several hours. We compared randomly selected profiles measured at 0000UTC and 1200UTC for each season at each site, and found while this assumption can largely hold, there are a number of cases with obvious differences. Second, as pointed out by previous study (Seidel et al 2010, Collaud Coen et al 2014, different observation and calculation techniques can yield different MLH results. However, selecting the 'best' method is case by case and sometimes not possible. On the other hand, given the limitation of the radiosonde launch frequency, our approach is currently the most feasible. Additionally, we noticed at many places, MMLH time series exhibit considerable interannual and decadal variability beyond the linear trend. These phenomena might be related to climate oscillations and deserves further study. Finally, the causes of q  v and DTR trends, likely associated with global warming but also modulated by regional scale climate variability, also needs investigation. At the current level, our work serves as the first to report a global MLH trend entirely based on observation, which is meaningful for explaining long term pollution trends and boundary layer parameterization in models.