Madhuchhanda Brahma
email: madhu.brahma@gmail.com

Marc Bescond
email: marc.bescond@im2np.fr

Demetrio Logoteta
email: demetrio.logoteta@im2np.fr

Ram Krishna Ghosh
email: rghosh1@nd.edu

Santanu Mahapatra
email: san-tanu@cedt.iisc.ernet.in

Germanane MOSFET for sub-Deca Nanometer High Performance Technology nodes

Keywords: Germanane, 2D materials, ballistic transport, ITRS, Non-Equillibrium Greens Function (NEGF) formalism

Ballistic transport in monolayer Germanane metal oxide field effect transistors (MOSFETs) are investigated for high performance(HP) applications. Characteristics of both n-and ptype transistors having channel lengths of 7 nm, 5 nm and 3 nm are studied and compared against the International Technology Roadmap for Semiconductor (ITRS) target of 2028. For this purpose we conceive single band effective mass Hamiltonian and double-band k.p Hamiltonian for n-and p-type channels respectively. The quantum transport model is based on the Non-Equillibrium Greens Function (NEGF) formalism self consistently coupled with Poisson equation. The ON current (I ON ) in 3 nm channel n-and p-MOSFETs, for a fixed OFF current I OFF =100 nA/µm, are found to be ∼ 890 µA/µm and 700 µA/µm respectively, which is indeed remarkable. We also observe that with the increase of channel lengths, the p-MOSFET starts to outperform the n-MOSFET in terms of I ON requirements as the direct source-to-drain tunneling gets suppressed. Other performance metrics like total gate capacitance, intrinsic switching delay and switching energy have also been calculated and found to be comparable to the ITRS 2028 HP technology requirements.

I. INTRODUCTION

I N recent times, atomically thin two dimensional (2D) materials are being extensively explored as an alternative channel material for transistors in an aspiration of CMOS technology growth beyond the Moores law. Most of such materials e.g. transition metal di-chalcogenides or TMD (MoS 2 , WS 2 , MoSe 2 , WSe 2 , MoTe 2 ), Phosphorene, Graphene, hBN etc. [START_REF] Radisavljevic | Integrated circuits and logic operations based on single-layer MoS 2[END_REF]- [START_REF] Novoselov | Two-dimensional atomic crystals[END_REF] are naturally available in layered form, where two-dimensional platelets are weakly bonded by vander Waals force to form three-dimensional crystals. In an alternate route, recent efforts are also observed to extract 2D analogues from conventional bulk semiconductors like Si, Ge, GaN etc. [START_REF] Kj Koski | The new skinny in two-dimensional nanomaterials[END_REF]- [START_REF] Bianco | Stability and exfoliation of Germanane: A germanium graphane analogue[END_REF]. Among these, hydrogenated monolayer of Germanium namely Germanane has attracted much attention. Unlike TMDs, which show higher band gap and higher effective mass, Germanane inherits a very interesting property of having low effective mass along with a relatively high band gap [START_REF] Bianco | Stability and exfoliation of Germanane: A germanium graphane analogue[END_REF] in the group of 2D materials. From recent first principle based assessment [START_REF] Krishna Ghosh | Germanane: A Low Effective Mass and High Bandgap 2-D Channel Material for Future FETs[END_REF], it was further concluded that chair morphology of Germanane has better potential for conventional and tunnel FET applications than its boat morphology. Stable structure of Germanane has been reported for the first time in 2013 [START_REF] Bianco | Stability and exfoliation of Germanane: A germanium graphane analogue[END_REF], and later used for fabrication of Schottky diodes and FETs [START_REF] Amamou | Large area epitaxial germanane for electronic devices[END_REF]- [START_REF] Gopal | Schottky diodes from 2D germanane[END_REF]. Since high performance switching devices demand channel materials with low effective mass to ensure high ON current, Germanane stands a very good chance in this respect. It is therefore important to theoretically estimate the performance limit of Germanane channel MOSFETs. However, previous studies in this context are either based on semi-classical approach [START_REF] Low Kain | Ballistic Transport Performance of Silicane and Germanane Transistors[END_REF] or limited only to the n-MOSFET [START_REF] Aziz | Performance Limit Projection of Germanane Field-Effect Transistors[END_REF]. In this work we focus on NEGF based ballistic transport in both n-and p-type Germanane chair FETs with channel lengths less than 10 nm to examine if such short Germanane FETs are able to meet the ITRS 2028 HP limit [21]. In section II of our paper, we describe the methodology used to simulate the transfer characteristics while in section III we present and analyse the results. From the transfer characteristics we find that the ON-state currents even in ultra-scaled Germanane nand p-MOSFETs of 3 nm show remarkably high values ,∼ 890 µA/µm and ∼ 700 µA/µm respectively. Both these values are comparable to the ITRS projected 2028 HP I ON value of 900 µA/µm. Finally, we have bench-marked the dynamic performance metrics of these sub-deca nanometer Germanane HP FETs with respect to the other state of the art 2D materials as well as the ITRS projection limits.

II. METHODOLOGY

In consideration with the ITRS-predicted physical gate length, effective oxide thickness, and power supply voltage for the 2028 HP devices, we conceive the device as shown in Fig. 1 with a channel length equal to the gate length and varying from 7 nm, 5 nm to 3 nm, along with an effective oxide thickness(EOT) of 0.47 nm, and a supply bias of V D = 0.6 V. The source(drain) lengths L S (L D ) are kept to be 30 nm each. Moreover, the gate work-function is also adjusted such that the OFF-state current (I OFF ) at zero gate bias can be fixed at 100 nA/µm. Monolayer Germanane chair has isotropic conduction band and valence bands (light hole and heavy hole) [START_REF] Krishna Ghosh | Germanane: A Low Effective Mass and High Bandgap 2-D Channel Material for Future FETs[END_REF]. Therefore, for n type MOSFET we have chosen a single band effective mass Hamiltonian with electron effective mass of 0.071m • . For the p -type MOSFET, we developed a twoband k.p Hamiltonian based on the approach described in [START_REF] Iv Tokatly | Orbital momentum Hall effect in p-doped graphane[END_REF] for its analogous Graphene counterpart and given as

H( k) ≡ γ1 2 ( kx 2 + ky 2 ) γ2 2 ( kx + i ky ) 2 γ1 2 ( kx -i ky ) 2 γ1 2 ( kx 2 + ky 2 ) , (1)
with the constants

γ 1 = 1 2 ( 1 m L + 1 m H ) γ 2 = 1 2 ( 1 m L -1 m H ) (2)
where m L and m H are the light hole and heavy hole effective masses. The band structure using k.p method is matched against the density functional theory (DFT) predicted values [START_REF] Krishna Ghosh | Germanane: A Low Effective Mass and High Bandgap 2-D Channel Material for Future FETs[END_REF] in the region close to the valence band maximum at the Γ point of the Brillouin zone by considering m L and m H equal to 0.071m • and 0.29m Since the main focus of this work is to meet the high performance on current limit, it is necessary to select an optimum source/drain doping for each length of the device. The source/drain doping for n-MOSFET and p-MOSFET with different channel lengths, has been chosen in the range varying from 3 × 10 12 /cm 2 to 5 × 10 12 /cm 2 , in order to ensure that maximum I ON is achieved for given I OFF as per the ITRS HP requirements. It has been found ( not shown here) that doping lower than 10 12 /cm 2 results in source exhaustion [START_REF] Kim | Source/drain doping effects and performance analysis of ballistic III-V n-MOSFETs[END_REF], whereas higher orders of doping result in degraded swing with the current rising rather slowly and not reaching the desired value at V D = 0.6 V. These doping concentrations result in the degeneracy of ∼ 0.03 eV and 0.14 eV in p-MOSFET and n-MOSFET respectively. The channel region is intrinsic. To model the ballistic transport, we solve Schrödinger's equation by using NEGF formalism self-consistently with the 2D-Poisson equation in the cross section of the device shown in Fig. 1. The retarded Green's function in the matrix expression is calculated as [START_REF] Datta | Quantum transport: atom to transistor[END_REF], [START_REF] Guo | Carbon nanotube electronics: Modeling, physics, and applications[END_REF]

] G r (k y , E) = [(E + iη + )I -H(k y ) -Σ S -Σ D ] -1 (3)
where k y is the transverse wave vector, H(k y ) is the Hamiltonian matrix discretized along the transport direction x, I is the identity matrix and Σ S /Σ D are self-energy matrices for the source/drain contacts. In this study, ballistic transport is assumed due to the relatively short channel lengths considered and hence no scattering self-energy is included in Eq.( 3). We consider a nano-sheet of Germanane and establish Bloch periodic boundary conditions along the transverse direction y. Moreover by neglecting mode mixing due to a smooth potential we implement uncoupled mode space approach for calculating NEGF equations, eventually resulting in computational savings. From the retarded Greens Function G r (k y , E) we calculate the electron and hole Greens Function G n (k y , E) and G p (k y , E) as follows:

G n (k y , E) = G r (k y , E)Σ n G r † (k y , E) G p (k y , E) = G r (k y , E)Σ p G r † (k y , E) (4)
where ,

Σ n ≡ (Σ S/D -Σ S/D † ) * f S/D Σ p ≡ (Σ S/D -Σ S/D † ) * (1 -f S/D ) (5)
with † denoting the transverse conjugate and f S/D denoting the Fermi-distribution function of source/drain. The carrier densities are calculated from the jth diagonal elements of G n (k y , E) and G p (k y , E) as

n j,j = ky 1 S +∞ -∞ G n (ky,E) 2π d(E) p j,j = ky 1 S +∞ -∞ G p (ky,E) 2π d(E) (6)
and the electron current density flowing from position j to j + 1 along the x direction is calculated from the off-diagonal elements (j, j + 1) of G n (k y , E) as

J n j,j+1 = ky 2q S +∞ -∞ d(E) 2π [H j,j+1 G n j+1,j (k y , E)- H j+1,j G n j,j+1 (k y , E)] (7)
where S is the device cross-section along the transport direction. Finite difference method has been used to discretize the 2D Poisson's equation in the cross-section of the device with the inclusion of the Dirichlet boundary condition at the metal gate electrodes and the Neumann boundary condition at the left, right, top and bottom edges. The Poisson's equation solved is:

∇ 2 U (z, x) = - ρ ( 8
)
where U is the electrostatic potential, is the dielectric constant and ρ is the charge density computed by the Recursive Green's Function Algorithm [START_REF] Lake | Single and multiband modeling of quantum electron transport through layered semiconductor devices[END_REF], [START_REF] Anantram | Modeling of Nanoscale Devices[END_REF]. The NEGF transport equation and the Poisson equation are solved iteratively until self-consistency is achieved. Mixed C and MATLAB R [28] programming was implemented to speed up the Recursive Green's Function algorithm [START_REF] Jiang | Quantum transport in RTD and atomistic modeling of nanostructures[END_REF]. Parallel computation for the transverse modes carried out on a 3.6 GHz, 6 core, CentOS workstation resulted in further reduction of the computational time.

III. RESULTS AND DISCUSSIONS

We first investigate the static current voltage characteristics of monolayer Germanane FETs followed by dynamic performance metrics. Transfer characteristics of both nMOSFET and p-MOSFET have been examined and the simulated I D -V G plots of p-MOSFET and n-MOSFET for varying channel lengths are shown in Fig. 3. The I ON in monolayer Germanane n-MOSFET for L CH =7 nm, 5 nm and 3 nm are 888.4 µA/µm, 1.7 mA/µm and 2.4 mA/µm respectively , whereas that for p-MOSFET are 697 nA/µm, 1.44 mA/µm and 4.16 mA/µm respectively as shown in Table I. Therefore it is quite evident that devices with channel length scaled upto 3 nm show comparable ON-current as that of ITRS 2028 HP metric. However, closer observation of the plots reveal an interesting trend. As the channel length increases the p-MOSFET starts performing better and exceeds the current in n-MOSFET. This phenomenon can be explained after studying the energy resolved current density spectrum plots in Fig. 4 to Fig. 6. The current spectrum plots for all the channel lengths of p-MOSFET and n-MOSFET have been calculated at that bias voltage where we see the current in p-MOSFET starting to cross over the current in n-MOSFET which is V G =0.4V, V D =0.6V. We also calculate the percentage of tunneling current I TUNNELING and current over the barrier I THERMIONIC for each of the channel lengths in both n-MOSFET and p-MOSFET at the said bias point as shown in Table II. We find that in shorter channel length devices the percentage of tunneling is more than the percentage of thermionic emission, thereby resulting in degraded slope of I D -V G curves. As the channel length increases, the percentage of tunneling reduces in both n and p-MOSFET, attributing to improvement in the slope of the transfer characteristics. However, we see that the percentage reduction in tunneling in p-MOSFET from 5 nm to 7 nm is more pronounced that in case of n-MOSFET. Due to this the improvement of slope from shorter to higher channel lengths, is better in p-MOSFET than its n-counterpart and hence the on-current starts to exceed at higher gate overdrives. Further, the modulation of barrier height in both the devices play an important role in determining the current cross-over phenomenon. Studying Fig. 4 to 6 we find that, in shorter channel length devices the height of the barrier is such that the current is mostly tunneling. Therefore, the n-MOSFET which has a lower effective mass results in a higher current value than the p-MOSFET. But, as the channel length increases, say for example 7 nm, when the barrier height is very low in both the devices at high gate overdrives, the current is mostly due to thermionic emission. Here, the p-MOSFET having a higher effective mass attributing to higher denisty of states has more number of modes to fill and hence exhibits a higher current than n-MOSFET.

Next, the dynamic HP parameters like intrinsic switching delay (τ ), total Gate to ground capacitance (C G ) and Power Delay Product (PDP) are calculated for all the channel lengths of n-MOSFET and p-MOSFET and compared. The values are shown in Table I. The intrinsic device delay, is computed following the method explained in [START_REF] Koswatta | Performance comparison between p-i-n tunneling transistors and conventional MOS-FETs[END_REF] as τ = (Q ON -Q OFF )/I ON ,where Q ON and Q OFF are the overall charges induced in the device at the ON-state and OFF-state, respectively. Total Gate to ground capacitance has been calculated by differentiating the total charge along the entire length of the device with respect to gate voltage as C G = δQ/δV GS . The PDP is calculated as V DD (Q ON -Q OFF ) [START_REF] Logoteta | Graphene-based lateral heterostructure transistors exhibit better intrinsic performance than graphene-based vertical transistors as post-CMOS devices[END_REF]. In Fig. 7 we show the C G values with respect to gate bias in case of both n-and p-MOSFETS. It is seen that in all cases C G is higher in case of p-MOSFET, the reason being that it has a higher density of states than n-MOSFET and hence a larger quantum capacitance translating to a higher C G . The ultrascaled n-and p-MOSFETs where the transport is mainly due to source to drain tunneling, do not exhibit much difference in the magnitudes of C G . However, as the channel length increases and the devices entire the regime of higher gate overdrives ( V D ∼= 0.3V to 0.4V ), there is a significant change in the magnitude of the C G values in case of ptype Germanane FETs. The reason is that, at such high gate overdrives, the barrier is very low in case of both n-and p-type FETs resulting mostly in thermionic emission, the p-MOSFET having higher density of states uncovers more number of modes thereby resulting in more modulation of charge with increment of the gate bias. The intrinsic switching delay, which is a manifestation of the gate capacitance, exhibits lower values in case of all the channel lengths of n-MOSFET than p-MOSFET. However, intrinsic switching delays for both p- and n-MOSFETs are comparable to the ITRS 2028 values. It should be noted , that the values of intrinsic switching delay in this sub-deca nanometer Germanane MOSFETS, tend to be increasing with the scaling of the channel length. The reason is that, though the difference in magnitudes of Q ON and Q OFF decreases as manifested by the reduction in values of C G , the on-current value also decreases, and hence the overall ratio of (Q ON -Q OFF ) to I ON increases thereby translating to an increment in the values of τ .

Next, we compare the energy efficiency and the switching capability of Germanane monolayer FETs with respect to other 2D material FETs. CMOS HP [START_REF] Anantram | Modeling of Nanoscale Devices[END_REF] BP AD 7.3nm [28] BP ZD 7.3nm [28] MoS 2 HP 5.1nm [29] MoS 2 HP 8.9nm [START_REF] Jiang | Quantum transport in RTD and atomistic modeling of nanostructures[END_REF] GeH nMOS GeH pMOS III-V HP [START_REF] Low Kain | Ballistic Transport Performance of Silicane and Germanane Transistors[END_REF] Ge HP CMOS HP [START_REF] Anantram | Modeling of Nanoscale Devices[END_REF] BP AD 7.3nm [28] BP ZD 7.3nm [28] MoS 2 HP 5.1nm [29] MoS 2 HP 8.9 nm [START_REF] Jiang | Quantum transport in RTD and atomistic modeling of nanostructures[END_REF] GeH nMOS GeH pMOS III-V HP [19] Ge HP [19] ITRS HP In Figs. 8 and9 we have plotted the switching energy versus the gate delay, and gate-delay versus I ON / I OFF of these devices respectively. In the plots, CMOS HP stands for 15nm technology node CMOS HP at 0.73V [START_REF] Majumdar | Benchmarking Transition Metal Dichalcogenide MOSFET in the Ultimate Physical Scaling Limit[END_REF]. BP AD and ZD stand for HP monolayer Black Phosphorus FETS in Armchair and Zigzag directions resp.,where the channel length is 7.3 nm with V D =0.69V [START_REF] Liu | Ballistic transport in monolayer black phosphorus transistors[END_REF]. Further, MoS 2 HP FETs for two different channel lengths are shown in the charts [START_REF] Liu | On Monolayer MoS2 Field-Effect Transistors at the Scaling Limit[END_REF]. The ITRS predicted trend lines till the end of the roadmap, as well as III-V/Ge HP technologies are also shown in the plots for reference. The straight lines in Fig. 8 and Fig. 9 denote constant energy-delay products and constant ratios of switching energy to I OFF resp. In both the charts, it can been seen that Germanane FETs meet the ITRS HP as well as III-V/Ge HP FET requirements very well. From Fig. 8 it is seen that the ultra scaled Germanane n-MOSFET has faster operation with a lower energy consumption with respect to the other 2D FETs. 7 nm Germanane n-MOSFET has almost similar performance as BP ZD FET , with lower energy consumption per switching event. On the other hand, Germanane p-MOSFETs has greater switching energy than their n-counterparts, but show comparable performance with respect to other 2D FETs. In Fig. 9 we see that the performance is best at higher channel lengths for both n-and p-Germanane FETs. Though they are not able to outperform the BP FETs in terms of faster switching and higher I ON /I OFF ratios, their performance are better compared to the MoS 2 HP FETs. Also, CMOS HP has almost the same I ON /I OFF ratio as the 7 nm Germanane FETs but at the cost of gate delay higher than an order of magnitude. Therefore as a whole, ultra-scaled Germanane FETs deliver a faster response at low switching energies, almost consistent energy-delay product, and practical I ON /I OFF ratios with uniform E/I OFF ratio -all of which combine to deliver a more energy efficient computation with a smaller footprint.

IV. CONCLUSION

The high performance characteristics of Germanane FETs in the sub-deca nanometer regime seem to be well comparable to the ITRS 2028 predicted values. Both Germanane n-MOSFET and p-MOSFET offer excellent I ON with large I ON /I OFF ratios even when scaled down to 3 nm channel length, which is beyond the ITRS projected limits. Dynamic parameters such as total capacitance, intrinsic switching delay and switching energy have also been calculated ans we show that n-MOSFET performs better in comparison to the p-MOSFET, though both the devices meet the ITRS projected limits. Our study therefore suggests that Germanane FETs possess benefits for high-performance device applications even in the sub-deca nm regime.

Fig. 1 :

 1 Fig. 1: Cross section of the device structure. t ox = 3 nm. L S (L D )= 30 nm in simulation. The length of the intrinsic channel is kept equal to the length of gate L G .

Fig. 4 :Fig. 5 :Fig. 6 :

 456 Fig. 4: Energy resolved current spectrum superimposed on the conduction band profiles for both (a) Germanane n-MOSFET and (b) p-MOSFET with channel length 3 nm

Fig. 3 :

 3 Fig. 3: Transfer characteristics of GermananeFETs for (a) L CH =3 nm (b)L CH =5 nm and (c)L CH =7 nm and for a drain bias of V D = 0.6 V.

Fig. 7 :

 7 Fig. 7: Gate capacitance as a function of gate bias for (a) L CH =3 nm, nMOSFET and p-MOSFET, (b)L CH =5 nm, nMOSFET and p-MOSFET and (c)L CH =7 nm, n-MOSFET and p-MOSFET.

Fig. 8 :

 8 Fig.8: Benchmarking the energy Vs delay with respect to the state of the art 2D materials and projected ITRS limits. E.τ denotes the energy delay product.

Fig. 9 :

 9 Fig. 9: Benchmarking the delay Vs I ON / I OFF with respect to the state of the art 2D materials and projected ITRS limits.

[k] -1.8 -1.6 -1.4 -1.2 -1 -0.8 Energy [eV] DFT k.p

		• respectively [18]. The k.p
	calculations are in agreement with DFT results for a range of
	∼ 0.2 eV, as can be seen in Fig. 2. It should be noted that
	since the Fermi level in our simulated p-MOSFET structure
	is below the valence band maximum level by 0.03 eV, the
	contributions from two valence bands can be readily captured
	by the proposed framework.
	-0.05	0	0.05
	Fig. 2: Band-structure close to the Γ point of Germanane chair
	valence band maximum. The 2-band k.p Hamiltonian shows
	a close agreement with the DFT calculations.

TABLE I

 I

	: Performance Metrics comparison of Germanane
	Monolayer FETs with ITRS 2028 Prediction
	Parameter	Germanane nMOS Germanane pMOS	2028 HP
	L CH (nm)	3	5	7	3	5	7 4.1
	I ON (mA/µm) 0.89 1.37 2.47 0.70 1.44 4.16 0.90
	ION IOFF (×10 3 )	8.90 13.70	25 7.00 14.40 41.6	9
	C G (fF/µm) 0.12 0.16 0.22 0.14 0.51 0.90 0.60
	τ (fs)	72	57	39 107	69	47 423
	P DP (aJ/µm) 38	47	58 45	82	118 600
	TABLE II: Percentage of tunneling component and thermionic
	component of the total current in Germanane Monolayer FETs
	for all channel lengths at V G =0.4V, V D =0.6V
	Parameter	Germanane nMOS Germanane pMOS
	L CH (nm)	3	5	7	3	5		7
	I TUNNELING (%) 60.80 27.93 14.30 76.89 49.68 28.15
	I THERMIONIC (%) 37.96 72.70 85.70 23.11 50.32 71.85

This work was supported by Science and Engineering Research Board (SERB), Department of Science and Technology (DST), Government of India, under Grant SB/S3/EECE/0209/2015, Robert Bosch Cyber Physical System, Indian Institute of Science (IISc) Bangalore, under grant no: RBC00011 and Raman-Charpak Fellowship, Centre Franco-Indien pour la Promotion de la Recherche Avancée (CEFIPRA) under grant no. RCF-IN-133. M. Brahma is with the Centre for Nano Science and Engineering (CeNSE) and S.


Alexandre Lhosmot

G Bertrand

M Steinmann

M L Toussaint

C Bertrand

L Collin

A Ponçot

G Magnon

H Jeanton

A V Walter

published or not. The documents may come   L'archive ouverte pluridisciplinaire

The CRITICAL PEAT Project : The importance of hydrology for carbon reactivity along with atmosphere -peatland interactions. Preliminary results from the Frasne peatland observatory (Jura Mountains, France).

Peatlands: a threatened 3 pillar equilibrium ~3 % of emerged surfaces (Joosten and Clarke, 2002) ~30 % of soil carbon (Parish et al, 2008), an underestimated stock? A diversified peatland patchwork: ecological mosaic (vegetation, peat fiber rate) (Briot, 2004;Grosvernier, 2005), variable peat thickness (Collin, 2016).

A large Jurassic syncline of the Jura arc: karstified limestones overlain by moraines. Mountainous area at 800 to 1100 m asl. Frasne peatland at 834 m asl.

Temperature increases significantly: + 1 °C in 10 years.

→ C sink function (NEE <0) partially dependent on the hydrological conditions. In particular, during the maximum metabolism period (ET>>0), great WTD (below -0,20 m) seems to be associated with a source function (R ECO >|GPP|). In contrast, when WTD is closer to the surface, |GPP|> R ECO ).

❑ Geological framework

❑ Pedological framework A Hydro-geochemical mosaic?

These heterogeneities could favor a patchwork-like hydro-ecological functioning.

1 LHOSMOT A, 1 BERTRAND G, 1 STEINMANN M, 2 JACOTOT A, 1 PONCOT A , 3 COLLIN L., 1 TOUSSAINT ML and the CRITICAL PEAT Team Local precipitation?

Evaporation Regional GW?

Seasonal and multi-year trend correlations

• Deepest water (< -50 cm) close to the LMWL (1200 m).

• Enriched superficial water (> -50 cm) suggesting evaporation.

Similar patterns for WTD and precipitation: the multiyear WTD follows precipitation. Buffered impact of the restoration at the center of the peat.

Precipitation (P) pattern: interannual variability and slight increase since 2008.

The Frasne peatland is representative of a sensitive compartment of the Critical Zone combining water, carbon and biodiversity stakes.

→ 2 hypotheses : (1) Water/moraine interaction may vary laterally due to variable residence times, leading to contrasted EC. (2) Alternatively, punctual injection of karstic water at the bottom of the peatland.

  (Nichols et al., 2019) ~10 % of global freshwater resources(Joosten and Clarke, 2002)

1•

  Bourgogne Franche-Comté University, Chrono Environnement UMR UFC-CNRS 6249, Besançon, France.; 2 ISTO, University of Orléans, France; 3 Communauté de communes Frasne Drugeon. Contact : alexandre.lhosmot@univ-fcomte.fr The peatland was a CO 2 sink in summer 2019, but a source of CO 2 during the rest of the time.• C balance of the peatland since July 2008: +167,8 gC/m².Outlet flow monitoring• Low EC during high flow (winter), high EC in summer. → Mobilization of strongly mineralized deep peat water or input of GW (moraines, karst) during dry periods.

Materials and methods

Bibliography


Elie Grinfeder

Sylvain Haupert

Manon Ducrettet

Julien Barlet

Marie-Pierre Reynet

Frédéric Sèbe

Jérôme Sueur
email: sueur@mnhn.fr

Soundscape dynamics of a cold protected forest: dominance of aircraft noise

Keywords: soundscape, cold forest, automatic classification, anthropophony, climate, protected area

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

Introduction

Due to anthropization, land use and climate change may alter all types of landscape through mixed effects on plant and animal population ecology [START_REF] Opdam | Climate change meets habitat fragmentation: linking landscape and biogeographical scale levels in research and conservation[END_REF][START_REF] Opdam | Changing landscapes to accommodate for climate change impacts: a call for landscape ecology[END_REF][START_REF] Jennings | Climate change and ecosystem composition across large landscapes[END_REF]. The analysis, prediction and possible mitigation of these effects require appropriate monitoring methods that can retrieve relevant ecological information on large spatial scales and over long time periods. Long-term research is a requisite to estimate past, current and future changes (Hugues et al. 2017) and to develop general theories in ecology and evolutionary biology [START_REF] Kuebbing | Long-term research in ecology and evolution: a survey of challenges and opportunities[END_REF]. Ecoacoustics, more specifically soundscape ecology which works at the landscape scale, offers valuable tools to retrieve ecological information with remote field-based sensors on large spatiotemporal scales [START_REF] Pijanowski | What is soundscape ecology? An introduction and overview of an emerging new science[END_REF][START_REF] Sueur | Ecoacoustics: the Ecological Investigation and Interpretation of Environmental Sound[END_REF]. The use of autonomous recorders led to a long series of recent papers which, among others, showed that ecoacoustics can be involved in a wide range of ecological and conservation applications [START_REF] Sugai | Terrestrial Passive Acoustic Monitoring: Review and Perspectives[END_REF][START_REF] Sugai | A roadmap for survey designs in terrestrial acoustic monitoring[END_REF]. Soundscape ecology and ecoacoustics have both been put forward as interesting approaches for the study of the impacts of climate change [START_REF] Krause | Using ecoacoustic methods to survey the impacts of climate change on biodiversity[END_REF][START_REF] Sueur | Climate change Is breaking Earth's beat[END_REF].

A mainstay of soundscape ecology is the understanding of the acoustic patterns that compose a soundscape. These patterns can be separated into three categories: the biophony (biotic sounds), the geophony (abiotic sounds), and the anthropophony (human related sounds) [START_REF] Krause | Anatomy of the soundscape: evolving perspectives[END_REF]. A sub-classification might be considered by including the technophony as a sub-component of anthropophony which refers to the sound of human technology [START_REF] Mullet | Temporal and spatial variation of a winter soundscape in south-central Alaska[END_REF]. Biophony, geophony and anthropophony are repeatedly referred to as major structuring components of soundscapes [START_REF] Krause | Anatomy of the soundscape: evolving perspectives[END_REF][START_REF] Farina | Soundscape ecology: principles, patterns, methods and applications[END_REF].

Biophony contains all sounds produced by living organisms, whether they are intentional (such as animal communication signals) or incidental (such as leaves rustling). Biophonic sounds can, therefore, work as proxies of life traits, ecological processes, and biodiversity [START_REF] Gibb | Emerging opportunities and challenges for passive acoustics in ecological assessment and monitoring[END_REF]. Acoustic biodiversity can be evaluated with two major types of indicators: acoustic richness and acoustic abundance (Farina 2018). Although geophonic sounds are less studied than anthropogenic and biophonic sounds, they are a major soundscape component, as found in any type of environment. Such sounds can be a source of noise or information that can alter species behavior in the short term [START_REF] Farji-Brener | Working in the rain? Why leaf-cutting ants stop foraging when it's raining[END_REF][START_REF] Geipel | Noise as an informational cue for decision-making: the sound of rain delays bat emergence[END_REF][START_REF] Lengagne | The effects of rain on acoustic communication : tawny owls have good reason for calling less in wet weather[END_REF][START_REF] Tishechkin | Vibrational background noise in herbaceous plants and its impact on acoustic communication of small Auchenorrhyncha and Psyllinea (Homoptera)[END_REF]) and the long term (Brumm 2004, Zollinger and[START_REF] Zollinger | Why birds sing loud songs and why they sometimes don't[END_REF]. Anthropogenic noise due to urbanization, transportation, industry, recreation, and military activity is invading all possible landscapes (marine, aquatic and terrestrial) even when they are protected by local or national legislation. If the negative effects of anthropogenic noise have been thoroughly studied at the individual level through animal behavior and animal physiology studies [START_REF] Brumm | Animal communication and noise[END_REF][START_REF] Slabbekoorn | Birdsong and anthropogenic noise: implications and applications for conservation[END_REF][START_REF] Barber | The costs of chronic noise exposure for terrestrial organisms[END_REF]Shannon 2015), less research has been conducted at larger ecological scales. Anthropophony has been shown to be a prevalent part of preserved soundscapes [START_REF] Mullet | Temporal and spatial variation of a winter soundscape in south-central Alaska[END_REF]Francis et al. 2011aFrancis et al. , 2011b)), and to have a role in the structure of soundscapes, in relation to massive noise generators, such as industrial sites [START_REF] Duarte | The impact of noise from open-cast mining on Atlantic forest biophony[END_REF][START_REF] Deichmann | Soundscape analysis and acoustic monitoring document impacts of natural gas exploration on biodiversity in a tropical forest[END_REF] or highways [START_REF] Munro | Traffic noise impacts on urban forest soundscapes in south-eastern Australia[END_REF][START_REF] Khanaposhtani | Effects of highways on bird distribution and soundscape diversity around Aldo Leopold's shack in Baraboo[END_REF].

Anthropogenic noise can alter the benefits that humans may experience from natural soundscapes [START_REF] Francis | Acoustic environments matter: Synergistic benefits to humans and ecological communities[END_REF]. Anthropogenic noise can also, through cascading effects, induce changes in singing species communities by altering species interactions [START_REF] Francis | Noise pollution changes avian communities and species interactions[END_REF], disturbing trophic chains [START_REF] Hanache | Noiseinduced reduction in the attack rate of a planktivorous freshwater fish revealed by functional response analysis[END_REF], and disrupting pollination and seed dispersal by modifying animal communities [START_REF] Francis | Noise pollution alters ecological services: enhanced pollination and disrupted seed dispersal[END_REF]. Because biophony, geophony and anthropophony are part of a single ecological system, it appears necessary to understand their cross-correlated patterns over the long term.

Therefore, it seems essential, as a very first step, to automatically assess the presence of biophony, geophony, and anthropophony in long soundscape monitoring protocols. This decomposition appears as a challenge since biotic, abiotic and human sound sources can display great diversity regarding amplitude, time and frequency pattern, and moreover, they often overlap and interfere with each other. So far, such decomposition was operated either manually through listening and spectrographic visualization [START_REF] Matsinos | Spatio-temporal variability in human and natural sounds in a rural landscape[END_REF][START_REF] Liu | Spatiotemporal variability of soundscapes in a multiple functional urban area[END_REF][START_REF] Duarte | The impact of noise from open-cast mining on Atlantic forest biophony[END_REF][START_REF] Mullet | Temporal and spatial variation of a winter soundscape in south-central Alaska[END_REF][START_REF] Gasc | Soundscapes reveal disturbance impacts: biophonic response to wildfire in the Sonoran Desert Sky Islands[END_REF][START_REF] Rountree | Temperate freshwater soundscapes: a cacophony of undescribed biological sounds now threatened by anthropogenic noise[END_REF] or on the basis of frequency delimitation, anthropophony mainly occurring between 1 and 2 kHz [START_REF] Joo | Analysis and interpretation of variability in soundscapes along an urban rural gradient[END_REF][START_REF] Gage | Visualization of temporal change in soundscape power of a Michigan lake habitat over a 4-year period[END_REF][START_REF] Fairbrass | Biases of acoustic indices measuring biodiversity in urban areas[END_REF][START_REF] Ross | Listening to ecosystems: data-rich acoustic monitoring through landscape-scale sensor networks[END_REF][START_REF] Doser | Characterizing functional relationships between anthropogenic and biological sounds: a western New York state soundscape case study[END_REF].

However, the choice of frequency limits is questionable as animals can produce sound at a lower or higher frequency range [START_REF] Kasten | The remote environmental assessment laboratory's acoustic library: An archive for studying soundscape ecology[END_REF]. The biophony/anthropophony ratio has been quantified through spectrographic decomposition but neither geophony nor silence were estimated [START_REF] Dein | Landscape influence on biophony in an urban environment in the European Alps[END_REF]. Eventually, artificial intelligence appears as a promising solution to tease apart the soundscape components. Convolutional neural networks (CNN) have indeed been involved in acoustic biodiversity monitoring [START_REF] Sethi | Characterizing soundscapes across diverse ecosystems using a universal acoustic feature set[END_REF], soundscape labelling [START_REF] Bellisario | Contributions of MIR to soundscape ecology. Part 3: Tagging and classifying audio features using a multi-labeling k-nearest neighbor approach[END_REF]), species identification (LeBien et al. 2020) and sound source separation [START_REF] Lin | Source separation in ecoacoustics: A roadmap towards versatile soundscape information retrieval[END_REF]). Yet, most of these recent tools require important efforts as most of them are based on complex algorithms, large training datasets and costly hardware components.

In Western Europe, mountain landscapes may suffer important land use and climate change at the same time. The Haut-Jura region, a middle altitude and cold area across France and Switzerland, is a typical example of a mountain landscape that could be substantially altered by climate change. The Haut-Jura region is dominated by a typical coniferous forest landscape inhabited by a very rich diversity of resident and migrating bird species. The region also includes the city of Geneva, the second largest urban area of Switzerland, which is served by one of Europe's major airports. The region is therefore a mix of cold, preserved natural environment and urbanized zones.

To estimate the possible changes in soundscape composition over the long-term, in relation to land use and climate change, we started a 15-year monitoring program in one of the largest continuous coniferous forests in Europe, localized in the Parc naturel Régional du Haut-Jura and protected by several conservation tools. Using the first year of audio data, encompassing 2,336 hours of sound files, we first aim at understanding the dynamics of the soundscape categories over the seasons. We trained an artificial neural network, based on cepstral features, to automatically classify biophony (birds, mammals and insects), geophony (wind and rain), anthropophony (aircraft) and silence. The annual and diel phenology of each class could then be estimated, serving as a baseline for further research in the long-term and revealing for the first time the predominance of aircraft noise competing with biophony and limiting silence in a protected natural area.

Material and Methods

Study site

The site consisted of a temperate cold climax forest locally named 'Risoux'. The forest is located in the East of France in the Jura Mountains at the Swiss border and the closest city is Morez (46°31'22'' N, 6°01'23'' E). The forest consists of a 4,400 ha West-East flat anticline crossed by 26 km of roads, hiking and cross-country skiing trails. The climate is semicontinental. The anticline and the average altitude of 1,230 m a.s.l. fosters cold air conservation, inducing low annual average temperature (5.5°C), important snowfalls (> 2 m), and an extended winter period with a very narrow vegetation period between April and October. Vegetation is mainly composed of European spruce (Picea abies), European beech (Fagus sylvatica), and European silver fir (Abies alba), their relative density depending on the altitude. Due to rich native flora and fauna, this old forest is protected by a natural park (Parc Naturel Régional du Haut-Jura) and a European Natura 2000 zone network, since 2003.

Conservation measures and legal protection aim at ensuring the long-term protection of Europe's most valuable and threatened species and habitats, listed under both the European Birds Directive and the Habitats Directive. However, these conservation tools do not exclude human activities so that the Risoux Forest is still exploited for recreational hunting and commercial logging.

Recording

The forest soundscape was recorded using four automatic SongMeter 4 recorders (Wildlife Acoustics Inc, Concord, MA, USA). These recorders were installed across the Risoux Forest central zone, along a West-East axis. The distance between the recorders was 1.00 ± 0.10 km.

The position of the recorders was chosen to (1) cover the forest area, (2) sample a single habitat defined by the dominance of the European spruce, (3) avoid pseudo-replication between neighboring recorders, and (4) avoid anthropogenic noise due to car traffic and local activities (Figure 1).

The recorders were attached to trees at a height of 2.5 m and oriented on a 45-315° axis so that the microphones were parallel to the main circulation axes. The microphones were protected with a plastic roof to reduce rain noise and damage and were facing south to limit ice and moisture. The recorders were powered with an external 12 V battery and were programmed to record for 1 min every 15 min (1' on, 14' off, 96 recordings per day) all year round, in accordance with reviewed literature [START_REF] Depraetere | Monitoring animal diversity using acoustic indices: Implementation in a temperate woodland[END_REF][START_REF] Bradfer-Lawrence | Guidelines for the use of acoustic indices in environmental research[END_REF]. The recordings started on July 13 2018 and were programmed to work over 15 years. Sounds were recorded in lossless stereo .wav format with a 44.1 kHz sampling frequency, 16bit digitization depth, no analog filtering and a total gain of 42 dB. Mono and right channels were considered as duplicate sampling units, so that only the left channel was considered.

The sound database recording, covering the first monitoring year, began on August 1st and stopped on July 31st 2019. This recording method resulted in a yearly database of 140,160 sound files (96 x 365 days x 4 sites) for a total of 2,336 hours of recording. Due to some recorder failures, the effective database contained 137,087 sound files for a total of 2,285 hours of effective recording.

Supervised sound classification

In order to assess variations in the main forest soundscape components over space and time, a supervised sound classification was achieved. Five sound types were classified: aircraft (anthropophony), wind (geophony), rain (geophony), biophony, and silence (Figure 2).

Anthropophony was dominated by aircraft, as the recording sites were chosen to be distant from other human activities (ground based transportation, hiking, skiing, hunting and logging). Biophony was defined as any sound attributed to animal vocalization or movement.

Silence was defined as the absence of any emergent sound. Classification was based on an artificial neural network (ANN) trained on a subsample of manually labelled sounds. The subsample contained 1% (i.e. 1,402 files) of all recordings and was built to ensure a balanced dataset that is representative of (1) the four recording sites, (2) the 12 months of the year, (3) the seven week days and (4) the 24 hours of a day-and-night cycle. . Each recording was then listened to and manually classified as one of the five sound types by one of us (EG) using Audacity v2.3.2 (Audacity team 2019). The files were listened to through high-quality headphones with a noise reduction option (Bose Quiet Comfort 35 II, Bose, USA) and visualized with a grey-scale spectrogram with a Fourier window of 1024 samples with no overlap and a dynamic range of 60 dB.

The recordings were parameterized with Mel frequency cepstral coefficients (MFCCs), obtained with the R package tuneR v1.3.3 [START_REF] Ligges | tuneR: analysis of music and speech[END_REF]. Successive MFCCs were computed over time, using a moving window of 512 samples with no overlap between successive windows. The time precision (δ t = 0.0116 s) was enough to track fast sound events. This computation resulted in a matrix of 86 columns (windows) and 26 rows (MFCCs) for each 1 min audio recording. The matrix was then compressed by computing the average according to the rows and the 26 MFCCs so that the initial recording was parameterized with a vector of 26 coefficients [START_REF] Sueur | Sound analysis and synthesis with R[END_REF]. These coefficients were the input data used by the ANN multi-label classification algorithm using the R package neuralnet v1.44.2 [START_REF] Fritsch | neuralnet: training of neural networks[END_REF]. The settings of the ANN algorithm were: (1) an input layer of 26 neurons corresponding to the 26 MFCCs, (2) a hidden layer of five hidden neurons, (3) an output layer of five neurons corresponding to the five sound types, and (4) a stopping criteria (threshold for the partial derivatives of the error function) of 0.4. The algorithm used was a resilient backpropagation with weight backtracking and the activation function was a logistic function. Due to the random start of the algorithm, the ANN procedure was repeated so that true positive rate (TPR) was maximized and false positive rate (FPR) minimized. The output of the ANN was the probability of the presence for each sound type. The optimized decision threshold for each sound type was automatically determined using the inflection point of the respective receiver operating characteristic (ROC) curve. The final selected ANN was then used to classify the remaining dataset (140,160-1,402 = 138,758 files). Consequently, the final output was a presence/absence (0/1) decision for each sound type and for each recording, so that a presence/absence matrix of 140,160 rows (recordings) and five columns (sound types) was obtained.

Statistical analysis

The presence/absence matrix obtained with the ANN was first treated with an ordination analysis to assess the role of temporal and spatial factors on sound type classification. A Multiple Correspondence Analysis (MCA) was applied to the matrix using sound type absence/presence as the variables to be explained and (1) recording sites, (2) months, (3) week days and (4) hours as explanatory variables. MCA results were visualized by plotting the projections of both explained and explanatory variables on the two first axes. The analysis was run with the R package FactoMineR v1.34 [START_REF] Lê | FactoMineR: An R package for multivariate analysis[END_REF].

Spatio-temporal trends of sound types were also assessed with generalized additive models and generalized additive mixed models. These models accounted for nonlinear relationships between variables to be explained (number of detections per hour and per month) and the explanatory variables, i.e., the spatio-temporal covariates (hour and month as continuous variables) and the sites (categorical variable) [START_REF] Hastie | Generalized additive Models[END_REF]. The nonlinear relationship was modelled with smooth functions of covariates which can be isotropic and/or tensor product smooths [START_REF] Wood | Generalized additive models: an introduction with R, Second Edition[END_REF]. Tensor product smooth modelled interactions between covariates with different units. The tensor product could include both the effect of the two variables and the interaction effect [START_REF] Fewster | Analysis of population trends for farmland birds using generalized additive models[END_REF]. Model selection was based on the Akaike's Information Criterion (dAIC < 2) [START_REF] Burnham | Model selection and multimodel inference: a practical information-theoretic approach[END_REF][START_REF] Shadish | Using Generalized Additive (Mixed) Models to Analyze Single Case Designs[END_REF]. Seasonal patterns, their differences between sites and common temporal nonlinear trends were assessed with boxplot, quantile-quantile, and frequency plots [START_REF] Zuur | A protocol for data exploration to avoid common statistical problems[END_REF][START_REF] Zuur | A Protocol for Conducting and Presenting Results of Regression-Type Analyses[END_REF]. Using the function gam from the R package mgcv [START_REF] Wood | Generalized additive models: an introduction with R, Second Edition[END_REF], a full model was then fitted on data exploration with a Poisson error distribution and a log link function for each of the five variables to be explained. Collinearity among explanatory variables was assessed on a standard regression model with the variance inflated factor (VIF), applying the function vif from the R packages car [START_REF] Fox | An R companion to applied regression[END_REF].

Concurvity among smooth functions of covariates was checked using the function concurvity from the R package mgcv [START_REF] Ramsay | The effect of concurvity in generalized additive models linking mortality to ambient particulate matter[END_REF][START_REF] Wood | Generalized additive models: an introduction with R, Second Edition[END_REF]. The residual temporal dependences of each of the five models were visualized with auto-correlation functions (ACF). Since GAM models do not integrate autocorrelation error structures, models were adjusted as GAM mixed model with an adaptive autoregressive error structure using the function gamm from R package mgcv. When over-dispersion occurred, models were re-fitted with a Tweedie distribution. Pertinence of the distribution correction was evaluated through AIC comparison. Models were validated after checking residual assumptions with graphical diagnostic tools using the function gam.check from the R package mgcv and non-parametric tests compiled in the function check.residuals from the R package forecast [START_REF] Hyndman | Automatic time series forecasting: the forecast package for R[END_REF]. All statistical analyses were performed with R software version 3.6.1 (R Core Team, 2019).

Models for biophony and silence sound types were designed to establish whether the number of detections followed a site-specific pattern over time, or whether sites shared a common pattern [START_REF] Pedersen | Hierarchical generalized additive models in ecology: an introduction with mgcv[END_REF]. Three models were built. Model 1 (M1) tested hourly and monthly effects with a common trend over time for the four sites, using a tensor product smooth function with a cyclic cubic spline (hour) and a thin plate regression spline (month).

Model 2 (M2) tested hours, months and sites with a common trend over time and speciesspecific random deviations around the common trend. It is specified with the same structure as M1 plus another three dimensional tensor product, which integrates sites as a random smooth. Model 3 (M3) tested hours, months and sites using a common trend and a sitespecific trend over time, using M1 structure plus a spatio-temporal two dimensional tensor product for each site. The mathematical expression of each is detailed in Supplementary Information.

Regarding biophony, model M1 fitted with a second-order autoregressive structure did not show any residuals issue. Model M2 violated the homoscedasticity hypothesis, even after re-fitting the model with a second-order autoregressive structure in link with the low number of levels of the variable site. Model M3, applied to biophony, did not converge.

Regarding silence, model M1 was adjusted for temporal dependence with a first-order autoregressive structure and fitted with a Tweedie distribution with a square root link and a p parameter of 1.01. Model M3 did not converge, whereas model M1 had a lower AIC than model M2.

Regarding wind, rain and aircraft sound types, the site effect was not considered due to a risk of pseudo-replication. Model M1 was fitted with a first-order autoregressive structure for aircraft and wind. A Tweedie distribution with a logarithm link and a p parameter of 1.06 was specified for the wind sound type model. For rain sound type, model M1 showed residual patterns, as well as a model with an hourly smooth by monthly trend.

Results

Automatic sound type classification

The area under the curve (AUC) of the ROC curve, when training the artificial neural, network was 0.92 for aircraft, 0.93 for wind, 0.88 for rain, 0.85 for biophony, and 0.95 for silence. The respective false positive rate (FPR) and true positive rate (TPR) were 3.75% and 71.66% for aircraft, 13.56% and 72.51% for biophony, 12.99% and 87.05% for wind, 11.11% and 78.26% for rain, and 10.65% and 94.12% for silence. The application of the trained artificial neural network on the total database led to a detection of 102,195 files with anthropophony (75%), 58,348 files with biophony (43%), 49,287 files with geophony (57%), 49,287 files with wind (36%), 28,546 files with rain (21%)), and 18,890 files with silence (14%).

The percentage of the number of detections per file per month and per hour reveals specific patterns for each sound type (Figures 34). Aircraft was nearly constant across months and less abundant during the night than during the day, with a maximum around noon and a minimum at 3 am. Wind was almost evenly distributed across the year at the exception of a peak in March and showed a broad peak around 4 pm. Rain showed a maximum value in December and a minimum value in January, and was evenly distributed across day and hours.

Biophony followed a clear seasonal pattern with a peak in June corresponding to spring and a trough in January corresponding to winter. Nevertheless, silence did not show any particular seasonal pattern but a clear pattern according to days and hours with less silence during the day than during the night and a peak at 3 am. This confirmed that silence patterns were opposite to aircraft sounds as suggested by the MCA.

Ordination analysis

The first two axes of the Multiple Correspondence Analysis (MCA) explained 65.13% of the total variance. The third and fourth axes explained 28.78% additional variance but did not reveal interesting patterns. The projection of the variables to be explained (aircraft, biophony, wind, rain, silence) and the explanatory variables (sites, months, week days and hours) along these two axes revealed that the variables sites and week days had a weight close to zero (Figure 5a). A second MCA was therefore calculated excluding these variables. The first two axes of this second MCA explained 65.40% of the total variance.

The projection of the variables to be explained along the two first axes showed that silence and aircraft contributed the most to the MCA first dimension axis, whereas rain and wind strongly contributed to the MCA second dimension axis (Figure 5b). Silence and aircraft were negatively correlated, meaning there was silence when there were no aircraft.

Rain and wind were positively correlated, meaning there was rain when there was wind. The projection of the explanatory variables revealed that the hours categories could be separated into two groups: (1) the night hours between 22:00 and 06:00 were located close to silence presence and aircraft absence, and (2) the day hours between 07:00 and 21:00 were located close to silence absence and aircraft, wind and rain presence. The months categories could be split into three groups: (1) March and December were located close to wind and rain presence, (2) August and September were close to wind and rain absence, and (3) other months were close to the origin.

Regression analysis

The phenology of each sound type shows a clear pattern according to months and hours (Figure 6). The model for aircraft sound type explained 93.7% of the variance, with a tensor product term significant at 169.3 estimated degrees of freedom, suggesting a high nonlinear relation between temporal covariates and the number of detections (Figure 6). Aircraft were characterized by a clear diel pattern with a morning peak activity followed by a plateau until 22 and then a sudden decline for the four sites. The morning peak was reached at 8 am from November to March and at 7 am for the rest of the year.

The model for wind explained 72 % of the variance with a tensor product term significant at 108.6 estimated degrees of freedom. Wind followed three patterns along the day depending on the month of the year: a constant and low number of detections from December to January, and a bell shape in the afternoon during the rest of the year, from 8 am to 7 pm in February and from 7 am to 9 pm in August. The number of detections per hour is the lowest around at 40 in October, November and January; the highest being around 80 in March and July. For rain sound type, the model did not converge and could not be analyzed.

The model for biophony explained 89,7% of the variance with a tensor product term significant at 121.3 estimated degrees of freedom, suggesting a high nonlinear relation between temporal covariates and the number of detections for the four sites. Biophony detections followed two diel patterns changing monthly: (1) from March to August, the highest peak was reached in the morning, followed by a downward trend, a second peak and a sudden decline the next hour, and (2) from September to December a similar pattern was observed but without the second afternoon peak. A slight peak in the morning in November was also observed, probably because of specific autumn sounds due to birds and deer vocalizations. The time of the first peak differed per month with the earliest activity at 6 am from May to July, progressively shifting to 9 am from December to January. Biophony activity was also characterized by a yearly trend with increasing activity from February to June, and a two-fold decrease in activity from September to January.

The model for silence explained 93.6 % of the variance, with a tensor product term significant at 130.5 estimated degrees of freedom (Figure 6). Silence occurred during the night following a bell shape and then was absent through the day. The presence of silence was prevalent from 10 pm to 9 am during cold months and in July and August, during warm months.

Discussion

Monitoring and understanding natural soundscapes inthe long-term necessitates distinguishing biophony, geophony and anthropophony [START_REF] Krause | Anatomy of the soundscape: evolving perspectives[END_REF][START_REF] Farina | Soundscape ecology: principles, patterns, methods and applications[END_REF]). This step is a prerequisite to (1) track the relative importance of each component across space and time and ( 2) study individually (that is to filter out) the dynamic and composition of each component. However, this preliminary task is challenging due to the variability of each component, so much that a separation based on simple sound parameterization is not optimal.

Convolutional neural networks (CNN) offer attractive solutions but are difficult to enforce due to a need for large training datasets and long computing time. Unsupervised clustering methods solve some of the CNN limitations. They do not require annotations and specific hardware, but the classification performance depends on the right choice of the acoustics features, such as acoustics indices (Sueur et al. 2014;[START_REF] Buxton | Efficacy of extracting indices from large-scale acoustic recordings to monitor biodiversity[END_REF], in order to be able to disentangle acoustics components and create meaningful clusters [START_REF] Phillips | Revealing the ecological content of long-duration audio-recordings of the environment through clustering and visualisation[END_REF]. Here, we used a simple artificial neural network (ANN) with only a single hidden layer applied on a vector of MFCC features to decompose the soundscape of a cold protected forest over the first year of a 15-year monitoring program. The training dataset contained only a small fraction (1%) of the complete dataset for which labelling was completed in less than 24 working hours by a single person (EG). The ANN algorithm was trained in a few minutes on a standard computer and returned robust results with an area under the curve between 0.85 and 0.95. Therefore, workflow appears reliable and easy to implement on local computers using standard and open-source scientific languages (e.g. R, Python and Julia). The main drawback of such an approach was that the output returned by the ANN consisted of presence/absence data so that no abundance information of the soundscape components was available. If this lack of abundance information could limit fine scale interpretation, data absence/presence could still reveal meaningful soundscape phenology on annual and diel scales.

The annual and diel phenologies of each soundscape component showed patterns that were, indeed, in agreement with the forest structure, surrounding human activities, and local weather. In particular, the annual phenology of biophony followed vegetation phenology, increasing at the time of plant growth and flowering between March and August, i.e. in spring and summer. During this period, the diel phenology of biophony followed the activity of the species that were acoustically the most active. In particular, biophony started with sunrise and decreased with sunset at the time of dawn and dusk bird choruses [START_REF] Gil | The bird dawn chorus revisited[END_REF].

It should be noted that anthropophony, which was dominated by aircraft noise, and biophony showed convergent daily patterns all year long as aircraft followed clear diel cyclicity, corresponding to human diel activities from the nearest airport (Geneva, Switzerland). Aircraft were detected in 75% of the files, revealing noise exceedance in the soundscape of a protected high-diversity forest. Aircraft noise has already been shown to be pervasive at a large scale in the USA's preserved or remote areas [START_REF] Dumyahn | Beyond noise mitigation: managing soundscapes as common-pool resources[END_REF][START_REF] Barber | Anthropogenic noise exposure in protected natural areas: estimating the scale of ecological consequences[END_REF][START_REF] Buxton | Noise pollution is pervasive in U.S. protected Areas[END_REF][START_REF] Stinchcomb | Extensive aircraft activity impacts subsistence areas: acoustic evidence from Arctic Alaska[END_REF]. Aircraft noise impacts bird distribution, behaviour and ecology [START_REF] Sierro | European blackbirds exposed to aircraft noise advance their chorus, modify their song and spend more time singing[END_REF][START_REF] Gil | Birds living near airports advance their dawn chorus and reduce overlap with aircraft noise[END_REF][START_REF] Dominoni | Airport noise predicts song timing of European birds[END_REF][START_REF] Wolfenden | Aircraft sound exposure leads to song frequency decline and elevated aggression in wild chiffchaffs[END_REF]) and might induce cascading effects that can affect ecosystems and landscapes [START_REF] Francis | Noise pollution alters ecological services: enhanced pollination and disrupted seed dispersal[END_REF]. The growth of air transport seems to be inescapable, French civil aerial transport is estimated as doubling between 2030 and 2050 [START_REF] Courteau | Sur les perspectives d'évolution de l'aviation civile à l'horizon 2040 : préserver l'avance de la France et de l'Europe[END_REF]). The Haut-Jura is particularly characterized by a rising trend of flight activity due to a century of aircraft development. In 2019, Geneva airport welcomed 17.9 millions passengers with a 58% growth since 2009 (Genève Aéroport, 2020). All year long, daily flight activity for noncargo airports (such as Geneva airport) starts around 5 am and ends around 10 pm (Wilkend et al. 2011), at the time of bird chorus. The evaluation of local aircraft-borne anthropogenic noise and its impact on local forest ecosystems is therefore a major issue for landscape conservation. The reduction of aircraft traffic since March 2020 due to the COVID-19 pandemic drastically changed the soundscape, so much that the effect on bird activity has already been documented [START_REF] Derryberry | Singing in a silent spring: Birds respond to a half-century soundscape reversion during the COVID-19 shutdown[END_REF]. Further studies should estimate the possible acoustic interference of aircraft and biophony inside and outside lockdown periods. Such a study should be achieved by estimating sound overlap at a fine time scale and, more importantly, along the frequency spectrum. Referring to civil aircraft traffic data, it would also be possible to compare the acoustic impact of aircraft traffic according to aircraft altitude, trajectory and national and international destinations.

Wind showed a clear acoustic diel cyclicity which corresponded to the effect of daily air temperature change on wind force. On a yearly scale, wind was more present in spring during the day, with a peak before sunset, and was most present in March. This time pattern corresponded to the seasonality of 'La Bise', the local wind coming from the North or North-East. Rain was correlated with wind, indicating bad weather and, therefore, strong geophony. Unfortunately, it was not possible to fit a confident linear model with the temporal explaining variables. Rainfall modelling probably requires the addition of other environmental variables such as topography [START_REF] Ranhao | A multivariate regression model for predicting precipitation in the Daqing mountains[END_REF], air temperature, atmospheric pressure or humidity [START_REF] Jennings | Spatial variation of the rain-snow temperature threshold across the Northern Hemisphere[END_REF]. The roles of wind and rainfall in animal acoustic communication have rarely been studied but they are major components of soundscape, as they could impair animal acoustic communication [START_REF] Lengagne | The effects of rain on acoustic communication : tawny owls have good reason for calling less in wet weather[END_REF], affect sound propagation [START_REF] Priyadarshani | The impact of environmental factors in birdsong acquisition using automated recorders[END_REF], and reduce biophony. Climate change might change rain and wind regimes. In particular, there is some evidence of an emerging consensus from the different regional climate models that annual wind density might slightly increase in northern Europe by the end of this century [START_REF] Pryor | Climate change impacts on wind power generation[END_REF]. In the Jura mountains, intensity and frequency of daily precipitation is already perceptible (Sherrer et al 2016), with high rainfall during late summer afternoons [START_REF] Barton | A multi-year assessment of sub-hourly gridded precipitation for Switzerland based on a blended radar-Rain-gauge dataset[END_REF]. Moreover, winter snowfalls, which are standard in Jura, are expected to gradually be replaced by rainfalls, increasing the number of rain events in the Risoux Forest [START_REF] Beniston | Climatic change in mountain regions: a review of possible impacts[END_REF]. All of these meteorological modifications will affect the soundscape composition with more wind and rain geophony, in particular during spring when the biophony reaches a peak.

As expected, silence was occurring when biophony and aircraft were diminishing, and overall followed a diel pattern that was opposed to the other sound types. Silence was more prevalent during winter nights, when there was neither biophony nor anthropophony, and during summer nights, probably because of a decrease of the local wind 'La Bise', inducing less geophony. Silent times in protected areas appeared as rather rare events which were fully constrained by anthropophony. Although the Risoux Forest benefits from important protection and conservation measures, animals and visitors can only benefit from very short time windows of quietness.

The estimated trends of each soundscape component illustrate a yearly acoustic profile with a winter soundscape dominated by aircraft, progressively evolving towards a more diverse soundscape, including biophony from early spring until the end of summer. The Risoux Forest soundscape is, therefore, closely tied to geophony and anthropophony. Here, we only monitored aircraft sounds that were dominating the anthropophony, but other sounds, even as minor as those of gunshots or chain saws due to hunting and logging, should not be totally neglected and could also be detected through acoustic methods [START_REF] Hill | AudioMoth: evaluation of a smart open acoustic device for monitoring biodiversity and the environment[END_REF][START_REF] Sethi | Characterizing soundscapes across diverse ecosystems using a universal acoustic feature set[END_REF].

The recording sites were close to each other (~1 km) and were all localized in the same subalpine spruce forest habitat. This spatial design was intentionally selected so that each recording site could be considered as a replicate sample. Indeed, our analysis could not reveal differences between the sites, confirming the replication role of each recorder.

However, we cannot rule out that further analyses over a longer observation time, that is over several years and dealing with more specific sounds, as birdsongs, mammal vocalizations or insect sounds, might reveal inter-site differences.

Using an automatic acoustic analysis, it was therefore possible to provide a first glimpse of the Risoux Forest soundscape structure over a complete year. This first analysis constitutes a baseline for future soundscape monitoring. The model generated by the ANN is efficient enough to analyze forthcoming datasets and display a reliable overview of the soundscape dynamics over years. Referring to this first year, it will be possible to estimate the relative importance of the biophony, geophony and anthropophony variation in the long term and to assess any possible shifts in diel and season patterns that could be linked to climate change. In addition, biophony automatic identification opens up the possibility to select recordings with biotic sounds and undertake more specific analyses about acoustic biodiversity. These analyses include the study of the composition and dynamics of acoustic communities in the long term, with a focus on the decline of specialist species in mountain environments [START_REF] Lehikoinen | Common montane birds are declining in northern Europe[END_REF][START_REF] Scridel | A review and metaanalysis of the effects of climate change on Holarctic mountain and upland bird populations[END_REF], an increase of generalist species, shifts in migration and nesting dates inducing new biophony phenologies. The development of sound classification tools should allow for more detailed long-term analysis about the evolution and impact of anthropic constraints over protected areas.

Figure 1 .

 1 Figure 1. Map of the study site. The recording site was localized in the Haut-Jura, in France at the border between France and Switzerland (inset). The four recorders were evenly distributed along the central zone of the Risoux forest which is crossed by roads, ski trails (left). The forest mainly consists of European spruce (Picea abies) (right).

Figure 2 .

 2 Figure 2. Sound types monitored. Examples of each sound type automatically recorded and detected. Short-time Fourier transform parameter: Hanning window of 512 samples with no overlap, 50 dB dynamic with a relative maximum defined for each sample. Obtained with seewave R package (Sueur 2018). For a sake of clarity the duration of the examples was limited to 10 s and with only one type of sound per example but the ANN could detect several sound types for each recording.

Figure 3 .

 3 Figure3. Sound type presence according to month. The percentage of each sound type was obtained by dividing the number of files containing the sound type by the total number of files. Note that the year starts with the beginning of the recordings, that is in August.

Figure 4 .Figure 5 .Figure 6 .

 456 Figure 4. Sound type presence according to hour. The percentage of each sound type was obtained by dividing the number of files containing the sound type by the total number of files. Note that the day starts at 1 am and finishes at midnight.

Acknowledgements

This work was supported by the Parc naturel du Haut-Jura who received funding from the Région Bourgogne-Franche-Comté and the Région Auvergne-Rhône-Alpes.