A Method for Deriving the Boundary Layer Mixing Height from MODIS Atmospheric Profile Data

The planetary boundary layer is the medium of energy, moisture, momentum and pollutant exchange between the surface and the atmosphere. In this paper, a method to derive the boundary layer mixing height (MH) was introduced and applied over the Heihe river basin. Atmospheric profiles from the MODerate Resolution Imaging Sepctroradiometer (MODIS) instrument onboard the NASA-Aqua satellite were used for the high spatial resolution of this method. A gap-filling method was used to replace missing MODIS data. In situ MH data were also calculated from HIWATER (Heihe Watershed Allied Telemetry Experimental Research) and WATER (Watershed Allied Telemetry Experimental Research) observational radiosonde sounding data from 2008 and 2012 using the Richardson number method combined with a subjective method. The MH occurs where there is an abrupt decrease in the MR (water vapor mixing ratio). The minimum vertical gradient of the MR is used to determine the MH. The method has an average RMSE of 370 m under clear skies and convective conditions. The seasonal variation in the MH at the Gaoya radiosonde station is also presented. The study demonstrates that remote sensing methodologies can successfully estimate the MH without the help of field measurements.


Introduction
The planetary boundary layer (PBL) or atmospheric boundary layer refers to the part of the atmosphere that is directly influenced by the Earth's surface and responds to surface forces within a timescale of an hour or less [1,2].Because substances emitted into the PBL are gradually dispersed horizontally and vertically through the action of turbulence, if sufficient time elapses and there are no significant sinks, the PBL will become completely mixed; hence, the height of this layer is often called the mixing height.The mixing height (MH) is not only a key parameter for air pollution models, but is also an important parameter in the surface energy balance system (SEBS) for the estimation of turbulent heat flux.In the application of pollution models, the MH determines the volume available for the dispersion of pollutants and is involved in many predictive and diagnostic methods and/or models to assess pollutant concentrations [3].Similarly, evapotranspiration (ET) estimation models used for SEBS analysis require atmospheric data.Normally, these variables are estimated either from a large-scale meteorological model or from Obukhov similarity theory using data such as air pressure, temperature, humidity and wind speed at a reference height [4].
Two basic possibilities for the practical determination of the MH are its derivation from profile data (measurements or numerical model output) and its parameterization using simple equations or models, which only need a few measured input values.The traditional approach uses the Parcel and Richardson number method to analyze the air temperature profile from radiosonde data [3].Currently, however, an increasing number of various remote sensing instruments and technologies are used to estimate MH.Compared with the free atmosphere, the MH has a much higher concentration of water vapor and aerosols.The abrupt decrease in water vapor and aerosols can cause a sharp Lidar signal change around the top of the boundary layer.Thus, Lidar can be used to monitor the temporal and spatial variations in the MH over urban and other areas [5][6][7].The use of SODAR, WTR and ceilometer instruments to estimate the MH has also been discussed extensively in the literature [8].These three instruments deliver five types of vertical profiles: wind and acoustic backscatter from the SODAR (Sound detection and ranging) instrument, wind, temperature and electromagnetic backscatter from the WTR (Wind-temperature-radar) instrument, and optical backscatter from the ceilometer.These can be combined to derive the MH [9].
However, radiosondes, being launched only several times each day, with limited height resolution or routine ascents, do not provide continuous vertical profiles over time.Thus, it is difficult to study the time evolution of the MH [3].Furthermore, radiosonde data do not have good spatial representation since balloons released at one site drift with the wind.Since landuse is a crucial factor for MH evolution, this is a significant problem in urban areas, given their complex spatial variation in land use, as well as in basin-scale studies [10].Relying on ground-based remote sensing also limits the choice of the parameters that can be monitored.The classical models of turbulent transport in the atmospheric boundary layer require other input variables, such as wind speed, which are very difficult to measure accurately using land-based remote sensing tools.Hence, although the spatial structure of the MH is an important characteristic of the atmosphere, especially when the research area is at the basin-scale or larger, current methodologies cannot provide a clear indication of its spatio-temporal variation.
To date, only a few studies of limited scale have examined PBL processes using space-based remote sensing [11].However, as remote sensing has evolved, satellites providing atmospheric temperature and water-vapor profile data are increasingly available.Among them, MODIS has accumulated more than 10 years of global data (MOD07 and AIRS product).Some geostationary meteorology satellites, such as the Fengyun series, have a visitation frequency of half an hour [12], but their vertical resolution is limited.If space-based remote sensing data can detect the abrupt decrease of water vapor at the top of the PBL, and thus diagnose the MH, it will alleviate the lack of detailed spatio-temporal information for the MH.
The basic advantages of space-based systems are their continuous operation and the absence of any impact on the investigated flow.The objective of the present study is to develop a method for estimating a high spatial-resolution MH and associated variables as input for a 1-km resolution SEBS model for the Heihe River Basin using only atmospheric variables from space-based remote sensing platforms without any reliance on ground measurements.For the future, we hope our results, when combined with other similar spaced-based data, can provide a global data set for the application of pollutant models and evapotranspiration estimation models.
In this paper, Section 2 contains the introduction of the study area and data; Section 3 contains the method of processing radiosonde data for the validation of the MH, the gap filling of the MODIS profile data and the main retrieval method of the MH; Section 4 is comprised of the results from Section 3 and some validation; and Section 5 and Section 6 contains the discussion and results of the paper.

Study Area and Data
The Heihe river basin (HRB) is a typical arid region of north western China (Figure 1).It is located between 96°42′-102°00′E and 37°41′-42°42′N with a drainage area of 130,000 km 2 , and has a desert-oasis-river landscape structure [13].It is an ideal field laboratory to pursue integrated PBL studies for several reasons.Firstly, although human activity, especially urban heat islands and skyscrapers, can have a large impact on the structure of the PBL, the economy here is relatively backward and there are few metropolises, so any effect is limited.Additionally, comprehensive experiments with radiosonde observations, such as WATER [14] and HiWATER [15] have been conducted here, and radiosondes have been released at Wuxing, Gaoya, Biandukou and Arou.Further, three national radiosonde stations in the region (at Zhangye, Jiuquan, and Ejinaqi) have accumulated standard sounding data for more than 20 years (Figure 1).
Radiosondes are the most common source of data for operational determination of the MH.They are widely distributed spatially, and the data are subjected to quality control on a continual basis.The three national stations Zhangye, Jiuquan and Ejinaqi within the HRB, operated by the China Meteorological Administration, launch radiosondes twice daily at 08 Beijing Time (BJT = UTC/GMT + 8 hours), and 20 BJT.These standard sounding data include regular pressure layers (1000, 920, 850, 700, 500, 400, 300, 250, 200 hPa, etc.) and were obtained through the Natural Science Foundation of China (NSFC).High vertical resolution sounding data for 2008 and 2012 were obtained from the WATER and HiWATER projects, respectively, which collected GPS-located radiosonde sounding observations for the middle and upper stream sections of the HRB.Both sets of sounding data were used for validation.
The MODIS MOD07 L2 product was chosen as our primary data source for deriving the MH, mainly because of its high spatial resolution and the advancement of its data collection process.Although the MODIS is not a sounding instrument, it does have many of the spectral bands found on the High Resolution Infrared Radiation Sounder (HIRS) currently in service on the polar orbiting NOAA TIROS Operational Vertical Sounder (TOVS).Thus, it is possible to generate the vertical profiles of temperature and moisture as well as total column estimates of precipitable water vapor.The MODIS algorithms were adapted from the operational HIRS and GOES algorithms, with adjustments to accommodate the absence of stratospheric sounding spectral bands and to realize the advantage of greatly increased spatial resolution (1 km MODIS vs. 17 km HIRS) with good radiometric signal to noise ratio (better than 0.35 C for typical scene temperatures in all spectral bands) [16].MODIS MOD07 atmosphere profile data include 20 specific pressure layers; each layer contains elevation, air temperature and dewpoint temperature information.The air pressure layers from the surface to upper air are at 1000, 950, 920, 850, 780, 700, 620, 500, 400, 300, 250, 200, 150, 100, 70, 50, 30, 20, 10 and 5 hPa.For most situations, the MH is at an elevation below 600 hPa; therefore, our analysis emphasizes the data from 600 hPa to 1000 hPa.

In Situ Radiosonde Data Processing
The most common method to determine the MH from radiosonde data is based on the Richardson number method, but it is usually used under convective conditions.The Richardson number Rib is defined as [1]: (1) where g is the gravity acceleration, z0 is the altitude of the measurement location above sea level, θ is the potential temperature, and μ(z) and v(z) the wind zonal and meridional components, respectively.The MH is defined as the height where the Richardson number equals the critical Richardson number of 0.21 [17,18].At heights where Rib is higher than the critical value, the atmosphere is considered to be free of turbulence (free troposphere).Unfortunately, the MH estimates based on standard radiosonde data can result in quite high uncertainty [19,20].Particular problems occur under stable atmospheric conditions since no universal relationships seem to exist between the profiles of temperature, humidity or wind.Consequently, other methods must be used in combination [21,22]: (1) The height at which the virtual potential temperature matches the surface value; (2) The level of the maximum vertical potential temperature gradient; (3) The base of an elevated temperature inversion.
In our paper, the Richardson number method is used to define the radiosonde MH under convective conditions and the combination method is used under stable conditions.

MODIS Data Processing
The satellite data suffer from data gaps.Some gaps are internal to the satellite, primarily related to its orbits or to instrument failures, while others are external such as cloud contamination [23].Several solutions have emerged in recent years [24].In the HRB, the presence of cloud cover is the main factor resulting in spatial gaps in the data.In order to provide complete spatial coverage, gap-filling methodology was devised, which has both spatial and temporal components.First, the spatial component investigates a 5 × 5 pixel area.If enough data are available (more than half of the total), the inverse distance weighted average method is used to fill the empty pixel; alternatively the window is expanded to 10 × 10.The expansion process is repeated until a specified maximum window size (usually, 15 × 15 pixels) is reached.Within the maximum window, if there is still too little available data (less than half of the total), then the average value at the same location from the most recent clear day before and after is substituted.
In addition to the primary air temperature and dewpoint temperature profiles, other variables were also derived, including saturation vapor pressure, actual vapor pressure, relative humidity, mixing ratio of water vapor, and virtual potential temperatures, using the conversion equations defined by the Food and Agriculture Organization (FAO) of the United Nations [25].

MH Retrieval
Under clear conditions within the entrainment zone, turbulent boundary layer eddies mix aerosol and water vapor-laden air with cleaner free tropospheric air.This mixture of clear and dirty air together in the entrainment zone prevents aerosols and water vapor mixing with free air.Thus, other researchers have suggested several methods to retrieve the MH based on the assumption that water vapor and aerosols are much more abundant within the PBL than in the free troposphere above it [26].After examining the ability of all the primary and conserved variables (Figure 2) to identify MH, only the mixing ratio of the water vapor (MR) profile was able to provide evidence of the MH in the Heihe river basin: There was an abrupt decrease for MR profile.Accordingly, an automated method based on the level of the minimum vertical gradient of the MR profile was proposed for identifying the MH.
As an example, Figure 3 shows two profiles of vertical MR gradients on 10 August 2012, for Jiuquan and Zhangye.The height of the minimum vertical gradient of the MR is coincident with an abrupt decrease in the MR and can be considered as the top of the mixed layer.Based on this, the height (location) of the MH is retrieved pixel by pixel.

Radiosonde Data
4 July and 8 July 2012 were used for analysis of the in situ radiosonde data from the HIWATER project.Because of different air conditions on these days, different methods were used for each day.Although we lack direct observations of the convection and buoyancy flux, we can infer the stability of air conditions from the structure of the potential temperature profile [1]. 4 July was stable and 8 July was convective.In Figure 4, the vertical dashed line represents Rib = 0.21, and the horizontal dashed line indicates the height of the MH.In the upper figure, there two dashed lines; the lower one is the height of Rib = 0.21, and the upper one is the level of the maximum vertical potential temperature gradient.For 4 July, the height of the stable boundary layer was identified as the level of the maximum vertical potential temperature gradient instead of the height where the Richardson number equaled 0.21, as shown in Figure 4.For 8 July, the height of the convection boundary layer was identified using the Richardson number method, i.e., as the height where the Richardson number equals its critical value of 0.21 (Figure 4, lower panels).On both days, there was an abrupt decrease in the MR at the height of the MH, which justifies our method of using the MR to derive the MH from MODIS profile data (Figure 4).

MODIS Data Processing
As previously discussed, there were missing data in the MODIS atmospheric profiles due mainly to cloud coverage.A spatial-based gap-filling method was used for continuity when the missing data were limited.Where large areas were missing, data from the most recent available clear day were substituted (Figure 5).
The gap filling and its impact were assessed by using a situation where all MODIS data were available.Using 25 March 2012 as an example, certain areas were screened out randomly, and gap filling was performed for these areas.Then, the gap-filled data were compared to the original data on a pixel basis.As shown in Figure 6, the left picture is the original MODIS dewpoint temperature data from a clear sky, the middle one is the MODIS dewpoint temperature data after screening some areas, and the right picture is gap-filled from the middle one.From Figures 6 and 7, it is clear that our gap filling method worked well when the area of missing data was limited and data from the most recent available clear day were used.The coefficient of determination (R 2 ) was 0.82 and the RMSE was 1.04 for the original and gap-filled data; if large missing areas were considered, a poorer relationship was obtained: the coefficient of determination (R 2 ) was 0.31 and the RMSE was 2.13.The moisture structure of the boundary layer has large spatial variations and can be totally different between cloudy and cloud-free conditions and between different days.Therefore, it is not a good idea to fill sections with large gaps.In particular, using the average value from the most recent clear day in the gap-filling process will produce problematic and unrealistic spatial changes of moisture data.Thus, we prefer gap filling using data from areas with few missing values and avoiding areas with extensive missing data.After gap filling, comparisons between the temperature and dewpoint temperature profiles from in situ radiosonde data and those from MODIS MYD07 for the Heihe river basin under clear and cloudy sky conditions showed a strong correlation.In Figure 8, validation of the MODIS data from clear and cloudy sky conditions is indicated in black and red colors, respectively.The coefficient of determination (R 2 ) between the satellite-retrieved moisture profiles and the in situ measured temperature profiles of the black points for the Zhangye, Jiuquan, and Ejinaqi stations during 2012 for each regular pressure layer is 0.93, 0.94 and 0.92, respectively.This indicates that the MODIS atmospheric profile data correlate well with measured moisture profiles from the radiosondes.The RMSE for the dewpoint temperature is approximately 6.5 K and 0.6 g/kg for the MR.Compared with clear sky, the data from gap filling under cloudy sky conditions indicate a weaker correlation; as shown in Figure 8, the R 2 and RMSE for the three stations are 0.70, 0.82, 0.69 and 13.82, 13.60, 14.66, respectively.These errors may partially be a consequence of the difference between the radiosonde sending time (8 and 20 BJT) and the satellite overpass time (about 13-15 BJT), arising mainly from the uncertainty originating from the gap filling method.

MH Validation
Most of the optical-thermal satellites, including the two MODIS satellites (Terra and Aqua), have overpass times ranging from 9:00 to 15:00 local time.The greatest depth of the PBL usually occurs in the afternoon between 13:00 and 17:00 local time.This most closely corresponds to Aqua's overpass time (13:00-15:00); therefore, MODIS MOD07 data from Aqua were used.Table 1 compares MH estimates from in situ radiosonde data to those from the MODIS data on selected days in 2008 and 2012 when the radiosonde balloon release time was close to the satellite overpass time.The data in Table 1 indicate that the MH calculated from MODIS is largely coincident with that from in situ radiosondes.The absolute error is below 500 meters, which is the resolution of the MODIS moisture profile at its base.However, there is a significant error for 4 July.As discussed previously and demonstrated in Figure 4, the potential temperature profile for 4 July suggests stable atmospheric conditions.For our method, determination of the MH from the MODIS profile data under stable conditions can result in substantial error.This appears to be the case for 4 July.Unfortunately, the lack of in situ data under stable air conditions, which would enable additional comparisons, inhibits further verification of our speculation.
The assumption that stable air conditions may cause substantial errors needs further examination.The effect of this error is apparent from Figure 9.By omitting the data for 4 July, the R 2 between the two estimations increases from 0.59 to 0.82, the RMSE decreases from 550 m to 370 m and relative error changes from 44% to 28%.If there is a regular diurnal cycle of the MH in the Heihe river basin, the residual layer of boundary layer will also be moist and it is likely that the top of the residual layer instead the top of the actual ABL will be found.Therefore, from Table 1 we can find that most of the daily MH values from MODIS data are larger than those from radiosonde data.

MH Spatial and Temporal Patterns
After data pre-processing and model calculation, the peak MH for each day was calculated for selected days in 2008 and 2012.The results from clear days in winter and summer are mapped in Figure 10 for two example days.The mixing layer height is primarily driven by variation in the vertical flux of energy and momentum caused by the balance between the heating of the earth's surface by incoming solar radiation and radiative cooling of the surface [27][28][29].Accordingly, as Figure 10 shows, the MH is lower in the winter and relatively higher in the summer, demonstrating that net radiation is the main factor driving the development of the PBL.Especially in the western portion of the downstream area, where desert land use dominates, the growth amount of the MH extends higher than that in the middle and upper stream regions where oasis landscapes occur.Forest and grass vegetation may have a buffer for temperature variation; therefore, the oases area results in more consistent land surface temperatures, which keeps sensitive heat more stable and lower than in the desert area, eventually leading to a similarly stable and lower change in the MH.For every single day, the MH showed a stronger correlation with the actual land conditions, i.e., higher for desert than grass/forest during both winter and summer.
The annual variation in the MH was analyzed at the Gaoya station.Figure 11 shows the daily variation and monthly average MH at this station in 2012.There is a clear annual cycle with the MH growing from about 2500 meters in January to 3000-3500 meters in the summer and autumn, and descending to about 2500 meters again at the end of the year.

Discussion
Although the average trend of the MH is reasonable at the Gaoya station, there is a large variation from day to day.Several reasons may account for this daily variation: the MH is strongly related to net radiation, which has a large daily variation; there is a difference between the overpass times of the satellite on each day, and the diurnal variation of the MH can be as much as several kilometers; the presence of cloud cover, and the consequent use of gap-filling to estimate missing data may generate large errors in the moisture profile; under stable air conditions in particular, the automated MH estimation method may not work well, resulting in erroneous MH estimates; and the elevation of each layer derived from MODIS is not always accurate, especially on cloudy days when the estimated MH can be even lower than the surface elevation.All of the above require consideration and further research if the results generated by our method are to be more broadly utilized.
Cloud is a major impediment for the accuracy of our method.With cloudy skies, underlying surface air conditions and other atmospheric variables are uncertain or unknown.A gap-filling method was proposed for filling these data gaps, but this method itself can generate errors.The vertical resolution of the atmospheric profile data is another impediment to our approach.The vertical resolution of the MODIS profile data is about 500 m below 780 mb, and about 1000 m above 700 mb.For most days, the MH is below 700, and the RMSE of our method was close to the vertical resolution of the MODIS profile data.Although our method can detect the correct layer of the MH, there is a possible several hundred meter bias, corresponding to the vertical resolution of the atmospheric profile data at its base.As discussed previously and demonstrated in Figure 4, for our method, determination of the MH from MODIS profile data under stable conditions can result in substantial error.This appears to be the case for 4 July.The lack of in situ data under stable air conditions in the Heihe river basin inhibits further verification of our speculation.The radiosonde data from citation [30], as Figure 12 showed, indicate the reason for the large difference in the MH from the radiosonde and the MODIS profile data.The MH method applied to radiosonde data in this paper is based on the potential temperature profile; for MODIS data, the MH is derived from the mixing ratio profile; under convective conditions (Figure 12a), the MH diagnosed from the radiosonde potential temperature profile is consistent with the height of minimum vertical gradient of MR.Consequently, our method for MODIS data based on the MR profile works well under convective conditions; it is different for the stable condition (Figure 12b), where the MH diagnosed from the radiosonde potential temperature profile is below 100 m, but the height of the abrupt decrease in the MR is about 1500 m, which corresponds to the top of the residual layer of previous days, similar to 4 July, i.e., 627 m MH for radiosonde and 2156 m for MODIS.
There are caveats associated with using the water vapor profile to diagnose the mixing height: (1) If there is a regular diurnal cycle, especially under stable conditions, the residual layer will also be moist and it is likely that the top of the residual layer is identified instead of the top of the actual ABL; (2) Moist layers in the lower atmosphere do not necessarily indicate local mixing, they could also be advected.Considering the mountainous terrain in the south of the study area, and the fact that it is not as arid as the northern parts of the Heihe river basin, this is a real possibility in the case of southerly winds.This should also receive attention when the method of deriving the MH in the paper applied to other areas, especially where there is a regular diurnal MH cycle or near the coast where moist layers in the lower atmosphere may not indicate the MH, but rather sea-land advection layer during the daytime.
Despite weaknesses, initial estimates of the PBL depth using the atmospheric moisture profile from space-based satellite seem reasonable although more validation is needed.With enough computer resources, this analysis could be extended to provide a global PBL data set using MODIS MOD07 data or relatively high accuracy but low spatial resolution data such as the AIRS product.Our result is also a potential data source for input into global pollutant models.New sources of high accuracy atmospheric moisture profile data with higher vertical resolutions, from future satellites, for example, will promote our method since observation of the atmospheric boundary layer is difficult and expensive using traditional, ground-based methods.If high temporal resolution data from geostationary meteorology satellite show a similar ability as MODIS data for MH calculation, the diurnal variation of MH could also be derived.

Conclusion
In this paper, the mixing ratio was selected from primary and conserved variables to identify the MH.The data suggest that the mixing ratio profile can provide evidence of the MH.The MH occurs when there is an abrupt decrease in the mixing ratio (MR).An automated method was developed to obtain the mixing height of the atmospheric boundary layer from MODIS moisture profile data.This study has demonstrated the feasibility of space-based remote sensing of the mixing layer height (MH) without the aid of ground measurements.The methodology for estimating the MH works well on clear-sky days (no clouds) with convective air conditions.However, cloudy skies and stable air conditions may result in large estimation errors.Consequently, cloud and vertical resolution of the atmospheric profile data are two major impediments to the accuracy of our method.Nevertheless, the preliminary results are encouraging both for spatial and temporal variation of the MH.In addition, the MH can be estimated more accurately if higher vertical-resolution data are available.

Figure 1 .
Figure 1.Heihe river basin location and landuse; location of radiosonde stations.

Figure 2 .
Figure 2. Profile of primary and conserved variables retrieved from MODIS atmospheric profile data for Jiuquan and Zhangye radiosonde stations on 10 August 2012.

Figure 3 .
Figure 3. Profile of the MR and vertical gradient of the MR (in red) retrieved from MODIS atmospheric profile data for Jiuquan and Zhangye radiosonde stations on 10 August 2012.

Figure 4 .
Figure 4. Profile of potential temperature, MR and Richardson number calculated from radiosonde data collected at Gaoya station on 4 July and Wuxing station on 8 July, 2012.

Figure 5 .
Figure 5.The moisture data for the Heihe river basin at 700 hPa before and after application of the gap-filling process for 2 January 2012.

Figure 6 .
Figure 6.The validation of moisture data for the Heihe river basin at 700 hPa: original data, screening area and data after gap filling for 25 March 2012.

Figure 7 .
Figure 7. MODIS original moisture data vs. those from the gap-filling method for 25 March 2012 at 700 hPa pressure within the Heihe River Basin.

Figure 8 .
Figure 8. Validation between MODIS moisture data and radiosonde data for 2012 at Zhangye, Jiuquan and Ejinaqi stations for regular pressure layers.

Figure 9 .
Figure 9. R 2 and RMSE for the MH values retrieved using our method and radiosonde data for all days (a) and for convective days (b) at the Heihe river basin in 2008 and 2012.

Figure 10 .
Figure 10.Spatial variation of the MH derived from the MODIS MOD07 product for the clear days 8 January and 6 July 2012.

Figure 11 .
Figure 11.Daily and monthly variation in the MH calculated from the MODIS profile data from the Gaoya station in 2012.

Figure 12 .
Figure 12.Case of water vapor mixing ratio and potential temperature profile under convective (a) and stable air conditions (b).The horizontal line depicts the position of the determined planetary boundary layer height (PBLH) [30].

Table 1 .
Comparison between the MH calculated from radiosonde data and that from MODIS data for 13 days in 2008 and 2012 for the Heihe river basin.