Diurnal and seasonal variation of planetary boundary layer height over East Asia and its climatic change as seen in the ERA-5 reanalysis data

The diurnal/seasonal structure of the boundary layer height (BLH) is investigated over East Asia by using the hourly synoptic monthly ERA5 reanalysis variables from 1979 to 2019. Sensible heat flux (SHF) is the major factor in the temporal and spatial variation of the BLH. Although BLH, in general, is positively correlated with SHF throughout the year, BLH-SHF relationship varies significantly based on the surface type, latitude and time of the year. Analysis also reveals that stability is an important parameter controlling the diurnal maximum BLH. The growth of BLH is strongly limited by the presence of a stable layer. On the other hand, BLH increases abruptly in the presence of a weakly stratified residual layer. In addition, regional warming tends to increase the BLH in the mid- to high-latitude continental area. In the low-latitude continental area, the sign of anomalous SHF varies seasonally and regionally. Stability plays only a minor role in the BLH change except over the Tibetan Plateau, where the increased stability at the top of boundary layer due to warming reduces BLH rather significantly.


Introduction
The planetary boundary layer (PBL) is the lower portion of the atmosphere in which the flow field is strongly influenced by direct interaction with the Earth's surface. The portion above the PBL is called the free atmosphere, and a clear difference between the two layers is the timescale of response to surface forcing [1][2][3]. Turbulence generated from the surface, by definition, is trapped within the boundary layer due to the stably stratified inversion layer, also called the entrainment layer, at the top of the PBL. Thus, the diurnal variation that comes from the daily insolation occurs primarily below the inversion height. In addition, the exchange of heat, momentum, water, and chemical constituents between the surface and the free atmosphere occurs in the PBL, and specific features of this exchange are highly dependent on the turbulent motions in the PBL [4]. In this respect, the dispersion of chemical contaminants produced by human activities is closely associated with the physics of the PBL.
This study focuses on analyzing the temporal and spatial variability of the boundary layer height (BLH) over East Asia. Height of the convective boundary layer (CBL) reflects the intensity of turbulent convection in the boundary layer. In fact, turbulent mixing is not readily obtainable from observations, but it is often parameterized in terms of the BLH [5][6][7], which can be obtained from the vertical profiles of virtual potential temperature, water vapor mixing ratio, relative humidity, and/or wind speed. In particular, Deardorff [5] used the vertical buoyancy flux and the BLH to calculate the velocity scale called the convective velocity, which characterizes the turbulent motions by free convection. The convective velocity and the BLH itself are the main scaling parameters in some CBL similarity hypothesis.
The BLH is also an important variable in the study of air quality [8][9][10][11]. A taller PBL implies more vigorous mixing of chemical constituents in the boundary layer, so that they disperse to higher elevations and surface concentrations become lower. In this respect, the volume of the PBL can be thought of as the holding capacity of the chemical constituents emitted from the surface. Likewise, Davy and Esau [12] considered the PBL as a reservoir for excessive heat generated at the surface and explained the shallow PBL as a major reason for the amplified warming in the Arctic. As exemplified here, the BLH is an important physical factor in elucidating diverse atmospheric phenomena.
In order to investigate the BLH from a climatological perspective, it should be ensured if the algorithm of determining the BLH is reasonable in terms of the maximum and minimum heights of the PBL as well as their diurnal and seasonal variations. Several methods of defining the BLH have been introduced and compared (e.g., [13][14][15], Sawyer and Li 2013; [16]). For instance, Seidel et al. [15] computed the BLH by using the profiles of different variables such as temperature, potential temperature, virtual potential temperature, relative humidity, specific humidity, and refractivity. They addressed that uncertainty is inherent in defining the BLH in a rigorous manner regardless of the season, surface condition, and thermal stratification; each method has its own merits and demerits depending on specific situations. Inter-comparison of different methods, however, is beyond the scope and will not be dealt with in this study.
In the present study, the ERA5 reanalysis dataset (https:// www. ecmwf. int/ en/ forecasts/ datasets/reanalysis-datasets/era5) is employed. The ERA5 reanalysis model defines the BLH as the minimum height for the bulk Richardson number reaching the value of 0.25. This algorithm was proposed by Vogelezang and Holtslag [17] and was used in recent studies [18][19][20][21]. According to the "Integrated Forecast System (IFS) documentation" of European Center for Medium-Range Weather Forecasts (ECMWF), this method is little influenced by the dataset employed (such as reanalysis, radiosonde sounding, and lidar observations) and returns appropriate values both for convective and stable boundary layers.
Medeiros et al. [22], Seidel et al. [15], Liu and Liang [23], Seidel et al. [20], Engeln and Teixeira [13], and Guo et al. [18] described the mean state of the seasonally repeating patterns of the BLH without referring to any climatological shift. Examples of BLH studies from a climatological perspective include Yang et al. [24] and Pal and Haeffelin [25], who examined the long-term and inter-annual variability of the BLH at specific sites. Zhang et al. [21], Guo et al. [19], and Zhang and Li [26] investigated the observed trends of the BLH over Europe, China, and Japan, respectively. They suggested several variables that are positively or negatively correlated with the climatological shift of the BLH by comparing their trends with that of the BLH. The diurnal variation of the climatological shift, however, was not examined.
Variation of the BLH exhibits two primary periodicities: diurnal and seasonal cycles. It should be noted that change in the PBL depth depends significantly on the temporal scales (diurnal versus seasonal) as well as the surface types. Furthermore, BLH is spatially inhomogeneous since it is strongly influenced by the surface types and flux. Therefore, it is imperative to understand how the climatology and the impact of climate change differ depending on the time scales and the surface types.
While there are many studies based on observational data across the analysis domain providing important insights into the variability of the BLH, they are typically limited in terms of spatial and/or temporal scope. Many studies are based on observational data at a few sampling stations without sufficient spatial coverage. Some studies use specific seasons or annual averages in addressing the inter-annual variability of the BLH and the impact of climate changes on the BLH. Not many studies investigate how climate changes affect the diurnal variation of the BLH in terms of maximum and minimum BLH. It is important to understand the variability of the BLH in the context of the diurnal and seasonal variation of the BLH.
In this study, therefore, detailed and comprehensive analysis of the BLH is conducted over East Asia together with the impact of regional warming by applying the cyclostationary empirical orthogonal function (CSEOF) technique (see [27][28][29] to the ERA5 reanalysis datasets). The main goals of the present study are (1) to investigate the variability of the BLH together with the major controlling factors over the entire domain and (2) to delineate the impact of regional warming on the maximum and minimum BLH on the seasonal basis. This study aims to depict the full spatiotemporal structure of BLH variability throughout the year, thereby complementing the observational studies that generally lack spatial and/or temporal details. The use of the ERA5 reanalysis product, instead of observational datasets, is an important caveat since uncertainty is inherent in the reanalysis product. Section 2 describes the data employed in this study and the methods used to analyze the data. In Sect. 3, the fundamental features of the diurnal and seasonal evolution of the BLH are investigated. The impact of regional warming on the BLH is addressed in Sect. 4, followed by summary in Sect. 5 and concluding remarks in Sect. 6.

Data
Data used in this study include monthly averaged variables at every hour of the day (also called synoptic monthly variables) for the period of 1979-2019 (41 years). Variables are analyzed over the domain (70.5°-160.5°E × 0°-60°N) at a 1.5° × 1.5° resolution (60 × 40 grid points). Thus, there are total of 11,808 (= 24 h × 12 months × 41 years) time steps. Based on the findings in earlier studies (e.g., Avissar and Schmidt 1988; [30][31][32][33][34]), such surface variables as sensible heat flux, latent heat flux, sea surface temperature, surface air temperature, and frictional velocity are analyzed to investigate the effect of surface forcing on the BLH. In addition, pressure (1000-200 hPa; 23 vertical levels) variables including air temperatures, geopotential height, and vertical (pressure) velocities are used to calculate the stability of the atmospheric column.

Method of analysis
In the present study, CSEOF analysis [27][28][29] is conducted in order to extract the diurnal/seasonal cycle of the BLH and the regional warming mode along with the associated physical variables. Unlike EOF analysis, each CSEOF loading vector describes a temporal evolution of a physically distinct mechanism (see supplementary Movies S1 and S2). Further, regression analysis in CSEOF space allows physical consistency among the evolutions of different variables. Thus, it is ideal for describing the diurnal/seasonal variation of the BLH and the impact of regional warming on it together with the responsible physical factors.
Each variable is subject to CSEOF analysis, which decomposes data as where B n (r, t) are called the cyclostationary loading vectors (CSLVs; see Movie S1 in the supplementary information), T n (t) are the principal component (PC) time series (see Figs. 1e and 8e), r is location, and t is time. Owing to the cyclostationarity assumption [29,35], CSLVs are periodic, i.e., where d is called the nested period, which is 288 (= 24 h × 12 months) in the present study. As in EOF analysis, CSLVs are mutually orthogonal and PC time series are mutually uncorrelated, i.e., where n are eigenvalues (variance of T n (t) ), and nm is the Kronecker delta. Thus, CSLVs are mutually uncorrelated and orthogonal evolutions extracted from the given variable.
In a typical setting, many variables are analyzed in order to understand the physical mechanism behind each CSEOF mode. For example, another variable is decomposed as (2) B n (r, t) = B n (r, t + d), P(r, t) = ∑ n C n (r, t)P n (t), where C n (r, t) are the CSLVs and P n (t) are the PC time series of variable P(r, t) . The two sets of decomposition in (1) and (5) are not physically consistent in that T n (t) ≠ P n (t) in general. In order to make the two sets physically consistent, regression analysis is conducted in CSEOF space between the so-called target variable T (r, t) and the predictor variable P(r, t) . It is a two-step process: where K ( = 20 in this study) is the number of CSEOF modes used for regression. Then, the two variables can be written, aside from a small regression error, as This procedure can be repeated for as many predictor variables as needed. Then, the entire dataset can be written as where the terms in curly braces represent CSLVs for different variables. The "regressed" CSLVs are physically consistent with B n (r, t) in that they all share the same PC time series T n (t) [see supplementary Movie S2 for an example].

Seasonal patterns of maximum and minimum boundary layer height
The first CSEOF mode represents the diurnal/seasonal cycle of the BLH. Figure 1 shows the daily maximum BLH (shade) and the minimum BLH (contour) derived from the loading vector of the first CSEOF mode (see Movie. S1 for the full loading vector). For the sake of brevity, loading vector is shown only in four months (January, April, July, and October) representing the four seasons. It should be noted that these patterns are nearly identical with the seasonal average patterns (figure not shown). In January, both the maximum and the minimum BLHs over the continent are generally lower than those in the other months; typical minimum values are found in January or December (Fig. 2b, d). Over the ocean, on the contrary, both the maximum and the minimum BLHs are Step higher than those in the other months; maximum values are typically found in January or February (Fig. 2a, c; see also [36,37]). In April, maximum BLH increases significantly over the continent reaching the annual maximum values in May or June (Fig. 2a). Over the ocean both the maximum and the minimum BLHs decrease slightly. In July, the maximum and the minimum BLHs slightly decrease over the continent compared with the spring months. Over the ocean, both the daily maximum and the minimum BLHs exhibit the annual minimum values in summer (Fig. 2b, d). In October, maximum BLH begins to decrease over the continent, and the minimum BLH begins to increase over the ocean. For each month, the PC time series of the seasonal cycle exhibits small inter-annual fluctuations of about ± 5% around its mean (Fig. 1e). The amplitude of the seasonal cycle shows a strong monthly variation with a maximum in May and a minimum in November. This peculiar behavior indicates that the inter-annual modulation effect on the amplitude of the seasonal cycle differs from one month to another, with a stronger modulation in May and a weaker one in November.
The structure of the seasonal evolution of the BLH over the domain is generally similar to that over Europe with a seasonal maximum around July [25]. The seasonal patterns of the daily maximum BLH over China are reasonably similar to those estimated from the CALIPSO (Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observation, [37]). As noted in Liu et al. [37], however, some differences both in terms of the amplitude and the phase of the seasonal BLH are seen between the CALIPSO and the ECMWF datasets. The annual and seasonal patterns of the BLH at specific hours (see Fig. S1 and the discussion therein) are also similar to the estimates based on the sounding data in China Research Article [18,38]. Contrast in the seasonal evolution patterns of the maximum BLH over the ocean and the continent is also similar to that in Chan and Wood [39] and Medeiros et al. [22]. Compared with the eastern North Pacific [40], western North Pacific exhibits a stronger seasonality in the variation of the BLH (Fig. 1). Figure 3 shows the pattern of frictional velocity together with the minimum BLH. Frictional velocity is given by which represents the horizontal eddy momentum carried by vertical eddies [3,41]. Eddy momentum flux is much higher over the ocean than over the continent, which implies that the mechanical production of eddy kinetic energy tends to be stronger over the ocean than over the continent given similar magnitudes of the vertical wind shear. There is only small seasonal variation of the frictional velocity over the continent, while the frictional velocity is much higher in cold months than in warm months over the ocean. This is primarily due to the increased lower tropospheric stability and wind [42] and is partly due to the increased vertical velocity associated with the seasonal variation of the temperature difference between the sea surface and surface air (see Fig. S2). As can be seen in the figure, frictional velocity is a good indication of the

Physical factors governing maximum and minimum boundary layer height
minimum BLH both over the continent and over the ocean, and pattern correlation is 0.64 over the continent and is 0.66 over the ocean. Figure 4 shows the pattern of the SHF over land and that of the LHF over the ocean together with the maximum BLH. Energy budget in the PBL is described as [3,43,44] where ⟨ ⟩ is the mean potential temperature in the PBL, ′ w ′ is the heat flux at the surface (S) and the top (T) of the PBL, and ΔF z is the longwave radiative forcing. An increase in the surface heat flux, then, increases the maximum BLH assuming that the heat flux at the PBL top is small compared to the surface heat flux. Radiative forcing is neglected in the present study. As can be seen, the pattern of SHF (LHF) matches reasonably with that of maximum BLH over land (ocean). Note that latent heat flux, when it is converted into heat upon condensation, also serves the same role as sensible heat flux. It appears that the maximum BLH increases as SHF (LHF) increases over the continent (ocean) [22,25]. There is a substantial difference between the continent and the ocean in terms of the relative magnitude of the maximum BLH to that of the heat flux. It should be noted that SHF (LHF) is mainly responsible for the variation of the maximum BLH over land (ocean), and pattern correlation is 0.75 (0.56) over land (ocean). Fig. 3 The patterns of frictional velocity (shade; cm s -1 ) together with the minimum boundary layer height (blue contours < 1.0  In January, stability quickly rises from the surface over the continent, limiting turbulent mixing close to the surface in the morning (see January 22 UTC in Fig. 5a, b, c). During the daytime, BLH increases by ~ 1-1.5 km up to the base of the entrainment layer with stability greater than 1.5 K per 25 hPa (see January 07-09 UTC in Fig. 5a, b, c). Over the ocean, BLH is fairly high throughout the day. As already discussed in conjunction with Fig. 4, LHF is highly correlated in pattern with the maximum BLH over the ocean. The release of a large amount of oceanic LHF is manifested as weak stability in the low troposphere (see January in Fig. 5). As a result, the maximum BLH over the ocean is relatively high.
In April, the morning BLH is seen near the level of maximum stability over the continent (see   Fig. 5a, b, c). Above the layer of maximum stability, a thick layer of weak stability is found. This relatively well-mixed layer, which may be a leftover from the daytime PBL, is known as the residual layer [3]. During the daytime, stability decreases near the surface so that turbulent mixing penetrates the layer of weak stability and BLH jumps to a much higher elevation (see April 08-10 UTC in Fig. 5). Over the ocean, BLH is relatively low as the release of LHF/SHF is reduced.
In July, situation is similar to that in April with a large diurnal excursion of the BLH over the continent and a shallow boundary layer over the ocean (see July in Fig. 5). It is typically between April and July when the BLH reaches a maximum in mid-to high-latitude land areas (Fig. 2). As can be seen in Fig. 6, the maximum BLH is observed earlier on the eastern side of the continent with 06-07 UTC (14-15 Beijing time) being the time of maximum near Beijing (116°E).
In October, a morning layer of high stability near the surface results in a shallow PBL over the continent (see right panels of October in Fig. 5). During the daytime, stability decreases appreciably, resulting in an increased BLH (left panels of October in Fig. 5). Over the ocean, BLH begins to increase because of the increased LHF. The vertical structure of stability (shade; S p × 25 K hPa −1 ) and its change due to regional warming (blue < 0 < red contours; S p × 2500 K hPa −1 ), and the boundary layer height (thick black line) at the approximate time of the maximum height (left column) and the minimum height (right column): a 40°N  It appears that the diurnal cycle of the BLH is strongly controlled by the SHF over the continent. Figure 6 shows the diurnal structure of the SHF, BLH, and the stability at the top of the BLH (see also supplementary Figs. S3, S5, and S7). During the cold months (January and October), BLH over the continent increases as SHF increases during the daytime. It appears that stability poses a strong restriction as to the maximum height of BLH; BLH rises until it encounters a relative stable layer (red dots). During the warm months, when the layer of relatively low stability is fairly thick, BLH increases significantly as the daytime SHF increases. In particular, BLH jumps to ~ 2.5 km above the surface in April-July along 40°N (see also [23,25], Chu et al. 2019). The maximum rise of the BLH seems to be controlled primarily by the SHF but is also limited by the stability. BLH rises typically until it encounters a stable layer, although the top of the PBL tends to reach the elevation of higher stability as the magnitude of SHF increases. In the late afternoon, BLH begins to decrease as SHF diminishes [45]. A sharp drop in the BLH is observed around the midnight and stays low until the next morning (see also Fig. 4 in Hegarty et al. [46]). Over the Indian subcontinent and the surrounding oceans, the diurnal variation of the BLH exhibits a structure similar to that in Sathyanadh et al. [47]. Over the ocean, the BLH exhibits a much flatter diurnal cycle [23,40].
A relationship between the SHF, BLH, and stability is examined in Fig. 7. In general, SHF and BLH are positively correlated. The bins with the red dots (higher occurrence frequencies) near the origin largely represent the events after the sunset until the sunrise. As should be expected, there is a strong seasonal variation in the maximum BLH and maximum SHF. The spring season exhibits the largest range of the BLH and SHF, whereas the fall season exhibits the smallest range. In general, the upper limb (with higher ratios of BLH to SHF) of the spread is associated with lower stability than the lower limb with lower ratios of BLH to SHF. On the other hand, the far end of the spread (strong SHF and/or large BLH) exhibits a wider range of stability from 1 to 2. The bins with tall boundary layers in DJF and MAM consist largely of the cases from the Tibetan Plateau (see Fig. S9). This indicates that the BLH responds more strongly to SHF over the Tibetan Plateau (see also Fig. 4). It seems that a thick residual layer is formed over the top of the Tibetan Plateau in spring and winter. As SHF increases during the day, stability near the top of the Tibetan Plateau decreases, and BLH jumps to the top of the residual layer. The anomalously high PBL over the Tibetan Plateau region in winter and spring is also discussed in Chen et al. [8]. They asserted that the upper-level PV structures and the jet position as well as the surface condition and heat flux are important factors for the BLH because of the high orography of the Tibetan Plateau. It should be noted that the BLH climatology over the Tibetan Plateau differs significantly among different reanalysis datasets; it is cautioned that physical relationship may vary depending on a particular reanalysis dataset employed because of the model biases and errors in these extreme conditions. Over the ocean, the role of the stability is unclear (see Fig. S10), although there is a general linear relationship between the BLH and LHF (and SHF).
4 The effect of regional warming Figure 8 shows the second CSEOF mode, which depicts the effect of regional warming on the BLH. The corresponding PC time series shows a positive trend (Fig. 8e), which clearly indicates that this mode is associated with regional warming (see also the associated surface air temperature in Fig. S11). The net effect of regional warming during the 1979-2019 period is determined to be approximately 1.6 (change in the PC amplitude between 1979 and 2019) times the loading patterns in Fig. 8 (see also Eq. (1)). The PC time series exhibits strong seasonal variations, suggesting that the effect of warming differs from one month to another. On top of the trend, the inter-annual fluctuations in the PC time series indicate that warming also arises due to natural sources. While a linear fit is used in Fig. 8e, the rate of increase seems to vary on decadal time scales. In particular, warming trend seems to be much lower during 2004-2019 compared with that during 1989-2003 (see also Fig. S12 and discussion therein). This result seems to be roughly consistent with the finding in 19, see their Fig. 2.
In January, maximum BLH has decreased over the midlatitude region of the Asian inland, while it has increased over the ocean. The minimum BLH has generally increased over the whole domain except over the mid-latitude inland China. In April, the maximum BLH has increased by 50 m or more over the continent except over the Tibetan Plateau and India. The minimum BLH has decreased over the northern part and increased over the southern part of the domain (mainly over the ocean). In July, fairly significant increases (> 100 m) in the maximum BLH are seen over much of East Asia with a reduction in the minimum BLH to the north of ~ 20°N. In October, a slight increase in the maximum BLH is seen over the continent with a reduction in the minimum BLH in the mid-and low-latitude land areas. As in earlier studies ( [21] in Europe [19], in China), BLH averaged over the entire domain shows a positive trend in association with warming. On the other hand, trend is not uniform spatially with even negative values around the Indian subcontinent particularly in winter and spring (see also [48]. Fig. 6 The diurnal structure of boundary layer height (contours; 0.2 km interval up to 1.0 km and 0.5 km interval thereafter) together with sensible heat flux (shade; W m −2 ) over the continent, latent heat flux (shade; W m −2 ) over the ocean, and stability at the top of the boundary layer (color dots; S p × 25 K hPa −1 ) for the first CSEOF mode (diurnal/ seasonal cycle). The upper color bar represents the shading levels for heat flux and the lower color bar the dot color levels for stability: a 40°N,  As shown in Sect. 3, SHF plays an important role in the growth of the boundary layer. Figure 9 shows surface heat flux anomalies in comparison with the maximum BLH anomalies for the regional warming mode. Change in the maximum BLH due to warming is generally consistent with the change in SHF (LHF) over land (ocean); pattern correlation is 0.74 (0.56) over land (ocean). The maximum BLH tends to increase as SHF increases and the opposite is true. SHF has increased in the northern part of the inland Asia, but decreased over India, Indochina, and Tibet. Positive correlation between the SHF and the maximum BLH is more robust in the warm season (May, July, September) than in the cold season (January, March, November).
Over the ocean, LHF has increased due to warming except over the Northern Pacific (see Fig. 9). In March and November, significantly increased heat flux is seen near the Kuroshio Current and its extension in the latitude band of 20°-35°N. This leads to increased maximum BLH in March and November (see also Fig. S13b). On the other hand, decreased BLH is seen in March and November north of ~ 40°N (see Fig. S13a, c).
It should be noted that change in the maximum BLH is not quite proportional to that of the SHF (see Fig. 9). As we have already discussed in conjunction with Fig. 5, static stability may be an important parameter controlling the PBL. The colored contours in Fig. 5 exhibit the changes in static stability under regional warming. It should be noted that there have been noticeable changes in the vertical structure of the atmosphere as well as the surface condition due to warming. Since turbulence generated from the surface is trapped in the lower troposphere below a strongly stratified layer, change in the static stability may alter the height that turbulent eddies can rise. Guo et al. 120°E 1 50°E [19] also investigated the stability change in comparison with the inter-annual trend of the BLH. They used the lower tropospheric stability defined as the potential temperature difference between the 700-hPa level and the surface. In the present study, static stability ( S p , Eq. 12) is calculated at the top of boundary layer where turbulence that originated from below encounters strong resistance. Figure 10 displays the anomalous stability at the top of the climatological PBL in the SHF-BLH space for the different seasons. Again, there is, in general, a positive correlation between SHF and BLH changes (see also Fig. S13). SHF has generally increased due to regional warming as indicated by a larger number of dots on the positive side of ΔSHF except in DJF. In the winter and spring seasons, there is a large number of bins scattered in the third quadrant, many of which represent infrequent events. The pattern of scattering indicates that the degree of BLH reduction to decreased SHF varies widely in winter and spring. As can be seen, stability at the top of the climatological BLH has generally decreased throughout the year. This suggests that change in SHF is the primary factor for the BLH change in conjunction with regional warming, and change in stability generally does not play a major role. There are, however, some notable exceptions. For example, many bins in the second quadrant in winter and spring exhibit significantly lowered stability. These bins represent the high-latitude regions of East Asia. It appears that the BLH has increased even though the SHF has decreased, because stability has decreased significantly at the top of the climatological boundary layer.
It is known that turbulence can develop even in stably stratified upper troposphere near the jet stream because of the strong vertical wind shear [49,50]. Thus, the jet stream can contribute substantially to the growth of the boundary layer over the Tibetan Plateau. In response to regional warming, the BLH over the Tibetan Plateau decreased significantly in winter and spring due to the decreased SHF and increased stability at the top of PBL (see Fig. S14). Judging from the higher frequency of bins along the y-axis, it appears that stability is the main reason for the reduction of BLH in these seasons (see also Movie S1). The dramatic change of stability over this region could be attributed to the change in the position of the jet stream in response to warming. A shift of the jet stream may be explained by the meridionally asymmetric warming pattern (see Fig. S11) in the context of the thermal wind balance.

Summary
In this study, the daily and seasonal cycles of the BLH were investigated in East Asia together with the impact of regional warming by using the hourly synoptic monthly ERA5 reanalysis variables. The diurnal structure of the BLH changes appreciably with the season in terms of the minimum/maximum height of the boundary layer and exhibits strong season-dependent spatial variation related to the surface type, elevation, and latitude (Fig. 1).
Over the continent, SHF is a major factor determining the maximum BLH throughout the day (Figs. 4 and  6). There is a strong positive correlation between the two variables, although the slope of the linear relationship varies significantly from one month to another and from one location to another (Fig. 7); this seems to be associated with several different factors including the surface type and the vertical structure of stability of the atmospheric column. The maximum BLH is achieved ~ 2-3 h after the maximum SHF in the warm months (Fig. 6). The static stability of the atmospheric column, to a lesser extent, also exerts influence on the maximum BLH, particularly in the cold months (Fig. 5). A strong stability above the cold continental surface in winter often limits the growth of the BLH. In the absence of substantial surface heat flux, the minimum BLH is largely controlled by mechanical production of eddies measured by the frictional velocity (Fig. 3).
The seasonal maximum BLH is observed typically between April and July north of ~ 40°N, and August-September and March-April to the south of ~ 40°N (Fig. 2). The Tibetan Plateau is an exception with the maximum BLH being observed in February. The seasonal minimum BLH is seen in December over much of the continent. The low-latitude regions including India, the Indochina Peninsula, and the Tibetan Plateau, on the other hand, exhibit the seasonal minimum BLH in the warm months. It should be mentioned that the residual layer is an important seasonal feature controlling the maximum BLH. In the presence of a weakly stratified residual layer, turbulent motion easily reaches the top of the residual layer with the aid of surface forcing. The seasonal patterns of the minimum BLH show complex patterns with a strong latitudinal and regional dependency ( Fig. 2c  and d).
Over the ocean, LHF rather than SHF seems more strongly related to the maximum BLH. While the diurnal variation of the BLH over the ocean is not significant, seasonal variation is relatively strong. In particular, the maximum BLH varies from ~ 100 m in July to ~ 1.5 km in January over the northern North Pacific (Fig. 4). Over the low-latitude ocean, seasonal variation of the BLH is not so dramatic. For the regional warming mode, changes in SHF and BLH, in general, display positive correlations with changes in BLH (Fig. 10). SHF has generally increased over much of the midand high-latitude continental domain due to regional warming except perhaps in spring (Fig. 9). Static stability at the top of the boundary layer has also decreased, in general, due to warming (Fig. S13a and c). This implies that turbulence can reach a higher elevation, thereby increasing the BLH under regional warming over the mid-and high-latitude continental areas. Over the low-latitude continental area, on the other hand, SHF has decreased (increased) on the western (eastern) side of the domain throughout the year. The daytime static stability is generally consistent in sign with the SHF. This implies that BLH reduction (augmentation) is observed on the western (eastern) side of the continent (Fig. 9). A clear exception is July, when the BLH has increased over the western part of the low-latitude continental area even though SHF has decreased (Fig. S13b). The effect of regional warming is not quite steady, since temperature fluctuations can come by several different sources such as the regional monsoons, sea ice reduction in the Arctic, and El Niño in addition to a steady increase in the concentration of greenhouse gases (Fig. 8e). As a result, the slope of the anomalous BLH can vary significantly from one period to another.

Concluding remarks
This is the first time that the diurnal/seasonal cycle of BLH is analyzed and discussed in its full detail over East Asia together with the impact of regional warming. This study aims to complement the observational studies with somewhat localized scope in space and time by providing the full spatiotemporal structure of BLH variability throughout the year. We hope that this study serves as (1) useful guideline for improving the boundary layer physics in computer models, (2) basic information on the spatiotemporal structure of boundary layer in East Asia, and (3) a benchmark solution aiding more localized and detailed analysis of the BLH and its physics. It should be mentioned that the use of reanalysis products poses crucial limitations because of uncertainty in estimating the BLH. On the other hand, BLH retrieved from the measurements from the CALISO differs from that obtained from the ground-based micropulse lidar (MPL) or the radiosonde (RS) data [51]. Further, the retrieved BLH varies according to the technique employed. Thus, the results in this study should be interpreted in the context of the inherent uncertainty in the BLH. In closing, it should be noted that the BLH climatology from the ERA5 reanalysis product in the present study may be different from those of the other reanalysis datasets and observations particularly over the Tibetan plateau.
Authors contribution KYK conceived the idea, conducted analysis, interpreted the results, and wrote the manuscript.
Funding This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea Government (NRF-2020R1A2B5B01002600).
Availability of data and materials Data used in this study are the ERA5 reanalysis products which are freely available from the ECMWF website (ecmwf.int). Results of analysis are available upon request.
Code availability CSEOF code (Fortran) is available upon request from the author.

Conflict of interest KYK declares no competing interests.
Ethical approval This study does not involve humans/animals (waiver from ethics approval).

Consent to publish Not applicable.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.