# Historical reconstruction of background air pollution over France for 2000–2015

Elsa Real, Florian Couvidat, Anthony Ung, Laure Malherbe, Blandine Raux,

Alicia Gressent, Augustin Colette

► To cite this version:

Elsa Real, Florian Couvidat, Anthony Ung, Laure Malherbe, Blandine Raux, et al. Historical reconstruction of background air pollution over France for 2000–2015. Earth System Science Data, 2022, 14 (5), pp.2419-2443. 10.5194/essd-14-2419-2022. ineris-03827420

HAL Id: ineris-03827420

https://ineris.hal.science/ineris-03827420

Submitted on 7 Nov 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Historical reconstruction of background air pollution over France for 2000–2015

Elsa Real, Florian Couvidat, Anthony Ung, Laure Malherbe, Blandine Raux, Alicia Gressent, and Augustin Colette

Atmospheric Modelling and Environmental Mapping Unit, Environmental Characterisation and Impacts on Human and Biodiversity Department, INERIS, France

**Correspondence:** Elsa Real (elsa.real@ineris.fr)

Received: 31 May 2021 – Discussion started: 13 July 2021

Revised: 17 March 2022 – Accepted: 18 March 2022 – Published: 24 May 2022

**Abstract.** This paper describes a 16-year dataset of air pollution concentrations and air quality indicators over France. Using a kriging method that combines background air quality measurements and modeling with the CHIMERE chemistry transport model, hourly concentrations of NO₂, O₃, PM₁₀ and PM₂.₅ are produced with a spatial resolution of about 4 km. Regulatory indicators (annual average, SOMO35 (sum of ozone means over 35 ppb), AOT40 (accumulated ozone exposure over a threshold of 40 ppb), etc.) are also calculated from these hourly data. The NO₂ and O₃ datasets cover the period 2000–2015, as well as the annual PM₁₀ data. Hourly PM₁₀ concentrations are not available from 2000 to 2007 due to known artifacts in PM₁₀ measurements. PM₂.₅ data are only available from 2009 onwards due to the limited number of measuring stations available before this date. The overall dataset was evaluated over all years by a cross-validation process against background stations (rural, sub-urban and urban) to take into account the data fusion between measurement and models in the method. The results are very good for PM₁₀, PM₂.₅ and O₃. They show an overestimation of NO₂ concentrations in rural areas, while NO₂ background values in urban areas are well represented. Maps of the main indicators are presented over several years, and trends are calculated. Finally, exposure and trends are calculated for the three main health-related indicators: annual averages of PM₂.₅, NO₂ and SOMO35. The DOI link for the dataset is https://doi.org/10.5281/zenodo.5043645 (Real et al., 2021). We hope that the publication of this open dataset will facilitate further studies on the impacts of air pollution.

## 1 Introduction

Air pollution is a major environmental risk for human health and ecosystems in Europe. Over the past decades the European Union (EU) has put in place several measures to reduce anthropogenic emissions of pollutants. In response to emissions reductions, concentrations of SO₂, NO₂ and particles measured over Europe have shown a clear decrease since 1990 (EEA, 2018; EMEP, 2016).

The evolution of O₃ trends is less clear despite the decrease in its precursors. The magnitude of high ozone episodes decreased, while annual average ozone levels measured at EMEP (European Monitoring and Evaluation Programme) stations were increasing in the 1990s and show a limited negative trend from 2002 onwards. As shown in

the Tropospheric Ozone Assessment Report (Tarasick et al., 2019), this feature is generally attributed to the changing global tropospheric ozone baseline for which further hemispheric control strategies are needed. The same conclusions could be drawn from the Malherbe et al. (2017) study, which focused on France, with significant reductions in NO₂ and particle concentrations and an increase in average O₃ offset by a slight decrease in peak O₃. Despite these reductions in emissions and pollutant concentrations (with the exception of the annual average O₃), a proportion of French citizens are still exposed to concentrations above the EU limit and target value, and air quality in EU remains one of the main reasons for premature deaths (IHME, 2013).
As a complement to observations (which provide only partial spatial information), accurate, high-spatial-resolution and up-to-date air pollution maps are important information for assessing air pollution trends and exposure. They should provide geographically detailed information on the concentrations of air pollutants over the whole territory. These maps serve as a basis for informing citizens, designing and stratifying monitoring networks, supporting policy strategies, and measuring their impact. They are also used to estimate population exposure to air pollutants, which is essential for epidemiological studies.

On a European scale, different mapping approaches have been used to produce maps of pollutant concentrations. These maps can be obtained by modeling using a regional chemistry transport model (CTM) that simulates the concentration of pollutants over Europe. However, these models cannot always be used over the whole of Europe with a high resolution and have some biases and limitations in spatial representativeness. Regression methods (Briggs et al., 2000; Beelen et al., 2009) are also used at different scales. These stochastic modeling techniques develop statistical associations between potential “predictor variables” (land use, emission sources, topography) and measured pollutant concentrations in order to predict concentrations at an unsampled site. Other frequently used techniques are kriging techniques. These geostatistical techniques are based on the assumption that the data are spatially autocorrelated and therefore take into account the distances between measurements and the spatial structure of the variable. Different types of kriging are used to map the concentrations of air pollutants. Over France, kriging methods combining information from a regional CTM (CHIMERE; Mailler et al., 2017) and observations are produced daily by the Prev'air operational forecasting and mapping system for air quality (Rouil et al., 2009). Since 2003 (for ozone) and 2005 (for PM₁₀), the concentration maps simulated for the day before in Prev'air have been corrected each morning using observations. The kriging technique used in Prev'Air has evolved over time, and PM₂.₅ and NO₂ concentrations are now also corrected for the day before. Today, a kriging of hourly observations with CHIMERE as external drift is applied to map NO₂ and O₃ concentrations. Since 2017, for the mapping of PM₁₀ and PM₂.₅ concentrations, the method used is an hourly kriging of PM₁₀ and PM₂.₅ data with CHIMERE as external drift. These choices are the results of successive studies that compared different kriging techniques (Malherbe and Ung, 2009, Beauchamp, 2015). A similar methodology was implemented for an earlier reconstruction of outdoor air pollution in Europe for the period 1989–2008 in Bentayeb et al. (2014). There are also ambient air pollution maps produced at European scale at 1 km resolution by the European Environment Agency but only for selected annual indicators and without consistency for multi-year reconstructions (Horálek et al., 2012, 2020). The Copernicus Atmosphere Monitoring Service has also produced European analyses since 2015, but

again there is no multi-year consistency as these European maps are produced on an annual basis with gradually improving methodologies (Marécal et al., 2015). At the global scale, the Global Burden of Disease also makes available air pollution exposure maps; a recent update of the methodology was presented in (Shaddick et al., 2017), but the resolution is 0.1° or about 10 km.

The purpose of this paper and the associated datasets is to present and provide mapped data of O₃, NO₂, PM₁₀ and PM₂.₅ concentrations at high spatial and temporal resolution and associated regulatory indicators covering the French metropolitan territory for the period 2000–2015 (2007–2015 and 2009–2015 for hourly concentrations of PM₁₀ and PM₂.₅). The same kriging technique as in the Prev'air system is used to combine modeled and observed concentrations. Hourly concentrations of PM₁₀, PM₂.₅, NO₂ and O₃ are produced and mapped over France, and these hourly data are then used to calculate and map European and French air quality standards.

## 2 Methods

Model outputs and measurements from the permanent monitoring network were combined by external drift kriging (Malherbe and Ung, 2009; Benmerad et al., 2017) to construct hourly concentration maps over France for a long period: 2000 to 2015. Details on the input data and methods used are described in the following paragraphs. From these corrected hourly concentration data, annual regulatory air quality maps are then constructed over France.

### 2.1 Monitoring data

Hourly measurements are extracted from validated reference datasets. For France, observations are extracted from the national air quality databases: BDQA (Base de Données de Qualité de l'Air) before 2013 and GEODAIR (https://www.lcsqa.org/fr/les-donnees-nationales-de-qualite-de-lair, last access: 19 May 2022) after 2013, as well as from the Airbase database (https://www.eea.europa.eu/themes/air/air-quality/map/airbase, last access: 19 May 2022) for other European countries from 2000 to 2012 and from air quality e-reporting (https://www.eea.europa.eu/data-and-maps/data/aqereporting-9, last access: 19 May 2022) from 2013 to 2015. All background monitoring data over the spatial domain are used in the kriging procedure except for stations with measurements above the 95 percentiles. This includes rural, suburban and urban stations but excludes industrial and traffic stations that are representative of very local concentrations and are difficult to reproduce in a national-scale mapping system. The number of background monitoring sites for each type of station and for each year is summarized in Table 1.

Until 1 January 2007, operational monitoring of PM₁₀ and PM₂.₅ was carried out in France by automatic measur-
<table><caption>**Table 1.** Number of background French monitoring sites for the years 2000 to 2015.</caption><thead><tr><th></th><th>2000</th><th>2001</th><th>2002</th><th>2003</th><th>2004</th><th>2005</th><th>2006</th><th>2007</th><th>2008</th><th>2009</th><th>2010</th><th>2011</th><th>2012</th><th>2013</th><th>2014</th><th>2015</th></tr></thead><tbody><tr><th>O<sub>3</sub></th><td>284</td><td>310</td><td>337</td><td>362</td><td>378</td><td>396</td><td>404</td><td>405</td><td>399</td><td>385</td><td>376</td><td>360</td><td>347</td><td>318</td><td>319</td><td>331</td></tr><tr><th>NO<sub>2</sub></th><td>274</td><td>290</td><td>299</td><td>322</td><td>337</td><td>353</td><td>353</td><td>350</td><td>352</td><td>337</td><td>334</td><td>316</td><td>299</td><td>284</td><td>282</td><td>300</td></tr><tr><th>PM<sub>10</sub></th><td>119</td><td>125</td><td>171</td><td>212</td><td>219</td><td>238</td><td>126</td><td>219</td><td>252</td><td>241</td><td>249</td><td>245</td><td>240</td><td>218</td><td>173</td><td>251</td></tr><tr><th>PM<sub>2.5</sub></th><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>62</td><td>69</td><td>74</td><td>84</td><td>89</td><td>90</td><td>105</td></tr></tbody></table>

ing systems of the TEOM (tapered element oscillating mi-crobalance; PM₁₀, PM₂.₅) or Beta (PM₁₀) type. However, compared to the reference method EN 12341 (gravimetry), these systems underestimate the concentrations of particles. This is a known artifact related to the loss of semi-volatile compounds. To correct PM₁₀ concentrations measured be-for 2007, a simple approach consists of applying a uni-form correcting factor over France. This method is not suit-able for correcting hourly or daily concentrations, but it has been shown to work well for annual average PM₁₀ concen-trations (Malherbe et al., 2017, Bessagnet et al., 2008). The factor (1.36) is a median value calculated on the PM₁₀ data from “reference” sites (Bessagnet et al., 2008). As a conse-quence, for the period 2000 to 2006, the only PM₁₀ indica-tor available is the annual average concentration. Concern-ing PM₂.₅, given the few reference measurements available before 2009, the reliability of even annual measurements is low. It was therefore decided to apply the kriging method-ology only from the year 2009 onwards, which is when the change in measurement method had become widespread.

## 2.2 CHIMERE simulations

The CHIMERE chemistry transport model (Couvidat et al., 2018) is used to estimate air pollution levels for metropoli-tan France with a resolution of about 4 km (0.06° × 0.03°) over the years 2000 to 2015. This model has long been im-plemented and assessed in France as the main component of the national air quality forecasting and monitoring system PREV’AIR (Honoré et al., 2007). Two types of input data are used to simulate concentrations.

Prior to 2010, a configuration similar to the one used in the EURODELTA-Trends project (Colette et al., 2017) is used. The methodology of Colette et al. (2017) is used to recon-struct the emissions of the main air pollutants (non-methane volatile organic compound (NMVOC), NOₓ, CO, SO₂, NH₃ and primary PM): the annual emissions of each country, bro-ken down by SNAP (Selected Nomenclature for reporting of Air Pollutants) sectors, are estimated using the GAINS (Greenhouse gases and Air pollution Interactions and Syner-gies) model (Amann et al., 2011) for the years 2000, 2005 and 2010. To derive emissions for intervening years, secto-rial results for 5-year periods are linearly interpolated. Mete-orological data are simulated with the Weather Research and Forecast Model (WRF version 3.3.1; Skamarock et al., 2008) from 2000 to 2010.

For the period 2011 to 2015, year-to-year emissions of the main pollutants are taken from the cooperative program for the monitoring and evaluation of the long-range transmission of air pollutants in Europe (EMEP) avail-able at https://www.ceip.at/webdab-emission-database/ emissions-as-used-in-emep-models (last access: 19 May 2022). Annual meteorological data were pro-vided by ECMWF with the Integrated Forecasting System (IFS) model with data assimilation.

For these two datasets, the spatialization of emissions over France is performed with a 1 km proxy based on the national bottom-up emission inventory (available at http: //emissions-air.developpement-durable.gouv.fr/, last access: 19 May 2022) which feeds the CHIMERE emission pre-processor described in Mailler et al. (2017). Furthermore, Denier van der Gon et al. (2015) showed that primary PM emissions from residential wood burning can be underesti-mated by up to a factor of 2–3 over Europe because the emis-sions largely lack semi-volatile compounds. To compensate for this underestimation, a country correction factor deter-mined from Denier van der Gon et al. (2015) is applied over the whole period.

## 2.3 Kriging

Hourly atmospheric concentration fields are estimated by universal kriging, a geostatistical method. Kriging aims to es-timate the value of a random variable (random process which describes the observations) at locations from the measure-ments. Kriging relies on the concept of spatial continuity, which implies that measurements that are close to each other will be more similar than distant measurements. In addition, kriging requires a good knowledge of the spatial structure of the interpolation domain which is represented by the var-iogram or co-variogram (second-order properties) of a ran-dom function (Goovaerts, 1997; Wackernagel, 1996; Chiles and Delfiner, 2012; Lichtenstern, 2013). Kriging involves de-riving the linear combination of the observations which en-sures the minimal estimation variance under a non-bias con-dition. At a point s₀, the concentration estimate ̂y(s₀) is given by Eq. (1).

$$ \hat{y}(s_0) = \sum_{i=1}^{N} \lambda_i y(s_i), \qquad (1) $$
<figure><img src="image_3.png" /><figcaption>Figure 1. PM₁₀: statistical indicators calculated using cross-validation technique on daily mean PM₁₀ values measured and estimated over rural background stations for the years 2007 to 2015. (a) Number of rural stations for each year. (b) Mean bias (black circles), RMSE (colored rectangles), correlation (grey crosses and the associated dashed lines) and mean observation (horizontal lines).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 2. PM₁₀: statistical indicators calculated using cross-validation technique on daily mean PM₁₀ values measured and estimated over URBAN background stations for the years 2007 to 2015. (a) Number of rural stations for each year. (b) Bias (black circles), RMSE (colored rectangles), correlation (grey crosses and the associated dashed lines) and mean observation (horizontal lines).</figcaption></figure>

where $y(s_i)$, $i = 1...N$, are the observed concentrations at sampling locations through the entire domain (unique neighborhood) or within a limited neighborhood of $s_0$ (moving neighborhood), and $\lambda_i$, $i = 1...N$, are the kriging weights.

Among the kriging methods, the universal kriging (especially external drift kriging) allows us to consider additional information to make estimates more accurate. This approach is based on a linear regression with auxiliary variables and a spatial correlation of the residuals, and it allows us to combine simultaneously observations and additional information. The main hypothesis is that the global mean of the random variable is not constant through the domain, and it relies on

explanatory variables. This kriging technique has been used for several years in monitoring air quality systems for spatial interpolation at the regional scale (PREV’AIR; Malherbe and Ung, 2009). For $y(s_0)$, which is the pollutant concentration to be estimated at a location $s_0$, the hypothesis is a linear relation between $y(s_0)$ and the considered auxiliary variables as explained by Eqs. (2) and (3).

$$y(s_0) = m(s_0) + \varepsilon(s_0), \qquad (2)$$

$$m(s_0) = b_0 + b_1x_1(s_0) + b_2x_2(s_0) + \dots + b_px_p(s_0), \qquad (3)$$

where $m(s_0)$ is the drift of the mean, $b_0, b_1, \dots, b_p$ are the coefficients of the linear regression, and $x_0, x_1, \dots, x_p$ are the
<figure><img src="image_5.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 3.** PM<sub>2.5</sub>: statistical indicators calculated using cross-validation technique on daily mean PM<sub>2.5</sub> values measured and estimated over rural background stations for the years 2009 to 2015. (a) Number of rural stations for each year. (b) Bias (black circles), RMSE (colored rectangles), correlation (grey crosses and the associated dashed lines) and mean observation (horizontal lines).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 4.** PM<sub>2.5</sub>: statistical indicators calculated using cross-validation technique on daily mean PM<sub>2.5</sub> values measured and estimated over URBAN background stations for the years 2009 to 2015. (a) Number of rural stations for each year. (b) Bias (black circles), RMSE (colored rectangles), correlation (grey crosses and the associated dashed lines) and mean observation (dotted horizontal lines).</figcaption></figure>

auxiliary variables. ε corresponds to the stationary random process which is associated with a semi-variogram. In addition, the kriging weights must satisfy the drift condition described in Eq. 4.

$$ \forall x_p : x_p(s_0) = \sum_{i=1}^{N} \lambda_i x_p(s_i) \quad (4) $$

In this work, kriging is performed with surface monitoring observations, and the drift is described by the outputs from

the CHIMERE chemistry transport model. European stations located outside the French domain are included in the kriging to increase accuracy at the borders. The kriging is performed using a moving neighborhood as this allows for local adjustment of the relationship between the measurements and CHIMERE. The concentration at each grid point is estimated within a window of 80 monitoring sites. This number has been adjusted in previous studies by sensitivity tests (Benmerad et al., 2017; Beauchamp et al., 2017). In addi-
<figure><img src="image_9.png" /><figcaption>**Figure 5.** O₃: statistical indicators calculated using cross-validation technique on daily mean O₃ values measured and estimated over rural background stations for the years 2000 to 2015. (a) Number of rural stations for each year. (b) Bias (black circles), RMSE (colored rectangles), correlation (grey crosses and the associated dashed lines) and mean observation (horizontal lines).</figcaption></figure>

tion, smoothing is applied to avoid discontinuities in the map (Beauchamp et al., 2015); the smoothing methodology was adapted from Rivoirard and Romary (2011). The final output resolution is the same as for the CHIMERE model: approximately 4 km resolution (0.06° × 0.03°).

For PM₁₀ (particles with a radius < 10 µm) and PM₂.₅ (particles with a radius < 2.5 µm) a co-kriging with external drift is applied. Co-kriging is an extension of kriging to the multivariate case. It allows the estimate of PM₁₀ or PM₂.₅ concentrations by a linear combination of the two-variable data. The particularity of co-kriging is the use of the cross variance or semi-variance between the principal variable and the secondary variable. In the case of co-kriging with external drift, the simple and cross variograms are built based on

residuals (De Fouquet et al., 2007). Co-kriging allows us to take into account the correlation between PM₁₀ and PM₂.₅ and to improve consistency between PM₁₀ and PM₂.₅ estimates (Beauchamp et al., 2015). This co-kriging also allows the PM₂.₅ estimate to benefit from the higher density of PM₁₀ monitoring stations.

## 2.4 Output: regulatory air quality indicators

From the hourly kriged concentrations, several air quality indicators (regulatory and used in health impact assessments) are calculated and mapped over France. The complete list and definition of these indicators are given in Table 2.
<figure><img src="image_10.png" /><figcaption>**Figure 6.** O₃: statistical indicators calculated using cross-validation technique on daily mean O₃ values measured and estimated over URBAN background stations for the years 2000 to 2015. (a) Number of urban stations for each year. (b) Bias (black circles), RMSE (colored rectangles), correlation (grey crosses and the associated dashed lines) and mean observation (horizontal lines).</figcaption></figure>

## 3 Data validation

Usually the quality of the estimated concentration maps is assessed using statistical indicators that compare observations and estimated concentrations at the monitoring stations in the domain. Here, information of all background stations in the domain is already used to produce the maps. Therefore, for a fair comparison, the cross-validation method is used. The cross-validation method calculates the quality of the spatial interpolation for each measurement station point from all available information except the selected station point; i.e., it retains one data point and then makes a prediction at the spatial location of this point. This procedure is repeated for all measurement points in the available set, thus allowing the quality of the predicted values to be assessed at lo-

cations without measurements (provided they are within the area covered by the measurements).

It was noticed that the scores are systematically different at rural and urban stations (even though the kriging technique used here is not differentiated by the type of station). This is why the results of the cross-validation are described by pollutant and differentiated by station type (rural and urban types are presented here). Three statistical indicators are calculated on the basis of the daily average concentration: the mean bias, the root mean squared error (RMSE) and the Pearson correlation (r²). For each year, they are first calculated for the “left out” station, and then the median values over all stations are calculated.

Leave-one-out validation is a commonly used method in the air quality community (see, for example, ETC reports on air quality mapping (Horálek et al. 2017)) which is
<table><caption>**Table 2.** Yearly regulatory air quality indicators from EU legislation or French legislation and usual indicators.</caption><thead><tr><th>ID</th><th>Pollutant</th><th>Statistics</th><th>Threshold</th><th>Threshold origin</th><th>Target to protect</th></tr></thead><tbody><tr><td>NO2_avgannual</td><td>NO<sub>2</sub></td><td>Yearly average</td><td>40 µg m<sup>-3</sup></td><td>Limit value (EU)</td><td>Human health</td></tr><tr><td>O3_avgannual</td><td>O<sub>3</sub></td><td>Yearly average</td><td></td><td></td><td></td></tr><tr><td>O3_AOT40</td><td>O<sub>3</sub></td><td>AOT40* from May to July</td><td>6000 µg m<sup>-3</sup></td><td>Long-term objective</td><td>Vegetation</td></tr><tr><td>O3_AOT40_5years</td><td>O<sub>3</sub></td><td>AOT40* from May to July (5-year average)</td><td>18 000 µg m<sup>-3</sup></td><td>Target value (EU)</td><td>Vegetation</td></tr><tr><td>O3_SOMO35</td><td>O<sub>3</sub></td><td>Sum of excess of max daily 8 h averages over 35 ppb (= 70 µg m<sup>-3</sup>) calculated for all days in a year; SOMO35 (sum of means over 35 ppb)</td><td></td><td>Health impact assessment</td><td>Human health</td></tr><tr><td>O3_T120</td><td>O<sub>3</sub></td><td>Number of days for which the running average over 8 h exceeds 120 µg m<sup>-3</sup></td><td></td><td>Quality objective (EU)</td><td>Human health</td></tr><tr><td>O3_T120_3years</td><td>O<sub>3</sub></td><td>Number of days for which the running 8 h average exceeds 120 µg m<sup>-3</sup> (averaged over 3 years)</td><td>Not to exceed more than 25 d a year</td><td>Target value (EU)</td><td>Human health</td></tr><tr><td>O3_T180</td><td>O<sub>3</sub></td><td>Number of hours exceeding the average value of 180 µg m<sup>-3</sup></td><td></td><td>Recommendation and information threshold (France)</td><td>Human health</td></tr><tr><td>O3_T240</td><td>O<sub>3</sub></td><td>Number of hours exceeding the average value of 240 µg m<sup>-3</sup></td><td></td><td>Alert threshold (France)</td><td>Human health</td></tr><tr><td>PM10_avgannual</td><td>PM<sub>10</sub></td><td>Yearly average</td><td>40 µg m<sup>-3</sup></td><td>Limit value (EU)</td><td>Human health</td></tr><tr><td>PM10_t50</td><td>PM<sub>10</sub></td><td>Number of days exceeding the average value of 50 µg m<sup>-3</sup></td><td>Not to exceed more than 35 d a year</td><td>Limit value (EU)</td><td>Human health</td></tr><tr><td>PM10_t80</td><td>PM<sub>10</sub></td><td>Number of days exceeding the average value of 80 µg m<sup>-3</sup></td><td></td><td>Alert threshold (France)</td><td>Human health</td></tr><tr><td>PM25_avgannual</td><td>PM<sub>25</sub></td><td>Yearly average</td><td>25 µg m<sup>-3</sup></td><td>Limit value (EU)</td><td>Human health</td></tr></tbody></table>

presently recommended by FAIRMODE (Janssen and Thunis, 2020). However, scores derived from the results of the leave-one-out validation might be influenced by areas where the density of sampling points is highest. For this reason, during the FAIRMODE project (Riviere et al., 2019), for which a kriging method similar to the one conducted here was conducted, a comparison has been performed between cross-validation results obtained by the leave-one-out cross-validation and cross-validation results obtained by the 5-fold cross-validation (cross-validation leaving 20% of stations out). Results and related scores were very similar. We therefore decided to continue with the leave-one-out cross-validation process for the validation of our kriging results.

## 3.1 PM<sub>10</sub>

The scores show a good representation of the observations by the kriged data with correlations between 0.77 and 0.86 and RMSE of about 7 µg m⁻³, i.e., between 30% and 50% of the annual mean PM<sub>10</sub> concentration. The mean biases are particularly low for urban stations with values below −1%. For rural stations the average bias is less than +3 µg m⁻³, i.e., less than +15%. The proportion between rural and urban

stations varies between 1/3 and 1/10. The larger number of urban stations allows a better capture of the spatial variability in concentrations in urban environments.

Looking at the evolution of the scores over the years for rural stations, the number of stations available first increases from 2009 to 2012 before decreasing until 2014. In 2015 a new increase in the number of stations in France begins. For urban stations, the decrease starts earlier (2010), but the evolution is the same. The temporal evolution of the scores generally follows the number of stations with higher correlations and smaller relative mean biases and RMSE when more stations are available. Indeed, the greater the number of stations there are, the more representative the kriging technique will be of the real spatial variability. There are exceptions, however, like in 2015 for rural stations, which had the second worst scores even though that year has the largest number of stations.

## 3.2 PM<sub>2.5</sub>

There are between a half and a third as many PM<sub>2.5</sub> stations as PM<sub>10</sub> stations. However, by using a co-kriging technique, the PM<sub>2.5</sub> mapping also benefits from PM<sub>10</sub> information so

<footnote>* AOT40 (expressed in µg m⁻³ h⁻¹) means the sum of differences between hourly concentrations greater than 80 µg m⁻³ (= 40 ppb or parts per billion) and 80 µg m⁻³ for a given period using only the 1 h values measured daily between 08:00 and 20:00 CET.</footnote>
<figure><img src="image_11.png" /><figcaption>**Figure 7.** NO₂: statistical indicators calculated using cross-validation technique on daily mean NO₂ values measured and estimated over rural background stations for the years 2000 to 2015. (a) Number of rural stations for each year. (b) Bias (black circles), RMSE (colored rectangles), correlation (grey crosses and the associated dashed lines) and mean observation (horizontal lines).</figcaption></figure>

that the correlations, mean bias and RMSE are almost similar to the PM₁₀ scores. The mean biases for rural stations do not exceed 20% of the mean concentrations and are very low for urban stations (between 0 and -3%). As for PM₁₀, this bias is systematically positive for rural stations (overestimation) and slightly negative over urban stations (underestimation). This is mainly related to the resolution of the data which smoothes the concentration gradients, giving a unique value for each grid (about 4 km horizontal resolution). For urban station, located close to PM₂.₅ precursor emissions and generally having high concentration values, this smoothing effect leads to an underestimation. For rural areas far from emission precursors, the opposite is observed.

The correlation is generally higher than 0.8, and the RMSE does not exceed 7 µg m⁻³ (at most 50% of the annual mean concentration).

### 3.3 O₃

The comparison between estimated and observed ozone at rural stations shows good correlations (0.8 to 0.87), small relative mean negative biases (−4% to −8%) and low RMSE (around 20% of the annual mean concentration). Between 2000 and 2007, the number of rural stations increased, resulting in improved modeled concentration maps. The small decrease in the number of stations after 2007 does not penalize the scores for these years.
<figure><img src="image_12.png" /><figcaption>**Figure 8.** NO₂: statistical indicators calculated using cross-validation technique on daily mean NO₂ values measured and estimated over URBAN background stations for the years 2000 to 2015. (a) Number of urban stations for each year. (b) Bias (black circles), RMSE (colored rectangles), correlation (grey crosses and the associated dashed lines) and mean observation (horizontal lines).</figcaption></figure>

The same conclusions can be drawn for the urban ozone scores. The higher number of urban stations even leads to slightly better scores, with correlations above 0.9 for all years and relative mean positive biases not exceeding 5 %. A satisfactory RMSE is also obtained for all years with values around 20 % of the annual mean concentration. It can be seen that the positive and negative biases are reverse with respect to the PM scores. Indeed, the highest O₃ values are generally observed in rural areas where precursors have had time to produce O₃ and where O₃ destruction is lower than in urban areas. Therefore, the smoothing effect has the opposite effect to that of PM.

### 3.4 NO₂

Rural scores for NO₂ are worse for particles and O₃. The correlations are between 0.55 and 0.7, but above all, strong positive biases are observed for all years with an overestimation of the observations of 60 % to 80 %. This also affects RMSE scores that can exceed 100 % of the annual mean concentration. This poor performance can be explained by the strong spatial gradients in NO₂ concentrations due to its shorter atmospheric lifetime than O₃ or particles. There are too few rural stations to properly capture this variability in the kriging technique used here, so the urban stations have too much weight, and the raw model concentrations also overestimate rural concentrations.

The urban scores for NO₂ are much better than the rural scores. The correlations are around 0.8, the biases do not ex-
<figure><img src="image_13.png" /><figcaption>Figure 9.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure 9.** PM₁₀ annual mean concentrations from 2000 to 2015. Concentrations are obtained through a combination of regional modeling and observations.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Figure 10.** Trends in PM₁₀ annual mean concentration. Sen slope coefficient (µg m⁻³ yr⁻¹) calculated over the period 2000–2015.</figcaption></figure>

ceed -3.5%, and the RMSE is between 10 and 12 µg m⁻³ (less than 25% of the annual mean concentration). The high number of urban background stations seems satisfactory to allow the kriging technique to correctly reproduce the spatial variability in NO₂ in urban background environments. It should be noted, however, that traffic stations are not used in the present analysis (whether as observational data to be compared with or included in kriging).

## 3.5 Comparison with other scores

In order to evaluate the added value of the kriging technique compared to the raw CHIMERE model simulations, the cross-validation scores can be compared to the raw model scores. Table 3 shows the scores averaged over all years and all observations without distinction of typology.

All scores are strongly improved by the kriging method of observations with CHIMERE in external drift. However, as can be seen in the previous figures, this improvement is more pronounced in urban areas than in rural areas due to the much larger number of stations in urban areas, which makes the kriging more representative of these areas.

The cross-validation scores can also be compared with those obtained in Europe with other mapping methods. Chen et al. (2019) compared 16 algorithms to develop Europe-wide spatial models of PM₂.₅ and NO₂, included linear stepwise regression, regularization techniques and machine learning methods. Those models were developed based on the 2010 routine monitoring data from the Airbase dataset, satellite observations, dispersion model estimates and land-use variables as predictors. De Hoogh et al. (2018) also performed

<table><caption>**Table 3.** Validation scores for the raw data and the kriged concentrations (cross-validation). Annual scores (bias, RMSE and the Pearson correlation coefficient r²) are calculated over France for all years and all stations and are averaged.</caption><thead><tr><th></th><th>NO<sub>2</sub></th><th>O<sub>3</sub></th><th>PM<sub>10</sub></th><th>PM<sub>2.5</sub></th></tr></thead><tbody><tr><th colspan="5">Raw</th></tr><tr><th>Bias</th><td>-3.51</td><td>3.46</td><td>-8.91</td><td>-4.02</td></tr><tr><th>RMSE</th><td>12.97</td><td>17.26</td><td>12.63</td><td>8.73</td></tr><tr><th>r<sup>2</sup></th><td>0.55</td><td>0.73</td><td>0.71</td><td>0.75</td></tr><tr><th colspan="5">Kriged concentration</th></tr><tr><th>Bias</th><td>-0.51</td><td>-0.07</td><td>-0.04</td><td>-0.15</td></tr><tr><th>RMSE</th><td>10.41</td><td>12.54</td><td>7.64</td><td>5.83</td></tr><tr><th>r<sup>2</sup></th><td>0.81</td><td>0.92</td><td>0.85</td><td>0.87</td></tr></tbody></table>

cross-validation of their fine-spatial-scale land-use regression models (also based on Airbase dataset, satellite observations, dispersion model estimates and land-use variables as predictors) used in Europe for the year 2010. Results from their cross-validation are compared to our own cross-validation results in Table 4.

The comparison of performance in these three studies is of course limited by the fact that the spatial coverage differs: in De Hoogh et al. (2018) and Chen et al. (2019), the cross-validation is computed over the whole of Europe. In this study, the performances are assessed over France.

For all pollutants the spatial correlation (r²) is better in our study. At the same time, higher RMSEs are also found for our study. This may be due to a larger bias, but we also demonstrated in our paper that the bias was very small, except at rural NO₂ stations. Since the RMSE score also depends on the absolute concentrations, the different spatial coverage may also play a role. The lower RMSE over Europe could be an artifact of including areas where absolute concentrations of NO₂, PM₂.₅ or O₃ are lower than over France.

The validation scores obtained, as well as the comparison with raw data and with other mapping methods, allow us to be confident about the validity of the concentrations obtained and their good representativeness of background concentrations, in particular in urban areas. A crucial step appears, however, when it comes to the representativeness of rural NO₂ concentrations which are overestimated in our results.

## 4 Results and discussion

After ensuring the validation of the kriged concentration data, yearly indicators, trend over years and human exposition are calculated. Hourly concentrations fields are produced from 2000 to 2015 for NO₂, O₃ and PM₁₀; however, as explained in Sect. 2, for PM₁₀ only annual mean indicator maps are produced before 2007. PM₂.₅ hourly concentrations
<figure><img src="image_16.png" /><figcaption>**Figure 11.** PM<sub>2.5</sub> annual mean concentrations from 2009 to 2015. Concentrations are obtained through a combination (kriging) of regional modeling and observations.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Figure 12.** Trends in PM<sub>2.5</sub> annual mean concentration. Sen slope coefficients (µg m<sup>-3</sup> yr<sup>-1</sup>) calculated over the period 2009–2015.</figcaption></figure>

are calculated for the years 2009 to 2015 due to the limited number of background stations available before 2009.

## 4.1 Concentration maps and trends

All the indicators presented in Sect. 2 are calculated, but the following sections focus on the annual averaged concentrations of PM₁₀, PM₂.₅, NO₂ and O₃, as well as SOMO35 and AOT (two indicators associated with O₃), for which mapped data are presented. These indicators are presented in this paper and available on a Zenodo repository and on an online map library (see Sect. 5). A total of 34 trend analyses over the period are performed by calculating the Sen–Theil regression slope for each grid point on the domain. To characterize the significance of these trend slopes, the 95 % confidence interval is calculated. This confidence interval represents the lower and upper values above or below which there is 95 % confidence that the trends will occur. The smaller the confidence interval is, the more statistically significant the trend will be. Large confidence intervals are considered as unrepresentative, especially those containing 0. Trend slopes and confidence intervals are calculated for each grid point in the domain, and country-averaged values are also given in Table 5.

### 4.1.1 PM₁₀

Maps of annual average PM₁₀ concentration maps are presented in Fig. 9 for the period 2000–2015. The resolution of the grid (around 4 km) allows us to see patterns such as interconnected cities, especially in the latest years during

which the patterns of large inter-regional concentrations are decreasing. The impact of meteorological conditions is also visible through the interannual variability. For example, the 2003 heat wave year is associated with higher PM₁₀ levels due to the increased formation of secondary aerosols.

Figure 10 shows the mapped trends in annual average PM₁₀ expressed as Sen–Theil regression slope in micrograms per cubic meter per year and calculated over the period 2000–2015.

There is a downward trend in PM₁₀ annual mean concentrations everywhere in France and in particular in the regions with the highest PM₁₀ concentrations at the beginning of the period: the south of France (east and west), the Auvergne–Rhône–Alpes region, the east (Grand Est) and the extreme north of France. A country-averaged downward trend in PM₁₀ concentrations of −0.8 µg m⁻³ per year is estimated over the period 2000–2015 (spatial average of the trends calculated on each grid point). This trend is statistically significant on average over France with a narrow 95 % confidence interval – [−0.50, −1.09] – that does not include zero (see Table 5) and applies to almost all grid points (maps of confidence interval not shown here). Taking the year 2000 as the base year, this amounts to a 39 % reduction. In a study conducted for France over the period 2000–2010, Malherbe et al. (2017) estimated a downward trend that was twice as small (0.4). This reflects the accelerated decline in concentrations in France in recent years.

This significant downward trend is the result of the decrease in primary pollutant emissions over these 16 years in response to emission reduction measures. From 2000 to 2015, primary PM₁₀ emissions over France were reduced by 39 %, as well as emissions of PM₁₀ precursors such as NOₓ emissions (−56 %) and SOₓ emissions (−87 %) (data calculated by CITEPA (Interprofessional Technical Centre for Studies on Air Pollution) and extracted from the 2015 French national air quality report: https://www.statistiques.developpement-durable.gouv.fr/sites/default/files/2018-10/datalab-bilan-de-la-qualite-de-l-air-en-france-en-2015-octobre-2016-c.pdf, last access: 19 May 2022).

### 4.1.2 PM₂.₅

The highest PM₂.₅ values are observed at the beginning of the period and are more concentrated in the main source regions than PM₁₀. Significant reductions in annual average background concentrations are observed over the years. The Sen slope coefficients calculated for the annual average PM₂.₅ (Fig. 12) over the period show negative trends over the whole territory and more pronounced ones over the southeast region, the Auvergne–Rhône–Alpes region, northern France and Brittany. A downward trend of −0.87 µg m⁻³ per year from a national average is calculated, again with statistical significance (95 % interval of [−0.48, −1.41] which does not contain zero). Taking 2009 as a reference year, this
<figure><img src="image_18.png" /><figcaption>Figure 13.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**Figure 13.** SOMO35 indicator for the period 2000 to 2015. Ozone concentrations are obtained through a combination (kriging) of regional modeling and observations.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>**Figure 14.** Trends in annual mean O₃ concentrations (in µg m⁻³ yr⁻¹) (**a**), as well as SOMO35 (in µg m⁻³ d⁻¹ yr⁻¹) (**b**) and AOT40 (in µg m⁻³ h⁻¹ yr⁻¹) (**c**) indicators. Sen slope are calculated over the period 2000–2015.</figcaption></figure>

amounts to a 35% decrease in 7 years. As for PM₁₀, this negative trend is associated with the reduction in primary PM₂.₅ emissions and in PM₂.₅ precursor emissions (SOₓ, NOₓ and volatile organic compounds (VOCs)).

### 4.1.3 Ozone

The SOMO35 indicator shows strong interannual variability. O₃ is a photochemical pollutant produced by secondary reactions in the presence of NOₓ, VOCs and sunlight. The hot year 2003 is distinguished by a very high SOMO35 over almost the entire territory. For each year, the highest SOMO35 is found in southeastern France and to a lesser extent in the Alsace region. The trends in SOMO35, annual average O₃

and AOT40 over the years are shown in Fig. 14 for the period 2000–2015.

For the O₃ average annual concentration, small positive trends are found over France. Two exceptions are the southeast (Provence–Alpes–Côte d’Azur, or PACA, region) and the Grand Est region (east of France), i.e., the regions with the highest O₃ concentrations, showing negative trends. Averaging over France, this leads to a positive trend of 0.32 µg m⁻³ yr⁻¹ which corresponds to an increase of 6.5% over 16 years. The same order of magnitude was found for the period 2000–2010 by Malherbe et al. (2017). Both negative (in the south of France) and positive trends are significant according to the mapped 95% confidence interval (not shown). SOMO35 and AOT40 indicators, which are indica-
<figure><img src="image_21.png" /><figcaption>Figure 15.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>**Figure 15.** NO₂ annual mean concentrations for the period 2000 to 2015. NO₂ concentrations are obtained through a combination of regional modeling and observations.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>**Figure 16.** Trends in yearly mean NO₂ concentrations. Sen slope coefficients (µg m⁻³ yr⁻¹) are calculated over the period 2000–2015.</figcaption></figure>

<table><caption>**Table 4.** Validation scores for De Hoogh et al. (2018), Chen et al. (2019), and this study. The following scores are calculated by cross-validation for the three studies: Pearson correlation coefficient r², the bias, and the root mean square error (RMSE).</caption><thead><tr><th></th><th></th><th>De Hoogh et al. (2018)</th><th>Chen et al. (2019)</th><th>This study</th></tr></thead><tbody><tr><td rowspan="3">NO₂</td><td>r²</td><td>0.57</td><td>0.57–0.62</td><td>0.81</td></tr><tr><td>RMSE</td><td>9.51</td><td>9–9.5</td><td>10.41</td></tr><tr><td>Bias</td><td></td><td></td><td>-0.51</td></tr><tr><td rowspan="3">PM₂.₅</td><td>r²</td><td>0.58–0.68</td><td>0.48–0.63</td><td>0.87</td></tr><tr><td>RMSE</td><td>2.97–3.3</td><td>3.1–3.9</td><td>5.83</td></tr><tr><td>Bias</td><td></td><td></td><td>-0.15</td></tr><tr><td rowspan="3">O₃</td><td>r²</td><td>0.63</td><td></td><td>0.92</td></tr><tr><td>RMSE</td><td>6.87</td><td></td><td>12.54</td></tr><tr><td>Bias</td><td></td><td></td><td>-0.07</td></tr></tbody></table>

tors with a threshold value below which concentrations are not taken into account, show mostly negative trends. However, according to the value of the mapped 95% confidence interval (not shown here) on most grid points, the confidence interval is wide and contains zero, indicating a lack of significance of the calculated trends. These results are consistent with other European studies (EMEP, 2016; Malherbe et al., 2017) that show an increase in background concentrations and a decrease in O₃ peaks.

<table><caption>**Table 5.** Country-averaged slope and its 95% confidence interval.</caption><thead><tr><th>Indicator</th><th>Mean tendency slope<br>(or mean trend; in<br>µg m<sup>-3</sup> yr<sup>-1</sup>)</th><th>Mean 95%<br>confidence interval<br>(in µg m<sup>-3</sup> yr<sup>-1</sup>)</th></tr></thead><tbody><tr><td>PM<sub>10</sub> – avg annual</td><td>-0.8</td><td>[-0.5, -1.09]</td></tr><tr><td>PM<sub>2.5</sub> – avg annual</td><td>-0.87</td><td>[-0.48, -1.41]</td></tr><tr><td>O<sub>3</sub> – avg annual</td><td>0.32</td><td>[0.005, 0.59]</td></tr><tr><td>O<sub>3</sub> – SOMO35</td><td>-5.52</td><td>[-102.7, 76.7]</td></tr><tr><td>O<sub>3</sub> – AOT</td><td>-142</td><td>[-641, 315]</td></tr><tr><td>NO<sub>2</sub> – avg annual</td><td>-0.32</td><td>[-0.3, -0.63]</td></tr></tbody></table>

### 4.1.4 NO₂

NO₂ is mainly emitted by road transport. All maps show the same pattern, with cities and interconnected major roads showing the highest NO₂ concentrations. Trends over the period 2000–2015 are shown in Fig. 15. Decreases in NO₂ concentrations are observed in both rural and urban areas throughout the country. However, we recall that rural levels were found to be overestimated with our approach (see Sect. 3.4). The decrease is more important when NO₂ concentrations are high. As with PM₂.₅, these results highlight the combined benefit of large-scale emission management policies that target emission sectors and locally oriented policies.

On average, a significant negative trend of -0.46 µg m⁻³ is calculated over France with a narrow 95% confidence interval (see Table 3). This downward trend is slightly stronger than that calculated in Malherbe et al. (2017) over the period 2000–2010 over France (-0.37 µg m⁻³ yr⁻¹) and corresponds to a reduction of about 30% (taking 2020 as the base year).

## 4.2 Exposure trends

Population-weighted annual average concentrations are good estimates of population exposure as they give more weight to the air pollution where people mainly live. Here, the country-averaged population-weighted concentrations of NO₂, PM₂.₅ and SOMO35, which are the three main indicators used to calculate health impacts, are calculated for each year from the hourly kriged mapped data over France. For one pollutant, it is obtained by adding the result of multiplying the concentration by the population on all the country's grids and then dividing by the total population of the country. The population database used in this study is the LCSQA (L'expertise au service de la qualité de l'air) national population database (Létinois et al., 2014) established for the year 2015. It is based on detailed files from the French Ministry of Finance with information at building level. It is important to note that the French population used here has not varied over the years. The French population increased by about 10% between 2000 and 2015. However, if we considered that the demographic evolution is homogeneous over the country
<figure><img src="image_24.png" /><figcaption>**Figure 17.** Yearly evolution of the country-averaged population-weighted (a) NO₂ concentration, (b) PM₂.₅ concentration and (c) SOMO35. Trends and 95 % confidence intervals are calculated (d).</figcaption></figure>

(the urban/rural ratio has only increased by about 2.5 % in France over the same period), the population-weighted concentration on national average should be the same whatever the year of the population database.

As for the concentrations, a very clear downward trend is observed for population-weighted NO₂ with a trend of -0.5 µg m⁻³ yr⁻¹ and a narrow 95 % confidence interval: [-0.4, -0.6], i.e., a reduction of about 30 % in 16 years. A downward trend of -1 µg m⁻³ yr⁻¹ is also clearly calculated for PM₂.₅ (95 % confidence interval: [-0.6, -1.5]) over the period 2009–2015, i.e., a reduction of about 31 % in 7 years. In contrast, there is no clear trend for the SOMO35 indicator over the period 2000–2015.

When the abovementioned indicators are multiplied by the total population (to obtain the total exposure, i.e., the sum of the population-weighted concentrations over a country), the outcome indicators are those used to calculate the health impact assessment based on dose-response functions, as suggested by the WHO review “Health Risks of Air Pollution in Europe”, described in Holland (2014a, b). Exposure to

SOMO35, anthropic PM₂.₅ and NO₂ (with or without threshold depending on the health impact indicator) contributes to both morbidity and mortality impacts. For example in France, they were used in the PREPA (National Air Pollutant Emissions Reduction Plan) evaluation study for which about 50 political measures to be implemented in France were evaluated and ranked on different criteria, such as air quality impact, health impact and cost-benefit assessment (Schucht et al., 2018). At constant population evolution, the trends are similar between both indicators (total exposure and population-weighted average concentration). However, the evolution in population (even if it is homogeneous over the territory) has an impact on the total exposure of the population. Therefore, we expected a reduced impact on health compared to those on population-weighted concentrations.

## 5 Data availability

Mapped regulatory indicators and exposure data for all 15 years and the four pollutants described here are
available on a Zenodo repository in the Netcdf format (version no. 4) and csv format for data at the municipal or regional level. The DOI link for the dataset is https://doi.org/10.5281/zenodo.5043645 (Real et al., 2021). It is also available through a web-based map library (https://www.ineris.fr/fr/recherche-appui/risques-chroniques/mesure-prevision-qualite-air/20-ans-evolution-qualite-air, Real et al., 2021). The web-based map library is intended to be updated annually. Those data have been provided to several research teams with different fields of expertise ranging from epidemiology to environmental economics and atmospheric science. Most of this work is still in progress, but other papers on the subject have been submitted or are being submitted (Favet et al., 2020; Mink, 2022; Cantrell and Michoud, 2022).

# 6 Conclusions

A 16-year datasets of mapped air pollution concentrations and indicators over France was constructed using a data fusion technique (kriging) that combines measurements from background surface monitoring stations and modeling from the regional model CHIMERE. The resulting data are hourly concentrations at a resolution of about 4 km over France for the period 2000–2015 (shorter period for PM<sub>2.5</sub> and PM<sub>10</sub> hourly indicators).

The kriging technique implemented combines kriging with external drift for NO<sub>2</sub> and O<sub>3</sub> and co-kriging with external drift for particulate matter, allowing the PM<sub>2.5</sub> estimation to benefit from the highest density of PM<sub>10</sub> monitoring stations. These datasets have been evaluated over several years using a cross-validation process that takes into account the incorporation of measurements in the correction process by retaining a data point before calculating the score. The kriging technique significantly improves the validation scores, especially in urban areas with very low biases and high correlations. However, a crucial point appears concerning the representativeness of NO<sub>2</sub> concentrations in rural areas which are overestimated by the model. A new methodology is being developed to better map NO<sub>2</sub> concentrations in these rural areas. It should be noted that the performance increases with the number of measurements taken into account until a threshold is reached at which the addition of stations no longer seems to improve performance. This threshold depends on the pollutant and is higher for pollutants with a strong spatial gradient (i.e., NO<sub>2</sub> which has a shorter lifetime).

The main annual indicators (mean NO<sub>2</sub>, PM<sub>10</sub>, PM<sub>2.5</sub>, O<sub>3</sub>, SOMO35 and AOT40) are analyzed in this article, and the annual trends are calculated. Significant downward trends are calculated over the whole period for annual average concentrations of PM<sub>10</sub>, PM<sub>2.5</sub> and NO<sub>2</sub>. They reflect the reductions in precursor emissions that have taken place in Europe since the 1990s. The trends for O<sub>3</sub> over the 16 years are less

significant. In general, the background O<sub>3</sub> level is increasing, mainly due to large-scale pollution, and high (peaks) O<sub>3</sub> levels are decreasing due to reductions in local O<sub>3</sub> precursor emissions. This results in a positive trend for the annual average O<sub>3</sub> concentration over most of France, but a small downward trend is also observed in the regions with the highest O<sub>3</sub> levels (southeast and east). No significant trend is calculated for the two O<sub>3</sub> indicators detailed here (SOMO35 and AOT40). Population exposure is also calculated over France. The average weight of NO<sub>2</sub> and PM<sub>2.5</sub> in the population of the country decreases respectively by 30% in 16 years and 31% in 7 years. No clear trend was found for the population weight of SOMO35.

**Author contributions.** Data kriging, results evaluation by cross-validation process, and maps and graph production for the paper were performed by ER. The CHIMERE modeling concentration data over the period 2000–2015 were produced by FC. Software developments for the kriging and cross-validation methods were provided by AU, LM and AG. The web-based map library used to store and visualize the data was developed by BR. All work has been supervised and conceptualized by AC. The manuscript draft was mainly written by ER with contributions from all co-authors.

**Competing interests.** The contact author has declared that neither they nor their co-authors have any competing interests.

**Disclaimer.** Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Acknowledgements.** This work was supported by the French ministry in charge of ecology. Part of the simulations was carried out in the XENAIR project funded by the ARC Foundation for Cancer Research.

**Financial support.** This research has been supported by the Ministère de l'Écologie, du Développement Durable et de l'Énergie (grant no. program DRC16).

**Review statement.** This paper was edited by Bo Zheng and reviewed by two anonymous referees.

# References

Amann, M., Bertok, I., Borken-Kleefeld, J., Cofala, J., Heyes, C., Höglund-Isaksson, L., Klimont, Z., Nguyen, B., Posch, M., Rafaj, P., Sandler, R., Schöpp, W., Wagner, F., and Winiwarter, W.: Cost-effective control of air quality and greenhouse gases
in Europe: Modeling and policy applications, Environ. Modell. Softw., 26, 1489–1501, 2011.

Beelen, R., Hoek, G., Pebesma, E., Vienneau, D., de Hoogh, K., and Briggs, D. J.: Mapping of background air pollution at a fine spatial scale across the European Union, Sci. Total Environ., 407, 1852–1867, 2009.

Beauchamp, M.: LCSQA notes, https://www.lcsqa.org/system/files/media/documents/lcsqa2015-note_cartes_analysees-drc-16-152350-11840a_0.pdf (last access: 19 May 2022), 2015.

Beauchamp, M., de Fouquet, C., and Malherbe, L.: Deal-
ing with non-stationarity through explanatory variables in
kriging-based air quality maps, Spatial statistics, 22, 18–46,
https://doi.org/10.1016/j.spasta.2017.08.003, 2017.

Benmerad, M., Slama, R., Botturi, K., Claustre, J., Roux, A., Sage, E., Reynaud-Gaubert, M., Gomez, M., Kessler, R., Brugière, O., Mornex, J.-F., Mussot, S., Dahan, M., Boussaud, V., Danner-Boucher, I., Dromer, C., Knoop, C., Auffray, A., Lepeule, J., Malherbe, L., Meleux, F., Nicod, L., Magnan, A., Pison, C., and Siroux, V.: Chronic effects of air pollution on lung function after lung transplantation in the Systems prediction of Chronic Lung Allograft Dysfunction (SysCLAD) study, Euro. Respir. J., 49, 1600206, https://doi.org/10.1183/13993003.00206-2016, 2017.

Bentayeb, M., Stempfelet, M., Wagner, V., Zins, M., Bonenfant, S., Songeur, C., Sanchez, O., Rosso, A., Brulfert, G., Rios, I., Chaxel, E., Virga, J., Armengaud, A., Rossello, P., Rivière, E., Bernard, M., Vasbien, F., and Deprost, R.: Retrospective model-
ing outdoor air pollution at a fine spatial scale in France, 1989–
2008, Atmos. Environ., 92, 267–279, 2014.

Bessagnet, B., Malherbe, L., and Aymoz, G. : Bilan de la première année de mesure des PM10 ajustées en France et évaluation des outils de modélisation, Rapport LCSQA, https://www.lcsqa.org/system/files/media/documents/DRC_08_94306_15151A_Bilan_PM10_corrige_vf.pdf (last access : 22 May 2022), 2008.

Briggs, D. J., de Hoogh, C., Gulliver, J., Wills, J., Elliott, P., King-
ham, S., and Smallbone, K.: A regression-based method for map-
ping traffic-related air pollution: application and testing in four
contrasting urban environments, Sci. Total Environ., 253, 151–
167, 2000.

Cantrell, C., Michoud, V., Formenti, P., Doussin, J.-F., Alhajj Moussa, S., Cirtog, M., Gratien, A., and Picquet-Varrault, B. and the ACROSS Team: A Future Multi-Platform Atmospheric Chemistry Measurement Campaign to Study Oxidation in Mixed Anthropogenic-Biogenic Air Masses, EGU General Assembly 2021, online, 19–30 Apr 2021, EGU21-15448, https://doi.org/10.5194/egusphere-egu21-15448, 2021.

Chen, J., de Hoogh, K., Gulliver, J., Hoffmann, B., Hertel, O., Ketzel, M., Bauwelinck, M., Donkelaar, A., Hvidtfeldt, H. A., Katsouyanni, K., Janssen, N. A. H., Martin, R. V., Samoli, E., Schwartz, P. E., Stafoggia, M., Bellander, T., Strak, M., Wolf, K., Vienneau, D., Vermeulen, R., Brunekreef, B., and Hoek, G.: A comparison of linear regression, regularization, and machine learning algorithms to develop Europe-wide spatial models of fine particles and nitrogen dioxide, Environ. Int., 130, 104934, https://doi.org/10.1016/j.envint.2019.104934, 2019.

Chiles, J.-P. and Delfiner, P. P.: Geostatistics: Modeling Spatial Un-
certainty, Vol. 497, John Wiley & Sons, 726 pp., 2012.

Colette, A., Andersson, C., Manders, A., Mar, K., Mircea, M., Pay, M.-T., Raffort, V., Tsyro, S., Cuvelier, C., Adani, M., Bessagnet, B., Bergström, R., Briganti, G., Butler, T., Cappelletti, A., Couvidat, F., D'Isidoro, M., Doumbia, T., Fagerli, H., Granier, C., Heyes, C., Klimont, Z., Ojha, N., Otero, N., Schaap, M., Sindelarova, K., Stegehuis, A. I., Roustan, Y., Vautard, R., van Meijgaard, E., Vivanco, M. G., and Wind, P.: EURODELTA-Trends, a multi-model experiment of air quality hindcast in Europe over 1990–2010, Geosci. Model Dev., 10, 3255–3276, https://doi.org/10.5194/gmd-10-3255-2017, 2017.

Couvidat, F., Bessagnet, B., Garcia-Vivanco, M., Real, E., Menut, L., and Colette, A.: Development of an inorganic and organic aerosol model (CHIMERE 2017β v1.0): seasonal and spatial evaluation over Europe, Geosci. Model Dev., 11, 165–194, https://doi.org/10.5194/gmd-11-165-2018, 2018.

De Fouquet, C., Gallois, D., and Perron, G.: Geostatistical charac-
terization of the nitrogen dioxide concentration in an urban area : Part I : Spatial variability and cartography of the annual concen-
tration, Atmos. Environ., 41, 6701–6714, 2007.

De Hoogh, K., Chen, J., Gulliver, J., Hoffmann, B., Hertel, O., Ketzel, M., Bauwelinck, M., Donkelaar, A., Hvidtfeldt, U. A., Katsouyanni, K., Klompmaker, J., Martin, R. V., Samoli, E., Schwartz, P. E., Stafoggia, M., Bellander, T., Strak, M., Wolf, K., Vienneau, D., Brunekreef, B., and Hoek, G.: Spatial PM2.5, NO2, O3 and BC models for Western Europe–Evaluation of spatiotemporal stability, Environ. Int., 120, 81–92, 2018.

Denier van der Gon, H. A. C., Bergström, R., Fountoukis, C., Johansson, C., Pandis, S. N., Simpson, D., and Visschedijk, A. J. H.: Particulate emissions from residential wood combustion in Europe – revised estimates and an evaluation, Atmos. Chem. Phys., 15, 6503–6519, https://doi.org/10.5194/acp-15-6503-2015, 2015.

European Environment Agency (EEA): Air Quality in Europe – 2018 Report, Copenhagen, Report No. 12/2018, 2018.

EMEP: Air pollution trends in the EMEP region between 1990 and 2012, joint report of TFMM/CCC/MSC-E and MSC-W, European Monitoring and Evaluation Programme, Norway, CCC-Report 1/2016, http://www-nilu.no/projects/ccc/reports/cccr1-2016.pdf (last access: 19 July 2019), 2016.

Fayet, Y., Praud, D., Fervers, B., Ray-Coquard, I., Blay, J.-Y., Ducimetiere, F., Fagherazzi, G., and Faure, E.: Beyond the map: evidencing the spatial dimension of health inequalities, Int. J. Health. Geogr., 19, 1–11, https://doi.org/10.1186/s12942-020-00242-0, 2020.

Goovaerts, P. (Eds.): Geostatistics for Natural Resources Evalua-
tion, Oxford University Press, New York, 483 pp., 1997.

Holland, M.: Implementation of the HRAPIE Recommendations for European Air Pollution CBA work, Health Impact Assessment and Cost Benefit Analysis, EMRC, https://ec.europa.eu/environment/air/pdf/CBA%20HRAPIE%20implement.pdf (last access: 19 May 2022), 2014a.

Holland, M.: Cost-benefit Analysis of Final Policy Scenarios for the EU Clean Air Package, Version 2, Corresponding to IIASA TSAP Report 11, Version 1, EMRC, https://ec.europa.eu/environment/air/pdf/TSAP%20CBA.pdf (last access: 19 May 2022), 2014b.

Honoré, C., Menut, L., Bessagnet, B., Meleux, F., Rouïl, L., Vau-
tard, R., Poisson, N., and Peuch, V. H.: PREV’AIR: A platform
for air quality monitoring and forecasting, Chapter 3.4, Dev. En-
vironn. Sci., 6, 293–300, 2007.

Horálek, J., de Smet, P., Corbet, L., Kurfürst, P., and de Leeuw, F.: European air quality maps of PM and ozone for 2010 and their uncertainty, ETC/ACM Technical Paper 2012/12, 74 pp., 2012.

Horálek, J., de Smet, P., de Leeuw, F., Kurfürst, P. and Benešová,
N.: European air quality maps for 2015, ETC/ACM Technical
Paper, vol. 7., 2017.

Horálek, J., Schreiberova, M., Vlasakova, L., Markova, J., Tognet, F., Schneider, P., Kurfürst, P., and Schovankova, J.: European air quality maps for 2018, Eionet Report, ETC/ATNI 2020/10, 2020.

Institute for Health Metrics and Evaluation (IHME):
The Global Burden of Disease: Generating Evidence, Guiding Policy – European Union and European Free Trade Association Regional Edition, Seattle, WA, IHME, https://www.healthdata.org/policy-report/global-burden-disease-generating-evidence-guiding-policy (last access: 20 May 2022), 2013.

Janssen, S. and Thunis, P.: FAIRMODE Guidance Document on Modelling Quality Objectives and Benchmarking, version 3.2, Publications Office of the European Union, Luxembourg, EUR 30264 EN, JRC120649, 54 pp., https://doi.org/10.2760/30226, 2020.

Létinois, L.: Méthodologie de répartition spatiale de la pop-
ulation, LCSQA report, https://www.lcsqa.org/system/files/
media/documents/drc-15-144366-00427a_modelisation_
methodologieglobal-burden-disease-generating-evidence-guiding-policy
(last access: 20 May 2022), 2014.

Lichternstern, A.: Kriging methods in spatial statistics, BA thesis,
Technische Universität München, Department of Mathematics,
https://mediatum.ub.tum.de/doc/1173364/file.pdf (last access: 19
May 2022), 2013.

Mailler, S., Menut, L., Khvorostyanov, D., Valari, M., Couvidat,
F., Siour, G., Turquety, S., Briant, R., Tuccella, P., Bessag-
net, B., Colette, A., Létinois, L., Markakis, K., and Meleux,
F.: CHIMERE-2017: from urban to hemispheric chemistry-
transport modeling, Geosci. Model Dev., 10, 2397–2423,
https://doi.org/10.5194/gmd-10-2397-2017, 2017.

Malherbe, L. and Ung, A.: Travaux relatifs à la plate-
forme nationale de modélisation PREV'AIR: Réalisa-
tion de cartes analysées d'ozone (2/2), rapport LC-
SQA, https://www.lcsqa.org/system/files/media/documents/
drc-10-103351-01139a_prevair-2_vf.pdf (last access: 19 May
2022), 2009.

Malherbe, L., Beauchamp, M., Bourin, A, and Sauvage, S.: Analyse de tendances nationales en matière de qualité de l'air, Rapport final LCSQA, https://www.lcsqa.org/system/files/media/documents/lcsqa2016-tendances_nationales_qa_vf.pdf (last access: 19 May 2022), 2017.

Marécal, V., Peuch, V.-H., Andersson, C., Andersson, S., Arteta,
J., Beekmann, M., Benedictow, A., Bergström, R., Bessagnet,
B., Cansado, A., Chéroux, F., Colette, A., Coman, A., Curier,
R. L., Denier van der Gon, H. A. C., Drouin, A., Elbern, H.,
Emili, E., Engelen, R. J., Eskes, H. J., Foret, G., Friese, E.,
Gauss, M., Giannaros, C., Guth, J., Joly, M., Jaumouillé, E.,
Josse, B., Kadygrov, N., Kaiser, J. W., Krajsek, K., Kuenen,
J., Kumar, U., Liora, N., Lopez, E., Malherbe, L., Martinez, I.,
Melas, D., Meleux, F., Menut, L., Moinat, P., Morales, T., Par-
mentier, J., Piacentini, A., Plu, M., Poupkou, A., Queguiner, S.,

Robertson, L., Rouil, L., Schaap, M., Segers, A., Sofiev, M., Tarasson, L., Thomas, M., Timmermans, R., Valdebenito, Á., van Velthoven, P., van Versendaal, R., Vira, J., and Ung, A.: A regional air quality forecasting system over Europe: the MACC-II daily ensemble production, Geosci. Model Dev., 8, 2777–2813, https://doi.org/10.5194/gmd-8-2777-2015, 2015.

Mink, J.: Putting a price tag on air pollution: the social healthcare costs of air pollution in France, in prep, 2022.

Real, E., Couvidat, F., Ung, A., Malherbe, L., Raux, B., Gressent, A., and Colette, A.: Historical reconstruction of background air pollution over France for 2000–2015, Zenodo [data set], https://doi.org/10.5281/zenodo.5043645, 2021 (data available at: https://www.ineris.fr/fr/recherche-appui/risques-chroniques/mesure prevision-qualite-air/20-ans-evolution-qualite-air, last access: 19 May 2022).

Riviere, E., Bernard, J., Hulin, A., Virga, J., Dugay, F., Charles, M. A., Cheminat, M., Cortinovis, J., Ducroz, F., Laborie, A., Malherbe, L., Piga, D., Real, E., Robic, P. Y., Zaros, C., Seyve, E., and Lepeule, J.: Air pollution modeling and exposure assessment during pregnancy in the French Longitudinal Study of Children (ELFE), Atmos. Environ., 205, 103–114, 2019.

Rivoirard, J. and Romary, T.: Continuity for kriging with moving neighborhood, Mathematical Geosciences, Springer Verlag, 43, 469–481 pp., 2011.

Rouil, L., Honore, C., Vautard, R., Beekmann, M., Bessagnet, B., Malherbe, L., Beekmann, L., Bessagnet, B., Malherbe, L., Meleux, F., Dufour, A., Elichegaray, C., Flaud, J.-M., Menut, L., Martin, D., Peuch, A., Peuch, V.-H., and Poisson, N.: PREV'AIR: an operational forecasting and mapping system for air quality in Europe, B. Am. Meteorol. Soc., 90, 73–84, 2009.

Schucht, S., Real, E., Couvidat, F., Rouil, L., Brignon, J. M., Alle-
mand, N., Le Clercq, G., and Fayolle, D.: Economic analysis of
health impacts in the National Air Pollution Control Programme,
Environnement, Risques & Santé, 17, 393–400, 2018.

Skamarock, W. C., Klemp, J. B., Dudhia, J., Gill, D. O., Barker, D. M., Wang, W., and Powers, J. G.: A description of the Advanced Research WRF version 3, NCAR Technical note-475+ STR, University Corporation for Atmospheric Research, https://doi.org/10.5065/D68S4MVH, 2008.

Shaddick, G., Salter, J. M., Peuch, V.-H., Ruggeri, G., Thomas, M. L., Mudu, P., Tarasova, O., Baklanov, A., and Gumy, S.: Global Air Quality: An Inter-Disciplinary Approach to Exposure Assessment for Burden of Disease Analyses, Atmosphere, 12, 48, https://doi.org/10.3390/atmos12010048, 2021.

Tarasick, D., Galbally, I. E., Cooper, O. R., Schultz, M. G., Ancel-
let, G., Leblanc, T., Wallington, T. J., Ziemke, J., Liu, X., Stein-
bacher, M., Staehelin, J., Vigouroux, C., Hannigan, J. W., Gar-
cía, O., Foret, G., Zanis, P., Weatherhead, E., Petropavlovskikh,
I., Worden, H., Osman, M., Liu, J., Chang, K.-L., Gaudel,
A., Lin, M., Granados-Muñoz, M., Thompson, A. M., Olt-
mans, S. J., Cuesta, J., Dufour, G., Thouret, V., Hassler,
B., Trickl, T., and Neu, J. L.: Tropospheric Ozone Assess-
ment Report: Tropospheric ozone from 1877 to 2016, ob-
served levels, trends and uncertainties, Elem. Sci. Anth., 7, 39,
https://doi.org/10.1525/elementa.376, 2019.

Wackernagel, H.: Multivariate Geostatistics, Springer Berlin Heidelberg, 387 pp., https://doi.org/10.1002/bimj.4710380409, 1996.