A Satellite-Based Spatio-Temporal Machine Learning Model to Reconstruct Daily PM2.5 Concentrations across Great Britain

Epidemiological studies on the health effects of air pollution usually rely on measurements from fixed ground monitors, which provide limited spatio-temporal coverage. Data from satellites, reanalysis, and chemical transport models offer additional information used to reconstruct pollution concentrations at high spatio-temporal resolutions. This study aims to develop a multi-stage satellite-based machine learning model to estimate daily fine particulate matter (PM2.5) levels across Great Britain between 2008–2018. This high-resolution model consists of random forest (RF) algorithms applied in four stages. Stage-1 augments monitor-PM2.5 series using co-located PM10 measures. Stage-2 imputes missing satellite aerosol optical depth observations using atmospheric reanalysis models. Stage-3 integrates the output from previous stages with spatial and spatio-temporal variables to build a prediction model for PM2.5. Stage-4 applies Stage-3 models to estimate daily PM2.5 concentrations over a 1 km grid. The RF architecture performed well in all stages, with results from Stage-3 showing an average cross-validated R2 of 0.767 and minimal bias. The model performed better over the temporal scale when compared to the spatial component, but both presented good accuracy with an R2 of 0.795 and 0.658, respectively. These findings indicate that direct satellite observations must be integrated with other satellite-based products and geospatial variables to derive reliable estimates of air pollution exposure. The high spatio-temporal resolution and the relatively high precision allow these estimates (approximately 950 million points) to be used in epidemiological analyses to assess health risks associated with both short- and long-term exposure to PM2.5.


PM2.5 and PM10 observed data
Daily PM10 and PM2.5 (µg/m 3 ) measurements in the study period were obtained from five monitoring network sources through the R package openair [29]: Automatic Urban and Rural Network (AURN), Air Quality England (AQE), King College London (KCL), Scotland Air Quality Network (SAQN), and Wales Air Quality Network (WAQN). When monitors from different sources showed the same temporal distributions (i.e., a correlation equal to 1) and were located at approximately the same coordinates, only the AURN monitor was kept. Monitors with less than 18 hours of PM2.5 records per day as well as less than 30 days by year were removed. The final set includes 581 and 183 monitors measuring PM10 and PM2.5 along the study period, respectively. Figure  1 shows the locations of these monitors across Great Britain, illustrating how the network coverage is densely located in major cities, leaving several small cities and rural areas with only a few or no AQ records. Each monitor was indexed using the cell-ID of the 1 km grid cell that contained it.

Spatially-Lagged and Nearest Monitor PM 2.5 Variables
Four spatio-temporal variables were generated from the monitor series of PM 2.5 to represent spatially-lagged annual average concentrations. Monitor types were grouped into two classes (background (urban, suburban, and rural) and hotspots (traffic and industrial)) to compute the annual average values of nearby monitors by class using an inverse-distance weighted leave-one-out cross-validated (IDW-LOOCV) approach. Two different weights were applied, namely the inverse distance and the inverse squared distance (in km). The former assigns relatively more weight to distant monitors and therefore can represent a regional background, while the latter captures local differences. Therefore, these four IDW-LOOCV variables were named as follows: (i) Spatially-lagged hotspot-PM 2.5 regional, (ii) Spatially-lagged background-PM 2.5 regional, (iii) Spatially-lagged hotspot-PM 2.5 local, and (iv) Spatially-lagged background-PM 2.5 local. To improve the model performance in the spatial domain, two additional spatial variables were generated based on the closest Euclidean distance for each monitor class, named as: (i) Nearest hotspot monitor distance and (ii) Nearest background monitor distance. These six variables were used as additional predictors to capture the heterogeneity across monitors and exploit their spatial autocorrelation, and thus help the model to better categorise the differences in the spatial patterns of measured-PM 2.5 series.

AOD Data: Satellite and Atmospheric Reanalysis Models
Daily satellite-AOD was obtained from the Collection 6 Level-2 gridded product (MCD19A2). These data are generated at a 1 km grid through the Multi-angle Implementation of Atmospheric Correction (MAIAC) algorithm using data from a Moderate Resolution Imaging Spectroradiometer (MODIS) sensor on board both Terra and Aqua Earth Observation satellites [30]. Four layers from the MCD19A2 product were extracted: (i) AOD Blue band (0.47 µm), (ii) AOD Green band (0.55 µm), (iii) AOD uncertainty (i.e., the level of uncertainty based on blue-band surface brightness (reflectance)), (iv) AOD_QA (quality assurance flags to retrieve only the best quality AOD). These layers are generated for each passing time of Terra and Aqua satellites over the area of study and combined by day. The layers AOD-0.47 µm and AOD-0.55 µm were used as the outcome variables after their values were filtered using AOD uncertainty and AOD_QA layers to guarantee high product quality.
This calibration process, together with pixels covered by clouds, removed a large sample of AOD grid cells over Great Britain. To fill the satellite-AOD gaps, the modelled-AOD total column was used from Copernicus Atmosphere Monitoring Service (CAMS) reanalysis provided by the European Centre for Medium-Range Weather Forecasts (ECMWF) [31]. CAMS reanalysis provides every 3-hourly modelled-AOD at five different wavelengths (0.47 µm, 0.55 µm, 0.67 µm, 0.87 µm, and 1.24 µm) with a spatial resolution of approximately 80 km, but the data were downloaded at 10 km, based on an interpolation performed through the ECMWF's API request. The satellite-AOD and CAMS modelled-AOD were indexed to the closest 1 km grid cell from their pixel centroid.

Modelled PM 2.5 from Chemical Transport Models
Atmospheric chemistry transport models (ACTMs) incorporate anthropogenic and natural sources of emission, land use, and meteorological conditions to simulate the atmospheric compositions and deposition of various air pollutants (trace gases and particles). Based on the European Modelling and Evaluation Programme (EMEP) ACTM, EMEP4UK, has been developed to represent the UK hourly atmospheric composition at a spatial resolution of approximately 5 km [32]. The description of the EMEP4UK model framework and setup can be found elsewhere [33,34]. Daily EMEP4UK simulations of PM 2.5 (µg/m 3 ) concentrations at the surface-level were included to represent ground-level contributions, in contrast to AOD products that refer to the total column of aerosol concentration. Each EMEP4UK 5 km pixel was linked to the closest 1 km grid cell centroid.

Meteorological Variables from Climate Reanalysis Models
Meteorological variables were retrieved from the ECMWF's climate reanalysis models with the highest spatial resolution available during 2008-2018 and at two sub-day times (0:00 and 12:00). Sea-level pressure and the boundary layer height (BLH) were downloaded from the ERA 5 global reanalysis with a spatial resolution of approximately 30 km [35]. Air temperature at 2 m height and total precipitation were obtained from the ERA 5 Land global reanalysis with a spatial resolution of approximately 9 km [36]. Relative humidity, wind direction, and wind speed were downloaded from the UERRA regional reanalysis at 5.5 km for the MESCAN-SURFEX system [37]. All meteorological variables were indexed to the closest 1 km grid cell to their centroid.

Normalized Difference Vegetation Index
Monthly Normalized difference vegetation index (NDVI) is used to quantify vegetation presence and it ranges from −1 to 1. The NDVI 1 km grid was obtained from MOD13A3 Version 6 Level 3, Terra-MODIS product [38]. Each NDVI pixel was indexed to the closest 1 km grid cell to its pixel centroid and the NDVI values repeated for the days inside the corresponded month.

Land Variables and Night-Time Light Data from Earth Observation Satellites
Three land predictors were collected from the Copernicus Land Monitoring Service (CLMS) [39] database: elevation, land cover, and impervious surfaces. Elevation data were obtained from the 2011 European Digital Elevation Model (EU-DEM) version 1.1 with a spatial resolution of 25 m. The elevation values were obtained from the mean of all 25 m-pixel values located inside each 1 km grid cell. Land cover data were obtained from the 2012 CORINE Land Cover (CLC) inventory. It was derived from high-resolution ortho-rectified satellites images that mapped all land elements at a spatial resolution ranging from 5 m to 60 m and aggregated into 100 m. Nine predictors were defined by grouping the original 44 CLC classes and each predictor represents its group proportion inside each 1 km Great Britain Grid cell. The imperviousness degree is a binary raster product at a spatial resolution of 100 m, where the value 0 represents natural land cover or water surface and value 1 represents entirely artificial surfaces (i.e., built-up areas). The amount of impervious surfaces was estimated by the proportion of artificial surfaces inside each 1 km Great Britain Grid cell.
Night-time lights data were provided by the visible Infrared Imaging Radiometer Suite (VIIRS) sensor aboard the Suomi-National Polar-orbiting Partnership (Suomi-NPP) satellite. The VIIRS Day/Night band collects cloud-free average radiance values at annual and monthly composites in a spatial resolution of 750 m [40]. The 2015 night-time lights annual mean was computed based on the weighted average of VIIRS pixels inside each 1 km Great Britain Grid cell.

Population Density
The resident population counts data from 2011 that were collected from the Office for National Statistics (England and Wales) and the National Records of Scotland. The smallest geographic unit of the UK census is output areas (OA), with a total of 227.769 OAs polygons for Great Britain. The proportion of each OA's area inside each 1 km Great Britain Grid cell was extracted to estimate the weighted average population by 1 km grid cell.

Road Density and Distance
Road density and length predictors were derived from the Ordnance Survey (OS) [41] Open Roads product, which offers a geospatial representation of Great Britain's Road network. Three density predictors were defined by grouping the original eight OS road types, where each predictor was computed as the sum of all roads length inside each 1 km Great Britain Grid cell by group (highway, secondary, and local). Three distance predictors were defined by computing the inverse distance of each 1 km grid centroid from the closest road group (highway, secondary, and local).

Inverse Distance from Airports and Seashore
Information on location and size of airports was derived from the Civil Aviation Authority (CAA) [42], which collects monthly statistics about air traffic movements for more than 60 UK Airports, including aviation activities for terminal passengers, commercial flights, and cargo tonnage. A total of 19 airports across Great Britain were selected based on a minimum of 1% for the annual percentage of passengers at the airport from 2015-2018. For each cell of the 1 km Great Britain grid, the inverse distance from the closest airport was calculated.
The inverse distance from the seashore was computed for each 1 km grid cell using the geographical information about the boundaries of England, Wales, and Scotland provided by the UK Data Service [43].

Statistical Methods
A four-stage model was developed to obtain daily PM 2.5 concentrations for all 234.429 grid cells covering Great Britain. Each stage is described in detail below. Briefly, Stage-1 applies a random forest (RF) algorithm to predict PM 2.5 concentrations in monitors with only PM 10 records. Stage-2 uses RF models to impute missing satellite-AOD from Terra and Aqua satellites, using modelled-AOD from CAMS. Stage-3 combines the output from Stage-1 and Stage-2 with a list of spatial and spatio-temporal synchronised predictors to estimate PM 2.5 concentrations at the locations of the monitors. Stage-4 uses the Stage-3 model to predict daily PM 2.5 across the whole of Great Britain.

Random Forest Algorithm
RF is a supervised tree-based design ML algorithm that trains an ensemble of independent decision trees (or forests) in parallel. The final model accuracy is estimated by the performance average of all decision trees. There are two main advantages of the RF architecture (known as the bagging ensemble method): first, it controls the bias-variance trade-off by feeding each tree model with two-thirds of the training set while one-third is left out for validation (i.e., out-of-bag (OOB)) [44]. When all decision trees receive the same amount and list of predictors, they become highly correlated, not solving the variance problem. Therefore, the second positive aspect of RF is that the number of predictors on each tree is less than the full list available and they are randomly selected. This approach changes the predictor placed on the top of the tree, generating different splits and internal nodes. The algorithm is then able to estimate an importance ranking by quantifying the amount of error decreased due to a split of a specific predictor [45].
The performance of the models in each stage was assessed using statistics based on OOB samples and then from a 10-fold cross-validation (CV) procedure based on monitors. In the latter, ten random groups of monitors were defined, and the complete outcome series in each group were predicted using a model fitted in the other nine. This procedure offers a measure of the true predictive ability of the RF model in locations where no monitor is available. Measures of performance were generated by regressing the OOB or cross-validated predicted values on the observed series, and computing the R 2 , the root mean square error (RMSE), and the intercept and slope of the prediction. These statistics were computed overall using the whole series and then separated into spatial and temporal contributions. The former was computed using the averages of predicted and observed values across the series, and it offers a measure of performance in capturing long-term average PM 2.5 concentrations. The latter was computed as daily deviations from the averages, and it quantifies the temporal variability explained by the model.

Stage-1: Increasing PM 2.5 Measurements Using Co-located PM 10 Monitors
The number of PM 2.5 monitors across Great Britain at the beginning of the study period was relatively low, and even if the quantity has increased substantially after 2010, most of the monitors were mostly installed in major cities. Stage-1 aims to increase the number of spatio-temporal ground-level PM 2.5 references using observations from co-located monitors measuring PM 10 , which are available in a higher number of locations and are better distributed across Great Britain. Specifically, in this stage, we fitted an RF model in locations with both PM 2.5 and PM 10 measurements, using only data from the specific year-model (2008-2018). The RF model for each year y is defined as: where: PM y 2.5 i,t and PM 10 i,t are the target variable and main predictor, respectively, measured in year y at monitor i on day t; mon.type i is a categorical variable classifying monitor i (traffic, industrial, urban, suburban, and rural); month t and dow t are categorical variables representing the months and day of the week of day t; and lat i and lon i define the coordinates of monitor i. The RF model was defined based on the best parameter setting. The optimised parameters were 500 decision trees as the RF ensemble (Ntree = 500) and 4 variables randomly selected to be used on each tree (mtry = 4). This model was eventually used to predict PM 2.5 measurements in locations/days in which only PM 10 was measured.

Stage-2: Imputing Missing Satellite-AOD from CAMS Modelled-AOD
The percentage of missing satellite-AOD measurements in Great Britain ranged between 87% and 94% during 2008-2018 with the greatest portion during autumn and winter, near to the coast, and in the North of Scotland. Stage-2 imputes satellite-AOD missing for every day and 1 km grid based on an optimised RF model (Ntree = 50 and mtry = 20) and satellite-AOD wavelength (0.47 µm and 0.55 µm), separately in each year within the study period. These RF models were built for each year y as follows: where: Satellite-AOD z, y i,t is the target variable representing satellite-AOD (wavelength z, year y) estimates at grid cell i on day t; CAMS.AOD is the main predictor, representing CAMS modelled-AOD estimates at grid cell i, on day t, at five wavelengths (0.47 µm, 0.55 µm, 0.67 µm, 0.865 µm, and 1.24 µm), and at seven sub-day times (3 h, 6 h, 9 h, 12 h, 15 h, 18 h, and 21 h); yday t defines the sequence of days in a year from 1 to 365 (366, for leap years); lat i and lon i represent the coordinates of grid cell centroid i. This model was eventually used to predict missing satellite-AOD measurements.

Stage-3: Estimating PM 2.5 Concentrations Using Spatial and Spatio-Temporal Variables
Stage-3 aims to build a predictive model for daily PM 2.5 concentrations, using as target variable the combined set of PM 2.5 directly measured from monitors or predicted from Stage-1, AOD measured from satellite instruments or predicted from Stage-2, together with all the spatially and spatio-temporally synchronized predictors described in the previous section. The RF models were fit separately by year y. The optimization process for this stage selected Ntree = 500 and mtry = 20 as the best set of parameter values. In this stage, PM 2.5 concentrations were log-transformed to ensure the prediction of non-negative values. The Stage-3 model for each year y is defined as: where: log PM 2.5 y i,t is the target variable representing the log-PM 2.5 concentrations in year y at the monitor located in grid cell i on day t, while SPT i,t and SP i represent the spatio-temporal and spatial predictors, respectively. Specifically, SPT i,t included: surface-level log-PM 2.5 values from the EMEP4UK model; Stage-2 AOD (at 0.47 µm, 0.55 µm); meteorological variables (air temperature, sea pressure, relative humidity, total precipitation, wind speed, and direction); BLH (at 12:00 and 24:00 hours); monthly-averaged NDVI; and day of the year, month, and week. SP i included: spatially-lagged regional and local log-PM 2.5 average concentrations from background and hotspot groups; nearest background and hotspot monitor distance; land variables (elevation, land cover in 9 groups, and % of impervious surface); night-time light; road density and inverse distance (each for highway, secondary, and local); population density; and inverse distance from closest airport and seashore.   Table 1 shows the Stage-1 model performance (i.e., a separate RF model for each year) reported by two CV methods: (i) OOB, the original RF CV algorithm, which presented better accuracy (R 2 average of 0.932), and (ii) 10-Fold CV, a targeted sampling, which leaves out monitors with the full set of observations (R 2 average of 0.855). As expected, the 10-Fold CV approach resulted in a slightly lower predictive accuracy although it still showed a good performance. There is some indication of bias in the 10-Fold CV, with a slope lower than the expected value of 1 and a slightly negative intercept. The results for 10-Fold CV spatial and temporal domains are in Table S1 in Supplementary Materials. Table 1. Predicted-PM2.5 concentrations obtained from Stage-1 RF models were regressed against measured-PM2.5 concentrations in a linear regression model. The performance was evaluated using two CV methods (OOB and 10-Fold) together with RMSE (a measure of the model error, µg/m3), intercept (µg/m3), and slope (µg/m3).

Stage-2 Results
The Stage-2 procedure is illustrated in Figure 2, showing missing satellite-AOD, imputed through a combination of multiple modelled-AOD wavelengths and sub-day times. Table 2 shows the performance of Stage-2 models validated using the OOB CV method. The results presented consistently high R 2 (ranging from 0.963 to 0.988), low RMSE (varying between 0.007 to 0.010), and almost no bias (intercept zero and slope close to 1 in almost all years and wavelengths).

Stage-3 Results
Stage-3, the main step of the satellite-based machine learning framework, combines the output

Stage-3 Results
Stage-3, the main step of the satellite-based machine learning framework, combines the output of Stage-1 and Stage-2 with a list of spatial and spatio-temporal predictors to estimate PM 2.5 at the locations of the monitors. The relative importance of the predictors for the Stage-3 RF models is ranked in Table 3. The list with the top-15 predictors demonstrates the larger contribution of the more informative spatio-temporal variables (EMEP4UK PM 2.5 , meteorological parameters, BLH, and sea-level pressure). All the proposed spatially-lagged PM 2.5 variables were classified as highly important and their ranking positions varied slightly across the displayed years. This suggests the presence of spatial correlations in PM 2.5 values that are not entirely captured by the other variables.  Table 4 shows the results of the 10-Fold CV Stage-3 RF models by year. The results indicate a good predictive performance of the model throughout the study period. Overall cross-validated R 2 ranged from 0.704 (2008) to 0.821 (2011), with an average of 0.767. The average prediction error is 4.042 µg/m 3 , with a negligible bias in intercept and slope. The inspection of the spatial and temporal contributions shows that the model performs well generally across the two components, displaying a spatial performance drop only for 2008 (0.486) and 2015 (0.579). The cross-validated spatial R 2 ranges from 0.486 (2008) to 0.746 (2017), while the cross-validated temporal R 2 from 0.738 (2010) to 0.843 (2011). The two components have an average of 0.658 and 0.795, respectively. The high spatial R 2 performance across the years demonstrates that the Stage-3 RF models were able to predict the spatial variation of long-term PM 2.5 across Great Britain with good accuracy. In Supplementary Materials, Table S4 shows the Stage-3 results by season, demonstrating that the seasonal patterns were well described by RF models, although with a lower accuracy for the temporal domain in summer, characterized by a lower 10-Fold R 2 . This drop in performance during summer was also seen by other studies [17,20]. Table 4. Predicted-PM 2.5 concentrations obtained from Stage-3 RF models were regressed against Stage-1 measured/predicted-PM 2.5 concentrations in a linear regression model. The CV-R 2 (how well the model described the PM 2.5 variability in new locations) described in three different patterns (overall, spatial, and temporal), RMSE (µg/m 3 ), intercept (µg/m 3 ), and slope (µg/m 3 ).

Stage-3
Overall   Figure 3, revealing a decrease of pollution levels in recent years across the whole territory, although slightly stronger in England. Table S5 in Supplementary Materials provides the same figures for all the years, confirming the decreasing trend. The spatial comparison suggests that PM 2.5 concentrations are lower in Scotland and Wales compared to the more populated southern regions of England, with hotspots located in urban areas such as Liverpool, Manchester, Birmingham, and Greater London. At the bottom of Figure 3, the maps display the corresponding annual average of PM 2.5 levels in London, demonstrating the precision of the multi-stage ML model in reconstructing PM 2.5 concentrations in 1 km grid cells within urban areas. The maps show local hotspots of high pollution, with a spatial distribution that, however, changes along the study period.

Stage
The greatest contribution of this study is in the ability of the satellite-based ML models to reconstruct daily levels of PM 2.5 over a 1 km grid across a wide geographical domain. Figure 4 displays the time series plots of observed and predicted PM 2.5 concentrations at three monitoring sites and their locations within the geographical domain of Great Britain. While not necessarily representative, the results show that the error varies depending on location, type of monitor, and period. However, generally the plots indicate a very good performance of the ML algorithms in recovering the observed temporal variation in PM 2.5 levels, capturing peaks and periods of stable low concentrations.
years, confirming the decreasing trend. The spatial comparison suggests that PM2.5 concentrations are lower in Scotland and Wales compared to the more populated southern regions of England, with hotspots located in urban areas such as Liverpool, Manchester, Birmingham, and Greater London. At the bottom of Figure 3, the maps display the corresponding annual average of PM2.5 levels in London, demonstrating the precision of the multi-stage ML model in reconstructing PM2.5 concentrations in 1 km grid cells within urban areas. The maps show local hotspots of high pollution, with a spatial distribution that, however, changes along the study period. The greatest contribution of this study is in the ability of the satellite-based ML models to reconstruct daily levels of PM2.5 over a 1 km grid across a wide geographical domain. Figure 4 displays the time series plots of observed and predicted PM2.5 concentrations at three monitoring sites and their locations within the geographical domain of Great Britain. While not necessarily representative, the results show that the error varies depending on location, type of monitor, and period. However, generally the plots indicate a very good performance of the ML algorithms in recovering the observed temporal variation in PM2.5 levels, capturing peaks and periods of stable low concentrations. Figure 5 complements the analysis of daily variations by displaying the spatial distribution of PM2.5 estimations across Great Britain (Top) and London (Bottom) for specific days within the study period. It is interesting to note the wide variation in PM2.5 concentrations between days in the same area and between areas on different days. For instance, the maps in the left panels represent a day with almost no variation and generally low concentrations; the maps in the mid panels display a strong north-south split likely linked to weather conditions; the right panels show a more complex pattern with a wider range of PM2.5 values, a large area with very high pollution concentrations located in east London, and hotspots in highly populated urban areas across England.   Figure 5 complements the analysis of daily variations by displaying the spatial distribution of PM 2.5 estimations across Great Britain (Top) and London (Bottom) for specific days within the study period. It is interesting to note the wide variation in PM 2.5 concentrations between days in the same area and between areas on different days. For instance, the maps in the left panels represent a day with almost no variation and generally low concentrations; the maps in the mid panels display a strong north-south split likely linked to weather conditions; the right panels show a more complex pattern with a wider range of PM 2.5 values, a large area with very high pollution concentrations located in east London, and hotspots in highly populated urban areas across England.

Discussion
This study presents the first application of satellite-based spatio-temporal ML methods to reconstruct levels of pollution across Great Britain, providing estimates of daily PM2.5 concentrations over a 1 km grid during 2008-2018. The multi-stage ML framework provided significant advantages, allowing the combination of information from multiple data sources, such as air quality monitoring networks, remote sensing satellite products, chemical dispersion models, reanalysis databases, and administrative census data, among others.
The beginning of the study period (2008 and 2009) had a lower quantity of monitors measuring PM2.5 across Great Britain, but this number increased considerably from 2010 after a new measuring network had been established in 2009 [46]. Therefore, Stage-1 was an extremely important step to extend the number of PM2.5 measurements. The implementation of Stage-2 was also relevant to fill the gaps in MAIAC AOD retrievals from satellites, thus maximizing the available information. Stage-3 produced an ML prediction model based on a long list of informative predictors, accounting for potentially complex inter-relationships and functional forms. The Stage-3 ML algorithms offered excellent performance, showing an average cross-validated R 2 of 0.767 across the period, with an increased predictive ability in the last years. Stage-4 provided a single PM2.5 estimation for each of the 234.429 1 km grid cells in each of the 4018 days, totalling around 950 million data points. The methodology provides complete spatial coverage, high resolution, and a relative small error of the predictions, and the ability to capture variations in PM2.5 concentrations across both spatial and

Discussion
This study presents the first application of satellite-based spatio-temporal ML methods to reconstruct levels of pollution across Great Britain, providing estimates of daily PM 2.5 concentrations over a 1 km grid during 2008-2018. The multi-stage ML framework provided significant advantages, allowing the combination of information from multiple data sources, such as air quality monitoring networks, remote sensing satellite products, chemical dispersion models, reanalysis databases, and administrative census data, among others.
The beginning of the study period (2008 and 2009) had a lower quantity of monitors measuring PM 2.5 across Great Britain, but this number increased considerably from 2010 after a new measuring network had been established in 2009 [46]. Therefore, Stage-1 was an extremely important step to extend the number of PM 2.5 measurements. The implementation of Stage-2 was also relevant to fill the gaps in MAIAC AOD retrievals from satellites, thus maximizing the available information. Stage-3 produced an ML prediction model based on a long list of informative predictors, accounting for potentially complex inter-relationships and functional forms. The Stage-3 ML algorithms offered excellent performance, showing an average cross-validated R 2 of 0.767 across the period, with an increased predictive ability in the last years. Stage-4 provided a single PM 2.5 estimation for each of the 234.429 1 km grid cells in each of the 4018 days, totalling around 950 million data points. The methodology provides complete spatial coverage, high resolution, and a relative small error of the predictions, and the ability to capture variations in PM 2.5 concentrations across both spatial and temporal domains. The model offers a prediction accuracy that makes the output suitable for application in epidemiological studies on the short-and long-term health effects of air pollution. As demonstrated in Figure 4, the methodological approach presented in this study was able to capture the daily PM 2.5 variability across different years, locations, and monitor types.
The output from the empirical ML model developed in this study complements existing databases of modelled PM 2.5 in the United Kingdom, with some advantages for applications in epidemiological studies. Country-wide maps generated by emission-dispersion models are usually available at a coarser spatial or temporal resolution [32,47,48], and they generally show lower small-scale accuracy when tested against observed monitoring data [49,50]. The spatio-temporal ML models presented here demonstrated comparable predictive performance to similar methods applied in other countries, based either on single-learner ML models [13,20], ensemble ML models [18,19,21], or generalised additive models (GAM) [15]. Ref [21] developed an ensemble ML model (composed by RF, Deep Neural Networks, GAM, Gradient Boosting, K-nearest Neighbour) to estimate PM 2.5 for Greater London, reaching a mean 2005-2013 CV spatial-R 2 of 0.396. Using the 2008-2013 period, this study reached a CV spatial-R 2 of 0.637 for the whole of Great Britain. Modelling a larger area might have provided more information and a higher spatial variability, improving substantially the CV spatial-R 2 . The'application of a single learner to model air pollution in both spatial and temporal domains for Great Britain achieved a satisfactory performance. Nonetheless, ensemble model formats can be used in future applications using the same area of study and variables to assess how much performance gain is reached compared to the RF-only learner.
In 2008, a new Air Quality Directive came into force, bringing considerable changes to the following UK annual air quality assessments in 2010, setting an annual mean target of 25 µg/m 3 [51]. Several policy controls and emission reductions had been put in place aiming to reduce from 2010 the road traffic PM emission by 83%, off-road mobile machinery by 54%, and energy production by 32% until 2020 [52]. Independently of the temporal aggregation (daily or annual), all maps shown in this study detected this considerable drop in the PM 2.5 concentrations from 2010 across Great Britain.
Some limitations of this study must be acknowledged. First, the multi-stage model relies on the extension of the observed series of PM 2.5 by predicting values based on co-located PM 10 in Stage-1. This step was necessary for the application of the method in the early years characterized by sparse PM 2.5 monitoring, which is likely to have contributed to the lower predictive performance in this period. Morover, the cross-validation procedures revealed the presence of some bias in the Stage-1 predictions, particularly in the spatial domain, which was probably linked to limitations in modelling the relative distribution of the two PM components. Second, while the model displayed a good performance throughout the period, the accuracy is worse in the temporal domain in the summer (Table S4), suggesting limitations in capturing the higher temporal variation in this season. Third, the generalization of the prediction model was dependent on the selected locations of the monitors, which may be not representative of the study domain. This can result in an underestimation of the error and potential biases in the predictions in more remote and less represented areas if structural spatial differences are not entirely captured by the model covariates.
The findings of this study also highlight additional limitations not yet discussed in previous studies related to the use of remote sensing satellite products for reconstructing air pollution exposure for epidemiological purposes. The RF importance ranking order showed that Stage-3 was mainly informed by EMEP4UK PM 2.5 , regional and local interpolated PM 2.5 estimations, and meteorological variables, while the contribution from Stage-2 outputs (i.e., predicted-AOD 0.47 µm and 0.55 µm) was very limited. Therefore, direct satellite observations (i.e., AOD, NDVI, elevation, land cover, impervious surfaces, and night-time light) offered a relatively minor contribution to represent ground PM 2.5 exposures (Stage-3), while satellite-based products (i.e., indirect satellite observations) such as EMEP4UK and climate reanalysis products played a much more relevant role to derive reliable PM 2.5 estimates across Great Britain.
Future directions in remote sensing are pointing to new satellite instruments developed for air pollution monitoring which are likely to provide better resolution and reliability, thereby improving the predictive performance. As examples, the recently launched Copernicus Sentinel-5 Precursor [53] and the two future earth observation missions from European Space Agency (e.g., Copernicus Sentinel-4 [54] and Sentinel-5 [55]) were developed to provide information on atmospheric variables, such as air quality parameters (nitrogen dioxide, ozone, and aerosols). Finally, future research developments related to advanced statistical methodology, as demonstrated by the ML framework proposed here, include the application of geostatistical techniques, the use of alternative single-learner or ensemble ML algorithms, and statistical downscaling methods to increase further the resolution of the predictions.

Conclusions
This study developed and applied a multi-stage ML model, combining data for multiple sources, including remote sensing satellite products, climate and atmospheric reanalysis models, chemical transport models, and geospatial features, to generate a complete map of daily PM 2.5 concentrations in a 1 km grid across Great Britain between 2008 and 2018. The model showed good performance overall and in both spatial and temporal domains, with an accuracy that is compatible with the use of such reconstructed values as a proxy for PM 2.5 exposures in epidemiological studies. In particular, the availability of high-resolution measures that can be linked as such or aggregated at different spatial and/or temporal scales makes the output suitable for investigations on both transient and chronic health risks associated to short and long-term exposures to PM 2.5 , respectively.