Estimating the urban atmospheric boundary layer height from remote sensing applying machine learning techniques

Estimating the urban atmospheric boundary layer height from 1 remote sensing applying machine learning techniques 2 Gregori de Arruda Moreira 1,2,3, Guadalupe Sánchez-Hernández 1,4, Juan Luis Guerrero3 Rascado 1,2 , Alberto Cazorla 1,2, Lucas Alados-Arboledas 1,2 4 1 Andalusian Institute for Earth System Research (IISTA-CEAMA), Granada, Spain 5 2 Dpt. Applied Physics, University of Granada, Granada, Spain 6 3 Federal Institute of São Paulo (IFSP), São Paulo, Brazil 7 4 Department of Physic, University of Jaén, Jaén, 2371, Spain 8 Correspondence to: Guadalupe Sánchez-Hernández ([email protected]) 9 Keywords Atmospheric Boundary Layer height, ceilometer, gradient boosting regression tree 10 Abstract 11 This study proposes a new methodology to estimate the Atmospheric Boundary Layer Height (ABLH), 12 discriminating between Convective Boundary Layer and Stable Boundary Layer heights, based on the machine 13 learning algorithm known as Gradient Boosting Regression Tree. The algorithm proposed here uses a first 14 estimation of the ABLH derived applying the gradient method to a ceilometer signal and several meteorological 15 variables to obtain ABLH values comparable to those derived from a microwave radiometer. A deep analysis of 16 the model configuration and its inputs has been performed in order to avoid the model overfitting and ensure its 17 applicability. The hourly and seasonal values and variability of the ABLH values obtained with the new algorithm 18 have been analyzed and compared with the initial estimations obtained using only the ceilometer signal. Mean 19 Relative Errors (MRE) between the ABLH estimated with the new algorithm and microwave radiometer show a 20 daily pattern with their highest values during the night-time (stable situations) and their lowest values along the 21 day-time (convective situations). This pattern has been observed for all the seasons with MRE ranging between - 22 5% and 35%. This result notably improves those ABLH values derived by applying the gradient method to 23 ceilometer data during convective situations and enables the Stable Boundary Layer height detection at night and 24 early morning, instead of only Residual Layer top height. Finally, the model performance has been directly 25 validated in three particular cases: clear-sky day, presence of low-clouds and dust outbreak event. In these three 26 particular situations, ABLH values obtained with the new algorithm follow the pattern obtained with the 27 microwave radiometer presenting very similar values, thus confirming the good model performance. In this way 28 it is feasible by the combination of the proposed method with gradient method, to estimate Convective, Stable and 29 Residual Boundary Layer height from ceilometer data and surface meteorological data in extended network that 30 include ceilometer profiling. 31 1 Introduction 32 The Atmospheric Boundary Layer (ABL) is defined as the part of the troposphere that is directly influenced by 33 the presence of the Earth’s surface, and responds to surface forcings with a timescale of about an hour or less 34 (Stull, 1988). The ABL is the atmospheric region directly affected by turbulent and evapotranspiration processes 35 and where the air pollutants are dispersed (Stull, 1988). The characteristics of the ABL, and particularly, the ABL 36 height (ABLH), play a fundamental role in numerous atmospheric areas such as weather forecasting, air quality 37 and/or numerical modeling (e.g. Cheng et al., 2011). 38 The estimation of the ABLH with high temporal resolution is not an easy task, due mainly to its high variability 39 throughout its daily cycle. Thus, based on an ideal scenario, some instants after the sunrise, the ground surface 40 temperature begins to increase, due to the positive net radiative fluxes. Such a phenomenon causes the warming 41 of air masses located at low heights favoring the convective process, and the heat transfer from the surface to 42 upper atmospheric layers in the troposphere. This process generates a layer known as Convective Boundary Layer 43 (CBL). Just before the sunset, the CBL becomes a layer called Residual Layer (RL), which is stably stratified and 44 contains the characteristics from the previous CBL. In conjunction with this process arises from the ground a 45 thermally stratified layer and endowed of lower heights (in comparison with RL and CBL), denominated Stable 46 Boundary Layer (SBL). 47 In the last years remote sensing systems, such as elastic lidars (e.g. Toledo et al., 2017; Bravo-Aranda et al., 2017; 48 Moreira et al., 2019; Vivone et al., 2021), ceilometers (e.g. Haeffelin et al., 2012; Caicedo et al., 2017; Lee et al., 49 2019; Uzan et al., 2020; Moreira et al., 2020a; Jiang et al., 2021), Doppler lidars (e.g. Manninen et al., 2018; 50 Marques et al., 2018; Moreira et al., 2019) and microwave radiometers (e.g. Cimini et al., 2013; Bravo-Aranda et 51 al., 2017; Moreira et al., 2020a; Jiang et al., 2021) have been widely used to characterize the ABLH. Among these 52 remote sensing systems, ceilometers have the advantage to be a low-cost and low-maintenance system that 53 monitors aerosol and clouds layers (Lee et al., 2019). Such characteristics have favored the creation of national 54 (e.g., Automated LIdar-CEilometer network – ALICEnet (Haefele et al., 2016); Unified Ceilometer Network - 55 UCN (National Research Council, 2009)) and international networks (e.g., EUMENET-Profiling Program 56 [https://www.eumetnet.eu/]; E-PROFILE [https://e-profile.eu]; Iberian Ceilometer Network - ICENET) (Cazorla 57 et al., 2017), which have been dedicated to standardize and expand the activities of ABL monitoring by ceilometer 58 data. 59 Ceilometers have been applied in many previous works related to ABLH detection, which vary from short-term 60 (e.g. Helmis et al., 2012; Bruine et al., 2017; Caicedo et al., 2017) to long-term studies (e.g. Stachlewska et al., 61 2012; Schween et al., 2014; Moreira et al., 2020), and various mathematical algorithms such as, vertical gradients 62 (Emeis et al., 2008), wavelet covariance transform (Baars et al., 2008; Granados-Muñoz et al., 2012), STRAT 63 (STRucture of the ATmosphere) [application of first derivative of the Gaussian filter on Range Corrected Signal 64 (RCS) profile] (Morille et al., 2007), STRAT-2D [it has same structure of STRAT and includes an edge detection 65 method based on both vertical and temporal gradients of RCS] (Haeffelin et al., 2012), STRAT+ [combination 66 of radiosoundings information and Canny edge detection applied to gradient and variance profiles of RCS] (Pal 67 et al., 2013), PathfinderTURB [it combines the strength points of gradient and variance of RCS methods and 68 addresses the layer attribution problem by adopting a geodesic approach] (Poltera et al., 2017), or COBOLT 69 (COntinuous BOundary Layer Tracing) [a time-height tracking procedure] (Geiß et al., 2017) have been 70 developed in order to improve the ABLH values derived from them. In spite of this, there are still limitations in 71 the application of ceilometers for ABLH monitoring. Special difficulties occur in cases considered as complex 72 such as rainy situations, presence of low clouds, and dust outbreaks. In these situations, the abrupt changes of the 73 aerosol vertical profile notably differ from the idealized profile on which most of the ABLH detection methods 74 are based. Although some methods have been proposed to improve the ABLH detection from lidar data in the 75 situations mentioned above (Bravo-Aranda et al., 2017; Liu et al., 2018), the use of only one wavelength and/or 76 low signal-to-noise ratio make it difficult to apply such techniques in ceilometers. Another weakness in the 77 application of ceilometers to obtain ABLH is the difficulty for discriminating between the RL top height (RLH) 78 and SBL height (SBLH) during stable periods (Moreira et al., 2020a). This limitation comes from the basis of the 79 detection procedure applied to ceilometers, based on the vertical profile of the atmospheric aerosol that prevents 80 the detection of the top of the thermal inversion (SBLH) in this situation. Having in mind these facts, there is still 81 room for some improvements in the ABLH retrieval with ceilometers on the basis of alternative data processing 82 of the ceilometer’s output. 83 Machine learning techniques have been widely applied in the environmental sciences during the last years (e.g. 84 Cadeddu et al., 2009; McGovern et al., 2017; Bonnin et al., 2018; Vassalo et al., 2020; Moreira et al., 2021). These 85 ML techniques can account for complex relationships on atmospheric processes and have been successfully 86 applied in several atmospheric areas, ranging from the estimation of atmospheric parameters such as the mixing 87 layer height (Bonin et al., 2018) or the analysis of sky-camera images to characterize the aerosol layer (Cazorla 88 et al., 2008; Cazorla et al., 2009) to predict pollutants concentration (Moreira et al., 2021). Particularly for the 89 estimation of the ABLH, Jiang et al. (2021) applied machine learning combined GPS radio occultation technology 90 to build a simulation model to estimate the ABLH, providing reliable results for several months. Krishnamurthy 91 et al. (2021) proposes a method based on Random Forest algorithm to estimate the ABLH from meteorological 92 and Doppler lidar data. Such an algorithm provides an improvement of 50%, in the CBLH detection during clear 93 sky or cloudy conditions, in comparison with a method based on vertical wind speed profiles. 94 Thus, the main objective of this study is to propose a machine learning algorithm to improve the ABLH 95 estimations obtained from a ceilometer. To this aim, a Gradient Boosting Regression Trees (GBRT) algorithm 96 has been trained using as input the ABLH values derived from a ceilometer, surface meteorological data and 97 ABLH values derived from a co-located microwave radiometer (Moreira et al., 2018; Moreira et al., 2020a). From 98 this algorithm it is possible to estimate the height of CBL and SBL combining ceilometer and meteorological 99 surface data. Therefore, application of such methodology can expand ceilometer data applicability, so that it is 100 possible to discriminate the three main ABL sublayers (CBL, SBL and RL) without acquiring expensive 101 instruments. These ABLH values derived from the microwave radiometer have been the reference dataset for both 102 the fitting and validation analysis. The model performance has been assessed analyzing its temporary and seasonal 103 variability as well as the dependence of its residuals against several meteorological variables. Finally, the model 104 has been directly validated in three particular cases (clear-sky day, presence of low-clouds and dust outbreak 105 event). 106 The paper is organized as follows. First, the experimental site and the instrumentation used in this study have been 107 described in Section 2. Then, in Section 3 the development and set-up of the machine learning algorithm here 108 proposed are presented. Next the results of a seasonal analysis as well as the analysis of particular cases are shown 109 in Section 4, while the main conclusions of this study are discussed in Section 5. 110 2 Experimental site and instrumentation 111 The measurements analyzed in this study were recorded at the University of Granada (UGR) station located on 112 the roof of the Andalusian Institute of Earth System Research (IISTA-CEAMA) at Granada (37.164º N, 3.605º 113 W, 680 m a.s.l.). These facilities are managed by the Atmospheric Physic Research Group (GFAT) and they are 114 part of the observatory AGORA (Andalusian Global ObservatoRy of the Atmosphere) in the framework of 115 ACTRIS (Aerosol, Clouds and Trace Gases Research Infrastructure) and of the Iberian Ceilometer Network 116 (ICENET) (Cazorla et al., 2017). 117 Granada is a medium sized non-industrialized city in the West Mediterranean region, at the Southeast of Spain, 118 and presents a large seasonal temperature range associated with its Mediterranean-continental conditions. The city 119 is characterized by cool winters and hot summers with the most humid period from late autumn to early spring 120 (AEMET, 2015). This region is usually affected by mineral dust outbreaks from the Sahara desert in summer and 121 spring (e.g. Guerrero-Rascado et al., 2008, 2009; Bravo-Aranda et al., 2015), and some extreme events have also 122 occurred in winter (Cazorla et al., 2017; Fernández et al., 2019), but also for more local sources of aerosol particles 123 such as traffic, domestic-heating or biomass burning in winter time (Titos et al., 2012, 2017). From roughly 124 February to July, primary biological aerosol particles (pollen-type) are present in the region (e.g., Cariñanos et 125 al., 2021). To a lesser extent, the area is also affected by advected fresh and aged smoke mainly from the Iberian 126 Peninsula (Alados-Arboledas et al., 2011) and from North America (Ortiz-Amezcua et al., 2017), respectively. 127 The combination of all these factors highly affect the local meteorology, and, therefore the ABL detection (Stull, 128 1988). 129 A ceilometer Jenoptik model CHM15k was operated at the UGR station. This instrument measures the 130 backscattered signal of a pulsed Nd:YAG laser emitting at 1064 nm, with an energy per pulse of 8.4 μJ, a repetition 131 frequency in the range of 5–7 kHz and a laser beam divergence less than 0.3 mrad. The backscattered signal is 132 received by a telescope with a field of view of 0.45 mrad. The spatial and temporal resolution used were 15 m and 133 15 s, respectively. The complete overlap of the instrument is found around 1500 m above a.g.l. and its overlap is 134 90% at 555 m a.g.l., in accordance with the overlap function provided by the manufacturer (Cazorla et al., 2017). 135 This equipment has been operating continuously since December 2012 and it is part of the Iberian Ceilometer 136 Network (ICENET), an initiative of the Atmospheric Physics Group of the University of Granada (Cazorla et al., 137 2017). Measurements recorded with this instrument were employed to derive initial estimations of the ABLH, 138 ABLHCEIL , as input (feature) of the machine learning algorithm. The range of values to ABLHCEIL [200 – 4500 139 m] (Table 1) is based on a previous long-term study (2012-2016) presented in Moreira et al. (2020a). 140 The surface meteorological dataset consists of 1-min data of air temperature (T), relative humidity (RH), 141 atmospheric pressure (P), and wind speed (WS) measured at the roof of the UGR station facilities and covering 142 the whole analyzed period 2015-2017. T and RH at this station were monitored by a HMP60 probe manufactured 143 by Vaisala. This probe has an accuracy of ±0.6 ºC and 2% for T and RH measurements, respectively. In this station 144 WS was measured by an anemometer model 05103, manufactured by Campbell Scientific, with an accuracy of 145 ±0.3 m/s. Simultaneously, P was monitored by a Vaisala PTB110 barometer with a silicon capacitive sensor 146 specially designed to guarantee accurate (±0.3 hPa at +20 °C) and stable (±0.1 hPa/year) measurements. Global 147 horizontal irradiance (G) in the range 280-2800 nm was measured by a CM-11 pyranometer manufactured by 148 Kipp & Zonen while downward infrared irradiance (IR) in the range of 4000-50000 nm is measured by a Precision 149 Infrared Radiometer (PIR) manufactured by EPPLEY. Both instruments comply with the specifications for the 150 first-class WMO pyranometer classification with an accuracy below ± 5 W/m2 for daily measurements. All these 151 sensors operated following the WMO standard protocols and procedures (WMO, 2013). These measurements and 152 magnitudes derived from them were employed as input for the machine learning algorithm developed in this study 153 as both input (features). 154 Co-located to the ceilometer, a ground-based passive microwave radiometer (MWR), model RPG-HATPRO G2 155 (Radiometer Physics GmbH) was operating in the scanning mode in automatic and continuous mode since 156 November 2011 as part of MWRnet [http://cetemps.aquila.infn.it/mwrnet/] (Rose et al., 2005; Caumont et al., 157 2016). This instrument measures the sky brightness temperature with a radiometric resolution between 0.3 and 158 0.4 K root mean square error at 1 s integration time. The MWR uses direct detection receivers within two bands, 159 22-31 GHz (water vapor - K band) and 51-58 GHz (oxygen - V band), for deriving RH and T profiles, respectively, 160 by inversion algorithms described in Rose et al. (2005). Both profiles have a range resolution varying between 10 161 and 200 m in the first 2 km and varying between 200 and 1000 m up to 10 km (Navas-Guzmán et al., 2014). This 162 change in the profile resolution is associated with an exponential decrease with height of the MWR weighting 163 functions (Spänkuch et al., 1996). The measurements recorded with this instrument were employed to derive the 164 reference values (target) of the ABLH, ABLHMRW. In the same way of the ABLHCEIL, the range of values to 165 ABLHMWR [200 – 4500 m] (Table 1) is based on Moreira et al. (2020a). The application of MWR data to ABLH 166 detection have been extensively validated with other instruments such as: Doppler lidar (Moreira et al., 2018; 167 Moreira et al., 2020b), elastic lidar (Granados-Muñoz et al., 2012; Bravo-Aranda et al., 2017; Moreira et al., 2018; 168 Moreira et al. 2020a) and radiosoundings (Bedoya-Velásquez et al., 2019). Particularly, ABLHMRW has shown to 169 be less influenced by presence of clouds (Moreira et al., 2020a) and decoupled aerosol layers (Bravo-Aranda et 170 al., 2017) compared with other devices. Moreover, the MWR temporal resolution of 2 min guarantees the volume 171 of data required for the development of the machine learning algorithm. 172 Finally, a database of hourly values of all these variables, listed in Table 1, were built for the period analyzed in 173 this study, which encompasses three entire years from 2015 to 2017. This final dataset was split in two subsets: 174 (1) a training subset, formed by measurements recorded along 2015 and 2016, and (2) a validation subset, 175 composed by measurements taken along 2017. 176 3 Methodology 177 3.1 Gradient Boosted Regression Trees 178 The Gradient Boosting Regression Trees (GBRT) is a supervised non-parametric machine learning technique 179 widely applied in classification and regression problems (e.g., Friedman 2001; Li et al. 2008; Ye et al. 2009; Chen 180 et al. 2015; Baturynska et al. 2020). The idea behind boosting is to sequentially fit multiple ‘weak learners’, that 181 is, simple models that perform relatively poorly with low accuracy (Friedman, 2001). In each iteration, a new 182 model is proposed using information from the previous model trying to learn from its mistakes and improving 183 iteration by iteration. In the case of GBRT, the ‘weak learners’ are decision tree models with very few branches. 184 Because GBRT operates with small models training sequentially, it is a faster process and requires lower memory 185 consumption than other machine learning techniques such as Random Forests. Additionally, GBRT does not 186 require the application of advanced normalization techniques on its inputs and enables a combination of different 187 numerical and categorical data as input. Similarly to the development of other types of atmospheric models, this 188 machine learning technique requires both independent variables and reference data to build the model. Particularly 189 in the machine learning vocabulary independent variables are named as features while the reference data is 190 denoted as target. The variables used in this study will be detailed described in the next sections. 191 Figure 1 shows a flowchart that briefly describes the main steps of how an ensemble of trees is created by the 192 GBRT algorithm. The ensemble consists of M trees built one-by-one. Thus, a first decision tree (Tree1) is trained 193 using the feature matrix, X (ABLHCEIL and meteorological variables), and the target variable y (ABLHMRW). The 194 predictions of Tree1 (F1(X)) are used to determine the pseudo-residual errors (r1) of the training set applying the 195 loss function L (Root Mean Square Error in our case). After that, a second decision tree (Tree2) is trained using X 196 and r1, which is the new target variable, as inputs. From the Tree2 predictions (F2(X)) the new pseudo-residual r2 197 are computed and used as input to build an improved third tree. Such a process is repeated M times until the 198 residuals are minimized and the improvement between consecutive trees is negligible. Finally, FM (X) is obtained 199 as the combination of the predicted values provided by each m-tree: 200 𝐹𝑀(𝑋)=𝐹0(𝑋)+𝐹1(𝑋)+...+𝐹𝑚−1(𝑋)+𝐹𝑚(𝑋) (1) 201 From Fig. 1 it is possible to observe that as more trees are added to the model, there is a progressive tendency to 202 reduce errors in predictions. A more detailed description of this process, including the most relevant mathematical 203 aspects, is given in Appendix 1. 204 3.2.- GBRT set up 205 3.2.1.- Inputs for the GBRT: features and target variables 206 The features initially selected to build the GBRT algorithm have been the ABLH obtained from the aerosol vertical 207 profiles measured with the ceilometer, ABLHCEIL, and set of near-surface meteorological variables which 208 influence on the ABLH as reported in previous works (e.g. Stull, 1988; Georgoulias et al., 2009; Granados-Muñoz 209 et al., 2012; Haeffelin et al., 2012 Allabakash and Lim, 2020; Rey-Sanchez et al., 2021)). ABLHCEIL have been 210 obtained applying the gradient method (Flamant et al., 1997) on the 1-hour averaged range corrected signal (𝑅𝐶𝑆 ´). 211 Considering the existence of an intense reduction in the aerosol load in the transition region between the ABL and 212 the Free Troposphere (FT), this methodology estimates the ABLH as the height (z) where the minimum in the 213 gradient of the 𝑅𝐶𝑆 ´ profile is detected (Moreira et al., 2020a). Rain cases were flagged by an empirical threshold 214 and removed (Moreira et al., 2020a). This methodology has been widely applied to numerous ceilometers 215 belonging to national or international networks such as E-PROFILE (Haefele et al., 2016) or ICENET (Cazorla 216 et al., 2017). 217 The initial near-surface meteorological dataset is composed of WS, T, P, RH, G and NR. On one hand, several 218 authors have reported significant correlations between ABLH and near-surface values of WS, T, P and RH 219 (Georgoulias et al., 2009; Wang et al., 2009; Allabakash and Lim, 2020; Krishnamurthy et al., 2021). On the other 220 hand, G accounts for the total energy reaching the surface while NR is a proxy of the brightness temperature of 221 the atmosphere, highly correlated with its composition. The solar zenith angle (SZA), the hour of the day (H) and 222 the season (S), have been also included as inputs in order to account for the Sun position and possible daily and 223 seasonal dependencies. Additionally, the clearness index (kt), estimated as the ratio between the solar radiation at 224 the top of the atmosphere and the global solar irradiance on the Earth’s surface, has been also initially considered 225 as a proxy of atmospheric transmissivity and cloudiness, respectively. In addition to their influence of these 226 variables on the ABLH they have been chosen because of their wide availability through national and international 227 meteorological and radiation networks as well as from reanalysis and satellite databases. 228 In this study, the reference values or target, also included as input in the GBRT algorithm, are the ABLH values 229 obtained from the MWR (ABLHMWR). Such an ABLH is calculated from the potential temperature profile in an 230 algorithm that combines gradient and parcel methods, for stable and convective situations, respectively. This 231 technique has been previously validated with respect to co-located elastic lidar (Granados-Muñoz et al., 2012; 232 Bravo-Aranda et al., 2017; Moreira et al. 2018) and Doppler lidar (Moreira et al. 2018; Moreira et al., 2020a) 233 presenting in both comparisons reasonable correlations with a coefficient of determination, R², above 0.7. In a 234 recent study, Bedoya-Velásquez et al. (2019) performed a validation of MWR data comparing them with 5 years 235 of radiosonde data at Granada-Spain. Such analysis demonstrated a very low bias in MWR profiles respects 236 radiosoundings, being this bias from 1.8 to −0.4 K with and standard deviation of 1.1 K for the temperature profiles 237 and from 3.0 to −4.0% with and standard deviation around 135 for the humidity profiles, under all-weather 238 conditions and below 2 km a.g.l.. 239 Additionally, from the MWR potential temperature (θ) profiles, a feature to describe the atmospheric stability 240 (AtSt) has been defined. Using the comparison criterion presented in Moreira et al. (2020), where each θ profile 241 is classified as convective, the AtSt categorical feature has been obtained being AtSt = 0 for convective situations 242 and AtSt = 1 for stable cases. 243 The initial Dataset is presented in Table 1. Hourly averages of all the relevant variables for the period 2015-2017 244 have been obtained from their original database, except for the values of H and S which were included as 245 categorical variables. Additionally, continuous variables have been normalized with respect to their mean values, 246 in order to homogenize their ranges of variability. Although this is not a required process in GBRT, Krishnamurthy 247 et al. (2021) have pointed out slight improvements in ABLH detection, mainly at nighttime, when this 248 normalization is applied. This final dataset was splitted in two subsets: (1) a subset with measurements recorded 249 in the period 2015-2016 that will be used for the model set-up and training, and (2) a validation subset, composed 250 by measurements taken along 2017. 251 3.2.2 Feature selection 252 In order to verify the relevance of each feature and to avoid data redundancy, as well as excessive complexity in 253 the model, a selection of the most relevant features from the initial dataset has been performed (Guyon and 254 Elisseeff, 2003). To this aim, the importance of each feature has been analyzed from two criteria. A first criterion, 255 namely the Boruta algorithm, estimates the importance of each feature by comparing its influence on the predicted 256 value with that of its randomly shuffled copies (Kursa et al, 2010). The second criterion, known as Recursive 257 Feature Elimination (RFE), trains a predetermined model starting with all features in the training dataset, and after 258 each iteration discards the least important features and refits the model (Yu and Liu, 2003). In this study, the 259 variables that after being discarted did not cause a 2% reduction in coefficient of determination (R²) were removed. 260 Both criteria have been applied on the entire database but also a specific feature importance analysis has been 261 performed in order to account for possible differences in the feature relevance between dayand night-time 262 situations. Figure 2 shows the relative importance of each feature, during day (a) and night (b), so that as higher 263 the value obtained, greater is the influence of this variable on the results provided by the ML model. For daytime 264 ABLHCEIL and G appear as the most relevant features while T, RH, NR, WS, P and WS show a lower relevance 265 and are sorted differently by each criteria. In the case of nighttime data, the most relevant feature is the hour (H), 266 which explains how the model can identify nighttime situations, while the ABLHCEIL takes the second position 267 and remains as one of the most important features for the model. On the opposite extreme AtSt, S, SZA, kt have 268 been classified as irrelevant features. This result can be explained by the correlation of these variables with some 269 of the features classified as relevant. Thus, i.e., all the near-surface meteorological features selected as relevant 270 present some seasonal dependence making the use of the parameter S redundant. Similarly, AtSt, appears in both 271 cases, nighttime and daytime, as one of the less relevant features. In the case of our location, this is explained 272 because nighttime/daytime classification is mostly equivalent to a stable/convective classification, making the 273 variable AtSt a redundant input. Thus, in a deep analysis of the entire database no stable cases during the daytime 274 while the 95.5% of nighttime cases are convective. These irrelevant (AtSt, S, SZA, kt) features have not been 275 included as input in the final GBRT model in order to avoid redundancy in the dataset. 276 3.2.3 Hyperparameters 277 GBRT algorithm requires a thorough setup of the so-called hyperparameters (parameters that cannot be updated 278 during the training process) in order to avoid overfitting in the training dataset. The most relevant hyperparameters 279 involved in the GBRT proposed in this study are: (1) the maximum depth of each tree, which represents the 280 maximum number of leaves in each tree, (2) the maximum number of features, which indicate the maximum 281 number of features inputted in each tree, (3) the learning rate, which indicates the influence of the previous 282 decision-trees on its successors, and (4) the minimum sample leaf, which represents the minimum number of 283 samples required to be at a leaf node in the tree. 284 In this study, the hyperparameters of the baseline model have been obtained from a large group of values randomly 285 selected over our setup-training subset, over which a cross validation and Bayesian optimization processes 286 (Frazier, 2018) have been applied using the Python library Scikit-learn (Pedregosa et al., 2011). Then, an empirical 287 fine-tuning was performed in order to detect the values that provide the best results. From this analysis, the most 288 suitable value for the maximum depth of each tree has been estimated as 5 while for the maximum number of 289 features a value of 4 has been selected. These low values of the hyperparameters contribute to reducing the 290 potential overfitting. A low value has been also obtained for the learning rate (0.0573), which ensures the 291 improvement of the correction under ceilometer data during stable periods. The optimal minimum sample leaf 292 value was indicated as 3, avoiding higher values of this parameter that can generate greater smoothing in the 293 predicted values. 294 3.3 Model training 295 Once the inputs and hyperparameters have been determined, the GBRT algorithm has been trained (stage where 296 the model is fitted) and tested (stage where the model performance is analyzed in terms of accuracy/precision). 297 As indicated in Section 2, this training has been performed using a two-year dataset (2015-2016) with 5.153 cases. 298 In order to reduce possible bias, the k-fold cross-validation methodology (James et al. 2013) has been applied. In 299 this methodology, the dataset is randomly shuffled and divided into k parts, approximately equal. Then, k 300 iterations are performed and, in each one of them, one group is selected as a test while the others k-1 are used for 301 training. After k iterations, the chosen performance parameters obtained from each iteration and mean absolute 302 error are averaged, and such values are considered as the performance parameters of the model. In this work k is 303 5 and, consequently, in each iteration an 80% and 20% of the data subset were employed for training and testing, 304 respectively. Figure 3 illustrates this process. 305 In the training stage, the model reached a R² of 0.97, which indicates a satisfactory performance and that the 306 overfitting was avoided. The Mean Absolute Error (MAE) obtained was 127 m. During the test stage, although a 307 reduction of around 20% in R² (0.76) was observed, the variation of MAE was lower than -2%, resulting in 129 308 m. 309 3.4 Analysis 310 The GBRT algorithm proposed in this study has been validated using data recorded in our station along the entire 311 2017. Thanks to that, different aspects have been analyzed. On one hand, the general performance of the algorithm 312 has been assessed analyzing the temporal and seasonal variability of the Mean Relative Error (MRE) among the 313 ABLHGBRT and ABLHMRW values. This statistic quantifies the mean relative deviation between the target value 314 (ABLHMRW) and that one provided by the model (ABLHGBRT). The MRE has been estimated by the following 315 equation: 316 𝑀𝑅𝐸𝐺𝐵𝑅𝑇(%)=100·∑(𝐴𝐵𝐿𝐻𝐺𝐵𝑅𝑇−𝐴𝐵𝐿𝐻𝑀𝑅𝑊 𝐴𝐵𝐿𝐻𝑀𝑅𝑊 )(2) 317 This statistic has been also calculated for the ABLHCEIL values in order to assess the improvement of the algorithm 318 proposed in this study with respect to the use of the ceilometer alone. 319 The statistical analysis has been completed with the estimations of the relative Root Mean Squared Error (rRMSE) 320 defined as: 321 𝑟𝑅𝑀𝑆𝐸𝐺𝐵𝑅𝑇(%)=100·√1 𝑁∑(𝐴𝐵𝐿𝐻𝐺𝐵𝑅𝑇−𝐴𝐵𝐿𝐻𝑀𝑅𝑊 𝐴𝐵𝐿𝐻𝑀𝑅𝑊 )2(3) 322 where n is the number of samples. 323 In order to identify possible limitations of the proposed algorithm under different atmospheric conditions, 324 cloudless, stable and convective situations have been differentiated and the MRE values for these situations have 325 been analyzed. Dayand nighttime have been separated in terms of the solar zenith angle values (SZA), with SZA 326 < 80º for daytime and SZA > 100º for nighttime. As mentioned above, because of the results of the 327 convective/stable analysis performed from the AtSt feature, in this study nighttime is equivalent to stable and 328 daytime is equivalent to convective situations. Additionally, cloudy and cloudless conditions have differentiated. 329 In this study, clouds have been detected from the intensity of the RCS measured by the ceilometer, which notably 330 increases in presence of clouds over the instrument. Clouds are detected when the RCS reaches values above 107, 331 which is the empirical threshold estimated for our station as representative of cloud presence (Moreira et al., 332 2020a). Dayand nighttime have been separated in terms of the solar zenith angle values (SZA), with SZA < 80º 333 for daytime and SZA > 100º for nighttime. As mentioned above, because of the results of the convective/stable 334 analysis performed from the AtSt feature, in this study nighttime is equivalent to stable and daytime is equivalent 335 to convective situations. 336 well as by the COSMO-2model. Atmos. Chem. Phys. 14, 13205–13221, https://doi.org/10.5194/acp-14-13205557 2014. 558 Emeis, S., Schäfer, K., Münkel, C., 2008. Surface-based remote sensing of the mixing-layer height – a review. 559 Meteorologische Zeitschrift, 17, 621–630.https://doi.org/10.1127/0941-2948/2008/031. 560 Eresmaa, N., Karppinen, A., Joffre, S. M., Räsänen, J., Talvitie, H., 2006. Mixing height determination by 561 ceilometer, Atmos. Chem. Phys., 6, 1485–1493, https://doi.org/10.5194/acp-6-1485-2006. 562 Fernández, A.J., Sicard, M., Costa, M.J., Guerrero-Rascado, J.L., Gómez-Amo, J.L., Molero, F., Barragán, R., 563 Basart, S., Bortoli, D., Bedoya-Velásquez, A.E., Utrillas, M.P., Salvador, P., Granados-Muñoz, M.J., Potes, M., 564 Ortiz-Amezcua, P., Martínez-Lozano, J.A., Artíñano, B., Muñoz-Porcar, C., Salgado, R., Román, R., 565 Rocadenbosch, F., Salgueiro, V., Benavent-Oltra, J.A., Rodríguez-Gómez, A., Alados-Arboledas, L., Comerón, 566 A., Pujadas, M., 2019. Extreme, wintertime Saharan dust intrusion in the Iberian Peninsula: Lidar monitoring and 567 evaluation of dust forecast models during the February 2017 event. Atmospheric Research, 228, pp. 223-241. 568 DOI: 10.1016/j.atmosres.2019.06.007. 569 Flamant, C., Pelon, J., Flamant, P.H., Durand, P., 1997. Lidar determination of the entrainment zone thickness at 570 the top of the unstable marine atmospheric boundary layer. Boundary-Layer Meteorol, 83, 247–284. 571 Frazier, P.I., 2018. A Tutorial on Bayesian Optimization. arXiv:1807.02811v1. 572 Friedman, J., 2001. Greedy boosting approximation: a gradient boosting machine. Ann. Stat. 29, 1189–1232. 573 https://doi.org/10.1214/aos/1013203451. 574 Geiß, A., Wiegner, M., Bonn, B., Schäfer, K., Forkel, R., von Schneidemesser,E., Münkel, C., Chan, K.L. and 575 Nothard, R., 2017. Mixing layer height as an indicator for urban air quality? Atmospheric Measurement 576 Techniques, 10,2969–2988. https://doi.org/10.5194/amt-10-2969-2017. 577 Georgoulias, A.K., Papanastasiou, D.K., Melas, D., Amiridis, V., Alexandri, G., 2009. Statistical analysis of 578 boundary layer heights in a suburban environment. Meteorol Atmos Phys 104, 103–111. 579 https://doi.org/10.1007/s00703-009-0021-z. 580 Granados-Muñoz, M. J., Navas-Guzmán, F., Bravo-Aranda, J. A., Guerrero-Rascado, J. L., Lyamani, H., 581 Fernández-Gálvez, J., and Alados-Arboledas, L., 2012. Automatic determination of the planetary boundary layer 582 height using lidar: One-year analysis over southeastern Spain, J. Geophys. Res.-Atmos., 117, D18208, 583 https://doi.org/10.1029/2012JD017524. 584 Guerrero-Rascado, J.L., Ruiz, B., Alados-Arboledas, L., 2008. Multi-spectral Lidar characterization of the vertical 585 structure of Saharan dust aerosol over southern Spain. Atmospheric Environment, 42 (11), pp. 2668-2681. DOI: 586 10.1016/j.atmosenv.2007.12.062. 587 Guerrero-Rascado, J. L., Olmo, F. J., Avilés-Rodríguez, I., Navas-Guzmán, F., Pérez-Ramírez, D., Lyamani, H., 588 and Alados Arboledas, L., 2009. Extreme Saharan dust event over the southern Iberian Peninsula in september 589 2007: active and passive remote sensing from surface and satellite, Atmos. Chem. Phys., 9, 8453–8469, 590 https://doi.org/10.5194/acp-9-8453-2009. 591 Guyon, I., Elisseeff, A., 2003. An introduction to variable and feature selection. Journal of Machine Learning 592 Research, 3, 1157-1182. 593 Haefele, A., Hervo M., Turp, M., Lampin, J.-L., Haeffelin, M., Lehmann, V., 2016. E-PROFILE team, TOPROF 594 team. The E-PROFILE network for the operational measurement of wind and aerosol profiles over Europe. Teco 595 2016 - Madri (Spain). Available at: https://www.eumetnet.eu/wp-content/uploads/2016/10/E596 PROFILE_TECO_Madrid_2016.pdf. 597 Haeffelin, M., Angelini, F., Morille, Y., Martucci, G., Frey, S., Gobbi, G.P.,Lolli, S., O’Dowd, C.D., Sauvage, L., 598 Xueref-Rémy, I., Wastine, B., Feist, D.G., 2012. Evaluation of mixing-height retrievals from automatic profil-ing 599 lidars and ceilometers in view of future integrated networks in Europe.Boundary-Layer Meteorology, 143, 49– 600 75. https://doi.org/10.1007/s10546-011-9643-z. 601 Helmis, C.G.; Sgouros, G.; Tombrou, M.; Schäfer, K.; Münkel, C.; Bossioli, E.; Dandou, A., 2012. A Comparative 602 Study and Evaluation of Mixing-Height Estimation Based on Sodar-RASS, Ceilometer Data and Numerical 603 Model Simulations. Bound.-Layer Meteorol., 145, 507–526. 604 Iqbal M., 1983. An Introduction to Solar Radiation, Academic Press, New York. 605 James, G., Witten, D., Hastie, T., Tibshirani, R., 2013. An Introduction to Statistical Learning with Application 606 in R (Springer Texts in Statistics) 7th ed. Springer New York Heidelberg Dordrecht London. 607 Jiang, R., Zhao, K., 2021. Using machine learning method on calculation of boundary layer height. Neural Comput 608 & Applic. https://doi.org/10.1007/s00521-021-05865-3. 609 Jiang, Y., Xin, J., Zhao, D., Jia, D., Tang, G., Quan, J., Wang, M., Dai, L., 2021. Analysis of differences between 610 thermodynamic and material boundary layer structure: Comparison of detection by ceilometer and microwave 611 radiometer. Atmospheric Research, 248, 105179, https://doi.org/10.1016/j.atmosres.2020.105179. 612 Krishnamurthy, R., Newsom, R. K., Berg, L. K., Xiao, H., Ma, P.-L., Turner, D. D., 2021. On the estimation of 613 boundary layer heights: A machine learning approach, Atmos. Meas. Tech. Discuss. [preprint], 614 https://doi.org/10.5194/amt-2020-439. 615 Kursa, M.B., Jankowski, A., Rudnicki, W.R., 2010. Boruta - A System for Feature Selection. Fundamenta 616 Informaticae, 101, 271–285. https://doi.org/10.3233/FI-2010-288. 617 Lee, J., Hong, J.W., Lee, K., Hong, J., Velasco, E., Lim, Y.J., Lee, J.B., Nam, K., Park, J., 2019. Ceilometer 618 Monitoring of Boundary-Layer Height and Its Application in Evaluating the Dilution Effect on Air Pollution. 619 Boundary-Layer Meteorol 172, 435–455. https://doi.org/10.1007/s10546-019-00452-5. 620 Li, P., Wu, Q., Burges, C. J., 2008. Mcrank: Learning to rank using multiple classification and gradient boosting. 621 In Advances in Neural Information Processing Systems 20, pages 897–904. 622 Liu, B., Ma, Y., Gong, W., Yang, J., Zhang, M., 2018. Two-wavelength Lidar inversion algorithm for determining 623 planetary boundary layer height. J. Quant. Spectrosc. Radiat. Transf., 206, 117–124. 624 Manninen, A. J., Marke, T., Tuononen, M. J., O'Connor, E. J., 2018. Atmospheric boundary layer classification 625 with Doppler lidar. Journal of Geophysical Research: Atmospheres, 123, 8172– 8189. 626 https://doi.org/10.1029/2017JD028169. 627 Marques, M. T. A., Moreira, G. de A., Pinero, M., Oliveira, A. P., Landulfo, E., 2018. Estimating the planetary 628 boundary layer height from radiosonde and doppler lidar measurements in the city of São Paulo - Brazil. EPJ 629 WEB OF CONFERENCES, v. 176, p. 06015. https://doi.org/10.1051/epjconf/201817606015. 630 McGovern, A., Elmore, K.L., Gagne, D.J., Haupt, S.E., Karstens, C.D., Lagerquist, R., Smith, T., Williams, J.K., 631 2017. Using artificial intelligence to improve real-time decision-making for high-impact weather. Bulletin of the 632 American Meteorological Society, 98(10), pp.2073-2090. 633 Moreira, G.A., Guerrero-Rascado, J.L., Bravo-Aranda, J.A., Benavent-Oltra, J.A., Ortiz-Amezcua, P., Róman, 634 R., Bedoya-Velásquez, A.E., Landulfo, E., Alados-Arboledas, L., 2018. Study of the planetary boundary layer by 635 microwave radiometer, elastic lidar and Doppler lidar estimations in Southern Iberian Peninsula. Atmospheric 636 Research, 213, 185-195. https://doi.org/10.1016/j.atmosres.2018.06.007. 637 Moreira, G.A., Guerrero-Rascado, J.L., Benavent-Oltra, J.A., Ortiz-Amezcua, P., Román, R., Bedoya-Velásquez, 638 A.E., Bravo-Aranda, J.A., Olmo-Reyes, F.J., Landulfo, E., Alados-Arboledas, L., 2019. Analyzing the turbulent 639 planetary boundary layer by remote sensing systems: the Doppler wind lidar, aerosol elastic lidar and microwave 640 radiometer. ATMOSPHERIC CHEMISTRY AND PHYSICS (ONLINE), v. 19, p. 1263-1280. 641 Moreira, G.A., Guerrero-Rascado, J.L., Bravo-Aranda, J.A., Foyo-Moreno, I., Cazorla, A., Alados, I., Lyamani, 642 H., Landulfo, E., Alados-Arboledas, L., 2020a. Study of the planetary boundary layer height in an urban 643 environment using a combination of microwave radiometer and ceilometer. Atmospheric Research, 240, 104932, 644 https://doi.org/10.1016/j.atmosres.2020.104932. 645 Moreira, G. A., Da Silva Lopes, F.J., Guerrero-Rascado, J.L., Ortiz-Amezcua, P., Cazorla, A., De Oliveira, A.P., 646 Landulfo, E., Alados-Arboledas, L., 2020b. Comparison Among the Atmospheric Boundary Layer Height 647 Estimated From Three Different Tracers. EPJ WEB OF CONFERENCES, v. 237, p. 03009, 648 https://doi.org/10.1051/epjconf/202023703009. 649 Moreira, G.A., Andrade, I.S., Cacheffo, A., Yoshida, A.C., Gomes, A.A., Silva, J.J., Lopes, F.J.S., Landulfo, E., 650 2021. COVID-19 outbreak and air quality: Analyzing the influence of physical distancing and the resumption of 651 activities in São Paulo municipality, Urban Climate, Volume 37, 100813, ISSN 2212-0955, 652 https://doi.org/10.1016/j.uclim.2021.100813. 653 Morille, Y., Haeffelin, M., Drobinski, P., Pelon, J., 2007. STRAT: An Automated Algorithm to Retrieve the 654 Vertical Structure of the Atmosphere from Single-Channel Lidar Data, Journal of Atmospheric and Oceanic 655 Technology 24, 5 , 761-775, https://doi.org/10.1175/JTECH2008.1. 656 Müller, A.C., Guido, S., 2016. Introduction to machine learning with Python. O’Reilly Media, Inc., Sebastopol. 657 National Research Council. 2009. Observing Weather and Climate from the Ground Up: A Nationwide Network 658 of Networks. Washington, DC: The National Academies Press. https://doi.org/10.17226/12540. 659 Ortiz-Amezcua, P., Guerrero-Rascado, J. L., Granados-Muñoz, M. J., Benavent-Oltra, J. A., Böckmann, C., 660 Samaras, S., Stachlewska, I. S., Janicka, Ł., Baars, H., Bohlmann, S., Alados-Arboledas, L., 2017. Microphysical 661 characterization of long-range transported biomass burning particles from North America at three EARLINET 662 stations, Atmos. Chem. Phys., 17, 5931–5946, https://doi.org/10.5194/acp-17-5931-2017. 663 Pal, S., Haeffelin, M., Batchvarova, E., 2013. Exploring a geophysical process-based attribution technique for the 664 determination of the atmospheric boundary layer depth using aerosol lidar and near-surface meteorological 665 measurements. Journal of Geophysical Research: Atmospheres, 118, 9277– 666 9295.https://doi.org/10.1002/jgrd.50710. 667 Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B. Grisel, O., Blondel, M., Prettenhofer, P., 668 Weiss, R., Dubourg, V., Vanderplas, J., Passos, A, Cournapeau, D., Brucher, M., Perrot, M., Duchesnay, E., 2011. 669 Scikit-learn: machine learning in Python. Journal of Machine Learning Research, 12, pp. 2825-2830. 670 Poltera, Y., Martucci, G., Collaud Coen, M., Hervo, M., Emmenegger, L., Henne, S., Brunner, D., Haefele, A., 671 2017. PathfinderTURB: an auto-matic boundary layer algorithm. Development, validation and application to 672 study the impact on in situ measurements at the Jungfraujoch. Atmospheric Chemistry and Physics, 17, 10051– 673 10070. https://doi.org/10.5194/acp-17-10051-2017. 674 Rey-Sanchez, C., Wharton, S., Vilà-Guerau de Arellano, J., Paw U, K. T., Hemes, K. S., Fuentes, J. D., Osuna, 675 J., Szutu, D., Ribeiro, J. V., Verfaillie, J., Baldocchi, D., 2021. Evaluation of atmospheric boundary layer height 676 from wind profiling radar and slab models and its responses to seasonality of land cover, subsidence, and 677 advection. Journal of Geophysical Research: Atmospheres, 126, e2020JD033775. 678 https://doi.org/10.1029/2020JD033775. 679 Schween, J.H., Hirsikko, A., Löhnert, U., Crewell, S., 2014. Mixing-layer height retrieval with ceilometer and 680 Doppler lidar: From case studies to long-term assessment. Atmos. Meas. Tech., 7, 4275–4319. 681 Stachlewska, I.S., Migacz, S., Szkop, A., Zielínska, A.J., Swaczyna, P.L., 2012. Ceilometer observations of the 682 boundary layer over Warsaw, Poland. Acta Geophys. 60, 1386–1412. 683 Stull, R. B., 1988. An Introduction to Boundary Layer Meteorology, 666 pp., Kluwer Acad., Dordrecht, 684 Netherlands. 685 Toledo, D., Córdoba-Jabonero, C., Adame, J.A., Benito, D.L.M., Gil-Ojeda, M., 2017. Estimation of the 686 atmospheric boundary layer height during different atmospheric conditions: A comparison on reliability of several 687 methods applied to lidar measurements. Int. J. Remote Sens., 38, 3203–3218. 688 Uzan, L., Egert, S., Khain, P., Levi, Y., Vadislavsky, E., Alpert, P., 2020. Ceilometers as planetary boundary layer 689 height detectors and a corrective tool for COSMO and IFS models, Atmos. Chem. Phys., 20, 12177–12192, 690 https://doi.org/10.5194/acp-20-12177-2020. 691 Vassallo, D., Krishnamurthy, R., Fernando, H. J. S., 2021. Utilizing physics-based input features within a machine 692 learning model to predict wind speed forecasting error, Wind Energ. Sci., 6, 295–309, 693 https://doi.org/10.5194/wes-6-295-2021. 694 Vivone, G., D'Amico, G., Summa, D., Lolli, S., Amodeo, A., Bortoli, D., Pappalardo, G., 2021. Atmospheric 695 boundary layer height estimation from aerosol lidar: a new approach based on morphological image processing 696 techniques, Atmos. Chem. Phys., 21, 4249–4265, https://doi.org/10.5194/acp-21-4249-2021. 697 Wang, W., Mao, F., Gong, W., Pan, Z., Du, L., 2016. Evaluating the Governing Factors of Variability in Nocturnal 698 Boundary Layer Height Based on Elastic Lidar in Wuhan. International Journal of Environmental Research and 699 Public Health. 13(11):1071. https://doi.org/10.3390/ijerph13111071. 700 Ye, J., Chow, J.H., Chen, J., Zheng, Z., 2009. Stochastic gradient boosted distributed decision trees. In 701 Proceedings of the 18th ACM conference on information and knowledge management (pp. 2061–2064), ACM. 702 https://doi.org/10.1145/1645953.1646301. 703 Yu, L., Liu, H., 2003. Feature Selection for High-Dimensional Data: A Fast Correlation-Based Filter Solution. 704 Proceedings of the Twentieth International Conference on Machine Learning (ICML-2003), Washington DC. 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 Table 1. Group of input variables initially considered along with the instrument used to measured them and the 735 range of variation of each variable during the period of study. 736 Instrument/Algorithm Variable Range Ceilometer/Gradient Method ABLHCEIL 200 m - 4500 m MWR/Gradient and Parcel Method ABLHMWR 200 m - 4500 m HMP60 Temperature (T) 0 - 42 ºC Hour (H) Categorical Variable Season (S) Categorical Variable Stability (AtSt) Categorical Variable (0-1) HMP60 Relative Humidity (RH) 4.7 - 91 % Barometer PTB110 Pressure (P) 920 - 952 hPa Anemometer 05103 Wind Speed (WS) 0 - 5 m/s Pyranometer CM-11 Global Radiation (G) 0 - 1016 W/m² Pyrgeometer PIR Net Radiation (NR) -167 - (-1) W/m² Blanco-Muriel et al. (2001) Solar Zenith Angle (SZA) 14 - 165 º Iqbal (1983) Clearness Index (kt) 0 - 1 737 738 Table2. rRMSE of the GBRT for all cases, as well as for each season, for both dayand nigh-time 739 (convective/stable) situations. 740 All cases Winter Spring Summer Autumn Day Night Day Night Day Night Day Night Day Night rRMSEGBRT (%) 15 25 18 24 14 25 10 20 15 26 741 742 Table 3. MRE and R² of the GBRT algorithm under cloudless and all-cloud-type conditions. Additionally, for 743 each category, stable, convective and all-stability conditions have been differentiated. Number of cases on each 744 category have been included in order to prove their representativeness. 745 All Cases Cloudless Cases Stable Convective All cases Stable Convective All cases Number of cases 1284 1579 2863 398 600 998 R² 0.75 0.89 0.90 0.75 0.89 0.91 MREGBRT (%) 28.0 11.1 20.9 27.5 10.2 19.1 746 747 Figure 1 - GBRT flowchart. X, y, and L represent the feature matrix, target variable, and loss-function, 748 respectively. rN and FN indicate de n-nth pseudo-residual error and prediction. 749 750 751 Figure 2 - Feature relative importance classification/ranking applying the Recursive Feature Elimination (RFE) 752 method, for day (a) and night (b) situations. 753 754 755 756 757 Figure 3 - Scheme of input dataset (left) and k-fold cross-validation methodology (right). 758 759 760 761 Figure 4 – Hourly Mean Relative Error for all the analyzed cases applied in the GBRT algorithm (black) and the 762 gradient method to the ceilometer data (pink). It should be highlighted the important difference between the scales 763 required for each methodology. 764 765 766 767 Figure 5 - Comparison between the hourly ABLH average measured (red line) and those predicted by the GBRT 768 algorithm (black line) and the ceilometer (pink) for (a) winter, (b) spring, (c) summer, and (d) autumn. The dark 769 shadow represents the GBRT model standard deviation. 770 771 772 Figure 6 - Hourly Mean Relative Error for all the analyzed cases applied in the GBRT algorithm (black) and the 773 gradient method to the ceilometer data (pink) during (a) winter, (b) spring, (c) summer, (d) autumn. 774 775 776 777 778 Figure 8 - Comparison among the hourly values of ABLHGBRT (black stars), ABLHMWR (red stars) and of 779 ABLHCEIL (pink stars) at January 24, 2017. 780 781 782 783 Figure 9 - Comparison among the hourly values of ABLHGBRT (black stars), ABLHMWR (red stars) and ABLHCEIL 784 (pink stars) at February 7, 2017. 785 786