François Jouanne

Ananta Gajurel

Jean-Louis Mugnier

Lok Laurent Bollinger

Bijaya Adhikari

Bharat Koirala

Nathalie Cotte

Roshanraj Bhattarai

Pascale Bascou

Pascale Huyghe

Bijaya Lok

Adhikari

Arnaud Pecher

Postseismic deformation following the April 25, 2015 Gorkha earthquake (Nepal): Afterslip versus viscous relaxation

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

southern part of the Main Himalayan Thrust (MHT). Three hypotheses for the mechanisms controlling postseismic deformation are tested through numerical simulations of the postseismic time series: (i) viscous relaxation, (ii) afterslip, or (iii) a combination of these two mechanisms. We can exclude postseismic deformation controlled by viscous relaxation of a thick deformation zone along the northern and lower flat of the MHT. However, it is impossible to discriminate between postseismic deformation controlled by either afterslip along the MHT (northern part of the rupture zone, crustal ramp, and lower flat of the MHT) or a combination of afterslip along the MHT (northern part of the rupture zone, crustal ramp) and viscous relaxation controlled by a thin (~3-4 km thick) low-viscosity body centered on the lower flat of the MHT. The occurrence of afterslip along the northern part of the upper flat of the MHT and its longitudinal variations have been established thanks to the densification of GNSS network by our team presented in this paper.

Introduction.

The large number of geophysical instruments which recorded the April 25, 2015 Mw 7.8 Gorkha earthquake provides a unique opportunity to document better than ever a large Himalayan earthquake. The main instrumental contributors to our understanding of the mainshock include Interferometric synthetic aperture radar (InSAR) investigations as well as teleseismic data (e.g. Avouac et al., 2015, Fan andShearer, 2015), with also the continuous Global Navigation Satellite System (GNSS) time series [START_REF] Grandin | Rupture process of the Mw = 7.9 2015 Gorkha earthquake[END_REF] and regional seismological records (e.g. [START_REF] Adhikari | The aftershock sequence of the April 25 2015 Gorkha-Nepal earthquake[END_REF][START_REF] Baillard | Automatic analysis of the Gorkha earthquake aftershock sequence: evidences of structurally segmented seismicity[END_REF]. The deployment of the Nepal GPS Geodetic Network, prior to the earthquake, and its maintenance during the seismic crisis has also allowed to quantifying the early postseismic surface deformation (KKN4 and CHLM stations) consecutive to this large earthquake [START_REF] Mencin | Himalayan strain reservoir inferred from limited afterslip following the Gorkha earthquake[END_REF][START_REF] Gualandi | Pre-and postseismic deformation related to the 2015, Mw 7.8 Gorkha earthquake, Nepal[END_REF].

In addition to being one of the few instrumental Himalayan earthquake, the setting of the 2015 earthquake (Fig. 1&2) is appropriate to study the complexity of the seismic cycle in a continental subduction zone. Indeed, the slip deficit along the main Himalayan thrust (MHT), on which the mainshock occurred, is accommodated during various types of earthquakes [START_REF] Mugnier | Structural interpretation of the great earthquakes of the last millennium in Central Himalaya[END_REF], that partially contributes to release the elastic strain accumulation [START_REF] Avouac | Lower edge of locked Main Himalayan Thrust unzipped by the 2015 Gorkha earthquake[END_REF]. This partial release has already been instrumentally measured in several oceanic subduction environments; for example, in the Andean convergent margin, the main slip patch of the 2010 Mw 8.8 earthquake occurred in an area that was highly coupled and had already released slip in 1960 [START_REF] Melnick | Darwin seismic gap closed by the 2010 Maule earthquake[END_REF][START_REF] Moreno | Toward understanding tectonic control on the Mw 8.8 2010 Maule Chile earthquake[END_REF]. However, contrarily to earthquakes ruptures at oceanic subduction zones rarely captured by GNSS network in their whole extend due to their frequent extension offshore, the 2015 Gorkha rupture was surrounded by 12 stations well distributed above and around it. In addition, this permanent network was complemented during the weeks following the main shock, by the deployment of additional permanent or semi-permanent GNSS stations (this paper), this network dedicated to the postseismic deformation has been designed to test the occurrence of afterslip along the upper flat of the MHT in the rupture zone and south of it.

These records add some details to the postseismic deformation pattern and give a unique opportunity to study in greater details the mechanism controlling the postseismic deformation associated with the 2015 Gorkha event.

Different mechanisms are usually inferred for the postseismic deformation: afterslip or viscous relaxation. In this paper, the Relax software [START_REF] Barbot | Postseismic deformation due to the Mw 6.0 2004 Parkfield earthquake: stress-driven creep on a fault with spatially variable rate-and-state friction parameters[END_REF] was used to model the postseismic deformation, taking into account the Main Himalayan Thrust geometry, rather well imaged in the region (Nabelek et al., 2009;[START_REF] Duputel | The 2015 Gorkha earthquake: a large event illuminating the Main Himalayan Thrust fault[END_REF] and the coseismic displacement field associated to it. Three postseismic mechanisms were tested: afterslip, viscous relaxation or a combination of these two phenomena. We used a flap-ramp-flat geometry [START_REF] Elliott | Himalayan megathrust geometry and relation to topography revealed by the Gorkha earthquake[END_REF] with a southern and upper flat corresponding to the geometry adopted by [START_REF] Grandin | Rupture process of the Mw = 7.9 2015 Gorkha earthquake[END_REF] to describe the coseismic slip distribution of the Gorkha earthquake along this flat.

Geological setting.

The structure of the Himalaya results from the underthrusting of the Indian lithosphere along the Main Himalayan Thrust (MHT) beneath the Tibetan Plateau [START_REF] Argand | La tectonique de l'Asie. 104 p. Compte rendu du congrès géologique international[END_REF]. Great earthquakes have episodically ruptured segments of the brittle upper part of the MHT (e.g. [START_REF] Avouac | Le cycle sismique en Himalaya[END_REF][START_REF] Mugnier | Structural interpretation of the great earthquakes of the last millennium in Central Himalaya[END_REF][START_REF] Bollinger | Estimating the return times of great Himalayan earthquakes in eastern Nepal: Evidence from the Patu and Bardibas strands of the Main Frontal Thrust[END_REF]. The Nepal earthquake of April 25, 2015 followed a series of great earthquakes in the Central Himalaya [START_REF] Mugnier | Seismites in the Kathmandu valley and seismic hazard in central Himalaya[END_REF], Bollinger et al., 2016): over the last two centuries, the city of Kathmandu was severely damaged in 1833, 1866, and 1934. The April 25, 2015 event (Fig. 1) was the first earthquake on the MHT to be simultaneously recorded by numerous instruments: high-rate GPS (e.g. [START_REF] Avouac | Lower edge of locked Main Himalayan Thrust unzipped by the 2015 Gorkha earthquake[END_REF], teleseismic waves (e.g. [START_REF] Fan | Detailed rupture imaging of the 25 April 2015 Nepal earthquake using teleseismic P waves[END_REF], SAR imaging (e.g. [START_REF] Lindsey | Line-of-sight displacement from ALOS-2 interferometry: Mw 7.8 Gorkha Earthquake and Mw 7.3 aftershock[END_REF], InSAR and teleseismic waves [START_REF] Grandin | Rupture process of the Mw = 7.9 2015 Gorkha earthquake[END_REF], strong-motion recordings (e.g. [START_REF] Bhattarai | Overview of the large 25 April 2015 Gorkha, Nepal, earthquake from accelerometric perspectives[END_REF], and a local seismometer network [START_REF] Adhikari | The aftershock sequence of the April 25 2015 Gorkha-Nepal earthquake[END_REF].

The rupture was mainly located to the NW of Kathmandu, at a depth of 13-15 km on a flat portion of the Main Himalayan Thrust (MHT) that dips N-NE at 7-10°. The northern boundary of the main rupture corresponds to the transition towards a steeper crustal ramp. This ramp, which is partly coupled during the interseismic period [START_REF] Jouanne | Estimation of coupling along the Main Himalayan Thrust in the central Himalaya[END_REF], was only locally affected by the earthquake [START_REF] Elliott | Himalayan megathrust geometry and relation to topography revealed by the Gorkha earthquake[END_REF]. The southern boundary of the rupture was near the leading edge of the Lesser Himalaya antiformal duplex [START_REF] Mugnier | Segmentation of the Himalayan megathrust around the Gorkha earthquake (25 April 2015) in Nepal[END_REF] and possibly near the frontal footwall ramp of the upper Nawakot duplex [START_REF] Hubbard | Structural segmentation controlled the 2015 Mw 7.8 Gorkha earthquake rupture in Nepal[END_REF]. This rupture was affected by transversal structures: on the western side, a ~20° dipping lateral ramp (Kumar et al., 2017), located beneath the Judi lineament [START_REF] Mugnier | Segmentation of the Himalayan megathrust around the Gorkha earthquake (25 April 2015) in Nepal[END_REF] and already defined by [START_REF] Kayal | Microearthquake Seismology and Seismotectonics of South Asia[END_REF], separates the main rupture zone from the nucleation area. On the eastern side, a lateral ramp [START_REF] Mugnier | Segmentation of the Himalayan megathrust around the Gorkha earthquake (25 April 2015) in Nepal[END_REF] located beneath the Gaurishankar lineament [START_REF] Kayal | Microearthquake Seismology and Seismotectonics of South Asia[END_REF] illustrated by its seismic signature at depth [START_REF] Jean | Teleseismic depth estimation of the 2015 Gorkha-Nepal aftershocks[END_REF][START_REF] Baillard | Automatic analysis of the Gorkha earthquake aftershock sequence: evidences of structurally segmented seismicity[END_REF] separates the April 25, 2015 rupture from the May 12, 2015 (Mw 7.2) rupture.

Data acquisition and analysis.

After the April 25 Gorkha earthquake, we installed a GNSS network designed to record postseismic deformation induced by this Mw 7.9 earthquake. We installed four permanent GNSS stations in June 2015 and 11 campaign GNSS stations between June and December 2015.

Seven of these stations were measured again in November 2016 (Figs 3 and4). We analyzed these data together with the continuous GNSS networks established by LDG-DASE, Caltech (www.unavco.org), Central Washington University (http://www.geodesy.cwu.edu) and with data from the GNSS sites defined in the ITRF2014 reference frame [START_REF] Altamimi | ITRF2014: A new release of the International Terrestrial Reference Frame modeling nonlinear station motions[END_REF] (ARTU, BADG,BAKU,BJFS, BJNM, CHUM, CUSV, GUAO, HYDE, IISC, IRKJ, IRKM, KIT3, KUWT, LHAZ, MDVJ, POL2, SGOC, SHAO, TALA, TASH, TCMS, TEHN, TNML, URUM, ZECK) obtained from unavco.

The data were analyzed using Bernese 5.2 software (Dach et al., 2016) with absolute antenna phase center offset models, together with precise IGS orbits, Earth rotation parameters, and ocean tidal loading and atmospheric tidal loading estimates. Velocities and time series were estimated in the ITRF2014 reference frame [START_REF] Altamimi | ITRF2014: A new release of the International Terrestrial Reference Frame modeling nonlinear station motions[END_REF] with discontinuities associated with this reference frame and expressed in terms of the India fixed reference frame by the use of the rotation pole proposed by [START_REF] Ader | Convergence rate across the Nepal Himalaya and interseismic coupling on the Main Himalayan Thrust: implications[END_REF]. This rotation pole, originally proposed to describe the motion of India relative to ITRF2005, accurately describes the India plate motion in the ITRF2014 reference frame, as shown by [START_REF] Jouanne | Estimation of coupling along the Main Himalayan Thrust in the central Himalaya[END_REF] and as exemplified by the lack of velocity of the cGPS SIM4 station located on the India plate (Fig. 3). In any case, any velocity changes introduced by the ITRF2005 / ITRF2014 change are negligible compared to the large postseismic displacements.

We followed the resolution strategy with (1) an initial ionosphere-free analysis with calculation of the residuals; (2) a residual analysis; (3) code-based wide-lane ambiguity resolution for all baselines [START_REF] Melbourne | The case for ranging in GPS based geodetic systems[END_REF][START_REF] Wübbena | Software developments for geodetic positioning with GPS using TI4100 code and carrier measurements[END_REF], using differential code bias (DCB) files when available and calculation of the ionosphere-free solution with the introduction of resolved Melbourne-Wübbena linear combination ambiguities; (4) phase-based wide-lane (L5) ambiguity resolution for baselines < 200 km and computation of the ionosphere-free solution with the introduction of resolved ambiguities; (5) resolution of the previously unresolved ambiguities for baselines < 2000 km using the quasi ionosphere-free strategy of resolution; (6) direct L1/L2 ambiguity resolution for baselines < 20 km with the introduction of an ionosphere model; (7) calculation of the normal equations; (8) a compatibility test between the daily free solution and ITRF2014 solution, selection of compatible ITRF2014 stations, and (9) transformation of the daily normal equation in the ITRF2014 reference frame with a sixparameter Helmert solution (three translation parameters and three rotation parameters) using the ITRF2014 selected stations. During these steps, site-specific troposphere parameters were estimated every two hours.

Normal equations were analyzed simultaneously to determine accurate velocities in the ITRF2014 reference frame with the introduction of ITRF2014 coordinates and velocities.

Outliers and new discontinuities were detected using the "Find Outliers and Discontinuities in Time Series" tool in Bernese 5.2 software which estimates annual seasonal fluctuations and reduces, step by step, the discrepancy between the functional model and the time series due to statistical adjustment [START_REF] Ostini | FODITS: a new tool of the Bernese GPS software[END_REF]. Bernese 5.2 software underestimates the daily errors because systematic errors or mismodeled parameters are not included in the formal error. We therefore rescaled the formal errors by multiplying them by a factor of 10 to obtain a more realistic estimated error [START_REF] Hugentobler | Bernese GPS Software Version 4.2[END_REF].

Results: postseismic displacements at the surface.

We consider that displacements contain three superimposed signals: interseismic displacements, postseismic displacements, and the seasonal signal associated with loading of the Ganga Plain during the monsoon and its unloading during the dry season [START_REF] Bettinelli | Seasonal variations of seismicity and geodetic strain in the Himalaya induced by surface hydrology[END_REF][START_REF] Fu | Seasonal and long-term vertical deformation in the Nepal Himalaya constrained by GPS and GRACE measurements[END_REF][START_REF] Gualandi | Pre-and postseismic deformation related to the 2015, Mw 7.8 Gorkha earthquake, Nepal[END_REF]. In order to analyze postseismic displacements, we removed the interseismic displacements that were either measured before the Gorkha earthquake for the sites that already existed before this event (BESI, BRN2, CHLM, DAMA, GUMB, KKN4, KLDN, LMJG, NAST, ODRE, RMJT, RMTE, SNDL, SYBC, TPLJ) (Fig. 2; [START_REF] Jouanne | Estimation of coupling along the Main Himalayan Thrust in the central Himalaya[END_REF] or were estimated using interseismic deformation modeling by [START_REF] Jouanne | Estimation of coupling along the Main Himalayan Thrust in the central Himalaya[END_REF] for the new GNSS sites (Fig. 3).

In order to evaluate the seasonal components, we make the hypothesis that time series are nearly linear during the period beginning 300 days after the 25 April 2015 earthquake. For the period beginning 300 days after the main shock, we have then modeled time series as the consequence of a linear velocity and of seasonal component (eq. 1).

U(t) = U*(t-tref) + a cos (2π*(t-tref)/365) + b sin (2π*(t-tref)/365) + c cos (2π*(t-tref)/182.5) + d sin (2π*(t-tref)/182.5) (eq. 1)
Where tref is the first date after the day D0+300 existing in the time series, D0 being the date of the main shock, U the velocity, and a, b, c and d are the amplitude of annual and semi-annual signals. In a second step, whole time series are corrected for seasonal signals using these amplitudes.

The postseismic displacements for the period when most of the GNSS stations were available (between June 24, 2015 and late 2016) are shown (Fig. 4) to compare site displacement measured over the same time span. Postseismic displacements are characterized by the southward displacements of the points north of the Gorkha rupture (CHLM, GUMB, TSM1). The station that recorded the greatest postseismic displacement is CHLM (Fig. The slip deficit could then possibly be resorbed by a future large earthquake between the southern end of the rupture near Kathmandu and the Main Frontal Thrust near the Ganga Plain, as probably happened in the 1866 earthquake [START_REF] Oldham | A catalogue of Indian earthquakes from the earliest time to the end of 1869 AD[END_REF][START_REF] Szeliga | Intensity, magnitude, location, and attenuation in India for felt earthquakes since 1762[END_REF], or by a larger earthquake rupturing the locked fault zone from the brittle-ductile transition zone to the surface, in a similar way as in the subduction zones [START_REF] Melnick | Darwin seismic gap closed by the 2010 Maule earthquake[END_REF][START_REF] Moreno | Toward understanding tectonic control on the Mw 8.8 2010 Maule Chile earthquake[END_REF] or to the 1934 Nepal-Bihar earthquake (e.g. [START_REF] Mugnier | Structural interpretation of the great earthquakes of the last millennium in Central Himalaya[END_REF]. Time series of southern component present a progressive decrease in postseismic displacement with time, suggesting progressive dissipation of the stress perturbation induced by the earthquake (Fig. 5). It can be observed that the time series of the eastern components (Fig. 6) presents spatial variations: BALE, CHLM and DNC4 time series indicate significant westward components of postseismic displacements, whereas time series of BHAZ, CHAR, GUMB, JIR2, KKN4 and XBAR do not present this component. This is why we will consider only the southern and vertical time series to simulate postseismic displacements. Time series of vertical component present a progressive decrease in postseismic displacement with time, suggesting progressive dissipation of the stress perturbation induced by the earthquake (Fig. 7). The CHLM time series of the south components presents a jump corresponding to the Mw 7.2 main aftershocks of May 12, whereas the other stations were not significantly affected or were installed too late to record this coseismic displacement. It can be also observed that this second earthquake does not induced a change in the evolution of the postseismic displacement of the CHLM and KKN4 station. This is why we have not consider this earthquake in the numerical simulations.

Modeling postseismic-displacement.

We have chosen to test the origin of the progressive release using three hypotheses with Relax 1.0.7 software [START_REF] Barbot | Space geodetic investigation of the coseismic and postseismic deformation due to the 2003 Mw 7.2 Altai earthquake: Implications for the local lithospheric rheology[END_REF][START_REF] Barbot | Postseismic deformation due to the Mw 6.0 2004 Parkfield earthquake: stress-driven creep on a fault with spatially variable rate-and-state friction parameters[END_REF]Barbot and Fialko 2010;[START_REF] Bruhat | Evidence for postseismic deformation of the lower crust following the[END_REF][START_REF] Rousset | Postseismic deformation following the 1999 Chi-Chi earthquake, Taiwan: Implication for lower-crust rheology[END_REF] For this modeling, we use the same MHT geometry as for the pre-earthquake deformation [START_REF] Jouanne | Estimation of coupling along the Main Himalayan Thrust in the central Himalaya[END_REF]; this geometry also fits with the location of the rupture zone proposed by [START_REF] Grandin | Rupture process of the Mw = 7.9 2015 Gorkha earthquake[END_REF] indicated in Figure 2. The geometry of the MHT is formed by a succession of flats and ramps characterized by a large southern flat (95-100 km) dipping 7° northward, a crustal ramp dipping 27° northward between 15 and 25 km, and a northern flat dipping 7° northward as illustrated by receiver-function data [START_REF] Duputel | The 2015 Gorkha earthquake: a large event illuminating the Main Himalayan Thrust fault[END_REF].

As postseismic displacements only affected stations north of Kathmandu, if postseismic displacements are controlled by afterslip along the MHT, afterslip did not affect all of the upper flat of the MHT but only its northern parts. If we consider the points BALE, KKN4, and GURJ (Fig. 4 and5) that are nearly the same distance from the upper flat-ramp connection, it appears that the eastern one (BALE) is affected by greater postseismic displacements than the other two located in the western part of the upper flat of the MHT. We then considered that these surficial variations in displacement reflect afterslip distribution along the MHT. These spatial variations in the time series led us to consider a 25 km-wide western segment and a 40 km-wide eastern segment and the boundary between these two segments coincides with the western edge of the May 12, 2015 rupture and with the Gaurishankar lateral ramp as illustrated in Figure 1. This spatial variation also reflects interseismic coupling variations between the western and eastern parts: the western part was completely locked whereas the eastern part was only partially locked, with a coupling of 0.7 (Fig. 2a) [START_REF] Jouanne | Estimation of coupling along the Main Himalayan Thrust in the central Himalaya[END_REF].

To test the models, we used a time series that was corrected for interseismic velocities and seasonal components. As postseismic displacements are southward with a lack of E-W components, we have only used the N-S and vertical components of time series.

The preferred solutions were selected using the comparison between the observed and simulated time series with the estimation of weighted root mean square WRMS for each simulation (Fig. 9).

Viscous relaxation model.

The estimation of coupling during the interseismic period (Betinelli et al., 2006;[START_REF] Ader | Convergence rate across the Nepal Himalaya and interseismic coupling on the Main Himalayan Thrust: implications[END_REF][START_REF] Jouanne | Estimation of coupling along the Main Himalayan Thrust in the central Himalaya[END_REF] suggests ductile behavior along the northern flat of the MHT.

We therefore test the hypothesis that postseismic displacements reveal viscous relaxation of the stress perturbation induced by the main shock. Ductile behavior is temperature dependent [START_REF] Scholz | Earthquakes and friction laws[END_REF] and does not occur in the brittle crust. Viscous relaxation is therefore excluded for the upper flat and the ramp characterized by a temperature of less than 375 °C [START_REF] Robert | Control of detachment geometry on lateral variations in exhumation rates in the Himalaya: insights from lowtemperature thermochronology and numerical modeling[END_REF]. Hereafter, we assume that the lower flat accommodates viscous relaxation only.

In this simulation, a low viscous body with a linear Maxwell rheology corresponding to the deeper part of the Main Himalayan Thrust is embedded in an elastic half space (Fig. 8a). The parameters to be determined are the viscosity and the thickness of this body. As the indian lithosphere is very thick, 135 km in eastern Nepal [START_REF] Uma | Imaging the Indian lithosphere beneath the Eastern Himalayan region[END_REF], we do not take into account the asthenosphere viscosity.

We chose to explore the following ranges of values: viscosity of the body from 10 16 to 10 19 Pa*s and thickness of the body from 0 to 14 km. The WRMS is invariably greater than 17 mm and no solution is clearly favored among the 84 numerical runs performed to determine parameter influence and summarized in Figure 9a.

Afterslip model.

If the postseismic displacements are controlled by afterslip induced by the Mw 7.9 earthquake of April 25 2015, these displacements can be simulated by progressive release of stress perturbation due to the main shock. We assume that the afterslip occurs (Fig. 8a) along the northern part of the upper flat of the MHT affected by the April 25, 2015 earthquake, along the crustal ramp north of the rupture zone, and along the lower flat north of this ramp.

We assign rate-strengthening properties derived from laboratory experiments [START_REF] Dieterich | Modeling of rock friction: 1. Experimental results and constitutive equations[END_REF][START_REF] Ruina | Slip instability and state variable friction laws[END_REF][START_REF] Barbot | Postseismic deformation due to the Mw 6.0 2004 Parkfield earthquake: stress-driven creep on a fault with spatially variable rate-and-state friction parameters[END_REF][START_REF] Lapusta | Elastodynamic analysis for slow tectonic loading with spontaneous rupture episodes on faults with rate-and state-dependent friction[END_REF][START_REF] Rice | Slip complexity in earthquake fault models[END_REF] to these three fault segments with:

(eq. 2)

Where V is the slip velocity on the fault, ∆τ is the shear stress perturbation induced by the earthquake that is progressively released by the afterslip, is the initial slip velocity and (ab) σ is a constitutive parameters with a value of the order of the stress drop during the earthquake, σ being the effective normal stress and (a -b) a frictional parameter.

The initial stress perturbation associated with the coseismic slip distributions of April 25, 2015 estimated by [START_REF] Grandin | Rupture process of the Mw = 7.9 2015 Gorkha earthquake[END_REF], with InSAR and teleseismic waves data (Fig. 2b) is calculated with Relax 1.0.7 software.

We explore the parameters space (Fig. 9b,c mm/year. The simulation of vertical time series allows to confirm that must be larger than 150 mm/year (Fig. 9d).

We decide to adopt as the better model the solution obtained for 0.7 frictional parameters and a value of 200 mm/year.

Afterslip and viscous relaxation model.

This simulation takes into account the hypothesis that the northern flat may be a zone affected by ductile deformation [START_REF] Beaumont | Himalayan tectonics explained by extrusion of a low-viscosity crustal channel coupled to focused surface denudation[END_REF], while the upper part of the MHT is probably in a brittle condition as shown by the high coupling of the upper flat and crustal ramp during the interseismic period [START_REF] Jouanne | Current shortening across the Himalayas of Nepal[END_REF][START_REF] Ader | Convergence rate across the Nepal Himalaya and interseismic coupling on the Main Himalayan Thrust: implications[END_REF][START_REF] Jouanne | Estimation of coupling along the Main Himalayan Thrust in the central Himalaya[END_REF].

We test the hypothesis of a mixed origin for postseismic deformation (Fig. 8), with viscous relaxation controlled by a deep low-viscosity body centered on the northern MHT flat, as done above, whereas afterslip occurs along the crustal ramp and the northern part of the southern flat.

We adopt, as an a priori model for the afterslip component, the above model that minimizes the WRMS among the models based solely on afterslip and we explore the following parameters: viscosity of the low-viscosity body (from 10 16 to 10 19 Pa.s) and the thickness of this body (1 to 6 km).

Among 50 numerical runs performed to determine the parameter influence with the southern time series (Figure 9e), a minimum of 3.2 mm is found for the WRMS (Fig. 9e) linked to a 3-4 km thickness for the viscous body and a 10 16 Pa*s viscosity. The used of vertical component to determine the best solution (Fig. 9 e) is not determining.

Discussion.

Comparison between the different models.

Comparison of the WRMS for the models performed with the three hypotheses tested (Fig. 9) indicates that the models with only viscous relaxation on the lower flat of the MHT have much higher WRMS and can be excluded (Fig. 9a).

The other two models, afterslip alone or afterslip associated with viscous relaxation on the lower flat of the MHT, have respectively WRMS for the southern time series of 3.6 mm for the afterslip model and of 3.2 mm for the mixed model. The use of vertical time series (Fig. 9d and9f) do not help to choose between these two models.

The mixed model seems to better simulate the observations than the afterslip only model, but if the Relax software allowed the use of variable friction coefficients along the MHT, it would probably be possible to find a more satisfactory solution with the afterslip model, with probably a lower friction coefficient along the northern flat of the MHT than along the ramp and the southern or upper flat of the MHT.

Therefore the numerical modeling approach does not enable us to choose between the two hypotheses -afterslip alone or afterslip combined with viscous relaxation on the lower flat (Fig. 9c and 9e and Fig. 10). Nevertheless, it must be underlined that afterslip combined with viscous relaxation on the lower flat allows to simulate the early postseismic displacements as shown on A friction coefficient of 0.7 is somewhat higher than the values usually inferred. For example, the friction angle was estimated at 20° for the Indian crust in modeling by Berger et al. (2003). Such high values have nonetheless been found experimentally for numerous rocks (Kulhawi, 1975).

The afterslip evolution over time and space of our preferred model indicates that (Fig. and then by local high deviatoric stress to be released. Locally, slip also occurred southward, in areas near the southern end of the mesh, the upper flat segment of the MHT affected by afterslip.

Our modeling suggests changes of width of the upper flat of the MHT affected by afterslip (Fig. 11). The distribution of the afterslip reflects the distribution of the coseismic displacement model, the roughness of the afterslip distribution reflects the roughness of the co-seimic displacement model considered at the input of the model. It can therefore be assumed that with a smoother coseismic displacement model, we would also obtain a smoother afterslip model.

But we assume that the model roughness considered as input does not affect the conclusions of our models. Viscous relaxation controlled by a low viscosity zone along the lower flat of the MHT.

The time series for the stations located above the rupture zone, or to the south of it, are not correctly simulated by only viscous relaxation controlled by a ductile body centered on the northern flat. In the numerical models where afterslip along the upper flat of the MHT and along the crustal ramp has been introduced, the low-viscosity body (10 16 Pa.s) has to be thin (close to 3-4 km) to minimize the discard between model and data (Fig. 9e). We therefore believe that viscous relaxation controlled by a low-viscosity body cannot alone explain the recorded postseismic deformation and another source of deformation located to the south of the lower flat (i.e. afterslip) is needed. Models combining afterslip and viscous relaxation controlled by a thin, viscous body centered on the lower flat of the MHT better predict early post-seismic displacement (CHLM and KKN4 time series) than models considering only afterslip along the MHT (Fig. 10).

The low-viscosity body is a reasonable assumption as it may explain the low coupling along the lower flat documented for the interseismic period (Fig. 2a and12). Furthermore, the low-viscosity level has been suggested from the interpretation of deep images of the Himalayan crust (e.g. [START_REF] Zhao | Deep seismic reflection evidence for continental underthrusting beneath southern Tibet[END_REF][START_REF] Nelson | Partially molten middle crust beneath southern Tibet: synthesis of project INDEPTH results[END_REF]. Since these articles, new seismic experiments demonstrated that the low velocity zones in Tibet is discontinuous, restricted to surfaces with maximum vertical extension of 10 km and maximum horizontal length about 50 km (Hetenyi et al., 2011).

A channel flow model [START_REF] Beaumont | Himalayan tectonics explained by extrusion of a low-viscosity crustal channel coupled to focused surface denudation[END_REF] was even proposed to describe the evolution of Himalayan tectonics and is based on a very thick ductile zone (more than 10 km thick); so thick that the horizontal velocity depicts a peculiar gradient, with a maximum velocity in the center of the ductile zone. Our results nonetheless suggest that the ductile zone is 3 -4 km thick beneath the extreme southern part of the Tibetan plateau. Therefore, we suggest that the recent tectonism of the Himalaya is not controlled by channel flow deformation at depth, although our results does not preclude a thicker ductile zone to the north.

If we assume that postseismic deformation ends after 6 years, our results (Fig. 12) show ( 1 Balakot-Bagh earthquake [START_REF] Jouanne | Postseismic deformation in Pakistan after the 8 October 2005 earthquake: evidence of afterslip along a flat north of the Balakot-Bagh Thrust[END_REF][START_REF] Wang | Space geodetic observations and models of postseismic deformation due to the 2005 M7.6 Kashmir (Pakistan) earthquake[END_REF], for the deformation consecutive to the 1995 Chi-Chi earthquake [START_REF] Yu | GPS measurement of postseismic deformation following the 1999 Chi-Chi, Taiwan, earthquake[END_REF][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, Taiwan[END_REF], and for the 1995 Jalisco earthquake [START_REF] Hutton | Slip dynamics during and after the 9 October 1995 Mw = 8.0 Colima-Jalisco, Mexico[END_REF].

Importance of afterslip in the seismic cycle.

As illustrated in Figure 12, the effects of afterslip consecutive to the Gorkha earthquake remain confined to the northern part of the rupture area, the crustal ramp, and the northern flat of the MHT. The afterslip does not, then, homogenize coseismic slips throughout the rupture area but The rupture of the Gorkha earthquake is indicated by a blue line and the rupture of its main aftershocks, on May 12, 2016, is drawn with an orange line. The star indicates the nucleation of the Gorkha earthquake to the west of its rupture. The dashed green lines indicate the location of the Judi lineament [START_REF] Kayal | Microearthquake Seismology and Seismotectonics of South Asia[END_REF] that separates the main rupture zone from the nucleation area (Kumar et al., 2017) and, on the eastern side, the Gaurishankar lineament [START_REF] Mugnier | Segmentation of the Himalayan megathrust around the Gorkha earthquake (25 April 2015) in Nepal[END_REF]. Errors ellipses are drawn for a 95% confidence level. GNSS stations installed by our team for this study are indicated in blue and permanent GNSS stations of LDG are indicated in orange. See Fig. 9 for location of the cross section.

  4 and 5) located to the north of the rupture, with 40 mm of southward displacement between April 25 and June 24, 2015 and 100 mm of southward displacement between June 24, 2015 and late 2016. The points located in the northern part of the rupture zone (KKN4, GURJ, BALE) are affected by moderate southward postseismic displacements, whereas the points located at the southern end of the rupture, near to and south of Kathmandu (CHIR, DAMA,NAST and XYAK), are subject to small postseismic displacement. It should be noted that the points at the western part of the rupture (LMJG, GORK and ANBU) are not affected by postseismic displacements whereas the points located several kilometers to the east of the April 25, 2015 rupture (CHAR, JIR2, MALU) are affected by significant postseismic displacements reaching 35 mm (Fig.4). Given that these points were installed after the Mw 7.2 main aftershock of May 12, 2015, these postseismic displacements could be related to the postseismic deformation induced by this earthquake, but the lack of change in the cumulative postseismic displacements in the CHLM time series indicates, before and after the May 12, 2015 earthquake indicates that this second large earthquake does not induce a significant postseismic deformation.Significant postseismic displacements are therefore restricted to the northern part of the rupture zone and to an area to the north of it(Figs 4, 5, 6 and 7). It is important to note that there are no clear postseismic displacements to the south of the Gorkha earthquake rupture and, like[START_REF] Mencin | Himalayan strain reservoir inferred from limited afterslip following the Gorkha earthquake[END_REF], we exclude the occurrence of significant afterslip along the southern part of the MHT. Therefore, the slip deficit along the MHT south of Kathmandu is not resorbed aseismically by slip along the MHT from the rupture zone to the MFT where the MHT emerges.

Figures 5 ,

 5 Figures 5, 6 and 7 illustrate the postseismic deformation with the time series for the southern,

  : a postseismic deformation controlled by afterslip, viscous relaxation of the stress perturbation with a linear rheology, and a combination of these two mechanisms. This software computes the displacement and stress due to dislocations and the nonlinear postseismic time-dependent deformation, the latter either controlled by viscous bodies or by rate-strengthening friction faults. This software takes gravity into account and is based on a Fourier-domain elastic Green's function and an equivalent body-force representation of deformation mechanisms. We consider a 512*512*512 grid size model around the relevant faults with 0.7*0.7*0.7 km sampling which corresponds at a 358*358*358 km grid.

  , d) formed by (from 50 to 500 mm/year), and the friction along the relevant planes (from 0.2 to 0.7) for (a -b) σ values of 10 and 22 MPa (22 MPa being the stress drop estimated by Denolle et al., 2015). The tests of the (a -b) σ values of 22 and 10 MPa (Fig. 9 c and d) allow us to determine that the best fit fort the southern component time series is obtained for the (a -b) σ value of 22 MPa, which is the value of the stress drop estimated by Denolle et al. (2015) , a friction ratio of 0.7 and a = 200

Figure 10

 10 Figure10for the CHLM and KKN4 time series which is not the case if we consider the afterslip

  11): the lower flat of the MHT is affected by slow large-scale afterslip decreasing over time, whereas the crustal ramp is affected by negligible afterslip without significant lateral variations.The northern part of the upper flat is characterized by local patches of high afterslip; the latter reached 1 m after two years with significant lateral variability and a significant decrease of afterslip evolution between D + 30 days and D + 6 months. The areas of the upper flat affected by significant slip are located close to the flat-ramp transition corresponding to the northern boundary of the rupture. Areas of the upper flat locally affected by significant slip are characterized by a local high coseismic displacement gradient inGrandin et al.'s model (2016)

  Our modelling clearly differs from the model proposed by[START_REF] Wang | Observations and modeling of coseismic and postseismic deformation due to the 2015 Mw 7.8 Gorkha (Nepal) earthquake[END_REF], we allow afterslip along the northern part of the 2015 rupture whereas these authors supposed that no afterslip occurred along the upper flat of the MHT. Our model allows good simulation of time series of points located in the southern part of the area affected by postseismic deformation BALE, BHAZ, GURJ, KKN4, JIR2 and CHAR whereas the model of Wang and Fialko does not allow a good simulation of the time series of JIR2, the points BALE, BHAZ, GURJ and CHAR, installed by our team, being not included in their study. On the contrary, the model of Wang and Fialko presents a better simulation of the CHLM and DCN4 time series. Our model failed to predict accurately the early postseismic deformation of CHLM and failed to well fit the time series of DCN4 and GUMB. We can then propose that rheological properties are different between the upper flat of the MHT and the crustal ramp and the northern flat. As in this study we present a larger number of time series sensitive to the afterslip along the upper flat, we can suppose that our solution reflects mainly the rheological properties of the upper flat whereas the results presented by Wang and Fialko, 2018 reflect mainly the rheological properties along the lower part of the MHT (ramp and northern flat).

  ) afterslip or viscous relaxation of a thin body (~3-4 km thick) of low viscosity along the lower flat of the MHT (up to 0.5 m after 6 years) and (2) moderate afterslip along the crustal ramp that extends down-dip from the Gorkha rupture (up to 0.5 m after 6 years) as well as (3) a highvelocity small-scale afterslip (up to 5 m after 6 years for the western segment) along the southern boundary of the upper flat affected by afterslip, and along the upper flat, high-velocity afterslip in areas characterized by a significant coseismic slip gradient during the earthquake (reaching 5.6 m after 6 years in its western part and 3.4 m in its eastern part). This pattern, in which moderate afterslip occurs along the crustal ramp and extends down-dip from the rupture, has also been shown for the postseismic deformation following the 2005

  only and partially in its northern part. Moreover, if we consider that the 1833 and 2015 earthquakes have the same rupture area, and if we take the cumulated slip along the MHT during the interseismic period as estimated by[START_REF] Jouanne | Estimation of coupling along the Main Himalayan Thrust in the central Himalaya[END_REF], then it appears that postseismic displacement along the MHT does not homogenize slip along the MHT between the ductile slip along the northern flat (4 m for 200 years) and the coseismic slip along the northern part of the upper flat. Conclusion. The cumulative postseismic velocity field clearly shows a lack of displacements to the south of Kathmandu. Afterslip along the southern part of the rupture area of the Gorkha earthquake and along the southern segment of the MHT, that was unaffected by the main shock, is then clearly excluded. The slip deficit along the southern part of the upper flat of the MHT, between Kathmandu and the Ganga Plain, would probably be transferred along the MHT by an earthquake of large magnitude, as may have been the case after the previous 1833 large earthquake along this segment of the MHT, or by a very large magnitude earthquake affecting the MHT from the ductile zone to the MFT. To the west and east of the Gorkha rupture, postseismic displacements are almost negligible, and afterslip does not affect the upper flat areas of the MHT to the west of the rupture zone. The simulation of the time series using Relax software suggests that the postseismic deformation consecutive to the April 25, 2015 earthquake is controlled by: heterogeneous afterslip along the northern part of the upper flat of the MHT corresponding to the northern part of the rupture; afterslip along the crustal ramp reaching 0.17m after two years; and afterslip reaching 0.20 m after two years or viscous relaxation of a thin (~3-4 km thick) ductile zone along the lower flat of the MHT. Viscous relaxation along a ductile body centered on the lower flat of the MHT alone cannot explain alone the postseismic deformation. Captions Graphic abstract Example of observed and simulated postseismic time series and test of mechanisms controlling the postseismic deformation, afterslip along the northern part of the upper flat, the ramp and the flower flat of the Main Himalayan Thrust (blue) or a combination of afterslip along the northern part of the upper flat and the ramp of the MHT and a viscous relaxation controlled by a low viscosity body centered on the lower flat of the MHT (red). The simulated time series are respectively drawn in blue and red.

Figure 1 .

 1 Figure1. Location of the Gorkha earthquake along the upper flat of the MHT. Dark purple for ramp geometry; pink, light purple, and green for fully locked, partial coupled, and ductile MHT, respectively; red for rupture zones. The rupture area is bounded by lateral ramps(Mugnier et al., 2017).

Figure 2 .

 2 Figure 2. Description of inter and co-seismic deformation used for modeling the post-seismic deformation with Relax software. (a) Coupling along the Main Himalayan Thrust during the

Figure 3 .

 3 Figure 3. Observed[START_REF] Jouanne | Estimation of coupling along the Main Himalayan Thrust in the central Himalaya[END_REF] and simulated interseismic velocities using the model proposed by[START_REF] Jouanne | Estimation of coupling along the Main Himalayan Thrust in the central Himalaya[END_REF]. Velocities are expressed in the India fixed reference frame.

Figure 4 .

 4 Figure 4. Postseismic displacements following the Gorkha earthquake for the June 24, 2015-December 2016 period. The recorded displacements are corrected for interseismic displacements.

Figure 5 .

 5 Figure 5. Southern components for the cGNSS station displacements, corrected from interseismic velocities and seasonal signals estimated during the period before the main shock or extrapolated for the new stations. The locations of cGNSS stations are indicated in Figure 2.

Figure 6 .

 6 Figure 6. Eastern components for the cGNSS station displacements, corrected from interseismic velocities and seasonal signals estimated during the period before the main shock or extrapolated for the new stations. The locations of cGNSS stations are indicated in Figure 2.

Figure 7 .

 7 Figure 7. Vertical components for the cGNSS station displacements, corrected from interseismic velocities and seasonal signals estimated during the period before the main shock or extrapolated for the new stations. We present only the time series that seem to be correctly correct for seasonal components. The locations of cGNSS stations are indicated in Figure 2.

Figure 8 .

 8 Figure 8. Geometry and rheological properties considered in the numerical simulations of postseismic deformation. (a) Viscous relaxation controlled by a high-fluidity body embedded in a half-space characterized by low fluidity; (b) afterslip along the northern part of the upper flat of the MHT, along the ramp, and along the lower flat of the MHT; and (c) viscous relaxation controlled by a high-fluidity body embedded in a half-space characterized by low fluidity and afterslip along the crustal ramp and along the northern part of the upper flat.

Figure 9 .

 9 Figure 9. Search for the best solution, with the plot of the WRMS, for (a) the viscous relaxation hypothesis controlled by a viscous body along the lower flat of the MHT; (b) the afterslip hypothesis with (a-b)σ= 10 MPA; (c) the afterslip hypothesis with (a-b)σ= 22MPA with the southern time series, (d) the afterslip hypothesis with (a-b)σ= 22MPA with the vertical time series; (e) afterslip and viscous relaxation controlled by a viscous body along the lower flat of the MHT with the southern time series and (f) afterslip and viscous relaxation controlled by a viscous body along the lower flat of the MHT with the vertical time series .

Figure 10 .

 10 Figure 10. Observed and simulated postseismic time series with our preferred afterslip model, the afterslip & viscous relaxation model, and the viscous relaxation model. The simulated time series with our preferred afterslip and afterslip & viscous relaxation models almost always become conflated. Blue curve: afterslip hypothesis, red curve: afterslip and viscous relaxation controlled by a viscous body along the lower flat of the MHT.

Figure 11 .

 11 Figure 11. (a) Evolution of slip over time and space, predicted by our preferred afterslip model, expressed for selected dates. The locations of the permanent GNSS stations used to constrain our modeling are indicated by circles, (b) zoom of afterslip distribution after two years.

Fig. 12 .

 12 Fig. 12. Comparison between interseismic slip rate, seismic slip, and after-slip predicted after 6 years, supposed to be the end of postseismic deformation, along the Main Himalayan Thrust.

Acknowledgments.

The figures and maps were prepared using Generic Mapping Tools software (Wessel and Smith, 1995) and Paraview software (Ahrens et al., 2005;Ayachit, 2015). This material is partially based on data provided by the Nepal permanent GNSS network (DMG, Caltech, DASE) and data services provided by the UNAVCO Facility with support from the National Science Foundation (NSF) and National Aeronautics and Space Administration (NASA) under NSF Cooperative Agreement No. EAR-0735156. Funding was from the Labex OSUG@2020, INSU-CNRS, the AAP of Savoie Mont Blanc University, the IRD, and the ANR Bhoutanepal.

We thank two anonymous reviewers for their detailed and constructive reviews.


Nadia El Alami

El Hassani

Abdoullatif Baraket

Selim Boudjaoui

Ernandes Taveira

Tenório Neto

Joan Bausells

Nezha El

Benachir Bouchikhi

Abdelhamid Elaissari

Abdelhamid Errachid

Nadia Zine
email: nadia.zine@univ-lyon1.fr

Nadia El

Alami El Hassani

Nezha El Bari

Development and application of a novel electrochemical immunosensor for tetracycline screening in honey using a fully integrated electrochemical Bio-MEMS

Keywords: Bio-MEMS, Magnetic nanoparticles (MNPs), Immunosensor, Soft lithography, micro contact printing, Food control

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

Introduction

Tetracycline (TC) is a broad-spectrum antibiotic used as a veterinary drug to prevent and to treat diseases caused by Gram-positive and Gram-negative bacteria, atypical organisms such as Chlamydiae, Mycoplasmas, Rickettsia, and protozoan parasites [START_REF] Eliopoulos | Tetracycline therapy: update[END_REF]. The favorable antimicrobial properties of this drug and the absence of significant adverse side effects have led to their extensive use in the therapy of both human and animal infections [START_REF] Chopra | Tetracycline antibiotics: mode of action, applications, molecular biology, and epidemiology of bacterial resistance[END_REF]. Furthermore, it has been reported as the most intensely growth-promoting agent widely used in animal husbandry [START_REF] Economou | Agriculture and food animals as a source of antimicrobialresistant bacteria[END_REF]. Consequently, despite their proven efficacy, the intensive use of this drug is no longer encouraged since its remain residues can be accumulated in the animal bones, as well as in fish, meat, eggs, and milk, on account of the complexation of TC with Ca 2+ [START_REF] Nelson | Brief communication: mass spectroscopic characterization of tetracycline in the skeletal remains of an ancient population from Sudanese Nubia 350-550 CE[END_REF]. To ensure human food safety, the European Union (EU) has established the maximum residue limits (MRLs) for TC in edible products and milk at 100 µg.kg -1 , in eggs at 200 µg.kg -1 and in the liver at 300 µg.kg -1 (European Commission, 1999). However, the MRLs for TCs in honey has not yet been established by Codex, since the treatment of honeybee diseases with antibiotics is unallowable in the EU [START_REF] Kivrak | Development of a rapid method for the determination of antibiotic residues in honey using UPLC-ESI-MS/MS[END_REF]. Despite banning, some countries, such as Switzerland, UK, and Belgium still use TC treatment for hive diseases and have established, therefore, the MRLs for TC ranging from 10 to 50 µg.kg -1 [START_REF] Bargańska | Determination of antibiotic residues in honey[END_REF].

However, TC residues can be successfully determined in various food products requiring highly performing sample preparation steps. Several analytical methods currently exist for TC detection and quantification. High performance liquid chromatography coupled to mass spectrometry [START_REF] Pokrant | In-House Validation of HPLC-MS/MS Methods for Detection and Quantification of Tetracyclines in Edible Tissues and Feathers of Broiler Chickens[END_REF], capillary electrophoresis coupled with electrochemiluminescence [START_REF] Deng | Pharmacokinetics and residues of tetracycline in crucian carp muscle using capillary electrophoresis on-line coupled with electrochemiluminescence detection[END_REF], enzyme immunoassay [START_REF] Gaurav | ELISA based monitoring and analysis of tetracycline residues in cattle milk in various districts of Punjab[END_REF], liquid chromatography-mass spectroscopy (LC-MS) [START_REF] Desmarchelier | Determination of five tetracyclines and their epimers by LC-MS/MS based on a liquid-liquid extraction with low temperature partitioning[END_REF] have been described for confirmatory analysis. Other methodologies based on spectroscopy analysis are described [START_REF] Qin | A high-sensitivity terahertz spectroscopy technology for tetracycline hydrochloride detection using metamaterials[END_REF]; unfortunately, they suffer from a lack of sensitivity compared to chromatographic techniques [START_REF] Mitra | Sample Preparation Techniques in Analytical Chemistry[END_REF]. Subsequently, Most of these approaches often require the use of extraction techniques, such as liquid-liquid [START_REF] Desmarchelier | Determination of five tetracyclines and their epimers by LC-MS/MS based on a liquid-liquid extraction with low temperature partitioning[END_REF] and solid phase extraction [START_REF] Shi | Group-selective molecularly imprinted polymer solid-phase extraction for the simultaneous determination of six sulfonamides in aquaculture products[END_REF], in order to accomplish both the analytes preconcentration and sample cleanup [START_REF] O'connor | Analysis of tetracycline antibiotics in soil: advances in extraction, clean-up, and quantification[END_REF]. However, the principal limitations of these steps lie in the high costs, the long time of achievement and the requirement of advanced technical skills [START_REF] Joshi | Recent advances of ionic liquids in separation science and mass spectrometry[END_REF]. For this reason, simple, rapid and accurate methods are recommended for the on-site screening of low TC residues without needing extraction or clean-up steps. Due to their advantages of high selectivity, rapid detection, and in-situ applications, several optical and electrochemical techniques based on biosensors have been investigated [START_REF] Lan | Recent advances in nanomaterial-based biosensors for antibiotics detection[END_REF]. Therefore, aptamer-based sensing techniques were widely used for the food safety determination. Notably, there is a growing rise in the aptasensors fabrication for the TC detection [START_REF] Han | Label-free photoelectrochemical aptasensor for tetracycline detection based on cerium doped CdS sensitized BiYWO 6[END_REF][START_REF] Tang | The development of a graphene oxide-based aptasensor used for the detection of tetracycline in honey[END_REF] with a few applications in honey samples [START_REF] Wang | A direct competitive assay-based aptasensor for sensitive determination of tetracycline residue in honey[END_REF][START_REF] Wang | Development of an indirect competitive assay-based aptasensor for highly sensitive detection of tetracycline residue in honey[END_REF]. Nevertheless, the main limitations of these systems are related to their relative low detection signals [START_REF] Yi | A universal electrochemical sensing system for small biomolecules using target-mediated sticky ends-based ligation-rolling circle amplification[END_REF]. Furthermore, in some recent works, the molecularly imprinted sensors had been successfully applied to the analysis of antibiotic residues in honey samples [START_REF] Bougrini | Development of a novel sensitive molecularly imprinted polymer sensor based on electropolymerization of a microporous-metal-organic framework for tetracycline detection in honey[END_REF][START_REF] Lv | Adsorption-controlled preparation of molecularly imprinted hybrid composites for selective extraction of tetracycline residues from honey and milk[END_REF]. This type of sensors, however, often show relatively low sensitivity and specificity when the imprinted membrane is fragile. Hence, on the last few decades, the development of immunosensors grew tremendously to the detection of drug residues in food. The use of these analytical systems has attracted considerable attention because of the low detection limits and the high selectivity for analyzing complex samples [START_REF] Majdinasab | An Overview on Recent Progress in Electrochemical Biosensors for Antimicrobial Drug Residues in Animal-Derived Food[END_REF]. In this regard, some immunosensors were performed for TC detection in milk samples [START_REF] Conzuelo | Integrated disposable electrochemical immunosensors for the simultaneous determination of sulfonamide and tetracycline antibiotics residues in milk[END_REF][START_REF] Liu | Electrochemical immunosensor based on the chitosan-magnetic nanoparticles for detection of tetracycline[END_REF] with the limits of detection being 0.032, 3.9, and 0.85 ng.mL -1 , respectively.

Herein, a novel and highly sensitive approach is proposed for the detection of remaining TC residues in honey samples. This immunosensor was conducted on bio-microelectro-mechanical systems (Bio-MEMS) transducers based on gold micro working electrodes (µWEs) with fully integrated, reference and counter electrodes. The novelty of this work is related to the functionalization process of µWEs, which was achieved through two-and threedimensional (2D and 3D) shapes. These functionalizations were conducted, in one hand, by electroaddressing diazonium salt and, in the other hand, by electrodepositing of a new structure of magnetic nanoparticles coated poly(pyrrole-co-carboxylic acid) (Py/Py-COOH/MNPs), while combining the two methods with a preconcentration technique. The Py/Py-COOH/MNPs as solid supports of biomolecules, is the pioneering aspect of this work dedicated to the combination of submicron, magnetic conducting particles. The quantification of this analyte was performed through competitive detection procedure with TC molecules immobilized on the µWEs surface toward polyclonal TC antibody. The prepared immunosensor was successfully applied to enhance the detection limit of TC in complex matrices such as honey.

Experimental

Reagents and solutions

Tetracycline (TC), Doxycycline (DX), Oxytetracycline (OX), Chlortetracycline (CT), 1-Ethyl-3-(3-dimethylaminopropyl)carbodiimide (EDC), Ethanolamine, Sulfuric acid (H 2 SO 4 ) (30 wt.% in H 2 O), 4-aminophenylacetic acid (CMA), sodium nitrite (NaNO 2 ), Sodium cyanoborohydride (95%), hydrochloric acid (37%, HCl), Carbon tetrachloride, heptane, octadecyltrichlorosilane (OTS), Pure ethanol (99.8%), Acetone (99.9%), Potassium hexacyanoferrate II, Potassium hexacyanoferrate III, Phosphate buffer saline (PBS) with pH 7.4 were all purchased from Sigma-Aldrich, France. N-Hydroxysuccinimide (NHS) and Hydrogen peroxide (H 2 O 2 ) (35%, v/v) were purchased from Acros Organics, France. 11-(triethoxysilyl)undecanal (TESUD) was purchased from abcr GmbH, Germany.

Polydimethylsiloxane (PDMS: Sylgard 184) was purchased from Dow Corning, France. The sheep polyclonal antibody raised against TC (Ab-TC) was supplied by Abcam, France. A fluorescent labeled antibody was generated by conjugation of Ab-TC with rhodamine using the rhodamine fast conjugation kit (ab188286) by Abcam, France. Rosemary honey sample, guaranteed tetracycline-free was supplied by Secrets d'Apiculteur, Lyon, France. Honey sample with an unknown concentration of TC was commercially provided. Milli-Q ultrapure water was used in all experiments.

Procedure for immunosensor elaboration

The development of this immunosensor was conducted by three methods (Fig. 1).

Before any functionalization, the Bio-MEMS device surface was cleaned with ethanol in an ultrasonic bath for 10 min, dried under nitrogen flow, then placed under UV-Ozone Cleaner-ProCleaner™ Plus from BioForce for 30 min to remove all organic contaminants.

Diazonium immobilization onto gold microelectrodes

The cleaned microelectrodes (µWEs) were functionalized by CMA solution (3 mM) previously diazotized in an aqueous solution of HCl (20 mM) and NaNO 2 (20 mM) for 15 min at 4 °C. Subsequently, nine-repetitive cyclic voltammograms were applied from 0.3 V to -1.0 V with a scan rate of 80 mV.s -1 . This was done sequentially until all desired µWEs were modified with CMA. The carboxylic groups formed on the electrode surfaces were activated by a mixture of 0.4 M EDC and 0.1 M NHS, prepared in ethanol for one hour at room temperature. Afterwards, the activated electrode surfaces were incubated in 40 µL of TC at 100 µg.mL -1 for 30 min at 4 °C. Then, the immunosensor was treated with ethanolamine (1% v/v) in PBS buffer (pH 7.4) for 10 min at room temperature. This step is crucial to prevent nonspecific bonding phenomenon at the detection stage of TC (Fig. 1A).

Preconcentrating technique

The core-shell magnetic nanoparticles coated with poly(pyrrole-co-pyrrole carboxylic acid) (Py/Py-COOH/MNPs) were exclusively synthesized by [START_REF] Tenório-Neto | Submicron magnetic core conducting polypyrrole polymer shell: Preparation and characterization[END_REF] using a seeded-polymerization technique. The morphology of these nanoparticles was already presented in our previous work [START_REF] Hassani | Novel strategy for sulfapyridine detection using a fully integrated electrochemical Bio-MEMS: Application to honey analysis[END_REF] using transmission electron microscopy (TEM). Here, the magnetic nanoparticles stock solution (100 µL) was washed three times with PBS buffer (pH 7.4) in a magnetic field to separate the nanoparticles from the storage solution. EDC and NHS solutions (each at a concentration of 100 mM) were prepared immediately before use. A volume of 250 µL each of EDC and NHS solutions were added to the washed magnetic nanoparticles and incubated with slow rotation at room temperature for 90 min. Then, the reaction mixture was rinsed three times with HCl (1 mM) at 4 °C in a magnetic field. Subsequently, 100 µL of purified Ab-TC polyclonal antibody (100 µg.mL -1 ) was added to the activated nanoparticles and incubated for 2 h at 4 °C. The antibody-coated nanoparticles were washed three times with PBS buffer (pH 7.4). The non-reacted activated carboxylic acid groups were blocked with ethanolamine (1% v/v) in PBS buffer for 10 min.

The antibody-coated magnetic nanoparticles were then separated from the mixture, resuspended in 500 µL of PBS buffer and incubated with different concentrations of TC (Fig. 1B). Simultaneously, the gold µWEs were functionalized by CMA and TC as described in paragraph 2.2.1.

Functionalization with Py/Py-COOH/MNPs

This functionalization was performed by electrodepositing Py/Py-COOH/MNPs onto microelectrodes by using pulsed chronoamperometric technique (PCA). Briefly, two continuous repetition potentials were applied: the first potential E 1 =-0.2 V for t 1 =0.1 s, induces a cathodic electrodeposition process, while the second potential E 2 =0.2 V for t 2 =0.1 s oxidizes the electrode surface, and in particular, cleans it from eventual extraneous deposited species. The cycling process was iterated 30 times. Then, the gold µWEs were rinsed with distilled water. Besides, it is interesting to point out that the amplitude of the charge-discharge process remains practically unchanged during the prolonged cycling process, indicating that the electrode surface exhibits a homogenous Py/Py-COOH/MNPs deposition. The activation of carboxylic groups, present in the Py/Py-COOH/MNPs, was carried out by the mixture of EDC/NHS (0.4 M/0.1 M) prepared as described above. The biofunctionalization and blocking steps were performed by incubating the sensor in TC solution (100 µg.mL -1 ) and ethanolamine/PBS (1% v/v), respectively. Finally, the immunosensor was rinsed with PBS and used for TC detection (Fig. 1C). The homogeneity and the roughness of the resulted µWEs surface were already examined in our preceding study [START_REF] Hassani | Novel strategy for sulfapyridine detection using a fully integrated electrochemical Bio-MEMS: Application to honey analysis[END_REF]. This was conducted by two-, and three-dimensional surface topography using atomic force microscopy (AFM). The interest of using this technique is to ensure that the µWE surfaces were entirely homogeneous with the presence of spherical shapes indicating the presence of MNPs coated with Py/Py-COOH.

Micro-contact printing

The soft-lithographical technique called microcontact printing (µCP) of selfassembled monolayers (SAMs) was employed to check the immune TC/Ab-TC reactivity. For this purpose, an elastomeric stamp based on polydimethylsiloxane (PDMS) was fabricated by replica molding (RM) as already described in previous work [START_REF] Baraket | Diazonium modified gold microelectrodes onto polyimide substrates for impedimetric cytokine detection with an integrated Ag/AgCl reference electrode[END_REF] . A mixture of pre-polymer PDMS and curing agent (10:1 w/w) was poured onto a silanized silicon-mold which contains micro-pillars on relief of its surface (Fig. 2a). Subsequently, complete degassing of PDMS/silanized silicon-mold was carried out to ensure that all air bubbles which may create defects on the PDMS stamp surface have been removed (Fig. 2b). After the polymerization step at 52 °C for 1 hour, the PDMS stamp was peeled-off from the silicon mold bearing the micro holes on its surface (Fig. 2c). The stamp was then inked by immersion for 1 min in heptane solution containing octadecyltrichlorosilane (OTS) (0.1%, v/v) and carbon tetrachloride (2%, v/v), then dried with a stream of nitrogen (Fig. 2d). Simultaneously, glass substrates were cleaned with acetone then ethanol, rinsed with distilled water and then dried with nitrogen. The substrates were activated by submersion in piranha solution (3:1 v/v of H 2 SO 4 : H 2 O 2 ) for 30 min then thoroughly rinsed with distilled water and dried with nitrogen. Afterwards, the µCP was realized by putting the prepared PDMS stamp in immediate and conformal contact with a freshly activated glass substrate (Fig. 2e). The PDMS stamp was peeled-off from the substrate, leaving then the SAMs of OTS on the glass surface which was incubated afterwards at 100 °C for 45 min to enhance OTS adhesion.

The freshly prepared substrates were submerged in an ethanol solution containing 1% of TESUD, rinsed with ethanol, dried with a stream of nitrogen and then heated in a vacuum stove at 100 °C for 60 min to complete the silanization process (Fig. 2f). Subsequently, the functionalized substrate surface (FSS) was incubated in 250 µL of Ab-TC antibody (10 µg.mL -1 ) already diluted in 4 mM of sodium cyanoborohydride. The use of this reagent was conducted to avoid adversely reducing aldehydes to nonreactive hydroxyls, allowing covalent bonding of Ab-TC with the amine of TESUD (Fig. 2g). After that, the FSS were incubated in TC solution (0.1 pg.ml -1 ) for 1 h to allow the antibody-antigen interaction followed by rinsing in PBS and drying under nitrogen (Fig. 2h). Finally, the rhodamine-labeled polyclonal Ab-TC* (10 µg.ml -1 ) were dropped on the FSS for 1 h (Fig. 2i). The samples were then rinsed in PBS, dried under nitrogen and observed by fluorescence microscopy.

Fluorescence microscopy characterization

Fluorescence images, as a rapid tool showing the recognition of Ab-TC for its corresponding antigen (TC), were taken using a fluorescence microscope (Zeiss Axioplan 2 Imaging), equipped with 10× and 40× lenses and a monochrome camera. Samples were observed by fluorescent light: TC sample was excited with a 550 (±25) nm band-pass filter and fluorescence from the sample was observed with a 605 (±70) nm band-pass filter.

Electrochemical measurements

All steps of gold µWEs modification were characterized by cyclic voltammetry (CV) and electrochemical impedance spectroscopy (EIS) techniques using a multi-channel potentiostat (VMP3 Bio-logic-Science Instrumentation, France). The electrochemical configuration was performed by using integrated µWEs, RE and CE of the Bio-MEMS. All the measurements were performed in a solution of 5.0 mM [Fe(CN) 6 ] 3/4-prepared in PBS buffer (pH 7.4) and in a Faraday cage. For the cyclic voltammetric measurements, the potential was scanned from 0.3 V to 0.6 V with a scan rate of 100 mV.s -1 . While the electrochemical impedance was measured in the 150 mHz-200 kHz frequency range with a DC potential of -0.3 V/ref and an AC potential of 90 mV.

Results and discussion

Fluorescence analysis

Fluorescent imaging is considered as an effective tool for bio-functionalization analysis. It's used to certify that the corresponding immune detection process is taking place based upon the recognition of the antibody for its corresponding antigen. In this study, fluorescent antibodies were immobilized (chemically) onto a glass substrate for the detection of TC antigens. This immuno-fluorescent test has reported the strong reactivity of Ab-TC against TC which can, therefore, be applied for immunosensor development.

Electrochemical characterization of modified gold µWEs

The behaviors of CMA and Py/Py-COOH/MNPs films were checked for each step of the Bio-MEMS functionalization by CV analysis. The process of CMA electro-addressing is shown in Fig. S1A (Supplementary information). The initial cycle shows a large and irreversible cathodic wave with a peak potential at -0.8 V revealing diazotated CMA onto the gold µWEs surface by diazonium salt reduction. However, the reduction wave has been kept unchanged upon the deposition of CMA. This functionalization was checked by variation of the CV before and after gold µWEs modification (Fig. S1B) (Supplementary information). In this figure, a significant decrease in peak-to-peak of the CV cycles was noticed after CMA electro-deposition against the bare gold µWEs. This was referred to the broad passivation region of this microelectrodes.

While Fig. S1C (Supplementary information) shows the chrono-amperometric process for Py/Py-COOH/MNPs electrodeposition. During this step it was interesting to specify that the amplitude of the charge-discharge process remains nearly unchanged, indicating that the electrode surface displays homogenous Py/Py-COOH/MNPs layer. CV analysis were also investigated to examine the behavior of the µWEs surfaces after Py/Py-COOH/MNPs electrodeposition. Fig. S1D (Supplementary information) shows the considerable increase in the peak current of cyclic voltammograms after Py/Py-COOH/MNPs deposition compared to bare gold. This was explained by the conductive behavior of Py/Py-COOH coated on MNPs.

At these steps, the modified gold µWEs were ready for bio-functionalization with the target analyte.

Optimization of conditions for detection

The amount of TC immobilized on the surface of gold µWEs as well as the antibody concentration were subject to optimization described in (Supplementary information).

Competitive detection of TC

In this study, the quantification of the analyte has been made by the competitive approach to enhance the detection signal of TC with a small molecular weight (444 Da) compared to Ab-TC (160 kDa). Furthermore, this approach was advantaged in the present work to make pre-concentration allowing to detect only TC and eliminate other irrelevant molecules. This step was performed by incubating the µWE surfaces of the Bio-MEMS in a mixture of a fixed concentration of Ab-TC previously optimized (10 µg.mL -1 ) and decreasing levels of standard solutions of TC from 1000 pg.mL -1 to 0.1 pg.mL -1 during 30 min at 4 °C for each concentration. At least three independents replicated experiments were carried out.

Ultimately, the electrode surface was rinsed with PBST (10 mM PBS with 0.05% Tween 20) and analyzed with EIS. Fig. 4A-C show the obtained Nyquist plot semi-circles for the three functionalizations of µWE surfaces. It can be noticed that for the three approaches, the Nyquist plot semi-circles increase by decreasing TC concentrations. As awaited, this level decreases produce an increase in the charge transfer resistance owed to the presence of free antibodies reacting with TC immobilized on the µWEs surface.

An excellent fitting between the simulated and experimental spectra was obtained for each TC concentration by using the Randles equivalent circuit shown in the inset of Fig. 4D (where R s is the solution resistance, R ct is the charge-transfer resistance, W is the Warburg impedance and Q is the constant phase element). For each TC concentration, the value of ∆R/R 0 was calculated when R is the charge transfer resistance. The normalized data show three linear equations: y a = -0.208 x + 0.960, y b = -0.148 x + 0.578 and y c = -0.109 x + 0.474 with the determination coefficients of 0.962, 0.986 and 0.946, respectively. As can be seen from this figure, the detection sensitivity was highly improved for the Py/Py-COOH/MNPs gold surface modification. For this approach, the configured impedimetric immunosensor have provided the more sensitivity of 0.208 mL.pg -1 with a limit of detection (LOD) of 1.2 pg.mL -1 . It is appealing to notice that this LOD value is nearly 83000 times lower than the maximum residue limits of TC in honey fixed by European regulation. As listed in Table 1, the proposed system provides relatively wide linear ranges and low LODs in comparison with other published studies.

Interferences study

The selectivity of the proposed system was investigated in the presence of various antibiotics from the same family of cyclines that can have cross-reactivity with Ab-TC antibody. This study was performed by comparing the responses for TC detection with three analogous molecules namely: chlortetracycline (CT), doxycycline (DX), and oxytetracycline (OX) (Fig. S3) (Supplementary information). Consequently, this immunosensor was found to be very selective for TC and can provide reliable results regardless of the presence of the interfering molecules.

Analysis of spiked honey samples

The TC determination in real samples was performed by spiking analyte in blank honey matrices which indicate that the calibration curve is suitable to determine the TC content in real unknown samples (Supplementary information). The results of the added spiked samples are shown in Table S1 (Supplementary information). Good recoveries from 80% to 98% were obtained with the RSD ranging between 1% and 6%, which indicate that the calibration curve is suitable to determine the TC content in real unknown samples.

Detection of TC in unknown honey sample

The standard addition method was further employed to determine the TC concentration in unknown honey sample. The competitive detection of TC was performed by a mixture of different standard solutions (3 pg.mL -1 , 5 pg.mL -1 , 10 pg.mL -1 , 15 pg.mL -1 )

with the presence of Ab-TC (10 µg.mL -1 ). For all TC concentrations, the normalized charge transfer resistance ∆R/R 0 were measured. The extrapolation of regression line predicts the TC concentration in the diluted honey sample at 2.52 pg.mL -1 . Consequently, the concentration of TC in the unknown sample was 25.2 pg.mL -1 according to the dilution factor of 10 times (Fig. 5). This level of contamination remains lower than the MRLs of 10 000 pg.mL -1 established by Switzerland, UK, and Belgium.

Conclusions

In the present work, we have discussed the development process of a novel sensitive and highly selective immunosensor based on 2D and 3D network of immobilization with the diazonium salt, and a new structure of MNPs coated with Py-COOH for specific detection of TC. Moreover, the electrochemical analysis performed by these materials allow high reproducibility and rapid response. The high sensitivity of this immunosensor was evaluated for the three functionalization methods and was highly improved for the Py/Py-COOH/MNPs gold surface modification. For this approach, the configured impedimetric immunosensor have provided a sensitivity of 0.208 mL.pg -1 with a limit of detection of 1.2 pg.mL -1 . The selectivity of the present immunosensor platform toward TC was higher when compared to analogous molecules, namely DX, CT and OX. We have presented a novel and less costly fabrication process to develop a highly sensing device able to detect low levels of TC residues in honey samples without any complicated pretreatment. The miniaturization of this immunosensor will be very beneficial to control the quality of food products on an industrial scale. However, there are still some limitations of this work related to the construction time which is a little long. Therefore, our further investigations in this direction are to shorten the development time of this immunosensor, to meet the requirements of a large scientific community which can extend its use for other applications. Future works will also be done, to assess the multi-functionalization of the same device allowing a simultaneous detection of several analytes.

Fig. 1 .

 1 Fig.1. Schematic illustrations of competitive immunoassay of TC on modified µWEs by functionalization with: (A) CMA electroaddressing, (B) CMA electroaddressing followed by preconcentration of immunomagnetic nanoparticles cross-linked with Ab-TC antibody, (C) Py/Py-COOH/MNPs.

Fig. 2 .

 2 Fig. 2. Microcontact printing process to immobilize Ab-TC onto TESUD patterns and sandwich immunoassay with TC and rhodamine-labeled Ab-TC.

  Fig. 3A-B shows the fluorescent pattern of Ab-TC *R . The red fluorescent regions indicate tags in which the antibodies are specifically immobilized, while the extinct areas are those blocked with OTS. In the Fig.3C, we notice an absence of fluorescence due to the non-cross-reactivity of the doxycycline (DX), as interfering molecules, towards Ab-TC.

Fig. 3 .

 3 Fig. 3. Fluorescent images of (Ab-TC *R ) chemically immobilized on glass substrates after recognition of TC, (A) Positive homogenous pattern with the blackened regions corresponding to the OTS blocked glass surface. (B) Magnification of the positive pattern showing specific detection between Ab-TC *R and TC. (C) Fluorescent image of the biorecognition assay between doxycycline and Ab-TC *R .

Fig. 4 .

 4 Fig. 4. Nyquist impedance plots in 5 mM of [Fe(CN) 6 ] 3-/4-for competitive detection of TC on the modified µWEs by: (A) Functionalization with CMA, (B) Functionalization with CMA followed by preconcentration step, (C) Functionalization with Py/Py-COOH/MNPs, where (a) Bare gold, (b) After TC immobilization at 100 µg.mL -1 , (c-g) After competitive detection of TC at 1000, 100, 10, 1 and 0,1 pg.mL -1 , respectively. (D) The sensitivity of the immunosensor by normalization of Nyquist plot data, inset: the equivalent circuit used for EIS fitting.

Fig. 5 .

 5 Fig. 5. Standard addition plots of different TC spiking in diluted unknown honey sample, using (a) Ab-TC and (b) Ab-155 antibodies.

Table 1 .

 1 Comparison of different techniques for TC detection.

	Methods	Linear range	LOD (ng.mL -1 )	Real	References
		(ng.mL -1 )		sample
	High-performance liquid-	15-5000	5	Milk	(Kargin et
	chromatography with				al., 2016)
	fluorescence detection
	Photoelectrochemical	0.2-1000.0	0.01	Drug	(Han et al.,
	aptasensor			sample	2018)
	Terahertz spectroscopy	0-20×106	0.45	Water	(Qin et al.,
				and	2017)

Acknowledgements

The authors acknowledge the financial support from PHC-Toubkal under the project No.32567YD, and funding from the European Union's Horizon 2020 research and innovation programme entitled HEARTEN under grant agreement No 643694" and microMole grant agreement No 643694 and No 653626, respectively.