key: cord-0888085-5rwmul08
authors: Fioravanti, Guido; Cameletti, Michela; Martino, Sara; Cattani, Giorgio; Pisoni, Enrico
title: A spatiotemporal analysis of NO [Formula: see text] concentrations during the Italian 2020 COVID‐19 lockdown
date: 2022-03-12
journal: Environmetrics
DOI: 10.1002/env.2723
sha: f8620052c9b45d4257d68bb4497492bf352c943d
doc_id: 888085
cord_uid: 5rwmul08

When a new environmental policy or a specific intervention is taken in order to improve air quality, it is paramount to assess and quantify—in space and time—the effectiveness of the adopted strategy. The lockdown measures taken worldwide in 2020 to reduce the spread of the SARS‐CoV‐2 virus can be envisioned as a policy intervention with an indirect effect on air quality. In this paper we propose a statistical spatiotemporal model as a tool for intervention analysis, able to take into account the effect of weather and other confounding factor, as well as the spatial and temporal correlation existing in the data. In particular, we focus here on the 2019/2020 relative change in nitrogen dioxide (NO [Formula: see text]) concentrations in the north of Italy, for the period of March and April during which the lockdown measure was in force. We found that during March and April 2020 most of the studied area is characterized by negative relative changes (median values around [Formula: see text] 25%), with the exception of the first week of March and the fourth week of April (median values around 5%). As these changes cannot be attributed to a weather effect, it is likely that they are a byproduct of the lockdown measures. There are two aspects of our research that are equally interesting. First, we provide a unique statistical perspective for calculating the relative change in the NO [Formula: see text] by jointly modeling pollutant concentrations time series. Second, as an output we provide a collection of weekly continuous maps, describing the spatial pattern of the NO [Formula: see text] 2019/2020 relative changes.

In the last decades a large number of initiatives have been taken in European countries in order to reduce air pollution and its adverse effects on human health. This is for example the case of low emission zones or congestion charge introduced in urban areas (see e.g., Fassò, 2013; Holman et al., 2015; Maranzano et al., 2020) , of more stringent limits on the content of sulfur in marine fuels (e.g., Grange & Carslaw, 2019) or of new air pollution control regulations (e.g., Font & Fuller, 2016; Han et al., 2021) , among others (see Burns et al., 2020 , for a review). A crucial issue in air quality management is the assessment of the efficacy of a specific intervention or a new environmental policy. In particular, quantifying the effective changes in pollutant concentrations due to the adopted strategy can be difficult because of the complexity of air pollution dynamics, strongly depending on weather conditions. Moreover, when the intervention is time-delimited and characterizes different areas, it may be interesting to evaluate the variability in space and time of its effect, if any. The lockdown measures adopted by many countries in 2020 in order to to prevent the spread of the SARS-CoV-2 virus can be considered as a policy intervention, with a possible indirect effect on air quality. In this regard, the main focus is on assessing how the lockdown affected air pollution levels, in particular after controlling for meteorology, long-term trend and other confounding factors. The literature on this topic is obviously quite recent but already large, as discussed in the two recent systematic reviews by Gkatzelis et al. (2021) and Rana et al. (2021) .

In this paper, we propose a new statistical modeling approach for assessing and quantifying the effectiveness of a policy intervention on air quality. In particular, our statistical model is defined for daily differences of pollutant concentrations and has a spatiotemporal specification which gives us the possibility to estimate in space and time the relative change in air pollution levels. We show here an application of our modeling strategy for nitrogen dioxide (NO 2 ) concentrations in northern Italy during the months of the first Italian lockdown (March and April 2020) as compared to 2019. In early 2020, Italy was the first European country to adopt strict lockdown measures (Remuzzi & Remuzzi, 2020) . Starting from late February 2020, people were banned from traveling and all the nonessential socioeconomic activities were stopped. These restrictive measures were initially adopted in some limited municipalities in northern Italy and then were extended to the entire country with the national lockdown in force since March 8, 2020 (Malpede & Percoco, 2020; Sanfelici, 2020) . The restrictions caused, incidentally, a strong decrease in anthropogenic emissions of the main air pollutants, especially for some sectors such as road transport and aviation.

Several environmental studies revealed that restrictions on mobility during the lockdown had a positive effect on NO 2 levels, even if not homogeneously across the considered spatial domains. For example, the broadly publicized data from the Copernicus Sentinel-5P satellite Ali, Abbas, Qamer, Wong, et al., 2021; Bar et al., 2021; Dutta et al., 2021; Muhammad et al., 2020) recorded for NO 2 a sharp drop, in the range 20%-55%, during January-April 2020 compared to 2019 in many cities in China, India, Pakistan, Western Europe, and United States. With regards to Italy, Bauwens et al. (2020) found the average TROPOMI NO 2 column during the lockdown period in 2020 to be between 38 (± 10%) lower than during the same period in 2019 in Milan. Similarly, Cersosimo et al. (2020) reported a general decrease in NO 2 levels over the Po Valley during the lockdown using both satellite and in situ measurements. For the city of Rome and its surroundings, Bassani et al. (2021) documented a 2019 versus 2020 reduction of 50% in NO 2 concentrations using surface measurements from urban traffic stations. The NO 2 decrease for the city of Rome was also documented in Kumari and Toshniwal (2020) . Finally, a general decrease in the NO 2 levels for three cities in Tuscany region was described by Donzelli et al. (2020) , for Reggio Emilia by Marinello et al. (2021) and for Naples by Sannino et al. (2021) . The main weakness of the aforementioned studies consists in their descriptive approach, based on the direct comparison of pollutant concentrations before and after the lockdown or between the lockdown period and the corresponding period of the previous year(s). In other words, these studies make no attempt to adjust for the effect of meteorological conditions which can be adverse or favorable to pollutant dispersion. In this regard, it is worth to note that the first quarter of 2020 experienced exceptional weather conditions with also stronger positive temperatures anomalies over Europe (Barré et al., 2020; van Heerwaarden et al., 2021) .

A number of studies approached the problem of assessing the lockdown effects on air quality using chemical transport models (CTM) in order to compare model forecasts under the business-as-usual (BAU) emission scenario with the observed ground-level measurements or with the expected concentrations computed under a COVID-19 pandemic scenario with reduced emissions (see e.g., Barré et al., 2020; Menut et al., 2020; Piccoli et al., 2020; Putaud et al., 2020; Zhe et al., 2021) . In any case, the use of deterministic models poses a number of practical and conceptual difficulties. First, collecting the needed input data (e.g., emission inventories and meteorology data) is far from straightforward. Second, deterministic models are complex to be run and are not able to properly assess the uncertainty of the results.

Another substantial research line is represented by model-based studies. In this context historical measurement data from previous years (or from the pre-lockdown period) are used to run machine learning algorithms (see e.g., Barré et al., 2020; Diémoz et al., 2021; Granella et al., 2021; Grange et al., 2021; Keller et al., 2021; Kim et al., 2021; Petetin et al., 2020) or to estimate statistical models, as multiple linear regression models (e.g., Bao & Zhang, 2020; Dacre et al., 2020; Hoermann et al., 2021) , Generalized Additive Models (e.g., EEA, 2020; Ordóñez et al., 2020; Solberg et al., 2021) or Autoregressive Integrated Moving Average models (e.g., Tyagi et al., 2020) . The fitted model is then used to predict concentrations for 2020 (or for the post-lockdown period) under the BAU (or counterfactual) scenario, that is, assuming that the lockdown did not take place. In order to correct for the weather effect and temporal trends, meteorological and time variables are included in the model as linear or nonlinear effects. The differences between the predictions, derived from the estimated model, and the observed concentrations (i.e., the out-of-sample prediction errors) are then used to evaluate the effect of the lockdown restrictions. Other papers adopt a different modeling approach and define the counterfactual by using data from cities not subject to lockdown restrictions; then difference-in-differences (DID) models are used to estimate the relative change in pollutant concentrations in the treatment group (locked-down cities) compared with the control group (nonlocked-down cities) (see e.g., He et al., 2020; Wang et al., 2021) . An extension of the DID method is proposed in Zheng et al. (2021) with the aim of computing the would-be average concentrations without the COVID-19 pandemic by removing the meteorological confounding and accounting for the temporal trend. Another modeling strategy makes use of time series model or dynamic panel data model, where pollutant concentrations, measured in 2020 and possibly in previous years, represent the response variable. The lockdown effect is included as a time-dependent dummy variable in the set of regressors, together with meteorological and time variables (see e.g., Bao & Zhang, 2020; Beloconi et al., 2021; Cameletti, 2020) . In this case the effectiveness of the lockdown is evaluated by means of the lockdown dummy coefficient and its interaction with time or other regressors. Whatever the adopted model-based strategy is, pollutant concentrations time series can be analyzed separately for each monitoring station or jointly. The second solution leads to more efficient parameter estimates and a better predictive capability because of the larger amount of available data. More importantly, when measurements come from several spatial sites it is convenient to account also for spatial correlation, besides the temporal one, for explaining any residual variability. This is a standard and well-established option in air pollution modeling (see e.g., Finazzi et al., 2013; Lee et al., 2016; Sahu et al., 2006) given that nearby monitoring stations are expected to show similar pollutant concentrations values. As far as we know, the study of Beloconi et al. (2021) is the only study which implements a spatiotemporal model for evaluating the effect of COVID-19 lockdown on air quality, while all the other papers fail to consider the spatial correlation. However, in Beloconi et al. (2021) the analysis of the impact of the lockdown is limited at the single point stations and no continuous maps are provided for the entire region of interest.

In this paper, using a spatiotemporal statistical model, we aim to produce spatiotemporal maps showing continuously in space and across time the impact of the lockdown on air pollution. We expect these highly spatially resolved maps to help in assessing if the lockdown effect was homogeneous in the considered area or was more consistent in particular zones. Having this goal in mind, we modeled the 2019/2020 daily differences of NO 2 concentrations (in March and April), rather than jointly modeling the data available for the 2 years. This makes it possible to focus on the change occurred between 2019 and 2020 in the NO 2 levels, while still accounting for the effect of meteorology. In particular, by including a spatial stochastic component, our model is able to take care of the spatial correlation between observations and generate spatially continuous prediction surfaces, also where no ground-monitoring stations are available and in remote or mountainous areas (Diémoz et al., 2021) . The model we propose is applied here to the Italian COVID-19 case study but it is a general modeling approach that can be implemented anytime it is necessary to evaluate the effect of an intervention on an output variable with a spatiotemporal dimension.

The paper is structured as follows: in Section 2 we introduce the study domain and the input data. Then, in Section 3 the proposed spatiotemporal regression model is described; details are also provided with respect to the prediction phase of the analysis. Finally, in Section 4 we describe the main results of our analysis with particular attention to the prediction maps. Section 5 concludes the paper.

In this study we focused on the Italian COVID-19 lockdown period (from March 1 to April 30, 2020), corresponding to 10 weeks of daily observations. The 2019 and 2020 raw NO 2 hourly data (expressed as μg∕m 3 ), together with monitoring stations metadata, were extracted from the national database used to store and process the Italian Air quality measurements (SNPA, 2020). The hourly data were measured at stations operated by the local Regional Environmental Protection Agencies, following the European standards EN 14211:2012 for NO, NO 2 , and NO X (CEN-TC 264, 2012). All ground concentrations were fully validated following internal quality assurance and quality control standard procedures.

Daily averages were computed provided that a daily 75% data coverages was achieved (i.e., at least 18 valid hourly records out of 24). The input dataset for our analysis regards 200 monitoring sites with a low proportion of missing data F I G U R E 1 Spatial distribution of the monitoring sites for NO 2 concentrations (red points) and mesh adopted for the implementation of the spatiotemporal model described in Section 3 (< 25% per station) and distributed across eight regions or autonomous provinces in the north of Italy (Valle d'Aosta 4 stations, Piedmont 17, Veneto 32, Lombardy 55, Autonomous province of Trento 5, Friuli Venezia Giulia 12, Autonomous province of Bozen 5, Emilia Romagna 36) plus Tuscany region in the center of Italy (with 34 stations). The spatial distribution of the selected stations is illustrated in Figure 1 . The monitoring stations cover urban (123 stations), suburban (39) and rural (38) areas. At the time of the analysis, we were not able to access to the NO 2 data from Liguria region (the western region filled in gray in the map in Figure 1 ).

The investigated area is characterized: in the north, by remote and relatively pristine areas with the high peaks of the Alps mountain system; in the southern part, by the mountains and hills of the Apennines range (the "spine" of the Italian peninsula); in the center, by the Po Valley. The latter is a large area which includes Lombardy, Emilia-Romagna, Veneto, and Piedmont, regions which account for 48.2% of Italy's GDP (OECD, 2020). Interestingly, in early March 2020, these were the regions with the highest incidence of COVID-19 cases. 1 The Po Valley exhibits intense human pressure and severe pollution levels, with high population density, urban sprawl, intensive agricultural practice and livestock management (Pezzagno et al., 2020; Romano et al., 2020) . This leads to large amounts of NO X emissions from vehicles, methane and ammonia from agricultural activities and particulate matter (PM) from residential heating. The Alps and the Apennines often limit the air flows between the Po valley and the rest of continental Europe causing air pollution stagnation. Table 1 provides the summary statistics for the 2019 and 2020 NO 2 daily concentrations. A certain variability in the NO 2 values is observed across areas, month and year. Nonetheless, it is apparent that lower NO 2 measurements characterize the 2020 months compared with 2019. For example, the monthly NO 2 mean values are in the range 13-35 μg∕m 3 in March 2019 and 12-26 μg∕m 3 in April 2019, while in March and April 2020 the same values vary between 11-24 μg∕m 3 and 7-17 μg∕m 3 , respectively. This means that for the lockdown months of 2020 the range of the monthly mean concentrations is approximately halved compared with the same months of the previous year. The same situation can be observed also in the median and maximum values.

To match the 2019 and 2020 daily NO 2 concentrations and compute the corresponding daily differences, we excluded the use of the calendar date. The reason for this choice is that we do not want to match business days in 2020 with weekend TA B L E 1 Summary statistics -mean, standard deviation (SD), minimum (Min), median, and maximum (Max) -for NO 2 concentrations (in μg∕m 3 ), for the considered region or autonomous province (AP) and for March and April of years 2019 and 2020 days in 2019, and vice versa, given the well-known NO 2 weekly cycle (characterized by lower values during the weekend). Rather, we aligned the 2019 daily measurements of each monitoring site by week number (according to the ISO-8601 standard) and day of the week (Monday, Tuesday, … , Sunday) taking as a reference those of the 2020 daily measurements. A similar approach is described in Ruan et al. (2020) to analyze the impact of COVID-19 on electricity consumption in the United States. By doing so, we were able to compare the first Monday of the first week in March 2020 with the first Monday of the first week in March 2019, the first Tuesday of the first week in March 2020 with the first Tuesday of the first week in March 2019, and so on. However, even using this approach it can happen that a public national holiday (e.g., Easter Monday or April 25) in 2020 is matched with a working day in 2019, and vice versa. The parallel plot in Figure 2 graphically shows the 2019 and 2020 NO 2 aligned datasets. Here, the daily data are visually aggregated at the weekly level. The plot nicely reveals the decrease in the NO 2 concentrations during March and April 2020, compared with 2019, especially during the last 2 weeks of March. At the same time, the positive increments occurring across the weeks in some of the stations suggest that it would be wrong to think that, if the lockdown affected the NO 2 concentrations levels, such an effect was a homogeneous phenomenon both in time and space.

We conclude this section observing that a small fraction (about 3.5%) of aligned daily measurements exhibits unusual combinations of extremely high values in 1 year and extremely low values in the other. These combinations correspond to 2019/2020 relative changes greater than 100% in absolute value, which abnormally lie outside the general pattern seen in our input dataset. These data have no meaning for the purposes of our analysis and were discarded.

For each monitoring site, we retrieved its meteorological conditions and spatial characteristics for a total number of 13 independent explanatory variables (see Table 2 ), whose selection was driven by expert knowledge (Fioravanti et al., 2021) . To describe part of the spatial variation of the NO 2 daily differences, three geographical variables (constant in time) were considered: the linear distance from the the nearest major road, the altitude and the percentage of agricultural and arable lands. Yeganeh et al. (2018) found that traffic volume and congestion data for all of the individual roads can effectively improve the spatiotemporal modeling of NO 2 concentrations. Unfortunately, this information is not available for the whole investigated domain.

How weather affects NO 2 concentrations is assessed using 10 variables (such as total precipitation, wind speed, relative humidity, surface pressure) retrieved from the ERA5-Land dataset from the Copernicus Climate Data Store (https:// climate.copernicus.eu/). ERA5-Land is a regridded version (∼9-km grid spacing) of the ERA5 climate reanalysis (Hersbach et al., 2020) of the European Centre for Medium-Range Weather Forecasts. The use of ERA5 data is documented, among others, by Chan et al. (2021) for investigating the effect of COVID-19 pandemic on NO 2 concentrations over Germany.

At this point it is important to stress that the meteorological variables enter the model as 2020-2019 daily differences. More precisely, for the meteorological daily data we applied the same alignment procedure used for the NO 2 concentrations and described in Section 2.1.1. Afterwards, we calculated the daily differences which we used as model regressors, TA B L E 2 Description of the spatial (altitude, linear distance to the nearest primary road, agricultural and arable lands) and 

Altitude m 1 km

Linear distance to the nearest primary road m 1 km

Agricultural and arable lands % 1 km as explained in Section 3, where the details of the statistical model are given. It is worth to note that the inclusion in the model of the differenced meteorological variables alleviates potential multicollinearity problems and makes it possible to include in the same model regressors like temperature and surface pressure, which otherwise would be highly correlated.

Consider a couple of NO 2 concentrations y m 2019 (t, s i ) and y m 2020 (t, s i ) temporally aligned according to the procedure described in Section 2.1.1, where s i denotes the location (with i = 1, … , 200) and t the day (t = 1, … , T m ) of month m = 1, 2 (where m = 1 denotes March) of year 2019 and 2020, respectively.

As the objective of this study is to evaluate the change in NO 2 levels between 2019 and 2020, we first of all defined the daily differences of the log-transformed NO 2 concentrations:

.

(1)

The logarithmic transformation is a common choice in air pollution analysis in order to reduce the typical positive skewness observed in concentrations distributions (Ott, 1990) .

As the Italian lockdown fell almost at the end of the winter season, when usually the meteorological conditions favor the dispersion of the pollutants, a downward trend in NO 2 concentrations is expected across March and April. This could have a confounding effect when trying to isolate the impact of COVID-19 lockdown measures on air quality. In order to control for this long-term trend, we decided to model the daily differences separately for the 2 months (March and April).

The measurement equation is given by

where (t, s i ) ∼ N(0, 2 ) is the Gaussian measurement error assumed to be independent in space and time. For the sake of simplicity in Equation (2) we omitted the index m given that the model structure is the same for the considered months. The mean level (t, s i ) is then defined as the sum of fixed and random effects as follows:

where 0 is the intercept and 1 t the linear temporal trend which should account for the short-term variation across days in the month. The term z(s i ) is the p z -dimensional vector of the purely spatial regressors, while x (s i , t) is the p x -dimensional vector of the differences computed using the values of the meteorological regressors, as described in Section 2.2. We assume a linear effect for the spatial and spatiotemporal covariates by means of the parameters' vector z and x , respectively. In order to test whether or not there was a variation in the Sunday effect during the lockdown restrictions, we included also the dummy variable I S t which is equal to 1 when t is Sunday and zero otherwise. The corresponding coefficient represents the additional expected change in the mean difference (t, s i ) during Sunday. The term v(s i ) ∼ N(0, 2 v ) represents a temporally and spatially uncorrelated Gaussian random effect which captures some of the small-scale spatial variability. Finally, u(t, s i ) is the residual space-time correlation for which the following first-order autoregressive dynamics was specified:

for t = 2, … , T and |a| < 1. The innovations (t, s i ) have a zero-mean Gaussian distribution, are uncorrelated in time (i.e., Cov

where h is the Euclidean distance between site i and j and (h) is the Matérn covariance function with variance 2 , range and smoothness parameter = 1. For more details about this separable spatiotemporal covariance structure see for example (Cameletti et al., 2011 (Cameletti et al., , 2013 .

The model described in Section 3 is developed in the Bayesian framework and is fully specified once prior distributions are set. Vague Gaussian priors were used for 0 , 1 and for the elements of z and x . For the remaining parameters Penalized Complexity (PC) priors (Simpson et al., 2017) were used. PC priors are a relatively new approach to specify weakly informative priors (Lemoine, 2019) in realistically complex statistical models with the twofold purpose of penalizing model complexity and avoiding overfitting.

For the standard deviation parameters (here and v ) PC priors are generally defined as Prob( > u ) = , where u > 0 is a quantile of the prior and 0 ≤ ≤ 1 is a probability value. In our analysis we set u = 1 and = 0.1 for both and v . This choice can be motivated considering that the total standard deviation of the observed daily differences of the log-transformed NO 2 concentrations is ∼ 0.5 in each month, so it is very likely that the variance of each component is lower than 1.

The joint PC prior suggested in Fuglstad et al. (2019) was used for and . This can be specified through

where we set u = 150, = 0.8, u = 1, = 0.01. Finally, for the autocorrelation parameter a we used the PC prior proposed in Sørbye and Rue (2017) . This can be specified through Prob(a > u a ) = a , where we set u a = 0.8 and a = 0.3. All these choices reflect both previous findings (see e.g., Cameletti et al., 2013; Fioravanti et al., 2021) and restrictions to the possible values of u a and a . Inference was carried out by using the INLA-SPDE approach (Lindgren & Rue, 2015; Martino & Riebler, 2022) , which has been proved to be computationally faster than the Markov chain Monte Carlo (MCMC) approach commonly used for Bayesian inference.

Once the model has been fitted to the observed data, we used Monte Carlo (MC) simulation to generate a large number of samples (say 1000) from the posterior predictive distribution of Δ(t, s) for any location s in the study area and day t:

is the vector of all the model parameters, whose uncertainty is averaged out given all the observed data . For simplicity the generic sampled value will be denoted byΔ(t, s).

To generate prediction maps, we used a raster grid with a spatial resolution of 1 km. For each pixel of this grid, corresponding to the generic location s, we simulated 1000 valuesΔ(t, s) from Equation (4). While computing the predictions, we set the values of x (s, t) equal to zero for each s and t. This corresponds to assume that the meteorological conditions are equal in 2019 and 2020 for each location and time point.

We derived the posterior distribution of the relative change of NO 2 concentrations between 2019 and 2020 by using the following transformation:Δ

whereΔ(t, s) takes negative (positive) values if lower (higher) NO 2 concentrations are expected in 2020 compared to 2019, and equal to zero in case of no change. Finally, we averaged the 1000 MC samples from the daily distributions ofΔ(t, s) at the week temporal resolution. The result is a collection of 1000 predicted raster maps for each week, from which it was possible to obtain the maps of the posterior mean and of the 2.5% and 97.5% quantiles. These latter two identify the bulk of the posterior distribution of each grid cell and were used to determine the 95% credible interval. When a credible interval does not include the zero, the corresponding relative change statistically differs from zero at the significance level of 0.05.

Finally, it is noteworthy to say that to take into account the weekly cycle of NO 2 concentrations, we generated both mean and quantile maps distinguishing between working days (Monday-Saturday) and weekends (Sunday).

R software was used for the model implementation by means of the R-INLA package (https://www.r-inla.org). For the manipulation of the raster maps, we used the R raster package (https://cran.r-project.org/web/packages/raster/ index.html) and the CDO software (https://code.mpimet.mpg.de/projects/cdo). Input data and part of the code used for this study are available through a dedicated GITHUB repository (https://github.com/progettopulvirus/A_spatiotemporal_ analysis_of_NO2_concentrations).

An application of the model to the entire Italian territory is available as an interactive web application (Chang et al., 2021) at the following link: https://guidoispra.shinyapps.io/pulvirus_no2/. Figure 3 shows the monthly posterior distributions for the fixed effects' coefficients ( 0 , 1 , , z and x ). Generally speaking, we observe that a significant effect, invariant in sign across the months of March and April, characterizes most F I G U R E 3 Posterior distribution of 0 (Intercept), 1 (Day), (Sunday) and of the covariate coefficients x and z . When the bulk of the posterior distribution lies far from zero (dashed line) it can be concluded that the corresponding parameter is significantly different from zero of the selected regressors. The shape of the posterior distributions is rather stable in the models for the 2 months, but some exceptions are apparent. The distribution of the linear trend coefficient 1 is narrower in March than in April, while the opposite happens for the posterior distribution of the intercept 0 . This is in line with the large variability which characterizes the weekly prediction maps of April (see Section 4.3). Finally, the surface pressure coefficient exhibits two almost flat posterior distributions, but this result is not of easy interpretation. A significant negative linear trend coefficient 1 was found in March, which suggests a decreasing trend for the daily differences of the log-transformed NO 2 concentrations in the month when the lockdown restrictions occurred. A significant negative effect was also found for the Sunday coefficient . We can interpret this result saying that a weekly cycle TA B L E 3 Posterior mean and standard deviation (in parentheses) of the parameters in the models for March and April. a: AR(1) coefficient; and : range (in km) and standard deviation of the Matérn spatial covariance function; v : standard deviation of the small-scale spatial random effect; : standard deviation of the Gaussian measurement error persists even when we consider the difference of the log-transformed NO 2 concentrations. Surprisingly, neither the linear trend nor the Sunday dummy variable coefficients are significant in April. With regards to the meteorological parameters, we distinguish those with a positive significant effect (the diurnal temperature range, the average temperature at 2 m, the relative humidity and surface pressure) and those with a negative significant effect (the min and max planet boundary layer, the net irradiance, the wind speed and the total precipitation). Conversely, all the spatial regressors (elevation, % of agricultural/arable land and distance to major roads) show a positive posterior mean. This could suggest that far from the urbanized centers and the road network the level of the NO 2 daily concentrations in 2020 tends to be equal or greater than the ones in 2019. Table 3 provides the posterior summary statistics for the remaining model parameters. We observe that the posterior mean of the AR(1) autocorrelation coefficient a, the spatial range and the standard deviation v have a larger posterior mean in April than in March.

For validation purposes, we stratified the input monitoring sites according to their type classification (urban, suburban, and rural stations). Within each of these groups, we sampled 10% of the stations in order to define a validation dataset. The remainder of the stations (training dataset) was used to fit the model and the fitted model was used to predict the daily differences of the log-transformed NO 2 concentrations (see Equation 1) on the validation dataset. Both the sampling and the estimation process were repeated three times for each month. Comparing predicted and observed values allows to investigate how the model performs on unseen data, specifically if the model generalizes well or suffers from overfitting/underfitting. Plots and summary statistics are shown as percentage relative changesΔ(t, s)% (see Equation 5 ) for ease of interpretation of the validation results. The scatterplots in Figure 4 shows the distribution of the fitted values versus the observed values. The points spread uniformly along the diagonal line, showing good agreement between observed and modeled data both in the training and validation stage. As expected, modeled relative change for the validation stage shows greater variability than in the training stage. Also the Pearson correlation coefficients support the good model performance with a value of 0.9 for the training stage and of ∼ 0.7 for the validation stage. Finally, the Root Mean Squared Error is ∼ 10% and ∼ 20% for the training and validation stage, respectively.

Using the predictors described in Section 2.2 as individual raster files, the procedure described in Section 3.2 was implemented to obtain maps of the relative change of NO 2 concentrations between 2019 and 2020, separately for the weeks of March and April. The maps in Figures 5 and 6 show the spatial pattern of the 2019/2020 relative change of NO 2 concentration in the north/center of Italy for March and April, respectively. These maps refer to the weekly averaged estimates obtained from working days (Monday-Saturday), while the weekly Sunday estimates are available in the Appendix (see Figures A1-A3 ). This distinction between working days and Sunday maps reflects our choice of including a dummy regressor for the Sunday effect in the model and its statistical significance observed in March. For visualization purposes, we limited the data range between −100% and +100% and used a scientifically derived color map (Crameri et al., 2020) . Furthermore, the islands from the Tuscan Archipelago are not displayed and Liguria region, where no daily NO 2 concentrations were available, is represented with a white background. It is worthwhile to stress that the maps presented in this section must be intended as meteorology-normalized maps. This means that they were generated assuming that, in each cell of the raster grid, the daily meteorological conditions are exactly the same in 2019 and 2020. Mathematically, this assumption is equivalent to setting to zero all the meteorological terms in our spatiotemporal model (see Section 3.2).

Figures 5 and 6 reveal a substantial decrease in NO 2 concentrations during March and April 2020 as compared to 2019. A statistically significant reduction persists during the third and fourth working weeks of March across the whole study area: we quantified the interquartile range of the corresponding relative changes distribution to be between −40% and −20%, as shown in the boxplot in the left panel of Figure 7 . Notably, the third and fourth week of March correspond to the stringent phase of the COVID-19 lockdown. The barplots of Figure 8 indicate that Lombardy, In April the predicted surfaces show more spatial variability than in March and all the weekly distributions of the relative changes (right panel of Figure 7 ) exhibit positive increments up to 100%. These results suggest that the NO 2 concentrations began to recover in April. Surprisingly, the fourth week of April is overwhelmed by a significant increment in the agricultural/arable area of the Po Valley (median increment around 10%), while a spatially homogeneous decrease, like the one observed in the second half of March, occurs only during the third week of the month (interquartile range between −40% and −20%). Note that during most of April 2020 the lockdown restrictions were the same of the third and fourth week of March. Looking at the maps of Figure 6 , it is apparent that the urbanized belt of the Po valley shows a persistent significant negative change during almost all the weeks of April. This is confirmed by the barplots of Figure 8 . Despite this variegated situation, the weekly distributions of the relative changes still are dominated by negative values with median values around −30%, with the exception of the fourth week.

In general, the estimates from the Sunday maps (see Appendix) are consistent with those obtained for working days. The most pronounced decline of the NO 2 concentrations occurs during the last two Sundays of March: we find a −40% median relative change in 2020 as compared to the same period in 2019. The same drop in the NO 2 concentrations is quantified for the second week of April (when Easter 2020 occurred).

Capturing the spatiotemporal variation in pollutant concentrations is of keen interest for the air pollution management's community. One of the reasons is that the effectiveness of environmental policies or interventions, or even exceptional events as the COVID-19 lockdown, must be evaluated across the whole study domain, in order to understand if the change is homogeneous in space or peculiar spatial patterns occur. This analysis is made even more difficult by the fact that air pollution dynamics strongly depends on weather conditions; thus, it is necessary to disentangle the effect of the adopted intervention by controlling for the known confounding factors. CTM can be used for this purpose. They can assess, with a very detailed spatial and temporal resolution, variations due to changes in the emission burden arising from interventions, like the implementation of air quality plans. However, some drawbacks exist: CTMs need very detailed input data (e.g., emission inventories, meteorology data), are complex to be run and rely on advanced IT infrastructure for their simulations. More importantly, they do not provide an estimate of the uncertainty of the final predictions which are also strongly dependent on the initial and boundary conditions. Moving to statistical models, it is straightforward to deal with uncertainty and spatiotemporal correlation structures, while keeping the computational costs at a reasonable level (especially if computationally efficient estimation methods, like the INLA-SPDE approach, are used).

In this paper we propose the use of a statistical Bayesian spatiotemporal model as a novel tool for assessing-in time but more importantly in space-the effectiveness of air quality policy interventions. As a case study, we focused on how the COVID-19 lockdown affected air pollution levels in the north/center of Italy, by means of the 2019/2020 NO 2 concentrations relative change, the main output of our modeling strategy. Even if our proposal does not explicitly simulate physical and chemical reactions in the atmosphere (like CTMs do), it has two key advantages: it is parsimonious in terms of data needs and, possibly more important, it naturally manages the uncertainty associated to the parameters and map outputs. The proposed method, applied to a real-world experiment, demonstrated its capability to provide a reliable picture of the temporal and spatial patterns of the NO 2 variations (compared to 2019), while accounting for the effect of meteorology.

In this study, the spatiotemporal variation of the NO 2 concentrations in March and April 2020, as compared with the same period in 2019, is illustrated through meteorology-normalized spatially continuous maps at weekly intervals, both for working days (Monday-Saturday) and weekends (Sunday). Our results show that during March and April 2020 the study domain was generally characterized by negative relative changes in the NO 2 concentrations with median values around −25%. Such estimate seems to be reasonably consistent with previous findings about NO 2 levels during the lockdown in Europe. However, the message from our output maps is richer than what a simple number can describe. First, a visual inspection of the maps shows that statistically significant reductions mainly occurred in the urban areas of the Po valley and Tuscany, while not significant variations persisted over the mountainous regions of Valle d'Aosta and the Autonoumous Provinces of Bolzano and Trento. Second, and more interestingly, our weekly analysis of NO 2 highlight that in March 2020 such reduction was synchronous with the lockdown measures adopted, while in April 2020 the concentrations, as compared to 2019, started to recover in some parts of the investigated area. To the authors of this study, this result comes not unexpected as it reflects the behavior of the input data seen in the parallel plot of Figure 2 . With the exception of Fassò et al. (2021) -who estimated for Lombardy region during the lockdown period a general decrease in NO 2 concentrations in metropolitan areas or for traffic-close stations along with a slight increase in NO 2 concentrations in rural and industrial sites-it comes as a surprise that none of the statistical studies we examined documented an increase in the NO 2 concentrations in April in the studied Italian domain. However, we are aware that previous findings are based on different models, data, and methodologies, so our results are not directly comparable, especially if we consider that all these studies do not provide continuous maps as an output.

It is worthwhile to observe that our modeling approach does not allow us to draw causal inference conclusions. However, given that the maps account for the weather effect, we can conclude that the reductions we observe across regions and weeks can be attributed to a factor different from the meteorology. Given that they occur in the same period of the restrictive measures, it is likely that the COVID-19 lockdown had an effect in the reduction of air pollution. Nevertheless, we can not exclude that other factors, not considered in the model, could have had an active role in the dynamics of NO 2 .

The model described in this paper takes inspiration from the model introduced in Fioravanti et al. (2021) for the spatiotemporal assessment of the daily log-transformed PM 10 concentrations in Italy. As here our objective was to describe the relative changes across the two considered years (2019 vs. 2020), in this study the target variable is the daily differences of the log-transformed NO 2 concentrations. As far as we know this modeling strategy has never been used before in the context of spatiotemporal models for intervention analysis. Other approaches could have been adopted, still in the context of spatiotemporal modeling. For example, the same problem could have also been tackled by, first, jointly modeling the 2019 and 2020 NO 2 concentrations and then computing the differences between the daily gridded predictions. However, this alternative solution has two important shortcomings: (1) it is not consistent with the AR(1) assumption, as our observations are not consecutive in time but are 1 year apart (March-April 2019 and 2020); 2) it is more time consuming as the size of the input dataset would be as twice as the size of the input dataset when daily differences are considered. Furthermore, the use of daily differences alleviates multicollinearity among regressors and allows to include more parameters in the regression equation. For this reason we believe that the spatiotemporal model we propose represents a valid solution to analyze pollutant concentrations measured by a network of monitoring stations before and after a given event of interest (e.g., the lockdown intervention). Interestingly, our model allows a straightforward spatial assessment of NO 2 variations through high-resolution estimates, since the posterior predictive distributions of NO 2 differences can virtually be derived for every pixel within the study domain. This is of keen interest for policymakers involved in intervention analysis or epidemiologist assessing population exposure change and gradients, since the assessment at urban hot spot, the between-city variability as well as the comparison among rural and urban context can be achieved at the finest resolution, with a reasonable computational effort. In this regard, our approach could be usefully replicated whenever the effects of an intervention aimed to reduce the pollutant emissions at local or regional scale have to be assessed. Last but not least, the proposed model is simple and can be easily extended to a larger spatial domain. This means that it could be also applied to other pollutants (e.g, PM, NO x , O 3 ) or, more generally, spatiotemporal phenomenon which are continuous in space and for which it is interesting to evaluate the change in consecutive years.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

ORCID Sara Martino https://orcid.org/0000-0003-4326-9029 Giorgio Cattani https://orcid.org/0000-0001- How to cite this article: Fioravanti, G., Cameletti, M., Martino, S., Cattani, G., & Pisoni, E. (2022) . A spatiotemporal analysis of NO 2 concentrations during the Italian 2020 COVID-19 lockdown. Environmetrics, e2723. https://doi.org/10.1002/env.2723 APPENDIX F I G U R E A1 Weekly maps of the NO 2 2019/2020 relative change (Sunday). Contour lines mark those areas where the relative change is statistically significant at the significance level of 0.05 F I G U R E A2 Weekly maps of the NO 2 2019/2020 relative change (Sunday). Contour lines mark those areas where the relative change is statistically significant at the significance level of 0.05

Environmental spatial heterogeneity of the impacts of covi-19 on the top-20 metropolitan cities of Asia-pacific

Environmental impacts of shifts in energy, emissions, and urban heat island during the COVID-19 lockdown across Pakistan

Does lockdown reduce air pollution? Evidence from 44 cities in northern China

Impacts of partial to complete COVID-19 lockdown on NO2 and PM2.5 levels in major urban cities of Europe and USA

Estimating lockdown induced European NO2 changes

Nitrogen dioxide reductions from satellite and surface observations during COVID-19 mitigation in Rome (Italy)

Impact of coronavirus outbreak on NO2 pollution assessed using TROPOMI and OMI observations

Spatio-temporal modelling of changes in air pollution exposure associated to the COVID-19 lockdown measures across Europe

Interventions to reduce ambient air pollution and their effects on health: An abridged Cochrane systematic review

The effect of corona virus lockdown on air pollution: Evidence from the city of Brescia in Lombardia Region (Italy)

Comparing spatio-temporal models for particulate matter in Piemonte

Spatio-temporal modeling of particulate matter concentration through the SPDE approach

European standard en 14211:2012 ambient air quality -Standard method for the measurement of the concentration of nitrogen dioxide and nitrogen monoxide by chemiluminescence

TROPOMI NO2 tropospheric column data: Regridding to 1 km grid-resolution and assessment of their consistency with in situ surface observations

Estimation of surface NO2 concentrations over Germany from TROPOMI satellite observations using a machine learning method

shiny: Web application framework for R

The misuse of colour in science communication

How have surface NO2 concentrations changed as a result of the UK's COVID-19 travel restrictions?

Air Quality in the Italian Northwestern Alps during Year 2020: Assessment of the COVID-19 lockdown effect from multi-technique observations and models

The effect of the COVID-19 lockdown on air quality in three Italian medium-sized cities

Recent advances in satellite mapping of global air quality: Evidences during COVID-19 pandemic

Air quality in Europe -2020 report

Statistical assessment of air quality interventions

Spatiotemporal variable selection and air quality impact assessment of COVID-19 lockdown

A model-based framework for air quality indices and population risk evaluation, with an application to the analysis of Scottish air quality data

Spatio-temporal modelling of PM10 daily concentrations in Italy using the SPDE approach

Did policies to abate atmospheric emissions from traffic have a positive effect in London? Environmental Pollution

Constructing priors that penalize the complexity of Gaussian random fields

The global impacts of COVID-19 lockdowns on urban air pollution: A critical review and recommendations

COVID-19 lockdown only partially alleviates health impacts of air pollution in Northern Italy

Using meteorological normalisation to detect interventions in air quality time series

COVID-19 lockdowns highlight a risk of increasing ozone pollution in European urban areas

A Bayesian LSTM model to evaluate the effects of air pollution control regulations in Beijing

The short-term impacts of COVID-19 lockdown on urban air pollution in China

The ERA5 global reanalysis

Separating the impact of gradual lockdown measures on air pollutants from seasonal variability

Review of the efficacy of low emission zones to improve urban air quality in European cities

Global impact of COVID-19 restrictions on the surface concentrations of nitrogen dioxide and ozone

Importance of satellite observations for high-resolution mapping of near-surface no2 by machine learning

Impact of lockdown on air quality over major cities across the globe during COVID-19 pandemic

A rigorous statistical framework for spatio-temporal pollution prediction and estimation of its long-term impact on health

Moving beyond noninformative priors: Why and how to choose weakly informative priors in Bayesian analyses

Bayesian spatial modelling with R-INLA

Lockdown measures and air quality: Evidence from Italian provinces

Statistical modeling of the early-stage

The impact of the COVID-19 emergency on local vehicular traffic and its consequences for the environment: The case of the city of Reggio Emilia (Italy)

Integrated Nested Laplace Approximations (INLA)

Impact of lockdown measures to combat COVID-19 on air quality over Western Europe

COVID-19 pandemic and environmental pollutionA blessing in disguise?

Italian regional SME policy responses

Early spring near-surface ozone in Europe during the COVID-19 shutdown: Meteorological effects outweigh emission changes

A physical explanation of the lognormality of pollutant concentrations

Meteorology-normalized impact of the COVID-19 lockdown upon NO 2 pollution in Spain

Spatial planning policy for sustainability: Analysis connecting land use and GHG emission in rural areas

Modeling the effect of COVID-19 lockdown on mobility and NO2 concentration in the Lombardy region

Impacts of the COVID-19 lockdown on air pollution at regional and urban background sites in northern Italy

A systematic literature review of the impact of COVID-19 lockdowns on air quality in China

COVID-19 and Italy: What next?

The urbanization run-up in Italy: From a qualitative goal in the boom decades to the present and future unsustainability

Spatio-temporal modeling of fine particulate matter

The Italian response to the COVID-19 crisis: Lessons learned and future direction in social development

Analysis of air quality during the COVID-19 pandemic lockdown in Naples (Italy)

Penalising model component complexity: A principled, practical approach to constructing priors

Quantifying the impact of the COVID-19 lockdown measures on nitrogen dioxide levels throughout

Penalised complexity priors for stationary autoregressive processes

Short-term change in air pollution following the COVID-19 state of emergency: A national analysis for the United States

Record high solar irradiance in Western Europe during first COVID-19 lockdown largely due to unusual weather

Air quality improvement from COVID-19 lockdown: Evidence from china

Estimating the spatiotemporal variation of NO2 concentration using an adaptive neuro-fuzzy inference system

Impacts of COVID-19 lockdown, spring festival and meteorology on the NO2 variations in early 2020 over China based on in-situ observations, satellite retrievals and model simulations

Effects of corona virus disease-19 control measures on air quality in North China

This research has been carried out within the framework of the PULVIRUS partnership agreement between ENEA, ISPRA, SNPA, and ISS. We wish to thank for their helpful contribution to the air quality and meteorological data pre-processing Raffaele Morelli, Walter Perconti and Arianna Trentini. In addition, we wish to thank all the Regional Environmental Protection Agencies (ARPA Valle d'Aosta, ARPA Piemonte, ARPA Lombardia, ARPA Veneto, APPA Bolzano, APPA Trento, ARPA Friuli, ARPAE Emilia Romagna, ARPA Toscana) that manage the air quality networks and provide fully validated NO 2 hourly data. This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

The authors declare no potential conflict of interests.