{'fraction_non_alphanumeric': 0.03953272895736061, 'fraction_numerical': 0.028577657467797904, 'mean_word_length': 4.48301623667696, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Coupled multifractal methods to reveal changes in nitrogen dioxide and tropospheric ozone concentrations during the COVID-19 lockdown

Due to COVID-19 pandemic, the lockdown effects on air pollution level are undeniable. Several studies around the world have detected an uneven behaviour in tropospheric ozone (O3) concentrations. In this work, Seville (Spain) is used as example of faced to traffic place in which the nitrogen dioxide (NO2) is drastically reduced (41%) while O3 has no significant changes. In order to evaluate the existence of differences in O3 behaviour that is not detected by statistical procedures, a multifractal approach was used to assess the coupled scale relationship between NO2 and O3 during the 2020 lockdown against a period reference (2017–2019). For this purpose, the two main coupled multifractal method were employed: multifractal detrended cross-correlation and joint multifractal analysis. While cross-correlation analysis did not detect differences between the cross-correlated fluctuations of NO2 and O3 in the periods analysed, the joint multifractal analysis, based on the partition function and the method of moments, found a loss of variability in O3 during the lockdown. This leads to a loss of multifractal characteristic of O3 time series. The drastically reduction of primary pollutants during the lockdown might be the responsible of the tendency to monofractality in O3 time series. These differences were found for a wide temporal extent ranging from 80 min to ~28 days.


Introduction
In order to analyse the stochastics properties of air pollutant concentration time series, (Lee, 2002) explored the long term memory, selfsimilar and invariant patterns at time scales for these time series using a multifractal approach. In this pioneer study, multifractal characteristics were found in several time series, such as tropospheric ozone (O 3 ) and nitrogen dioxide (NO 2 ). The existence of multifractal features in pollution time series means that, on the one hand, data scale differently in different regions, which can be explained from completely different stochastics processes such as air movement and pollution transportation (Lee, 2002). On the other, multifractality is related to the random dilution process, that can be described through a random multiplicative process (Lee et al., 2003). These pioneers' works validates multifractal approach for enhancing insight into the complex nature of pollution time series.
O 3 dynamics have attracted the attention of researches due to the multiple processes, precursors and factors involved in its formation, transport and destruction. Indeed, O 3 could be the result of a deterministic system influenced by chaotic interdependent variables, which could be also capture by means multifractal analysis (Lee et al., 2006). Thus, O 3 time series have been studied from different multifractal methods, by considering different time periods and resolution. Several works focused on the study of O 3 as a single variable, while others used joint techniques for coupling with other factor and precursors. Table 1 gives a brief description of the works concerning O 3 pollution time series from a multifractal approach, encompassing pollutants involved, time resolution, period and multifractal method.
Regarding O 3 as a single variable, it is well known that this pollutant exhibits a temporal scaling strongly conditioned by the season of the year, i.e. climatological parameters influence O 3 formation (Jiménez-Hornero et al., 2010b;Plocoste et al., 2019). These results are also supported in the contextual multifractal analysis of visibility graphs (Carmona-Cabezas et al., 2019). Indeed, the degree of multifractality is conditioned by the year studied and prevail over the measurement location within the city (Diosdado et al., 2013). According to (Pavón-Domínguez et al., 2013), chemical dynamics highlight the differences in multifractal features among different monitoring sites when weather conditions for O 3 formation are not optimum (cold seasons), and also where anthropogenic primary emissions are not significant (rural places). In addition, there is a weaker degree of multifractality when pollution levels are low (Plocoste et al., 2017). O 3 , among other pollutant time series, was also analysed as a single variable (Dong et al., 2017) by means the Multifractal detrended fluctuation analysis (MF-DFA) (Kantelhardt et al., 2002).
There are several multifractal methods for jointly evaluate O 3 time series and its precursors or factors. The most employed coupled methods for studying tropospheric ozone behaviour have been the Joint Multifractal Analysis based on the Partition Function (MF-X-PF) (Meneveau et al., 1990) and Multifractal Detrended Cross-correlation Analysis (MF-DXA) (Zhou, 2008). Regarding MF-X-PF method, seasonal effects were also evident in ground-level relationships between O 3 and NO 2 (Jiménez-Hornero et al., 2010a). When the effect of air temperature is considered in the MF-X-PF, the higher temperature and NO 2 concentration, the greater degree of multifractality in O 3 time series (Pavón-Domínguez et al., 2015). With regard to MF-DXA approach, crosscorrelation patterns between NO 2 and O 3 exhibits a long-term crosscorrelations behaviour in rural (Xu et al., 2018) and traffic sites . According to , multifractal degree, shape and length of the tails in the cross-correlated multifractal spectra between both pollutants depends on the time resolution of the data and the season of the year.
Recently, new multifractal approaches for long-term cross-correlation ozone description along with its chemical precursors have been tested. The Coupled Detrended Fluctuation Analysis (CDFA) was employed for assessing the long term coupling correlation between O 3 and two main precursors: NO 2 and NMHC  and for calculating the proportion of coupling sensitivity degree of each precursor (Liu et al., 2021). CDFA method is an extension of the MF-DXA developed by (Hedayatifar et al., 2011) for multifractal analysis when more than two time series are correlated to each other. Indeed, a CDFA for six pollutant time series (CO,NO 2 , O 3 , PM 2.5 , PM 10 , SO 2 ) was developed by (Wang, 2019). On the other hand, (Stan et al., 2020) obtained multifractal properties and persistent long-range cross-correlations among O 3 and several precursors such as NO x , CH 4 , NMHC and THC using the Multifractal Multiscale Detrended Cross-correlation Analysis (MM-DXA). The MM-DXA (Shi et al., 2014) allows to visualize the cross-correlation by a Hurst surface, where the points on the surface represent the generalized dependence H xy (q, s). A further step jointly considering 6 pollutant time series consist of the representation of the Hurst surface in a compact format called multiscale multivariate multifractal detrended fluctuation analysis method (MMV-MFDFA) (Fan et al., 2019).
In literature, several epidemiological studies have already demonstrated the harmful effect of O 3 on human health (Chen et al., 2007;Gryparis et al., 2004;Lippmann, 1989;Lippmann, 1991;Nuvolone et al., 2018), to cite a few. It is a respiratory irritant that may cause decreased lung function and has been linked to other important respiratory health effects. This is the reason why during the COVID-19 pandemic, several studies around the world have focused on studying the behaviour of O 3 following the lockdown or the curfew. Due to the complexity of O 3 formation, its behaviour varies from region to region. Indeed, many studies have found an expected increase in O 3 concentration (Anil and Alagha, 2021;Briz-Redón et al., 2021;Donzelli et al., 2021;Rojas et al., 2021;Salazar et al., 2020;Venter et al., 2020;Zambrano-Monserrate and Ruano, 2020;Zhang et al., 2020) while others showed a decrease (Adams, 2020;Donzelli et al., 2021;Kumari and Toshniwal, 2020). In some cities, there is no significant change in O 3 concentration before and during the lockdown (Aljahdali et al., 2021;Donzelli et al., 2020;Varotsos et al., 2021). These studies are mainly based on statistical analysis of ground-based measurements. However, the statistical analysis of time series is only based on a single scale analysis and might not necessarily reflect the features of the time series over other scales (Zeleke and Si, 2006). This approach may not be appropriate for nonlinear and nonstationary time series as O 3 (Chelani, 2010).
To authors knowledge, only two studies have applied multifractal method for characterizing the lockdown effects on air pollution due to COVID-19 pandemic (Li, 2021;Sipra et al., 2020). The first one studied PM 2.5 behaviour while the second focused on air quality index time series before and during the lockdown. The only one work concerning multifractality in O 3 during a road blockage is found during the Hong Kong protest in 2014 . In the three aforementioned works, MF-DFA is employed for assessing air pollutants time series as single variables. Therefore, in this work is conducted a first study on the joint scaling behaviour of O 3 and NO 2 at ground level during a strong reduction of traffic emissions, i.e. during the lockdown due to COVID-19 pandemic. For that purpose, the two main coupled multifractal algorithms (MF-X-PF and MF-DFA) are applied.
In order to carry out this study, the remainder of the paper is organized as follows. Section 2 presents the study area and the data. Section 3 describes the theoretical framework used. Section 4 comments on the results obtained and provides a discussion. Finally, a conclusion and an outlook for future studies are given in Section 5.

Site description and data collection
This analysis is developed from ground-level air pollution data recording in a monitoring station faced to traffic emissions in Seville (Spain). Seville (37 ∘ 24 ′ N, 5 ∘ 60 ′ W and 20 masl) is the Andalusia's capital and the largest city in the southern region of Iberian Peninsula. Its metropolitan area has a population around 1,500,000 inhabitants (Feria-Toribio, 2015) and is comprised by several localities surrounding the main city of Seville. Mediterranean climate governs the weather in Seville, which is characterized by a high solar radiation (~3000 h year − 1 ) and an irregular rainfall, mainly in spring and autumn. Average temperatures range from 6.3 ∘ C in winter to 33.8 ∘ C in summer. A detailed description of O 3 behaviour and the study area can be found in (Massagué et al., 2021).
According to the ambient air quality report in 2019 elaborated by the Andalusian regional government (CAGPDS, 2019), days with poor air quality in Seville are mainly due to high levels of PM 10 and O 3 . Specifically, around 20 days in 2019 had ozone levels higher than 180 μgm − 3 for 1 h average concentration, which is known as population information threshold in the EU Directive (Directive 2008/50/EC, 2008. The Environmental Department of the Regional Government of Andalusia manages the air quality network in the region. In Seville's metropolitan area are located 10 monitoring stations that record a total of 5 pollutants: SO 2 , CO, NO x , O 3 and PM 10 . As the aim of this work is to assess the impact of lockdown on urban mobility emissions, Torneo monitoring station is selected for collecting the data. Torneo monitoring station is located next to a road with heavy traffic and close both exit towards roads and inner-city ring road. Thus, it is highly influenced by traffic emission (Adame et al., 2008).
The declaration of the state of alarm due to COVID-19 pandemic in Spain came into effect from March 16, 2020 to May 3, 2020. This period was associated to a very strict home confinement, drastic reduction of job activity and full closure of educational centers from kindergarten to university. During the strict lockdown, only basic consumer goods and medicines could be purchased and economic activity was limited to essential activities. The lockdown period was comprised by 7 weeks, in which 10-min NO 2 and O 3 concentrations were collected in Torneo monitoring station, this resulting in 7056 10-min data (CMAOT, 2020). For comparison purposes, the same dates in 2017, 2018 and 2019 are established as reference period. Missing values represent less than 4% of the data. They were replaced by the average value of the day after and day before at the same time.

Seasonal and Trend decomposition (STL)
NO 2 and O 3 time series are strongly linked to daily cycles (Adame et al., 2008;Massagué et al., 2021), these resulting in the presence of a seasonal component in time series. The presence of seasonality produce distortions impeding the understanding of their long-range correlations (Li et al., 2015). In order to analyse the inner fluctuations in pollution time series, the seasonal and trend components are removed, maintaining the remainder (residual) component. The study of the residuals produces more reliable results (Laib et al., 2018) since the analysis is focused on the inner dynamics of NO 2 and O 3 time series. Hence, a STL decomposition of the signal (Cleveland et al., 1990) is carried out for removing the deterministic component (seasonality and trend) and to extract the stochastic component (remainder). Seasonal-Trend de composition is a filtering method based on Loess smoother for decomposing time series (Y t ) into three components: seasonality (S t ), trend (T t ), and remainder (R t ). The sum of S t , T t , and R t is the original signal (Y t ): The STL parameters employed in this work are defined from (Cleveland et al., 1990) as follows: As the seasonal cycle is one day and the amount of 10-min data in a day is 144, the number of observations in each cycle of the seasonal component is n p = 144. Due to the absence of outlier in time series, the number of robustness iterations of the outer loop is n o = 0, and the number of passes through the inner loop is selected as n i = 2. The smoothing parameter for the low-pass filter, n l , is an odd integer greater or equal to n p , thus n l = 145. The smoothing parameter for seasonal component, n s , is an odd integer greater than 7. It was selected experimentally by developing STL with n s = 17, 27, 37 and 47. Decomposition results were quite similar with all the n s values tested, so thus the minimum one was selected (n s = 17). Finally, the smoothing parameter for the trend component is calculated from n p and n s as the lower odd satisfying: 1.5n p /(1 − 1.5n s − 1 ), so n t = 237. After applying the STL decomposition on the original NO 2 and O 3 time series, deterministic components were removed and the remainder components were used as input signals in subsequent multifractal analysis.

Multifractal detrended cross-correlation analysis (MF-DXA)
MF-DXA (Zhou, 2008) is based in the detection of long-range crosscorrelations and multifractal scaling properties in nonstationary time series for different qth order moments. Therefore, this method is capable of evaluating how two time series are cross-correlated in several time scales (Shadkhoo and Jafari, 2009). It is a generalization of the Multifractal Detrended Fluctuation Analysis (MF-DFA) for a single variable developed by (Kantelhardt et al., 2002) and the Detrended Crosscorrelation Analysis (DXA) devised by (Podobnik and Stanley, 2008). Here, the main interest is focused on the exploration of the multifractal cross-correlation properties between the remainder of NO 2 and O 3 concentration time series in each period analysed.
Initially, consider that x k and y k are time series of length N and of compact support, i.e., non-zero values or an insignificant fraction of zero values are present. Then, MF-DXA is conducted by following five steps: • Step 1: Determine the integrated time series X(i) and Y(i).
Commonly named"profile", the integrated time series are computed by subtracting the mean x and y to the accumulative sum of each signal respectively: (2) i taking values from 1 to the total length of the series (N). • Step 2: Divide the integrated time series X(i) and Y(i) into N s = int (N/s) non-overlapping segments of equal lengths s. Since the series length N could not be a multiple of the considered time scale s, the procedure is conducted twice. One, from the initial time series to the end, and secondly from the end to the beginning. Thereby, 2N s segments are obtained for each s value. In both cases, a short part at the end (or at the beginning) of the profiles may remain. This step exhibits an advantage of this method against the joint multifractal analysis based on the partition function (MF-X-PF), in which the total length of the time series must be a power of time scale.
• Step 3: Calculate the covariance by detrending the local trend for each of the 2N s segments. The local trend is evaluated in each segment v i (s) by means a polynomial fit x v (i) and y v (i) obtained from the least-squares method, respectively. In this work linear fits have been used, even if linear, quadratic or higher order polynomials can be taken in the fitting procedure. Later, considering the aforementioned fits, trends are removed by subtracting x v (i) from x(i) and y v (i) from y(i) for each segment v. Then, the detrended covariance of the residual, F 2 depending on s, is determined for each segment v as follows: Step 4: Determine qth order detrended covariance. From this step, the multifractality of the local covariances are evaluated through magnifying both higher and lower local fluctuations by employing different q order moments. Local variances are assembled in the denominated q-order fluctuation function for two variables, F xy (q, s), by averaging over all the fluctuations of the segments, defined as: when q∕ = 0. For q= 0 the fluctuation function diverges and needs to be estimated by applying a logarithmic averaging procedure: In this method we are interested in the dependence of the generalized fluctuation function, F xy (q, s), versus different time scales. Hence, repeating steps 2-4 for several time scales s, it is apparent that F xy (q, s) will increase with increasing s.
Whether x k and y k time series are long term power-law crosscorrelated, the fluctuation function F xy (q, s) scales with s, according to a power law: where h xy (q) is the generalized cross-correlation exponent. In the case that h xy (q) is independent of q, the generalized cross-correlation exponent is a horizontal straight line. By contrast, multifractality exists when h xy (q) is a decreasing function depending on q. The greater dependence of h xy (q) of q, the greater multifractality. The strength of the multifractality is measured through Δh xy (q). Specifically, when q = 2, the generalized cross-correlation exponent exhibits the same properties as a univariate Hurst exponent. h xy (2) is named the second cross-correlation exponent and provides useful information about the long-term powerlaw cross-correlation between both time series. h xy (2) > 0.5 indicates persistency, that is, large (or small) fluctuations in one variable are more likely followed by large (or small) values in another variable. Conversely, h xy (2) < 0.5 indicates anti-persistency, meaning that large (or small) fluctuations in one variable are more likely followed by small (or large) values in another variable. When h xy (2) = 0.5, both time series are no long term cross-correlated.
In spite of the above method is equivalent to the single MF-DFA when x k = y k , in this study the MF-DFA is directly conducted through the method proposed by (Kantelhardt et al., 2002). For more details of this method we refer the reader to (Kantelhardt et al., 2002;Peng et al., 1994). A brief description of the 5 steps that comprised the MF-DFA are as follows: Initially, consider that x k is a time series of length N and of compact support.
• Step 1: The integrated time series X(i) is computed by subtracting the mean (x) to the accumulative sum of the time series. • Step 2: Divide the integrated time series X(i) into N s = int (N/s) nonoverlapping segments of length s.
The integrated time series X(i) is divided into N s = int (N/s) nonoverlapping and identical segments of length s. In order to consider the whole dataset, the procedure is carried out from the beginning to the end and from the end to the beginning of the integrated time series. Hence, 2N s segments are obtained for each s value.

•
Step 3: Local trends are calculated by least-square fits for each 2N s segments. Then, the variance of the residual, F 2 (s, v), is calculated: for • Step 4: Assemble the variances in the denominated generalized fluctuation function F q (s), by averaging over all the fluctuations of the segments and considering a set of qth order, according to: and for q = 0: • Step 5: Determine the scaling behaviour of the fluctuation functions from the log-log plots F q (s) versus s for each q value. If a power-law can be defined for a certain scaling range, h(q) can be obtained from: Whether small and large fluctuations scale differently, h(q) depends on q and the time series is considered as a multifractal signal. By contrast, when the contribution of both large and small fluctuations are similar, h(q) is independent of q and the time series is monofractal. Additionally, the value of h(q) for q = 2 gives useful information about the time series. For non-stationary signals, h(2) is greater than 1 and the Hurst exponent is calculated as H = h(2) − 1. For stationary signals, h(2) is identical to the Hurst exponent H. The Hurst exponent can take values greater or smaller than 0.5 when the series exhibits persistency or antipersistency, respectively. Time series is uncorrelated when H is equal to 0.5.

Joint multifractal analysis based on the partition function (MF-X-PF)
Multifractal spectrum f(α) formalism is related to the set of generalized dimensions D(q) (Grassberger, 1983;Hentschel and Procaccia, 1983). The determination of multifractal spectrum f(α) based on the strange attractor formalism (Halsey et al., 1987) was extended by (Meneveau et al., 1990) for describing the degree of correlation among two or more multivariate measures coexisting in the same geometric support. In this work, we employed the method for the direct estimation of the multifractal spectrum ) and the methods of moments (Evertsz and Mandelbrot, 1992) for computing the joint multifractal spectra between NO 2 and O 3 for the study period.
By considering that x k and y k are time series of length N coexisting in the same domain, the joint multifractal analysis based on partition function for two multifractal measures (Meneveau et al., 1990) can be summarized into five steps: • Step 1: Divide the domain. For time series, the domain is a temporal line that is successively divided into n identical and non-overlapping intervals of increasing length ranging from δ ini to δ final time resolutions. The length of the intervals must be power of 2, this resulting in a shortcoming of this method when the total length N of the time series is not a power of 2. The values [x k ] ini and [y k ] ini in any initial interval j are set to be equal to the sample measurement and thus, δ ini is defined as the minimum time resolution that allows to assign a non-zero measure for each initial interval. Increasing the time resolution (δ), the sum of the values [x k ] ini and [y k ] ini which are contained in each interval of resolution (δ) represent the value of this interval. Finally, the total sum of the time series of each variable is assigned for the interval δ final . • Step 2: Calculate de probability mass functions. The probability mass function c i , of each interval δ, are defined in terms of a mass density for each time series: where [x k ] i and [y k ] i are the sum of the values of each time series which comprise the interval i for a specific time resolution δ. Thus, for each time resolution, probability mass functions c i [x k ](δ) and c i [y k ](δ) represent the measure of each interval i respect to the total sum of the values of each time series. • Step 3: Obtain the joint partition function. The joint study of the probability mass functions through different time scales is conducted by the method of moments (Evertsz and Mandelbrot, 1992). Joint partition function allows to evaluate the existence of scaling properties in the joint distribution of the time series. It is obtained by: where n is the amount of intervals for a specific time resolution δ; and q x and q y are the q-moments for variables x k and y k , respectively. As noted, joint partition function exclusively depends on the time resolution δ and the values of the two q-moments. The higher positive value of q x (q y ), the greater emphasis on high values of x k (y k ). Conversely, the higher negative value of q x (q y ) the greater emphasis on low values of variable x k (y k ). For (multi-)fractal signals, the joint partition function has the following scaling property: • Step 4: Estimate the mass exponent function τ(q x , q y ). The mass exponent function, which depends only on the q-moments, is estimated from the slope of the linear fits of the logâ€"log plot of the joint partition function χ(q x , q y , δ) versus the time resolution δ. From this logâ€"log plot it is determine the possible existence of multifractality between both time series, as well as the range of scales for which multifractality exists. • Step 5: Obtain the joint multifractal spectrum. As multifractal behaviour can be represented by the multifractal spectrum (Halsey et al., 1987), the joint multifractal behaviour for two variables can be analysed by means the joint multifractal spectrum, f(α x , α y ) versus α x and α y . f(α x , α y ) can be considered as the fractal dimension of the set of intervals that corresponds to a combination of the singularities α x and α y (Meneveau et al., 1990). The singularity exponents α x and α y , also named Lipschitzâ€"HÃ ¶lder exponents, quantify the singularity strength of the measures. The singularity exponents are determined by the Legendre transformation of the mass exponent for each qmoment (Evertsz and Mandelbrot, 1992): Let N(α x , α y , δ) be the number of intervals of time resolution δ where a given combination of α x and α y values are found, and define f(α x , α y ) from the scaling relationship: Hence, f(α x , α y ) is calculated from Meneveau et al., 1990): Multifractal spectrum for time series reaches a maximum value for f (α x , α y ) equal to 1, which corresponds to the capacity dimension of the geometric support for time series (one-dimensional). Whether high values of q x are considered, the single multifractal spectrum f(α y ) represents the multifractal behaviour of time series y k when high values of x k are emphasized. Conversely, when low values of q x are considered, the single multifractal spectrum f(α y ) represents the multifractal behaviour of time series y k when low values of x k are emphasized. When q x is taken as zero, the joint multifractal spectrum is identical to the single multifractal spectrum of f(α y ), i.e. the multifractal behaviour of y k when x k is not considered. For monofractal time series, δ is the same for all the intervals of identical time resolution and the multifractal spectrum is a single point (Kravchenko et al., 1999). Table 2 summarizes the descriptive statistics for NO 2 and O 3 time series during the lockdown and the reference period. As seen, whereas changes in NO 2 descriptive statistics are evident, O 3 remains equal in 2020 respect to the reference period. During the lockdown, the mean [NO 2 ] is notably reduced, as well as the maximum value. Indeed, compared to the reference period, a reduction of 41% is observed for [NO 2 ]. This behaviour has been observed in many countries around the world (Kumari and Toshniwal, 2020;Venter et al., 2020). Skewness is greater, which evinces that the frequency of low data in 2020 is greater than in the reference period. By contrast, the differences in [O 3 ] are slight and not conclusive. Compared to [NO 2 ], [O 3 ] during the lockdown remain equal to the reference period. Overall, the maximum [O 3 ] in 2020 is slightly higher and the dispersion lower. Moreover, despite the reduction of [NO 2 ] in 2020, the mean value of [O 3 ] is not the highest of the studied years. It is important to recall that O 3 formation is a complex nonlinear process involving the interaction of sunlight with nitrogen oxides (NO x = NO 2 + NO) and Volatile organic compounds (VOC) (Duan et al., 2008;Krupa and Manning, 1988). Usually, for high VOC/NO x ratios (r> 12), the chemistry of the air masses is NO x -limited, and NO x control is the more effective way of reducing ozone levels. These situations are specific of suburban and rural areas. In urban centers where VOC/NO x ratios is traditionally low (r≤ 6), these systems are VOCcontrolled, i.e. reducing NO x at constant VOC leads to an increase in ozone concentrations (Da Silva et al., 2018;Dantas et al., 2020;Siciliano et al., 2020). A previous study made by (Adame et al., 2008) already showed that Seville is a city under VOC-limited regime as VOC/NO x ratios is equal to 2.4. This may be one of the reasons why [O 3 ] is not sensitive to the reduction of [NO 2 ]. In addition, NO 2 which has an atmospheric lifetime of about a day are clearly discernible locally while O 3 which has a lifetime of several weeks are affected by long-distance transport associated with specific weather patterns (Venter et al., 2020). Fig. 1 exhibits the typical diurnal pattern for these air pollutants. In the reference period, two evident peaks in the hourly average NO 2 concentration are found at 8.00 and 21.00 h linked to heavy traffic emissions due to rush hour ( Fig. 1(a)). However, during the lockdown, the mean NO 2 values are much lower than those of the reference period. In addition, the evening peak is notably mitigated due to very strict traffic movement. As consequence of the drastic traffic reduction, a significant decrease in pollutants, especially NO 2 , was detected in the  most populated cities in Spain (Cárcel-Carrasco et al., 2021). Whereas a significant NO 2 reduction was achieved in most Spanish cities, other pollutants such as CO, PM 10 and SO 2 were only reduced in some locations, as well as O 3 concentration, that only increased in some cities (Briz-Redón et al., 2021). Regarding O 3 in Seville (see Fig. 1(b)), the daily cycle follows the expected spring behaviour in Torneo, which is in agreement with (Adame et al., 2008). The minimum value is found at 8.00 h, when traffic emissions in the early morning reach the maximum. O 3 concentration increase from early morning to 16.00 h due to optimum mechanisms for photochemical formation. In addition, recirculation of air masses promotes horizontal and vertical transport of O 3 from other places to Seville, which is added to the photochemical O 3 formation in situ (Adame et al., 2008). From afternoon, there is a decreasing trend in O 3 concentration until it takes stable values during the night. According to Fig. 1(b), lockdown period resembles as a common year with very similar mean hourly values than reference period. Even if no increase in O 3 concentration during the lockdown is found, O 3 destruction from 21.00 h differs from the reference period. Indeed, the reduction of nitrogen oxides emissions in 2020 lead to reduce the local titration of O 3 , i.e. the reaction of NO with O 3 (Venter et al., 2020). Consequently, the chemical destruction of O 3 is more limited in late evening.

Meteorological parameter
Beside traffic and industry emissions during the lockdown, meteorological conditions migth also affect air pollutant concentration in cities (Anil and Alagha, 2021;Ginzburg et al., 2020;Sulaymon et al., 2021). In several places in the world, there has been no statistically significant reduction in O 3 concentration during the shutdown period (Adams, 2020;Aljahdali et al., 2021;Varotsos et al., 2021), this being attributed to a prevalence of the meteorological conditions (Varotsos et al., 2021). Meteorological conditions may even mask the consequences of the reduction in human activity (Sipra et al., 2020). Meteorological data during the study periods in Seville were provided by the (CAGPDS, 2020) of the regional government of Andalusia. The agrometeorological monitoring station of La Rinconada (located at 8 km from Seville) provided daily measurements of the main meteorological variables involved in O 3 behaviour, such as temperature, air humidity, solar radiation and rain fall (see Table 3). Wind speed and direction are depicted in Fig. 2.
According to Table 3, meteorology during the lockdown period in Seville is comparable with the reference period. Mean temperature and air humidity are in agreement with the reference period, but mean solar radiation is lower than in the reference period. In addition, 2020 was a rainy year in comparison to 2019 and 2017, but rather similar to 2018. Fig. 2 depicts the wind roses for the study periods, which are involved in air pollutant movement from surroundings areas to the city. As seen, wind charts indicate a SW-NE regime for the lockdown (2020) and reference periods in 2017 and 2019, which matches with the typical spring regime in Seville. In spring, a dual wind regime is usual in Seville across the Guadalquivir Valley, with nocturnal air masses coming from NE and diurnal Atlantic air masses from SW. In 2018 wind regime is more similar to a winter scenario (Adame et al., 2008).

Multifractal analysis
Before performing the multifractal analysis, STL decomposition described in subsection 3.1 is conducted in order to remove seasonal components in time series. Thus, multifractal algorithms are applied for the remainder components of each air pollutant. The analysis of the remainder components provides more reliable results and information on the inner dynamics of the time series (Laib et al., 2018). As example, Fig. 3 shows STL decomposition of NO 2 and O 3 for the lockdown period (2020). From the top to the bottom, panels depict the 10-min original data, seasonal, trend and remainder components.

Multifractal detrended cross-correlation analysis (MF-DXA)
After applying the MF-DXA described in section 3.2 to the remainder of both pollutants for the lockdown and the reference periods, the loglog plot of the cross-correlated fluctuation function F NO2, O3 (q, s) is depicted in Fig. 4. The procedure was conducted by using a set of q moments ranging from − 3.5 to 3.5 at 0.5 increments and 80 timescales (s) logarithmically spaced from 10 to N. For determining the scaling behaviour among these variables we focus on the linear regions. As seen in all the study years, a linear region for small scales can be detected for scales lower than a crossover (s x ). Analogously, a second linear region is found for larger scales from the crossover (s x ) to higher scales. Lower and upper limits for linear fits has been made according to (Kantelhardt, 2012). Thus, scales lower than s < 10 have been discarded because they produce systematic deviations from the scale behaviour in Eq. 8. On the other hand, scales higher than s > N/4 are excluded because fluctuation function becomes statistically unreliable since the smaller amount of intervals N s in the averaging procedure for determining the qth order detrended covariance. For comparison purposes, the same cut-off have been selected for the four analysed periods, this resulting in two regimens with different scaling exponents. Whereas small-scale regime is established from s min =10 (~2 h) to s x =142 (~1 day), long-scale regime is defined from s x =142 (~1 day) to s max =N/4= (~12 days). Fig. 5 shows the q-dependence of the generalized cross-correlation exponent h xy (q) for x = NO 2 and y = O 3 determined by linear fits for (a) small-scale and (b) large-scale regimes. h NO 2 , O 3 (q) are decreasing functions with the increase of q, which confirms the multifractal behaviour in the cross-correlation between NO 2 and O 3 in each study period. (Xu et al., 2018) also confirmed the multifractal nature of the power-law cross-correlations between NO 2 and O 3 in three traffic sites in Hong Kong. The goodness of the linear fitting is measured from the coefficient of determination (R 2 ), whose minimum value for the most unfavourable q moment is R 2 > 0.95 for small-scales and R 2 > 0.91 for large-scales. Additionally, slope estimation errors are also depicted as vertical error bars in Fig. 5. As seen, greater errors are found for the estimation of h NO 2 , O 3 (q) for negative q-values, but these being lower than 0.045 for the most unfavourable case. Differences in multifractal results for small and large scales confirms that the choice of the crossover was suitable.
Focusing on small-scales ( Fig. 5(a)), h NO 2 , O 3 (q) decrease stronger for negative moments than for positive ones. In fact, functions tend to resemble flatted for q positive values, which highlight an insensitive multifractal structure of NO 2 and O 3 to local cross-correlation fluctuations with a large magnitude (Ihlen, 2012). Second cross-correlation exponents, h NO 2 , O 3 (2), are quite similar, ranging between 0.565 and 0.591, the h NO 2 , O 3 (2) for 2020 (lockdown period) being within this interval. Thus, all of them are greater than 0.5, but close to the uncorrelated situation h xy (2)=0.5, indicating that there are a slightly persistent cross-correlations between NO 2 and O 3 for small-scales. On the other hand, multifractal cross-correlation degree, which is measured as Δh NO 2 , O 3 (q), indicates that multifractality strength in 2020 (0.283) is greater Table 3 Daily average values of several meteorological variables during the lockdown period (2020) and the reference period (2018, 2018 and 2017)    Beside similar multifractal parameters, generalized cross-correlation exponent functions overlapped, and thus, no differences between the lockdown and the reference period is found in multifractal cross-correlations between NO 2 and O 3 for small scales (lower than 1 day). As regards to large-scales ( Fig. 5  (b)), h NO 2 , O 3 (q) are decreasing functions for all q moments considered, which confirm the existence the multifractality in the cross-correlations between NO 2 and O 3 . Contrary to small-scales, second cross-correlation exponents, h NO 2 , O 3 (2), are lower than 0.5, which indicates an antipersistent cross-correlations between NO 2 and O 3 . h NO 2 , O 3 (2) values range between 0.378 and 0.485, with years 2017 and 2018 close to the uncorrelated case. Whereas in 2020 and 2019 the strength of multifractality is lower than for small-scales, it is greater for 2018 and 2017 when timescales are greater than 1 day. Therefore, differences among lockdown and reference years are not conclusive, while MF-DXA results in 2020 and 2019 are quite similar. By comparing to previous works, persistency was found between NO 2 and O 3 after applying the MF-DXA at traffic sites in Hong Kong (Xu et al., 2018). However, time series in the current work are shorter than in the aforementioned one. Other differences in the dataset, such as time resolution, range of scales for linear fits and the existence of crossovers might could hinder the comparison of the results. Given that the evident NO 2 reduction during the lockdown (see Fig. 1 and Table 2) is not detected by MF-XDA, a simple MF-DFA analysis is then carried out with the purpose of finding possible differences between the periods analysed when the variables are studied independently. Multifractal parameters are the same that those employed in the MF-DXA with identical scales and q moments. Fig. 6 depicts the generalized Hurst exponent, h(q), for NO 2 and O 3 in the analysed periods. It should be noted that a crossover at 1 day is also found, and thus, the same timescales are established for linear fits. Goodness of the fitting is ensured by coefficient of determination greater  Fig. 6(a) and Fig. 6(b), respectively. The main difference in temporal scale regime is found in the value of the Hurst exponent, h(q = 2), which ranges from 1.098 to 1.155 for small scales and from 0.742 to 0.880 for large scales. This indicates that NO 2 remainder is close to a pink noise (H = 1) for scales lower than 1 day and exhibits persistency for those scales greater than 1 day. Surprisingly, intrinsic fluctuation across time for NO 2 are kept in 2020, given the similarities displayed with respect to the reference period. On the other hand, as same as before, O 3 fluctuation across time are also kept in 2020. According to results, O 3 for both small and large scales are slightly multifractals with a multifractal degree (Δh) in the range 0.155-0.211 and 0.051-0.283, according to Fig. 6(c) and Fig. 6(d), respectively. For small scales (lower than 1 day), O 3 remainder are non-stationary signal. For non-stationary signals the Hurst exponent is calculated as H = h(q = 2) − 1, thus H are in the range of 0.235-0.277, this evincing that O 3 exhibits anti-persistency for small scales. Finally, for large scales, Hurst exponents range from 0.737 to 0.892, which denote the persistency of O 3 for timescales longer than 1 day. Results obtained in Seville follows the tendency described in Hong Kong during the road blockage in 2014 for NO 2 , but differs for O 3 . According to these author, VOC-limited conditions in Hong Kong (Xing et al., 2011) may raise the risk of increasing urban ozone levels when local NO x emissions are reduced.

Joint multifractal analysis based on the partition function (MF-X-PF)
Before applying the method described in section 3.3 to NO 2 and O 3 remainder time series, one main issue concerning the MF-X-PF is highlighted. This method has an important shortcoming related to the length of the time series, which must be a power of 2. Since time series are comprised by 7056 data, we have selected the maximum power of 2 amount of values from the end part of the time series, resulting in 4096 data = 2 12 . The analysed period correspond from April 16 at 13:30 h to Mai 3 at 00:00 h for the lockdown and the reference years. The end part of the time series have been selected as the most representative period for evaluating the effects of the decrease in mobility and traffic, weeks after the lockdown started.
We have selected a time resolution from δ ini = 2 0 to δ max = 2 12 and a set of q-moments order ranging from q = − 3.5 to 3.5 at 0.5 increments. From the log-log plot of the joint partition function versus time resolution, linear regions are evident from time resolutions greater than the minimum δ ini = 20 = 10 min. For comparison purposes, the same time scaling range for linear fits was selected for the lockdown and the reference years. Thus, linear fits are established from δ ini = 2 3 (80 min) to δ max = 2 12 (~28 days) with a mean value of R 2 of 0.965, 0.984, 0.952 and 0.972 for 2020, 2019, 2018 and 2017, respectively. This broad temporal scaling range agrees to similar multifractal studies involving O 3 and NO 2 during a month with 10-min time resolution (Jiménez-Hornero et al., 2010a;Jiménez-Hornero et al., 2010b;Pavón-Domínguez et al., 2013), in which δ ini is in the range 2 2 -2 4 (20-160 min) and δ max is the total length of the time series (~1 month).
Joint multifractal spectra are convex surfaces that reach f(α NO 2 , α O 3 )= 1 when both q-moments are equal to zero, as expected for the onedimensional temporal data. Fig. 7, depicts a projection of these surfaces on the plane α NO 2 − α O 3 . It is evinced the joint multifractality between NO 2 and O 3 for each year, by means the relationships of a certain fractal dimension f(α NO 2 , α O 3 ) for each combination of the singularity exponents α NO 2 and α O 3 . All of them exhibit triangular shapes oriented from the top left region to the bottom right with no high values of singularity exponents α NO 2 and α O 3 . This lack of values in the top right region indicates that no combinations of low NO 2 and O 3 concentrations are found. By contrast, the orientation top leftbottom right evince the inversely proportional relationship of these variables (Clapp and Jenkin, 2001;Han et al., 2011;Plocoste et al., 2018;Pudasainee et al., 2006), in which, higher values of NO 2 are likely related to lower O 3 concentrations and vice-versa. In addition, both regions are associated to lower fractal dimensions f(α x , α y ), this indicating that these extreme events are infrequent.
Even so, several differences in the shape of multifractal spectra can be detected. Focusing on years 2019 and 2017 (see Fig. 7(b) and (d)), they exhibit triangular shapes with a broader range of singularity exponents of O 3 than NO 2 , similar to those described for the springsummer seasonality in (Jiménez-Hornero et al., 2010b). By contrast, for years 2020 and 2018 (see Fig. 7(a) and (c)) a broader range of singularity exponents for NO 2 is found. As seen, joint multifractal spectra for the reference period is more similar to the year 2018, but completely different to 2019 and 2017. Despite the greater similarity in their shapes, the main difference between 2018 and lockdown year is that a narrower range of singularity exponents for both NO 2 (Δα NO 2 =0.85) and O 3 (Δα O 3 =0.19) is found in 2020 (see Fig. 7(a)). This means that the lowest NO 2 (O 3 ) concentrations are related to a smaller variability of O 3 (NO 2 ) concentrations in 2020. In general, the triangular shape described above tends to be smaller with narrower range of singularity exponents for both NO 2 and O 3 in the lockdown period. This loss of multifractality is consequence of both a greater degree of association and a lower variability in both variables. Multifractal analysis is a very sensitive tool that can detect small variations in different periods. This is usual in joint multifractal literature, in which different joint multifractal shapes where detected in the study of crop yield and terrain slope during four growing seasons in the same plot (Kravchenko et al., 2000) and in the analysis of the relationship of solar radiation and PM 10 in different years (Plocoste and Pavón-Domínguez, 2020).
By considering NO 2 as independent variable, three scenarios depending on the NO 2 concentration are selected for describing the multifractal O 3 behaviour for each study period. These scenarios are single multifractal spectra f(α NO2 , α O3 ) vs α O3 and are obtained by selecting the corresponding values of q NO2 from the joint multifractal spectra. The first scenario consider high values of NO 2 by selecting q NO2 = +3, the second corresponds to the study of O 3 without considering NO 2 (q NO2 = 0) and the third one focuses on low NO 2 values by taking q NO2 = − 3. They are depicted by triangles with vertex up, squares and triangles with vertex down in Fig. 7, respectively.
According to Fig. 8(a), O 3 multifractal spectrum for high NO 2 concentrations (q NO 2 = +3) are very asymmetric, exhibiting long right tails with α max reaching values of 1.13 (2020), 1.39 (2019)  are infrequent. As seen in 2020, high NO 2 concentrations (q NO 2 = +3) influences O 3 behaviour reducing the maximum f(α) and narrowing the width of its right tail. A narrower spectrum indicates that O 3 variability during the lockdown period was lower than in the reference period when high NO 2 concentrations are considered.
When O 3 is studied as a single variable (see Fig. 8(b)), asymmetric spectra with long right tails are found, mainly for the reference year. However, the asymmetry is lower than in the case of q NO 2 = +3, with α max reaching values of 1.05 (2020), 1.32 (2019), 1.09 (2018) and 1.41 (2017). Therefore, multifractality in O 3 time series is mainly due to the greater variability in low O 3 concentrations (longer right tails), while there is an accumulation of points in the left tails of the spectra which brings the variability of high concentrations closer to the monofractal case. This also suits with the high NO 2 concentration case (q NO 2 = +3).
The asymmetry with respect α is usually found in air pollutant time series, indicating that the low and high concentration values are not distributed evenly (Lee, 2002). A greater multifractal heterogeneity in low ozone values has also been describe in cities, mainly for spring and summer months (Jiménez-Hornero et al., 2010b). Thus, multifractal analysis adds suitable insights that complement statistical. On the one hand, statistical analysis showed that the greatest variability in ozone concentrations occurred in 2020, and on the other, the values of standard deviation of the study years are very similar (see Table 2). However, multifractal analysis highlights that there is a loss of variability in O 3 concentration during the lockdwon with respect to the reference period across the different time scales studied. Indeed, if O 3 is studied as a single variable (q NO 2 = 0), the lockdown period tend to show a monofractal behaviour. According to (Pavón-Domínguez et al., 2013) rural stations show less multifractality than those close to urban areas, due to the absence of anthropogenic emissions in situ. According to this, the drastically reduction of primary pollutants during the lockdown might be the responsible of the loss of O 3 multifractality in Torneo station.
As seen in Fig. 8(c) multifractal behaviour of ozone is extremely sensitive for temporal intervals of low NO 2 (q NO2 = − 3). The main P. Pavón-Domínguez and T. Plocoste similarity is that the spectra are left-skewed, indicating that when NO 2 concentrations are low, O 3 concentrations are usually high. However, the values of f(α) are notably lower than in the cases of q = + 3 and q =0, which is related to very infrequent episodes. As observed in Fig. 8(c), maximum values of f(α) are less than 0.3 for the years 2020, 2018 and 2017. Looking for analogies among studied periods, it might be observed that spectra in years 2018 and 2020 completely lack a right tail, while this do not occur for 2017 and 2019. Years 2018 and 2020 show spectra with overlapping points, which confirms their tendency to be monofractal. When low NO 2 concentrations are considered (q NO2 = − 3), multifractal spectra show notably differences among the studied years and no significant differences can be extracted from the lockdown respect to the reference period.

Conclusions
Last decades, multifractal analysis has been carried out for studying air pollutant time series. In this work, novel contributions in the scaling behaviour of NO 2 and O 3 time series are shown. Firstly, this study employs the two main coupled multifractal approaches for assessing their efficiency to detect no significant changes in air pollutant time series. On the one hand, multifractal analysis based on cross-correlations (MF-DXA), which emphasizes the long-range relationships between the fluctuations of two variables and, on the other, the traditional joint multifractal analysis (MF-X-PF), based on the relationships between concentration data of each variable. Secondly, both scaling analysis techniques have been employed in order to detect possible changes in O 3 behaviour during the lockdwon due to the COVID-19 pandemic.
10-min time series of NO 2 and O 3 during the 7 weeks of strong lockdown have been analysed at an air quality monitoring station focused on traffic emissions in Seville (Spain). To evaluate the changes produced in the relationships of these pollutants in 2020, a reference period from 2017 to 2019 were selected. According to the preliminarily descriptive statistical analysis, Seville, as same as other locations in the world, experienced a drastic NO 2 reduction (41% compared to the reference period) but O 3 levels remained stable and kept the expected daily cycle.
For MF-DXA approach, two multifractal regimens were found for NO 2 and O 3 with a temporal crossover separating both scaling regimens located at 1 day. While for small temporal scales (<1 day), crosscorrelation between both variables is close to an uncorrelated situation, large scales (>1 day) exhibit anti-persistency. However, no conclusive differences between the lockdown and the reference period are found in multifractal cross-correlations between NO 2 and O 3 . Thus, MF-DXA is insensitive to changes produced in the NO 2 -O 3 relationship during the lockdown. This may be due to the fact that the study period is 7 weeks, which could be a short period for a multifractal analysis based on the detection of long-range cross-correlations.
As regards MF-X-PF approach, multifractal features are detected for a broad temporal scaling range defined from 80 min to ~28 days. Joint multifractal spectra exhibit triangular shapes oriented from the top left region to the bottom right, which indicates an inversely proportional relationship between NO 2 and O 3 with no combination of low values for both variables. During 2020, a narrower range of singularity exponents for both pollutants is found which meaning that the lowest NO 2 (O 3 ) concentrations are related to a loss of variability in O 3 (NO 2 ) concentrations during the lockdown. From joint multifractal spectra are derived different scenarios of O 3 as a single variable by taking certain q moments of NO 2 . When high NO 2 concentrations are considered (q NO 2 = +3), O 3 single spectra are asymmetric with longer right tails. Specifically, in 2020 the spectrum is narrower which indicates that O 3 variability during the lockdown period was lower than in the previous years. When O 3 is studied removing the influence of NO 2 (q NO 2 = 0), the lockdown period is characterized by a quasi-monofractal behaviour as it is deduced from the accumulation of points of the spectrum, which differs from the asymmetric shape and long right tails of the reference period. When low NO 2 concentrations are considered, no differences between the lockdown and the reference period can be established.
To conclude, MF-X-PF is sensitive to O 3 variations under VOC-limited regime when descriptive statistics and MF-DXA are not able to detect them. Temporal scaling variability of O 3 has decreased during the lockdown period, mainly when high concentrations of NO 2 (q NO 2 = +3) are considered. Furthermore, when NO 2 variable is removed from the analysis (q NO 2 = 0), O 3 tends towards a monofractal behaviour. In order to better understand the inversely proportional relationship between NO 2 and O 3 , the fundamental mechanism of this trend should be investigated in a future study.

Disclosure statement
No potential conflict of interest was reported by the authors.

Funding
This work was funded by the Spanish Ministry of Science, Innovation and Universities and ERDF under contract RTI2018-100754-B-I00 (iSUN project).

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.