H Agurto-Detzel

Y Font

P Charvis

M Régnier

A Rietbrock

D Ambrois

M Paulatto

A Alvarado

S Beck

F Courboulex

L De Barros

A Deschamps

M J Hernandez

S Hernandez

M Hoskins

S León-Ríos

C Lynner

A Meltzer

D Mercerat

F Michaud

J M Nocquet

F Rolandone

M Ruiz

L Soto-Cordero

Ridge Subduction and Afterslip Control Aftershock 1 Distribution of the 2016 Mw 7.8 Ecuador Earthquake

We characterise the aftershock sequence following the 2016 Mw=7.8 Pedernales earthquake.

13 More than 10,000 events were detected and located, with magnitudes up to 6.9. Most of the aftershock seismicity results from interplate thrust faulting, but we also observe a few normal and strike-slip mechanisms. Seismicity extends for more than 300 km along strike, and is16 constrained between the trench and the maximum depth of the coseismic rupture. The most striking feature is the presence of three seismicity bands, perpendicular to the trench, which are also observed during the interseismic period.

Additionally, we observe a linear dependency between the temporal evolution of afterslip and aftershocks. We also find a temporal semi-logarithmic expansion of aftershock seismicity along strike and dip directions,further indicating that their occurrence is modulated by afterslip. Lastly, we observe that the spatial distribution of seismic and aseismic slip processes is correlated to the distribution of bathymetric anomalies associated with the northern flank of the Carnegie Ridge, suggesting that slip in the area could be influenced by the relief of the subducting seafloor. To explain our observations, we propose a conceptual model in which the Ecuadorian margin is subject to a bimodal slip mode, with distributed seismic and aseismic slip mechanically controlled by the subduction of a rough oceanic relief. Our study sheds new light on the mechanics of subduction, relevant for convergent margins with a complex and heterogeneous structure such as the Ecuadorian margin.

Introduction

The largest earthquakes on Earth occur in subduction zones, which also host a diversity of processes including seismic and aseismic slip along the subduction interface (e.g. Bilek and Lay, 2018, and references therein). What controls the occurrence and distribution of these phenomena remains an outstanding problem in Earth sciences. One way to gain a better insight into the nature of the subduction mechanism and the physical medium that host them, is by studying the aftershocks sequence that follows a large megathrust earthquake.

Moreover, the high rate of seismicity during aftershock sequences, combined with recent technological and logistical improvements in seismological network deployments and data processing (e.g. [START_REF] Beck | Advancing Subduction Zone Science After a Big Quake Access and Use of Data Lead to Innovative Ways[END_REF], allows us to collect and analyse vast amounts of data with increased spatio-temporal resolution.

Aftershocks occur either because of the release of residual stresses on the mainshock fault and surrounding medium, or as a result of static or dynamic stress perturbations due to the coseismic rupture and subsequent aftershocks (e.g. [START_REF] Das | Spatial relation between main earthquake slip and its aftershock distribution[END_REF][START_REF] Freed | Earthquake Triggering By Static, Dynamic, and Postseismic Stress Transfer[END_REF].

Consequently, aftershocks can provide an independent constraint in the shape and extension of the rupture area and interface heterogeneities, as well as help us identify areas of partially released and/or accumulated stress over the megathrust interface following the mainshock, thus delineating potential source areas for future earthquakes.

The often intricate distribution of aftershocks accounts for a complex distribution of remaining stresses and interface heterogeneities following the mainshock. For instance, after the 2005 Mw=8.7 Nias-Simeulue earthquake in Sumatra, [START_REF] Hsu | Frictional afterslip following the 2005 Nias-Simeulue earthquake[END_REF] found that aftershocks clustered in the boundary area between the coseismic rupture and the afterslip area, with afterslip concentrated mostly up-dip of the coseismic rupture. Furthermore, it is often observed that regions of large co-seismic slip tend to have little seismicity after the mainshock rupture, whilst the largest aftershocks concentrate around the patches of large coseismic slip (e.g. [START_REF] Das | Spatial relation between main earthquake slip and its aftershock distribution[END_REF][START_REF] Rietbrock | Aftershock seismicity of the 2010 Maule Mw=8.8, Chile, earthquake: Correlation between co-seismic slip models and aftershock distribution?[END_REF][START_REF] Agurto | Seismic-afterslip characterization of the[END_REF][START_REF] Wetzler | Systematic deficiency of aftershocks in areas of high coseismic slip for large subduction zone earthquakes[END_REF]. On the other hand, aftershock activity is not only limited to the megathrust interface, but also to the surrounding seismogenic volume, often showing a diversity of focal mechanisms and complex interactions between activity in the slab and in the overriding plate (e.g. [START_REF] Asano | Spatial distribution and focal mechanisms of aftershocks of the 2011 off the Pacific coast of Tohoku Earthquake[END_REF]. Lastly, for some subduction earthquakes, such as the 2011 Tohoku, Japan earthquake, the reduction of shear stresses after the mainshock is such that it produces a rotation of the deviatoric stress field, potentially causing extensional earthquakes in a previously compressional setting (e.g. [START_REF] Ryder | Large extensional aftershocks in the continental forearc triggered by the 2010 Maule earthquake, Chile[END_REF][START_REF] Hardebeck | Coseismic and postseismic stress rotations due to great subduction zone earthquakes[END_REF].

Moreover, the physics behind aftershock generation is still not fully understood. Aftershocks were first described and used as a proxy for the mainshock rupture extension, and subsequently explained as ruptures on surrounding faults due to the re-distribution of strain energy following the mainshock. Consequently, aftershocks triggering mechanism would be related to dynamic and/or static stress transfers, following the mainshock and subsequent aftershocks [START_REF] Stein | The role of stress transfer in earthquake occurrence!: Abstract!: Nature[END_REF]. More recently, observational and theoretical studies have proposed that afterslip plays an important role in the occurrence and distribution of aftershocks (e.g. Henry and Das, 2001;[START_REF] Perfettini | A Model of Aftershock Migration Driven by Afterslip[END_REF]. For example, following the 2005 Mw=8.7 Nias-Simeulue earthquake in Sumatra, [START_REF] Hsu | Frictional afterslip following the 2005 Nias-Simeulue earthquake[END_REF] found that the cumulative number of aftershocks increased linearly with the postseismic displacement, suggesting that the temporal evolution of aftershocks is governed by afterslip.

On the 16 of April 2016, a Mw=7.8 earthquake struck the coast of northern Ecuador rupturing a ~100 km-long asperity of the interface between the Nazca plate and South America [START_REF] Nocquet | Supercycle at the Ecuadorian subduction zone revealed after the 2016 Pedernales earthquake[END_REF]. Shortly after the mainshock, we deployed an amphibious temporary network of seismic stations to monitor the evolution of the seismic activity. In this paper, we benefit from the continuous seismic waveform dataset acquired during one year of the aftershock deployment to explore the distribution of hypocentral locations and magnitudes for the Pedernales sequence. We also use full waveform inversions to compute moment tensors for a selection of events, providing a seismotectonic constraint to the characterization of the sequence. We discuss our results in the light of the earthquake cycle, exploring the relations between seismic and aseismic processes within the context of a subduction zone with highly heterogeneous frictional properties. Finally, we present a conceptual model in which we explain the distribution and diversity of slip processes in the Ecuadorian margin, and the control factors that affect them.

Seismotectonic context and previous studies

The Ecuador-Colombia subduction margin has generated four large tsunamigenic megathrust earthquakes (Mw > 7.5) in the 20 th century. In 1906, an Mw ~8.8 event (the largest thus far documented offshore Ecuador) ruptured a roughly 500 km-long segment of the margin, causing widespread damage and tsunami waves [START_REF] Kanamori | Variable rupture mode of the subduction zone along the Ecuador-Colombia coast[END_REF]. Subsequent events occurred in 1942 (Mw 7.8), 1958 (Mw 7.7) and1979 (Mw 8.2;[START_REF] Kanamori | Variable rupture mode of the subduction zone along the Ecuador-Colombia coast[END_REF][START_REF] Beck | The rupture process of the great 1979 Colombia earthquake: evidence for the asperity model[END_REF], partially overlapping the rupture area of the 1906 event. This sequence of three earthquakes presented a northward migration pattern (Fig. 1), and the sum of their combined seismic moments accounts for only a fifth of the moment released by the 1906 event (Keller, 1972;[START_REF] Kanamori | Variable rupture mode of the subduction zone along the Ecuador-Colombia coast[END_REF]. This would imply that the 1906 event not only ruptured the other three isolated asperities simultaneously, but also broke the adjacent subduction interface which otherwise creeps during the interseismic period.

The area that ruptured in 2016 had already been identified as a highly coupled region [START_REF] Chlieh | Distribution of discrete seismic asperities and aseismic slip along the Ecuadorian megathrust[END_REF][START_REF] Nocquet | Motion of continental slivers and creeping subduction in the northern Andes[END_REF], and the same asperity had allegedly been ruptured by the earthquake of 1942 [START_REF] Nocquet | Supercycle at the Ecuadorian subduction zone revealed after the 2016 Pedernales earthquake[END_REF]. In this region, the convergence rate between Nazca and South America is 58 mm/yr, which is partially accommodated by the north-eastern motion of the North-Andean sliver, resulting in a slip rate of 46 mm yr -1 at the megathrust [START_REF] Chlieh | Distribution of discrete seismic asperities and aseismic slip along the Ecuadorian megathrust[END_REF][START_REF] Nocquet | Motion of continental slivers and creeping subduction in the northern Andes[END_REF]. Also, this area is located within the northern flank of the aseismic Carnegie Ridge (hereafter CR), which currently subducts beneath South America between 0° to 2.5° lat. S.

To date, several co-seismic slip models of the 2016 earthquake have been published based on a complete or partial use of teleseismic, tsunami, GPS, InSAR and regional accelerometric data (e.g. [START_REF] Ye | The 16 April 2016, MW7.8 (MS7.5) Ecuador earthquake: A quasi-repeat of the 1942 MS7.5 earthquake and partial re-rupture of the 1906 MS8.6 Colombia-Ecuador earthquake[END_REF][START_REF] Nocquet | Supercycle at the Ecuadorian subduction zone revealed after the 2016 Pedernales earthquake[END_REF][START_REF] Yoshimoto | Depth-dependent rupture mode along the Ecuador-Colombia subduction zone[END_REF][START_REF] Gombert | Strain budget of the Ecuador-Colombia subduction zone: A stochastic view[END_REF] and references within). All models have in common an extension of the rupture area of roughly 100 km along strike, a southward propagation rupture, and the presence of two patches of high coseismic slip with no shallow slip near the trench. They differ, however, in the maximum and average amount of slip, with maximum slip ranging from 2 m [START_REF] Yoshimoto | Depth-dependent rupture mode along the Ecuador-Colombia subduction zone[END_REF] to 6-7 m [START_REF] Nocquet | Supercycle at the Ecuadorian subduction zone revealed after the 2016 Pedernales earthquake[END_REF][START_REF] Gombert | Strain budget of the Ecuador-Colombia subduction zone: A stochastic view[END_REF]. These last two models are very similar regarding magnitude and distribution of the co-seismic slip, and are the most comprehensive up to date in terms of diversity of used datasets and methodology.

Previous studies using geodetic and seismological data highlight the diverse nature of slip processes in the interseismic period. [START_REF] Font | Seismicity patterns along the ecuadorian subduction zone: New constraints from earthquake location in a 3-D a priori velocity model[END_REF] produced a seismicity catalogue for a 13-yr period based on locations in a 3-D a priori velocity model. [START_REF] Vallée | Intense interface seismicity triggered by a shallow slow slip event in the Central Ecuador subduction zone[END_REF] characterized a one-week-long slow slip event (SSE), accompanied by a seismic swarm, that occurred in August 2010 below La Plata Island (hereafter LPI), south of the 2016 rupture.

Similarly, [START_REF] Vaca | Recurrent slow slip events as a barrier to the northward rupture propagation of the 2016 Pedernales earthquake (Central Ecuador)[END_REF] described a six-week-long SSE accompanied by a seismic swarm that occurred between December 2013 and January 2014 at the northern limit of the 2016 rupture, arguing that this area acted as a barrier for the 2016 rupture propagation northwards. Finally, [START_REF] Segovia | Seismicity Distribution Near a Subducting Seamount in the Central Ecuadorian Subduction Zone, Space-Time Relation to a Slow-Slip Event[END_REF] studied the seismicity distribution during a twoyear experiment in the south of the region, describing the interface geometry, and associating swarm-like activity to a SSE below LPI.

Data and Methods

Earthquake rapid response deployment

Following the Pedernales earthquake, an international effort involving institutions from Ecuador (IG-EPN), France (Géoazur, Cerema, IRD and CNRS), the UK (U. of Liverpool) and the USA (IRIS, U. of Lehigh, U. of Arizona) rapidly installed a network of 50 inland stations and 10 ocean-bottom seismometers (OBS) to record for one year after the mainshock (Fig. 2; [START_REF] Meltzer | The 2016 Mw 7.8 Pedernales Earthquake, Ecuador: RAPID Response Deployment[END_REF]. This temporary deployment complemented the permanent Ecuadorian network [START_REF] Alvarado | Seismic, Volcanic, and Geodetic Networks in Ecuador: Building Capacity for Monitoring and Research[END_REF]. Instruments included broadband, intermediate and short period stations, in addition to some accelerometers from the Ecuadorian network, all recording at a sampling rate of 100 Hz or higher.

Data processing

The continuous waveforms were collected and archived in mini-seed format. They were processed using the software package SEISCOMP3 (SC3; https://www.seiscomp3.org) which provides in-built capacity to detect, associate and locate seismic events including the calculation of magnitudes. Although SC3 is primarily designed for real-time monitoring with continuous injection of data, it can also be used in 'playback mode', that is, injecting and processing the whole of the collected data at once. Parameterization of the different SC3 modules is critical, and therefore we adopted an empirical approach in which several tests were systematically performed looking for the best set of parameters that would maximize the number of real events while minimizing the number of false detections. Control days, for which we manually detected events, were used to asses this fine-tuning process. Additionally, we visually inspected the detected events and discarded false detections as well as classified real events into first and second quality events according to the number and accuracy of their automatic picks (see Sup. Mat.).

The workflow was as follows: after injection of the continuous waveform dataset, detection of arrival times was performed using a standard STA/LTA algorithm for P-phases and the AIC picker implemented in SC3 for S-phases, on band-pass filtered waveforms (1-10 Hz for seismometers; 1-8 Hz for accelerometers and OBS). Subsequently, we used the SC3 module SCANLOC, which is based on the cluster-search algorithm DBSCAN (Easter et al., 1996), to associate picks and locate events. Relocation of these initial events was performed using the NonLinLoc (NLL) algorithm [START_REF] Lomax | Probabilistic Earthquake Location in 3D and Layered Models[END_REF] configured in standard global mode. The visual quality-inspection described above was carried out on these preliminary locations. Finally, the whole set of events was relocated outside SC3 using NLL configured in regional mode (Cartesian coordinates) and a simplified velocity model taken from a newly derived 1-D velocity structure for the region (León-Ríos et al., 2017; see Sup. Mat.).

Initially, a total of 15,233 aftershocks were detected and located for the period between April 16 2016 to April 30 2017. Visual analysis of seismic sections was performed to identify and discard false detections, spurious events, and to assess the pick quality to assemble a highquality subset of events. After this visual inspection, a total of 4,963 (33%) events were discarded as false detections or spurious events (poor signal and/or too noisy). The 10,270 events left were classified into two categories according to their picks and location quality:

1 st quality: events with at least four P-phases and clear arrival picks -7326 events 2 nd quality: events with pick residuals larger than ~2 s and greater location errors -2944 events

Moment tensor inversions

We selected aftershocks with M L >4.5 to compute moment tensors from full waveform inversions, including both body and surface waves. For this we used the software ISOLA (Sokos and Zaharadnik, 2008) which can handle inversions of local to regional waveforms.

Green functions were computed using the 1-D model produced by León-Ríos et al. (2017) and waveforms were inverted in the 10 -25 s period range. Horizontal centroid position was kept fixed to the epicentral position from the earthquake locations, but a grid-search was performed to obtain optimal centroid depth and time. Examples of the inversion and fitting are provided in the Supplementary Material.

Magnitudes

Local magnitudes (M L ) were calculated from maximum P-wave amplitudes on vertical components. The obtained magnitudes vary between 0.7 and 6.9, with a magnitude of completeness Mc=2.5 (Sup. Mat.). In general, there is a good agreement between the calculated local magnitudes (M L ) and the moment magnitudes (Mw) obtained from our moment tensor inversions and those from the GCMT catalogue (Sup. Mat.). Nonetheless, we observe that for Mw>5.6 there is an underestimation of local magnitudes, probably due to saturation of the M L scale. On the other hand, for Mw<5.6 we observe an overestimation of M L by ~0.3 units. These differences are commonly observed when comparing local magnitudes with moment magnitudes (e.g. [START_REF] Deichmann | Local Magnitude, a Moment Revisited[END_REF].

Residual bathymetry

In order to compare the distribution of seismicity with the distribution of the incoming oceanic relief, we produced a residual bathymetry grid for the Ecuadorian margin following the ensemble averaging approach of [START_REF] Basset | Gravity anomalies, crustal structure, and seismicity at subduction zones: 1. Seafloor roughness and subducting relief[END_REF], and using the higher resolution GEBCO2014 grid. We calculated the average topography for a series of trench-normal profiles. Then we subtracted this averaged topography from the original grid to produce an elevation map where large-amplitude trench-normal variations associated with the subduction zone have been removed and short-wavelength/lower amplitude structures are preserved and highlighted.

Spatio-temporal distribution of aftershocks

Along strike, the aftershock seismicity extends beyond the coseismic rupture, over 300+ km, from latitude 1°N to at least 1.5°S (Fig. 3). Along the dip direction, the seismicity seems to be constrained by the coseismic rupture maximum depth, with most of the aftershocks located in the upper 30 km and no aftershock seismicity locates deeper than the coseismic rupture termination.

The most striking feature is the presence of three bands of seismicity perpendicular to the trench, and located up-dip west of the mainshock rupture area (profiles BB', CC' and DD' in Fig. 3; see also [START_REF] Soto-Cordero | New insights on co-and post-seismic deformation and slip behavior associated with the Mw7.8 2016 Pedernales, Ecuador earthquake and its aftershock sequence[END_REF]. Interestingly, this seismicity pattern is also observed in the background seismicity during the interseismic period [START_REF] Font | Seismicity patterns along the ecuadorian subduction zone: New constraints from earthquake location in a 3-D a priori velocity model[END_REF]. The northern band (BB') extends for about 40 km up-dip with a width of about 10 km. The central band (CC') is more diffuse, starting at the upper termination of the rupture area and extending 40 km up-dip with a width of around 20 km. The southern band also starts at the upper termination of the coseismic rupture and extends up-dip 60 km with a width of about 25 km. Both, the southern and central bands reach the trench, whilst seismicity is more diffuse close to the trench for the northern band. Although we do observe seismicity near the trench, we do not observe any extensional focal mechanism in this area that could be related to outer rise seismicity following the mainshock (e.g. [START_REF] Sladen | Shallow megathrust earthquake ruptures betrayed by their outer-trench aftershocks signature[END_REF].

Considering the location uncertainties, most of the seismicity in these three alignments occurs at the interface or within 10 km from it. Additionally, all large aftershocks (M>=5) occur outside the mainshock rupture and mostly along bands BB' and DD', located up-dip at the northern and southern limits of the co-seismic rupture, respectively. Inside the mainshock rupture area, seismicity occurs mostly between the two patches of maximum coseismic slip (Figure 5, see Section 5).

To the north (0.9° N), we observe a cluster of seismicity within the subduction interface below the coastline (cluster G1 in Fig. 3). Further to the east, a cluster of crustal seismicity (G2, hereafter called Esmeraldas sequence) is observed at 10-20 km depth. This group of shallow seismicity started to develop at the end of June 2016, with a burst of seismicity during July 5-8 and its largest earthquake, normal faulting Mw=4.9, occurring on July 6, 2016 (see details in Section 4).

South of the mainshock rupture area, we observe three separate groups of seismicity. The first one is a cluster of events occurring beneath the coastline, at around latitude 0.9°S (G3).

This cluster seems to occur on the megathrust interface, and as seen in Section 4, presents thrust focal mechanisms compatible with subduction earthquakes. The second group corresponds to the seismicity observed inland at around latitude 1.3°S (G4) which also occur at the interface. The third group (G5) is located offshore, nearby LPI. This seismicity is sparsely distributed, and because of its location offshore at the southern end of the network, it is difficult to assess hypocentral depths with certainty. Nevertheless, a clue regarding the origin of this seismicity comes from previous studies which have found swarm-like seismicity and SSEs in this area (Vallee et al., 2013;[START_REF] Segovia | Seismicity Distribution Near a Subducting Seamount in the Central Ecuadorian Subduction Zone, Space-Time Relation to a Slow-Slip Event[END_REF], as well as a SSE during the early postseismic period of the 2016 mainshock [START_REF] Rolandone | Areas prone to slow slip events impede earthquake rupture propagation and promote afterslip[END_REF]. Like the trench-normal bands, these three seismicity groups had also been observed during the interseismic period [START_REF] Segovia | Seismicity Distribution Near a Subducting Seamount in the Central Ecuadorian Subduction Zone, Space-Time Relation to a Slow-Slip Event[END_REF].

The spatio-temporal analysis of the aftershock sequence (Sup. Video) shows that during the first 24 hours after the mainshock, aftershocks start to nucleate mostly along profiles DD' and CC', and in particular between the two patches of maximum co-seismic slip. The aftershocks then extend along profiles BB' and EE'. Seismicity around LPI starts on the third day with peaks of activity on the 11 th and 12 th days after the mainshock. A last burst of seismicity in this area occurs between 1 st and 3 rd December 2016. As stated above, the shallow clustered seismicity of the Esmeraldas sequence occurs mostly during early July 2016. Finally, the seismicity observed at the interface along the profile AA' develops during December 2016.

Seismotectonics and moment tensor inversions

For the 12-month period following the Pedernales mainshock (April 16 2016 -April 30 2017) there are 32 moment tensors with Mw between 4.8 and 6.9 available in the GCMT catalogue (http://www.globalcmt.org/). We complemented these with 29 additional events with Mw"between"4.1"and"5.0,"for"a"total"of"61"moment"tensors"(Fig. "4"and"Sup."Mat.)."Most" of"the"moment"tensors"indicate"thrust"faulting"at"the"subduction"interface."No"large" aftershocks"(Mw>5)"occur"inside"the"coseismic"rupture"area."The"largest"thrust"aftershocks" occur"along"the"seismicity"bands"located"at"the"northern"and"southern"termination"of"the" mainshock"rupture."Besides"these"two"bands"dominated"by"thrust"faulting"at"the"interface," we"also"observe"subduction"earthquakes"to"the"south,"around"latitude"1°S,"and"towards"the" north"by"the"coastline"up"to"1°N."The"geometry"of"the"reverse"faulting"focal"mechanisms"is" similar"to"that"of"the"mainshock,"with"an"average"rotational"angle"(Kagan"angle)"of"22°" relative"to"the"mainshock's"focal"mechanism"(inset"Fig. "4)."""" " We"also"observe"a"few"normal"and"strikeNslip"events."StrikeNslip"events"seem"to"be"sparsely" located"and"within"the"subducting"slab."A"possible"explanation"for"this"activity"could"be"the" presence"of"preNexisting"structures"in"the"subducting"CR,"reactivated"by"the"mainshock."On" the"other"hand,"two"similar"normal"fault"earthquakes,"of"Mw"5.1"and"4.9"respectively," occurred"in"the"marine"forearc"around"latitude"0.3°N,"on"June"1 st "2016,"separated"by"5" hours."The"GCMT"centroid"depths"for"these"earthquakes"(12"and"17"km"depth)"place"them" close"to"the"subduction"interface,"but"on"our"own"regional"moment"tensor"inversions"we" found"the"lowest"waveform"misfit"at"5"km"depth."Despite"the"depth"uncertainties,"a"possible" explanation"for"this"faulting"could"be"given"by"the"existence"at"this"location"of"a"subducted" seamount,"previously"imaged"using"multiNchannel"seismic"reflection"data" [START_REF] Marcaillou | Seamount subduction at the North-Ecuadorian convergent margin: Effects on structures, inter-seismic coupling and seismogenesis[END_REF]." LeónNRíos"et"al."(2017)"hypothesize"that"the"subduction"of"this"structure"produces"an" anomalous"extensional"stress"field"parallel"to"the"convergence"vector,"which"in"turn"could" have"been"affected"by"the"2016"mainshock."In"fact," [START_REF] Marcaillou | Seamount subduction at the North-Ecuadorian convergent margin: Effects on structures, inter-seismic coupling and seismogenesis[END_REF]"observed"a" complex"and"highly"fractured"margin"structure"in"this"region,"and"argued"that"the"absence" of"background"seismicity"and"low"interseismic"coupling"here"suggest"that"this"area"is" incapable"of"storing"sufficient"elastic"strain"to"produce"large"thrust"earthquakes"and" tsunamis.""" """""""" Two"additional"normal"fault"events"are"observed"in"our"dataset."One"is"a"Mw=4.4" intermediateNdepth"event,"most"likely"intraNslab,"located"at"0.6°N,"200"km"east"of"the" trench."The"other"is"a"Mw=4.9,"crustal"normal"fault"event"with"a"strikeNslip"component," belonging"to"the"Esmeraldas"sequence."Unfortunately,"the"uncertainties"of"our"hypocentral" locations"in"this"area"do"not"allow"us"to"distinguish"the"fault"plane"from"the"two"nodal" planes."On"the"other"hand,"the"geological"map"for"this"area" [START_REF] Reyes | Mapa Geológico de la Margen Costera Ecuatoriana (1:50000[END_REF];"Sup." Mat.)"shows"a"set"of"normal"faults"striking"ESE"and"dipping"to"the"S,"which"coincide"with"one" of"the"nodal"planes"of"this"event"(strike"103,"dip"42,"rake"N29)."We"suggest"that"crustal" activity"on"one"of"these"faults"might"be"responsible"for"the"seismicity"observed"during"the" Esmeraldas"sequence"(see"also "Hoskins"et"al.,"2018)."Some"previous"large"megathrust" aftershock"sequences,"such"as"Maule"2010"and"Tohoku"2011,"have"shown"similar"shallow" normal"faulting"at"the"edges"of"the"coseismic"rupture"area"(e.g. [START_REF] Kato | A normal-faulting seismic sequence triggered by the 2011 off the Pacific coast of Tohoku Earthquake: Wholesale stress regime changes in the upper plate[END_REF]"Ryder"et" al.,"2012)."A"similar"tectonic"configuration"could"be"responsible"for"our"normal"event"in"the" Esmeraldas"area,"which"indicates"horizontal"extension"in"the"overriding"plate"following"the" mainshock."Since"these"events"are"shallow,"near"the"coast,"and"can"produce"considerable" vertical"displacement,"they"are"important"to"consider"when"estimating"earthquake"and" tsunami"hazard"at"a"local"scale.""""

The April 16 2016 Mw=4.9 foreshock Nearly 11 minutes before the Pedernales earthquake, an event M W =4.9 nucleated about 14 km ESE of the mainshock's epicentre. We also obtained the moment tensor for this event, which indicates a thrust faulting mechanism, likely on the subduction interface (Fig. 4). The possibility of this earthquake to have triggered the Mw=7.8 mainshock is worth exploring, although a dynamic or static triggering would be difficult to reconcile with the time and distance between the two events. More accurate relocations of both the foreshock and main event hypocentres, and a detailed analysis of the Coulomb stress change field, beyond the scope of this study, would be necessary to resolve this issue.

Relation between coseismic rupture and aftershock distribution

As a first order feature, we observe an inverse correlation between the number of aftershocks and co-seismic slip, with highs in slip associated to lows in seismicity and vice versa (e.g. at 20, 45 and 60 km south of the mainshock in profile N-S of Fig. 3). Figure 5 shows in detail the distribution of aftershocks and co-seismic slip. We observe that most of the large aftershocks occur outside the mainshock rupture area (defined as the 1 m slip contour area).

When we consider all magnitudes, 28% of the aftershocks occur inside the mainshock rupture, but when we consider only events with M L ≥3.5, only 14% of aftershocks nucleate inside and, moreover, no aftershock larger than M L =5 nucleated inside the mainshock rupture area.

Additionally, the histograms in Fig. 5 show the normalized areal distribution of co-seismic slip together with the normalized frequency distribution of aftershocks inside the coseismic rupture. Accordingly, if the aftershocks occurrence were randomly distributed, the aftershock frequency curve would resemble the slip frequency distribution. Instead, we observe that aftershocks tend to concentrate at intermediate levels of coseismic slip (2 -3.5 m), particularly in areas of large slip gradient, such as in between the two patches of coseismic slip maxima. On the other hand, areas of low coseismic slip (< 2 m) present less seismicity than expected, whilst areas of high coseismic slip (> 4.5 m) seem to present a random distribution of aftershocks (histogram Fig. 5a), although when we consider only events with M L ≥3.5, there is a lack of aftershocks compared to a random distribution (histogram Fig. 5b).

If we look at the aftershock density, we observe that in terms of number of events, the highest density is located inside the mainshock rupture area, in between the two patches of maximum coseismic slip (Fig. 6a). If instead we look at the seismic moment density (Fig. 6b), we observe that inside the mainshock rupture area the moment density is relatively low (< 1e 17 N m / 0.1°x0.1°). On the other hand, high moment density (> 1e 18 N m / 0.1°x0.1°) is observed outside the mainshock rupture, along the three trench-normal seismicity bands and particularly nearby the coastline at latitude 0.5°N, due to the occurrence here of the largest aftershock of the sequence (Mw=6.9, thrust faulting).

Relation between seismic and aseismic processes

We compare the temporal evolution of the aftershock sequence with that of the geodetic afterslip during the first 30 days following the mainshock. Following [START_REF] Rolandone | Areas prone to slow slip events impede earthquake rupture propagation and promote afterslip[END_REF], we consider the afterslip and aftershocks as three discrete patches (North, South and LPI; see Fig. 3 and Sup. Mat) according to their spatial distribution, and analyse them separately (Fig. 7). Cumulative seismicity (red curve) for the northern and southern patches show an Omoritype decay in which a steep slope is observed immediately after the mainshock, followed by a deceleration after the first week of aftershocks. On the other hand, the LPI patch shows a rather slow start in aftershocks generation, and then an increase from day 8 until day ~20 when it decreases again. The different behaviour in the LPI patch could be explained because this area hosted a slow slip event associated to seismicity during this period [START_REF] Rolandone | Areas prone to slow slip events impede earthquake rupture propagation and promote afterslip[END_REF].

We observe for all three patches that the curve for cumulative number of earthquakes closely follows that of the afterslip cumulative moment release, implying a linear relationship between both processes. In fact, if we assume that both afterslip and aftershocks cumulative distributions present an exponential behaviour, their curves should resemble a straight-line in a semi-logarithmic plot, as seen in the right panels of Fig. 7, which also show both curves present similar slopes (segmented lines). Leaving the LPI patch aside, the linear relation between cumulative aftershocks and afterslip release is remarkable.

Furthermore, for each of the patches we observe that after 30 days of postseismic activity, the total cumulative moment released by the aftershocks represents about 10% of the cumulative moment released by the afterslip, indicating that most of the postseismic deformation is aseismic (Sup. Mat).

Additionally, we explore the spatial dependency between afterslip and aftershocks. As seen from the geographic distribution of seismicity outside the mainshock rupture area, aftershocks are spatially associated with afterslip, particularly in the area of the trench-normal bands and around LPI (Fig. 3). Figure 8 shows the temporal evolution of seismicity as a function of along-strike distance from the mainshock epicentre, clearly showing a log-time expansion of the aftershocks. A similar behaviour is seen for the along-dip direction (Sup. Mat). These observations are consistent with previous studies (e.g. [START_REF] Frank | Mapping the rheology of the Central Chile subduction zone with aftershocks[END_REF], and numerical modelling (e.g. [START_REF] Ariyoshi | The key frictional parameters controlling spatial variations in the speed of postseismic-slip propagation on a subduction plate boundary[END_REF][START_REF] Perfettini | A Model of Aftershock Migration Driven by Afterslip[END_REF] which suggest that this type of semilogarithmic migration is indicative of afterslip-driven aftershock activity.

Discussion

Where do aftershocks occur?

The largest aftershocks occur outside the mainshock rupture area. This finding is in agreement with previous studies which have found that regions of high coseismic slip are mostly devoid of large aftershocks, whilst post-seismic seismicity tends to concentrate at the edges of the coseismic rupture (e.g. [START_REF] Das | Spatial relation between main earthquake slip and its aftershock distribution[END_REF][START_REF] Asano | Spatial distribution and focal mechanisms of aftershocks of the 2011 off the Pacific coast of Tohoku Earthquake[END_REF]Rietbrock et al., 2012;[START_REF] Agurto | Seismic-afterslip characterization of the[END_REF][START_REF] Frank | Mapping the rheology of the Central Chile subduction zone with aftershocks[END_REF][START_REF] Wetzler | Systematic deficiency of aftershocks in areas of high coseismic slip for large subduction zone earthquakes[END_REF]. For the 2010 Mw=8.8

Maule earthquake, Agurto et al. ( 2012) also found that aftershocks concentrated at intermediate levels of coseismic slip, with areas of low and large coseismic slip lacking in aftershocks. Therefore, this could be a common feature for large megathrust earthquakes with a heterogeneous distribution of coseismic slip.

Additionally, a large number of aftershocks do occur within the co-seismic rupture area, although presenting low magnitudes. The fact that aftershocks nucleate inside the mainshock rupture area indicates that the accumulated strain energy within the fault is not totally released during the mainshock, or at least that this release is not homogeneously distributed along the megathrust rupture. Attempting to investigate this issue, [START_REF] Yabe | Why Do Aftershocks Occur Within the Rupture Area of a Large Earthquake?[END_REF] produced quasi-dynamic numerical simulations in which they replicate several megathrust frictional scenarios and mainshock ruptures with their respective aftershock sequences. They observed aftershocks around and within the mainshock rupture area for cases in which frictional heterogeneity varies significantly along the fault. On the other hand, aftershocks were not produced when frictional heterogeneities along the fault were small. Similarly, the fact that for the Pedernales sequence we observe the highest density of aftershocks within the mainshock rupture area, might be indicative of the highly heterogeneous distribution of frictional properties along the northern Ecuador megathrust.

When we account for location uncertainties, the low-magnitude seismicity located within the co-seismic rupture area seems to occur distributed within the seismogenic volume and not only at the megathrust interface (Fig. 3). This volume represents the off-fault damage zone produced by successive megathrust ruptures, and it usually concentrates a diversity of aftershocks focal mechanisms in structures re-activated by the mainshock (e.g. [START_REF] Asano | Spatial distribution and focal mechanisms of aftershocks of the 2011 off the Pacific coast of Tohoku Earthquake[END_REF][START_REF] Agurto | Seismic-afterslip characterization of the[END_REF].

What controls the evolution of the aftershock seismicity?

The temporal linear dependency between afterslip and aftershocks shown here (Fig. 7) suggests a causative time-based relationship between these two processes, and therefore the temporal distribution of aftershocks associated to patches of afterslip would be modulated by the stressing rate associated with afterslip (e.g. [START_REF] Perfettini | Postseismic relaxation driven by brittle creep: A possible mechanism to reconcile geodetic measurements and the decay rate of aftershocks, application to the Chi-Chi earthquake[END_REF][START_REF] Hsu | Frictional afterslip following the 2005 Nias-Simeulue earthquake[END_REF].

Additionally, the semi-logarithmic migration of aftershocks both along strike and dip (Fig. 8) suggests that afterslip also controls the spatial extension and migration speed of aftershocks (e.g. [START_REF] Frank | Mapping the rheology of the Central Chile subduction zone with aftershocks[END_REF][START_REF] Perfettini | A Model of Aftershock Migration Driven by Afterslip[END_REF]. We notice that the origin of the two red lines indicating the propagation front in Fig. 8 is not located at the epicentre but approximately 40 km south of it, in the area where most of the aftershock seismicity take place during the first 24 hours following the mainshock (Section 3). This corresponds to the centre of the coseismic rupture, and therefore we hypothesize that the expansion of aftershocks is initiated at this point, subsequently propagating outwards.

Another explanation for the observed aftershocks expansion could be related to fluid diffusion. Nevertheless, in such a case we would observe that the distance D associated with the migration front of the seismicity is related to time t as D~sqrt(t)c, where c is the hydraulic diffusivity coefficient [START_REF] Wang | Theory of Linear Poroelasticity[END_REF]. This is unlike our observations, in which we see that D~log(t).

Finally, we notice that a similar relationship between seismic and aseismic processes in our study area has also been described during the interseismic period [START_REF] Vallée | Intense interface seismicity triggered by a shallow slow slip event in the Central Ecuador subduction zone[END_REF][START_REF] Rolandone | Areas prone to slow slip events impede earthquake rupture propagation and promote afterslip[END_REF][START_REF] Segovia | Seismicity Distribution Near a Subducting Seamount in the Central Ecuadorian Subduction Zone, Space-Time Relation to a Slow-Slip Event[END_REF][START_REF] Vaca | Recurrent slow slip events as a barrier to the northward rupture propagation of the 2016 Pedernales earthquake (Central Ecuador)[END_REF]. These previous studies describe seismic swarms associated to SSEs in the offshore area in front of Punta Galera (lat. ~0.7°N; [START_REF] Vaca | Recurrent slow slip events as a barrier to the northward rupture propagation of the 2016 Pedernales earthquake (Central Ecuador)[END_REF], and around LPI [START_REF] Vallée | Intense interface seismicity triggered by a shallow slow slip event in the Central Ecuador subduction zone[END_REF][START_REF] Segovia | Seismicity Distribution Near a Subducting Seamount in the Central Ecuadorian Subduction Zone, Space-Time Relation to a Slow-Slip Event[END_REF]. A similar SSE around LPI occurred during the postseismic period of the 2016 Pedernales earthquake, also associated with seismicity [START_REF] Rolandone | Areas prone to slow slip events impede earthquake rupture propagation and promote afterslip[END_REF]. Therefore, it seems that the close spatio-temporal correlation between seismic and aseismic processes in this region is persistent during the whole of the earthquake cycle.

Persistent seismicity patterns over the earthquake cycle

Aseismic slip seems to modulate the rate and spatio-temporal expansion of the aftershock seismicity. But why do these slip processes occur where they occur in the first place? In our study area, the presence of persistent spatial seismicity patterns over the earthquake cycle, such as the three trench-normal bands and the seismicity south of the mainshock rupture area [START_REF] Font | Seismicity patterns along the ecuadorian subduction zone: New constraints from earthquake location in a 3-D a priori velocity model[END_REF], suggest that earthquake nucleation in these areas is somehow controlled by long-lived structural features. We also notice that the bands are dominated by thrust events (Fig. 4), and oriented perpendicular to the trench, similar to the slip vector of the mainshock, as opposed to parallel to the convergence vector.

To our knowledge, no other subduction zone presents this type of permanent seismicity pattern transcending the earthquake cycle. Observations in other tectonic settings such as Parkfield, in the San Andreas fault, show sub-horizontal alignments of seismicity along the fault plane that also persist through many seismic cycles. Because of its geometry and the motion of the fault, it has been proposed that this seismic activity is related to rheological transitions within the fault zone and/or stress concentrations between locked and creeping areas (e.g. [START_REF] Waldhauser | Streaks, multiplets, and holes: High-resolution spatio-temporal behavior of Parkfield seismicity[END_REF]. Nonetheless, invoking rheological transitions in our area is a less plausible hypothesis to explain our observations, mainly because the seismicity within the bands lacks any clear depth-dependency.

One additional hypothesis is that the interface frictional properties in these regions of high seismicity are somehow different than in the rest of the area. In this sense, the interseismic coupling map for our study region (Fig. 1) shows that the general area of the bands is only slightly coupled (<40%), but the model lacks the resolution to see any difference along strike, between areas with seismicity (bands) and areas without.

Influence of the subducting seafloor relief

Previous studies have proposed an along-strike segmentation of the Ecuadorian margin in which large subduction earthquakes only occur north of the CR, which acts as a barrier to the southward propagation of megathrust ruptures (e.g. [START_REF] Collot | Are rupture zone limits of great subduction earthquakes controlled by upper plate structures? Evidence from multichannel seismic reflection data acquired across the northern Ecuador-southwest Colombia margin[END_REF]. More generally, it has been proposed that rugged subducting seafloor, as in the case of seamounts and ridges, give rise to heterogeneous stresses, promoting creep as expressed in transient events of various spatial and temporal scales, accompanied with small and medium-sized earthquakes (Wang and Bilek, 2014). Bassett and Watts (2015) produced a compilation of residual bathymetric anomalies for several subduction zones of the world, and found that regions with subducted seamounts were correlated to reduced levels of megathrust activity, suggesting that these areas are mostly associated with small earthquakes and creep rather than with large megathrust events. Furthermore, they argue that larger bathymetric features, such as aseismic ridges, exhibit seafloor roughness over a larger scale than subducted seamounts, presenting widths comparable to the rupture length of large (Mw~7) megathrust earthquakes. They observe that the maximum roughness is located at the flanks of the ridges, which often serve as barriers of rupture propagation. For the Ecuador subduction zone, some authors observed that the northern flank of the CR has acted as a barrier against the southward propagation of the 1906 and 1942 earthquakes [START_REF] Kelleher | Rupture zones of large South American earthquakes and some predictions[END_REF][START_REF] Collot | Are rupture zone limits of great subduction earthquakes controlled by upper plate structures? Evidence from multichannel seismic reflection data acquired across the northern Ecuador-southwest Colombia margin[END_REF].

Following the ensemble averaging approach of Bassett and Watts (2015), and benefiting from combined high resolution datasets, we produced improved maps of residual bathymetry for the Ecuadorian margin. We compared the spatial distribution of these anomalies with the distribution of the seismic and aseismic processes before and after the Pedernales earthquake (Fig. 9). Landward from the trench, the down-dip limit of the area with high residual bathymetry (>2 km) coincides with the up-dip limit of the Pedernales earthquake rupture area. Bassett and Watts (2015) notice that this limit coincide with the continental slope break, and suggest that the slope break corresponds to the updip limit of the seismogenic zone, and that the outer portion of the plate interface, below the steep continental slope, is weak/conditionally stable and would slip aseismically. Furthermore, we notice that both the 1942 and the 2016 epicentres are located nearby this limit, with the 2016 mainshock rupture area itself extending down-dip from this limit, within an area of smoother residual bathymetry. We also notice that the trench-normal bands of seismicity observed during the interseismic and post-seismic periods occur in areas of higher gradient and residual bathymetry. In particular, the seismicity band DD', which marks the southern boundary of the Pedernales rupture zone, is in front of the highest bathymetric and gravity anomaly, which correspond to the thickest part of the CR crust (~20 km; [START_REF] Collot | Are rupture zone limits of great subduction earthquakes controlled by upper plate structures? Evidence from multichannel seismic reflection data acquired across the northern Ecuador-southwest Colombia margin[END_REF][START_REF] Sallarès | Seismic structure of the Carnegie ridge and the nature of the Galápagos ho[END_REF]. Lastly, both the SSEs observed during the interseismic period, and the afterslip patches observed during the post-seismic period occur in areas dominated by high residual bathymetry due to the subduction of the CR (Fig. 9). We summarize our observations in an interpretative figure (Fig. 10) in which we suggest that the Ecuadorian margin hosts a bimodal slip mode mechanically controlled by the distribution of the subducting oceanic relief. The bimodal slip mode produces seismic and aseismic slip processes, and is present both along-strike and along-dip. In the area where the CR subducts beneath the margin (latitude 0° to ~2.5°S), particularly in the region containing a high residual bathymetry (>2 km, from the trench until ~90 km landward; Zone A in Fig. 10), the overall ISC is low (<40%), and the subduction slip mode is dominated by creep and small to medium-sized earthquakes (Mw<6), swarm-like seismicity and SSEs during the interseismic phase, and aseismic afterslip during the postseismic period. Down-dip of this limit (i.e. over 90 km horizontally from the trench, down to the maximum seismogenic depth; Zone B), the ISC is higher (>40%) and the slip mode is dominated by large subduction earthquakes (Mw>7) as in the case of the 2016 Pedernales Earthquake and similar past ruptures. Along strike to the north of the ridge flank, away from the area of influence of the CR (Zone C), the overall ISC is high up to the trench, and megathrust earthquake ruptures could reach the trench, as allegedly was the case for the 1906 earthquake and possibly the 1979 earthquake. Therefore, Zone A presents an overall stable regime (velocity-strengthening) whilst Zones B and C are unstable/conditionally stable (velocity-weakening). Thus, the area of high residual bathymetry (> 2 km) would act as a barrier to up-dip (trench-normal) propagation of megathrust ruptures, whilst the lateral flanks of the ridge would act as barriers to along-strike (trench-parallel) rupture propagation.

Conclusion

We characterised the aftershock seismicity occurring in the Ecuadorian margin over one year following the 2016 Mw=7.8 Pedernales earthquake. More than 10,000 events were detected and located, with magnitudes up to 6.9. Most of the seismicity results from interplate thrust faulting but we also observe a few normal and strike-slip mechanisms. Within the mainshock rupture area, seismicity concentrates in regions of intermediate coseismic slip, particularly in between the two patches of slip maxima. Outside the rupture area, seismicity extends for more than 300 km along strike. The most striking feature is the presence of three seismicity bands, perpendicular to the trench, which are also observed during the interseismic period.

We observe a linear dependency between the temporal evolution of afterslip and number of aftershocks, confirming previous results [START_REF] Rolandone | Areas prone to slow slip events impede earthquake rupture propagation and promote afterslip[END_REF]. Additionally, aftershocks present a temporal semi-logarithmic expansion along the strike and dip directions, which further suggest their spatio-temporal occurrence is regulated by afterslip. A comparison of the distribution of seismic and aseismic slip processes with the distribution of bathymetric anomalies reveals that slip in the area seems to be controlled by the subduction of oceanic plate roughness. To explain our observations, we propose a conceptual model in which the Ecuadorian margin presents a bimodal slip mode mechanically controlled by the subduction of a rough oceanic relief. In this sense, the flanks of the CR act as a barrier to the propagation of megathrust ruptures, both up-dip and along-strike. On the other hand, the area of maximum influence of the CR (residual bathymetry > 2 km) is characterized by small magnitude earthquakes (Mw<6), aseismic slip, repeating events and earthquake swarms. Figure 1. Interseismic coupling [START_REF] Nocquet | Motion of continental slivers and creeping subduction in the northern Andes[END_REF] and main seismotectonic features. White stars and solid white lines show epicentres and approximate rupture areas of past megathrust earthquakes respectively [START_REF] Kanamori | Variable rupture mode of the subduction zone along the Ecuador-Colombia coast[END_REF]Mendoza and Dewey, 1984). Yellow star shows epicentre of 2016 mainshock together with its GCMT focal mechanism. Blue contour shows rupture area of 2016 event [START_REF] Nocquet | Supercycle at the Ecuadorian subduction zone revealed after the 2016 Pedernales earthquake[END_REF]. Black contours show depth of subduction interface every 10 km (Hayes, 2012). Segmented black line indicates the Dolores-Guayaquil Fault Zone (Collot et al., 2002). Convergence NAZ/NAS from [START_REF] Chlieh | Distribution of discrete seismic asperities and aseismic slip along the Ecuadorian megathrust[END_REF]. NAZ Nazca Plate, NAS North Andean Sliver, SAM South American Plate. We propose that the area influenced by the subduction of the CR, as shown by the residual bathymetry contours, delimits the slip mode along dip and along strike in the Ecuadorian margin. Along dip, Zone A presents a rough and highly heterogeneous interface with the presence of fractures, possible fluids and overall low coupling. The interface at Zone A would be weak and seismically stable (velocity strengthening), and its slip mode is dominated by creeping, and includes SSE, repeating earthquakes, small to medium size (M<6) earthquakes and swarm activity, including the permanent bands of seismicity. Down dip, Zone B is less influenced by the CR, presenting an overall high coupling and a smoother interface allowing for large megathrust ruptures, although contained within ~15 to 40 km depth as in the case of the 1942 and 2016 ruptures. North of the CR along strike, Zone C is out of the influence of the CR and presents overall high ISC and large (M>7.5) megathrust ruptures that occasionally can reach the trench as in the case of the 1906 earthquake. The interfaces of both Zones B and C therefore would be unstable/conditionally stable (velocity weakening).

  Figure2. Seismic network and average location errors (68% confidence) for events with confidence ellipse semi-axes less than 20 km (50% of total events). Left: average horizontal error; Right: average vertical error.

Figure 3

 3 Figure3(previous page). Aftershock locations in map view and depth sections. Light blue circles show all first quality locations; dark blue circles show high accuracy locations with ellipse semi-axis errors less than 5 km. Coseismic rupture model is shown as white contours every 1 m slip[START_REF] Nocquet | Supercycle at the Ecuadorian subduction zone revealed after the 2016 Pedernales earthquake[END_REF]. Red stars are aftershocks with M L >=5. Pink contours show afterslip every 10 cm[START_REF] Rolandone | Areas prone to slow slip events impede earthquake rupture propagation and promote afterslip[END_REF]. Clusters (G1-G5) indicate seismicity groups described in Section 3. Slab depth model (white lines in map, black line in depth sections) from Hayes, 2012. Histograms with blue bars show number of earthquakes for each profile. Histogram with orange line in N-S profile show distribution of coseismic slip along strike.

Figure 4 .

 4 Figure4. GCMT mechanisms and regional moment tensors obtained in this work. Distribution shows epicentral location from this study for all events, and depth from computed centroid depth. Inset (top) rose histogram showing strike of nodal planes for all reverse fault mechanisms. Blue segmented line shows strikes of nodal planes for mainshock. Inset (bottom) shows histogram of rotational angle relative to mainshock mechanism for all reverse fault events. For details, see also Table2in Supplementary Material.

Figure 5 .

 5 Figure 5. Distribution of aftershocks and coseismic rupture (Nocquet et al., 2017). (a) all magnitudes; (b) magnitudes equal or greater than 3.5. Histograms show normalized frequency distribution of coseismic slip (colour bars) and aftershocks (blue line).

Figure 6 .

 6 Figure 6. Density plots for (a) number of earthquakes, (b) seismic moment. Other features same as in Fig. 2.

Figure 7 .

 7 Figure 7. Temporal evolution of afterslip and aftershocks for the three different afterslip patches during the first 30 days following the mainshock. Released afterslip distribution after Rolandone et al. (2018). Left panels: cumulative distribution as a function of day. Right panels: cumulative distribution as a function of logarithm of day. Segmented line is bestfitted straight-line.

Figure 8 .

 8 Figure 8. Expansion of earthquakes along strike in function of time since the mainshock. Red line indicates semi-logarithmic migration velocity of events (drawn by hand).

Figure 10 .

 10 Figure 9. (a) Residual bathymetry and slip processes in the Ecuadorian margin. Blue contours every 500 m above 2000 m of residual bathymetry. Black box shows zoomed area in right-side panels. (b) interseismic (1943-2016) slip processes over residual bathymetry in grey scale. Seismicity from ISC catalogue. (c) postseismic slip processes (after 2016 mainshock) over residual bathymetry in grey scale.

  Table 2 in Supplementary Material.

Acknowledgments

The rapid deployment of the seismic array after the Pedernales earthquake was supported by IG-EPN and IRD in the frame of the International Joint Laboratory "Earthquakes and Volcanoes in the Northern Andes", INSU-CNRS, the University of Liverpool, Lehigh University and the University of Arizona. The US Seismic Rapid Response deployment was supported by the NSF RAPID Program Award EAR-1642498 and by the PASSCAL facility of the Incorporated Research Institutions for Seismology (IRIS) through the PASSCAL Instrument Center at New Mexico Tech with support from the National Science Foundation under Cooperative Agreement EAR-1261681 and by the Department of Energy National Nuclear Security Administration. The OBS were deployed from R/V "ORION" provided by INOCAR (May 2016) and recovered from Coast Guard vessel "LG-52" in harsh sea conditions thanks to Davide Oregioni and Deny Malengros (Géoazur). A special thank to Comandante Andres Pazmiño (INOCAR) and Esmeraldas Coast-Guard Captain Patricio Estupinian for securing shiptime and providing invaluable help to prepare OBS on-land. We are also grateful to IRD for supporting the postdoc grant of HAD and project REMAKE (ANR-15-CE04-0004) for supporting the processing and interpretation of the data. HAD acknowledges fruitful discussions with I. Manighetti and C. Twardzik. Figures were compiled using GMT software (Wessel and Smith, 1998).


Lifeng Xu
email: lfxu@zjut.edu.cn

Zhongzhu Yang
email: zhongzhu_yang@126.com

Weilong Ding
email: wlding@zjut.edu.cn

PhD Gerhard Buck-Sorlin
email: gerhard.buck-sorlin@agrocampus-ouest.fr.

Physics-based algorithm to simulate tree dynamics under wind load

Keywords: physics-based algorithm, tree model, L-system, Cantilever beam virtual plant modeling

Rapid development in the different computer science fields during the recent decades has facilitated the creation of new applications in the area of dynamic simulation of plant development. Among these new applications, simulation of trees swaying in the wind is of great importance, as those computer graphics related areas, e.g., computer games, tree cultivation and forest management simulations, help a lot in revealing the mechanisms of tree dynamics under wind load. However, it is a big challenge to balance the effect of visualization in real time and calculation efficiency for any simulation algorithm. A physics-based algorithm to simulate tree dynamics under wind load was proposed in this study. A mechanistic model simulating the bending of a cantilever beam was used within the algorithm to simulate deformation of stems, and the algorithm was integrated with a landscape model in which different types of trees were constructed with an L-system-based formalism. Simulation results show that realistic dynamic effects can be achieved with reasonably high computational efficiency.

Introduction 

Trees are an essential component of the natural environment, and virtual tree models have been widely used in many aspects such as digital games, movies, urban modeling, plant biology, forestry, and agriculture. Given such complex structures, and the large amounts of data coming along with them, it is difficult to construct a model with instant rendering and interaction between trees and the environment. With the continuous development of virtual reality theory and technology, more and more researchers are focusing on the dynamic simulation of trees in the wind [START_REF] Ancelin | Development of an individual tree-based mechanical model to predict wind damage within forest stands[END_REF] . As a non-linear flexible system, it is relatively difficult to simulate the deformation of a tree in the wind accurately and quickly, thus making it a challenging topic in computer graphics. Previous methods [START_REF] Sakaguchi | Modeling and animation of botanical trees for interactive virtual environments[END_REF] can only be used within relatively simple tree models, in which the control nodes and shape fitting patches are relatively small, mostly counted in hundreds of control nodes and thousands of shape fitting patches. Obviously, the use of such a simple model to generate the natural sway of complex trees is difficult; the dynamic deformation of the calculation accuracy and imaging of the realism will be less precise.

Virtual plant technology uses computers to simulate the growth process of plants, and represents the morphological structure of plants with a 3D model. Virtual plants can help us to research plants with a new perspective. It can not only reveal the complex morphological structure and life characteristics of plants but also promote the development of computer modeling and animation theory. Virtual plant modeling approaches include two major categories, i.e. graphical-based approaches [START_REF] Cannon | The fractal geometry of nature[END_REF] and image-based approaches [START_REF] Teng | tree modeling from a few images with very narrow viewing range[END_REF] . Aono and Kunii [START_REF] Aono | Botanical tree image generation[END_REF] made a detailed estimate of the distribution of the stem and the leaves on the parent by using the tree as biological object model, and developed an interactive synthetic tree operating system called A-system. Oppenheimer [START_REF] Oppenheimer | Real time design and animation of fractal plants and trees[END_REF] used fractals to construct a tree model, adjusting the parameters to change the geometric and topological shapes of the simulated object, considering the random effects of the environment to avoid the self-similarity of the fractal, resulting in greater diversity and realism. In order to solve the problems inherent to plant image processing, Shlyakher et al. [START_REF] Shlyakhter | Reconstructing 3D tree models from instrumented photographs[END_REF] used image segmentation and image-based three-dimensional reconstruction techniques to produce the skeleton (trunk and branch) of the plant and L systems [START_REF] Lindenmayer | Mathematical models for cellular interactions in development I. Filaments with one-sided inputs[END_REF] to generate the 3D model of the plant. Based on image rendering technology (IBR, Image-Based Rendering), the rapid visual modeling of plant images was well realized. However, the reconstruction of skeletons mainly depended on the approximate 3D shape of the plant, and the result was different from the actual skeletons of the plants.

Okabe et al. [START_REF] Okabe | Interactive design of botanical trees using freehand sketches and example-based editing[END_REF][START_REF] Ijiri | Floral diagrams and inflorescences: interactive flower modeling using botanical structural constraints[END_REF] designed and implemented a 3D model system using hand-drawn sketches and instance-based editing operations that could be used to generate 3D tree graphics quickly and easily based on the 2D tree sketches, which were also manually drawn by the user.

A number of authors have been working on the simulation of swaying of trees. Stam et al. [START_REF] Jos | Stochastic Dynamics: Simulating the effects of turbulence on flexible structures[END_REF] obtained the turbulence effect by filtering the white noise in the frequency domain space, and then used the finite-element model analysis method to calculate the natural frequency of the branches, to obtain the vibration effect of the branches under wind load. This method simulated various types of trees, but using the white noise filtering simulation model made it difficult to represent the natural wind characteristics faithfully. Shin et al. [START_REF] Ota | A hybrid method for real-time animation of trees swaying in wind fields[END_REF] simulated the wind field by using the 1/fβ noise function. The spring model was used to simulate branch motion, while the leaves were regarded as rigid bodies. The motion was the vector sum of translation and rotation. With this method, the wind field was simulated relatively faithfully, but there was not enough space for the spring model to move with the rigid body, rendering the simulated effect far from elaborated. Akagi et al. [START_REF] Akagi | Computer animation of swaying trees based on physical simulation[END_REF] developed a wind model using the Navier-Stokes equation, and used the "boundary condition map" method to analyze the occlusion of the foliage. Hu et al. [START_REF] Hu | Data-driven modeling and animation of outdoor trees through interactive approach[END_REF] presented an interactive approach for the modelling and animation of outdoor trees from photographs and videos. The fast Fourier transform was used in that study to estimate physical parameters of branches and leaves, which were then applied to a simplified physics-based model to generate animations of trees. In the simulation of a large-scale scene, Fan et al. [START_REF] Fan | The study of real-time animation of forest scene in wind projection[END_REF] used empirical equations from meteorology instead of physical equations to accelerate the computing process. Different kinds of trees could be simulated by changing coefficients in the equations. Ancelin et al. [START_REF] Ancelin | A population approach to study forest stand stability to wind: individual tree-based mechanical modeling[END_REF] introduced a model of wind damage using the Transfer Matrix Method of beams. This model could be used to explore how damage is related to wind speed. For the uncertainty about the appropriate exponent and reliability of the conventional form, Cullen [START_REF] Cullen | Trees and wind: A practical consideration of the drag equation velocity exponent for urban tree risk management[END_REF] reported on modeling of aerodynamic drag equation forms and concluded that conventional form could be appropriate for trees. Gaffrey et al. [START_REF] Gaffrey | The elasto-mechanical behaviour of Douglas fir, its sensitivity to tree-specific properties, wind and snow loads, and implications for stability -a simulation study[END_REF][START_REF] Gaffrey | Modifying the elastomechanics of the stem and the crown needle mass distribution to affect the diameter increment distribution: A field experiment on 20-year old Abies grandis trees[END_REF] developed a complete 3D model to simulate the elastic mechanical behavior of trees under wind and gravity. The stress and strain distribution at the surface of stem could be estimated using this model. Sellier et al. [START_REF] Sellier | A numerical model of tree aerodynamic response to a turbulent airflow[END_REF] presented a dynamic model to analyze the mechanical response of trees in turbulence. The 3D description of tree architecture was integrated in this finite-element model.

In this study, we proposed a new physics-based approach for calculating deformation of trees subject to different ranges of wind loads. Furthermore, a method to develop a 3D scene within the simulation was introduced, as well as a modeling method for different types of trees. The deformation algorithms proposed were implemented and realized together with the modeling approach in a landscaping context for validation.

Materials and methods

Simulation of the scene

The 3D virtual scene contains three modules: a sky model, a terrain model, and a radiation model.

For the sky model, SkyBox [START_REF] Ping | Research and realization of modeling method for virtual geographic scenes[END_REF] was used to construct the background module of the 3D virtual scene. Representing the 3D scenes as a cube, SkyBox does not have the ability to bend texture coordinates. Thus the standard UV coordinate mapping method was applied to handle the corresponding texture, together with Cubic texture technology which was applied to simulate the details of the sky. Unlike normal textures, the six sides (i.e. Back, Front, Up, Down, Left, and Right) of the skybox can be tiled into the entire space using cube texture technology. The skybox texture is generally produced by a panorama, making the whole scene look continuous and realistic.

A rectangle was used to represent the ground to improve fusion with the Terrain model into the background. Then texture-mapping technology was used to map the ground texture image to the Terrain model to enhance simulation of the entire 3D scene.

The virtual sun was simulated as a directional light source. The direction of the incident light was determined by the location of the sun, i.e. its altitude and azimuth [START_REF] Kolivand | ReVitAge: Realistic virtual heritage taking shadows and sky illumination into account[END_REF] .

Plant modeling

The L-system formalism was used to construct the tree model.

A ternary formula can be used to represent the syntax of an L-system [START_REF] Prusinkiewicz | The algorithmic beauty of plants[END_REF] , as follows: G = < V, ω, P > where, V is the alphabet; ω is the axiom or start word, i.e. a symbol or string from V + (the set of all not empty words that can be generated with the alphabet) which represents the initial state of the tree model; P represents the set of growth or production rules of trees, as in the following example:

V : {A, B, C} ω : A P 1 : A  BC P 2 : B  AC P 3 : C  B There are three symbols in the alphabet V. The object A is a start word. Three rules (P 1 , P 2 and P 3 ) are defined in this L-system rule set: A will be replaced by BC, B be replaced by AC, and C be replaced by B. The rule set will be applied repeatedly and in parallel to a series of strings (ω -σ 1 -σ 2 -…-σ n ), which can be further interpreted with Turtle geometry [START_REF] Abelson | Turtle geometry[END_REF] , thereby rendering a developmental series of a virtual tree.

Surface modeling

To enhance the realism of the tree model, Bezier surfaces [START_REF] Hu | Approximate degree reduction of triangular bezier surfaces[END_REF] were used to simulate the trunk surface. The Bezier surface was defined as: ,,

mn ij i m j n ij P u v P B u B v    , [0,1] uv 00 ( , ) ( ) ( )
i i m i i m m B u C u u   (1) , ( ) (1 )
j j n j j n n B v C v v   (2) , ( ) (1 )
where, P ij is the space point range; u and v are the parameters of the Bezier surface. Equations ( 2) and ( 3) are the Bernstein basis functions. The space grid (also called feature grid) is formed by connecting the two adjacent points in P ij . The matrix of the Bezier surface is calculated with Equation (4), as follows:

  0, 00 01 0 1, 10 11 1 0, 1, , , 01 () () ( , ) ( ), ( ), , ( ) ()
                                (4) 2.

Topology modeling

Ranks were used to specify the branches on an individual tree model (as shown in Figure 1). A frustum of a cone was used as the basic unit for tree branches. The branch unit model was defined with the following structure, B(r t , r b , l, r, o) (5) where, r t represents the top radius of the branch; r b is the bottom radius of the branch; l indicates the length of the branch; o denotes the sequence of the branch; r indicates the position of the branches with the same order. Since the growth direction of each branch is not necessarily the same from one branch to the next, parameters to control growth directions for the branches were needed. Three turtle commands, named RU(a1), RH(a2), RL(a3) applying geometry interpretation based on Turtle geometry [START_REF] Prusinkiewicz | Graphical applications of L-systems[END_REF] , were used to accomplish this task.

Branch modeling

The tree model established with L-systems has a modular structure that is recursively generated from simple structures according to the growth and development rules of the tree in question. Taking advantage of this, calculation and realization of deformation for the branches under wind load can be made with low cost by an algorithm representing the relationship between wind forces and branch bending. In this algorithm, a single branch was divided into multiple segments that were connected with several joints. Each segment was considered as a rigid body. Therefore, the bending of the branches can be obtained by calculating the deflection angle of the nodes. A schematic diagram of the branch model is shown in Figure 2.

Wind

The wind field was simulated as a uniform flow field. Movement of deformable solid with the effect of fluid loads could exert an effect on the uniform flow field to change the distribution of fluid loads. This phenomenon was named fluid-structure interaction [START_REF] Pei | Fluid-structure coupling effect on periodically transient flow of a single-blade sewage centrifugal pump[END_REF] . Interactions between wind load and trees were considered as a similar phenomenon. Thus, a similar approach [START_REF] Li | Realistic Simulation of three-dimensional trees swaying in the wind[END_REF] could be used to simulate the effect on trees from different wind directions and speeds. Two modules, branch and leaf, were designed in this study, both were considered within the calculation of interaction between tree modules and wind load.

The virtual 3D space was divided into many voxels to calculate the influence of tree modules on the wind field. Assuming the entering wind velocity is V, then it would decay to V i after having travelled through i continuous voxels (Figure 3).

Occlusion of branches

Only one single branch segment was considered per voxel (Figure 4a). The particles of the fluid around a branch, which latter was regarded as a cylinder, were disturbed. The input wind velocity V of this voxel was decomposed into two components (V s being parallel to the direction of that branch, V n perpendicular to the direction of that branch). Since V s was not obstructed by the branches, the component V s1 of the output wind velocity in this direction would be equal to the component V s of the input wind velocity in this direction. The output wind velocity V 1 could be obtained after calculating the changes of V n . Calculation of V n can be transformed into calculation of a cylinder with circular cross-section (Figure 4b). Ideal fluid around the cylinder is composed of uniform flow and doublet flow. The complex potentials of the uniform flow and doublet flow (W 1 and W 2 , respectively) could be calculated with Equations ( 6) and ( 7):

1 0 2 (cos sin ) () (cos sin ) 2 n W V r i rr M Wi r            (6) 2 0 2 n M r V   (7)
where, M is the dipole moment, and r represents the distance between P(x, y) and the origin O. The radius of the cylinder is r 0 . θ is the angle between OP and the X axis. According to the principle of superposition of potential flows, the complex potential W of a new flow formed by the superposition of uniform flow and doublet flow could be obtained as follows:

W = W 1 + W 2 (8)
Then the velocity potential function Φ and the stream function ψ can be calculated as:

               (9)
The velocity components at position P in the flow field are:

x y V x V y                (10)
V n1 is obtained from Equation ( 10) and then combined with V s1 to get the output wind velocity V 1 .

Occlusion of leaves

Compared to the rigidity of the branches, the rigidity of leaves is low and the variation in the wind is rather complicated. In this study, a simplified method was used to approximately simulate the effect of the blade on the wind.

As shown in Figure 5, the two-dimensional projection of the blade in the direction of the wind was calculated first; then the input wind velocity V of the voxel interacts with the blade. The output wind velocity V 1 was calculated as follows:

1 1 sin sin , 0 1 ll vv SS V V k V k SS          (11)
S l is the area of the leaf in the voxel in question. S v represents the projected area of the voxel in the direction of the wind velocity.

Application

Depending on the accuracy of voxel partitioning, the number of branches and leaves in a voxel were different. In actual calculation, there might be multiple branches and leaves in one voxel. The branches were discretized (Figure 6), i.e. each branch segment was considered as a node and then could be analyzed using the above-mentioned circular cylinder model. Similarly, the leaves were also discretized into nodes. Assuming that the number of branch segments is m and the number of leaves is n, the output wind velocity V 1 will be:

1 0 0 0 0 1 sin sin n m n li bi li ii i i i v v v m bi i v S S S V V k V S S S S V S                   (12)
where, S b is the projected area of branch segment in the direction of wind velocity. V′ represents the average wind velocity. Deformations can be observed on branches when the wind blows through the trees. In doing so, the deflection angle was found to be the key parameter during simulation of the processes.

Factors affecting the deflection angle of the branches were mainly of three aspects, the wind velocity, branch diameter, and branch rigidity. The rigidity of the branch is in the Young's modulus: the higher it is, the less deformation the branch would be subjected to. As a complex system, the values of Young's modulus for different tree parts vary. To simplify it in the simulation, the branches of a tree were divided into three parts: trunk, branch, and twig. The range of values of Young's modulus for the parts are given in Table 1 (according to US Department of Agriculture [START_REF] Concepts | Wood handbook -wood as an engineering material[END_REF] ).

Table 1 The value range of Young's modulus for a tree [START_REF] Concepts | Wood handbook -wood as an engineering material[END_REF] Different parts of the tree Young's modulus value/GPa trunk 5-12 branch 1-5 twig 0.5-1

Force analysis of branches

As shown in Figure 7, the branch segment was connected with its father branch at point A, and there might be a child branch at the end B. The branch segment was affected by the wind force F 1 and the force F 2 exerted by the rear segments at the end B if there was a child branch attached on end B. The point of action of the wind force F 1 is at the center of mass of the branch segment. A moment of couple m is attached to the branch according to the force translation theorem to move F 2 to O. m = F 2 ×L BO (13) where, L BO is the center of mass of the branch segment relative to the end point B of the branch segment. The relationship between wind velocity and wind force can be calculated with Equation ( 14). where, C d represents the resistance coefficient; C t , the lift coefficient; ρ, the air density; V, the velocity of wind; F d , the resistance of the branch, and F t the drag force of the branch.

Figure 7 Force analysis for a branch segment. The force F 2 was exerted by the rear segments at the B end in the case of a child branch attached to B

Calculation of deformation

The end of a branch which attached to its father branch was considered as fixed on the upper branch and the other end was free to move. Thus, deformation of the cantilever beam from the mechanics of materials can be applied on this matter. Since uniformly distributed loads were acting upon the beam, the bending moment on the full length of the beam could be described by a function. As shown in Figure 8, it was cut off from any cross section, with a uniformly distributed load permanently present on the left side of the section, and a reactive force F RO on the fulcrum generated by the wind load q, but only a uniformly distributed load existed on the right side of the section. The balance of all the forces on the right part was used to calculate the bending moment, as follows:

2 1 ( ) ( ) 2 M x q L x    (15)
where, q is uniformly distributed load; L is the length of the beam; M represents the bending moment. The deflection differential equation is:

2 2 () () d W x EI M x dx  (16)
where, E represents the Young's modulus of the branch, and I is the area moment of inertia. Equation (15) and Equation ( 16) are integrated in the following equation:

3 ( ) 1 () 6 dW x EI q L x C dx     (17)
where, C is an arbitrary constant. The constraints at the fixed end are as follows: 0 () 0

x dW x dx   (18)
Then C is calculated as Equation ( 19):

3 6 qL C  (19)
Combining ( 17) with ( 19), the following equation can be obtained:

33 () [ ( ) ] 6 dW x q L L x dx EI    (20)
Finally, the deflection angle can be calculated as:

() xL dW x dx    (21)
Figure 8 Diagram showing the stress analysis of cantilever beam.

Uniformly distributed loads q are applied to the beam. F RO is a reactive force on the fulcrum generated by the wind load q 3 Simulation Results

Using the method described above, a dynamic simulation system with 3D tree model was implemented in the Eclipse development environment using Java programming language and JOGL graphics engine. The computer hardware used in the simulation was configured with Intel (R) Core (TM) i7-6700HQ CPU, 8GB memory, and NVIDIA GeForce GTX 970M graphics card.

Since the entire 3D space was divided into multiple voxels, the division accuracy would have a direct influence on the simulation results. With higher division accuracy or more voxels designed, the results would be more accurate, but also be more time consuming. Therefore, it was necessary to find a division model with a relatively high computational efficiency and which at the same time would meet the requirement of simulation accuracy. For a tree model, the total number of discretization nodes was set to a fixed value, which meant the number of nodes within a certain voxel would be different when the total voxel number was changed. For a tree model with 23,684 triangles, Table 2 shows the relationship between the division accuracy and the calculation time for each frame. It can be clearly observed that the calculation took much longer when the space was divided into more voxels, among which the 80×80×80 voxels took the longest time to calculate, while the 40×40×40 only took about one third of the former, which, compared to those lower division settings, meant that it gained enough accuracy but with acceptable calculation time. The simulation results were verified on three aspects: 1) movement of a single branch under different wind velocities; 2) movement of the same tree under different wind velocities; 3) movement of different trees under the same wind velocity.

A branch from Koelreuteria paniculata was used for simulation to demonstrate the results of the algorithm described for a single branch, of which the bending was affected with four settings of wind velocity values, i.e. 0 m/s, 4 m/s, 8 m/s, and 12 m/s. The Young's modulus of the branch is 3.4 GPa, and the maximum curvatures of the branch bending under different wind velocities are shown in Figure 9. The dynamics of two tree species (Koelreuteria paniculata and Ulmus pumila) under four different wind loads, i.e. 0 m/s, 4 m/s, 8 m/s and 12 m/s, were simulated as the second step. The maximum curvatures were captured and compared (as shown in Figure 10). Different dynamics of bending and swaying of the branches between two trees and among different wind velocities on the same tree could be observed. Physical properties of the two trees are different, including tree architecture, Young's modulus, etc. This is the main reason why the two trees reacted differently under the same wind load. Concerning the same tree under different wind velocities, it was clear that it became more bent as the wind velocity increased. Finally, the method proposed by Akagi and Kitajima [START_REF] Akagi | Computer animation of swaying trees based on physical simulation[END_REF] was applied to the same tree models, as shown in Figure 11, and the result shows that there was no significant difference between dynamic reactions from their method and the algorithm proposed in this study. Figure 11 Snapshots of the maximum bending state of the two tree models, Koelreuteria paniculata (a1-d1) and Ulmus pumila (a2-d2), under four different wind loads with algorithm discussed in the paper [START_REF] Akagi | Computer animation of swaying trees based on physical simulation[END_REF] . Velocity of wind: a1, a2: 0 m/s; b1, b2: 4 m/s; c1, c2: 8 m/s; d1, d2: 12 m/s

Discussion

A physics-based algorithm to simulate tree dynamics of two different tree species under various wind loads was proposed in this study. By considering the physical properties of the trees, i.e. branches and leaves, the wind field and forces on trees can be simulated. Leaf shedding under high wind load has not yet been considered in our model, which should be considered in the simulation in further study.

One of the applications of the algorithm lies in tree cultivation and forest industry. In the forest, the strong wind can affect the growth process of the trees, resulting in structural deformation of the trees, even affecting the structural layout of the forests and the tree species diversity [START_REF] Dupont | Wind damage propagation in forests[END_REF] . In urban areas, trees on both sides of a road are always at risk to be ruptured by strong winds, often causing harm and damage to people, vehicles or houses [START_REF] Kennethr | Mechanical stability of trees under dynamic loads[END_REF] , all this representing substantial financial risks for insurance companies. In a forestry context, the simulation algorithm can used to predict storm damage or to optimize planting density and row orientation in a stand in certain geographical locations, and it can also be used to simulate and quantify the financial risks for urban landscaping. Currently, only the morphological effects of the trees under wind load were considered in the algorithm, while the impacts on physiological processes of trees were not. It is also essential to combine the effects of the wind load on biomechanics of the trees in order to better simulate its dynamics. We will do this in the near future.

Another application of our model could be in the context of ecological modelling, more particularly bioclimatology and photosynthesis models: it is well known that the rate of photosynthesis depends on the concentration of CO 2 inside the leaf, which in turn depends on stomatal conductance [START_REF] Nikolov | Coupling biochemical and biophysical processes at the leaf level: an equilibrium photosynthesis model for leaves of C3 plants[END_REF] . The latter process is controlled by a number of plant and biophysical processes, amongst which features wind velocity. Detailed information (at a fine spatial resolution) about the speed and direction of air movement around the leaf would greatly improve estimations of the CO 2 concentration and relative humidity in the leaf boundary layer, which latter serves as the base for estimating the amount of CO 2 available for photosynthesis, and which currently cannot be very well quantified.

Conclusions

A method for simulating the movement of trees in the wind based on physics was presented in the study. The L-system formalism was used to construct a tree model and a whole branch was divided into multiple segments. Considering the wind load acting on each branch, we calculated the deflection of each branch of the tree in the wind, and then superimposed the deflection of all the small segments to obtain the movement of a whole branch. The experimental results show that the method is more realistic and greatly reduces the computation time.

It is true that L-systems are not necessarily the only method to simulate trees with static architecture exposed to wind loads, as other methods of architecture reconstruction would be more efficient in some way. L-systems were used in this study, as it is only the first step to consider the static structure of the simulated trees. Our next step will include the effect of wind load on growth of a tree, and eventually consider the influence on the whole growth period and not only on structure, but also on functional physiological processes.

Admittedly, distribution of the wind field in a 3D space is relatively complicated, not to mention movement of a tree in the wind. To simplify the question, we assumed that the power of the wind within the space stayed the same, in an attempt to control calculation cost, however, simulation results were also less convincing. In order to get more realistic results, analysis of spatial distribution of wind field will be added in the future.

Figure 1

 1 Figure 1 Structure of branches with different ranks designed within the topological models

Figure 2

 2 Figure 2 Schematic diagram showing the structure of a branch and its reaction to the wind force in the model. Points refer to joints connecting different segments, while arrows represent the wind force, which has an impact on the joint to bend the branch

Figure 3

 3 Figure 3 Diagram of squares each represents a voxel. The wind velocity V decays to V 1 after having crossed the first voxel. As before, V 1 decays to V 2 after having crossed the second voxel 2.3.1 Occlusion of branchesOnly one single branch segment was considered per voxel (Figure4a). The particles of the fluid around a branch, which latter was regarded as a cylinder, were disturbed. The input wind velocity V of this voxel was decomposed into two components (V s being parallel to the direction of that branch, V n perpendicular to the direction of that branch). Since V s was not obstructed by the branches, the component V s1 of the output wind velocity in this direction would be equal to the component V s of the input wind

Figure 4

 4 Figure 4 Diagram of a voxel with branch segment and wind load upon the stem

Figure 5

 5 Figure 5 Schematic diagram of occlusion of leaves. The rectangle in the middle represents the projection of a leaf. θ is the angle between the direction of leaf projection and input wind velocity 2.3.3 ApplicationDepending on the accuracy of voxel partitioning, the number of branches and leaves in a voxel were different. In actual calculation, there might be multiple branches and leaves in one voxel. The branches were discretized (Figure6), i.e. each branch segment was considered as a node and then could be analyzed using the above-mentioned circular cylinder model. Similarly, the leaves were also discretized into nodes. Assuming that the number of branch segments is m and the number of leaves is n, the output wind velocity V 1 will be:

Figure 6 (

 6 Figure 6 (a) Branch segments and leaves within a voxel; (b) Discretized nodes within a voxel, a green node indicates the leaf while a black node indicates the branch 2.4 Dynamic simulationDeformations can be observed on branches when the wind blows through the trees. In doing so, the deflection angle was found to be the key parameter during simulation of the processes.Factors affecting the deflection angle of the branches were mainly of three aspects, the wind velocity, branch diameter, and branch rigidity. The rigidity of the branch is in the Young's modulus: the higher it is, the less deformation the branch would be subjected to. As a complex system, the values of Young's modulus for different tree parts vary. To simplify it in the simulation, the branches of a tree were divided into three parts: trunk, branch, and twig. The range of values of Young's modulus for the parts are given in Table1(according to US Department of Agriculture[START_REF] Concepts | Wood handbook -wood as an engineering material[END_REF] ).

Figure 9

 9 Figure 9 Maximum curvatures of the branch (Koelreuteria paniculata) bending under different wind velocities

Figure 10

 10 Figure 10 Snapshots of the maximum bending state of the two tree models, Koelreuteria paniculata (a1-d1) and Ulmus pumila (a2-d2), under four different wind loads with velocities of: a1, a2: 0 m/s; b1, b2: 4 m/s; c1, c2: 8 m/s; d1, d2: 12 m/s, with the physics-based algorithm proposed in this study

Table 2 Calculation time cost within the algorithm under different division accuracies

 2

	Division accuracy	Calculation time/ms
	5×5×5	5.6253
	10×10×10	8.3641
	20×20×20	13.8216
	40×40×40	34.2537
	80×80×80	113.2165

Open Access at https://www.ijabe.org Vol. 13 No. 2

Acknowledgements

This work was supported by the National Natural Science Foundation of China (Grant No. 31301230, 61571400, 31471416); and the Scientific Research Fund of Zhejiang Provincial Education Department (Grant No. Y201224867). The authors are grateful to the anonymous reviewers whose comments helped to improve this paper.