Yiming Liu

Tao Wang

Trissevgeni Stavrakou

Nellie Elguindi

Thierno Doumbia

Claire Granier

Idir Bouarar

Benjamin Gaubert

Guy P Brasseur

Diverse response of surface ozone to COVID-19 lockdown in China

Keywords: Surface ozone, meteorological condition, emission reduction, COVID-19

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The outbreak of coronavirus disease 2019 (COVID-19) has severely threatened public health worldwide, leading to millions of deaths [START_REF] Who | Evolution and Assessment of the Atmospheric Composition in Hangzhou and its Surrounding Areas during the G20 Summit[END_REF]. China, where the first case of COVID-19 was reported in the city of Wuhan, imposed country-wide measures from 23 January to 13 February 2020 to prevent the spread of the disease, including social distancing, teleworking, and closure of non-essential businesses [START_REF] Chinazzi | The effect of travel restrictions on the spread of the 2019 novel coronavirus (COVID-19) outbreak[END_REF][START_REF] Li | Substantial undocumented infection facilitates the rapid dissemination of novel coronavirus (SARS-CoV-2)[END_REF]. These restrictions drastically reduced anthropogenic activities, resulting in a sharp decrease in emissions of air pollutants [START_REF] Doumbia | Changes in global air pollutant emissions during the Covid-19 pandemic: a dataset for atmospheric chemistry modeling[END_REF][START_REF] Huang | Enhanced secondary pollution offset reduction of primary emissions during COVID-19 lockdown in China[END_REF]Wang et al., 2020a).

The huge and large-scale emission reductions during the COVID-19 lockdown can be treated as a natural outdoor experiment to improve our understanding of the air pollutant's response to emission control. According to satellite and surface observations, compared with the period before the lockdown, nitrogen dioxide (NO2) concentrations decreased by over 50% in China during the lockdown period [START_REF] Bauwens | Impact of Coronavirus Outbreak on NO2 Pollution Assessed Using TROPOMI and OMI Observations[END_REF]Liu et al., 2020;[START_REF] Shi | The Response in Air Quality to the Reduction of Chinese Economic Activities During the COVID-19 Outbreak[END_REF][START_REF] Zhang | NOx Emission Reduction and Recovery during COVID-19 in East China[END_REF]. The concentrations of other pollutants, including SO2, particulate matter with an aerodynamic diameter less than 2.5 μm (PM2.5), particulate matter with an aerodynamic diameter less than 10 μm (PM10), and carbon monoxide (CO), also declined in a large area of China [START_REF] Miyazaki | Air Quality Response in China Linked to the 2019 Novel Coronavirus (COVID-19) Lockdown[END_REF][START_REF] Liu | Worsening urban ozone pollution in China from 2013 to 2017 -Part 2: The effects of emission changes and implications for multi-pollutant control[END_REF]. However, surface ozone (O3) concentrations in northern and central China increased by over 100% [START_REF] Lian | Impact of city lockdown on the air quality of COVID-19hit of Wuhan city[END_REF][START_REF] Shi | The Response in Air Quality to the Reduction of Chinese Economic Activities During the COVID-19 Outbreak[END_REF]. Similar O3 increases have been reported in southern Europe, India, and Brazil despite the large decrease in other pollutants [START_REF] Sharma | Effect of restricted emissions during COVID-19 on air quality in India[END_REF][START_REF] Sicard | Amplified ozone pollution in cities during the COVID-19 lockdown[END_REF][START_REF] Siciliano | Increased ozone levels during the COVID-19 lockdown: Analysis for the city of Rio de Janeiro, Brazil[END_REF].

However, the underlying factors driving the O3 changes during the city lockdowns remain unclear.

Surface O3 is produced by photochemical reactions of ozone precursors, NOx, volatile organic compounds (VOCs), and carbon monoxide (CO) and can also be transported from higher levels of the atmosphere and from outside regions [START_REF] Akimoto | Analysis of monitoring data of ground-level ozone in Japan for long-term trend during 1990-2010: Causes of temporal and spatial variation[END_REF]Liu and Wang, 2020a;[START_REF] Roelofs | Model study of the influence of cross-tropopause O3 transports on tropospheric O3 levels[END_REF]. It is well known that O3 has a non-linear dependence on its precursors and that NOx can either decrease or increase O3 depending on the relative abundance of NOx to VOCs [START_REF] Atkinson | Atmospheric chemistry of VOCs and NOx[END_REF]Wang et al., 2017a). In general, the O3 production in urban areas with high NOx/VOCs ratios is VOCs limited, and reducing NOx emissions can increase O3 due to decreased titration of O3 and radicals. In addition to the two precursors, particulate matters can influence ozone by altering the solar irradiance and chemical reactions on aerosol surfaces (Li et al., 2019b;[START_REF] Liu | Worsening urban ozone pollution in China from 2013 to 2017 -Part 2: The effects of emission changes and implications for multi-pollutant control[END_REF][START_REF] Stadtler | Ozone impacts of gas-aerosol uptake in global chemistry transport models[END_REF]. Meteorological factors affect surface ozone by changing transport pattern, wet and dry depositions, chemical reaction rates, and natural emissions (Liu and Wang, 2020a;[START_REF] Lu | Exploring 2016-2017 surface ozone pollution over China: source contributions and meteorological influences[END_REF].

The responses of ozone (and other air pollutants) to short-term emission reductions have been previously studied for a number of public and political events in China, such as the Beijing Summer Olympic Games (August 2008), the Asia-Pacific Economic Cooperation (APEC) meeting in Beijing (November 2014), and the G20 summit in Hangzhou (September 2016). During these events, various emission-reducing measures were implemented in the cities concerned and their surrounding areas. Whereas atmospheric concentrations of primary air pollutants (NOx, CO, primary PM, and SO2) in the concerned cities generally decreased in response to the temporary control measures, the O3 concentrations showed mixed responses. O3 decreased after emission reductions for some events [START_REF] Huang | Impacts of emission reduction and meteorological conditions on air quality improvement during the 2014 Youth Olympic Games in Nanjing, China[END_REF]Wang et al., 2017b) but increased in others [START_REF] Wang | Quantifying the Air Pollutants Emission Reduction during the 2008 Olympic Games in Beijing[END_REF][START_REF] Wang | Changes in atmospheric composition during the 2014 APEC conference in Beijing[END_REF]Wu et al., 2019). The different O3 responses have been qualitatively attributed to differences in the meteorological conditions (including regional transport of air masses) and to different control measures implemented by the local governments.

Compared with the previously studied situations, the COVID-19 lockdown is unique in that emissions decreased across the whole country (and later worldwide) as opposed to a specific city or region, and the decreases were also much more drastic than those due to transportation restrictions alone. Moreover, the COVID-19 lockdown took place in winter, whereas the previous interventions occurred in summer and autumn, when meteorology and atmospheric chemistry are different from winter. The present study analyzes surface O3 data across China before and during the COVID-19 lockdown. We find that O3 decreased in southern China while increasing in most other regions during the lockdown. Using a regional chemistry transport model, we isolate the impacts of meteorological changes and anthropogenic emission reductions on O3. Our results highlight the importance of meteorological influences on the short-term O3 changes and the diverse response of O3 to the emission reductions of its precursors in different climate and emission-mix regions.

Materials and Methods

Surface measurement data

We obtained the observed concentrations of surface O3 and other pollutants (PM2.5, PM10, SO2, CO, NO2) at 1643 stations from the China National Environmental Monitoring Center (http://106. 37.208.233:20035/). Data quality control was conducted for the measurement data in accordance with previous studies [START_REF] Lu | Severe Surface Ozone Pollution in China: A Global Perspective[END_REF][START_REF] Song | Air pollution in China: Status and spatiotemporal variations[END_REF]. Fig. 1 shows the locations of these environmental monitoring sites.

The country-wide measures to control the spread of COVID-19 were implemented starting from 23 January 2020 (the exact date varies for different cities), just before the Chinese New Year. All enterprises remained closed until no earlier than 13 February, except those required for essential public services, epidemic prevention and control, and residential life needs. We focused on the period during the COVID-19 lockdown from 23 January to 12 February 2020 (hereafter referred to as the CLD period), 3 weeks in total. We derived the changes in O3 and other pollutants by comparing the CLD period with the 3 weeks before the COVID-19 outbreak, from 2 to 22 January 2020 (hereafter referred to as the pre-CLD period). We focused on three typical regions in China (Fig. 1): north China (NC,(113)(114)(115)(116)(117)(118)(119)including Beijing,Tianjin,Hebei,and western Shandong), central China (CC, 28.8-33°N, 108-117°E, including Hubei province, where Wuhan is situated, and the surrounding regions), and south China (SC, 21.5-24°N, 111-116°E, including the Pearl River Delta and the surrounding regions). The NC region is situated in the North China Plain, which is known to suffer from severe haze in winter; CC was the original epicenter of the COVID-19 outbreak in China and is an important economic hub for the central regions of China; the Pearl River Delta, where the megacities of Guangzhou and Shenzhen are situated, is the most developed region in southern China.

Numerical modeling

The CMAQ model (Community Multiscale Air Quality model, v5.2.1) was applied to simulate the O3 mixing ratios over China from 2 January to 12 February 2020. The WRF model (Weather Research and Forecasting model, v3.5.1) was driven by the dataset of the National Center for Environmental Prediction (NCEP) FNL Operational Model Global Tropospheric Analyses with a horizontal resolution of 1° × 1° and provided meteorological inputs for the CMAQ model. The CMAQ target domain covered the continental China at a horizontal resolution of 36 km × 36 km. SAPRC07TIC [START_REF] Carter | Development of the SAPRC-07 chemical mechanism[END_REF][START_REF] Hutzell | Interpreting predictions from the SAPRC07 mechanism based on regional and continental simulations[END_REF] and AERO6i [START_REF] Murphy | Semivolatile POA and parameterized total combustion SOA in CMAQv5.2: impacts on source strength and partitioning[END_REF][START_REF] Pye | On the implications of aerosol liquid water and phase separation for organic aerosol mass[END_REF] were adopted as the gas-phase chemical mechanism and aerosol mechanism, respectively. The CMAQ model has been improved with updated heterogeneous reactions to better predict the O3 concentration; details can be found in Liu and Wang (2020a). Although the WRF-CMAQ model was run in offline mode, the CMAQ model employs an in-line method that uses the concentrations of particles and O3 predicted within a simulation to calculate the solar radiation and photolysis rates [START_REF] Binkowski | Examining Photolysis Rates with a Prototype Online Photolysis Module in CMAQ[END_REF]. As a result, the effect of aerosol on O3 concentrations via changing the photolysis rates were also considered in the simulation. The chemical boundary conditions were provided by the results of Whole Atmosphere Community Climate Model (WACCM, https://www.acom.ucar.edu/waccm/download). The anthropogenic emissions in China were obtained from the Multi-resolution Emission Inventory for China (MEIC) in 2017 (http://www.meicmodel.org) with scaling factors to the year 2020 (Table S1, see text in Supplementary Information), which were estimated based on the Three-Year Action Plan (2018-2020) issued by the government and the changes in the multipollutant emissions of different sectors in recent years [START_REF] Zheng | Trends in China's anthropogenic emissions since 2010 as the consequence of clean air actions[END_REF]. The emission adjustments during the lockdown period are based on recent publications (see text in Supplementary Information). Emissions from the other countries were derived from the MIX emission inventory [START_REF] Li | MIX: a mosaic Asian anthropogenic emission inventory under the international collaboration framework of the MICS-Asia and HTAP[END_REF]. International shipping emissions were taken from the Hemispheric Transport Atmospheric Pollution (HTAP) emission version 2.2 dataset for 2010 [START_REF] Janssens-Maenhout | HTAP_v2.2: a mosaic of regional and global emission grid maps for 2008 and 2010 to study hemispheric transport of air pollution[END_REF]. Biogenic emissions were calculated by the Model of Emissions of Gas and Aerosols from Nature (MEGAN) version 2.1 [START_REF] Guenther | The Model of Emissions of Gases and Aerosols from Nature version 2.1 (MEGAN2.1): an extended and updated framework for modeling biogenic emissions[END_REF] with meteorological inputs from the WRF model.

Two experiments were conducted to investigate the impacts of meteorological changes and emission reductions on O3 during the CLD period. The first (baseline) used the same anthropogenic emissions for the pre-CLD and CLD periods, and the second (Reduction case) used emission reductions of 70%, 40% and 30% for transportation, industry and power generation, respectively, and a 10% increase of residential emission during the CLD period. These emission reductions for the whole country were estimated according to the previous literature [START_REF] Doumbia | Changes in global air pollutant emissions during the Covid-19 pandemic: a dataset for atmospheric chemistry modeling[END_REF][START_REF] Huang | Enhanced secondary pollution offset reduction of primary emissions during COVID-19 lockdown in China[END_REF]Wang et al., 2020a).

Comparing these two model simulations, the O3 changes during the CLD period relative to the pre-CLD period for the Reduction Case were considered to be entirely due to the meteorological changes and emission reductions. The impacts of the meteorological changes (including the changes in chemical boundary conditions) were quantified by subtracting the O3 mixing ratios of the pre-CLD period from those of the CLD period for the baseline experiment, while the impacts of emission reduction were estimated by comparing the O3 mixing ratio during the CLD period between the Reduction Case and the baseline experiment. Furthermore, we individually reduced the emissions of nitrogen oxide (NOx), VOCs, CO, PM (particulate matter, including PM10, PM2.5, black carbon, and organic carbon), and SO2 during the CLD period to elucidate the response of O3 to each pollutant reduction.

The performance of the CMAQ model in simulating the O3, NO2, PM2.5, SO2, and CO concentrations for the Reduction Case was evaluated (Fig. S1 and Table S2), showing reasonable agreements with the respective surface observations. Details of the emission estimation and the model evaluation are presented in Supplementary Information.

Results

Observed O3 changes in different parts of China

Figures 2 and 3 present the changes in observed concentrations of O3 and other pollutants during the CLD period compared with pre-CLD. The concentrations of most pollutants (SO2, CO, PM2.5, PM10) that partially or fully originated from the direct emissions declined in China during the lockdown. NO2, a precursor of O3, decreased by about 50% across the entire continental China, and by similar amounts in all regions (Fig. 3b). However, the O3 mixing ratio exhibited varying changes in different regions (Fig. 2a). In NC and CC, the daily average O3 increased significantly, by 112% and 73%, respectively (Fig. 2b); in contrast, it remained almost unchanged in SC. The O3 changes also varied between daytime (8:00-20:00 LST) and nighttime (20:00-8:00 LST). During daytime, the O3 increase in most parts of China was smaller than the daily average (Fig. 2c), 92% and 71% in the NC and CC regions (Fig. 2d). In the SC region, most stations displayed a decrease in O3 during daytime, leading to a regional average O3 drop of 14%. During nighttime, the O3 mixing ratio increased significantly across China (Fig. 2e), by 154%, 77%, and 18% in NC, CC, and SC, respectively (Fig. 2f). These results reveal the diverse response of O3 during the lockdown in different regions, especially for the daytime. The changes in the Ox (NO2+O3) concentration (Fig. 3a), which takes into account the NO titration, also varied in different regions. The daytime average Ox levels increased by 4% in NC and by 11% in CC, and decreased by 29% in SC. These results suggest that the NO titration effect was not the only cause of the O3 increase in northern and central China, as Ox would have decreased with sharply reduced NOx emissions.

Contribution of meteorological changes and emission reductions to O3

Ground-level O3 is influenced by both chemical reactions of O3 precursors and meteorology. In this study, we used the WRF-CMAQ model to separate the impacts of meteorological changes and emission reductions on the changes in O3 across China (Fig. S2), which reveals significant contributions of both meteoroglgy (over most of contrinental China) and emissions (mainly in populated areas of eastern China) . Fig. 4 shows the more detailed results for the NC, CC, and SC regions for both daytime and nighttime and Fig. S3 presents the meteorological impact and emission impact in terms of percent change. The observed O3 changes in these regions were reasonably captured by the simulations. For the daytime average, the O3 increase in NC was attributed to the comparable contributions from both meteorological changes (58%) and emission reductions (42%) (Fig. 4a).

In CC (Fig. 4b), the meteorological change (98%) was the primary cause of the O3 increase, whereas the contribution of emission reduction was much lower (2%). In SC (Fig. 4c), the meteorological changes (73%) and emission reductions (27%) both contributed to the O3 decrease. During nighttime, the emission reduction increased O3 in all three regions (including SC), and its impact was stronger; the effect of meteorological changes weakened at night (Fig. S3).

Impacts of meteorological changes on O3

The impacts of meteorological changes on O3 for the NC, CC, and SC regions can be explained by the changes in the weather pattern and specific meteorological factors. In winter, continental China is generally controlled by a cold high-pressure system (Fig. 5). During our study period, the center of this high-pressure system was located in northern China, moving southward from the pre-CLD to the CLD period, with weakening strength. The high-pressure system therefore became increasingly dominant in southern China, and the strengthened southward winds brought colder air masses from the north (Fig. 6c), which decreased the temperature locally (Fig. 6a). In contrast, in central and northern China, the winds shifted to a more northward direction, transporting warmer air masses from the south (Fig. 6c), which increased the temperature (Fig. 6a). During daytime, the decrease (increase) in temperature in the SC region (CC and NC regions) weakened (enhanced) the surface O3 chemical production. Biogenic emission is an important source of VOCs and thereby contributes to O3 formation in China [START_REF] Wu | Estimation of biogenic VOC emissions and their corresponding impact on ozone and secondary organic aerosol formation in China[END_REF]. The temperature changes led to an increase (decrease) of biogenic emissions in the CC (SC) region (Fig. S4). Thus, the temperature changes increased (decreased) O3 in the CC (SC) region by influencing chemical reaction rates directly [START_REF] Fu | Positive but variable sensitivity of August surface ozone to large-scale warming in the southeast United States[END_REF][START_REF] Steiner | Observed suppression of ozone formation at extremely high temperatures due to chemical and biophysical feedbacks[END_REF] and altering biogenic emissions indirectly [START_REF] Im | The impact of temperature changes on summer time ozone and its precursors in the Eastern Mediterranean[END_REF]Liu and Wang, 2020a).

The changes in the weather pattern also resulted in less clouds and precipitation in northern and central China, but more clouds and precipitation in southern China (Fig. 6e andf). Clouds can reduce the amount of solar radiation reaching the surface and thus the chemical production of O3 [START_REF] Lelieveld | Influences of Cloud Photochemical Processes on Tropospheric Ozone[END_REF], while precipitation can also reduce O3 through the scavenging of its precursors [START_REF] Seinfeld | Atmospheric Chemistry and Physics-from Air Pollution to Climate Change[END_REF][START_REF] Shan | Observational study of surface ozone at an urban site in East China[END_REF]. The cloud and precipitation patterns therefore contributed to O3 increases in CC and NC and decreases in SC. Furthermore, in NC and CC, the significant increase in the planetary boundary layer height during the lockdown (Fig. 6d) might promote the transport of O3 from the upper air to the surface, contributing to the O3 increase in these regions [START_REF] He | Air pollution characteristics and their relation to meteorological conditions during 2014-2015 in major Chinese cities[END_REF][START_REF] Sun | Vertical observations and analysis of PM2.5, O3, and NOx at Beijing and Tianjin from towers during summer and Autumn 2006[END_REF]. The increase (decrease) in specific humidity in NC and CC (SC) might also have contributed to the decrease (increase) in O3 mixing ratios in those regions (Li et al., 2019c;[START_REF] Ma | Air pollution characteristics and their relationship with emissions and meteorology in the Yangtze River Delta region during 2014-2016[END_REF] (Fig. 6b). During nighttime, the changes in meteorological factors were similar to those in daytime (Fig. S5) but exerted smaller impacts on O3 changes due to the decreasing effects of temperature and cloud cover (negligible biogenic emissions and solar radiation).

Response of O3 to emission reductions

We further investigated the impact of multi-pollutant reductions on the O3 changes. Because transportation and industrial activities were reduced significantly during the lockdown and they are the major sources of NOx (>80%) and VOCs (>60%) (Fig. 7), the estimated reductions of NOx and VOC emissions were more significant than those for CO, particulate matter (PM), and SO2 (Fig. 8a,c,e). The NOx emission reductions were 46%, 49%, and 55% in the NC, CC, and SC regions, respectively, while the respective reductions for VOC emissions were 32%, 24%, and 45%. The relationship between O3 and the emissions of its precursors is non-linear. We used the ratio of production rates between H2O2 and HNO3 (PH2O2/PHNO3) [START_REF] Gaubert | Global Changes in Secondary Atmospheric Pollutants During the 2020 COVID-19 Pandemic[END_REF][START_REF] Tonnesen | Analysis of radical propagation efficiency to assess ozone sensitivity to hydrocarbons and NOx : 1. Local indicators of instantaneous odd oxygen production sensitivity[END_REF] to identify the O3 formation regime in China for the periods before and during the lockdown (Fig. 9). PH2O2/PHNO3<0.06 is VOC-limited region; PH2O2/PHNO3≥0.2 is NOx-limited region, and 0.06≤PH2O2/PHNO3<0.2 is the transition zone. For the pre-CLD period, during daytime, the VOC-limited (or NOx-saturated) regions included North China Plain and other urban areas, while NOx-limited regions located in southern China and other rural areas (Fig. 9a). During nighttime, most regions are VOC-limited (Fig. 9c).

The O3 formation regime determines the response of O3 to the NOx reduction during the CLD period. During daytime, NOx reduction increased O3 in NOx-saturated regions, but decreased it in NOx-limited regions (Fig. 10b). We also found that although the daytime O3 formation regime in most regions shifted from the VOC-limited regime to the NOx-limited regime after the emission reductions during the CLD period, the daytime O3 formation in the North China Plain was still controlled by the VOC level (Fig. 9b), which suggests that the NOx level is still high in this region. During nighttime, the reduction of NOx emission contributed increased O3 due to the NO titration effect in a large areas (Fig. 10c). The reduction of VOC emission decreased O3 across China (Fig. 10d-f). As an O3 precursor, the reduction of CO emission also contributed to a small decrease in the O3 mixing ratio (Fig. 10g-i); in contrast, the PM and SO2 emissions reductions increased O3 (Fig. 10j-o) through the weakening of aerosol effects (Li et al., 2019a;[START_REF] Liu | Worsening urban ozone pollution in China from 2013 to 2017 -Part 2: The effects of emission changes and implications for multi-pollutant control[END_REF], but their impacts were much smaller and were insignificant (< 1 ppbv) due to the smaller reductions, compared with NOx and VOCs.

The response of O3 to the emission reductions in different regions depended on the levels of NOx and VOC reductions. For the daytime average, in the NOx-saturated NC region, the O3 increase by the NOx reduction counteracted the O3 decrease by the VOC emission reduction, leading to the decrease in increased O3 production rates (Fig. 11) and a substantial net O3 increase (Fig. 8b). In CC, the contributions of the NOx and VOC reductions were comparable in magnitude, and their opposing impacts resulted in only a slight change in O3 (Fig. 8d). In the NOx-limited SC region, the impact of the NOx reduction on O3 was smaller than that of the reduction of VOCs, leading to the decrease in O3 production rates (Fig. 11), and a net decrease in O3 (Fig. 8f). During nighttime, the effect of the VOC reduction was weakened due to the lower rate of degradation of VOCs by radicals compared with daytime, and the O3 level increased in all three regions due to decreases in the NO titration effect (Fig. 11). The impacts of emission reductions on whole-day average O3 were similar to those during daytime.

The above modeling results show that the contribution of NOx reductions (by 46%-55%) to the rise of O3 decreased from NC to CC and to SC, reflecting the decreasing level of NOx saturation from north to south. In contrast, the impact of the estimated VOC reduction on O3 increased from north to south, which can in part be attributed to the regional variation of VOC reductions.

In the SC region, transportation and industry are the predominant sources of VOCs (97%, compared with 85% and 60% in NC and CC, respectively) (Fig. 7). During the CLD period, the reduction of VOC emission in SC (45%) was significant and comparable with the NOx reduction (55%). In contrast, the VOC reductions in the NC (32%) and CC (24%) regions were much lower (Fig. 8a,c) and could not offset the impact of NOx reduction on O3. The residential sector (mainly household coal burning) is an important source of VOC emission in the NC and CC regions, whereas its contribution is smaller in SC. The residential emissions increased during the CLD period because many migrant workers came back for the Chinese New Year holiday and were stranded there due to the lockdown (Wang et al., 2020a;[START_REF] Liu | Worsening urban ozone pollution in China from 2013 to 2017 -Part 2: The effects of emission changes and implications for multi-pollutant control[END_REF].

Conclusion

The first country-wide lockdown during the COVID-19 outbreak in China drastically reduced transportation and industrial activities, leading to sharp declines in air pollutant emissions from these sectors. Surface O3 in urban areas of China responded differently in the northern (increase) and southern regions (decrease) compared to the three-week period before the lockdown, which can be explained by changes in meteorology and differences in the O3 chemistry regimes and the magnitudes of precursor reductions in these regions. The model simulated contributions of meteorology to daytime O3 changes were larger or comparable to most regions. The extent of VOC reduction, which suppressed O3 formation, was insufficient to offset the large NO titration effect during daytime in northern China, and that larger reductions of VOCs (e.g., from residential sectors) would have been needed to reduce the O3 concentration in the northern and central China. The rising O3 concentration in northern China during the COVID-19 lockdown and in recent winters should receive greater attention because O3 boosts the atmospheric oxidative capacity and therefore production of secondary aerosols [START_REF] Fu | Persistent Heavy Winter Nitrate Pollution Driven by Increased Photochemical Oxidants in Northern China[END_REF][START_REF] Huang | Enhanced secondary pollution offset reduction of primary emissions during COVID-19 lockdown in China[END_REF][START_REF] Zhu | The seesaw impacts between reduced emissions and enhanced AOC on O3 during the COVID-19[END_REF], which are important components of winter haze in northern China. Our findings in China are relevant to untangling the underlying factors driving the O3 changes in other parts of the world during their COVID-19 lockdowns.

Figure 1 :

 1 Figure 1: Location of 1643 environmental monitoring stations (red "+" symbols) operated by the Ministory of Ecology and Environmental

Figure 2 :

 2 Figure 2: Observed changes in O3 mixing ratios across mainland China before and during the COVID-19 lockdown period. (a, c, e) The

Figure 3 :Figure 4 :

 34 Figure 3: Percentage change of (a) observed daytime average Ox (NO2+O3), whole-day average (b) NO2, (c) CO, (d) SO2, (e) PM2.5, and (f)

Figure 5 :

 5 Figure 5: Averaged sea-level pressure during the pre-CLD and CLD periods. Data are from the National Center for Environmental Prediction

Figure 6 :

 6 Figure 6: Model simulated changes in daytime temperature at 2 m height, specific humidity at 2 m height, wind field at 10 m height, planetary

Figure 7 :

 7 Figure 7: Percentage contribution to NOx, VOCs, CO, PM, and SO2 emissions from industrial (IND), power plant (POW), residential (RES),

Figure 8 :

 8 Figure 8: The estimated reductions of multi-pollutant emissions due to the COVID-19 lockdown and their impacts on the O3 changes for

Figure 9 :

 9 Figure 9: Ozone formation regime in the daytime and nighttime before and during the lockdown periods estimated by the ratio of the

Figure 10 :

 10 Figure 10: Model simulated changes in O3 mixing ratios for all-day average, daytime average, and nighttime average due to the reductions

Figure 11 :

 11 Figure 11: O3 chemical production rates before and after the anthropogenic emission reductions and the changes during the COVID-19

Acknowledgements

This work was supported by the Hong Kong Research Grants Council (T24-504/17-N and A-PolyU502/16) and the National Natural Science Foundation of China (91844301). We would like to thank Prof. Qiang Zhang from Tsinghua University for providing the emission inventory.

Data availability

Materials & Correspondence

Correspondence and requests for materials should be addressed to T.W. or Y.M.L.

The codes and data used in this study are available upon request from Yiming Liu (liuym88@mail.sysu.edu.cn) and Tao Wang (cetwang@polyu.edu.hk).

Author

Competing interests

The authors declare that they have no conflict of interest.

Additional information

Supplementary information is available for this paper.

Reference


Alexander Heyde

Lijie Guo

Christian Jost

Guy Theraulaz

L Mahadevan
email: lmahadev@g.harvard.edu

Self-organized biotectonics of termite nests

Keywords: termite nests, stigmergy, partial differential equations, morphogenesis

The termite nest is one of the architectural wonders of the living world, built by the collective action of workers in a colony. Each nest has several characteristic structural motifs that allows for efficient ventilation, cooling, and traversal. We use tomography to quantify the nest architecture of the African termite Apicotermes lamani, consisting of regularly spaced floors connected by scattered linear and helicoidal ramps. To understand how these elaborate structures are built and arranged, we formulate a minimal model for the spatiotemporal evolution of three hydrodynamic fields-mud, termites, and pheromones-linking environmental physics to collective building behavior using simple local rules based on experimental observations. We find that floors and ramps emerge as solutions of the governing equations, with statistics consistent with our observations of A. lamani nests. Our study demonstrates how a local self-reinforcing biotectonic scheme is capable of generating a architecture that is simultaneously adaptable and functional, as well as relevant for a range of other animal-built structures.

collective animal behavior

T ermite nests are among the most complex and impressive structures produced by animal societies [START_REF] Hansell | Animal architecture[END_REF][START_REF] Perna | When social behaviour is moulded in clay: on growth and form of social insect nests[END_REF][START_REF] Lee | Caste composition and mound size of the subterranean termite macrotermes gilvus (isoptera: Termitidae: Macrotermitinae)[END_REF] and serve to create a microniche that allows for the controllable exchange of matter and energy with the environment [START_REF] Turner | On the mound of macrotermes michaelseni as an organ of respiratory gas exchange[END_REF][START_REF] King | Termite mounds harness diurnal temperature oscillations for ventilation[END_REF][START_REF] Sa Ocko | Solarpowered ventilation of african termite mounds[END_REF]. Even when compared with the many animal species that exhibit collective behaviors, including coordinated motion, active synchronization, or shared decision making [START_REF] Camazine | Selforganization in biological systems[END_REF][START_REF] Sumpter | Collective animal behavior[END_REF], the collective activities of social insects such as the highly eusocial termite are exceptional, in that they often lead to the formation of intricate physical structures: shelters [START_REF] Anderson | Intermediate-level parts in insect societies: adaptive structures that ants build away from the nest[END_REF], trail networks [START_REF] Tj Czaczkes | Trail pheromones: an integrative view of their role in social insect colony organization[END_REF], and, most prominently, nests [START_REF] Hansell | Animal architecture[END_REF][START_REF] Perna | When social behaviour is moulded in clay: on growth and form of social insect nests[END_REF]. These structures are not merely the byproduct of animal behavior, however, since they also play a central role in regulating the flow of information necessary for their own construction and function. Nest building in social insects is thus considered an exemplar of functional selforganization, and studying this process can inform us about how structure and function emerge on multiple length and time scales in ecophysiology [START_REF] Bonabeau | Self-organization in social insects[END_REF].

The eusocial termites of the genus Apicotermes typify the architectural complexity that can be produced by social insects. These termites, native to the savannahs and forests of Africa, construct small ovoid nests just 20-40 cm in diameter, located 5-50 cm underground [START_REF] Desneux | Les constructions hypogées des apicotermes termites de l'afrique tropicale[END_REF]. The nest structure must be efficiently constructed so as to be capable of passive ventilation and cooling [START_REF] Turner | On the mound of macrotermes michaelseni as an organ of respiratory gas exchange[END_REF][START_REF] King | Termite mounds harness diurnal temperature oscillations for ventilation[END_REF][START_REF] Sa Ocko | Solarpowered ventilation of african termite mounds[END_REF][START_REF] Schmidt | Functions of apicotermes nests[END_REF][START_REF] Korb | The architecture of termite mounds: a result of a trade-off between thermoregulation and gas exchange?[END_REF][START_REF] Singh | The architectural design of smart ventilation and drainage systems in termite nests[END_REF], while remaining habitable and traversable by the termites within. This is no small feat-given the scale and complexity of Apicotermes nests, it is not feasible that any single termite has a sufficiently broad cognitive map of the organization of the whole structure to coordinate the actions of the colony [START_REF] Theraulaz | A brief history of stigmergy[END_REF]. Yet even in the absence of centralized control, these nests clearly show co-herent structure at a global scale, pointing to the collective self-organization of several thousand termites acting only on their local conditions to produce global order [START_REF] Lee | Caste composition and mound size of the subterranean termite macrotermes gilvus (isoptera: Termitidae: Macrotermitinae)[END_REF][START_REF] Bonabeau | Self-organization in social insects[END_REF][START_REF] Theraulaz | A brief history of stigmergy[END_REF].

Rather than sharing information by direct communication (such as via antennal contact or visible gestures), it is instead likely that Apicotermes workers communicate indirectly by coupling their building actions with the deposition of a local stimulus (such as a secreted building pheromone) in nearby substrate [START_REF] Pp Grassé | La reconstruction du nid et les coordinations interindividuelles chezbellicositermes natalensis etcubitermes sp. la théorie de la stigmergie: Essai d'interprétation du comportement des termites constructeurs[END_REF]. This principle of behavioral and physical coupling during collective construction has already been implicated in the morphogenesis of the large termite mound structures which enclose the nest. Over the long length and time scales of the overall mound structure, the small internal tunnels and chambers within the mound influence the porosity of the mound material [START_REF] Perna | When social behaviour is moulded in clay: on growth and form of social insect nests[END_REF][START_REF] Singh | The architectural design of smart ventilation and drainage systems in termite nests[END_REF], mediate their ventilation [START_REF] King | Termite mounds harness diurnal temperature oscillations for ventilation[END_REF][START_REF] Sa Ocko | Solarpowered ventilation of african termite mounds[END_REF], and suffice to explain the coarse features of the mound, e.g. its shape, size, and temporal evolution [START_REF] Sa Ocko | Morphogenesis of termite mounds[END_REF]. However, this leaves open the question of how the internal structural motifs are formed and maintained at shorter length scales in the mound microstructure. More specifically, the architectural complexity of Apicotermes nests begs two important questions: First, what are the structural motifs within the nest, and how are they arranged? And second, how are they formed via the dynamic feedback loop between individual behavior and emerging structure in the nest construction process?

Digitization and structure of Apicotermes nests

To answer these questions, we first collected, scanned, and statistically analyzed nests of the species Apicotermes lamani.

Significance Statement

Termite nests are a remarkable example of functional selforganization that show how structure and function emerge on multiple length and time scales in ecophysiology. To understand the process by which this arises, we document the labyrinthine architecture within the subterranean nests of the African termite Apicotermes lamani and develop a simple mathematical model that relies on the physical and biological interactions between mud, termites, and pheromones in the nest. Our model explains the emergent formation of parallel floors connected by linear and helical ramps, consistent with observations of real nests. Overall the study of how physics and behavior couple together in architecture might serve as inspiration for the design and construction of sustainable human architectures.

D R A F T

Two nests (MeMo14, Fig. 1A, and MeMo13, Fig. 1B) were collected near Libreville, Gabon (fig. S1), although one of these (MeMo13) showed too much internal damage for a statistical analysis. A third nest (MeMo80, Fig. 1C) was collected in 2008 near Pointe Noire in the Republic of the Congo (fig. S1); in this nest, we measured that workers had length 4-5 mm and height 2-3 mm, while soldiers had length 8-10 mm with a height of 2-3 mm. We digitized and imaged the excavated nests with a medical X-ray computed tomography (CT) scanner (Materials and Methods), allowing for non-destructive exploration of the complete three-dimensional structure of each nest.

All of the nests showed a similar external appearance, and based on our X-ray CT scans, also shared a similar internal architecture consisting of many floor layers arranged in parallel, with several vertical or diagonal pillars connecting adjacent floors (Fig. 1A-B). We recorded a video fly-through of these complex, 3D internal structures (Movie S1). To investigate these structures that enable termites to efficiently traverse their nests, we systematically searched the digitized A. lamani nests for open corridors connecting adjacent floors (Materials and Methods), and we identified several occurrences of two types of ramps (Fig. 1C): linear ramps in which an excavated path connects several floors along a linear incline, and helicoidal ramps, in which the path spirals around a central pillar while connecting adjacent floors. Helicoidal ramps in the nest had no preferred chirality, although a pair of co-occurring left-handed and right-handed helicoidal ramps was identified.

We then examined the internal architecture of two A. lamani nests (MeMo80 and MeMo14, see also Table S1) in order to quantify the spacing and arrangement of floors and ramps. For each nest, we assembled a series of 359 and 177 vertical slices respectively (Fig. 2A) and measured the thickness of each floor and the vertical spacing between floors (Fig. 2B). The floor thickness was on average 1.7 mm (median = 1.6 mm) in nest MeMo80, and on average 1.8 mm (median = 1.2 mm) in nest MeMo14. The floor spacing was on average 4.6 mm (median = 4.7 mm) in nest MeMo80, and on average 7.2 mm (median = 7.7 mm) in nest MeMo14. For both nests, the distribution of floor thickness had low variance (coefficient of variation (CV) = 0.31 in MeMo80, CV = 0.59 in MeMo14), and the distribution of floor spacing was similarly narrow (CV = 0.19 in MeMo80, CV = 0.22 in MeMo14), demonstrating that the flooring of the termite nests are arranged at regular vertical frequencies, with a strong continuity between adjacent floors (Fig. S1B-C).

Given this consistent arrangement of floors, we next investigated the positioning of ramps within the nests by measuring the horizontal distance from each ramp to the nearest ramp on the same floor (Fig. 2C). In MeMo80, we found a mean samefloor distance of 27.4 mm (median = 25.4 mm, CV = 0.41), and in MeMo14 this mean distance was 25.7 mm (median = 27.1 mm, CV = 0.33). We compared these measurements with the horizontal distance between ramps on adjacent floors. Because these ramps were often directly connected, this distance was much smaller for both nests, with a mean of 12.9 mm (median = 9.2 mm, CV = 0.79) in MeMo80, and a mean of 7.1 mm (median = 6.5 mm, CV = 0.47) in MeMo14. Hence, while ramps tended to be spaced at fairly large, regular intervals across a floor, they tended to be positioned near to ramps on adjacent floors, allowing for efficient vertical traversal between nest levels.

A minimal model of termite tectonics

To address how these structures are formed, we turn to a series of observations, old and new, to motivate a minimal continuum theory of nest construction. With no evidence of design or a designer, termites build in response to local cues such as the nearby mound structure and secreted pheromones [START_REF] Pp Grassé | La reconstruction du nid et les coordinations interindividuelles chezbellicositermes natalensis etcubitermes sp. la théorie de la stigmergie: Essai d'interprétation du comportement des termites constructeurs[END_REF][START_REF] Oh Bruinsma | An analysis of building behaviour of the termite Macrotermes subhyalinus (Rambur)[END_REF]. The resulting architecture both enables and constrains the movement of pheromones and termites, and thereby modifies behavior, so that nest construction can be seen as a result of a feedback loop linking physical and behavioral dynamics (Fig. 3A), similar to a recent model for the macrostructural morphogenesis of the termite mound [START_REF] Sa Ocko | Morphogenesis of termite mounds[END_REF]. The architecture of the nest dictates which spaces are accessible to termite workers, the density of termite workers in turn controls the concentration of secreted pheromone, and the information carried by the pheromone profile in turn serves as a template for the ongoing remodeling of the nest architecture, thereby completing the feedback loop of nest morphogenesis.

We therefore model the spatiotemporal dynamics of three fields that depend on location denoted by the vector x and time t: the nest material density u(x, t), the termite worker density n(x, t), and the pheromone concentration ρ(x, t). These three fields jointly evolve according to a set of conservation laws, expressed as the following partial differential equations for the nest material, termite workers, and pheromone levels:

Nest material: ∂tu = ∂z • [g∂zu] + f+ -f-, [1]
Termite workers:

∂tn± = ∇ • [(1 -u)D∇n± + χn±∇u] ± (f--f+)k -1 , [ 2]
Pheromone: ∂tρ = ∇ • [δ∇ρ] + Hf+ -γρ, [START_REF] Lee | Caste composition and mound size of the subterranean termite macrotermes gilvus (isoptera: Termitidae: Macrotermitinae)[END_REF] where ∂t and ∂z are the differential operator in time and height respectively, ∇• is the three-dimensional divergence operator, and ∇ is the three-dimensional gradient operator. Here Eq. 1 reflects the addition and removal of the nest material, with f+ and f-denoting the building and removal rates of dirt (whose functional form is discussed below). Moreover, the poroelastic diffusivity of dirt g captures its capacity to settle under the influence of gravity [START_REF] Skotheim | Dynamics of poroelastic filaments[END_REF]; this term breaks the rotational symmetry of the model equations and establishes a defined vertical orientation for the mound. Eq. 2 reflects the dynamics of termites carrying dirt n+ or not carrying dirt n-, with a flux that has two components: an effective diffusivity that is proportional to the amount of open space, i.e. (1 -u)D, as well as a chemotactic term that drives termite workers into open and low-density nest regions with chemotactic coefficient χ. The final term represents switching between n+ and n-, with the pellet size k denoting the average amount of dirt transported by a termite worker. Lastly, Eq. 3 reflects the dynamics of secreted pheromones, driven by diffusion with a diffusion coefficient δ, production by termites at a rate H, and degradation at a rate γ. We note that the pheromone is relatively non-volatile and thus not advected by the ambient fluid, since the subterranean nests of Apicotermes do not typically have noticeable temperature gradients [START_REF] Schmidt | Apticotermes nests[END_REF].

To complete the formulation of our model, we require functional forms for the building rate f+ and removal rate f-. Although these rates have not been directly measured in A. lamani termites, the building rate has been investigated for the neotropical termite Procornitermes araujoi [START_REF] Fouquet | Coordination of construction behavior in the termite procornitermes araujoi: structure is a stronger stimulus than volatile marking[END_REF], and both rates have also been quantified in the ant Lasius niger [START_REF] Khuong | Stigmergic construction and topochemical information shape ant nest architecture[END_REF]. In all cases we find that the local nest density scales directly with the building rate and inversely with the removal rate, resulting in high-density regions becoming increasingly compacted while low-density regions are excavated. Moreover the amplification of building rate scales with a pheromone added by the workers to the building material (Fig. 3B), so that the termites are more likely to deposit dirt in regions marked by pheromone as being sites of active remodeling, and so that building does not take place either where there is no nest material (u = 0) or where the nest material is already fully compacted (u = 1). It is reasonable also that the removal rate scales with the density of termite workers not carrying dirt and with the density of dirt u available for removal. Given these basic considerations, the simplest possible functional forms for f+ and f-are Building rate: f+ = r+n+ρu(1 -u), [START_REF] Turner | On the mound of macrotermes michaelseni as an organ of respiratory gas exchange[END_REF] Removal rate: f-= r-n-u, [START_REF] King | Termite mounds harness diurnal temperature oscillations for ventilation[END_REF] where r± are rate constants. With these choices, the dynamics of nest morphogenesis as described by our model can be captured by only a small number of non-dimensional parameters: the scaled pheromone potency Hr+/r-, the scaled evaporative flux γ/r-, the pellet size k, and the scaled diffusivities δ/D, χ/D, and g/D (Materials and Methods). For simplicity, and to study the self-contained organization of this system, we model the nest as a closed three-dimensional domain, using no-flux boundary conditions, such that when integrated over the nest domain the total number of termites |n| and the total quantity of dirt u are conserved.

Natural patterns in simulated mounds

We conducted 500 independent simulations of nest construction by initializing our model each time with uniform nest density u(x, t) and pheromone concentration ρ(x, t), and with a uniform density of termite workers n(x, t), using measured and inferred parameter values from several studies on termite nest structure and behavior (Table S2) and a numerical implementation of an Euler differencing scheme (Materials and Methods). The model generated simulated nests (Fig. S2A-C) containing regularly spaced floors, with a mean planar orientation orthogonal to the direction of gravity and a floor thickness that scales as (δ/γ) 1/2 in some regimes (Materials and Methods, Fig. S2D), which gives an average prediction of ∼1.9 mm, directly in line with our observed measurements of 1.7 mm and 1.8 mm in MeMo80 and MeMo14, respectively. Moreover, the simulated nests contained surface edges (topological defects) spanning across floors and sometimes forming traversable connections (Fig. 3C) corresponding to linear ramps (edge dislocations)-in which one floor terminates to provide the vertical space necessary for a simple ramp between adjacent floors-and helicoidal ramps (screw dislocations), which wind about a vertical pillar serving as a dislocation line spanning multiple floors (24) (Fig. S3).

By identifying the eventual location of a ramp, we observe the dynamical process by which ramps emerge from these defects during the construction process (Fig. 4A,Movie S2). Varying our model parameters by an order of magnitude on either side of our estimates, we find that only two parameters control whether helicoidal ramps form: the scaled pheromone potency Hr+/r-and the scaled pheromone evaporation rate γ/r-, yielding a phase space for their emergence (Fig. 4B). This is consistent with previous agent-based simulations of ant nest construction that pointed to the central role of the pheromone evaporation rate in influencing nest structures [START_REF] Khuong | Stigmergic construction and topochemical information shape ant nest architecture[END_REF].

A fly-through of our simulated nests (Movie S3) reveals the structural similarity of the simulated nests produced by our model with the scanned A. lamani nests. To quantitatively compare the simulated and scanned nests, we gathered an identical set of statistics to describe the simulated nest topologies. We found that the floor thickness and spacing between floors were similarly consistent in both the simulated and scanned nests, indicated by distributions with a similarly narrow variance (Fig. 4C). Moreover, the horizontal spacing between neighboring ramps on the same floor was substantially greater than the horizontal spacing between ramps on adjacent floors that often connected directly, as observed in the scanned nests (Fig. 4D). To assess the similarity of the simulated and scanned nests on large scales, we carried out a spectral analysis of the average nest density as a function of height for both the scanned and simulated nests and found that these density profiles are very similar (Fig. 4E). Together, these metrics indicate that the simulated nests resemble natural nests across global and local scales.

Emergent biotectonics from physics and behavior

Our minimal theoretical framework links three spatiotemporal fields involved in the biotectonics of termite nests: dirt to constitute the nest, termite workers to shape the nest, and secreted pheromone to mark active regions of the nest. They allow us to capture two key geometrical and topological features of termite nests, namely the regular vertical spacing of floors and the horizontal spacing of ramps on the same or adjacent floors. This feedback loop involving architecture, behavior, and information that drives nest construction is likely to be quite general; it is implicated also in the macromorphogenesis of termite mounds [START_REF] Sa Ocko | Morphogenesis of termite mounds[END_REF] and might additionally provide insight into emergence of other animal architectures, such as ant nests, beehives or other termite subfamilies [START_REF] Hansell | Animal architecture[END_REF][START_REF] Khuong | Stigmergic construction and topochemical information shape ant nest architecture[END_REF][START_REF] Škarka | Mathematical model of building behavior of apis mellifera[END_REF][START_REF] Perna | Topological efficiency in three-dimensional gallery networks of termite nests[END_REF]. The finding that the nest structure emerges spontaneously for a range of parameter values points to it being a nonequilibrium steady state that is robust and adaptable-a direct consequence of the ability of a colony to both lay down and respond to cues in its environment.

An important open question is the role of heat gradients and air flow within the nest [START_REF] King | Termite mounds harness diurnal temperature oscillations for ventilation[END_REF][START_REF] Sa Ocko | Solarpowered ventilation of african termite mounds[END_REF][START_REF] Singh | The architectural design of smart ventilation and drainage systems in termite nests[END_REF], which are determined by the nest structure and in turn help to propagate pheromones and other signals, such as carbon dioxide levels [START_REF] Cox | Gaseous templates in ant nests[END_REF]. The pheromones we consider here have low diffusivities (19) so as to provide a localized signal useful in nest construction, such that heat-driven air flow results in a negligible change when incorporated into our model (SI Text, Fig. S4); however, other odors with a higher diffusivity, such as those secreted by the queen, would be more susceptible to internal nest flow, and hence more important to the larger-scale problem of the global templating and emergence of overall mound shape [START_REF] Sa Ocko | Morphogenesis of termite mounds[END_REF][START_REF] Oh Bruinsma | An analysis of building behaviour of the termite Macrotermes subhyalinus (Rambur)[END_REF]. In the subterranean nests of Apicotermes, we expect that the role of heat gradients and air flow is minor relative to above-ground nests which are directly subjected to sunlight and air currents.

The spontaneous emergence of these functional, palatial architectures from the construction dynamics of social insects interacting with the environment stands in stark contrast to planned human architectures that break down on moderate to large scales. Perhaps by studying the natural structures that termites and other insects build and rely upon, and by under-standing the coupling of physical and behavioral processes that underlie their morphogenesis, we can take inspiration for the design and construction of sustainable human architectures [START_REF] Turner | Beyond biomimicry: What termites can tell us about realizing the living building[END_REF]. After all, while we as humans have a mere ∼400,000 years of experience with constructing primitive structures (29), mound-building termites of the family Termitidae have spent close to 50 million years [START_REF] Da Arab | Parallel evolution of mound-building and grass-feeding in australian nasute termites[END_REF] developing one of the grandest examples of architecture in the natural world.

Materials and Methods

X-ray computed tomography (CT). Each nest was imaged using X-ray computed tomography with a medical scanner and reconstructed into a series of virtual cuts (512×512 pixels). MeMo13 and MeMo14 were scanned with a Somatom Sensation16 (Siemens, Erlangen, Germany) at CHU Toulouse Rangueil, France using exposure parameters of 120 kV and 150 mAs, slice thickness of 1 mm, interslice distance 0.5 mm. MeMo80 was scanned with a LightSpeed Ultra (GE Medical Systems, Buc, France) at CHU Dijon using exposure parameters 100 kV and 170 mAs, slice thickness of 0.625 mm, interslice distance 0.3 mm. Reconstruction of the virtual cuts (Fig. 1) was performed in the open-source software Horos (LGPL license).

Investigation of internal structures. Ramps (with connecting holes to the floor above) and pillars (supporting the above floor without a passway) were identified by visual inspection of three orthogonal planes at each voxel of the original CT in the Horos software. For helices, identification was confirmed by reconstructing a perspective view (3D Surface Rendering) that could be inspected from all sides (e.g. Fig. 1C). Mean coordinates for the base of each ramp or pillar were extracted (Fig. 2A) for the statistics in Fig. 2C and Table S1. For Fig. 2B, MeMo14 was rotated to align horizontal floors with the x-y plane (imageJ/Fiji, rotate tool with bilinear interpolation). Images were then binarized in imageJ with a threshold calibrated by visual comparison between the resulting structures, the original grey-scale tomograph images, and photos of cut-open same-species nests. The outer walls were removed by first dilating each slice five

D R A F T

times to close all outside openings, and then using the wand tool to define an outline ROI of the nest for each slice (also used to compute the total nest volume). All pixels outside the new ROIs were removed to obtain nests without outer walls as in Fig. 2A.

Floor thickness and floor spacing were measured by counting the number of consecutive black or white pixels in each vertical pixel line of the stacked images, with code written in R. The presence of pillars/ramps led to heavy right tails in both distributions, which were removed by cutting off all measures above twice the median floor height and alpha-trimming the resulting distributions with α = 0.02 on both sides. To check for the continuation of ramps across several floors, we selected for each floor the horizontal slice for which the horizontal pixel density was minimal (Fig. S1B), computed the correlation coefficient between the pixels of adjacent floors, and compared them to the correlation coefficients between slices at least three floors separated (Fig. S1C).

Non-dimensionalization of construction model. Using the dimensional parameters given in Table S2, we can construct five nondimensional parameters that allows us to rewrite Eqs. 1-5 in terms of a series of normalized variables T = r -t, X = x/ , ∇ = ∇, and P = ρr + /r -for time, position, gradient, and pheromone concentration respectively, where = (r -/D) 1/2 denotes the length scale used here for normalization. This gives the non-dimensional system

∂ T u = ∂ Z • [ g D ∂ Z u] + f + -f -, [6] ∂ T n ± = ∇ • [(1 -u) ∇n ± + χ D n ± ∇u] ± (f --f + )k -1 , [7] ∂ T P = ∇ • [ δ D ∇P ] + Hr + r - f + -γ r - P, [ 8]
with f + = n + P u(1 -u) and f -= n -u. These coupled equations contain a set of dimensionless parameters that are the scaled pheromone potency Hr + /r -, the scaled evaporative flux γ/r -, the pellet size k, and the scaled diffusivities δ/D, χ/D, and g/D that control the range of possible outputs from our model.

Analysis of construction model. Our model permits the derivation of a simple scaling prediction regarding the regularity of floor thickness. At steady state, the building and removal rates must be balanced:

f + -f -= r + n + ρu(1 -u) -r -n -u = 0. [9]
In terms of the steady state ratio η = n -/n + of termite workers without dirt, this balance leads to either a no-density tunnel or chamber (u = 0) or a high-density floor or wall with density value

u * = 1 - ηr - ρr + . [10]
To understand which of these steady states is stable we need to evaluate the stability condition ∂u(f + -f -) < 0 (see SI). This reveals that the tunnels (u = 0) occur when the pheromone level is low (ρ < ηr -/r + ); since in these regions f + = f -= 0, the pheromone profile converges to ρ(x) = 0. In contrast, structures (u = u * ) occur when the pheromone level is high (ρ > ηr -/r + ), and in the vicinity of a structure, building and removal are balanced at rate f + = f -= r -n + ηu * .

In the SI, we show that regularly spaced floors require that the condition 0 ≤ α ≤ 1 needs to be satisfied, where α = 1 -4γ/Hr + . This allows us to derive an expression for the thickness L of a floor which follows the approximate scaling law

L ∼ γ δ - Hr + 4δ (1 - √ α) -1/2
. [START_REF] Bonabeau | Self-organization in social insects[END_REF] In the regime α ≈ 1, for example, the floor thickness scales as L ∼ δ/γ, the prediction one would expect in the case where floors are fully packed, u * ≈ 1. This regime is characterized by slow pheromone evaporation relative to its high deposition rate and potency. For our parameter values, this thickness is a realistic 1.9 mm. For the alternative extremal case α ≈ 0, in which the pheromone evaporation is sufficiently high to result in more diffuse nest interiors, the scaling law depends also on H and r + , as plotted in Fig. S2D, consistent with our full numerical simulations.

Numerical solution of model equations.

To simulate our construction model, we implemented a finite difference solver to numerically integrate the system of equations [START_REF] Hansell | Animal architecture[END_REF][START_REF] Perna | When social behaviour is moulded in clay: on growth and form of social insect nests[END_REF][START_REF] Lee | Caste composition and mound size of the subterranean termite macrotermes gilvus (isoptera: Termitidae: Macrotermitinae)[END_REF][START_REF] Turner | On the mound of macrotermes michaelseni as an organ of respiratory gas exchange[END_REF][START_REF] King | Termite mounds harness diurnal temperature oscillations for ventilation[END_REF]. We initialized each simulation with uniform nest density u = 0.5 (halfway packed) and uniform pheromone concentration ρ = 0.1 ng/cm 3 , and with randomly scattered termite workers, such that the density of termites n at each grid location x was drawn from an exponential distribution with non-dimensional mean 0.1. At each time step, the gradient and Laplacian of each of the three fields was calculated using second-order differencing in space, and the nest interior was then updated with first-order differencing in time. For each set of model parameters, 500 independent simulations were conducted.

Fig. 1 .

 1 Fig. 1. Digital reconstruction of Apicotermes lamani nests reveals layered floors and chambers connected by linear and helicoidal ramps. (A-B) Cutaway view of two A. lamani nests, collected around Libreville (Gabon): MeMo14 (A) and MeMo13 (B). The nests were digitized with X-ray computer tomography and reconstructed in 3D. (C) A large A. lamani nest (left, MeMo80) collected in 2008 near Pointe Noire (Republic of the Congo) with examples of a linear ramp (i, red dots) and two helicoidal ramps with opposing chirality (ii, red and cyan dots).

Fig. 2 .

 2 Fig. 2. Statistical analysis of A. lamani nests shows consistent floor and ramp spacing. Top row = nest MeMo80, bottom row = nest MeMo14. (A) Representative vertical slices of nest structure with floors labeled in average height order. (B) Histograms for floor thickness (pink) and spacing between floors (cyan), as measured in mm, corresponding to the nest depicted at left. Both the thickness of and spacing between floors tend to fall within a consistent band in each nest. (C) Density plot for the horizontal distance from a ramp to the nearest other ramp on the same floor (red) or on an adjacent floor above or below (blue), corresponding to the nest depicted at left. Ramps on the same floor tend to be spaced out, while ramps on adjacent floors often connect directly, resulting in minimal spacing.

Fig. 3 .

 3 Fig. 3. Biotectonic model predicts floor spacing and ramp emergence in termite nests. (A) Model schematic of the feedback loop driving nest construction, highlighting the interactions between nest material u, termite workers n, and secreted pheromone ρ. (B) Illustration of a local region of a nest, showing the processes in our model. Termite workers migrate preferentially to low-density regions and cannot travel through very high-density regions. Workers remove dirt throughout the nest but are more likely to deposit dirt near pheromones which they release during deposition. Pheromones are assumed to have a low diffusivity and hence provide a local signal. (C) 3D reconstruction of a nest simulated according to our construction model, shown here at two angles. This simulated nest contains one linear ramp (i, blue) and two helicoidal ramps (ii-iii, red).

Fig. 4 .

 4 Fig. 4. Simulated nests resemble natural A. lamani nests and produce helicoidal ramps across a range of parameters. (A) Time snapshots of helicoidal ramp emergence during simulated nest construction, from early (top) to late (bottom) in the building period. (B) Heatplot for the frequency of helicoidal ramps in a simulated nest as a function of the two model parameters governing pheromone dynamics. (C) Histograms for floor thickness (pink) and spacing between floors (cyan), as measured in mm, averaged over simulated nests. The histograms show a pattern similar to the natural nests in Fig. 2B. (D) Density plot for the horizontal distance from a ramp to the nearest other ramp on the same floor (red) or on an adjacent floor above or below (blue), averaged over simulated nests and resembling the natural patterns shown in Fig. 2C. (E) The power spectrum of nest density, averaged over horizontal slices, peaks sharply at a period of two floor heights for both simulated (gray) and natural (purple) nests, indicating regularly spaced floor structures.

6

  | www.pnas.org/cgi/doi/10.1073/pnas.XXXXXXXXXX

	621
	622
	623
	624
	625
	626
	627
	628
	629
	630
	631
	632
	633
	634
	635
	636
	637
	638
	639
	640
	641
	642
	643
	644
	645
	646
	647
	648
	649
	650
	651
	652
	653
	654
	655
	656
	657
	658
	659
	660
	661
	662
	663
	664
	665
	666
	667
	668
	669
	670
	671
	672
	673
	674
	675
	676
	677
	678
	679
	680
	681
	682

Heyde et al.

Heyde et al.

ACKNOWLEDGMENTS.

We thank A. Robert for providing nest MeMo80 and the natural history museum in Paris for providing nests MeMo13 and MeMo14, both collected by Grassé. This work was supported by NSF Grant DGE-1144152 (to A.H.), ANR Grant ANR-06-BYOS-0008 (to G.T.), CSC PhD Grant (to L.G.), and the NSF Physics of Living Systems Grant PHY1606895 (to L.M.).