Ian Hough

Ron Sarafian

Alexandra Shtein

Bin Zhou

Johanna Lepeule

Itai Kloog

Gaussian Markov random fields improve ensemble predictions of daily 1 km PM 2.5 and PM 10 across France

Keywords: Particulate matter, Exposure assessment, Aerosol optical depth, Ensemble model, Epidemiology

Understanding the health impacts of particulate matter (PM) requires spatiotemporally continuous exposure estimates. We developed a multi-stage ensemble model that estimates daily mean PM2.5 and PM10 at 1 km spatial resolution across France from 2000 to 2019. First, we alleviated the sparsity of PM2.5 monitors by imputing PM2.5 at more common PM10 monitors. We also imputed missing satellite aerosol optical depth (AOD) based on modelled AOD from atmospheric reanalyses. Next, we trained three base learners (mixed models, Gaussian Markov random fields, and random forests) to predict daily PM concentrations based on AOD, meteorology, and other variables. Finally, we generated ensemble predictions using a generalized additive model with spatiotemporally varying weights that exploit the strengths and weaknesses of each base learner. The Gaussian Markov random field dominated the ensemble, outperforming mixed models and random forests at most locations on most days. Rigorous cross-validation showed that the ensemble predictions were quite accurate, with mean absolute error (MAE) of 2.72 µg/m 3 and R 2 of 0.76 for PM2.5; PM10 MAE was 4.26 µg/m 3 and R 2 0.71. Our predictions are available to improve epidemiological studies of acute and chronic PM exposure in urban and rural France.

Introduction

Ambient particulate matter (PM) air pollution is a leading environmental health risk and the 7 th overall risk factor for death and disease worldwide [START_REF] Murray | Global burden of 87 risk factors in 204 countries and territories, 1990-2019: a systematic analysis for the Global Burden of Disease Study[END_REF]. To protect public health, the World Health Organization sets guideline limits [START_REF] Who | Air quality guidelines[END_REF] for both chronic (e.g. annual mean) and acute (e.g. daily mean) exposure to PM with a diameter of 10 µm or less (PM10) and 2.5 µm or less (PM2.5), and many countries maintain air quality monitoring networks to ensure compliance with national limits. But these networks are not always sufficient for epidemiological studies: monitors tend to be clustered in cities, which makes it difficult to estimate exposure for suburban and rural populations, and are often too sparse to capture spatial variation in PM concentrations within a city.

Recently, methods have been developed to improve exposure estimates for epidemiological studies by modelling the continuous spatiotemporal distribution of PM. One approach is to use chemical transport models that simulate the formation, dispersion, and deposition of PM based on emissions, meteorology, and atmospheric chemistry. These are effective at large and small scales and can be used to forecast future PM concentrations (Zhang et al., 2012a), but are limited by the accuracy of the input emissions and meteorological data and the completeness of their representation of atmospheric physics and chemistry (Zhang et al., 2012b). They are also computationally intensive, although advances in computing have allowed increasing resolutions to be used over larger areas and timescales. For example, a recent study estimated PM concentration in France in 2010 and 2011 by combining a national model at 4 km-hourly resolution, 7 regional models at 3 to 4 km resolution, and 43 urban models at 10 to 200 m resolution [START_REF] Riviere | Air pollution modeling and exposure assessment during pregnancy in the French Longitudinal Study of Children (ELFE)[END_REF].

An alternative approach is to calibrate a statistical model that relates measured PM concentration to variables such as aerosol optical depth (AOD), a measure of the absorption and scattering of light by particles suspended in the atmosphere. AOD can be retrieved by satellite instruments across large areas at fairly high spatiotemporal resolution (e.g. 1 km-daily), making it a useful proxy for the spatial and temporal distribution of ground-level PM. However, satellite AOD is often missing due to cloud cover, glint on snow or water, or instrument malfunction. Early studies accommodated this by calibrating two relationships: one to estimate PM based on AOD and another to estimate PM when AOD was not available [START_REF] Hu | Estimating ground-level PM2.5 concentrations in the Southeastern United States using MAIAC AOD retrievals and a twostage model[END_REF][START_REF] Kloog | Assessing temporally and spatially resolved PM2.5 exposures for epidemiological studies using satellite aerosol optical depth measurements[END_REF]. Recently, methods have been developed to fill gaps in satellite AOD based on modelled AOD from atmospheric reanalyses, allowing AOD-based prediction for all days and locations [START_REF] Di | An ensemble-based model of PM2.5 concentration across the contiguous United States with high spatiotemporal resolution[END_REF][START_REF] Di | Assessing PM2.5 Exposures with High Spatiotemporal Resolution across the Continental United States[END_REF][START_REF] Stafoggia | Estimation of daily PM10 and PM2.5 concentrations in Italy, 2013-2015, using a spatiotemporal land-use random-forest model[END_REF].

Another challenge is that the relationship between PM and AOD varies over both time and space.

Studies in the United States [START_REF] Chudnovsky | Fine particulate matter predictions using high resolution Aerosol Optical Depth (AOD) retrievals[END_REF][START_REF] Chudnovsky | Prediction of daily fine particulate matter concentrations using aerosol optical depth retrievals from the Geostationary Operational Environmental Satellite (GOES)[END_REF][START_REF] Hu | Estimating ground-level PM2.5 concentrations in the Southeastern United States using MAIAC AOD retrievals and a twostage model[END_REF][START_REF] Kloog | A new hybrid spatio-temporal model for estimating daily multi-year PM2.5 concentrations across northeastern USA using high resolution aerosol optical depth data[END_REF][START_REF] Kloog | Assessing temporally and spatially resolved PM2.5 exposures for epidemiological studies using satellite aerosol optical depth measurements[END_REF][START_REF] Lee | A novel calibration approach of MODIS AOD data to predict PM2.5 concentrations[END_REF][START_REF] Lee | Spatiotemporal prediction of fine particulate matter using high-resolution satellite images in the Southeastern US 2003-2011[END_REF], Europe [START_REF] Beloconi | Estimating urban PM10 and PM2.5 concentrations, based on synergistic MERIS/AATSR aerosol observations, land cover and morphology data[END_REF][START_REF] De Hoogh | Modelling daily PM2.5 concentrations at high spatio-temporal resolution across Switzerland[END_REF][START_REF] Nordio | Estimating spatio-temporal resolved PM10 aerosol mass concentrations using MODIS satellite data and land use regression over Lombardy, Italy[END_REF][START_REF] Stafoggia | Estimation of daily PM10 concentrations in Italy (2006-2012) using finely resolved satellite data, land use variables and meteorology[END_REF], China [START_REF] Liang | MAIAC-based longterm spatiotemporal trends of PM2.5 in Beijing, China[END_REF][START_REF] Xiao | Full-coverage high-resolution daily PM2.5 estimation using MAIAC AOD in the Yangtze River Delta of China[END_REF][START_REF] Xie | Daily Estimation of Ground-Level PM2.5 Concentrations over Beijing Using 3 km Resolution MODIS AOD[END_REF][START_REF] Zhang | National scale spatiotemporal land-use regression model for PM2.5, PM10 and NO2 concentration in China[END_REF][START_REF] Zheng | Estimating ground-level PM2.5 concentrations over three megalopolises in China using satellite-derived aerosol optical depth measurements[END_REF], Mexico [START_REF] Just | Using High-Resolution Satellite Aerosol Optical Depth To Estimate Daily PM2.5 Geographical Distribution in Mexico City[END_REF], and Israel [START_REF] Kloog | Estimating daily PM2.5 and PM10 across the complex geo-climate region of Israel using MAIAC satellite-based AOD data[END_REF][START_REF] Shtein | Estimating daily and intra-daily PM10 and PM2.5 in Israel using a spatio-temporal hybrid modeling approach[END_REF] have used mixed models to allow the PM -AOD relationship to vary from day to day and between regions. This approach performs well and is computationally cheap but implies sharp changes in the daily PM -AOD relationship at the borders of predefined regions, which may be unrealistic. Smooth spatiotemporal variation is possible for large datasets with a Gaussian Markov random field (GMRF) solved via integrated nested Laplace approximations (INLA) [START_REF] Lindgren | An explicit link between gaussian fields and gaussian markov random fields: The stochastic partial differential equation approach[END_REF][START_REF] Rue | Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations[END_REF]. Recently, [START_REF] Sarafian | Gaussian Markov Random Fields versus Linear Mixed Models for satellite-based PM2.5 assessment: Evidence from the Northeastern USA[END_REF] showed that GMRFs with spatially smooth daily random effects predicted daily PM2.5 more accurately than mixed models both near to and far from monitors in the northeastern United States. GMRFs have also performed well predicting daily PM10 in northwest Italy [START_REF] Cameletti | Spatio-temporal modeling of particulate matter concentration through the SPDE approach[END_REF] and annual mean PM across Europe [START_REF] Beloconi | Bayesian geostatistical modelling of PM10 and PM2.5 surface level concentrations in Europe using high-resolution satellite-derived products[END_REF].

Other studies have used various statistical approaches, including geographically weighted regression [START_REF] Hu | Estimating ground-level PM2.5 concentrations in the southeastern U.S. using geographically weighted regression[END_REF][START_REF] Ma | Estimating Ground-Level PM2.5 in China Using Satellite Remote Sensing[END_REF][START_REF] Song | A satellite-based geographically weighted regression model for regional PM2.5 estimation over the Pearl River Delta region in China[END_REF][START_REF] Van Donkelaar | Global Estimates of Fine Particulate Matter using a Combined Geophysical-Statistical Method with Information from Satellites, Models, and Monitors[END_REF], geographically and temporally weighted regression [START_REF] Guo | Estimating ground-level PM2.5 concentrations in Beijing using a satellite-based geographically and temporally weighted regression model[END_REF][START_REF] He | Satellite-based mapping of daily high-resolution ground PM2.5 in China via space-time regression modeling[END_REF]Liu et al., 2020), and machine learning algorithms such as random forests [START_REF] Chen | A comparison of linear regression, regularization, and machine learning algorithms to develop Europe-wide spatial models of fine particles and nitrogen dioxide[END_REF][START_REF] Hu | Estimating PM2.5 Concentrations in the Conterminous United States Using the Random Forest Approach[END_REF][START_REF] Schneider | A satellite-based spatio-temporal machine learning model to reconstruct daily PM2.5 concentrations across great britain[END_REF][START_REF] Stafoggia | A random forest approach to estimate daily particulate matter, nitrogen dioxide, and ozone at fine spatial resolution in Sweden[END_REF][START_REF] Stafoggia | Estimation of daily PM10 and PM2.5 concentrations in Italy, 2013-2015, using a spatiotemporal land-use random-forest model[END_REF], gradient boosting [START_REF] Chen | A comparison of linear regression, regularization, and machine learning algorithms to develop Europe-wide spatial models of fine particles and nitrogen dioxide[END_REF][START_REF] Just | Advancing methodologies for applying machine learning and evaluating spatiotemporal models of fine particulate matter (PM2.5) using satellite data over large regions[END_REF], and neural networks [START_REF] Chen | A comparison of linear regression, regularization, and machine learning algorithms to develop Europe-wide spatial models of fine particles and nitrogen dioxide[END_REF][START_REF] Di | Assessing PM2.5 Exposures with High Spatiotemporal Resolution across the Continental United States[END_REF][START_REF] Park | Estimating PM2.5 concentration of the conterminous United States via interpretable convolutional neural networks[END_REF][START_REF] Yan | New interpretable deep learning model to monitor real-time PM2.5 concentrations from satellite data[END_REF]. Machine learning algorithms perform particularly well as they can capture complex nonlinear relationships, and recent work has shown that performance can be slightly improved by ensembling predictions from multiple base learners [START_REF] Di | An ensemble-based model of PM2.5 concentration across the contiguous United States with high spatiotemporal resolution[END_REF][START_REF] Murray | A Bayesian ensemble approach to combine PM2.5 estimates from statistical models using satellite imagery and numerical model simulation[END_REF]Shtein et al., 2019;[START_REF] Zhai | Development of a stacked ensemble model for forecasting and analyzing daily average PM2.5 concentrations in Beijing, China[END_REF]. However, the flexibility of machine learning algorithms makes them vulnerable to overfitting, so it is important to evaluate their accuracy on independent data to ensure they can generalize to unmonitored locations [START_REF] Just | Advancing methodologies for applying machine learning and evaluating spatiotemporal models of fine particulate matter (PM2.5) using satellite data over large regions[END_REF][START_REF] Sarafian | Gaussian Markov Random Fields versus Linear Mixed Models for satellite-based PM2.5 assessment: Evidence from the Northeastern USA[END_REF]. A common approach is cross-validation (CV): data are repeatedly split into training and test sets, the model is calibrated using only the training data, and its predictions are compared to the held-out test data. Since PM concentrations are often spatiotemporally autocorrelated, the splitting must be done in a way that ensures test data are far in space and time from training data. Recent studies have done this using spatial blocking, holding out all data from a group of monitors (Just et al., 2020;Meng et al., 2021;[START_REF] Murray | A Bayesian ensemble approach to combine PM2.5 estimates from statistical models using satellite imagery and numerical model simulation[END_REF][START_REF] Park | Estimating PM2.5 concentration of the conterminous United States via interpretable convolutional neural networks[END_REF]Pu and Yoo, 2021;[START_REF] Schneider | A satellite-based spatio-temporal machine learning model to reconstruct daily PM2.5 concentrations across great britain[END_REF][START_REF] Shtein | Estimating Daily PM2.5 and PM10 over Italy Using an Ensemble Model[END_REF][START_REF] Stafoggia | Estimation of daily PM10 and PM2.5 concentrations in Italy, 2013-2015, using a spatiotemporal land-use random-forest model[END_REF]Xiao et al., 2020) or temporal blocking, holding out all data from one year (He et al., 2021(He et al., , 2020;;Meng et al., 2021;Pu and Yoo, 2021;Xiao et al., 2020;[START_REF] Yan | New interpretable deep learning model to monitor real-time PM2.5 concentrations from satellite data[END_REF]. Ensemble models require special care: data held out to test the ensemble should not have been used to train the base learners, and ensembles should be calibrated using base learner CV predictions (predictions for held-out test data), as these reflect each base learner's ability to generalize to new areas [START_REF] Shtein | Estimating Daily PM2.5 and PM10 over Italy Using an Ensemble Model[END_REF].

The goal of this study was to estimate PM2.5 and PM10 concentrations across continental France from 2000 through 2019. To the best of our knowledge, this is the first study to estimate daily 1 km PM concentration in France over two decades using satellite data and geostatistical models.

We incorporated recent methodological advances to mitigate the sparsity of PM2.5 monitors and fill gaps in satellite AOD [START_REF] Stafoggia | Estimation of daily PM10 and PM2.5 concentrations in Italy, 2013-2015, using a spatiotemporal land-use random-forest model[END_REF], and developed the first ensemble model of daily PM concentration that incorporates a GMRF. We used a rigorous cross-validation scheme to estimate accuracy and provide the first evidence that GMRFs may predict daily PM concentration more accurately than random forests.

Materials

Study domain

Continental France covers a roughly hexagonal area of 542,973 km 2 in western Europe bounded by the Atlantic Ocean to the west and the Mediterranean Sea to the southeast (Figure 1). Most of the terrain is at low elevation, but the Pyrenees in the southwest rise to over 3000 m and in the southeast the Alps reach 4,809 m. Annual mean temperature ranges from about 0 °C at high elevations to about 17 °C in the Mediterranean southeast [START_REF] Hough | A multi-resolution air temperature model for France from MODIS and Landsat thermal data[END_REF]. The population is approximately 64.5 million, of which 12.5 million (20%) live in the greater Paris metropolis. About 20% of the population is rural, and 37% live in towns and small cities with fewer than 500,000 residents [START_REF] Insee | Estimation de la population au 1 er janvier 2020[END_REF]. For this study, we defined a grid of 632,571 approximately 1 km 2 cells covering continental France coincident with the pixels of the satellite AOD data (section 2.3). We considered the 7,245 days from 1 March 2000 through 31 December 2019, giving a total study domain of 4.58×10 9 cell-days.

Air quality monitoring data

We obtained hourly PM2.5 and PM10 (µg/m 3 ) measurements from 12 regional air quality monitoring networks (federated by ATMO France) through the French Central Air Quality Monitoring Laboratory. Monitors were mostly clustered in urban areas; the number of PM10 monitors increased from 222 to 330 and the number of PM2.5 monitors increased from 9 to 142 over the course of the study period (Figure 1). To limit the impact of instrument malfunctions and rare events, we excluded hourly PM2.5 concentrations > 200 µg/m 3 and hourly PM10 concentrations > 300 µg/m 3 (0.003% of all observations). We indexed each monitor to the containing 1 km grid cell and calculated daily mean PM for days with at least 18 hourly observations.

Aerosol optical depth

We obtained satellite-derived 0.469 µm AOD at approximately 1 km spatial resolution from the Moderate Resolution Imaging Spectroradiometer (MODIS) Multi-Angle Implementation of Atmospheric Correction (MAIAC) AOD product (MCD19A2v006) [START_REF] Lyapustin | MODIS Collection 6 MAIAC algorithm[END_REF].

MCD19A2v006 provides AOD up to four times per day (between 9:00 and 15:00 UTC in France).

We used the quality assurance band to identify all "best quality" observations; we also included "land; research quality" and "clear; within 2 km of coast" as these represent potentially usable observations over urban areas and coasts where there are few "best quality" observations. We indexed these observations to the 1 km grid (whose cells were defined to coincide with the MAIAC AOD pixels) and calculated daily mean AOD.

To fill gaps in MAIAC AOD (mostly due to cloud cover), we obtained modelled 3-hourly 0.469 µm AOD at approximately 80 km spatial resolution from the Copernicus Atmospheric Monitoring Service EAC4 Reanalysis [START_REF] Inness | The CAMS reanalysis of atmospheric composition[END_REF]. Since EAC4 begins on 1 March 2003, for 1 March 2000 to 28 February 2003 we obtained modelled hourly 0.55 µm AOD for 08:30 to 15:30 UTC at approximately 60 km spatial resolution from the MERRA2 reanalysis [START_REF] Randles | The MERRA-2 aerosol reanalysis, 1980 onward. Part I: System description and data assimilation evaluation[END_REF]. We bilinearly interpolated EAC4 and MERRA2 AOD to the 1 km grid, giving 8 values per cell-day (0 UTC, 3 UTC, …, 21 UTC for EAC4; 08:30 UTC, 09:30 UTC, …, 15:30 UTC for MERRA2).

Meteorology

Meteorological parameters such as wind, rain, temperature, and the height of the planetary boundary layer affect surface PM concentrations and indicate the extent to which AOD represents aerosols near the surface or higher in the atmosphere. We obtained hourly meteorological parameters at approximately 30 km spatial resolution from the Copernicus Climate Change Service ERA5 reanalysis [START_REF] Hersbach | The ERA5 global reanalysis[END_REF]. We bilinearly interpolated the parameters to the 1 km grid and calculated 10 daily values: boundary layer height at 0:00 and 12:00 UTC, total precipitation, mean and standard deviation of 2m air temperature, mean 2m dewpoint temperature, mean surface pressure, mean u-and v-components of 10m wind speed, and mean cloud cover.

Normalized difference vegetation index

Vegetation may influence PM dispersion and the density of PM sources. We obtained monthly composite normalized difference vegetation index (NDVI) at approximately 1 km 2 spatial resolution from the MODIS MOD13A3v006 product [START_REF] Didan | MODIS Vegetation Index User's Guide (MOD13 Series )[END_REF], which is spatially coincident with the MAIAC AOD data. We indexed NDVI to the 1 km grid, filled rare missing values with the Gaussian kernel mean of nearby cells, and used the same value for every day of each month.

Spatial predictors

In addition to the previous spatiotemporal predictors, we used impervious surfaces, land cover, road and railway density, elevation, population, climatic region, distance to coast, and PM2.5 and PM10 emissions as indicators of the typical spatial distribution of PM. Since these data are time invariant, we used the value from the closest reference year for every day of each year.

Supplemental Table S1 describes how we indexed these data to the 1 km grid and derived 19 spatial predictors.

Methods

We used a four-stage process to predict PM2.5 and PM10 for the 4.58×10 9 1 km grid cell-days in the study domain (Figure 2). Briefly, we: 1) alleviated the sparsity of PM2.5 monitors by training a random forest (RF) to impute daily PM2.5 at monitors that only measured PM10; 2) filled gaps in MAIAC AOD data by training monthly RFs to impute missing MAIAC AOD based on co-located EAC4 or MERRA2 AOD; 3) trained three base learners for each year (linear mixed models [LMM],

Gaussian Markov random fields [GMRF], and random forests [RF]) to predict daily 1 km PM based on gap-filled AOD, meteorology, NDVI, and spatial predictors; 4) increased accuracy by ensembling the base learner predictions with annual generalized additive models (GAM) that weight the base learners according to spatiotemporal variations in their performance. We performed all data processing and statistical analyses in R 3.6.3 (R Core Team, 2020) using the packages lme4 for LMM [START_REF] Bates | Fitting Linear Mixed-Effects Models Using lme4[END_REF], R-INLA for GMRF [START_REF] Bakka | Spatial modeling with R-INLA: A review[END_REF], ranger for RF [START_REF] Wright | Ranger: A fast implementation of random forests for high dimensional data in C++ and R[END_REF] with mlr and mlrMBO for tuning via model-based optimization [START_REF] Bischl | mlrMBO: A Modular Framework for Model-Based Optimization of Expensive Black-Box Functions[END_REF][START_REF] Bischl | Mlr: Machine learning in R[END_REF], and mgcv for GAM [START_REF] Wood | Generalized Additive Models: An Introduction with R[END_REF].

Figure 2. Flowchart of four-stage process to predict daily 1 km PM2.5.

Stage 1: imputing PM2.5 at PM10 monitors

Most monitors in France measured PM10 but not PM2.5. To mitigate the sparsity of PM2.5 monitors, we used all co-located daily measures of PM2.5 and PM10 (n = 474,761) to tune and train a RF of 500 trees to predict PM2.5 based on measured PM10 and monitor characteristics: is the error at monitor m on day t. To reduce bias in the variable importance estimates, we sampled 63.2% of observations without replacement for each tree and estimated importance by permutation. [START_REF] Strobl | Bias in random forest variable importance measures: Illustrations, sources and a solution[END_REF] We tuned mtry (the number of variables to consider at each split)

𝑃𝑀 2.5 𝑚𝑡 =
to minimize mean absolute error via model-based optimization and estimated accuracy using 5fold CV blocked by monitor (section 3.5). We used the RF to impute PM2.5 for the 1.71×10 6 monitor-days where only PM10 was measured.

Stage 2: filling gaps in MAIAC AOD

Clouds and snow cover often prevented MAIAC AOD retrieval over part of the study area. To fill these gaps, we trained RFs to predict MAIAC AOD based on co-located modelled AOD from atmospheric reanalysis. For computational reasons, we used 96 trees per forest, tuned using one spatiotemporally blocked 50% subsample of the data (section 3.5), and trained one RF for each month in the study period (mean observations per month ≈ 4.36×10 6 ): and εst is the error at cell s on day t. We estimated accuracy using 5-fold CV with spatiotemporal blocking (see section 3.5) and used the RFs to predict AOD for the 3.54×10 9 1 km grid cell-days without MAIAC AOD.

𝐴𝑂𝐷 𝑠𝑡 𝑀 = 𝑓(𝑅

Stage 3: predicting daily 1 km PM using three base learners

In stage 3, we trained LMMs, GMRFs, and RFs to predict PM2.5 (from stage 1) and PM10 based on gap-filled MAIAC AOD (from stage 2), 11 spatiotemporal predictors (sections 2.4 and 2.5), and 19 spatial predictors (section 2.6). We scaled the predictors to have similar range, and for the LMMs and GMRFs we log-transformed PM to approximate normality and prevent negative predictions.

We trained each base learner for each of PM2.5 and PM10 in each year, yielding 120 base models (mean observations per year = 111,610; range 54,353 to 126,544). We estimated the accuracy of each base model using multi-stage CV blocked by monitor (section 3.5) and used the base models to predict PM for the 4.58×10 9 1 km grid cell-days of the study domain.

Linear mixed models

For each PM size fraction in each year, we calibrated a LMM with a random effect that allowed the is the spatiotemporal random effect at cell s on day t. We assumed that the error was independent and identically distributed (i.i.d.) following 𝒩(0, 𝜎 𝜀 2 ) and the spatiotemporal random effect was temporally i.i.d. with Matérn spatial covariance, i.e.:

𝐶𝑜𝑣(𝜔 𝑠𝑡 , 𝜔 𝑠 ′ 𝑡 ′ ) = { 0, 𝑡 ≠ 𝑡′ 𝜎 𝜔 2 𝒞(𝑑 𝑠𝑠 ′ ; 𝜌 𝜔 ), 𝑡 = 𝑡′ (5)
where 𝜎 𝜔 2 is the variance of the spatiotemporal random effect, 𝒞 is the Matérn function, 𝑑 𝑠𝑠 ′ is the Euclidean distance between locations s and s', and 𝜌 𝜔 is a hyperparameter that governs the range (distance at which the correlation falls to less than about 10%). We assigned penalized complexity priors to 𝜎 𝜀 2 , 𝜎 𝜔 2 , and 𝜌 𝜔 that shrank the spatiotemporal random effect towards the null [START_REF] Fuglstad | Constructing Priors that Penalize the Complexity of Gaussian Random Fields[END_REF][START_REF] Simpson | Penalising model component complexity: A principled, practical approach to constructing priors[END_REF] and fit the model using INLA.

Random forests

For each PM size fraction in each year, we trained a RF with the equation: and ys are the spatial coordinates of cell s; and wdayt and ydayt are the day of week and day of year. We used 250 trees and fixed mtry at 5 because exploratory tuning suggested that mtry > 5 provided little benefit and risked overfitting.

𝑃𝑀 𝑠𝑡 𝐹𝑌 = 𝑓(𝐴𝑂𝐷 𝑠𝑡 , 𝑋

Stage 4: ensembling daily 1 km PM predictions to improve accuracy

In stage 4, we calibrated a GAM to ensemble the predictions of the stage 3 base learners. We used predictions for held-out monitors to calibrate the GAMs because these reflect accuracy at unmonitored locations (section 3.5). We fit a GAM for each PM size fraction in each year (20 GAMs total; mean observations per year = 111,610; range 54,353 to 126,544) using tensor product smooths that allowed the coefficient for each base learner's predictions to vary smoothly over space and time: data from the fold that contains cell s; and 𝜀 𝑠𝑡 𝐹𝑌 is the error at cell s on day t. We estimated accuracy using multi-stage CV blocked by monitor (section 3.5) and used the GAMs to predict AOD for the 4.58×10 9 1 km grid cell-days of the study domain.

𝑃𝑀 𝑠𝑡 𝐹𝑌 =

Cross-validation

To limit bias due to spatiotemporal autocorrelation and avoid information leakage between the base learners and the ensemble, we extended the blocked CV scheme described by [START_REF] Shtein | Estimating Daily PM2.5 and PM10 over Italy Using an Ensemble Model[END_REF].

In stage 1 (imputing PM2.5 at PM10 monitors), we estimated accuracy at monitors that never measured PM2.5 using 5-fold CV blocked by monitor: we randomly assigned each monitor that measured both PM2.5 and PM10 to one of 5 folds. This ensured that no observations from test monitors were in the training set.

In stage 2 (filling gaps in MAIAC AOD data), we estimated accuracy at locations far from same-day MAIAC AOD observations (because MAIAC AOD tends to be missing in spatial clumps). We used 5fold CV with spatiotemporal blocking: we split the study area into 50 regions and randomly assigned MAIAC AOD in each day-region to one of 5 folds (Figure S1). This ensured that no sameday observations from test regions were in the training set.

In stages 3 and 4 (predicting daily 1 km PM with three base learners and ensembling the predictions), we estimated accuracy at unmonitored locations. We also trained the stage 4 ensembles on base learner predictions for held-out monitors (to emulate performance at unmonitored locations) and ensured that data held out to test the ensemble had not been used to train the base learners. We used multi-stage CV blocked by monitor (Figure S2). In stage 3, we randomly assigned each monitor to one of 5 folds, and for each combination of three of the folds (10 combinations total) we trained the base learners and predicted for the two held-out test folds.

In stage 4, we used the test predictions from all base learners that held out the same one of the five folds to train an ensemble and predict for the held-out fold. This ensured that the ensemble was trained on base learner predictions for held-out monitors and evaluated at monitors that were held out from both the base learners and the ensemble.

Performance metrics

We evaluated the models using mean absolute error (MAE), which reflects the typical difference between a model's predictions and measured PM, R 2 , which reflects the fraction of spatiotemporal variation in PM captured by a model, and root mean squared error (RMSE), which can be compared with the standard deviation (SD) of measured PM to see by how much a model improves upon a naïve prediction of the mean. We also split each of these metrics into a spatial and temporal component as described by [START_REF] Kloog | Assessing temporally and spatially resolved PM2.5 exposures for epidemiological studies using satellite aerosol optical depth measurements[END_REF].

Results

4.1. Stage 1: imputing PM2.5 at PM10 monitors Mean PM2.5 was 13.7 µg/m 3 (standard deviation [SD] 9.9 µg/m 3 ) and mean PM10 was 21.6 µg/m 3 (SD 12.5 µg/m 3 ) across all monitors in continental France for 2000-2019 (Figure S3). PM concentrations declined over the study period and were generally highest in winter and lowest in summer (Figure S4). The cross-validated predictions of the stage 1 RF showed good correspondence with observed PM2.5 at monitors that were not used to train the RF (R 2 = 0.87, MAE = 2.48 µg/m 3 ) with little bias (mean error = -0.157) but a tendency to underestimate very high concentrations (Figure 3). Performance was good even in early years when there were few PM2.5 monitors (R 2 >= 0.82 in every year except 2008; MAE < 2.5 µg/m 3 in most years) (Table S2).

The drop in performance in 2007 and 2008 coincided with a change in monitor technology that increased measured PM10 concentrations and likely complicated the relationship between PM10 and PM2.5. PM10 concentration was by far the most important predictor of PM2.5 concentration (Figure S5).

Stage 2: filling gaps in MAIAC AOD

Mean MAIAC AOD over continental France for 2000-2019 was 0.126 (SD 0.084); MAIAC AOD was missing for 77% of the 1 km cell-days in the study domain (Figure S6), similar to other areas [START_REF] Di | An ensemble-based model of PM2.5 concentration across the contiguous United States with high spatiotemporal resolution[END_REF][START_REF] Schneider | A satellite-based spatio-temporal machine learning model to reconstruct daily PM2.5 concentrations across great britain[END_REF][START_REF] Stafoggia | Estimation of daily PM10 and PM2.5 concentrations in Italy, 2013-2015, using a spatiotemporal land-use random-forest model[END_REF]. Cross-validated R 2 for the stage 2 RFs typically ranged from about 0.55 in winter to about 0.78 in summer (Figure S7), coinciding with fewer MAIAC observations in winter and more in summer. MAE typically ranged from about 0.025 in fall to about 0.034 in summer, coinciding with lower AOD in fall and higher AOD in summer.

Performance was similar between periods that used modelled AOD from MERRA2 vs. EAC4 as predictors. There was a slight tendency to overestimate high AOD (slope = 0.94), but average performance was good (mean R 2 = 0.70; mean MAE = 0.030) (Table S3). Prior to 2003, the most important predictors of MAIAC AOD were modelled AOD at 12 UTC, the day of year, and the spatial y coordinate (Figure S8). From 2003 on, modelled AOD at 15 UTC was the second most important predictor. This likely related to the mid-2002 launch of the Aqua satellite, which passes over continental France around 13 UTC; previously, MAIAC AOD was only available around 11 UTC from the Terra satellite. Figure S9 shows an example of gapfilled AOD.

4.3. Stages 3 and 4: predicting daily 1 km PM10 and PM2.5 with three base learners and ensembling the predictions to increase accuracy

Table 1 shows the average cross-validated performance of the stage 3 base learners and stage 4 GAM ensemble. GMRF was the most accurate base learner (mean PM2.5 R 2 = 0.75, MAE = 2.72 µg/m 3 ; mean PM10 R 2 = 0.70, MAE = 4.26 µg/m 3 ), followed by RF, with LMM the least accurate. The stage 4 GAM ensembles slightly improved performance (mean PM2.5 R 2 = 0.76, MAE = 2.72 µg/m 3 ; mean PM10 R 2 = 0.71, MAE = 4.26 µg/m 3 ), almost eliminating the GMRFs' slight bias and increasing spatial R 2 compared to both the GMRFs and RFs. The relative importance of the base learners in the GAM ensemble varied over space and time (Figure S10), but GMRF predictions usually had the highest weight, consistent with their high cross-validated accuracy. S11). The base learners and GAM ensemble captured day-to-day variation in PM concentration at individual locations better than between-locations differences in annual mean PM concentration (GAM ensemble spatial R 2 ≈ 0.46, temporal R 2 ≈ 0.80). This is in part because PM concentration varies more over time than space; spatial MAE was lower than temporal MAE (PM2.5 GAM ensemble spatial MAE = 1.6, temporal MAE = 2.2), indicating that predicted annual mean PM concentrations were quite accurate. It may also reflect difficulty capturing spatial variation in PM concentrations in urban areas. The lowest spatial R 2 and highest spatial MAE were in Île-de-France, the densely populated region that contains Paris, which also had the highest and most variable PM concentrations (Table S4, Figure S12). There may not have been enough monitors for the model to capture complex spatial variation in PM concentration over greater Paris.

Since the majority of our PM2.5 data consisted of imputed PM2.5 concentration at PM10 monitors from the stage 1 RF, we performed a sensitivity analysis comparing the cross-validated predictions of the GAM ensemble to only observed PM2.5 concentration at PM2.5 monitors. Apart from 2000

(when there were only 5 PM2.5 monitors), performance was similar at PM2.5 monitors (mean R 2 = 0.77, mean MAE = 2.98) and across all monitors (mean R 2 = 0.77, mean MAE = 2.71). We also constructed an alternate model by retraining the base learners and ensemble using only PM2.5 monitors for all years except 2000. This alternate model was less accurate (mean R 2 = 0.66, mean MAE = 3.62) than the stage 4 GAM ensemble evaluated only at PM2.5 monitors (Figure 5), indicating that increasing the quantity of training data by imputing PM2.5 at PM10 monitors resulted in more accurate final PM2.5 predictions than if we had relied solely on PM2.5 monitors.

Discussion

Our finding that GMRFs predicted daily 1 km PM concentration more accurately than LMMs is consistent with results in the northeastern United States [START_REF] Sarafian | Gaussian Markov Random Fields versus Linear Mixed Models for satellite-based PM2.5 assessment: Evidence from the Northeastern USA[END_REF]. Our finding that GMRFs were also more accurate than RFs is novel and of note, as RFs and other tree-based machine learning algorithms performed well in several recent studies [START_REF] Di | An ensemble-based model of PM2.5 concentration across the contiguous United States with high spatiotemporal resolution[END_REF][START_REF] Just | Advancing methodologies for applying machine learning and evaluating spatiotemporal models of fine particulate matter (PM2.5) using satellite data over large regions[END_REF][START_REF] Schneider | A satellite-based spatio-temporal machine learning model to reconstruct daily PM2.5 concentrations across great britain[END_REF][START_REF] Stafoggia | A random forest approach to estimate daily particulate matter, nitrogen dioxide, and ozone at fine spatial resolution in Sweden[END_REF][START_REF] Stafoggia | Estimation of daily PM10 and PM2.5 concentrations in Italy, 2013-2015, using a spatiotemporal land-use random-forest model[END_REF]. We emphasize the importance of careful performance evaluation when using flexible machine learning algorithms: GMRFs had the best cross-validated performance (corresponding to accuracy at unmonitored locations), but RFs performed better than GMRF on non-independent training data (corresponding to accuracy at monitors). Evaluation methods that do not ensure independence between training and testing data risk mistaking good performance at monitors for good performance everywhere.

Our GAM ensemble captured temporal variation in PM concentration better than spatial variation.

We attempted to improve spatial performance in urban areas by downscaling the residuals of the GAM ensemble using RFs trained on the high spatial resolution predictors listed in Table S1, as was done in a few previous studies [START_REF] Di | An ensemble-based model of PM2.5 concentration across the contiguous United States with high spatiotemporal resolution[END_REF][START_REF] Kloog | A new hybrid spatio-temporal model for estimating daily multi-year PM2.5 concentrations across northeastern USA using high resolution aerosol optical depth data[END_REF][START_REF] Stafoggia | Estimation of daily PM10 and PM2.5 concentrations in Italy, 2013-2015, using a spatiotemporal land-use random-forest model[END_REF]. Unlike previous studies, we used 5-fold cross-validation blocked by monitor to assess whether the downscaling improved accuracy: the downscaled predictions were less accurate than the 1 km GAM ensemble predictions (higher MAE, lower R 2 ). Downscaling over cities is an area for further research, as epidemiological studies would benefit from better estimates of differences in PM exposure within a city.

Despite good overall performance, our approach has some limitations. First, the sparsity of the monitoring network limited performance in early years. Even in later years, most monitors were clustered in cities, making it difficult to evaluate accuracy in smaller towns and rural areas and risking overreliance on predictors that work well in urban areas but may not work well elsewhere.

The clustering of monitors near cities means our model is roughly weighted by population density, which may or may not be appropriate depending on the intended use for the predictions [START_REF] Sarafian | A Domain Adaptation Approach for Performance Estimation of Spatial Predictions[END_REF]. New low-cost PM sensors might complement the existing monitoring network, particularly since our model's weaker spatial performance suggests that a few PM observations at new locations might be more useful than a long timeseries of measurements at a single location.

Second, MAIAC AOD was the only predictor with both a high spatial (1 km) and temporal (daily)

resolution, but since it is based on a few daytime observations of the entire atmospheric column, it is both vertically and temporally misaligned with surface-level daily mean PM concentration. We included planetary boundary layer height to help distinguish between surface-level vs. highaltitude aerosols, but it had a much coarser spatial resolution than MAIAC AOD. Our model might have benefitted from AOD at coarser spatial but higher temporal resolution, such as from geostationary weather satellites, or from considering longer time periods and giving greater weight to rare observations when filling gaps in MAIAC AOD.

Despite these limitations, our multi-stage ensemble approach was able to predict daily 1 km PM2.5 and PM10 with low error across a large area over 20 years. To the best of our knowledge, this is the first work conducted in France with such a high spatiotemporal resolution (1 km-daily), large spatial extent (national) and long temporal coverage (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019). We increased accuracy by supplementing sparse PM2.5 observations with imputed data and by ensembling the predictions of three base learners. We confirmed that Gaussian Markov random fields predict daily PM concentration better than linear mixed models and provide the first evidence that that may also outperform random forests. Our dataset of daily 1 km PM2.5 and PM10 is available to health and ecosystems researchers in France and may inform policy makers on air quality issues.

Figure 1 .

 1 Figure 1. Spatial distribution of PM2.5 monitors (top row) and PM10 monitors (bottom row) in

Figure 3 .

 3 Figure 3. Cross-validation (CV) predicted vs. observed daily PM2.5 concentrations from the stage 1

Figure 4 .

 4 Figure 4. Annual cross-validated R 2 (top) and MAE (bottom; µg/m 3 ) of the stage 3 base learners

Figure 5 .

 5 Figure 5. Annual cross-validated performance at PM2.5 monitors of the stage 4 GAM ensemble (red

Figure 6

 6 Figure6shows the mean 2000-2019 PM2.5 and PM10 concentration predicted by each base learner

Figure 6 .

 6 Figure 6. Mean PM2.5 (top) and PM10 (bottom) concentration predicted by the stage 3 base

Figure 7

 7 Figure 7 shows PM2.5 and PM10 concentrations over greater Paris predicted by the GAM ensemble

Figure 7 .

 7 Figure 7. Mean 24-hour PM2.5 (top) and PM10 (bottom) concentration over greater Paris predicted

  𝑓(𝑃𝑀 10 𝑚𝑡 , 𝑣𝑜𝑙 𝑚𝑡 , 𝑙𝑜𝑐 𝑚𝑡 , 𝑖𝑛𝑓𝑙 𝑚𝑡 , 𝑙𝑎𝑡 𝑚 , 𝑙𝑜𝑛 𝑚 , 𝑤𝑑𝑎𝑦 𝑡 , 𝑦𝑑𝑎𝑦 𝑡 , 𝑑𝑎𝑡𝑒 𝑡 ) + 𝜀 𝑚𝑡[START_REF] Riviere | Air pollution modeling and exposure assessment during pregnancy in the French Longitudinal Study of Children (ELFE)[END_REF] where 𝑃𝑀 2.5 𝑚𝑡 and 𝑃𝑀 10 𝑚𝑡 are, respectively, the PM2.5 and PM10 measured by monitor m (1, …, 205) on day t(1, …, 7245); volmt indicates whether on day t monitor m excluded, included, or included an estimate of the semi-volatile fraction of PM10; locmt and inflmt are, respectively, the locale (rural, suburban, or urban) and predominant influence (traffic, industrial, or background) of monitor m on day t; latm and lonm are the latitude and longitude of monitor m; wdayt, ydayt, and datet are, respectively, the day of week (to capture trends related to commuting or business activity), day of year (to capture seasonal trends), and date (to capture long-term trends); and εmt

  1 𝑠𝑡 , … , 𝑅 8 𝑠𝑡 , 𝑥 𝑠 , 𝑦 𝑠 , 𝑤𝑑𝑎𝑦 𝑡 , 𝑦𝑑𝑎𝑦 𝑡 ) + 𝜀 𝑠𝑡 …, 15:30 UTC for MERRA2; 0 UTC, 3 UTC, …, 21 UTC for EAC4); xs and ys are the spatial coordinates of cell s; wdayt and ydayt are, respectively, the day of week and day of year;

	(2)
	where for each month M (1, …, 238), 𝐴𝑂𝐷 𝑠𝑡 𝑀 is the MAIAC AOD observed at 1 km grid cell s (1, …,

632571) on day t (1, …, number of days in month M); 𝑅 ℎ 𝑠𝑡 is the AOD from atmospheric reanalysis (MERRA2 before 1 January 2003; EAC4 otherwise) at cell s on day t at each of eight times (8:30 UTC, 9:30 UTC,

  𝛽 𝐴𝑂𝐷 𝐹𝑌 is the fixed slope of AOD and 𝜈 𝑡𝑟 𝐹𝑌 is the random slope of AOD on day t for the climatic region r that contains cell s; AODst is the AOD at cell s on day t. 𝛽 𝑝 𝐹𝑌 is the coefficient and 𝑋 𝑝 𝑠𝑡 the value at cell s on day t for each of the 11 spatiotemporal predictors, 19 spatial predictors, and sine and cosine transforms of the day of week; and 𝜀 𝑠𝑡 𝐹𝑌 is the error at cell s on day

	t.
	3.3.2. Gaussian Markov random fields
	For each PM size fraction in each year, we calibrated a GMRF with a spatiotemporal random effect
	that varied smoothly over space on each day:
	32
	log(𝑃𝑀 𝑠𝑡 𝐹𝑌 ) = 𝛼 𝐹𝑌 + 𝛽 𝐴𝑂𝐷 𝐹𝑌 𝐴𝑂𝐷 𝑠𝑡 + ∑ (𝛽 𝑝 𝐹𝑌 𝑋 𝑝 𝑠𝑡	)	+ 𝜔 𝑠𝑡 𝐹𝑌 + 𝜀 𝑠𝑡 𝐹𝑌	(4)
	𝑝=1
	PM-AOD relationship to vary daily for each of 8 climatic regions:
	32
	log(𝑃𝑀 𝑠𝑡 𝐹𝑌 ) = (𝛼 𝐹𝑌 + 𝜇 𝑡𝑟 𝐹𝑌 ) + (𝛽 𝐴𝑂𝐷 𝐹𝑌 + 𝜈 𝑡𝑟 𝐹𝑌 )𝐴𝑂𝐷 𝑠𝑡 + ∑ (𝛽 𝑝 𝐹𝑌 𝑋 𝑝 𝑠𝑡	)	+ 𝜀 𝑠𝑡 𝐹𝑌	(3)
	𝑝=1
	where for each PM size fraction F (2.5 or 10) and year Y (2000, …, 2019), log(𝑃𝑀 𝑠𝑡 𝐹𝑌 ) is the log-
	transformed PM concentration at 1 km grid cell s (1, …, 632571) on day t (1, …, number of days in

year Y); 𝛼 𝐹𝑌 is the fixed intercept and 𝜇 𝑡𝑟 𝐹𝑌 is the random intercept on day t for the climatic region r that contains cell s; where F, Y, s, t, log(𝑃𝑀 𝑠𝑡 𝐹𝑌 ), 𝛼 𝐹𝑌 , 𝛽 𝐴𝑂𝐷 𝐹𝑌 , 𝐴𝑂𝐷 𝑠𝑡 , 𝛽 𝑝 𝐹𝑌 , 𝑋 𝑝 𝑠𝑡 , and 𝜀 𝑠𝑡 𝐹𝑌 are as in equation (3), and 𝜔 𝑠𝑡 𝐹𝑌

  1 𝑠𝑡 , … , 𝑋 30 𝑠𝑡 , 𝑥 𝑠 , 𝑦 𝑠 , 𝑤𝑑𝑎𝑦 𝑡 , 𝑦𝑑𝑎𝑦 𝑡 ) + 𝜀 𝑠𝑡 𝑃𝑀 𝑠𝑡 𝐹𝑌 , 𝐴𝑂𝐷 𝑠𝑡 , and 𝜀 𝑠𝑡 𝐹𝑌 are as in equation (3); 𝑋 1 𝑠𝑡 , … , 𝑋 30 𝑠𝑡 are, respectively, the value for each of the 11 spatiotemporal predictors and 19 spatial predictors at cell s on day t; xs

	𝐹𝑌	(6)
	where F, Y, s, t,

  𝑡𝑒(𝑥 𝑠 , 𝑦 𝑠 , 𝑡) is the tensor product of penalized cubic regression splines of the spatial coordinates of cell s (xs and ys) and the temporal index t; 𝐿𝑀𝑀 𝑠𝑡 𝐹𝑌 , 𝐺𝑀𝑅𝐹 𝑠𝑡 𝐹𝑌 , and 𝑅𝐹 𝑠𝑡 𝐹𝑌 are, respectively, the CV prediction at cell s on day t from a LMM, GMRF, and RF that were trained while holding out all

	𝑡𝑒(𝑥 𝑠 , 𝑦 𝑠 , 𝑡)𝐿𝑀𝑀 𝑠𝑡 𝐹𝑌 + 𝑡𝑒(𝑥 𝑠 , 𝑦 𝑠 , 𝑡)𝐺𝑀𝑅𝐹 𝑠𝑡 𝐹𝑌 + 𝑡𝑒(𝑥 𝑠 , 𝑦 𝑠 , 𝑡)𝑅𝐹 𝑠𝑡 𝐹𝑌 + 𝜀 𝑠𝑡 𝐹𝑌	(7)

where for each PM size fraction F (2.5 or 10) and year Y (2000, …, 2019), 𝑃𝑀 𝑠𝑡 𝐹𝑌 is the PM concentration at 1 km grid cell s (1, …, 632571) on day t (1, …, number of days in year Y);

Table 1

 1 Slope of regression of CV predicted PM on observed PM R 2 for all models increased in early years with the number of monitors and remained high from 2009-2019; MAE covaried with mean observed PM (Figure 4). The sharp increase in PM10 MAE in 2007 coincided with a change in monitor technology: in 2007, all PM10 monitors were modified to measure semi-volatile particles in addition to non-volatile particles, increasing observed PM10 concentrations. PM2.5 monitors were modified in 2008 and 2009, corresponding to the increase in PM2.5 MAE in 2008 and 2009. R 2 was highest in winter and spring and lowest in summer; MAE was highest in winter and lowest in summer, corresponding to typically higher PM concentrations in winter and lower concentrations in summer (Figure

	. Cross-validated (CV) performance (averaged over 2000-2019) of the stage 3 base learners
	(LMM, GMRF, RF) and stage 4 ensemble (GAM) predicting daily 1 km PM (µg/m 3 ).
		Multi-stage CV Performance
	Observed PM	Total			Spatial	Temporal
	Mean SD a Model RMSE Bias b Slope c R²	MAE R²	MAE R²	MAE
	PM2.5 13.8 8.5 LMM 5.03	0.67 0.63 0.63 3.35	0.38 1.82	0.68 2.92
	GMRF 4.09	0.46 0.75 0.75 2.72	0.45 1.68	0.81 2.24
	RF	4.52 -0.14 0.63 0.70 3.18	0.47 1.69	0.74 2.68
	GAM 4.02 -0.01 0.76 0.76 2.72	0.49 1.63	0.81 2.23
	PM10 21.5 11.9 LMM 7.65	1.04 0.60 0.58 5.21	0.32 2.98	0.64 4.39
	GMRF 6.40	0.73 0.72 0.70 4.26	0.39 2.79	0.78 3.36
	RF	7.07 -0.19 0.57 0.64 5.00	0.41 2.80	0.70 4.11
	GAM 6.28 -0.02 0.72 0.71 4.26	0.43 2.71	0.78 3.34
	a Standard deviation
	b Mean error

c

Acknowledgements

Funding: This work was supported by the French National Centre for Scientific Research (ANR PRC 2018-2020), Fondation de France (no. 00081169), and Israeli Ministry of Science and Technology (PRC 2018(PRC -2020)). I.H. is supported by Univ. Grenoble Alpes (ANR-15-IDEX-02) and Ben Gurion University of the Negev. Most data processing and analyses were performed using the GRICAD infrastructure (https://gricad.univ-grenoble-alpes.fr), which is supported by Equip@Meso (ANR-10-EQPX-29-01). We thank the LCSQA for providing validated PM monitor data, INERIS for providing a spatial emissions inventory, and NASA, ECMWF, Copernicus, IGN, and

OpenTransportMap for their open datasets.