Understanding temporal and spatial changes of O3 or NO2 concentrations combining multivariate data analysis methods and air quality transport models

Understanding temporal and spatial changes of O 3 or NO 2 concentrations combining multivariate data analysis methods and air quality transport models Stefan Platikanov a , Marta Terrado b , María Teresa Pay c ,AlbertSoret b , Romà Tauler a, ⁎ a Department of Environmental Chemistry, IDAEA-CSIC, Jordi Girona, 18-26, 08034 Barcelona, Spain b Earth Sciences Department, Barcelona Supercomputing Center (BSC), Jordi Girona, 31, 08034 Barcelona, Spain c Department of Genetics, Microbiology and Statistics, University of Barcelona, Faculty of Biology, Diagonal, 643, 08028 Barcelona, Spain HIGHLIGHTS •CALIOPE O 3 and NO 2 concentration forecasts at different spatial resolutions in Catalonia (Spain). •Experimental and CALIOPE O 3 and NO 2 concentration forecasts are analyzed by the MCR-ALS method •The accuracy of CALIOPE forecasts depends on the spatial resolution, time and contaminant. •MCR–ALS describes the most meaningful contributions to NO 2 and O 3 concentration changes. •Sunlight, season, traffic and local station environment rule different patterns. GRAPHICAL ABSTRACT abstractarticle info Article history: Received 19 July 2021 Received in revised form 7 October 2021 Accepted 7 October 2021 Available online 13 October 2021 Editor: Philip K. Hopke The application of the multivariate curve resolution method to the analysis of temporal and spatial data variability of hourly measured O 3 and NO 2 concentrations at nineteen air quality monitoring stations across Catalonia, Spain, during 2015 is shown. Data analyzed included ground-based experimental measurements and predicted concentrations by the CALIOPE air quality modelling system at three horizontal resolutions (Europe at 12 × 12 km 2 , Iberian Peninsula at 4 × 4 km 2 and Catalonia at 1 × 1 km 2 ). Results obtained in the analysis of these different data sets allowed a better understanding of O 3 and NO 2 concentration changes as a sum of a small number of different contributions related to daily sunlight radiation, seasonal dynamics, traffic emission patterns, and local station environments (urban, suburban and rural). The evaluation of O 3 and NO 2 concentrations predicted by the CALIOPE system revealed some differences among data sets at different spatial resolutions. NO 2 predictions, showed in general a better performance than O 3 predictions for the three model resolutions, specially at urban stations. Our results confirmed that the application of the trilinearity constraint during the multivariate curve resolution factor analysis decomposition of the analyzed data sets is a useful tool to facilitate the understanding of the resolved variability sources. © 2021 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http:// creativecommons.org/licenses/by/4.0/). Keywords: Air quality modelling Chemometrics MCR-ALS Ozone Nitrogen dioxide CALIOPE system Science of the Total Environment 806 (2022) 150923 ⁎Corresponding author. E-mail address: [email protected] (R. Tauler). https://doi.org/10.1016/j.scitotenv.2021.150923 0048-9697/© 2021 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Contents lists available at ScienceDirect Science of the Total Environment journal homepage: www.elsevier.com/locate/scitotenv 1. Introduction Air pollution is a serious threat to both human and environmental health, being currently one of the most pressing challenges for cities in Europe and around the globe (EEA, 2019;WHO, 2016). In 2016, 91% of the world population was living in places where the World Health Organization (WHO) air quality guideline levels were not met. Having information about air pollution contributions and dynamics is essential to design adequate policies to improve air quality and reduce the negative impacts of pollution on health (European Environmental Agency (EEA), 2011;Nieuwenhuijsen, 2018). Monitoring and modelling of air pollutants are essential tasks to evaluate the impact from the continuous increase of human activities on the environment and public health. Therefore, current research efforts are mainly focused on the improvement of monitoring and modelling systems, and on accurately forecasting the behavior of hazardous air pollutants in order to understand their origin, transport, geographical distribution and time evolution (Schaap et al., 2015;Pay et al., 2019;Massagué et al., 2019). Ozone (O 3 ) and nitrogen dioxide (NO 2 ) are among the major air pollutants directly associated with negative effects on the human health (Bell et al., 2004;WHO, 2013;VicedoCabrera et al., 2020). Roughly 90% of O 3 is distributed in the stratosphere and around 10% in the troposphere. The origins of O 3 in the troposphere include its migration from the stratosphere and its formation from various photochemical reactions among carbon monoxide (CO), peroxy radicals (generated by the photochemical oxidation of volatile organic compounds, VOCs) and nitrogen oxides (NO x )(Crutzen, 1974). High temperatures, high solar radiation and low precipitation also favor positively the formation of tropospheric O 3 (Otero et al., 2016). Photochemical and oxidation reactions take place specially when the sunlight interacts with NO x and VOCs, mostly derived from anthropogenic emissions, thus converting the industrial and traffic emissions in additional contributions of O 3 formation. In these urban and industrialized areas, the concentrations of O 3 sharply peak at sunny middays, because of the rapid interconversion of O 3 and NO x . Conversely, when the impact of sunlight decreases in the evenings, destruction of O 3 occurs because of its reaction with NO from different activities, such as the evening rush hour car traffic. The circulation of air masses also controls the short and long-range transport of O 3 , affecting its lifetime in the atmosphere and its levels in rural and remote areas (Monks et al., 2015;Pay et al., 2019). For example, the transport of precursors emitted in urban and industrialized areas has been shown to cause O 3 production downwind (Holloway et al., 2003; Querol et al., 2017). NO 2 is also considered to be an important trace gas in the atmosphere. Apart from its participation in the photocatalytic formation of O 3 in the troposphere, NO 2 can also produce secondary aerosols and acid rain (Chan et al., 2015). Generally, NO 2 can have both natural and anthropogenic origins, and most of the anthropogenic contributions are related to fossil fuel combustion, biomass burning, and various types of industrial emissions (Tack et al., 2015;Gratsea et al., 2016;Chan and Chan, 2017). Across industrial and urban areas, the air pollution is largely due to emissions of pollutants from motor vehicles, industrial plants and generation of energy (i.e. power plants) (Khan et al., 2018). Sudden increases in pollutants’concentrations are correlated with the intensive morning and evening car traffic. The European Directive on Ambient Air Quality and Clean Air (European Commission, 2008) stands for air quality monitoring and modelling as useful tools to understand the dynamics of air pollutants, analyze and forecast air quality, develop plans to reduce emissions and alert the population when health-related episodes occur (EEA, 2011). In line with these objectives, the CALIOPE air quality forecasting system (Baldasano et al., 2008a, 2008b;http:// www.bsc.es/caliope/es, last access: 15 March 2021) has been developed to forecast the air quality in Spain. CALIOPE, which runs on high-performance computing platforms, uses multi-model inputs about meteorology, atmospheric composition and emissions, to predict the air quality at three different spatial resolutions: 12 × 12 km 2 for the European domain (EU12); 4 × 4 km 2 ,forthe Spanish Iberian Peninsula domain (IP4); and 1 × 1 km 2 for the Catalonian domain (CAT1); in all cases with 1 h time resolution for the next 48 h. Air quality modelling requires the use of physicochemical parametric models strongly dependent on accurate prior knowledge about the atmospheric processes that involve air pollutants’transport dispersion and reaction. In the CALIOPE system, all possible local contributions of pollution for a given territorial domain, together with the possible influence from more distant geographical zones - transboundary transport models –should be explicitly considered and, in general, need long computation time, especially when they are run at high spatial resolutions. Various studies have been conducted with the objective to compare the pollutants’concentrations predicted by CALIOPE with the real experimental concentrations measured by the ground-based monitoring stations (Pay et al., 2010, 2012;Baldasano et al., 2011;Aguilera et al., 2013). Satisfactory results have been observed in reproducing the temporal and spatial distribution of various pollutants, including O 3 and NO 2 , at different resolutions (Pay et al., 2014). In this work, a multivariate data analysis approach based on chemometrics soft modelling methods is used to evaluate the CALIOPE forecasts of O 3 and NO 2 concentrations. Unlike the more traditional parametric approaches, based on complex physicochemical and meteorological forecasting models, the proposed chemometrics approach offers the possibility to directly investigate the hidden information present in the experimental data without the assumption of a physical model. This is achieved by the application of a multilinearity constraint during the factor analysis decomposition of the experimental data arranged in multiway/multimode data structures (Norris et al., 2014;Belis et al., 2019;Malik and Tauler, 2013;Alier et al., 2009; Dadashi et al., 2020). In this work, we apply a chemometric analysis to both the experimental O 3 and NO 2 concentrations and to the concentrations predicted by the CALIOPE system at three different spatial resolutions (EU12, IP4 and CAT1) during the year 2015 at 19 air quality monitoring stationsacross Catalonia. The comparison of the results of the analysis of experimental and CALIOPE predicted data allowed us to identify the applications for which the CALIOPE system is already fit-for-use and the aspects that can be improved for a better performance. This study upgrades previous investigations for Europe (Schaap et al., 2015) and for the Iberian Peninsula (Pay et al., 2014), as it includes a full year simulated data at higher spatial resolution (1 × 1 km 2 ) for the whole territory of Catalonia. 2. Datasets 2.1. Experimental dataset Hourly O 3 and NO 2 concentration data were obtained from the Catalan air quality monitoring network (named Xarxa de Vigilància i Previsió de la Contaminació Atmosfèrica, XVPCA), which provides a relatively dense geographical coverage of the Catalan territory. The network contributes to the European Environment Information and Observation Network (EIONET; https://www.eionet.europa.eu/, last access: April 2019). For 2015, O 3 and NO 2 measurements were obtained from 19 monitoring stations across Catalonia with a temporal coverage higher than 85% on an hourly basis (Fig. 1). As defined by the European Environment Agency (EU, 2011), depending on the distribution/density of buildings, the area surrounding the station is classified as urban (6 stations, u1-u6 in Fig. 1), suburban (5 stations, s7-s11) and rural (8 stations, r12r19), from coastal to mountain sites. All measurements were recorded in UTC time. S. Platikanov, M. Terrado, M.T. Pay et al. Science of the Total Environment 806 (2022) 150923 2 2.2. CALIOPE predicted datasets The CALIOPE air quality forecasting system (http://www.bsc.es/ caliope/es, last access: 15 March 2021) was used to predict the NO 2 and O 3 concentrations over Catalonia during year 2015 using three spatial grid resolutions. The system, described in previous works (Baldasano et al., 2008a, 2011;Pay et al., 2010), is based in numerical modelling including a full and precise description of the processes involved for gas-phase and particulate matter modelling, both important factors to obtain accurate results of air pollutant concentrations in a complex region as Spain. CALIOPE is composed of the HERMES v2.0 emission model (Guevara et al., 2013), the WRF-ARW v3.6 meteorological model (Skamarock and Klemp, 2008), the CMAQ v5.0.2 chemical transport model (Byun and Schere, 2006) and the BSC-DREAM8bv2 mineral dust model (Basart et al., 2012). CALIOPE first runs over Europe at a 12 × 12 km 2 horizontal resolution (EU12 domain), then over the Iberian Peninsula at a 4 × 4 km 2 resolution (IP4 domain) and finally over Catalonia at a 1 × 1 km 2 resolution (CAT1 domain). The WRF-ARW system provides the meteorological fields required by the CMAQ chemical transport model. In the present work, this system has been configured with 38 sigma layers up to 50 hPa to resolve the troposphere-stratosphere exchanges properly. The planetary boundary layer (PBL) is characterized by approximately 11 layers, and the bottom layer's depth is 39 m. The EU12 domain uses meteorological initial and boundary conditions from the final analyses provided by the National Centers of Environmental Prediction (FNL/NCEP) at a 0.5° × 0.5° resolution. The first 12 h of each meteorological run are treated as cold start, and the next 23 h are provided to the chemical transport model. The CMAQ chemical transport model combines current knowledge in atmospheric science and air quality modelling with multiprocessor computing techniques in an open-source framework to deliver concentration estimates of the main air pollutants. CMAQ v5.0.2 uses the CB05 chemical mechanism (Yarwood et al., 2005), the AERO5 for aerosol modelling, and the in-line photolysis calculation. CALIOPE considers desert dust contribution by means of the BSC-DREAM8bv2, which runs offline at a 0.5° × 0.5° resolution covering Europe, North Africa and the Middle East. The emissions model HERMESv2.0 uses information and the state-ofthe-art methodologies for the estimation of emissions (Baldasano et al., 2008b;Guevara et al., 2013). For the non-Spanish countries included the IP4 domain, HERMESv2.0 processes the original annual European Monitoring and Evaluation Programme (EMEP) gridded emissions (50 × 50 km 2 ), performing a Selected Nomenclature for Air Pollution (SNAP) sector-dependent spatial and temporal (1 h) disaggregation as well as a speciation treatment (top-down approach). For the IP4 and CAT1 domains, HERMESv2.0 uses a combination of local information and state-of-the-art methodologies (mostly bottom-up approaches) to estimate anthropogenic and biogenic emissions at a high spatial (1 × 1 km 2 ) and temporal (1 h) resolution over the whole Spanish territory, which for the case of the IP4 are then aggregated to a 4 × 4 km 2 resolution (Guevara et al., 2013). Hourly simulations of NO 2 and O 3 concentrations for year 2015 covering the territory of Catalonia from the EU12, IP4 and CAT1 domains were used for this study. For each domain, the grid cells coinciding with the location of the monitoring stations were used for comparison of simulated results with observations. 2.3. Data arrangement During 2015, each monitoring station (k = 1,..,19) produced two long data vectors, (d k ), one for O 3 and another for NO 2 , with 24-h recordings for the 365 days (with dimensions 1 × 8760). These row vectors for every station were concatenated one on the top of the other giving the data matrix D* with size: 19 stations × 8760 recordings, as shown in Fig. 2a. In a second type of data arrangement (Fig. 2b), the O 3 or NO 2 data vectors (d k ) for every data station (k = 1,.0.19) were first individually arranged in a data matrix D k with 365 rows (days) and 24 columns (hours). In the simultaneous analysis of the data from the 19 stations, these individual matrices D k , were vertically concatenated in a columnwise augmented data matrix, D aug of size 6935 rows (365 days × 19 stations) and 24 columns (hours). These data arrangements were applied to the experimental and to the CALIOPE predicted data at the three spatial resolutions (EU12, IP4 and CAT1). Less than 5% of the O 3 and NO 2 measured concentrations were missing in the experimental Fig. 1. Location of the EIONET O 3 and NO 2 monitoring stations in Catalonia, Spain. Different icons represent the three environment defined by the EEA. Urban stations: u1 - Cuitadella district (Barcelona), u2 - Vall d’Hebron district(Barcelona), u3Tarragona, u4Mataró, u5 - Sant Cugat, u6Rubí.Suburban stations: s7 - Manlleu, s8 - Amposta, s 9 - Gavà, s10 - Vilafranca de Penedès, s11 - Viladecans. Rural stations: r12 - Berga, r13 - Tona, r14 - Juneda, r15 - Bellver, r16 - Cap de Creus, r17 - Els Torms, r18 - Montsec, r19 - Montseny. S. Platikanov, M. Terrado, M.T. Pay et al. Science of the Total Environment 806 (2022) 150923 3 data. For the CALIOPE data sets, missing values were below 1%. In this work, missing values in the individual data matrices D k were replaced by the mean (average) in the same column of the data matrix. Since this number was relatively small, the results will not be affected significantly by the imputation method and more rigorous techniques for data missing imputation (Stanimirova, 2013) were not required. Additionally, the Whittaker smoother (Eilers, 2003) was used for offset correction in O 3 and NO 2 data matrices. 3. Methods 3.1. PCA and bilinear MCR-ALS Different chemometric methods have been proposed in the literature for the analysis of environmental data. Principal Component Analysis (PCA) (Jolliffe, 2002) is the method most widely used to analyzethe variance sources and patterns of variation of environmental datasets. Moreover, the US Environmental Protection Agency (EPA) recommends the application of Positive Matrix Factorization (PMF) for factor analysis in air quality studies (Hopke, 2008;Norris et al., 2014). Similarly to PMF, Multivariate Curve Resolution Alternating Least Squares (MCR-ALS) is a frequently used method in spectrochemical mixture data analysis, which can also be easily extended to the analysis of environmental source apportionment data sets (Alier and Tauler, 2013;Alier et al., 2011). MCR-ALS is a flexible soft-modelling factor analysis tool which allows for the application of natural constraints (see below) and can be easily adapted to the analysis of complex multiway (multimode) data structures, such as three-way and four-way environmental data sets using trilinear and quadrilinear models (Tauler et al., 1998; Smilde et al., 2004;Malik and Tauler, 2013;Marín-García and Tauler, 2020). In this work, PCA and MCR-ALS methods are used to reveal and compare the information present in the experimental and CALIOPE predicted O 3 and NO 2 data sets. Both, PCA and MCR-ALS methods, are based on a bilinear model that performs the factor decomposition of a two-way data set (i.e. a data table or a data matrix). Eq. (1) summarizes this bilinear model in its element-wise form, while Eq. (2) presents the same model in a matrix linear algebra format: dij ¼∑N n¼1xi,nyj,n þei,j i¼1, ::IrowsðÞ,j¼1, ...J columnsðÞð1Þ D¼XY TþEð2Þ In Eq. (1), the individual data values, d i,j elements (in this case the O 3 or NO 2 concentration values) are decomposed as the sum of a number of components (contributions), n= 1,..N, each one of them defined by the product of two factors, x i,n (scores) and y j,n (loadings). In addition, the term e i,j is the residual part of d i,j ,which cannot be explained by these N components and accounts for the experimental noise and uncertainties. In Eq. (2), the data matrix, D,of dimensions IxJ is decomposed into the scores factor matrix X(IxN) and the loadings factor matrix, Y T (NxJ). The number of components, N, is selected to explain as much as possible the data variance, while the unexplained small contributions of data variance and experimental noise are in E. InthecaseofPCA,thecomponentsareobtainedundertheconstraints of orthogonality and maximum variance, and they are ranked in descending order of explained variance (see Jolliffe, 2002 for more details). Thus, a reduced set of N (N ≪JorI)principalcomponents (PCs) are selected, which contain most of the relevant variance (not noise) in the original data matrix. The constraints used in PCA provide mathematical unique solutions, but they imply negative values in the resolved component profiles, which for example in the case of O 3 and NO 2 concentrations have not a direct physical meaning. Fig. 2. Data arrangement and application of (a) PCA and (b) MCR-ALS bilinear and trilinear. Hourly O 3 and NO 2 concentrations experimentally observed during one year were arranged in long row vectors, d k for the different monitoring stations(k = 1,... 19), in the wide data matrix D* with 19 rows and 8760 columns, (as shown inFig. 2a), in multiple individual station data matrices D k with 365 rows and 24 columns and as column-wise augmented data matrix D aug with 365 × 19 = 8760 rows and 24 columns, as shown in Fig. 2b). PCA, MCR bilinear and MCR-ALS trilinear decompositions are also shown graphically (see Sections 3.1 and 3.2 and Eqs. (1)–(4)). S. Platikanov, M. Terrado, M.T. Pay et al. Science of the Total Environment 806 (2022) 150923 4 Multivariate Curve Resolution (Tauler, 1995) performs a similar bilinear model decomposition as shown in Eq. (2) but in this case different constraints are applied to provide more interpretable and meaningful factor profiles. The bilinear model factor decomposition is performed using an alternating least squares (ALS) optimization under constraints fulfilled naturally by the real data, for example, non-negativity constraints (Bro and de Jong, 1997;de Juan and Tauler, 2003). In the MCR-ALS method, an initial selection of the number of components is required. This can be gathered from PCA or singular value decomposition (Golub and van Loan, 1989) of the analyzed data set, by inspection of the sizes of the singular (or eigen) values of the analyzed data matrix. This selection can also be made by taking into account the amount of explained variance by the successive components, thus avoiding components associated with noise, usually with the smaller singular (eigen) values. Models with different numbers of components can be tested and a final decision is taken considering the data fit and the shapes of the resolved profiles. The ALS algorithm also needs initial estimates of either Xor Y T factor matrices. These initial estimates can be obtained from the ‘more ‘dissimilar’rows or columns of the original data matrix (Windig and Guilment, 1991). Eq. (2) for Dis solved iteratively using an ALS optimization, which updates the solutions (vector profiles in X and Y T matrices) until they fit the data optimally and fulfill the proposed constraints. The profiles in Xand Y T matrices are constrained to be nonnegative and they are directly interpretable in physical terms. The bilinear data modelling used in Eqs. (1) and (2) can be applied to the experimental and CALIOPE predicted data sets explained in the previous data section (data matrices D* and D aug )andinFig. 2. 3.2. Trilinear MCR-ALS The bilinear factor decomposition model given for a two-way data matrix can be extended to a three-way dataset, D, using a trilinear model. This model can be expressed individually for every data value as given by Eq. (3). dijk ¼∑N n¼1xi,nyj,nzk,n þei,j,k ð3Þ where d i,j,k are the individual data values (concentrations of O 3 or NO 2 ) in the three experimental data modes: the day of the year i=1,…,365, the hour of the day j = 1, …,24, and the station k = 1,..,19, (urban, u, suburban, s, and rural, r). Data are modeled as the sum of a number of components (contributions), n= 1,..N, defined by the product of three factors: x i,n ,y j,n , and z k,n, These factors are related with the three data modes respectively. e i,j,k is the part of d ijk not explained by the contribution of these N components. This trilinear model can be written in a matrix form using the decomposition of every individual D k data slice (every individual matrix D k )inEq.(4). Dk¼XZkYTþEkfor k ¼1, :::19 ð4Þ Under the trilinear model, all individual data matrices, D k (I,J) are simultaneously decomposed with the same number of components N and the same daily, X(I,N) and hourly Y T (N,J) profiles. Thus, they differ only in a diagonal matrix Z k (N,N) different for every one of the k = 1,... 19 stations (station profiles), which gives the relative amounts of the N components in every data matrix (station), D k . These N diagonal elements of the Z k canalsobegroupedinthe third factor matrix Z(K,N). Therefore, the proposed trilinear model takes advantage of the natural structure of the analyzed data sets, especially in relation to their different temporal modes (i.e. hourly, daily, seasonal) and to the different type of monitoring stations analyzed simultaneously. Fig. 2b shows the practical implementation of the trilinear model in the MCR-ALS analysis of a three-way data set. Nineteen individual data sets (one per station), were arranged in the column-wise augmented data matrix D aug of dimensions 365 days × 19 stations = 6935 rows and 24 hourly measures in columns. The application of the trilinearity constraint implies that the augmented profiles of every component n, x n aug , of dimensions 6935 × 1, is first refolded (step 1) in a data matrix X n aug of dimensions 365 days × 19 stations. Then this matrix is decomposed by SVD considering only the first component into the product of two vector profiles, x n (365 × 1) describing the 365 daily profiles, and z n (19 × 1) describing differences among the different stations. The application of this trilinearity constraint implies that for every component, the daily changes are described by the same single x n vector profile which changes station by station by the corresponding scalar value in z n . See previous works for a more detailed description of the algorithm used for the practical implementation of the trilinearity constraint in MCR-ALS (Tauler et al., 1998;de Juan and Tauler, 2001). MCR-ALS allows for the implementation of the trilinearity (or multililinearity) constraint separately for every component and it also allows for the relaxation of the synchronization of the component profiles, permitting therefore building mixed bilinear-trilinear models. It is worth also to mention that the application of the trilinearity constraint in MCR-ALS converges fast. MCR-ALS with the trilinearity constraint gives similar results to PARAFAC as we have shown in previous works (Tauler, 1995;Tauler et al., 1998;De Juan and Tauler, 2001), and solves the rotational ambiguity problem associated to MCR-ALS bilinear. The evaluation of the MCR-ALS fitting results is performed calculating the explained data variances (R 2 )usingEq.(5), R2¼100 1− ∑19 k¼1∑365 i¼1∑24 j¼1dijk −b dijk  2 ∑19 k∑365 i¼1∑24 j¼1d2 ijk 0 B @1 C Að5Þ where d ijk are the experimental or CALIOPE predicted O 3 (or NO 2 ) concentrations, and b dijk are the corresponding calculated values by MCR-ALS using either the bilinear (Eqs. (1)–(2)) or the trilinear model (Eqs. (3)–(4)). 3.3. Accuracy of CALIOPE predicted data The accuracy of the CALIOPE predicted data at three spatial resolutions can be assessed with the calculation of the Root Mean Square Error (RMSE) between the experimental and the CALIOPE model predicted data, averaging them according to Eqs. (6) and (7), respectively. RMSEdaily ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ∑19 k¼1∑24 j¼1dijk −b dijk  2 19x24 v u u tð6Þ RMSEhourly ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ∑19 k¼1∑365 i¼1dijk −b dijk  2 19 x 365 v u u tð7Þ where d ijk are the experimental or CALIOPE predicted O 3 or NO 2 concentrations at day i=1,…,365,hour j = 1, …,24 and station k=1,…,19 and b dijk are the corresponding experimental concentrations. RMSE daily and RMSE hourly profiles (Eqs. (6) and (7)) are the representative daily and hourly mean profiles averaged for all stations. To further compare and summarize the prediction accuracy, additional mean scalar values were calculated for all the data values averaging all data values per day, hour and station. RMSEavg ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ∑19 k¼1∑365 i¼1∑24 j¼1dijk −b dijk  2 19 x 365 x 24 v u u tð8Þ This gives a single scalar value RMSE avg for all stations or they can grouped also according to the urban, suburban and rural types for S. Platikanov, M. Terrado, M.T. Pay et al. Science of the Total Environment 806 (2022) 150923 5 comparison purposes (see experimental data set Section 2.1,Fig. 1 and Table 2 of results). 3.4. Chemometrics software All data were transferred to MATLAB computer and visualization environment (The Mathworks, Natick, USA) for their analysis. PCA was performed using PLS 8.5 Toolbox software (Eigenvector Research Ltd., Manson, WA, USA) and MCR-ALS analysis of the different data sets was performed using the MCR toolbox (Jaumot et al., 2015;http:// www.mcrals.info/). 4. Results and discussion In the Supplementary Material section, a brief description of the experimental O 3 and NO 2 concentration data is given. Supplementary Fig. S1 shows the average profiles of O 3 and NO 2 experimental concentrations measured over 365 days (a and b), over the 24 h of the day (c and d), and for the different stations (e and f). Values from urban, suburban and rural stations have been marked in different colors. Supplementary Fig. S2 shows the average concentration profiles of O 3 (a) and NO 2 (b) experimental (black) and predicted by CALIOPE (EU12 green, IP4 red and CAT1 blue) over the 24 h. 4.1. PCA and MCR-ALS results PCA results of the mean-centered O 3 and NO 2 column-wise augmented experimental data matrices (see Fig. 2a) of the 24 h × 365 days concentration values at the 19 stations in 2015, D* O3 and D* NO2, are given in Fig. 3. PC1 and PC2 of the O 3 matrix explain 52.2% and 9.9% of the total data variance, respectively. In the PC1vsPC2 scores plot (Fig. 3a), the 19 stations can be well distinguished in three groups according to the station type. The urban (blue triangles) and rural (green diamonds) type of stations are significantly separated one from another, while the suburban stations (red squares) occupy intermediate positions between them. This separation agrees well with the stations classification given by the local environmental agencies. In the loadings plot (Fig. 3c), the different symbols and colors refer to the different seasons. Loadings on PC1 are mostly positive in all seasons with a homogeneous distribution. In the case of NO 2 (Fig. 3b and d), PC1 and PC2 explain 67.8% and 5.7% of the total data variance. The 19 stations are again distinguished in three different groups for NO 2 concentrations (Fig. 3b). The dominating group in this plot corresponds to urban stations (blue triangles), showing a wider dispersion, while suburban and rural stations are distributed around the axis origin. No clear seasonal pattern can be observed in the loadings plot (see Fig. 3d). In summary, PCA results in Fig. 3 show some differences in the trends of O 3 and NO 2 concentrations according to the type of station whereas these differences are smaller among the different seasons. In a similar way, PCA was also applied to the three CALIOPE predicted data sets at different resolutions, D* EU12,O3 and D* EU12,NO2 ;D* IP4,O3 and D* IP4,NO2 and D* CAT1,O3 and D* CAT1,NO2 . The results of PCA analysis of CALIOPE predicted data are given in more detail in the Supplementary material. Table 1 summarizes the explained data variances obtained in the analysis of experimental and CALIOPE predicted data by MCR-ALS Fig. 3. PCA results of the mean-centered O 3 and NO 2 2015 experimental data matrix D*. Upper panels show PC1 vs PC2 scores for the analysis of O 3 (a) and NO 2 (b) colored by station type: urban, suburban and rural. Lower panels show PC1 vs PC2 loadings for the analysis of O 3 (c) and NO 2 (d), colored by season: autumn, winter, spring, summer. S. Platikanov, M. Terrado, M.T. Pay et al. Science of the Total Environment 806 (2022) 150923 6 bilinear and trilinear (i.e. with trilinearity constraint) using three components. In the MCR-ALS results, all three components explained relatively similar amounts of data variance for O 3 and NO 2. With the same number of components (three), the explained variances by the bilinear (98.4%, 93.6%) and trilinear (93.9%, 85.7%) MCR-ALS models for O 3 and NO 2 were rather similar. This suggests the convenience of applying the trilinearity constraint for this type of data, which has also the advantage of providing a more simplified interpretation of the results. MCR-ALS models for all O 3 datasets showed higher amounts of overlapping variance than the models for NO 2 data sets, suggesting a more complex behavior in the case of O 3 concentrations. Overlapping information for MCR-ALS is reduced with the application of the trilinearity constraint. The results of MCR-ALS analysis of the experimental data and of the data predicted by the EU12, IP4 and CAT1 CALIOPE models are given in Fig. 4 (O 3 )andFig. 5 (NO 2 ). According to the MCR-ALS first component (C1), O 3 daily profiles (Fig. 4a) have a similar shape in all cases, although for the experimental data analyzed by MCR-ALS (in black in Fig. 4a), it has a higher intensity in the warmer seasons. C1 O 3 hourly profiles (Fig. 4b) are also similar in shape but a time shift of two hours is observed for the experimental data profile (peak maximum around 20:00 h, in the evening) compared with the homologous profiles for the CALIOPE predicted data (peak maximum around 18:00 h). C1 O3 station profiles for IP4 and CAT1 (Fig. 4c) resemble very much the profile obtained in the analysis of the experimental data, showing that this C1 contribution is more important for rural stations. However, the C1 O 3 EU12 profile in Fig. 4c has a rather constant contribution for all the monitoring stations, without differentiating rural from the other types of stations. This indicates that the EU12 data had not enough spatial resolution to capture the local changes among different monitoring stations. C2 O 3 daily profiles in Fig. 4d have similar shapes for all data sets but with rather different intensities. In this case, the EU12 O 3 daily profile has the highest intensity, while the IP4 and CAT1 profiles have lower intensities than the experimental profile. C2 O 3 EU12 hourly profile (Fig. 4e) is different to the homologous experimental, IP4 and CAT1 profiles. While hourly profiles for IP4 and CAT1 have peak maxima at 02:00 h and 04:00 h respectively, the EU12 hourly profile keeps aflat shape. Again, in this case, a 2 h’time shift is observed for the experimental profile compared to IP4 and CAT1 profiles. C2 O 3 station profiles in Fig. 4f show no significant differences among the four different data sets. In all cases, the C2 O 3 station profile has higher values for the rural stations, with a more pronounced O 3 night peak. Finally, C3 O 3 daily profile (Fig. 4g) has a more pronounced peak shape in the experimental profile than in the CALIOPE profiles. C3 O 3 hourly profiles in Fig. 4h show a pronounced peak maximum around noon, shifted again by 2 h in the case of the experimental data. This O 3 profile clearly shows the ozone formation due to the increasing solar radiation at noon and afternoon time. C3 O 3 station profiles in Fig. 4i show similar trends for the different monitoring stations, decreasing from urban to rural stations. Fig. 5 shows the comparison of the daily, hourly and station MCRALS profiles obtained in the investigation of NO 2 data sets (experimental, EU12, IP4 and CAT1). C1 NO 2 daily profiles in Fig. 5a show some differences in their shape throughout the year. Whereas in the case of the experimental data (black line), the profile has a pronounced U-shape with their maxima in autumn and winter (increase of NO 2 ), the C1 NO 2 annual daily profiles for the CALIOPE predicted data are almost flat lines. The C1 hourly (Fig. 5b) and station (Fig. 5c) profiles are very similar in shape and intensity in all cases. In this case, no pronounced time shift was observed and the evening peak maximum of the C1 NO 2 hourly profile coincides around 21:00 h in all cases. The C1 NO 2 station profiles are also coincident, describing similar urban and suburban background NO 2 pollution. C2 NO 2 daily profiles (Fig. 5d) are also similar for all the investigated data sets. Only the C2 NO 2 daily profile obtained for the experimental data presents slightly more pronounced changes, mainly in spring. Hourly and station C2 NO 2 profiles (Fig. 5e and f) are practically identical in all cases, suggesting a great similarity in the results obtained in the MCR-ALS analysis of the experimental and CALIOPE predicted data. C3 NO 2 daily profiles in Fig. 5g are also coincident in all cases. As for C1 and C2 NO 2 daily profiles, a similar seasonal trend with a U shape was observed, with their peak minima in the summer and their peak maxima in autumn and winter. A small difference can be observed for the C3 NO 2 hourly CALIOPE profiles (Fig. 5h) compared to the experimental data. A double humped plateau with two peaks and a valley at 14:00 h is more pronounced for the profile obtained from experimental data, and it is not so well defined for the profiles obtained from the three CALIOPE models. In contrast, the latter have a more pronounced second peak maximum at 17:00 h, which is higher than the first peak at 09:00 h. A more detailed comparison of the hourly profiles for O 3 and NO 2 detected in the MCR-ALS analysis of the experimental data, shows that Table 1 Summary of the PCA and MCR-ALS (bilinear and trilinear) explained variances (in %) for the different data sets (experimental, EU12, IP4 and CAT1), using three components in all cases. Data MCR-ALS Bilinear MCR-ALS Trilinear O 3 Experimental Total a 98.4 Total a 93.9 C1 49.6 C1 59.6 C2 44.8 C2 42.1 C3 43.3 C3 33.6 Sum b 137.7 Sum b 135.3 O 3 EU12 model Total a 99.5 Total a 98.5 C1 61.3 C1 30.2 C2 48.9 C2 90.7 C3 36.9 C3 8.7 Sum b 147.1 Sum b 129.6 O 3 IP4 model Total a 99.3 Total a 97.7 C1 56.1 C1 63.1 C2 52.6 C2 56.1 C3 36.5 C3 26.2 Sum b 145.2 Sum b 145.4 O 3 CAT1 model Total a 99.2 Total a 97.3 C1 53.3 C1 68.6 C2 54.3 C2 47.3 C3 34.4 C3 14.3 Sum b 142.0 Sum b 130.2 NO 2 Experimental Total a 93.6 Total a 85.7 C1 43.8 C1 41.6 C2 44.9 C2 38.9 C3 37.1 C3 36.9 Sum b 125.8 Sum b 117.4 NO 2 EU12 model Total a 96.8 Total a 92.8 C1 37.3 C1 34.8 C2 46.3 C2 42.3 C3 41.9 C3 44.1 Sum b 125.5 Sum b 121.2 NO 2 IP4 model Total a 95.4 Total a 89.7 C1 50.3 C1 41.6 C2 34.0 C2 32.5 C3 39.4 C3 45.8 Sum b 123.7 Sum b 119.9 NO 2 CAT1 model Total a 94.4 Total a 86.7 C1 47.3 C1 40.6 C2 36.0 C2 34.1 C3 37.5 C3 37.4 Sum b 120.8 Sum b 112.1 a Total states for the percentage ofexplained variance by MCR-ALS bilinear and trilinear models with all components. b Sum states for the sum of the percentages of explained variances by the individual MCR-ALS components C1, C2 and C3. This Sum can be above 100% due to overlapping of information among the different MCR-ALS components. S. Platikanov, M. Terrado, M.T. Pay et al. Science of the Total Environment 806 (2022) 150923 7 the first maximum of NO 2 profile at 08:00 h (see C2 Fig. 5e) comes before the C3 O 3 hourly profile starts to increase in the morning (see Fig. 4h). This confirms that the sunlight was triggering the reaction between NO (from morning traffic) and O 3 to form NO 2 ,thusresulting in O 3 concentration depletion (Alier et al., 2011). After reaching its maximum at the end of the rush traffic hours at 8:00 AM, with the accumulation of solar irradiance hours, the photochemical catalytic reactions turn again in favor of O 3 formation, which reached its daily maximum between 14:00 and 16:00 h. This phenomenon can be observed in detail in Fig. 6a where the maximum of the O 3 hourly profile at 14:00 h coincides with the bottom of the valley in the humpback NO 2 profile. Fig. 6b and c provide additional information about the behavior of the O 3 -NO 2 interconnected patterns, e.g. the maximum of O 3 concentrations at 14:00 h occurred for all days, including working days and weekends (Fig. 6b), whereas the NO 2 humpback profile was more pronounced in working days, in agreement with the increased traffic and industrial activities (Fig. 6c). It is worth mentioning here that the simultaneous MCR-ALS analysis of both parameters, O 3 and NO 2 , have been also attempted. However, the results did not show any other interaction between them than the one already explained above for the independent analysis of both. 4.2. Comparison of prediction errors using the different CALIOPE model resolutions Comparison of the experimental data with the predicted data by CALIOPE at different model resolutions is performed using the RMSE daily Fig. 4. MCR-ALS results obtained in the analysis of O 3 experimental data (black) and of the three CALIOPE EU12 (green), IP4 (red) and CAT1 (blue) predicted data. C1, C2 and C3 daily (a, d and g), hourly (b, e and h) and station (c, f, i) profiles of the three MCR-ALS components. MCR-ALS explained variances are given in Table 1. S. Platikanov, M. Terrado, M.T. Pay et al. Science of the Total Environment 806 (2022) 150923 8 and RMSE hourly profiles calculated by Eqs. (6) and (7).Fig. 7a shows the RMSE daily O 3 profiles obtained for the CALIOPE predicted data. Results present some bias, with curved RMSE profiles in the case of CAT1 and IP4 models. The EU12 predictions gave a RMSE daily profile without a significant seasonal dependence. In summer, when the highest O 3 concentrations are reported, the CAT1 and IP4 prediction errors also increased. In winter, the CAT1 and IP4 models yield lower RMSE daily values than the EU12 model, meaning both models are working better in the lower range of O 3 concentrations. This can be related to the different spatial resolution adjustments used in the CALIOPE system. Fig. 7bshowsRMSE daily NO 2 profiles obtained for the CALIOPE predicted data. Slight increases of the RMSE daily rofile can be observed for the winter season, coinciding with the higher NO 2 concentrations reported during this season. In this case, the discrepancies among the CALIOPE predicted data and the experimental data were lower and more similar. Fig. 7c shows the RMSE hourly (Eq. (7))O 3 profiles obtained for the CALIOPE predicted data. Larger RMSE hourly peak shaped O 3 profiles were obtained for the EU12, IP4 and CAT1 data. The RMSE hourly O 3 profile for the EU12 data had a completely different shape to those obtained for the IP4 and CAT1 data. This suggests Fig. 5. MCR-ALS results obtained in the analysis of NO 2 experimental data (black) and of the three CALIOPE EU12 (green), IP4 (ref) and CAT1 (blue) predicted data. C1, C2 and C3 daily (a, d and g) hourly (b, e and h) and station (c, f, i) profiles of the three MCR-ALS components. MCR-ALS explained variances are given in Table 1. S. Platikanov, M. Terrado, M.T. Pay et al. Science of the Total Environment 806 (2022) 150923 9