Mengying Li

Shaocai Yu
email: shaocaiyu@zju.edu.cn

Xue Chen

Zhen Li

Yibo Zhang

Liqiang Wang

Weiping Liu

Pengfei Li

Eric Lichtfouse

Daniel Rosenfeld

John H Seinfeld

Large scale control of surface ozone by relative humidity observed during warm seasons in China

Keywords: Ozone, Relative humidity, Temperature, Warm season, China

Rising air pollution by surface ozone (O 3 ) in China has induced extensive efforts to control ozone generation in major urban and industrial areas, yet mechanisms ruling the ozone production and loss are not well understood. In particular, ozone levels are strongly influenced by meteorological factors such as relative humidity, but this has been explored only in local situations, and the effect of relative humidity on ozone levels in warm seasons on a large scale in China is still unknown. Here we studied surface ozone, relative humidity, temperature, and other meteorological variables in 74 major cities in China during 2017-2018, focusing on the warm seasons in seven regions. Results show that ozone levels decrease with increas-ing relative humidity in all cities, with an average correlation coefficient of -0.58, ranging from -0.17 in Zhangjiakou to -0.84 in Hengshui. At high relative humidity levels, above 75%, average ozone levels ranged from 44.6 to 122.5 μg m -3 , which is lower than Chinese quality threshold of hourly average ozone level of 200 μg m -3 . The decreases of ozone with relative humidity were more pronounced at high temperature, above 30 °C, than below 25 °C. The increases of ozone with temperature were more pronounced at low relative humidity, below 40%. Overall, our findings reveal that mechanisms rul-ing surface ozone levels are similar on a large scale. This is promising to design common methods of climate engineering to protect human health.

Introduction

Tropospheric ozone (O 3 ) is a major oxidant and greenhouse gas related partly to pollution. Ozone has negative impacts on human health, regional climates, agricultural yields, and forest productivity [START_REF] Yu | Performance and diagnostic evaluation of ozone predictions by the eta-community multiscale air quality forecast system during the 2002 new England air quality study[END_REF]Shao et al. 2006;[START_REF] Booker | The ozone component of global change: Potential effects on agricultural and horticultural plant yield, product quality and interactions with invasive species[END_REF][START_REF] Monks | Tropospheric ozone and its precursors from the urban to the global scale from air quality to short-lived climate forcer[END_REF]. Anthropogenic ozone is produced by photochemical reactions of nitrogen oxides (NO x ) and volatile organic compounds (VOCs) in the presence of sunlight [START_REF] Seinfeld | Atmospheric chemistry and physics: from air pollution to climate change, 2nd edn[END_REF]. High near-surface O 3 levels occur in many urban and industrial areas across China. The establishment of an atmospheric monitoring network since 2000 in China has allowed regular ozone monitoring. Hourly monitoring data in 74 major cities disclosed an increase of daily average maximum 8-h O 3 concentrations from about 69.5 ppbv in 2013 to 75 ppbv in 2015 (http:// www. mep. gov. cn). Moreover, a 4.2% increase of mean O 3 concentrations in 336 prefecture-level cities was observed from 2015 to 2016 [START_REF] Li | Air pollution characteristics in China during 2015-2016: spatiotemporal variations and key meteorological factors[END_REF]. Recently, maximum 8-h O 3 concentrations exceeding the China National Ambient Air Quality Standard of 160 μg/m 3 (CNAAQ, GB3095-2012) have been observed in many regions in China, notably in Beijing-Tianjin-Hebei, Yangtze River Delta, and Pearl River delta [START_REF] Wang | Ozone pollution in China: a review of concentrations, meteorological influences, chemical precursors, and effects[END_REF].

Extensive efforts have been made to control surface ozone by reducing ozone precursors and controlling meteorological conditions [START_REF] Ou | Ambient ozone control in a photochemically active region: short-term despiking or long-term attainment?[END_REF][START_REF] Wang | Ozone pollution in China: a review of concentrations, meteorological influences, chemical precursors, and effects[END_REF]Chen et al. 2019a;[START_REF] Yang | Quantitative impacts of meteorology and precursor emission changes on the long-term trend of ambient ozone over the Pearl River Delta, China, and implications for ozone control strategy[END_REF]. Parameters ruling ozone levels include emissions, solar radiation, precipitation, temperature (T), relative humidity (RH), wind speed and direction, and cloud cover [START_REF] Camalier | The effects of meteorology on ozone in urban areas and their use in assessing ozone trends[END_REF][START_REF] Davis | A comparison of CMAQbased and observation-based statistical models relating ozone to meteorological parameters[END_REF][START_REF] Reddy | Meteorological factors contributing to the interannual variability of midsummer surface ozone in Colorado, Utah, and other western U.S. states[END_REF][START_REF] Wang | Ozone pollution in China: a review of concentrations, meteorological influences, chemical precursors, and effects[END_REF]Chen et al. 2019b;[START_REF] Li | Air pollution characteristics in China during 2015-2016: spatiotemporal variations and key meteorological factors[END_REF]. [START_REF] Yu | Fog geoengineering to abate local ozone pollution at ground level by enhancing air moisture[END_REF] proposed a fog geoengineering method to increase atmospheric moisture and, in turn, abate near-surface O 3 levels. In this study, we focused on the effects of relative humidity on surface O 3 .

Relationships between O 3 and relative humidity have been investigated in modeling studies. For instance, maximum daily temperature and average relative humidity had the largest effects in 74 cities in the US according to a cityspecific generalized linear model [START_REF] Davis | A comparison of CMAQbased and observation-based statistical models relating ozone to meteorological parameters[END_REF]. The role of water vapor in the reduction of surface O 3 levels has also been demonstrated by simulations with Climatechemistry models (CCMs) [START_REF] Thompson | Sensitivity of tropospheric oxidants to global chemical and climate change[END_REF][START_REF] Johnson | Role of climate feedback on methane and ozone studied with a coupled ocean-atmosphere-chemistry model[END_REF][START_REF] Murazaki | How does climate change contribute to surface ozone change over the United States?[END_REF][START_REF] Doherty | Impacts of climate change on surface ozone and intercontinental ozone pollution: a multi-model study[END_REF]. A 19% increase of water vapor could lead to approximately 4-5% decrease of O 3 concentrations [START_REF] Doherty | Impacts of climate change on surface ozone and intercontinental ozone pollution: a multi-model study[END_REF]. Similarly, the increasing relative humidity from 9 to 87% can reduce O 3 concentrations by 30% [START_REF] Wang | Effect of particle water on ozone and secondary organic aerosol formation from benzene-NO 2 -NaCl irradiations[END_REF]. The decrease of ozone with relative humidity has been explained by the following processes (Table S1):

• O 3 photolysis produces O( 1 D). O( 1 D) then reacts with a water molecule (H 2 O) to generate two OH radicals, which favor ozone decomposition. • Relative humidity favors heterogeneous reactions of O 3 with particles in aerosols (He et al. 2017a;[START_REF] Zhao | A high-resolution emission inventory of primary pollutants for the Huabei region, China[END_REF]. For instance, humidity enhance ozone uptake by increasing aerosol size and decreasing aerosol viscosity [START_REF] Han | Quantifying the role of the relative humidity-dependent physical state of organic particulate matter in the uptake of semivolatile organic molecules[END_REF]). • Ozone decomposition is enhanced by heterogeneous bonding on aqueous aerosol surfaces [START_REF] Kotelnikov | Role of aqueous aerosols in ozone decomposition in the near-surface atmosphere[END_REF]). • Wet aerosols formed under high humidity reduce ultraviolet actinic flux and inhibit photochemical reaction rates [START_REF] Deng | Modeling heterogeneous chemical processes on aerosol surface[END_REF][START_REF] Liu | Hygroscopic properties of aerosol particles at high relative humidity and their diurnal variations in the north China plain[END_REF]). • Elevated relative humidity decreases temperature by water evaporation [START_REF] Olszyna | The correlation of temperature and rural ozone levels in southeastern U.S.A[END_REF]). • Relative humidity inhibits the chain length of hydroxyl radicals (OH and HO 2 ) and NO 2 , thus limiting O 3 generation [START_REF] Reichert | Investigation of the effect of water complexes in the determination of peroxy radical ambient concentrations: implications for the atmosphere[END_REF][START_REF] Wang | Effect of particle water on ozone and secondary organic aerosol formation from benzene-NO 2 -NaCl irradiations[END_REF]).

• High relative humidity promotes the reaction of gaseous N 2 O 5 with aqueous NaCl particles, thus greatly reducing O 3 generation [START_REF] Leu | Heterogeneous reactions of HNO3(g)+NaCl(s).fwdarw. HCl(g)+NaNO3(s) and N2O5(g)+NaCl(s).fwdarw. ClNO2(g)+NaNO3(s)[END_REF][START_REF] Jia | Ozone and secondary organic aerosol formation from ethylene-NO x -NaCl irradiations under different relative humidity conditions[END_REF].

Decreases of ozone with relative humidity have been reported in some regions and are commonly observed in summer seasons [START_REF] Olszyna | The correlation of temperature and rural ozone levels in southeastern U.S.A[END_REF][START_REF] Dueas | Assessment of ozone variations and meteorological effects in an urban area in the Mediterranean Coast[END_REF][START_REF] Elminir | Dependence of urban air pollutants on meteorology[END_REF][START_REF] Camalier | The effects of meteorology on ozone in urban areas and their use in assessing ozone trends[END_REF][START_REF] Tu | Temporal variations in surface ozone and its precursors and meteorological effects at an urban site in China[END_REF]. A study in a coastal Spanish site showed that, together with temperature and wind speed, relative humidity was one of three dominant meteorological factors influencing ozone levels [START_REF] Dueas | Assessment of ozone variations and meteorological effects in an urban area in the Mediterranean Coast[END_REF]. Pearson correlation and stepwise multiple regression analysis also showed negative correlations between ozone-RH [START_REF] Manju | Spatio-seasonal variation in ambient air pollutants and influence of meteorological factors in Coimbatore Southern India[END_REF]. In Nanjing, a report shows a maximum O 3 level of 43.4 ppbv below 40% relative humidity, versus the lowest value of 10.9 ppbv above 80% RH [START_REF] Tu | Temporal variations in surface ozone and its precursors and meteorological effects at an urban site in China[END_REF]. Relative humidity was also found as the strongest influencing factor in central and southern parts of eastern China [START_REF] Han | Local and synoptic meteorological influences on daily variability in summertime surface ozone in eastern China[END_REF]. Relative humidity appeared more important than radiation and temperature to assess the adverse effects of ozone on vegetation [START_REF] Gong | A humidity-based exposure index representing ozone damage effects on vegetation[END_REF]. Overall, previous investigations focused on local situations or several cities during the same time period. Therefore, here we present a study of ozone, relative humidity and temperature in different warm seasons at a large geographical scale in the major 74 cities in China.

Experimental

Monitoring data of O 3 and its precursors (NO 2 , CO) were collected from the website of the China National Environmental Monitoring Center (CNEMC, http:// www. mee. gov. cn/ hjzl/ dqhj/). We focused on hourly O 3 concentrations measured at all monitoring sites in 74 major cities in China during 2017-2018. The monitoring data in each city consist of urban and background sites. The average O 3 concentrations at all sites in each separate city were used and regarded as the city-level O 3 values. Meteorological variable data on an hourly basis including temperature (T, ℃), relative humidity (RH, %), precipitation (Prec, mm), and wind speed (WS, m/s) were provided by the China Meteorological Administration (http:// data. cma. cn/ site/ index. html). Only hours with both valid O 3 and RH values were used for analysis. The data of leaf area indices (LAI) in the summer seasons and topography height (HGT_M) in China were derived from MODIS products (https:// modis. gsfc. nasa. gov/ data/).

Under the comprehensive consideration of the spatial patterns of O 3 and national administrative divisions in China, the following 74 cities were geographically grouped into seven representative divisions, namely Northeast China (Dalian, Haerbin, Shenyang, Changchun), North China (Baoding, Beijing, Cangzhou, Chengde, Handan, Hengshui, Huhehaote, Langfang, Qinhuangdao, Shijiazhuang, Taiyuan, Tangshan, Tianjin, Xingtai, Zhangjiakou), Central China (Wuhan, Changsha, Zhengzhou), East China (Changzhou, Fuzhou, Hangzhou, Hefei, Huzhou, Huaian, Jinan, Jiaxing, Jinhua, Lishui, Lianyungang, Nanchang, Nanjing, Nantong, Ningbo, Qingdao, Quzhou, Xiamen, Shanghai, Shaoxing, Suzhou, Taizhouz, Taizhouj, Wenzhou, Wuxi, Suqian, Xuzhou, Yancheng, Yangzhou, Zhenjiang, Zhoushan), South China (Dongguan, Foshan, Guangzhou, Haikou, Huizhou, Jiangmen, Nanning, Shenzhen, Zhaoqing, Zhongshan, Zhuhai), Southwest China (Chengdu, Guiyang, Kunming, Chongqing, Lasa), and Northwest China (Lanzhou, Wulumuqi, Xian, Xining, Yinchuan), as seen in Fig. 1 and Table S2. Then, concentrations and meteorological data averaged across all cities located in each division were used. While conducting correlation analyses between O 3 and RH, we chose different study periods for the seven regions due to various warm seasons for each region, as based on monthly variations of temperature in Fig. 1. The national warm season typically spans April-September. But obviously the northern regions have lower temperatures than southern regions in May and September. So summer season (June-August) is regarded as the warm season in northern areas in China. Since O 3 concentrations in South China typically peak in October [START_REF] Lu | Severe surface ozone pollution in China: a global perspective[END_REF], also shown in Fig. 1, the study period of cities in South China was expanded to October. Therefore, the final chosen study periods were June to August for Northeast China (NE), North China (NC), and Northwest China (NW), May to September for East China (EC), Central China (CC), and Southwest China (SW), May to October for South China (SC). Due to geographical proximity of Fuzhou and Xiamen cities to South China, time periods of these two cities were kept in correspondence with the cities in South China.

Results and discussion

Seasonal variations of O 3 , NO 2 and CO

Average levels of ozone (O 3 ), nitrogen dioxide (NO 2 ), carbon monoxide (CO), relative humidity (RH), and temperature (T) in various Chinese regions in 2018 are shown in Fig. 1, with warm seasons shadowed in orange. Results show that annual mean O 3 concentrations are heterogeneous in 74 cities, with values ranging from 65.8 μg m -3 in Haikou, Hainan province, to 143.8 μg m -3 in Cangzhou, Hebei province. The highest O 3 levels were mainly observed in North and East China and have been explained by greater emissions of ozone precursors in industrialized and urbanized regions. Indeed, in North China, power plants, coal combustion industries, and biomass burning sources produce more NO x and VOCs emissions [START_REF] Zhao | A high-resolution emission inventory of primary pollutants for the Huabei region, China[END_REF][START_REF] Chai | Spatial and temporal variation of particulate matter and gaseous pollutants in 26 cities in China[END_REF]. Higher NO 2 levels are also observed in cities having more vehicle emission sources (Wang et al. 2010a, b).

NO 2 and CO levels in all Chinese regions display similar variations, with maximum values in winter and late autumn, and minimum values in summer, late spring, and early autumn (Fig. 1). The low concentrations of these ozone precursors in warmer seasons can be attributed to active photochemistry and stronger vertical atmospheric mixing. Whereas the high levels in colder seasons may be a result of larger consumption of heating fuels, as well as weak solar radiation, low wind speed and mixing layer height, enhancing NO 2 and CO concentrations (Wang et al. 2010a;[START_REF] Tai | Correlations between fine particulate matter (PM 2.5 ) and meteorological variables in the United States: implications for the sensitivity of PM 2.5 to climate change[END_REF].

Considerable seasonal variations of O 3 are observed in all Chinese regions, with high levels in summer, late spring, and early autumn, and low levels in winter (Fig. 1). Maximum ozone concentrations range from 73.1 μg m -3 in South China to 117.2 μg m -3 in North China. In general, monthly variations of O 3 can be explained by prevailing meteorological conditions and anthropogenic activities [START_REF] Reddy | Meteorological factors contributing to the interannual variability of midsummer surface ozone in Colorado, Utah, and other western U.S. states[END_REF]. Peak values of ground-level O 3 tend to appear in the warm seasons with intense solar radiation, high temperature, stagnant high-pressure conditions, as observed previously [START_REF] Ou | Speciated OVOC and VOC emission inventories and their implications for reactivity-based ozone control strategy in the pearl river delta region, China[END_REF]He et al. 2017b).

Nevertheless, we also observed inconsistent fluctuations in variation patterns. For example, monthly average O 3 concentrations displayed similar variations in six regions, peaking in May or June. Whereas South China displayed a seasonal ozone maximum in October when temperature was not maximum. This is likely due to the sharp decline of relative humidity from about 80% in August to 65% in October in South China. Similar findings have been shown in previous studies conducted in the Pearl River delta, with October being the most polluted month, different from other regions [START_REF] Bell | A meta-analysis of time-series studies of ozone and mortality with comparison to the national morbidity, mortality, and air pollution study[END_REF][START_REF] Zanobetti | Mortality displacement in the association of ozone with mortality: an analysis of 48 cities in the United States[END_REF]. This phenomenon has also been attributed to the adverse effect of the Asian summer monsoon on O 3 , with O 3 peaking before and after monsoon in the Yangtze River and Pearl River deltas [START_REF] Lu | Severe surface ozone pollution in China: a global perspective[END_REF][START_REF] Han | Local and synoptic meteorological influences on daily variability in summertime surface ozone in eastern China[END_REF]. Although higher temperature and longer sunshine duration in summer favors O 3 formation in South China, larger rainfall amounts in summer versus autumn can retard O 3 formation [START_REF] Shao | Ground-level ozone in the pearl river delta and the roles of VOC and NOx in its production[END_REF].

Overall, our findings show that high ground-level O 3 levels are generally accompanied by high temperature and low relative humidity in all regions, whereas high NO 2 and CO concentrations appear under low temperature and high relative humidity. Since ozone peaked in warm seasons, we focused further research on warm seasons highlighted in orange in Fig. 1 to investigate relationships between ozone and relative humidity.

Correlation between ozone and relative humidity

Figure 2 displays correlations of O 3 versus relative humidity (RH) during the noontimes, 11:00-16:00 local time, during warm seasons of 2017 and 2018 in the major 74 cities in China. We excluded data of rainfall time to avoid the influence of precipitation scavenging. We selected noontimes when O 3 formation is higher due to stronger solar radiation, according to [START_REF] Monks | Tropospheric ozone and its precursors from the urban to the global scale from air quality to short-lived climate forcer[END_REF]. Moreover, data with low relative humidity, below 25% and temperature, below 15 °C, were excluded owing to little effect of low humidity and temperature on ozone. Correlation coefficient (R) and slope of linear regression equation (S) were used to represent the relationships. Unlike prior studies choosing the same time periods for all cities, here we focused on dates of warm seasons, which are different in the seven regions.

Results show an overall decrease of ozone levels with relative humidity in all major cities, with R coefficient ranging from -0.84 in Hengshui city to -0.17 in Zhangjiakou city, and slope values varying widely from -2.57 to -0.48 µg m -3 /1%RH (Fig. 2). The average R coefficient for all cities is -0.58 ± 0.14. Moreover, hourly O 3 levels are strongly correlated with RH in most cities, with correlation coefficient below -0.5 in 57 cities. Also, we did not observe significant heterogeneity of correlations between Southern and Northern regions, contrary to previous research in China and abroad [START_REF] Camalier | The effects of meteorology on ozone in urban areas and their use in assessing ozone trends[END_REF][START_REF] Kavassalis | Understanding ozone-meteorology correlations: a role for dry deposition[END_REF][START_REF] Cheng | Regionalization based on spatial and seasonal variation in ground-level ozone concentrations across China[END_REF]). This discrepancy is explained by the fact that analyses of previous investigations compared data at the same dates for different locations, whereas here we used only data of warm seasons, which occur at different dates in different regions. For example, greater responses were observed in the southern versus northern areas in 39 US cities from May to September [START_REF] Camalier | The effects of meteorology on ozone in urban areas and their use in assessing ozone trends[END_REF]). We also studied the effect of other meteorological factors on ozone variations (Fig. S5). Results show that the relative humidity effect is comparable to or even stronger than temperature, and far exceeds wind speed. In summary, we found that almost all cities have large negative correlations between RH and ozone in warm seasons.

The different influences of relative humidity on ozone observed in different cities and regions (Fig. 2) can be explained by chemical and meteorological drivers. For instance, low temperatures in Southwestern areas throughout the year reduces O 3 formation, while high radiation intensity and long solar duration favors O 3 production via enhancing the photolysis rates of NO 2 and reaction rates of NO and O 3 [START_REF] Gaur | Four-year measurements of trace gases (SO 2 , NO x , CO, and O 3 ) at an urban location, Kanpur, in Northern India[END_REF]. Also, long-range transport of pollutants from high emission areas to low emission regions can change precursor concentrations and further affect ozone-RH correlations.

Moreover, rising relative humidity increases the opening of tree stomata and, in turn, should induce the removal of ozone from air as a dry ozone sink [START_REF] Kavassalis | Understanding ozone-meteorology correlations: a role for dry deposition[END_REF]. To check this hypothesis, we tested the correlation of leaf area indices (LAI) with correlation coefficients of ozone and RH (Figure S1, Table S2). Our results indicate that the dry deposition of ozone in trees do not play a significant influence on the strong negative correlations between ozone and RH, in our conditions.

Table S3 shows the effects of high RH on O 3 decreases in all cities. For RH over 75%, the average O 3 concentrations varied from 44.6 to 122.5 μg m -3 in these 74 cities, and the 90% percentile of hourly O 3 concentrations ranged from 75.6 to 168.7 μg m -3 , less than the quality standard of 200 μg m -3 . For O 3 concentrations exceeding 200 μg m -3 , the mean RH values ranged from 31 to 66%, with 90% percentile of hourly RH values varying from 36 to 81%. This shows that O 3 quality exceedance generally appeared under relatively low RH. In summary, our findings show that high RH exerted great influence on O 3 attainment in warm seasons for all 74 cities.

Overall, our findings show that ozone levels decreases with relative humidity in all Chinese regions, that there is no significant differences between southern and northern regions if only warm seasons are considered, and that dry ozone deposition on trees is not a major influencing factor in our conditions. Further effects of meteorological factors are studied below.

Meteorological and chemical factors

We studied the effect of temperature, wind speed, surface pressure, NO 2 , and particulate matter (PM 2.5 ) concentrations, on ozone levels (Fig. 3a). Results show that O 3 concentrations display larger variation ranges with RH at high temperature, low wind speed, high pressure, and high NO 2 and PM 2.5 levels. Ozone can exceed 160 μg m -3 below 40% low RH, with PM 2.5 over 35 μg m -3 and NO 2 over 20 μg m -3 . This suggests the inducing effect of PM 2.5 and NO 2 on ozone, and implies that adverse health effects should be reduced by decreasing particulate matter and NO 2 emissions. On the other hand, when relative humidity exceeds 60%, O 3 concentrations are lower than 130 μg m -3 under all conditions, except high PM 2.5 levels, above 75 μg m -3 . We conducted a stratified analysis at four representative cities, Baoding in North China, Zhengzhou in Central China, Hangzhou in East China, and Chengdu in Southwest China to gain more knowledge on ozone response to relative humidity and temperature. These cities show high O 3 pollution, high negative correlations between O 3 and RH, with R below -0.5, and a large distribution range of hourly temperature and relative humidity data. Correlations of O 3 versus RH at different temperature ranges are given in Figs. 3b,c, S4a, and S4b. Correlations of ozone versus temperature, grouped by relative humidity, are given in Figs. 3c,d,S4c, and S4d. Results show that the negative relationship between O 3 and RH became stronger under higher temperature in Baoding city, with an average decreasing rate of -1.73 μg m -3 /1%RH above 30 °C compared to -0.29 μgm -3 /1%RH below 25℃. Similar variations were also found in Zhengzhou, Hangzhou and Chengdu with larger decreasing rates of O 3 as RH rose at higher temperature. This result indicates that the increase of relative humidity is effective to decrease O 3 concentrations at high temperature. The choice of the temperature range is therefore important to control O 3 levels in mitigation practices.

We also observed that stronger increasing rates of O 3 with temperature occurred under low relative humidity, below 40%, with 11.6, 10.86, 2.78, and 9.19 μg m -3 /°C in Baoding, Zhengzhou, Hangzhou, and Chengdu, respectively, much higher than 0.85, 0.73, 0.04, and 1.70 μg m -3 /°C in corresponding cities at high relative humidity, over 80%. Although under intermediate relative humidity, from 40 to 80%, increasing rates of O 3 were comparable to and even greater in Hangzhou and Chengdu than under low RH, the weakening effects of high RH, above 80%, on O 3 increase were evident in these two cities. This means that relative humidity can not only decrease ozone levels but also alleviate the positive effect of temperature on ozone generation. In summary, we found that high temperature strengthens the effects of relative humidity on ozone decrease, and that high relative humidity weakens the positive correlations between O 3 and temperature.

Conclusion

On the basis of observed pollutant concentrations and meteorological variables for the major 74 cities in China, we first studied seasonal variation patterns of O 3 , then conducted correlation analyses of O 3 and relative humidity during different warm seasons in 74 cities. Furthermore, the study discussed the effects of relative humidity and temperature on O 3 using stratified analyses aimed at four cities, i.e., Baoding, Zhengzhou, Hangzhou, Chengdu. O 3 levels displayed spatial differences among different regions, with higher levels mainly distributed in North China and East China. A considerable seasonal variation of O 3 was recorded in each region, showing peak values in warm seasons with intense solar radiation, stagnant high-pressure conditions. We found the overall negative correlations between O 3 and relative humidity in all cities, indicating the explanatory power of relative humidity to hourly variability of midday O 3 in the warm seasons. Also, the role of high relative humidity in O 3 attainment can be seen at all the cities. Then the stratified analyses show that high temperature can strengthen O 3 dependence on relative humidity, and high relative humidity can weaken the positive impacts of temperature on O 3 . Our results may have an implication for the great need of controlling the rising O 3 pollution, especially in regions having drier climate conditions with low relative humidity and soil water content.

Fig. 1

 1 Fig. 1 Monthly variations of average concentrations of ozone (O 3 ), nitrogen dioxide (NO 2 ), and carbon monoxide (CO), and values of relative humidity (RH), and temperature (T) for Chinese areas in 2018. The orange-shaded areas denote chosen warm seasons: June to August for Northeast China (NE), North China (NC), and Northwest China (NW), May to September for East China (EC), Central China (CC), and Southwest China (SW), May to October for South China

Fig. 2

 2 Fig. 2 Correlations between noontime O 3 and relative humidity (RH), 11:00-16:00 local time, excluding data points with precipitation and low RH and T (RH ≤ 25%, T ≤ 15℃) in warm seasons of 2017 and 2018 for the major 74 cities in China. Results without exclusion are

Fig. 3 a

 3 Fig.3a The average O 3 concentrations for each relative humidity (RH) bin (5% intervals) in all the major cities grouped by temperature (T), wind speed (WS), pressure (P), NO 2 , and PM 2.5 concentrations (units: °c, m s -1 , hPa, μg m -3 , μg m -3 , respectively). Scatter plots of O 3 against RH grouped by T during the noontime (11:00-16:00, local time) excluding precipitation in different warm seasons during 2017-2018 in b Baoding and c Zhengzhou. Points represent the average values of O 3 concentration versus RH for each RH bin (5% intervals). Scatter plots of O 3 against T grouped by RH during the same period as above in d Baoding and e Zhengzhou. Points represent the average values of O 3 concentrations versus T for each T bin(2 °C intervals). The slopes of linear regression equations are also provided in the plots (α < 0.05). Results show that high temperature strengthened negative dependence of O 3 on RH, and high relative humidity not only lowered O 3 levels, but also weakened the positive effects of T on O 3 ◂

Fig

  Fig. S3 Correlations between noontime O3 and relative humidity (RH) (11:00 to 16:00, local time) by excluding data points

  TableS1Related water vapor-involved chemical reactions of O3 formation[START_REF] Yu | Fog geoengineering to abate local ozone pollution at ground level by enhancing air moisture[END_REF].

	O3 decomposition by water vapor hydroxyl radicals termination reactions	NO2 termination reactions
	O3 + hv → O( 1 D) + O2	HO2 + CO + nH2O → CO2 + OH	O( 3 P) + NO2 + M → NO3 + M
	O( 1 D) + H2O → 2OH	HO2 + NO + nH2O → nonradical products	NO2 + O3 → NO3 + O2
	O3 + OH → HO2 + O2	OH + CO + nH2O →nonradical products	NO2 + NO3 + M → N2O5 + M
	O3 + HO2 → OH + 2O2		2NO2 + H2O → HONO + HNO3
			N2O5 + H2O → 2HNO3
			N2O5(gas) + NaCl → ClNO2 +NaNO3

Acknowledgements

This work was supported in part by the Department of Science and Technology of China (No. 2016YFC0202702, 2018YFC0213506, and 2018YFC0213503), National Research Program for Key Issues in Air Pollution Control in China (No. DQGG0107), and National Natural Science Foundation of China (No. 21577126 and 41561144004). Part of this work was also supported by the "Zhejiang 1000 Talent Plan" and Research Center for Air Pollution and Health in Zhejiang University. Pengfei Li is supported by National Natural Science Foundation of China (No. 22006030), Initiation Fund for Introducing Talents of Hebei Agricultural University (412201904), and Hebei Youth Top Fund (BJ2020032).

Supplementary Information

The online version contains supplementary material available at https:// doi. org/ 10. 1007/ s10311-021-01265-0.

Supplementary Information


domain: stat sde math

text:
Philom Ène

Le Gall

Anne-Catherine Favre

Philippe Naveau

Alexandre Tuel

Non-parametric multimodel Regional Frequency Analysis applied to climate change detection and attribution

Keywords: CMIP, detection & attribution, extreme precipitation, F-madogram, Regional Frequency Analysis

A recurrent question in climate risk analysis is determining how climate change will affect heavy precipitation patterns. Dividing the globe into homogeneous sub-regions should improve the modelling of heavy precipitation by inferring common regional distributional parameters. In addition, in the detection and attribution (D&A) field, biases due to model errors in global climate models (GCMs) should be considered to attribute the anthropogenic forcing effect. Within this D&A context, we propose an efficient clustering algorithm that, compared to classical regional frequency analysis (RFA) techniques, is covariate-free and accounts for dependence. It is based on a new non-parametric dissimilarity that combines both the RFA constraint and the pairwise dependence. We derive asymptotic properties of our dissimilarity estimator, and we interpret it for generalised extreme value distributed pairs. As a D&A application, we cluster annual daily precipitation maxima of 16 GCMs from the coupled model intercomparison project. We combine the climatologically consistent subregions identified for all GCMs. This improves the spatial clusters coherence and outperforms methods either based on margins or on dependence. Finally, by comparing the natural forcings partition with the one with all forcings, we assess the impact of anthropogenic forcing on precipitation extreme patterns.

Introduction

Since the early 19 th century, fossil fuels-based human activities have become one of the major forces of ecosystem and climate change, defining a new geological era, called Anthropocene [START_REF] Crutzen | The "anthropocene[END_REF] or Capitalocene [START_REF] Malm | The geology of mankind? A critique of the Anthropocene narrative[END_REF][START_REF] Campagne | Le capitalocène: aux racines historiques du dérèglement climatique[END_REF]. The global warming caused by these activities induces important changes in the climate system (IPCC, 2021). Working Group I of the IPCC, which assesses the physical science of climate change, summarises the latest advances in climate science to understand the climate system and assess climate change, by combining data from paleoclimate, observations and global circulation model (GCM) simulations. The latter are based on differential equations linked to the fundamental laws of physics, thermodynamics and chemistry. GCMs simulate the evolution of various climate variables on discretised tridimensional meshes with a typical horizontal resolution of 100 [km] or more. The coupled model intercomparison project (CMIP) [START_REF] Meehl | The coupled model intercomparison project (CMIP)[END_REF][START_REF] Alexander | Historical and projected trends in temperature and precipitation extremes in Australia in observations and CMIP5[END_REF] aims at comparing the performances of several dozen of GCMs developed by different research centres, e.g. see Table 1 in Appendix. As numerical experiments and approximations of the true climate system, these GCMs can produce different climate responses to different given inputs, e.g. emission scenarios. To reduce model errors and gain robustness in signal detection, GCMs are often analysed jointly. In particular, CMIP models have been used in the field of detection and attribution that aims at finding causual links between the climate response and known external forcings (see, e.g. [START_REF] Ribes | Making climate projections conditional on historical observations[END_REF][START_REF] Van Oldenborgh | Pathways and pitfalls in extreme event attribution[END_REF][START_REF] Naveau | Statistical methods for extreme event attribution in climate science[END_REF]. As a yardstick, the so-called "natural forcings" runs have not been influenced by human activities and were only driven by external forcings, e.g. solar variations, explosive volcanic eruptions like Mont Pinatubo in 1991 (see, e.g. [START_REF] Ammann | Statistical volcanic forcing scenario generator for climate simulations[END_REF]. Such a numerical setup can be viewed as a thought-experiment and it corresponds to a counterfactual world, but not to the observed one. In contrast, a factual world is produced by integrating all forcings, including rising greenhouse gazes, and factual runs aims at reproducing the observed climatology over the last century. Future periods, say 2071-2100, can also be explored with GCMs but future forcing and emission scenarios need to be chosen. For example, RCP8.5 for CMIP5 (IPCC, 2013) and SSP5-8.5 for CMIP6 (IPCC, 2021) will be analysed in this paper. In this context, a natural question is to wonder how the climate system will change under these scenarios.

Due to their large societal and economical impacts, a vast literature has be dedicated to answering this question for extreme events. In particular, heavy rainfall and heatwaves have received a particular attention, see chapters 10 and 11 in the Working Group I contribution of IPCC (2021) report. In this paper, we focus on annual maxima of daily precipitation from 1850 to 2100 provided by the factual (all forcings) and counterfactual (natural forcings only) models listed in Table 1 of the Appendix. Note that our main climatological goal is not to directly assess changes in heavy rainfall intensities and frequencies, but rather to detect how spatial patterns (clusters) of yearly maxima of daily precipitation could be modified by anthropogenic forcing.

To model yearly block maxima, one classical statistical approach is to impose a parametric generalised extreme value (GEV) distributions (see e.g. [START_REF] Coles | An introduction to statistical modeling of extreme values[END_REF][START_REF] Davison | Statistical modeling of spatial extremes[END_REF]. For example, each grid point of each individual CMIP model could be fitted with a spatial structure embedded within the GEV parameters (see, e.g. [START_REF] Kharin | Changes in temperature and precipitation extremes in the cmip5 ensemble[END_REF]. However, the computational cost can be high (more than 200 years of precipitation data at thousands of grid points for 16 models), especially if the spatial dependence is included. Another aspect is the ease of interpretation. Well defined spatial patterns (clusters) in extreme precipitation are very useful for climatologists who can interpret them according to known physical phenomena (e.g., [START_REF] Pfahl | Understanding the regional pattern of projected future changes in extreme precipitation[END_REF][START_REF] Tandon | Understanding the dynamics of future changes in extreme precipitation intensity[END_REF][START_REF] Dong | Attribution of extreme precipitation with updated observations and cmip6 simulations[END_REF]. For example, the so-called regional frequency analysis (RFA) has been frequently used in hydrology, see [START_REF] Dalrymple | Flood-frequency analyses, manual of hydrology: Part 3[END_REF]; [START_REF] Hosking | Regional frequency analysis: an approach based on L-moments[END_REF], but it has been rarely implemented in a D&A context, especially within the CMIP repository. The main idea of RFA is to identity homogeneous regions with identical distributional features, up to normalising constants. More precisely two positive absolutely continuous random variables (r.v.) Y 1 and Y 2 are said to be homogeneous if there exists a positive constant λ such that

Y 2 d = λY 1 ,
where d = denotes equality in distribution. This condition can be reformulated in terms of their cumulative distribution functions (cdf)

F i (x) = P(Y i ≤ x) with i ∈ {1, 2} as F 2 (λx) = F 1 (x).
(1)

Hence, two climate model grid points are said to belong to the same homogeneous region if they satisfy (1). To visually understand this condition within the CMIP archive, three grid points, say A, B and C, from the CCSM4 counterfactual run are plotted in panel (a) of Figure 1. In panel (b), ranked annual precipitation maxima (rescaled by the empirical mean) of point A are compared to the ones from point B. Panel (d) provides the same information but between point A and point C. It appears that points A and B are likely to satisfy (1) and, consequently, could belong to the same homogeneous region. In contrast, the rescaled distribution at point A is much more heavy-tailed than at point C. This is not surprising because A and B are nearby and C far away from them. Still, panels (b) and (d) only rely on the marginal behaviours, and pairwise dependence information and/or covariates could help finding of homogeneous regions.

Various RFA techniques based on explanatory covariates (e.g., see [START_REF] Asadi | Optimal regionalization of extreme value distributions for flood estimation[END_REF]Fawad et al., 2018, for recent work) have been developed to identify homogeneous regions which rely on station location features and/or weather patterns to explain precipitation spatial distributions (see e.g. [START_REF] Burn | Evaluation of regional flood frequency analysis with a region of influence approach[END_REF][START_REF] Hosking | Regional frequency analysis: an approach based on L-moments[END_REF][START_REF] Evin | A regional model for extreme rainfall based on weather patterns subsampling[END_REF]. For example, [START_REF] Toreti | Precipitation extremes in the mediterranean region and associated upper-level synoptic-scale flow structures[END_REF] let scale parameters vary as a function of weather station locations. However, selecting relevant covariates is constrained by their availability, expert subjectivity and the scale of the problem. In particular, finding appropriate covariates for heavy rainfall patterns at the global scale is tedious. In addition, assessing the homogeneity of regions [START_REF] Hosking | Regional frequency analysis: an approach based on L-moments[END_REF] relies on specific moments like skewness and kurtosis that are not necessary robust (based on the spatial independence assumption). Other techniques bypass the use of covariates by only working with the data at hand, here precipitation [START_REF] Saf | Regional flood frequency analysis using L-moments for the West Mediterranean region of Turkey[END_REF]. For example, Le [START_REF] Gall | Improved regional frequency analysis of rainfall data[END_REF] considered a ratio of probability weighted moments, see [START_REF] Greenwood | Probability weighted moments: definition and relation to parameters of several distributions expressable in inverse form[END_REF] and applied a clustering algorithm on this ratio. More precisely, this ratio, denoted ω ∈ [0, 1], is mean and scale invariant, i.e. in compliance with (1), and it is a simple increasing function of ξ when rainfall extremes can be assumed to either follow a GEV or Pareto distribution with shape parameter ξ. To illustrate the spatial variability of CMIP rainfall tail index (i.e. of ω), panel (a) of Figure 2 displays the ratio ω at each grid point of a counterfactual CCSM4 annual maxima run. Note that grid points A and B exhibit similar ω estimates, while grid point C differs (lighter tail). All aforementioned RFA techniques has one major drawback. They rely on the assumption of pairwise independence or pairwise conditional independence (given the co-variates). Note that Eq.(1) also constraints the marginal behaviour, but does not take into account of any information about the spatial dependence strength. Still, precipitation series at two nearby grid points are likely to be dependent. To illustrate this point, we can go back to Figure 2. Panels (b) and (e) display the scatter plots (rescaled by their means) between points A and B, and between points A and C, respectively. As expected from their local proximity, not only A and B have same similar marginals, but annual maxima of daily precipipation appears to be strongly correlated. This information coupled with constraint (1) should play an important role in improving RFA methods.

Modelling the dependence structure in clustering algorithms can be handled in different ways depending on the assumptions one is ready to make. Fully non-parametric or parametric approaches can be developed. Explanatory covariates can be included or difficult to find. For example, [START_REF] Kim | Spatial cluster detection in mobility networks: a copula approach[END_REF] introduced a parametric approach based on copulas in the context of cluster detection in mobility networks. They grouped sites subject to intense traffic according to covariates (e.g. geographical), and checked the dependence strength within each cluster by fitting a multivariate Gumbel copula. [START_REF] Drees | Principal component analysis for multivariate extremes[END_REF] and [START_REF] Janßen | k-means clustering of extremes[END_REF] proposed approaches based on exceedances; after projecting observations onto the unit sphere, they reduced their dimension through K-means clustering [START_REF] Janßen | k-means clustering of extremes[END_REF] and principal component analysis [START_REF] Drees | Principal component analysis for multivariate extremes[END_REF]. Finally, [START_REF] Bernard | Clustering of maxima: Spatial dependencies among heavy rainfall in France[END_REF] applied a non-parametric approach based on the F-madogram to weekly precipitation maxima. The so-called F-madogram [START_REF] Cooley | Variograms for spatial max-stable random fields[END_REF] is defined by

d = 1 2 E |F 1 (Y 1 ) -F 2 (Y 2 )| , (2)
where Y i is the continuous r.v. with cdf F i . It is a distance which, by construction, is marginal-free because the r.v. F 1 (Y 1 ) and F 2 (Y 2 ) are both uniformly distributed on [0, 1]. Note that if Y 1 and Y 2 are equal in probability, the distance d = 0. Whenever the bivariate vector (Y 1 , Y 2 ) follows a bivariate GEV distribution (see e.g. [START_REF] Gumbel | Distributions des valeurs extremes en plusieurs dimensions[END_REF][START_REF] Tawn | Bivariate extreme value theory: models and estimation[END_REF], this distance can be interpreted as linear transformation of the extremal coefficient (see e.g. [START_REF] Cooley | Variograms for spatial max-stable random fields[END_REF][START_REF] Naveau | Modelling pairwise dependence of maxima in space[END_REF], and Section 2.2). [START_REF] Bernard | Clustering of maxima: Spatial dependencies among heavy rainfall in France[END_REF], [START_REF] Bador | Spatial clustering of summer temperature maxima from the CNRM-CM5 climate model ensembles & E-OBS over Europe[END_REF] and later [START_REF] Saunders | A regionalisation approach for rainfall based on extremal dependence[END_REF] computed this distance to build a pairwise dissimilarity matrix that was used as an input of a clustering algorithm. In these two former studies, a partitioning around medoids (PAM) algorithm [START_REF] Kaufman | Finding groups in data: an introduction to cluster analysis[END_REF] was applied whereas the latter used hierarchical clustering. But, the RFA requirement defined by (1) was not imposed, and so the marginal differences between Y 1 and Y 2 were not taken into account. To visualise this issue within the CMIP repository, it is simple to cluster a counterfactual CCSM4 annual maxima run with the PAM algorithm † based on the distance d. The resulting map displayed in panel (b) of Figure 2 shows a few spatially coherent structures, but, overall is very patchy. In †In all our CMIP analysis, PAM was applied separately to the southern and northern hemispheres. Global analysis (available upon request) were also made, but the climatological interpretation was not as clear as with the hemispheric scale. Also, different numbers of clusters were investigated and basic criteria like the silhouette coefficient were computed. No particular number could be clearly identified. But, in terms of interpretation, four clusters appear as a reasonable compromise between climate understanding, visual simplicity and statistical criteria.

addition, panel (a) related to the marginals behaviour appears to be unrelated to panel (b) that describes the spatial dependence. This was expected from the F-madogram distance, but it would make sense to cluster grid points that are both correlated but also the same type of marginal, see (1), the essence of the RFA.

To reach this goal, we propose the following work plan. In Section 2, we integrate the homogeneity condition (1) into a new definition of the F-madogram distance. The properties of this new dissimilarity, which we call RFA-madogram, is explained by analysing a special case: the logistic bivariate GEV model in Section 2.2. A non-parametic estimator of the RFA-madogram is proposed and its asymptotic consistency in law is detailed in Section 3. Concerning the CMIP database, we compute, in Section 4, a RFA-madogram dissimilarity matrix on annual maxima of daily precipitation for each CMIP models listed in Table 1, and then cluster them with the PAM algorithm. Finally, we propose a method to build a "central" partition that summarises the partitions obtained for each model and compare the spatial patterns obtained for counterfactual and factual (2071-2100) experiments. Section 5 concludes the paper by providing a short discussion.

Joint modelling of dependence and homogeneity

RFA-madogram

To introduce homogeneity criteria, see Eq.( 1), into distance defined in Eq.( 2), we propose to define and study the following expectation

D(c, Y 1 , Y 2 ) = 1 2 E F 2 (cY 1 ) -F 1 Y 2 c , (3)
where c > 0 is a normalising positive constant. The D(c, Y 1 , Y 2 ) is always non-negative and equal to zero for c = λ when Y 2 a.s = λY 1 . The homogeneous regions are not defined a priori, so the existence of λ and its value are not known. We denote 3) is that, under the homogeneity condition of Eq.( 1),

c * 12 = argmin{D(c, Y 1 , Y 2 ) : c > 0}. Note that D(c, Y 1 , Y 2 ) = D 1 c , Y 2 , Y 1 ,
D(λ, Y 1 , Y 2 ) = d(Y 1 , Y 2 ),
where d is the classical F-madogram, see Eq.(2). To simplify notations, D or D(c) will be a shortcut for

D(c, Y 1 , Y 2 ).
The key point from a RFA point of view is that, if Eq.( 1) is satisfied, D behaves as the classical F-madogram distance. Note that D is not a true distance, but a dissimilarity. The triangle inequality is satisfied under homogeneity condition but may not be valid in general. Still, D captures information about the extremal dependence like the F-madogram, and, in addition, it encapsulates marginal information concerning the departure from Eq.(1). More precisely, one can show (see Appendix A for the proof)

that 2 |d -D| ≤ E [∆(c, Y 1 )] + E [∆(c, Y 2 /c)] , (4)
where the function ∆(c, x) = |F 2 (cx) -F 1 (x)| measures the difference between the rescaled cdfs.

To deepen our understanding of D, we comment on the special case of a bivariate-GEV distributions.

RFA-madogram for bivariate GEVs

In this section, we suppose that the bivariate vector (Y 1 , Y 2 ) follows a max-stable distribution [START_REF] Coles | An introduction to statistical modeling of extreme values[END_REF][START_REF] Fougères | Multivariate extremes[END_REF][START_REF] Guillou | Madogram and asymptotic independence among maxima[END_REF] with dependence function V (., .)

P(Y 1 ≤ x; Y 2 ≤ y) = exp -V -1 log F 1 (x) , -1 log F 2 (y)
,

where

F i corresponds to a GEV marginal cdf. If F i (x) = exp - x σ i -1/ξi with ξ 1 = ξ 2 = ξ, then the equality Y 2 d = σ 2 σ 1
Y 1 holds and we are in the homogeneity case.The shape parameter ξ describes the common upper-tail behaviour. The larger ξ is, the heavier the upper-tail of the distribution. Although complex, Eq. ( 8) in Appendix C, summarises how D(c) can be expressed in function of V (., .) and the marginal parameters.

To simplify the dependence strength interpretation, it is common to focus on the extremal coefficient defined as the scalar θ 12 such that

P(Y 1 ≤ u, Y 2 ≤ u) = {P(Y 1 ≤ u)P(Y 2 ≤ u)} θ 12 2 .
It is equal to θ 12 = V (1, 1). If Y 1 and Y 2 are independent, then θ 12 = 2, while if they are fully dependent, then θ 12 = 1. Appendix D provides the mathematical details to link the extremal coefficient with D(c). It allows to find an optimal value for rescaling parameter c * 12 . For example, it is possible to show that c * 12 =

σ 2 σ 1 = λ. for the logistic GEV model, V (x, y) = x -1 α + y -1 α α
, with α > 0.

(5)

In particular, the value of the dissimilarity D(c * 12 ) can be plotted as a function of the logistic coefficient α and of the ratio ξ 1 /ξ 2 . From Figure 3, one can see that the full dependence case corresponds to α ≈ 0, and the independence case to α = 1. In addition, the ratio ξ 1 /ξ 2 varies between 1 (homogeneity case) and 10, i.e. cases with ξ 1 = 0.1 and ξ 2 = 0.01. The dissimilarity is small when both the dependence is strong and the marginals are homogeneous (leftmost corner). Large dissimilarities correspond to the opposite cases, a near independence and/or strong heterogeneity in the shape parameters (rightmost corner). Note also that, as the homogeneity and the dependence strength decrease jointly, dissimilarity increases (concavity of the surface). These features correspond to our goal that, given the same dependence strength, the price to pay is high when the RFA condition (1) does not hold. In other words, our aim to cluster grid points that are jointly strongly dependent and in compliance with (1) seems, at least conceptually, to have been reached. The remaining question is to know if this strategy works in practice with the CMIP archive. To answer this, we need to first check that a non-parametric estimator can be developed and its asymptotic properties can be well understood.

RFA-madogram Inference

Given X ⊂ R n and n ∈ N, let ∞ (X ) denote the spaces of bounded real-valued functions on X . For f :

X → R, let f ∞ = sup x∈X |f (x)|. The arrows " a.s.
-→", "⇒", and ";" denote almost sure convergence, convergence in distribution of random vectors (see Vaart, 1998, Ch. 2) and weak convergence of functions in ∞ (X ) (see Vaart, 1998, Ch. 18-19), respectively. Let L 2 (X ) denote the Hilbert space of square-integrable functions f: X → R, with X equipped with n-dimensional Lebesgue measure; the L 2 -norm is denoted by

f 2 = X f 2 (x) dx 1/2 .
In this section, given a sample of bivariate observations, say (Y 1 , . . . , Y n ) t , we focus on the asymptotic properties of two RFA-madogram estimators. Two cases can be studied: when the marginal distributions, F 1 and F 2 , are known or when we need to use their empirical estimator, say F 1 and F 2 . In both cases, the copula function of the bivariate vector (Y 1 , Y 2 ) t , say C(u 1 , u 2 ), that captures the dependence structure needs to be inferred. To derive our asymptotic results, we adapt the main ingredients of Theorem 2.4 from [START_REF] Marcon | Multivariate nonparametric estimation of the pickands dependence function using bernstein polynomials[END_REF] to our settings, see Appendix B for details. With the notation

a c (u) = F 2 {cF ← 1 (u)} , we can write D(c) = 1 2 E |a c (U 1 ) -a ← c (U 2 )| ,
where the bivariate vector U = (U 1 , U 2 ) t follows the copula C(u). This leads us to the estimators

D n (c) = 1 n n i=1 D c (U i ) , with U i = (F 1 (Y 1,i ), F 2 (Y 2,i )) t and D c (U i ) = |a c (U 1,i ) -a ← c (U 2,i )| .
If F 1 and F 2 are unknown and are replaced by their empirical estimators, we have, with âc

(u) = F2 c F ← 1 (u) , D n (c) = 1 n n i=1 Dc U i , with U i = F 1 (Y 1,i ), F 2 (Y 2,i ) t and D c U i = âc U 1,i -â← c U 2,i .
In practice, D n (c) is directly computed from the expression The normalising coefficient is chosen as the optimal one, c * . The x and y-axis indicate the dependency coefficient α in the logistic dependence, see Eq.( 5) and the ratio of tail parameters i.e. the homogeneity of the two r.v. A ratio equal to one corresponds to the homogeneous case.

D n (c) = 1 n n i=1 F 2 (cY 1,i ) -F 1 (Y 2,i /c) .
A ratio equal to 10 can be illustrated by the realistic case of ξ 1 = 0.1, ξ 2 = 0.01.

Still, the definition of D n (c) with U i facilitates the derivation of theoretical results by leveraging existing properties of the empirical copula

C n (u) = 1 n n i=1 I(U i ≤ u) and by writing D n (c) = [0,1] 2 D c (U) dC n (u) .
In particular, the following classical smoothness condition on copula C is needed, see Example 5.3 in Segers (2012) for details.

Condition (S)

For every i ∈ {1, 2}, the partial derivative of C with respect to u i exists and is continuous on the set {u ∈ [0, 1] 2 : 0 < u i < 1}.

Proposition 3.1. Let (Y 1 , . . . , Y n ) t be n independent and identically distributed random vectors whose common distribution has continuous margins and a copula function C that satisfies condition (S).

Let D be a C-Brownian bridge, that is, a zero-mean Gaussian process on [0, 1] 2 with continuous sample paths and with covariance function given by

Cov(D(u), D(v)) = C(u ∧ v) -C(u) C(v), u, v ∈ [0, 1] 2 . ( 6
)
Here u∧v denotes the vector of componentwise minima. We define the Gaussian process

D on [0, 1] 2 by D(u) = D(u) - ∂C ∂u 1 D(u 1 , 1) - ∂C ∂u 2 D(1, u 2 )
Then we can write that a) We have

D n (c) -D(c) ∞ → 0 almost surely as n → ∞. Moreover, as n → ∞, √ n {D n (c) -D(c)} ; {1 + D(c)} 2 2 1 0 {D(a ← c (x), 1) -D(a ← c (x), a c (x))} dx + 1 0 {D(1, a c (x)) -D(a ← c (x), a c (x))} dx b) We have D n (c) -D(c) ∞ → 0 almost surely as n → ∞, and as n → ∞, √ n D n (c) -D(c) ; -{1 + D(c)} 2 1 0 D {a ← c (x), a c (x)} dx c>0 .

Analysis of CMIP precipitation for 16 models under two experiments

We now apply the RFA-madogram to the problem of partitioning annual precipitation maxima from 16 CMIP GCMs (see Table 1 in Appendix) into homogeneous regions.

For each hemisphere of a given GCM run, we estimate the dissimilarity matrix D(c * ) (Eq.( 3)) between each pair of grid points. To cluster from a dissimilarity matrix, the PAM clustering algorithm is implemented as it is fast, adapted to max-stable distributions [START_REF] Bernard | Clustering of maxima: Spatial dependencies among heavy rainfall in France[END_REF], and it does not require the triangle inequality [START_REF] Schubert | Fast and eager k-medoids clustering: O (k) runtime improvement of the PAM, CLARA, and CLARANS algorithms[END_REF]. The counterfactual and factual (2071-2100) runs are analysed separately and later compared to identify possible differences. With 16 partitions in four clusters for each 16 counterfactual (factual) hemispheric runs, GCM in-between-model error becomes an issue in terms of interpretation. We therefore summarise them in one "central" partitions, which we obtain in two steps. First, partitions for each counterfactual hemispheric runs are relabelled so as to minimise the pairwise difference between two partitions by taking each (alternatively) as reference score. As an example with five grid points, the partitions {1 1 1 2 2 3} and {3 3 3 1 1 2} are equal up to the permutation (1, 3, 2). Then, we compute the probability of each grid point to belong to each of the clusters, and associate the corresponding grid point to the cluster of highest probability. For instance, grid point B is assigned to cluster 1 for 6 models out of 16, to cluster 2 for 9 models and to cluster 3 for only one model.

In the so-called central partition, B is then assigned to cluster 2 with probability 9/16. Partitions for the factual experiment are relabelled in order to minimise the difference with the counterfactual central partition.

For example, Figure 4 shows the central partitions in four clusters by hemisphere. Intense colours correspond to points that belong to the same cluster in most, if not all, model partitions. Beginning with the counterfactual experiment, we first note that the clusters are very coherent spatially, in stark contrast to marginal-(ω) and dependencebased (F-madogram) partitions (Figure 2), even though no geographical covariates were used in the clustering. The Northern Hemisphere is dominated by two clusters (pink and yellow), with two others (blue and turquoise) with limited spatial extent. The distribution is more even in the Southern Hemisphere, and also more zonally symmetric. These partitions, driven both by homogeneity and dependence, are generally consistent with precipitation climatology. In the Northern Hemisphere, the pink cluster extends over the storm track regions of the North Atlantic and Pacific Oceans, and over the Inter-Tropical Convergence Zone (ITCZ) around 10 • N. The blue cluster covers the dry subtropics above the Sahara, Southwest Asia and southwest of North America. The turquoise cluster is located in the dry zone above the cold Pacific tongue, while the yellow cluster includes most regions with semi-arid and continental climates. Still, it also includes monsoon-dominated regions (e.g., India) and the dry Arctic. In the Southern Hemisphere, arid regions in Antarctica and in the dry descent regions at the eastern edge of the subtropical anticyclones are grouped together in the purple cluster, while the red cluster covers much of the wet tropics. The orange and green clusters correspond to the Southern Hemisphere storm track.

Most of the clusters appear to be quite robust across GCMs. Notable exceptions are the ITCZ regions in the Northern Hemisphere, and the equatorial Pacific and the eastern Indian Ocean west of Australia in the Southern Hemisphere. This lack of robustness may be due to the choice of cluster number. In any case, some differences are expected across GCMs, as they differ in their representation of storm tracks, monsoons or ITCZ location and dynamics. At first order, it appears that homogeneity of the distributions plays the dominant role, with arid or wet regions grouped together in both hemispheres. Still, the clustering is by design not only based on marginal distributions but also on dependence strength. To measure the importance of dependence in the spatial structure, we apply our clustering algorithm to temporally shuffled annual maxima at each grid point. This removes any spatial dependence between variables while preserving their marginal distributions. The results of Figure 5 for the CCSM4 model show a much less spatially coherent partition for the shuffled data. The dependence thus plays an important role in the coherence of the partition. This role can be further quantified by computing the relative difference between RFA-madogram on shuffled and non-shuffled data (with respect to the medoids). For about 2/3 of the grid points, the RFA-madogram takes lower values on the nonshuffled data, in particular near the medoids. We now turn to the comparison of the central partitions between the factual and counterfactual experiments. The overall partition structure is very similar in both experiments (Figure 4). The clusters are better defined in the counterfactual experiment (i.e. cluster probabilities closer to 1) because the sample size is much larger than for the factual experiment (155 versus 30 years). Globally, differences between the two central partitions are not significant compared to variability of model partitions compared to the central partition for either the factual or the counterfactual experiment (not shown). Hence, we cannot conclude to more spatial pattern variability in the factual world. The most likely cluster changes for a number of grid points, however, as indicated by crosses on the bottom panel of Figure 4. In the Northern Hemisphere, the pink (humid) and blue (arid) clusters expand slightly Northwards. More specifically, the probability of a given grid point to belong to the pink cluster generally increases at high latitudes, while the probability to belong to the blue cluster increases around the 25 • N latitude. In the Southern Hemisphere the green cluster (humid) also expands Southwards. While the resolution of our analysis is rather low (5 • ), these differences are consistent with the expected polewards shift of major climate zones under climate change, particulary the arid subtropics and the storm track regions of both hemispheres [START_REF] Scheff | Robust future precipitation declines in cmip5 largely reflect the poleward expansion of model subtropical dry zones[END_REF].

Conclusion

When considering multivariate data, extreme value theory can be difficult to handle. Reducing the dimensionality of extreme precipitation data set is then a challenging task. Our main goal in this work was to show that a simple and fast clustering approach based on an interpretable dissimilarity could highlight climatologically coherent regions.

The proposed approach coupled the main RFA idea, i.e. a normalising factor, with the dependence structure via the F-madogram. The introduced dissimilarity has links with extreme value theory via the extremal coefficient and tail parameters. The RFAmadogram neither requires estimating any marginal parameters nor dependence parameters. It is fully data-driven and bypasses the need of selecting relevant covariates or dependence structure.

Our analysis of annual maxima of daily precipitation from each CMIP model provides more spatially coherent hemispheric regions than some other non-parametric methods focusing on only one aspect (either homogeneity or dependence). Another contribution of this work is the handling of multi-partitions as our selected CMIP set has 16 GCM runs. Our combining approach enables us to compare one multi-model partition of the factual (all forcings) world with another multi-model partition of counterfactual (natural forcings) world. It appears that spatial variability between all models for the factual (resp. counterfactual) experiment appears to be significantly higher than between the two factual and counterfactual experiments.

In this work, we focus on the spatial structure of annual maxima precipitation in CMIP models, and on the forcing impact. We did not directly study the changes in rainfall distributions and frequencies. One interesting perspective would be to model precipitation intensities and dependence structure within each cluster. This could be useful for the D&A community. Another aspect is that the statistical approach developed therein is easy-to-implement and flexible, e.g. it can be used on non-gridded products.

For example, it could be applied to large weather networks, reanalysis (ERA 5) and radar products. Such datasets have finer spatial resolution scales than GCMs, and the dependence structure could be stronger, and consequently the analysis of heavy rainfall spatial patterns at fine spatial scales improved.

Proof of Lemma B.1: Note that

δ (u) = max (a(u 1 ), a ← (u 2 )) - 1 2 (a(u 1 ) + a ← (u 2 )) .
For any u ∈ [0, 1] 2 , we have max (a(u 1 ), a ← (u 2 )) = 1 -

1 0 I (u 1 ≤ a ← (u), u 2 ≤ a(u)) du and 1 2 (a(u 1 ) + a ← (u 2 )) = 1 - 1 2 1 0 I (u 1 ≤ a ← (u)) du + 1 0 I (u 2 ≤ a(u)) du.
Subtracting both expressions and integrating over H implies

[0,1] 2 δ (u) dH (u) = 1 2 [0,1] 2 1 0 I (u 1 ≤ a ← (u)) dudH(u 1 , u 2 ) + [0,1] 2 1 0 I (u 2 ≤ a(u)) dudH(u 1 , u2) - [0,1] 2 1 0 I (a ← (u 1 ) ≤ u, a(u 2 ) ≤ u) dudH(u 1 , u 2 ).
The stated lemma can be deduced by applying Fubini's theorem on the three double integrals. By Lemma B.1, we obtain for a

c (u) = F 2 (cF ← 1 (u)) D n (a c ) = (φ(C n ))(a c ) and D(a c ) = (φ(C))(a c ). this leads to ||D n (a c ) -D(a c )|| ∞ ≤ 2||C n -C|| ∞ .
Classical results about empirical copulas gives uniform strong consistency, see Segers .... Similar arguments can be used for D n (â c ). Now, we can consider the empirical process

D n = √ n(C n -C), D n = √ n( C n -C).
and we can write

√ n D n (a c ) -D(a c ) = (φ(D n ))(a c ) and √ n D n (â c ) -D(â c ) = (φ( D n ))(â c ).
We recall now that in the space ∞ ([0, 1] d ) equipped with the supremum norm, D n ; D, as n → ∞, where D is a C-Brownian bridge, and, as condition (S) holds, then D n ; D, as n → ∞, where D is the Gaussian process defined in (3.1), see [START_REF] Segers | Asymptotics of empirical copula processes under non-restrictive smoothness assumptions[END_REF] for details.

In addition, âc converges in probability to a c . The continuous mapping theorem then implies, as n

→ ∞, √ n D n (a c ) -D(a c ) = φ(D n ) ; φ(D), √ n D n (â c ) -D(â c ) = (φ( D n )) ; φ( D), in ∞ ([0, 1]).
From the continuity of its sample paths and by the form of the covariance function ( 6), the Gaussian process D satisfies

P{∀ u ∈ [0, 1] : D(u, 1) = D(1, u) = 0} = 1.
Ths provides all the elements to conclude the proposition. .

C. Expression of D(c) in the bivariate GEV case

As |a -b| = 2 max(a, b) -a -b, we have

2D(c) = 2E [max (F 2 (cY 1 ) , F 1 (Y 2 /c))] -E [F 2 (cY 1 )] -E [F 1 (Y 2 /c)]
To deal with each term, we recall that the quantile function of F (x; ξ, σ) = exp -x σ

-1/ξ is F -1 (u; σ, ξ) = σ (-log u) -ξ = σz ξ , with z = -1/ log(u),

This implies that

Y i d = σ i Z ξi i ,
where Z i follows an unit Fréchet. If follows that, with a 12 = cσ By noticing that W i follows a Weibull distribution with P(W 1 > w) = exp(-w -ξ2/ξ1 ), the expectation E[F 2 (cY 1 )] can be linked as the Laplace transform of a Weibull r.v. For the bivariate structure, we can write that, for any u ∈ (0, 1), P [max (F 2 (cY 1 ) , F 1 (Y 2 /c)) ≤ u] = P max exp -a 12 Z -ξ1/ξ2

where W i follows a Weibull distribution with P(W 1 > w) = exp(-w ξ1/ξ2 ). Note that (a 12 )

ξ 2 ξ 1 = 1 a 21
Conversely, (a 21 )

ξ 1 ξ 2 = 1 a 12

D. Homogeneous case

In the special case where ξ 1 = ξ 2 = ξ, we denote θ c = V (a 12 , a 21 ), where a 12 = cσ 1 σ 2 -1/ξ = 1/a 21 . Then, we have (a12,a21) .

P [max (F 2 (cY 1 ) , F 1 (Y 2 /c)) ≤ u] = exp V σ 2 cσ 1 -1/ξ , cσ 1 σ 2 -1/ξ log u , = u V
We can write

D = 1 0 1 -u θc du - 1 2 E [exp (-a 12 W 1 )] - 1 2 E [exp (-a 21 W 2 )] (9)
where W i , i = 1, 2 has cdf equal to exp(-x).

Hence,

D = θ c θ c + 1 - 1 2(1 + a 12 )
-1 2(1 + a 21 ) .

To minimise D as a function of c, we study the variations of r : x -→ V x,

1 x 1 + V x, 1 x - 1 2(1 + x) - x 2(1 + x)
. We suppose that V is differentiable. If the previous function r admits a minimum, its derivative cancels in some c 0 . The r cancels if and only if the derivative of x -→ V x, 1

x 1 + V x, 1

x cancels, if and only if there exists x s.t. ∂V ∂x x, 1 x = 1 x 2

∂V ∂y

x, 1

x . In the special case where the dependence is logistic i.e.

V (x, y) = 1 x 1/α + 1 y 1/α α , we have ∂V ∂x x, 1 x = ∂V ∂y x, 1

x , for all positive x. Therefore, if r admits a minimum, it is for x = ±1. Eventually, for logistic dependence, D is minimal for c = σ 2 σ 1 .

Fig. 1 .

 1 Fig. 1. Localisation (a), QQ-plots (b) and (d) and scatter plots (c) and (e) of annual precipitation maxima at three grid points A,B and C in the counterfactual run of the CCSM4 model (1850-2005). Panels (b) and (d) show the QQ-plots of rescaled precipitation for pairs (A,B) (b) and (A,C) (d). Panels (c) and (e) display the (rescaled) scatter plots for the same pairs.

  for all positive c. Therefore, c * 12 case of equality in distribution, Y 1 d = Y 2 , corresponds to the case where c * 12 = c * 21 = 1. An important feature of Eq.(

Fig. 2 .

 2 Fig. 2. Two summaries of the structure of the precipitation annual maxima of counterfactual (1850-2005) CCSM4 model. (a) Pointwise ω ratio (Le Gall et al., 2021). High values of ω correspond to heavy tailed distributions. (b) Results of PAM clustering with the F-madogram distance[START_REF] Bernard | Clustering of maxima: Spatial dependencies among heavy rainfall in France[END_REF], with four clusters for each hemisphere separately. Each color corresponds to a cluster.

Fig. 3 .

 3 Fig. 3. Distance (z-axis) D defined in Eq.(3) in the logistic bivariate GEV model example.The normalising coefficient is chosen as the optimal one, c * . The x and y-axis indicate the dependency coefficient α in the logistic dependence, see Eq.(5) and the ratio of tail parameters i.e. the homogeneity of the two r.v. A ratio equal to one corresponds to the homogeneous case. A ratio equal to 10 can be illustrated by the realistic case of ξ 1 = 0.1, ξ 2 = 0.01.

Fig. 4 .

 4 Fig. 4. Central partitions of CMIP models (with four clusters for each hemisphere), for (top) the counterfactual experiment (1850-2005) and (bottom) the factual experiment (2071-2100).Each colour corresponds to a cluster, with the shade indicating the probability of belonging to that cluster. In the bottom map, brown crosses indicate points where the most likely cluster is different between the counterfactual and the factual experiments.

Fig. 5 .

 5 Fig. 5. Partition of CCSM4 model in the counterfactual experiment based on the RFA-madogram dissimilarity D(c * ) and PAM algorithm, for (top) original data, and (bottom) data randomly shuffled in time at each grid point. The clustering algorithm is applied to each hemisphere independently.

E

  [F 2 (cY 1 )] = E [exp (-a 12 W 1 )] and E [F 1 (Y 2 /c)] = E [exp (-a 21 W 2 )] .

1

  , exp -a 21 Z Since the r.v. max (F 2 (cY 1 ) , F 1 (Y 2 /c)) ≤ u is positive, in the general setup, we have

											-ξ2/ξ1 2	≤ u ,
					= P Z 1 ≤	-a 12 log u	ξ2/ξ1	, Z 2 ≤	-a 21 log u	ξ1/ξ2	,
					= exp -V	-a 12 log u	ξ2/ξ1	,	-a 21 log u	ξ1/ξ2	.
	D =	0	1	1 -exp -V	a 12 -log u	ξ2/ξ1	,	a 21 -log u	ξ1/ξ2	du
		-	1 2	E [exp (-a 12 W 1 )] -	1 2	E [exp (-a 21 W 2 )] ,

Acknowledgement

Within the CDP-Trajectories framework, this work is supported by the French National Research Agency in the framework of the "Investissements d'avenir" program (ANR-15-IDEX-02).

Part of this work was supported by the DAMOCLES-COST-ACTION on compound events, the French national program (FRAISE-LEFE/INSU and 80 PRIME CNRS-INSU), and the European H2020 XAIDA (Grant agreement ID: 101003469). The authors also acknowledge the support of the French Agence Nationale de la Recherche (ANR) under reference ANR-20-CE40-0025-01 (T-REX project), and the ANR-Melody.

A. Proof of Eq. (4)

We can write that

In the same way, we can write that

It follows that the inequality expressed in Eq.( 4) is valid.

B. Proof of Proposition 3.1

Let a(u) be any continuous non-decreasing function from [0, 1] to [0, 1] and denote its inverse by a ← (u). The map

is linear and bounded, and therefore continuous. To continue, we need the following lemma.

Lemma B.1. For any cumulative distribution function H on [0, 1] 2 and for any nondecreasing function a(.) on [0, 1], the function

CNRM-CM5

Centre National de Recherches Météorologiques France 1.4

CSIRO and Bureau of Meteorology Australia 1.9


Mengfan Cai

Christophe Guy

Martin Héroux

Eric Lichtfouse

Chunjiang An
email: chunjiang.an@concordia.ca

The impact of successive COVID-19 lockdowns on people mobility, lockdown efficiency, and municipal solid waste

Keywords: COVID-19, Lockdown, Resurgence, Municipal solid waste, Waste amount and composition

The COVID-19 pandemic has induced many issues for all societal sectors, in particular in the production and disposal of municipal solid waste. This may be because successive easing and reimposing of lockdown measures have deeply changed people's movements, consumers' behaviors and waste management. Previous studies have focused on the shortterm effects of lockdowns on waste changes, yet there is little knowledge on waste variations during successive lockdowns and unlock-ing of various lockdown intensities. Moreover, the efficiency of lockdown and its relation to people's mobility in different countries are still not clear. Here, we studied the variations of amount and composition of municipal solid waste before the pandemic in 2019 and during the pandemic in 2020-2021 in USA, Brazil, Canada, UK, France and Italy. We used a stringency index and a composite mobility index to assess the lockdown intensity and people's movements. Results show that the mobility index sharply decreased with lockdown intensity, and enforcing measures were more efficient in France and Italy. Compared to 2019, prolonged lockdowns caused larger decreases in the quantity of commercial and construction wastes versus household waste. The initial implementation of lockdown or unlocking measures promoted inhabitants' con-sumption, generally leading to the increased waste amount, by about 9% for Trento and 12% for Montreal at the beginning of lockdown, respectively. Moreover, larger variations in the waste amount in Trento, from -25.3 to 9.8%, were in line with higher lockdown intensity compared to those in Montreal, from -9.5 to 12.7%, affected by people's mobility, consumers' behaviors and waste management.

Introduction

The coronavirus disease 2019 (COVID-19) reported 172 million confirmed cases and 3.7 million deaths by June 3, 2021, globally (Google News 2021). The most affected countries included the United States (USA), Brazil, India, Mexico, the United Kingdom (UK), Italy, Russian, and France (Google News 2021; [START_REF] Huo | Groundwater contamination with the threat of COVID-19: insights into CSR theory of Carroll's pyramid[END_REF]. Protective measures such as mask wearing, traffic restriction, social distancing, and home isolation have been globally implemented to reduce the spread of coronavirus [START_REF] Wang | Disposable masks release microplastics to the aqueous environment with exacerbation by natural weathering[END_REF]). Many regions have experienced lockdowns since the beginning of this pandemic in March-April 2020, followed by gradual unlocking since May 2020. After this first pandemic wave, increasing economic and social activities have induced a second pandemic wave, as expected, in Euro-American countries since July 2020, and, in turn, new restrictions and lockdowns were imposed since September 2020. Lockdown periods have generally had positive impacts on environmental quality, with a notable decrease in urban air pollution (Wang et al. 2020a;[START_REF] Zambrano-Monserrate | Indirect effects of COVID-19 on the environment[END_REF]), yet there are concerns about the effect of lockdown on waste generation and recycling [START_REF] Kulkarni | Repercussions of COVID-19 pandemic on municipal solid waste management: challenges and opportunities[END_REF]. Municipal solid waste (MSW) is an unavoidable by-product of human activity, and its recycling and recovery can facilitate the progression of a circular economy [START_REF] Das | COVID-19 and municipal solid waste (MSW) management: a review[END_REF]. However, the COVID-19 pandemic has brought MSW systems severe challenges, such as the quantity variation and composition change of MSW [START_REF] Vaverková | Municipal solid waste management under COVID-19: challenges and recommendations[END_REF][START_REF] Vu | Modeling of municipal waste disposal rates during COVID-19 using separated waste fraction models[END_REF][START_REF] Yousefi | Municipal solid waste management during COVID-19 pandemic: effects and repercussions[END_REF]. For instance, a significant escalation of medical waste was observed in some Asian cities and countries, such as India, during the lockdown compared to levels before the pandemic [START_REF] Somani | Indirect implications of COVID-19 towards sustainable environment: an investigation in Indian context[END_REF][START_REF] You | COVID-19's unsustainable waste management[END_REF][START_REF] Zand | Emerging challenges in urban waste management in Tehran, Iran during the COVID-19 pandemic[END_REF]. Nonetheless, the lockdown caused only a small waste rise in Singapore and even substantial reductions of MSW in India and Shanghai [START_REF] Fan | An update of COVID-19 influence on waste management[END_REF][START_REF] Somani | Indirect implications of COVID-19 towards sustainable environment: an investigation in Indian context[END_REF]. Therefore, the effect of lockdown on waste dynamics can vary among regions due to distinct geographical, economic, and sociological characteristics. In particular, MSW production may be influenced by lockdown intensity and people's movements. Most studies mentioned above investigated the short-term effect of single lockdown on waste variations. The efficiency of lockdown and its relationship to people's mobility in different countries are still not clear. Here, we studied the effect of lockdown measures on people mobility, and the influence of successive lockdowns on waste variations in Euro-American regions.

Experimental

We selected six countries in Euro-America, including USA, Brazil, Canada, UK, France and Italy, in consideration of their large number of confirmed cases and deaths. Information on the timing and stringency of lockdown measures across the world was collected from the Oxford Coronavirus Government Response Tracker (OxCGRT 2021) and the COVID-19 Community Mobility Reports (2021). We also collected documents that provided the MSW amount in 2020-2021 or the variations in MSW amount between 2019 and 2020-2021 from news and municipal information. The City of Montreal in Canada and the City of Trento in Italy were selected as representative cities for America and Europe, respectively, to analyze the influence of varying pandemic lockdowns on waste variations during different lockdown periods (Trento Municipality 2021; City of Montreal 2020).

Lockdown intensity and people mobility

During the pandemic, the government lockdown intensity and people's movements changed with varying lockdowns, resulting in variations in the amount, composition and management of MSW. Here, we adopted a stringency index (S I , dimensionless) and a composite mobility index (M j , dimensionless) to represent the government lockdown intensity and people's movements, respectively. The S I was comprehensively calculated by nine response indicators (OxCGRT, 2021): where I j is the sub-index score for any given indicator (j), such as closures of school, workplace, and public transport; restrictions on public events, public gatherings, internal movements, international travel, and public information campaigns; and stay-at-home requirements. The detailed equation for the S I can be found at [START_REF] Oxcgrt | Oxford Covid-19 Government Response Tracker[END_REF]. In this study, the lockdown intensity was evenly classified into five levels based on the S I , that is, level 1, 2, 3, 4 and 5 for S I of 0-20, 20-40, 40-60, 60-80 and 80-100, respectively. If lockdown measures vary at the subnational level, S I is shown as the response level of the strictest subregion. This may result in the overestimation of the lockdown intensity of some cities. Mobility trends can show how peoples' movements have changed throughout the pandemic; thus, they were adopted to describe the specific lockdown intensity at different spatial scales. According to the COVID-19 Community Mobility Reports (2021), the mobility change (M i j , %) in visitor numbers (N i ) to specific categories of location between a COVID-19 pandemic period (j) versus the baseline period before the pandemic was measured using Eq. 2:

where i indicates one of six broad categories of location, including residential places, grocery and pharmacy stores, workplaces, parks, transit stations, and retail and recreation. In this study, principal component analysis (PCA) was conducted using SPSS 16.0 to generate a composite mobility index for the comparison of different regions (Wang et al. 2020b). All the statistical analyses and figures were conducted and prepared using Microsoft Excel 2020, SPSS 16.0, and Origin Pro 9.1.

PCA results are shown in Table 1 and Supplementary Material's Tables S.1-6. Table S.1 shows that there is a linear correlation between six variables because the correlation coefficient between each group of variables is generally greater than 0.3. The overall KMO test coefficient is 0.707; the KMO test coefficient of each variable is greater than 0.5; and Bartlett's test result is p < 0.001 (Tables S.2-S.3). These results mean that the data structure is reasonable, and the data can be extracted by principal components. Table 1 and Table S.4 suggest that only the eigenvalue of the top one principal component in this study is greater than 1, explaining 69.318% of the total data variation. Tables S.5-S.6 indicate the high explanatory ability of this principal component to six variables. Thus, the extracted principal component, (1) S I = ∑ 9 j=1 I j 9

(2)

M j i (%) = N COVID j i -N baseline j i N baseline j i × 100%
the composite mobility index (M j , dimensionless), could be directly used as a newly generated variable for further data analysis.

Amount of municipal solid waste

In this study, organic and recyclable wastes were the major compositions, accounting for around 32% and 36% of the total waste amount (tons) in Trento and 40% and 29% in Montreal in 2019, respectively (Table 2). The recyclable waste included paper and cardboard, glass, metals and plastics in Montreal. Besides these waste types, lightweight packaging was additionally classified as recyclable waste in Trento. For the two cities, the amount of organic waste was the sum of food and green wastes. The impacts of the COVID-19 lockdown in changing waste amount and composition were evaluated by waste variations. The variations (V i j , %) of any type of MSW (W i , tons) between a COVID-19 pandemic period (j) in 2020-2021 versus the same period in 2019 were calculated using Eq. 3:

Results and discussion

Effects of lockdown intensity on people mobility

We studied the effect of successive lockdowns and unlocking on the amount of municipal waste during the pandemic in 2020-2021, compared to 2019, in six countries. The (3)

V j i (%) = W COVID j i -W 2019 j i W 2019 j i × 100%
mobility index represented the mobility change in visitor numbers throughout the pandemic, and the rapid variations in mobility were generally seen on weekends or holidays. Figure 1 shows that these countries have generally experienced four lockdown periods in 2020: pre-lockdown in January-February, first lockdown in February-May, unlocking in May-October, and second lockdown in November-December. However, Italy, France and USA implemented unlocking measures in January 2021, while the other three countries imposed stricter restrictions. Italy and France implemented the lockdown measures about half a month earlier than Canada and Brazil (Fig. 1).

Results show that the lockdown intensity increased sharply at the onset of the first lockdown for all countries (Fig. 1). During the first lockdown, the intensity was slightly higher, at level 5, for France and Italy, whereas the intensity of other countries was at level 4. Then, during the first unlocking, the intensity decreased more for France and Italy, to level 3, than for other countries, with few or no decrease at level 4 until 2021. The second lockdown induced again a sharp increase of lockdown intensity for France and Italy. Overall, France and Italy experienced stronger variations in lockdown intensity.

The people mobility displayed in general reversed trends versus lockdown intensity, with more amplitude for France and Italy (Fig. 1). Mobility was slightly lower in major cities such as Trento and Montreal, versus the whole countries, which possibly indicates stricter enforcement in cities versus rural areas. Pearson correlation of stringency index and people mobility shows stronger relationships (r) for France (-0.817) and Italy (-0.803) than for UK (-0.758), USA (-0.698), Brazil (-0.600) and Canada (-0.597) (Table S.7). Therefore, the lockdown was more efficient in France and Italy, which both applied a 'rest' period during the first unlocking. In other words, a too long period of intense locking is less effective, as clearly observed for Brazil, and to a lesser extent in USA and Canada where people mobility increased steadily despite high lockdown intensity. A possible explanation is that during prolonged lockdown, people get stressed by immobility, adapt and find a way to move. Overall, higher variations in lockdown intensity in France and Italy were more efficient at reducing the people mobility, versus prolonged lockdowns in UK, USA, Brazil and Canada.

Waste amount dynamic during successive lockdowns

Figure 2 shows the waste variations (%) in the amount and composition of Montreal and Trento. Compared to the normal periods in 2019, significant decreases in total waste amount were observed in most of the months in 2020-2021 in Montreal (e.g., by 9.5% in May 2020) and Trento (e.g., by 13.7%, 25.3%, 14.7% from March to May 2020 and 16.5% in January 2021). For the major waste composition, both organic and recyclable wastes showed similar trends as the total waste amount, down by 8-19%. These results were observed when the national lockdown intensity is high, at levels 4-5, and the municipal mobility index is low, ranging from -2 to -1 (Fig. 1). Moreover, waste decreases in Trento were generally larger than those in Montreal because of the higher lockdown intensity and lower mobility index in Trento. This indicates that the significant reduction of people's movements could cause large decreases in waste amount. Notably, New York City and Stanford University communities in USA experienced significant decreases in commercial and construction wastes, nearly 50% and 60%, respectively, in March (The Conversation, 2020; WasteAdvantage, 2020). Brno's household waste also decreased by 2.6% and 6.9% in April and May, respectively, and its commercial and industrial wastes dropped nearly 40% in June [START_REF] Fan | An update of COVID-19 influence on waste management[END_REF]). UK's residual waste decreased by about 11%, with commercial and industrial waste down by about 40% during April-June [START_REF] Tolvik | Residual Waste during Q2 2020 lockdown[END_REF]. These results indicate that prolonged lockdowns caused larger decreases in commercial and construction wastes than household waste due to the drastic reduction in business, construction and tourist activities. There were about 12% increases in the total waste amount in both March and June for Montreal, while nearly 9% increases were all observed in February, June and September for Trento. Each type of Montreal's waste increased by about 11% in March and June. However, Trento's organic waste increased by 10.6% in February and its packaging waste increased by 12.1% and 16.4% in June and September, respectively (Fig. 2). Besides, February-April 2020 and November 2020-January 2021 experienced the first down and second lockdown periods, respectively. Thus, waste variations in Trento showed periodic changes under re-lockdown measures as that observed since the first lockdown period (Figs. 1 and2). In the first month of lockdowns, the panic shopping of necessary items, inappropriate cooking habits, such as overcooking, broken food supply chains, packaged take-out orders and home-delivered businesses, might result in increased food and packaging wastes [START_REF] Aldaco | Food waste management during the COVID-19 outbreak: a holistic climate, economic and nutritional approach[END_REF][START_REF] Jribi | COVID-19 virus outbreak lockdown: what impacts on household food wastage?[END_REF]). In addition, unlocking might cause many inhabitants to return to normal food habits and produce more food waste than the average (WRAP 2020). Thus, revenge consumption might explain June's increase in Montreal's MSW. Besides, the transition from in-store shopping to online shopping would increase the amount of in-home recyclable packaging waste [START_REF] Ikiz | Impact of COVID-19 on household waste flows, diversion and reuse: the case of multiresidential buildings in Toronto Canada[END_REF]. The visitors to retail & recreation places and grocery & pharmacy stores were still much fewer than those during the normal periods, even though the lockdown intensity was reduced to level 3 (Fig. 1; COVID-19 Community Mobility Reports 2021). This means that lockdown measures influenced consumers' behaviors during the pandemic, leading to an increase in online shopping and packaging wastes in Trento. These increases were observed in the first month of lockdowns and the second month of unlocking, likely because unlocking requires more response time. Overall, the initial implementation of lockdown or unlocking measures promoted inhabitants' consumption, generally leading to an increase in waste amount. Moreover, consumers' behaviors in different lockdown periods affect variations in waste amount and composition. Besides people's movements and consumers' behaviors, waste management affects waste amount and composition. For example, MSW was generally collected more on Mondays, Tuesdays and Thursdays than on Wednesdays and Fridays in Montreal. Therefore, the significant monthly variations in Montreal's MSW amount were likely due to more "major collection weekdays" in March and June while less in April, May and July in 2020 compared to 2019. Considering the spreading risks of contaminated wastes in recycling services, some Euro-American areas, including the six selected countries, restricted selective collections, e.g., a 35.5% reduction between March 23 and April 24 in Brazil, and recycling services, e.g., a 46% decrease in USA. Alternately, single-use plastic packaging was resurged, and recovery facilities were partially stopped in these countries. Consequently, the amount of recyclable waste generally decreased, but mixed and bulky wastes increased during March-April 2020 because of significant changes in waste management services during the early pandemic [START_REF] Fan | An update of COVID-19 influence on waste management[END_REF][START_REF] Gorrasi | Back to plastic pollution in COVID times[END_REF][START_REF] Ragazzi | Municipal solid waste management during the SARS-COV-2 outbreak and lockdown ease: Lessons from Italy[END_REF]. Overall, larger variations in the waste amount in Trento were in line with higher lockdown intensity compared to those in Montreal, affected by people mobility, consumers' behaviors and waste management.

Conclusion

Our results showed the long-term effect of varying pandemic lockdowns on waste variations during different lockdown periods in some Euro-American regions during 2020-2021. Overall, the government lockdown intensity and people's movements changed with the varying lockdowns, resulting in variations in the amount and composition of MSW during the pandemic. Generally, there have been mainly four successive lockdown periods in 2020 and people mobility appears better controlled in Italy and France versus USA, Brazil, Canada, UK. We found prolonged lockdowns caused larger decreases in commercial and construction wastes than household waste due to the drastic reduction in business, construction and tourist activities. The beginning of lockdown or unlocking measures promoted inhabitants' consumption, generally leading to an increase in waste amount. Moreover, larger variations in the waste amount in Trento were in line with higher lockdown intensity compared to those in Montreal, affected by people mobility, consumers' behaviors and waste management. Notably, the inadequacies of waste management systems in the early pandemic may lead to a decrease in waste collection of recyclable waste but an increase of bulky and mixed wastes. Thus, an appropriate strategy for the whole waste management system, including waste generation, collection, transport, recycling and disposal, should be considered to reduce waste for incineration and direct landfilling and increase waste recycling. Furthermore, the use of some new approached such as the automated waste management system through the Internet of Things can alleviate the risk of labor shortage and workers' infection. Lastly, the resurgence of single-use plastic packaging is a temporary alternative, although it may cause further microplastic pollution and affect public engagement in MSW management after the pandemic. The recovery of waste management practices and policies based on a circular economy can be a priority.

Fig. 1

 1 Fig. 1 Lockdown intensity and mobility index in 2020 in different cities/countries in Euro-America during the pandemic. The red lines show the lockdown intensities (I) on the primary vertical axis which are evenly classified into five levels based on the stringency index, that is, level 1, 2, 3, 4 and 5 for 0-20, 20-40, 40-60, 60-80 and 80-100, respectively. The green lines show the mobility indexes (M j ) on the secondary vertical axis based on the national mobility changes

Fig. 2 ×

 2 Fig.2Increases (red) and decreases (green) of municipal waste amount and composition in Montreal and Trento during the pandemic. The variations (V i j , %) of any type of waste (W i , tons) between a COVID-19 pandemic period (j) in 2020-2021 versus the same period in 2019 were calculated using Eq. 3:V j i (%) =

Table 1

 1 Total variance explained

			Component	Initial eigenvalues		Extraction sums of squared loadings
					Total	% of Variance	Cumulative %	Total	% of Variance	Cumulative %
			1		4.159	69.318	69.318	4.159	69.318	69.318
			2		0.981	16.355	85.673
			3		0.528	8.808	94.481
			4		0.201	3.346	97.827
			5		0.088	1.470	99.297
			6		0.042	0.703	100.000
			Extraction method: principal component analysis
	Table 2 MSW amount (tons) in Montreal and Trento in 2019
	City	Waste amount (tons)
		Organic	Recyclable	Others	Total
	Montreal	372, 168	266, 855	281, 865	920, 888
	Trento	17, 695	19, 546	17, 480	54, 721

Funding This research was supported by the Natural Sciences and Engineering Research Council of Canada.

Supplementary Information

The online version contains supplementary material available at https:// doi. org/ 10. 1007/ s10311-021-01290-z.

Authors' contribution Mengfan Cai conducted the study and wrote the article. Christophe Guy, Martin Héroux, Eric Lichtfouse, and Chunjiang An conceived the study and revised the article.

Declarations

Conflict of interest This is to certify that the authors of this manuscript have no affiliations with or involvement in any organization or entity with any financial interest, or non-financial interest in the subject matter or materials discussed in this manuscript.

SUPPLEMENTARY MATERIAL

The impact of successive COVID-19 lockdowns on people mobility, lockdown efficiency, and municipal solid waste Principal component analysis (PCA) was adopted to generate the composite mobility indexes from six variables of locations, including residential places, grocery & pharmacy stores, workplaces, parks, transit stations, and retail & recreation. The results are shown in Table S.1-S.6.