Maysa Yusef-Buey

Tzonka Mineva
email: tzonka.mineva@enscm.fr

Mathias Rapacioli
email: mathias.rapacioli@irsamc.ups-tlse.fr

Coupling Density Functional based Tight Binding with Class 1 Force Fields in a hybrid QM/MM scheme

Keywords:

QM/MM (Quantum Mechanics/Molecular Mechanics) hybrid methods have become very popular schemes to incorporate environmental effects in the calculation of molecular properties, when it is mandatory to have both a quantum description of electrons to compute these properties and an atomistic description of the environment. However, even Density Functional Theory/MM schemes may become timecosting when a large part of the system should be treated at the QM level or when plenty of single point energy calculations are intended to be done. We report a new implementation, within the deMonNano code, of a hybrid QM/MM scheme combining the density functional based tight binding with class 1 force fields. Two types of additive couplings can be chosen, namely the mechanical coupling, consisting of a Lennard-Jones potential and the electrostatic coupling, in which the MM part of the system is also polarizing the region described at the QM level. As first test-case application, the harmonic infrared spectra of simple molecules in the gas phase and in water clusters are computed and compared to those obtained at the DFT/MM level. Binding energies are also compared. Similar trends are obtained with the two levels of calculations and the main differences are discussed.

Introduction

The hybrid quantum mechanics/molecular mechanics (QM/MM) computational schemes owing their origin in the works of Warshel and Karplus [1] and Warshel and Levitt [2] incorporate simultaneously the efficiency of molecular mechanics algorithms and the accuracy of ab initio methods. The QM/MM approach is attractive with its simple idea to combine different levels of atomistic approaches and beneficial for large size systems requiring to explicitly consider the electrons only around a reactive center embedded in the electrostatic field of the environment, treated at molecular mechanics levels. The total energy of a QM/MM system (E QM/M M ) is obtained as sum of the energies of QM (E QM ) and MM (E M M ) subsystems, corrected by the energy term (E QM -M M ), accounting for the interaction between the two subsystems:

E QM/M M = E QM + E M M + E QM -M M (1)
The main question is how to compute the interaction energies between QM and MM regions. Among the different implementation schemes, two approaches emerged: the subtractive and the additive [START_REF] Groenhof | Introduction to QM/MM Simulations[END_REF]. In the former approach, the QM-MM interaction is obtained from the difference between the potential energy of the entire (QM + MM) system and the potential energy of the QM subsystem, both calculated at the MM level. An example of subtractive scheme is ONIOM [START_REF] Maseras | [END_REF]. The additive scheme is an embedding scheme considering the QM subsystem embedded by the MM subsystem and the QM-MM interactions are calculated explicitly at different levels of complexity by including mechanical, electrostatic and/or polarisation embedding techniques. Here, we have chosen to implement the additive QM/MM approach. The computations of energy gradients of MM atoms and their coupling to the energy gradients of QM atoms is straightforward. Geometrical optimisation and dynamics simulations based on both subtractive and additive DFT/MM schemes have been employed with success to study reactivity, spectroscopic features [5][6][7], to explore conformational space or to compute free energy profiles [8,9].

In DFT/MM schemes, most of the computational time is related to the calculation of the DFT energy, and its gradients when they are needed. When many single point energy calculations are intended to be done, like in global exploration searches or in long molecular dynamics (MD) simulations, or when the region treated at the QM level contains many atoms, it is appealing to replace the DFT QM scheme by less computationally expensive methods. In this domain, the Density Functional based Tight Binding (DFTB) is a well suited scheme [10][11][12][13][14]. DFTB can be seen as an approximated DFT, whose computational efficiency arises from the use of a minimal valence atomic basis (matrix diagonalisation scales as the cube of the number of basis functions) and the use of tabulated values (the so-called Slater Koster tables [15]) to avoid explicit calculations of integrals. With respect to force fields and empirical based tight binding schemes, a main advantage of DFTB is it's ab initio grounding as it is formally derived from DFT allowing for a greater transferability. As a consequence, it was shown to be a powerful method to simulate various processes, like chemical reactivity [16] or fragmentation [17], to compute different spectral properties, like electronic [18], harmonic [19,20] and anharmonic [21] IR and Raman [22] spectroscopy, to derive ionisation and electronic affinity energies [23], and to be valuable method to address various scientific fields ranging from material science [24] to biology [25].

Several implementations of DFTB/MM schemes within various MM packages such as AMBER [26,27], GROMACS [28], CHARMM [29] or the coupling of DFTB with the Universal Force Field (UFF) scheme in the deMonNano code [START_REF] Lima | Chapter 5 -Molecular Dynamics of Polypeptides and Their Inclusion Compounds with Cyclodextrin in Aqueous Solution Using DC-SCC-DFTB/UFF Approach[END_REF] have been previously reported. Let us also mention the coupling of DFTB with a polarisable force field to account for Argon matrix environment [START_REF] Iftner | [END_REF]. Other strategies have been used to consider the perturbation of a system described at the DFTB level by its environment, such as implicit models based on a continuous description of the medium [32,33] or models based on a coarse-graining strategy [34]. DFTB has also been included in QM/QM' schemes, acting as the lower level QM' method when combined with DFT [35].

In the present work, we report the coupling of DFTB with various force fields of class 1 (inter-molecular interactions are obtained by harmonic terms only and Lennard-Jones 6-12 potential), its implementation within the deMonNano code [START_REF] Heine | deMonNano[END_REF] and first applications to test cases. The present implementation is based on an additive coupling scheme, differing from the previous DFTB/UFF scheme within the same deMonNano code based on a substractive scheme [START_REF] Lima | Chapter 5 -Molecular Dynamics of Polypeptides and Their Inclusion Compounds with Cyclodextrin in Aqueous Solution Using DC-SCC-DFTB/UFF Approach[END_REF]. In addition, its implementation allows the use of any kind of parametrized class 1 force fields, as long as the user provides the appropriate parameters file. Furthermore, unlike some DFTB/MM coupling schemes combining different softwares [START_REF] Gillet | [END_REF]38], we report, a DFTB/MM hybrid implementation in a single code, as previously done in some other softwares [29,[START_REF] De M Seabra | The journal of physical chemistry[END_REF]40]. It avoids the use of any interface, reducing the loss of computational efficiency associated to communication processes. We compare harmonic vibrational spectra obtained at the DFT/MM and DFTB/MM levels for simple molecules, isolated or trapped in water clusters. In the next section, the basics of DFTB and MM methods are briefly reviewed before presenting the coupling strategy employed. Details about DFT/MM calculations are also given as well as the strategy used to compute energetic properties and harmonic spectra at the DFTB/MM level and to provide cross-comparisons with DFT/MM results. The laters are shown and discussed in section 5.

Incode additive DFTB/MM implementation

Density Functional based Tight Binding

An extended description of the DFTB method can be found in several review papers [13,14,41,42]. This theory was build in the 90's [10,11] following the pioneering works of Foulkes and Haydock [43]. In this scheme, the molecular orbitals (MOs) {φ i } are expanded on a non-orthogonal minimal atomic valence basis {ϕ µ } : φ i = µ c iµ ϕ µ , where c iµ are the MO's decomposition coefficients. The energy expression relies on a Taylor expansion of the effective Kohn-Sham potential around a reference density (ρ = ρ 0 + δρ) which can be performed up to the first, second or third order, leading to the various levels of DFTB, known as zeroth-order DFTB [10,11], self-consistent charges (SCC)-DFTB [44] and third order DFTB (also called DFTB3) [45,46]. In the widely used SCC-DFTB scheme, the energy is expressed as :

E = a,b E rep (R ab ) + i,µ,ν n i c iµ c iν H 0 µν + 1 2 a,b γ ab q a q b (2)
In this expression, E rep (R ab ) is a two-body repulsive contribution depending on the distance between atoms a and b H 0 µν is a matrix element of the Kohn-Sham operator taken at the reference density ρ 0 , n i is the occupation of the i th MO, q a is the atomic charge of atom a computed from the MO's coefficients within the Mulliken scheme and γ ab is a diatomic matrix accounting for the long range coulomb interaction between fluctuating charges plus some short range contributions of exchange correlation energy, the onsite term being equal to the atomic chemical hardness. The neglect of 3-centers integrals allows to parameterise the matrix elements of H 0 from atomic and diatomic DFT calculations and to write the second order contribution as a function of atomic charges. The electronic problem is then solved by applying the variational principle to this energy expression leading to the following secular equation :

(H 0 + H 1 (q))C i = i SC i ( 3
)
Where S is the atomic orbital overlap matrix and C i the vector containing the coefficients of the i th MO.

H 1 is the contribution to the operator arising from the second order energy contribution (last term in Eq. 2) and can be written as :

H 1 µν (q) = 1 2 S µν c q c (γ ac + γ bc ) (4)
with µ ∈ a and ν ∈ b. In the zeroth-order DFTB, the second order term is neglected and the operator is restricted to H 0 . In the DFTB3 scheme [45,46], including third order terms in the DFT Taylor expansion allows to consider the fluctuation of the Hubbard parameters with the electronic density expressed again as a function of the atomic charges q. This results in adding a new contribution to the DFTB energy expression (Eq. 2) and its corresponding perturbation H 2 (q) to the DFTB operator (Eq. 3). In the following, E DF T B and H DF T B can refer to either zeroth-order DFTB, SCC-DFTB or third order DFTB.

In practice, the results presented in section 5 were obtained from the third order DFTB combined with a correction for the γ matrix as recommended in reference [46]. We used the 3ob parameter set [47] downloaded from the www.dftb.org website.

Molecular Mechanics (MM)

In this section we briefly recall the fundamentals of the MM approach. The MM potential energy function is computed from the classical pairwise bonded and nonbonded inter-atomic interactions [START_REF] Monticelli | Force Fields for Classical Molecular Dynamics[END_REF]. The general ex-pression is given by:

E M M = N bond i V N bond i + N angles j V j + Ntors l V Ntors l (5) + Nimp l V N Imp l + N M M c N M M d>c V Coul cd + N M M c N M M d>c V LJ cd
The first four terms correspond to the bond stretching and angle bending contributions as well as to the dihedral and improper torsion (also known as local or shortrange) contributions. In class 1 FFs the bond stretching and angle bending terms are computed from harmonic functions [START_REF] González | [END_REF]:

V bond = 1 2 k b (r -r 0 ) 2 V angles = 1 2 k a (θ -θ 0 ) 2 (6)
where r -r 0 , θ -θ 0 are the deviations from the equilibrium bond distances and valence angles, and k b ,k a are the bond and angle force constants, respectively [50]. The torsion (dihedral and improper) potentials are modelled from periodic cosines functions :

V tors = V n 2 [1 + cos(nφ -δ)] (7)
where φ is the dihedral or torsional angle, n is an integer, non-negative constant, indicating the periodicity, δ is the phase and V n is the height of the rotational potential barrier [START_REF] González | [END_REF]. The last two terms in Eq.5 correspond to the Coulomb and the Van der Waals (vdW) interactions (non-bonding interactions). VdW are modelled by 12-6 Lennard-Jones potential as a combination of Pauli short-range repulsion and dispersion long-range attraction, given by the respective C cd 12 and C cd 6 parameters, and R cd the distance between atoms c and d. [START_REF] Groenhof | Introduction to QM/MM Simulations[END_REF][START_REF] González | [END_REF]:

V LJ cd = C cd 12 R cd 12 - C cd 6 R cd 6 (8)
An advantage of MM approach is that the bonding and non-bonding interactions are computed using different parameter sets, obtained from ab intio or semiempirical QM calculations, as well as from fitting experimental data. [START_REF] González | [END_REF]. Because we have coupled DFTB with class 1 FFs we presented above the energy potential terms only for this class of force fields. Some examples of class 1 FFs are OPLS, AMBER, CHARMM, GROMOS, etc [START_REF] Groenhof | Introduction to QM/MM Simulations[END_REF].

The general expression for the MM potential can be tuned, adding new terms in order to increase the level of accuracy of the molecular interactions. This is for example, the case of Urey-Bradly interactions [START_REF] González | [END_REF] V

U B = angles 1 2 k U B (s -s 0 ) 2 , ( 9
)
where s is the distance between two non-bonded atoms (1, 3) and k U B is the force constant parameter. As a final appreciation, in MM approaches, the electronic degrees of freedom are not considered explicitly. As a consequence, the impact of the E M M energy contribution on the secular equation can be, thus, ignored and its contribution to the total energy can be derived straightforwardly. [51] 2.3 DFTB/MM couplings In this section our implementations of DFTB/MM mechanical and electrostatic embedding schemes are presented.

Mechanical coupling

The simplest mechanical embedding comprises only the non-bonding interactions between QM and MM regions, thus the QM subsystem is not polarized by the MM subsystems and vice versa [START_REF] Groenhof | Introduction to QM/MM Simulations[END_REF]. The E QM -M M term in the DFTB/MM coupling, referred throughout this work as E QM -M M/mech is :

E QM -M M/mech = a∈QM c∈M M C ac 12 R ac 12 - C ac 6 R ac 6 (10)

Electrostatic coupling

The Coulomb interaction between the MM point charges and the DFTB density should in principle involve Coulomb integrals over the DFTB density described in the atomic basis. However, following the DFTB approximation used to compute the coulomb term in the second order energy correction, the DFTB density is approximated for this term by atomic charges using the Mulliken scheme formalism to replace the electronic density by atomic monopoles. It leads to the following expression for the QM/MM Coulomb potential energy :

E QM -M M/elec = a∈QM c∈M M
q a q c R ac (11) arising from the interaction between MM atomic charges (q c ) and DFTB atomic charges (q a ), derived from the electronic density in the DFTB region, is added to the Lennard-Jones potential (Sec. 2.3.1). As the atomic charges from the QM region are optimised self-consistently, the DFTB operator used to solve the secular Eq. 3 contains an additional term leading to the following expression for the matrix elements in the atomic basis :

H µν = H DF T B µν + H QM M M/elec µν with H QM M M/elec µ∈a;ν∈b = 1 2 S µν c∈M M q c ( 1 R ac + 1 R bc ) (12)
The analytical gradients have been derived and implemented for both mechanical and electrostatic embedding schemes. Benchmarks are reported in the supplementary material, comparing the calculation of the energy gradients from the analytical implementation and from numerical finite differences as well as the conservation of the total energy during MD micro-canonical dynamics. We note that the present implementation benefited in part from the former SCC-DFTB/UFF implementation within the deMonNano code [START_REF] Lima | Chapter 5 -Molecular Dynamics of Polypeptides and Their Inclusion Compounds with Cyclodextrin in Aqueous Solution Using DC-SCC-DFTB/UFF Approach[END_REF] and from the MM implementation within the deMon2k code [52] Actual version of DFTB QM/MM hybrid method implementation is already available in the deMonNano website1 [START_REF] Heine | deMonNano[END_REF] 3 DFT computational details

The DFT calculations were carried out with the de-Mon2k code [53] including the incode DFT/MM [52].

The PBE exchange correlation functional [54] was used in combination with the DZVP basis set [55]. Automatically generated auxiliary basis sets were employed [56]. An adaptive grid for numerical quadrature was used.

The SCF and gradient convergence tolerances were set to 10 -7 and 10 -5 a.u, respectively.

Computational Strategy

To validate the implemented scheme, we have compared DFT/MM and DFTB3/MM calculations, with the electrostatic coupling scheme. OPLS-AA (TIP3P for water molecules) has been used for the MM potential in both cases. The test systems consist of small molecules, namely ammonia, glycine, pentane and water, isolated (gas phase) or trapped in a water cluster ( between 205-215 water molecules). The goal here is to evaluate the effects of a water cluster environment, treated at MM level, on the infrared (IR) spectra of the solute molecules using either DFT/MM or DFTB3/MM schemes and to compute the corresponding binding energies.

A complete study of such systems would require global explorations of the complex potential energy surfaces, which is out of the scope of the present work with the focus on the in-code coupling between DFTB3 and MM schemes. Therefore, the procedure to generate the geometrical configurations prevents the use of the present results beyond methods comparison. In practice, our benchmark approach has two steps. First, the selected small molecules were optimized at both the DFT and DFTB3 levels in gas-phase. In the second step, each molecule was introduced in a water cluster. A short molecular dynamics (50 ps) was performed at the DFTB3/MM level in the canonical ensemble (temperature of 60 K maintained by a chain of 5 Nose Hoover thermostats [57][58][59] with energy exchange frequency of 800 cm -1 ). The final geometry was optimised at the DFTB3/MM level. Finally, these geometries were further optimised at the DFT/MM level. All the cluster's water molecules were treated at the MM level. For these optimised geometries in gas-phase and in water clusters, the IR spectra were computed in the harmonic approximation, involving the diagonalisation of the mass weighted Hessian matrix. In both, DFT/MM and DFTB3/MM implementations, only the dipole moments arising from the solute molecule, therefore the QM dipole moments, were considered in the calculations of the IR intensities. The relevant information is therefore the solvation effects on the IR spectra of the solute molecules and not the spectrum of the full cluster.

5 Applications to test cases 5.1 Ammonia Figure 1-top compares the ammonia IR spectra computed for the gas phase ammonia at the DFT and DFTB3 levels. At both the DFT and DFTB3 levels, the most intense band corresponds to N-H wagging vibrational mode and the frequencies are in good agreement (1022 cm -1 for DFTB3 vs 1030 cm -1 for DFT). A second band, by far less intense corresponds to H-N-H scissoring modes and the even less intense bands (almost no intensity at the DFT level) located around 3400 cm -1 result from N-H stretching modes. When NH 3 is inside the cluster, the main band (wagging mode) computed at the DFT/MM level is shifted towards higher frequencies (figure 1-bottom). This trend is captured by the DFTB3/MM spectra (figure 1-top), although the shift is smaller (93 cm -1 vs 187 cm -1 at the DFT/MM level). The intensity of this band is increased at both the DFT/MM and DFTB3/MM levels with respect to the isolated molecule spectrum. Very few changes are observed for the scissoring bands whatever the potential used. The intensities of the stretching modes are increased with either DFT/MM or DFTB3/MM. The last effects of the environment is the appearance of bands below 500 cm -1 visible in the DFT/MM spectrum and present, although much less intense, in the DFTB3/MM spectrum. These frequencies correspond to soft intermolecular modes. In some of them, the ammonia presents rotation and translation as a rigid body, whereas some others only correspond to intermolecular vibrations of water molecules without the involvement of displacements of ammonia atoms. In the latter modes, only the ammonia electronic density is perturbed. A little contribution from water environment is also observed at both the DFTB3/MM and DFT/MM levels around 3855 cm -1 and corresponds to combinations of the water stretching modes, inducing ammonia dipole fluctuations. Similarly, at both levels of theory, the small feature around 2100 cm -1 is associated to combinations of bending modes of many water molecules.

Neutral Glycine

The results obtained for the neutral glycine molecule are presented in Figure 2. Focusing first on the isolated molecule (top), it can be seen that the DFT spectrum presents an intense band at 3242.5 cm -1 assigned to the O-H stretching mode. This band is also among the most active ones at the DFTB3 level, although presenting a smaller intensity as well as a blueshift (3523.2 cm -1 ) with respect to the DFT results. Apart from this band, the DFT neutral glycine spectrum is dominated by two modes corresponding mostly to the displacement of a carbon atom from the carboxyle group at 1795 cm -1 and to a C-O-H bending mode at 1427 cm -1 . At the DFTB3 level, these two modes are the most active ones located at 1702 and 1252 cm -1 , respectively. Finally, several other modes are observed at lower frequencies in both spectra, the main difference being related to the modes around 900 cm -1 whose intensities are underestimated at the DFTB3 level in comparison to their DFT intensities.

When glycine is inside the water cluster, both methods lead to the same trends for the O-H streching mode, namely an increase of the intensity and a shift toward lower frequencies, the latter being stronger at the DFT/MM (395 cm -1 ) level than at the DFTB3/MM level (162 cm -1 ). The band at 1795 cm -1 (DFT value) is also redshifted in the cluster at both DFT/MM and DFT-gas DFT/MM-cluster DFT-gas Fig. 1 Harmonic IR spectra for ammonia. Top : gas phase molecule at the DFT and DFTB3 levels and solvated ammonia at DFTB3/MM levels; Bottom : gas phase and solvated ammonia at the DFT and DFT/MM levels. Harmonic stick spectra were convoluted with Lorentzian function (FWMH=20 cm -1 ) to produce the line spectra. For further analysis find in supplementary material comparison between DFTB3/MM and DFT/MM spectra for ammonia system DFTB3/MM levels, whereas a slight blue shift of the C-O-H bending mode is only observed at the DFTB3/MM level. Let us finally mention that, as reported in the case of ammonia, modes associated to either intermolecular vibrations, combinations of water bending or stretching modes are observed, in the low frequency domain and in the 2100 and 3850 cm -1 regions, respectively.

Pentane

The harmonic IR spectra obtained for the pentane molecule are reported on Figure 3. For the isolated DFT-gas DFT/MM-cluster DFT-gas Fig. 2 Harmonic IR spectra for neutral glycine. Top : gas phase molecule at the DFT and DFTB3 levels and solvated glycine at DFTB3/MM levels; Bottom : gas phase and solvated glycine at the DFT and DFT/MM levels.. Harmonic stick spectra were convoluted with Lorentzian function (FWMH=20 cm -1 ) to produce the line spectra.

molecule, the two levels of theory agree on the main spectral features, namely the presence of two distinct active regions : the [2900-3000] cm -1 region corresponding to various combinations of C-H streching modes and the [1300-1500] cm -1 region, corresponding to various combinations of C-H bending modes (essentially nodding and shearing modes). At lower frequencies, small contributions arise from modes associated to global deformations of the carbonaceous skeleton as well as some bending modes (torsion and swinging modes). From Figure 3, both top and bottom pannels, it can be seen that, in the opposite to what was observed in the case of ammonia and glycine, the main features of the har- DFT-gas DFT/MM-cluster DFT-gas Fig. 3 Harmonic IR spectra for pentane. Top : gas phase molecule at the DFT and DFTB3 levels and solvated pentane at DFTB3/MM levels; Bottom : gas phase and solvated pentane at the DFT and DFT/MM levels.. Harmonic stick spectra were convoluted with Lorentzian function (FWMH=20 cm -1 ) to produce the line spectra. monic spectrum do not present significant changes due to the cluster environment, whatever the level of theory. This is well in line with the known hydrophobic nature of pentane [60], which does not form strong Hbonds with the surrounding water molecules. The major observed changes are the appearance of frequencies in the 2100 and 3850 cm -1 regions, assigned to water H-O-H bending modes and O-H stretching vibrations, as well as in the low frequencies region, associated to intermolecular modes.

Water

Our last example is the analysis of water spectra which are reported on Figure 4. Focusing first on the isolated molecule, the DFT spectrum indicates a bending mode at 1635 cm -1 and two stretching modes at 3721 and 3854 cm -1 for the symetric and antisymetric modes, respectively. These bands are also observed at the DFTB3 level at 1461, 3604 and 3877 cm -1 , respectively. The intensities are in agreement with the exception of the symmetric stretching mode which is active in the DFTB3 spectrum but presents very low intensity in the DFT spectrum. Inside the water cluster, the investigated water molecule presents similar modifications when computed either at the DFT/MM or DFTB3/MM levels, namely an increase of the absorption intensities, in particular for the two stretching modes, a red-shift of the stretching modes and a blueshift of the bending mode. One can also observe, in both spectra, broad features below 1100 cm -1 which originate from global intermolecular modes. The main difference is that, in the DFT/MM spectrum, one of these modes, located at 1085 cm -1 presents a strong intensity, whereas in the DFTB3/MM spectrum no intermolecular mode is observed with such a strong intensity.

We finally discuss the influence of the size of the region treated at the QM level. To do so, we have optimized and computed the harmonic spectrum for the water cluster treating all the water molecules (i.e. those previously described either at the QM or at the MM level) at the DFTB3 level. To be consistent with the previous calculations, only the dipole moment of the molecule previously described at the QM level is included in the intensity calculation. This spectrum is compared with the one of the isolated water molecule on Figure 5. As in the calculation performed at the DFTB3/MM level, the O-H stretching modes are redshifted and their intensities are increased. Let us note that the intensity is now spread over many almost degenerated modes. We also recover the appearance of a broad pattern at low frequencies. The main difference relies in the bending modes, which was blue shifted in the DFTB3/MM spectrum, but is almost not impacted by the environment, apart from the fact that its intensity is now spread over many modes with very similar frequencies. This can clearly be related to the fact that the frequency of the water bending mode is 1471 cm -1 when computed at the DFTB3 level and 2028 cm -1 when computed at the MM level. As a consequence, at the DFTB3/MM level the coupling of the bending mode of the QM molecule with the bending modes of the MM molecules at higher frequency are responsible for the observed blue shift. To conclude on this point, we would like to stress that the size of the subsystem treated at the QM level should be large enough to limit such artifacts originating from the difference between the water molecules frequencies in the two levels of calculations and that force fields of higher accuracy (polarisable force fields or TIP4P, or SPC, or TIP3F) could provide better descriptions of water frequencies.

Binding energies

We finally compare the binding energies between the two subsystems, i.e. between the studied molecule (NH 3 , C 2 H 5 NO 2 , C 5 H 12 or H 2 O), treated at the QM level, and the surrounding water cluster, treated at the MM level). This binding energy is computed from the energies of the full system minus the sum of energies of the two subsystems :

E bind = E(QM + M M ) -E(QM ) -E(M M ) (13)
The geometries were not relaxed to compute the energies E(QM) and E(MM), i.e. the coordinates of the QM molecule for E(QM) (resp. of the water cluster for E(MM)) having the same values as in the full system. , the atomic charges of the QM and MM regions induce couplings between the orbitals of different atoms but not between orbitals of the same atom, due to the fact that the overlap matrices elements S µν are zero if µ and ν correspond to two different orbitals of the same atom. As a consequence, an MM region would not induce polarisation of an isolated QM atom. This neglect of atomic polarisation energy could be, at least partially, at the origin of an underestimation of the binding energy at the DFTB3/MM level. Let us note in particular that the atomic polarisabilities of some elements [START_REF] Miller | Atomic and Molecular Polarizabilities-A Review of Recent Advances[END_REF], like nitrogen (1.100 Å3 ) or oxygen (0.802 Å3 ), are significant with respect to molecular polarisabilities [START_REF] Olney | [END_REF] of ammonia (2.103 Å3 ) or water (1.501 Å3 ) .

Conclusion

We have reported a new implementation for an hybrid QM/MM scheme where the QM part of the system is treated at the DFTB level and the MM part of the system using class 1 force fields approaches. The DFTB level can vary from zeroth-order DFTB, SCC-DFTB to third order DFTB. Only the class 1 force fields AMBER and OPLS-AA have been used up to now but the implementation allows to use any class 1 Force Filed, as long as the user provides the corresponding parameters. The additive scheme has been implemented and two levels of coupling can be chosen by the user, namely the mechanical and the electrostatic coupling. In both cases, the interactions between the MM and the QM subparts of the system are computed from a Lenard-Jones potential. In the electrostatic coupling, an additional term describes the coulomb interactions between the atomic charges of the two regions. The atomic charges from the MM domain remain constant but those of the QM region are optimised self-consistently within the DFTB variational scheme, including the polarisation of the QM region by the charges of the MM region. The energy gradients, (mandatory to perform molecular dynamics ) of compute vibrational spectra , have also been implemented for this hybrid scheme. Benchmarks have been performed by comparing the latter computed either from the analytical formula or from finite differences and by checking that the total energy was conserved during micro-canonical simulations. For a set of selected small molecules in water clusters, harmonic spectra have been computed at the DFTB3/MM and DFT/MM levels, making use of the same force field. The geometrical configurations were arbitrarily generated, , they do not exemplify any real system. Any analysis of the properties and characteristic of this geometries can not be extracted. Nevertheless, these test cases showed that the main environment effects computed at the DFT/MM level are reproduced at the DFTB3/MM level.

In the future, the combination of the present implementation with schemes similar to the one developed by Iftner et al. [START_REF] Iftner | [END_REF] to couple polarisable Argon atoms description at the MM level with the DFTB scheme would allow us to mix DFTB with polarisable force fields similarily to what was done at the DFT level [52], [63] ents, namely a Table S1 comparing energy gradients computed for an arbitrarily chosen geometry of a cluster consisting of an ammonia molecule (QM level) and a water molecule (MM level) computed either from the analytical implementation and through finite differences and a Figure S1 showing the conservation of the total energy during an MD simulation in the microcanonical ensemble for an ammonia molecule embedded in a water cluster (215 molecules). Furthermore, a comparison of IR Harmonic spectra of soluted ammonia in water cluster at DFTB3/MM and DFT/MM levels of theory are included in the supplemental material (Figure S2). A brief analysis from both plots is extracted, showing that the effect of water environment on the system treated by DFT/MM method is reproduced by the DFTB3/MM level of theory.

Fig. 4

 4 Fig.4Harmonic IR spectra for water. Top : gas phase molecule at the DFT and DFTB3 levels and solvated water molecule at DFTB3/MM levels; Bottom : gas phase and solvated water molecule at the DFT and DFT/MM levels. Harmonic stick spectra were convoluted with Lorentzian function (FWMH=20 cm -1 ) to produce the line spectra.

Fig. 5

 5 Fig.5Harmonic spectra for water molecule in the gas phase and solvated in a water cluster computed at the DFTB3 level.

Table 1

 1

	reports the binding energies computed
	for the geometries optimized either at the DFTB3/MM
	or DFT/MM levels. Small differences (less than 0.13
	eV) can be observed between the binding energies ob-
	tained at the DFT/MM level for the geometries result-
	ing from local optimisation at the DFTB3/MM (3 rd col-
	umn) or DFT/MM (2nd column) levels. This indicates
	small structural relaxation when changing the compu-
	tational scheme. The binding energies computed at the
	DFT/MM and DFTB3/MM levels for the same geome-
	tries (optimised at the DFTB3/MM level, columns 3
	and 4) are in very good agreement for the glycine and
	pentane molecule, whereas for ammonia and water, the
	values obtained at the DFT/MM level (-0.649 and -
	1.380 eV, respectively) are larger than those obtained
	at the DFTB3/MM level (-0.394 and -1.001 eV, respec-
	tively). In order to track the origin of this difference, we
	report the contribution of the Lennard-Jones potential
	to the coupling energy, this value being the same in

Table 1

 1 Binding energies computed at the DFT/MM or DFTB3/MM levels from Eq. 13 for the four investigated systems at the geometries optimized either at the DFT/MM or DFTB3/MM levels of theory. Lennard-Jones contributions to the binding energies are also reported.both the DFT/MM and DFTB/MM schemes. It appears that, appart from the pentane case where it almost equals the binding energy, this contribution is repulsive. It means that most of the stabilising interaction arises from the interaction between the charges of the MM subsystem and electronic density of the QM subsystem. The difference could originate in a lack of polarisation at the DFTB/MM level. In the DFTB and DFTB/MM operator expressions (Equations 4 and 12)

The present version of the code corresponding to this work can be downloaded from deMonNano website (10.12.2021) and corresponds to git hash commit f3707842dd69b9629954a55a28693189d572f425

Acknowledgements The authors acknowledge the ANR agency (project ANR-19-CE29-0011-02 RUBI) for financial support and the supercomputing facilities of CALMIP for generous allocation of computer resources (project P0059 and P18009).

Supplemental Material

We provide, in Supplemental Material, benchmarks for the implementation of DFTB3/MM energy and gradi-


David Kaniewski
email: david.kaniewski@univ-tlse3.fr

Nick Marriner

Rachid Cheddadi

Christophe Morhange

Miguel Ángel

Cau Ontiveros

Joan J Fornós

Matthieu Giaime

Valérie Trichon

Thierry Otto

Frédéric Luce

Elise Van Campo

Affiliations

Recent anthropogenic climate change exceeds the rate and magnitude of natural Holocene variability on the Balearic Islands

Keywords: Ecosystem dynamics, Pollen-based climate reconstruction, climate change, Industrial-era warming, Aridification, island, Holocene, Mallorca, Mediterranean

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

The emergence of the industrial economy has triggered climate anomalies since the mid-late 19 th century [START_REF] Abram | Early onset of industrial-era warming across the oceans and continents[END_REF][START_REF] Gleckler | Industrial-era global ocean heat uptake doubles in recent decades[END_REF]. Although the resulting environmental stresses have affected many regions of the world, the International Panel on Climate Change (IPCC, 2001;IPCC, 2014;IPCC 2018) deem the Mediterranean to be "particularly sensitive" to the effects of global warming [START_REF] Fischer | Consistent geographical patterns of changes in high-impact European heatwaves[END_REF][START_REF] Homar | Recent trends in temperature and precipitation over the Balearic Islands (Spain)[END_REF][START_REF] López Mayol | PREGRIDBAL 1.0: towards a high-resolution rainfall atlas for the Balearic Islands (1950-2009)[END_REF][START_REF] Cramer | Climate change and interconnected risks to sustainable development in the Mediterranean[END_REF]. For the period 2081-2100, the IPCC's climate models project a 2.6-4.8°C rise in average surface temperatures (relative to 1986-2005), accompanied by a pronounced shift in rainfall regimes with anomalies of -4% to -27% in the Mediterranean (from Representative Concentration Pathways 4.5 to 8.5; IPCC, 2014). Such changes will generate exceptional stress on already limited water resources [START_REF] Guiot | Climate change: The 2015 Paris Agreement thresholds and Mediterranean basin ecosystems[END_REF], particularly in regions where anthropogenic demands already exceed supply. The projected warming of the Mediterranean is clearly attested on the Balearic Islands, which have witnessed an intensification of temperature anomalies during the last 40 years (0.6°C per decade) and a clear shortening of the cold season [START_REF] Cuxart | Observed trends in frost and hours of cold in Majorca[END_REF][START_REF] Homar | Recent trends in temperature and precipitation over the Balearic Islands (Spain)[END_REF][START_REF] López Mayol | PREGRIDBAL 1.0: towards a high-resolution rainfall atlas for the Balearic Islands (1950-2009)[END_REF]. According to the IPCC (2001), a rise in temperature (1.5 to 3.6 °C) and a drop in rainfall amounts (10% to 20%) are forecasted for 2080-2100 on the Balearics, with an increased probability of drought [START_REF] Kerr | Confronting the bogeyman of the climate system[END_REF] and extreme weather events [START_REF] Candela | Impact assessment of combined climate and management scenarios on groundwater resources and associated wetland (Majorca, Spain)[END_REF]. Even if the ocean's thermal inertia mitigates the temperature of the Mediterranean sea surface, delaying the accumulation of heat in coastal zones, notable increases in sea surface temperatures have already been recorded in the Balearic Sea (+1.1 °C during the last 40 years; [START_REF] Marbà | Mediterranean warming triggers seagrass (Posidonia oceanica) shoot mortality[END_REF][START_REF] Jordà | Mediterranean seagrass vulnerable to regional climate warming[END_REF][START_REF] Adloff | Mediterranean Sea response to climate change in an ensemble of twenty first century scenarios[END_REF].

Through the prism of Mallorca, this study aims to answer two main questions. First, to what extent does the industrial-era warming and aridification of Mediterranean islands deviate markedly from long-term Holocene trends? Second, how can an 8000-year palaeoecological archive from the Albufera Natural Park serve as a basis to quantify the current climate anomalies (1975( -2018 CE) CE)? As climate models project widespread anomalies across the Mediterranean Basin by the end of the 21 st century [START_REF] Cook | Old World megadroughts and pluvials during the Common Era[END_REF], and because warming and aridification can significantly impact the environment, economy, society [START_REF] Greve | Global assessment of trends in wetting and drying over land[END_REF] and human health [START_REF] Gale | Assessing the impact of climate change on vector-borne viruses in the EU through the elicitation of expert opinion[END_REF], quantifying the present-day climate trend in relation to the past is critical in improving climate-change projections. It is also important for developing robust hazard assessments and adaptation strategies.

Albufera Natural Park

The swamp/wetland of Albufera (elevation 1 m a.s.l.; surface of ~25 km 2 ; Fig. 1) is located in the eastern part of a Plio-Quaternary depression, between the northern and southern ridges of the island (dated from the Mesozoic to the Upper Miocene). The swamp is relict of a palaeolagoon formed following the Holocene marine transgression [START_REF] Fornós | Evolució neógena de la zona de s'Albulera d'Alcúdia (Mallorca, lIIes Balears): 1. Descripció de la serie estratigrafica (sondatge 1) i lacies sedimentaries[END_REF][START_REF] Goy | Evolución de la áreas margino-litorales de la costa de Mallorca (I. Baleares) durante el Último y Presente Interglaciar nivel del mar Holoceno y clima[END_REF][START_REF] Giaime | In search of Pollentia's southern harbour: Geoarchaeological evidence from the Bay of Alcúdia (Mallorca, Spain)[END_REF]. This swamp area is fed by the Sant Miquel and Muro streams and by groundwater inputs. We selected this protected natural reserve to limit any potential anthropogenic influences on the environmental-climate reconstructions. Today, the swamp is mostly populated by Phragmites australis (Cav. Trin. ex Steud.), Cladium mariscus (L. Pohl) and Typha latifolia (L.). Alisma plantago-aquatica (L.) and Sparganium erectum (L.) are also common around the edges, between the water and the marsh. In the water column, Potamogeton pectinatus (L.), Potamogeton crispus (L.), Riccia (L.), Myriophyllum spicatum (L.), Ceratophyllum submersum (L.) and Zannichellia pedunculata (Rchb. in Mössler) occur.

Around the park and closer to the dunes, the dominant plants are Pinus halepensis (Mill.), Pistacia lentiscus (L.), Rosmarinus officinalis (L.), Phillyrea angustifolia (L.), Cistus salviifolius (L.), Smilax aspera (L.), Asparagus acutifolius (L.), Halimium halimifolium (L.) and Brachypodium retusum (L.). On the seaward flank of the swamp, Juniperus oxycedrus subsp. macrocarpa (Sm. Ball) is the dominant tree species, with a cortege of Cakile maritima (Scop.), Eryngium maritimum (L.), Pancratium maritimum (L.), Matthiola sinuata (L.) and Helichrysum stoechas (L. Moench). The nomenclature derives from the Flora Iberica (http://www.floraiberica.es/).

Materials and methods

To answer the questions raised in the introduction, we used the instrumental climate dataset from Mallorca's airport weather station and produced a pollen-based climate reconstruction from a core drilled in the Albufera Natural Park (NE Mallorca). We used a 14 C chronology to calculate an age model for the core. We further tested the reliability of our climate reconstruction by estimating the role of anthropogenic impacts in modulating ecosystem dynamics. All data are available in tabulated form in the supplementary information.

Present-day climate data

The instrumental climate dataset was provided by Mallorca's airport weather station (Palma Aeropuerto, 39°33'39"N, 02°44'12"E, 3 m. a.s.l.). The dataset constitutes the Balearic Islands' longest temperature and precipitation records and covers the periods 1975-2018 (temperatures) and 1972-2018 (precipitation). The records from Palma Puerto (39°33'12"N, 02°37'31"E; 1 m a.s.l.; period 1978-2018) and Sa Pobla (39°44'57"N, 03°01'00"E; 27 m a.s.l.; period 2012-2018) were used as independent controls (for extreme values). Previous datasets were already available in the literature (e.g. [START_REF] Cuxart | Observed trends in frost and hours of cold in Majorca[END_REF][START_REF] Homar | Recent trends in temperature and precipitation over the Balearic Islands (Spain)[END_REF][START_REF] Jansà | Extension of summer climatic conditions into spring in the Western Mediterranean area[END_REF][START_REF] López Mayol | PREGRIDBAL 1.0: towards a high-resolution rainfall atlas for the Balearic Islands (1950-2009)[END_REF].

Core chronology and stratigraphy

We drilled an 8-m long core (39°47'44.43"N, 3°05'31.47"E, 1 m a.s.l.; termed POL12) at the heart of the Albufera Natural Park using a hydraulic rotary corer (Fig. 2). The sedimentary deposits mainly consist of clay with inclusions of silt and sand in the upper part. The chronology of the core is based on ten 14 C dates (Beta Analytic and Poznan Laboratory; see Table S1) calibrated to 1 and 2 sigmas (σ) using Calib Rev 7.0.4 (Fig. 2). The oldest level (790 cm depth) is dated to 7080 ± 40 14 C year BP, calibrated to 7905 year BP (7980-7830 calibrated 2σ BP, 7955-7865 calibrated 1σ BP). The age of each sample was calculated by interpolating between two adjacent 14 C dates. Our record thus covers ~8000 years of northeastern Mallorcan environmental and climate history.

Palaeoecological data

We prepared a total of 120 samples for pollen analysis using standard procedures for clay deposits [START_REF] Faegri | Textbook of pollen analysis[END_REF]. Pollen grains were counted under x400 and x1000 magnification using an Olympus Bh-2 BHS microscope. Pollen percentages are based on the terrestrial pollen sum, excluding local helophytes (Cyperaceae and Typha), macrophytes and spores of non-vascular cryptogams. Percentages of helophyte-macrophyte taxa were calculated using a sum that includes local freshwater plants and terrestrial vegetation. The mean pollen sum was 379 ± 53 pollen grains, with a minimum of 284 pollen grains. The median value was 374 pollen grains, with a 25th percentile of 336 pollen grains and a 75th percentile of 409 pollen grains. The mean number of taxa was 28 ± 3 taxa, with a minimum of 21 taxa. The median value was 27 taxa, with a 25th percentile of 26 taxa, and a 75th percentile of 29 taxa. The calculation of pollen concentrations follows the formula developed by [START_REF] Cour | Nouvelles techniques de détection des flux et des retombées polliniques : étude de la sédimentation des pollens et des spores à la surface du sol[END_REF], using 1 cm 3 of sediment. The mean concentration was 2728 ± 382 pollen grains per cm 3 with a median value of 2692 pollen grains per cm 3 , a 25th percentile of 2419 pollen grains per cm 3 , and a 75th percentile of 2948 pollen grains per cm 3 . Atlases [START_REF] Reille | Pollen et spores d'Europe et d'Afrique du nord[END_REF][START_REF] Reille | Pollen et spores d'Europe et d'Afrique du nord[END_REF][START_REF] Reille | Pollen et spores d'Europe et d'Afrique du nord[END_REF]), books (e.g. Faegri and Iversen, 1989;[START_REF] Moore | Pollen analysis[END_REF][START_REF] Punt | The northwest European pollen flora[END_REF], articles (e.g. [START_REF] Renault-Miskovsky | Observations de quelques pollens d'Oléacées au microscope électronique à balayage[END_REF][START_REF] Colombo | Pollen grain morphology supports the taxonomical discrimination of Mediterranean oaks (Quercus, Fagaceae)[END_REF][START_REF] Gortemaker | A method to identify pollen of some recent and fossil species of Fagus L. (Fagaceae)[END_REF], as well as local reference slides were used for pollen identification. A plate with pollen photos is provided in the supplementary information (see Fig. S1). The fire history of the coastal plain was reconstructed by counting charcoal particles (50-200 μm) on the pollen slides. The 50 μm size criterion was chosen to avoid confusion between microscopic charcoal fragments and opaque minerals, which are typically < 50 μm [START_REF] Parshall | Fire on the New England landscape: regional and temporal variation, cultural and environmental controls[END_REF][START_REF] Pederson | Medieval Warming, Little Ice Age, and European impact on the environment during the last millennium in the lower Hudson Valley[END_REF]. We followed the 200-300 items (sum of charcoal particles) defined by [START_REF] Finsinger | Minimum count sums for charcoal concentration estimates in pollen slides: accuracy and potential errors[END_REF]. A charcoal concentration curve (fragments per cm -3 ) was drawn on a linear age-scale (Marriner et al., 2019). The mean concentration was 77 ± 19 fragments including the standard error, with a maximum value of 1550 fragments and minimum of 0. The pollen data are detailed and presented as a pollen diagram (Fig. 3), with a focus on cereals, olive trees and fire activity (Fig. 4). All of these elements can be consulted in the supplementary information. We curated all the pollen slides at the Université Toulouse 3.

Statistical analyses

We applied a homogeneity test to fire activity, cereals and olive trees to detect potential shifts in their long-term dynamics (Fig. 4). The homogeneity test is based on a null hypothesis that a time series is homogenous. Possible alternative hypotheses such as change in distribution, changes in average or presence of trend can be detected and subdivided. Each change is highlighted and its average is indicated by mu. Before any modelling, a CABFAC analysis (Fig. 5) was used to test the potential weight of anthropogenic factors in mediating the ecological signal (e.g. [START_REF] Hammer | Paleontological data analysis[END_REF][START_REF] Kaniewski | Croatia's Mid-Late Holocene (5200-3200 BP) coastal vegetation shaped by human societies[END_REF]. The CABFAC factor analysis (Q-mode factor analysis) implements the classic method of factor analysis and environmental regression (CABFAC and REGRESS). We regressed selected environmental data on the CABFAC factors using the second-order (parabolic) method, with cross terms. The environmental regression model (RMA) reports the observed values against the values reconstructed from the factors [START_REF] Kaniewski | Croatia's Mid-Late Holocene (5200-3200 BP) coastal vegetation shaped by human societies[END_REF]. The selected environmental data for this study (wet/dry conditions, warm/cool conditions, anthropogenic activities, and fire) cover a wide range of natural and anthropogenic pressures. The consistency of each reconstructed model is indicated by the coefficient of determination R 2 . The R 2 values were ordered from highest to lowest, and used as a hierarchy of pressures on ecosystems.

Reconstructions of past climate variables

We did not consider local helophytes, macrophytes and spores of non-vascular cryptogams in the pollen-based climate reconstruction. The cultivated species (primary anthropogenic indicators: Poaceae cerealia) were also excluded from the pollen-based climate reconstruction in addition to the secondary anthropogenic indicators (e.g. Plantago, Rumex).

Pollen assemblages were used to reconstruct the palaeoclimate history of Mallorca. Our statistical approach presents two key caveats: (1) a fossil pollen assemblage may not have a modern analog sensu stricto; and (2) the climate range of a species/genus can shift in time [START_REF] Cheddadi | Thermal niches of three European tree species during the last millennia[END_REF]. Our new method is based on the assumption that the co-occurrence of plant taxa in an assemblage may only occur under a common climate range or an overlap of the climatic range of plant taxa comprising a given assemblage. Based on this assumption, one should not expect two plant species to occur in the same assemblage if their climate ranges do not intersect. We assigned the fossil pollen taxa to modern species, and subsequently calculated a weighted median value and its standard deviations of the climate range encompassed by all assigned modern species. Such an approach has been used to infer past climate variables using insects [START_REF] Elias | The mutual climatic range method of palaeoclimate reconstruction based on insect fossils: new applications and interhemispheric comparisons[END_REF], plant fossil remains [START_REF] Pross | Reconstructing paleotemperatures for the Early and Middle Pleistocene using the mutual climatic range method based on plant fossil[END_REF], mollusks [START_REF] Moine | Paleoclimatic reconstruction using mutual climatic range on terrestrial mollusks[END_REF], and ostracods [START_REF] Horne | A mutual temperature range method for Quaternary palaeoclimatic analysis using European nonmarine Ostracoda[END_REF] or to derive climate probability density functions [START_REF] Kühl | Probability density functions as botanicalclimatological transfer functions for climate reconstruction[END_REF] from fossil pollen data.

Our method [START_REF] Cheddadi | Microrefugia, climate change, and conservation of Cedrus atlantica in the Rif Mountains, Morocco[END_REF]) also performs a "leave-one-out calculation" to accommodate the temporal variability of pollen taxa in the fossil record. For each fossil sample, a taxon is removed and a weighted median (using the pollen percentages as a weight) is calculated using all of the remaining taxa. For each sample, we calculated the weighted median as many times as there were taxa in each fossil level. The final figure is the median value of all iterations. In addition, the standard deviations correspond to the median value of all the standard deviations.

This additional leave-one-out approach reduces the effects of either large variations in the fossil record or over-or under-representation of some pollen species. The code was written using the software R and exploits the available libraries (Akima and Gebhardt, 2016), RMySQL [START_REF] Ooms | RMySQL: database interface and[END_REF], and the statistical functions which are part of R. The modern plant database is based on georeferenced species from Flora Europaea [START_REF] Jalas | Committee for Mapping the Flora of Europe and Societas Biologica Fennica[END_REF][START_REF] Jalas | Committee for Mapping the Flora of Europe and Societas Biologica Fennica[END_REF][START_REF] Jalas | Committee for Mapping the Flora of Europe and Societas Biologica Fennica[END_REF] and from GBIF [START_REF] Heughebaert | Global Biodiversity Information Facility[END_REF]. We obtained present climate variables from the WORLDCLIM database [START_REF] Hijmans | Very high resolution interpolated climate surfaces for global land areas[END_REF].

Data availability

All data needed to evaluate the conclusions reported in the paper are present in the article and/or the supplementary information (excel files).

Results

The core

The bottom of the core comprises clay deposits (from 800 cm to 320 cm depth) consistent with a low-energy environment (Fig. 2). We identified a sand lens between 320 and 300 cm depth, which probably results from aeolian transport into the lagoon [START_REF] Giaime | In search of Pollentia's southern harbour: Geoarchaeological evidence from the Bay of Alcúdia (Mallorca, Spain)[END_REF]. The upper part (from 300 cm to 170 cm depth) consists of silt associated with freshwater inputs into the area. Some sand particles from the shallowing lagoon were mixed with the silt deposits during this period. The human-modified deposits correspond to the last 170 cm.

Pollen analysis

The 14 C dates (Fig. 2) frame the chronology of core POL12. The pollen data allowed us to reconstruct ecosystem dynamics during the past 8000 years at Albufera Natural Park (Fig. 3).

We compared and contrasted our data with earlier studies by [START_REF] Burjachs | Dinámica de la vegetación durante el Holoceno en la isla de Mallorca[END_REF][START_REF] Burjachs | Overview of environmental changes and human colonization in the Balearic Islands (Western Mediterranean) and their impacts on vegetation composition during the Holocene[END_REF] in the same area, and with cores drilled on Minorca and Ibiza [START_REF] Yll | Palynological evidence for climatic change and human activity during the Holocene on Minorca (Balearic Islands)[END_REF][START_REF] Pérez-Obiol | Similarities and dissimilarities, synchronisms and diachronisms in the Holocene vegetation history of the Balearic Islands and Sicily[END_REF][START_REF] Pérez-Obiol | Mid-Holocene vegetation and climatic history of the Iberian Peninsula[END_REF][START_REF] Carrión | Paleoflora y Paleovegetación de la Península Ibérica e Islas Baleares: Plioceno-Cuaternario[END_REF][START_REF] Burjachs | Overview of environmental changes and human colonization in the Balearic Islands (Western Mediterranean) and their impacts on vegetation composition during the Holocene[END_REF].

The date proposed for the introduction of agricultural practices and human-induced fires in the Albufera Natural Park is ~3600 BP (Fig. 4). According to the 2-sigma calibration, the full range is framed between 3650 and 3550 BP (3600 ± 50 BP). As stated by [START_REF] Burjachs | Overview of environmental changes and human colonization in the Balearic Islands (Western Mediterranean) and their impacts on vegetation composition during the Holocene[END_REF], cereals are relatively scarce in the area. Data from the core POL12 depict a similar trend, with a maximum value < 10%. In the core POL12, the earliest record of human impacts on the environment thus occurred after the shift from a mesophilous community to a maquis scrub (~5900 BP; Fig. 3).

The abrupt decline of Juniperus at ~5900 BP clearly illustrates this phenomenon (Fig. 3). The drop seems to result from environmental pressures because no strong anthropogenic footprints (fire activity or agriculture) were detected (Fig. 4). This decline is associated with greater aridity in the Mediterranean at ~5900 BP [START_REF] Jalut | Holocene climatic changes in the western Mediterranean, from south-east France to south-east Spain[END_REF][START_REF] Fletcher | Abrupt climate changes of the last deglaciation detected in a western Mediterranean forest record[END_REF]. As mentioned by [START_REF] Burjachs | Overview of environmental changes and human colonization in the Balearic Islands (Western Mediterranean) and their impacts on vegetation composition during the Holocene[END_REF], this climate event (along with others at 8.2 and 4.3 ka BP) directly influenced the composition and evolution of vegetation on the Balearic Islands. According to [START_REF] Pérez-Obiol | Similarities and dissimilarities, synchronisms and diachronisms in the Holocene vegetation history of the Balearic Islands and Sicily[END_REF], the important ecosystems changes ~5000-4000 BP seem to have been amplified, but not caused, by human action. Most south-western European pollen records show a general change from wetter to drier climatic conditions [START_REF] Jalut | Holocene climatic changes in the western Mediterranean, from south-east France to south-east Spain[END_REF][START_REF] Magny | Assessment of the impact of climate and anthropogenic factors on Holocene Mediterranean vegetation in Europe on the basis of palaeohydrological records[END_REF]. The core POL12 meshes perfectly with previously published data.

A homogeneity test was applied to three curves (fire activities, cereals and olive trees) to detect potential shifts in their long-term dynamics (averages are indicated by mu; Fig. 4). Results

showed that cereals and fire activity are well correlated, whereas olive trees are mediated by a different dynamic. To further establish whether olive trees were an integral part of early diversified agriculture, we calculated a cross-correlation (Pvalue = 0.05). We considered positive correlation coefficients, focusing on the Lag0 value (with +0.50 as the significant threshold).

The test, applied to Poaceae cerealia versus Olea europaea, showed no correlation between the two (Lag0 = 0.008, Pvalue < 0.001). A second test, based on fire activity versus Olea europaea, also demonstrated no correlation (Lag0 = -0.07, Pvalue < 0.001). According to [START_REF] Pérez-Obiol | Similarities and dissimilarities, synchronisms and diachronisms in the Holocene vegetation history of the Balearic Islands and Sicily[END_REF], from ~5000 to ~4000 BP, some taxa (Buxus and Corylus) were replaced by others such as Olea europaea, a process that reflects an ecological response to global climatic trends. A similar phenomenon was clearly recorded at Albufera (Fig. 3). At Cala'n Porter (Minorca), high occurrences of olive trees have been recorded since ~6500 BP. At Algendar (Minorca) and Prat de Vila (Ibiza), high values are recorded since ~5500 BP [START_REF] Burjachs | Overview of environmental changes and human colonization in the Balearic Islands (Western Mediterranean) and their impacts on vegetation composition during the Holocene[END_REF]). It appears that olive-tree dynamics are independent of human influence in the Balearics.

Nonetheless, some chronological discrepancies exist between the two studies in the Albufera Natural Park. These differences may result from the absence of dates between the 14 C ages 6270 ± 70 and 2440 ± 60 (respectively 7330-6985 and 2625-2350 calibrated 2σ BP) in the sequence published by [START_REF] Burjachs | Dinámica de la vegetación durante el Holoceno en la isla de Mallorca[END_REF][START_REF] Burjachs | Overview of environmental changes and human colonization in the Balearic Islands (Western Mediterranean) and their impacts on vegetation composition during the Holocene[END_REF]. Because the quality of an age model is dependent on the number of dated levels, more 14 C dates for the core POL12 have probably reduced the uncertainties for the older part of the sequence (Fig. 2).

Anthropogenic influence on the signal

Our pollen diagram from the Albufera Natural Park shows that, while human activities may have partially shaped the area since ~3600 BP (attested through an increase in cereal cultivation and fire activity), the CABFAC analyses indicate that the two main factors controlling the ecological signal are humidity (R 2 = 0.73, Pvalue < 0.001) and temperatures (R 2 = 0.65, Pvalue < 0.001). Human impacts (R 2 = 0.33, Pvalue < 0.001) and fire activities (R 2 = 0.31, Pvalue < 0.001) do not meet the significance threshold of 0.5 (Fig. 5). A further control is provided by the evolution of olive trees versus cereals (Fig. 4). The highest scores of Olea europaea are recorded long before the emergence of human impacts on the environment and also long before the oldest known date of human presence on the island [START_REF] Burjachs | Overview of environmental changes and human colonization in the Balearic Islands (Western Mediterranean) and their impacts on vegetation composition during the Holocene[END_REF]. Olive trees, therefore, are seemingly not a corollary of anthropogenic activities but rather of warmer conditions. The rank of pressures deriving from the CABFAC analyses clearly suggests the dominant role of environmental/climatic agents in mediating the island's palaeoecological dynamics (wet/dry and cool/warm).

Mallorca's mean seasonal temperatures

We first converted past and present-day datasets into mean seasonal temperatures to probe when climate change was most prominent (Fig. 6). According to the instrumental and reconstructed datasets, the sharpest "warming" contrasts between the last eight millennia and the period 1975-2018 occur in winter and autumn, with a respective difference in mean temperature of 4.8 ± 1.8°C and 3.9 ± 1.5°C. Spring and summer seasons display less pronounced shifts of 2.1 ± 1.3°C and 2.3 ± 1.4°C.

Mallorca's mean annual temperatures

We then calculated mean annual temperatures (MAT) for the instrumental and reconstructed datasets (Fig. 7A). The present-day MAT is 16.6 ± 0.8°C (Fig. 7A), with a mean high (Tmax) of 22.5 ± 0.7°C and a mean low (Tmin) of 10.6 ± 0.9°C. The warming trend is 0.55°C per decade for the period 1975-2018. A previous study based on the whole Balearic archipelago, for the period 1977-2006, showed that Tmin and Tmax are respectively 12.78°C and 21.76°C, with an estimated warming trend of 0.5°C per decade [START_REF] Homar | Recent trends in temperature and precipitation over the Balearic Islands (Spain)[END_REF]. The two warming estimates are similar, with a minor difference of 0.05°C. In comparison, warming on the island of Tenerife (Canary archipelago, North Atlantic Ocean) is much lower, with an average of 0.09 ± 0.04°C per decade since 1944 [START_REF] Martín | Assessment of global warming on the island of Tenerife, Canary Islands (Spain). Trends in minimum, maximum and mean temperatures since 1944[END_REF]. Focusing on the last eight millennia, the MAT is 13.3 ± 1.5°C, close to the present-day Tmin of the Balearic Islands (difference of ~0.5°C; [START_REF] Homar | Recent trends in temperature and precipitation over the Balearic Islands (Spain)[END_REF], but with significant temporal variations. The warmest period (15.6 ± 0.4°C, ~4000-3000 BP) is counterbalanced by the coolest phase (12.5 ± 1.1°C, ~8000-7000 BP), which yields a thermal adjustment that has not been observed since the 20th century. A comparison between present-day and past scores shows that only the ~4000-3000 BP values are similar to those recorded during the 1970s (MAT of 15.6 ± 0.2; Fig. 7A). Since the 1980s (MAT of 16.3 ± 0.5°C), according to our pollen-based model, the measured temperatures seem to have no analogue in the Holocene climate history of Mallorca. This finding is supported by the MAT for the last 400 years that clearly shows a major ecological break between the pre-industrial and industrial eras (Fig. 7B). The shift dates back to ~1900 ± 30 (according to our age-model) but changes have been particularly acute since the 1980s (instrumental data). It has been suggested that while the onset of the industrial-era warming can be dated to the 19 th century (~1830; [START_REF] Abram | Early onset of industrial-era warming across the oceans and continents[END_REF], the effects were more acute during the second half of the 20 th and 21 st centuries [START_REF] Brohan | Uncertainty estimates in regional and global observed temperature changes: A new data set from 1850[END_REF]IPCC, 2007;[START_REF] Abram | Early onset of industrial-era warming across the oceans and continents[END_REF]. In the Balearic Islands, scholars have demonstrated that climate-change anomalies have mainly increased since the 1980s [START_REF] Homar | Recent trends in temperature and precipitation over the Balearic Islands (Spain)[END_REF][START_REF] López Mayol | PREGRIDBAL 1.0: towards a high-resolution rainfall atlas for the Balearic Islands (1950-2009)[END_REF]. Our dataset (reconstructed and instrumental data) suggest a transition between the preand post-1980s.

Temperature anomalies

We subsequently calculated our reconstructed past climate data as anomalies compared to present-day conditions (Fig. 7C), using the period 1900-1980 (mean: 15 ± 0.6°C) as a baseline.

Two phases, centered on ~5800-5600 BP and ~4200-3000 BP, were warmer than the first period of industrial-era warming. By contrast, when the second period of instrumental data ; mean: 16.7 ± 0.6°C) is used as a baseline, all reconstructed temperatures from the preindustrial era were lower (i.e. cooler). Only two peaks at ~5665 ± 40 BP (with a range from -0.66 to 0.26°C) and ~4065 ± 30 BP (with a range from -0.39 to 0.31°C) are close to the 1981-2018 baseline (Fig. 8A). If the mean of the whole instrumental dataset is used as a baseline (Fig. 8B), the results illustrate the rapidity of the current changes compared to the climate history of Mallorca.

The warming pattern that emerges from our study suggests that the instrumental MAT is 3.3 ± 1.2°C warmer than the reconstructed MAT (for the last eight millennia). Furthermore, the reconstructed MAT is closer to the present-day Tmin (difference of 2.7 ± 1.2°C). With a mean warming rate of 0.55°C per decade, the present-day Tmin may become equivalent to the pre-/early industrial-era MAT in about five decades.

Mallorca's precipitation

We subsequently reconstructed the dynamics of the rainfall regime during the last 8000 years because variations in humidity were identified as the main forcing agent by the CABFAC analysis (Fig. 5). The MedCLIVAR consortium, and other scholars, have long insisted on the aridification of the Mediterranean area, with a likely increase in heat waves and droughts in the coming decades [START_REF] Homar | Recent trends in temperature and precipitation over the Balearic Islands (Spain)[END_REF][START_REF] Lionello | Introduction: Mediterranean climate -background information[END_REF][START_REF] Luterbacher | A review of 2000 years of paleoclimatic evidence in the Mediterranean[END_REF][START_REF] López Mayol | PREGRIDBAL 1.0: towards a high-resolution rainfall atlas for the Balearic Islands (1950-2009)[END_REF]. In effect, from 1900 to 2019, droughts and heat waves killed more than 62,000 people in the Mediterranean and affected a further > 9 million with total damages exceeding US$ 29 billion (EM-DAT, 2020).

We compared the trends highlighted in the reconstructed and instrumental data to assess whether major changes are detectable (Fig. 9). The two average values showed a significant shift, with a 33% drop in precipitation between the last 8000 years and the industrial era.

Focusing on the instrumental period, the driest years were recorded in 1983 (258.8 mm), 1995 (265.3 mm), 1999 (227.3 mm) and 2000 (255.9 mm). During the past 8,000 years, the driest millennium was ~4000-3000 BP (Fig. 9), with the lowest score recorded at ~3820 ± 30 BP (579.7 ± 11.8 mm). Although important variations were recorded during the Holocene, the highly variable annual precipitation patterns of the industrial era appear to be lower than reconstructed values for the last 8000 years.

Holocene climate dynamics

We finally compared our long-term Holocene reconstruction with other proxy records (Fig. 10).

The temperatures and precipitation data mimic variations in solar activity [START_REF] Solanski | Unusual activity of the Sun during recent decades compared to the previous 11,000 years[END_REF][START_REF] Vonmoos | Large variations in Holocene solar activity: Constraints from 10Be in the Greenland Ice Core Project ice core[END_REF][START_REF] Steinhilber | 9,400 years of cosmic radiation and solar activity from ice cores and tree rings[END_REF], particularly the warm and dry long interval that started ~5700 years ago and ended ~2400 years ago (Fig. 10).

Discussion

This paper reports a new climate reconstruction for the Balearics and underlines the fragility and vulnerability of these islands in the face of ongoing industrial-era warming. The long-term temperature and precipitation trend for Mallorca, based on both reconstructed and instrumental data, clearly depicts a marked offset between the present-day situation and the last 8000 years (Figs 6-9).

Caveats

Nonetheless, some caveats follow. A first limitation to the study results from missing data for a 35-yr interval spanning 1940-1975(Fig. 6) or 1940-1972 (Fig. 9) due to a gap in the instrumental records. A composite record, based on several Balearic weather stations, was not a viable option as no station has a complete time series for this interval. A second limitation relates to the 14 C chronology, which is based on calibrations with varying chronological ranges.

Although our reconstructions rely on these chronological data to provide a climate framework, the furnished climate signal remains unchanged. A third limitation relates to the influence of human activities on ecosystem dynamics. As suggested by the CABFAC analyses, human impacts are not the main pressures on ecosystem dynamics (for the pre-industrial era) but their influence is evidently recorded and must be accommodated in our reconstruction, because it is impossible to clearly separate natural forcing agents from anthropogenic activities. A final limitation concerns the number of endemic species in Mallorca and the potential non-analogue ecology when the present-day components are compared to the pre-anthropogenic era. Despite these limitations, our reconstructions are relevant because we focus on long-term trends, limiting over-interpretations resulting from a focus on particular events in the past. In effect, several similar climate reconstructions based on pollen (with the same limitations) have already been published and are relevant to ongoing global warming debates (e.g. [START_REF] Feurdean | Pollen-based quantitative reconstructions of Holocene climate variability in NW Romania[END_REF][START_REF] Guiot | The Mediterranean basin in a warmer and drier world: what can we learn from the past? Front[END_REF][START_REF] Mauri | The climate of Europe during the Holocene: a gridded pollen-based reconstruction and its multi-proxy evaluation[END_REF][START_REF] Guiot | Climate change: The 2015 Paris Agreement thresholds and Mediterranean basin ecosystems[END_REF][START_REF] Kaniewski | Cold and dry outbreaks in the eastern Mediterranean 3200 years ago[END_REF].

Local environmental significance

Our data provide a long-term perspective on current climate change on the Balearic Islands, as evoked by several authors [START_REF] Cuxart | Observed trends in frost and hours of cold in Majorca[END_REF][START_REF] Homar | Recent trends in temperature and precipitation over the Balearic Islands (Spain)[END_REF][START_REF] Jansà | Extension of summer climatic conditions into spring in the Western Mediterranean area[END_REF][START_REF] López Mayol | PREGRIDBAL 1.0: towards a high-resolution rainfall atlas for the Balearic Islands (1950-2009)[END_REF]. The rate and the magnitude of observed changes during the industrial era are significant when compared to the Holocene climate history of Mallorca (Figs 6789).

Moreover, the CABFAC analysis (Fig. 5) underscores that human activity was not the main pressure behind observed changes in the northeastern part of the island during the pre-industrial period. This argument was previously raised by [START_REF] Burjachs | Dinámica de la vegetación durante el Holoceno en la isla de Mallorca[END_REF]2017) for the same area.

Our study also stresses that, during the pre-industrial period, warmer and drier conditions developed from 5700 to 2400 BP (Figs 7 and 10). The same chronological boundary of 5700 BP was indicated by [START_REF] Bar-Matthews | Mid-Holocene climate variations revealed by highresolution speleothem records from Soreq Cave, Israel and their correlation with cultural changes[END_REF], with a major shift in aridity in the Levant, considered by the authors as the driest event of the mid-Holocene. [START_REF] Anderson | Postglacial history of alpine vegetation, fire, and climate from Laguna de Río Seco, Sierra Nevada, southern Spain[END_REF] also stressed this boundary for the Sierra Nevada (Laguna de Río Seco, southern Spain), where drier conditions commenced regionally after 5700 BP, marked by the onset of aridification of the flora of southern Spain. Our precipitation reconstruction tightly meshes with this observation (Fig. 9). According to [START_REF] Pérez-Sanz | Holocene climate variability, vegetation dynamics and fire regime in the central Pyrenees: the Basa de la Mora sequence (NE Spain)[END_REF], milder climate conditions characterized northern Spain from 8100 to 5700 BP. Mediterranean climate (warmer and drier) has spread since 5700 BP. [START_REF] Guiot | The Mediterranean basin in a warmer and drier world: what can we learn from the past? Front[END_REF] previously demonstrated that temperatures increased during the interval 5600-5700 BP. Although the period 5700-2400 BP was warm and dry, it never exceeded the values recorded since the end of the 20th century and those of the 21st century CE.

A major consequence of the ongoing warming is falling water availability in island contexts, whose environmental pressure absorption capacities are much lower than for mainland areas [START_REF] Marriner | Fire as a motor of rapid environmental degradation during the earliest peopling of Malta 7500 years ago[END_REF]. As stated by [START_REF] Kent | Tourism and sustainable water supply in Mallorca: a geographical analysis[END_REF], water shortages are a key environmental problem in Spanish islands, attracting both popular and political attention. Our pollen-based climate reconstruction supports these challenges, suggesting that current precipitation levels are substantially lower than reconstructed values for the Holocene (Fig. 9). For much of the Mediterranean, recent studies have indicated that, for each 1°C of global warming, mean rainfall will probably decrease by about 4% [START_REF] Cramer | Climate change and interconnected risks to sustainable development in the Mediterranean[END_REF][START_REF] Lionello | The relation between climate change in the Mediterranean region and global warming[END_REF].

Global significance

This study suggests two different climate pacemakers for the Balearic Islands. The climate oscillations playing out during the Holocene seem to result from an ocean-atmosphere coupling with a clear imprint of solar activity as the main driver. A connection between solar irradiance and water balance has previously been suggested for Spain (Morellón et al., 2011). On Mallorca, the warmest and driest interval (2400-5700 BP) during the pre-industrial period correlates with increases in solar activity (Fig. 10) as reconstructed from total solar irradiance (∆TSI; [START_REF] Steinhilber | 9,400 years of cosmic radiation and solar activity from ice cores and tree rings[END_REF], sunspot numbers [START_REF] Solanski | Unusual activity of the Sun during recent decades compared to the previous 11,000 years[END_REF] and from the solar modulation function Φ [MeV] based on cosmogenic radionuclides ( 10 Be production rates; [START_REF] Vonmoos | Large variations in Holocene solar activity: Constraints from 10Be in the Greenland Ice Core Project ice core[END_REF]. The industrial era shows a different trend with a rate and magnitude of changes that exceed the long-term natural variability (Fig. 7; IPCC, 2014IPCC, , 2018)). Changes in solar activity appear to have been a primary pacemaker of long-term climate variations in Mallorca. By contrast, anthropogenic change has been the primary pacemaker since the industrial era, with a clear acceleration of warming and aridification trends.

Social and political implications

In the coming decades, warming will have a significant impact on the local economy of the Balearics, which is heavily dependent on tourism (e.g. [START_REF] Kent | Tourism and sustainable water supply in Mallorca: a geographical analysis[END_REF][START_REF] Amelung | Mediterranean tourism: exploring the future with the Tourism Climatic Index[END_REF]. Since the 1960s, tourism has reached an apex on Mallorca and Ibiza, with approximately 11 million visitors in 2000-2001 and more than 12 million in 2011 [START_REF] Hof | The linkages between real estate tourism and urban sprawl in Majorca (Balearic Islands, Spain)[END_REF]. As shown by a recent tourism model, climate is an important variable in estimating the economic sustainability of these islands [START_REF] Álvarez-Díaz | Forecasting British tourist arrivals in the Balearic Islands using meteorological variables[END_REF]. Different climate change scenarios coupled with the Tourism Climatic Index suggest that, during the 21 st century, the Balearic Islands, in addition to the wider Mediterranean region, will probably become unsustainably hot during the summer months, while, by contrast, northern Europe will be climatically more attractive for tourism [START_REF] Amelung | Mediterranean tourism: exploring the future with the Tourism Climatic Index[END_REF]. Socially and economically, the elucidated changes could become detrimental to tourist-based economies throughout the Mediterranean [START_REF] Amelung | Mediterranean tourism: exploring the future with the Tourism Climatic Index[END_REF].

Water management on Mallorca is already under pressure from increasing variability in seasonal and inter-annual rainfall patterns [START_REF] Candela | Impact assessment of combined climate and management scenarios on groundwater resources and associated wetland (Majorca, Spain)[END_REF][START_REF] Homar | Recent trends in temperature and precipitation over the Balearic Islands (Spain)[END_REF][START_REF] López Mayol | PREGRIDBAL 1.0: towards a high-resolution rainfall atlas for the Balearic Islands (1950-2009)[END_REF], creating an unbalanced distribution of water/groundwater resources (domestic supply and irrigation-based agriculture versus natural ecosystem replenishment).

This situation is likely to present serious challenges to the future sustainability of the island and its economic model based on mass tourism [START_REF] Candela | Impact assessment of combined climate and management scenarios on groundwater resources and associated wetland (Majorca, Spain)[END_REF][START_REF] Hof | Urban and tourist land use patterns and water consumption: Evidence from Mallorca, Balearic Islands[END_REF].

Reducing the impact

Climate change concerns have led to major political changes on the Balearic Islands. Its parliament, for instance, has recently passed a law on energy transition and climate change that commits local tourism to draw all of its energy needs from renewable sources by 2050. This study reinforces the idea that the Balearics must rapidly adopt such measures to ensure the medium-to long-term viability of its economic model based on mass tourism.

Conclusions

In conclusion, the analyses presented above provide answers to the two research questions posed in this paper. First, through the prism of Mallorca, results have shown that the industrialera warming and aridification of Mediterranean islands deviates markedly from the long-term Holocene trends. The instrumental MAT appears 3.3 ± 1.2°C warmer than the reconstructed MAT for the last eight millennia, with a 33% drop in precipitation between the last 8000 years and the industrial era. Second, this study has also demonstrated the importance of palaeoecological archives in contextualizing current climate anomalies. It has revealed the slow pace of climate dynamics on Mallorca over the past 8000 years, both in rate and magnitude, compared to the speed of changes during the industrial era. This contrast also suggests that, despite human colonization of the Balearic Islands around 4400 years ago, the "Anthropocene" of accelerating human interactions with Earth systems closely intertwines with the industrial era. A human imprint on climate is clearly recorded in the data since 1900 CE, with a sharp acceleration in temperatures and aridity since the 1980s. When temperatures and precipitation during the last 40 years are used as benchmarks (period 1980-2018), the pre-industrial climate of Mallorca appears cooler and wetter than anything that has been recorded since. [START_REF] Solanski | Unusual activity of the Sun during recent decades compared to the previous 11,000 years[END_REF][START_REF] Steinhilber | 9,400 years of cosmic radiation and solar activity from ice cores and tree rings[END_REF]. The yellow shaded part underscores the warmer and driest phase during the last 8,000 years.

FiguresFigure 1 .

 1 FiguresFigure 1.Geographical location the Albufera Natural Park in Mallorca. The study area is

  Figure 1. Geographical location the Albufera Natural Park in Mallorca. The study area is denoted by a red star. The geomorphological map of the Bay of Alcúdia and the swamp area (Albufera Natural Park) is adapted from Giaime et al. (2017).

Figure 2 .

 2 Figure 2. Age-model for the core POL12. The lithology of the core is described in relation to depth. The radiocarbon dates are depicted as intercepts and 2-sigma calibrations (95% probability).

Figure 3 .

 3 Figure 3. Detailed pollen diagram (arboreal pollen -AP; non-arboreal pollen -NAP) of the core POL12. The arboreal and herbaceous taxa are plotted on a linear age-scale. The vertical lines indicate the main transitions.

Figure 4 .

 4 Figure 4. Fire activity, agriculture (Poaceae cerealia) and olive trees (Olea europaea) plotted on a linear age-scale. The vertical lines denote the results of the homogeneity test (each average is indicated by the letters mu).

Figure 5 .

 5 Figure 5. Environmental regression models from the Albufera Natural Park. The selected environmental data are indicated in each graph. The consistency of each reconstructed model is indicated by the R2. The lower graph corresponds to a ranking of the environmental data based on the R2 values.

Figure 6 .

 6 Figure 6. Long-term trends in Mallorca's mean seasonal temperatures. The chronological scales are expressed in years BP and Before Common Era -Common Era (BCE -CE). The modelled temperatures for each season are displayed, with the standard deviation for each data point. The univariate analyses (mean, 25 th percentile-Q1, median-Q2, 75 th percentile-Q3, minimum value, maximum value) are summarized using a boxplot. The period with no data (1940 to 1975 CE) is denoted by a vertical grey bar.

Figure 7 .

 7 Figure 7. Long-term evolution of mean annual temperatures on Mallorca. (A) The modelled and measured values are presented as curves and smoothing splines with the standard deviation for each reconstructed data point. The mean value for each dataset is represented in purple. (B) Focus on the period 1600-2018, showing the mean annual temperatures and associated standard errors. The curves are continuous despite the absence of data between 1940 and 1975. (C) Temperature anomalies for the pre-industrial era with the mean 1900-1980 used as the baseline.

Figure 8 .

 8 Figure 8. Long-term temperature anomalies on Mallorca. (A) Temperature anomalies for the pre-industrial era, with the mean 1981-2018 (16.71°C) used as the baseline. A smoothing spline and a polynomial curve were fitted to the sequence. (B) Temperature anomalies for the pre-industrial and industrial eras with a mean calculated on the whole dataset (14.12°C). This value was used as the baseline. A smoothing spline and a polynomial curve were added to the sequence. The age is expressed as b2k.

Figure 9 .

 9 Figure 9. Long-term trends in Mallorca's precipitation. The chronological scales are expressed in years BP and the Common Era (CE). The reconstructed annual precipitation values are displayed with the standard deviation. The univariate analyses (mean, 25th percentile-Q1, median-Q2, 75th percentile-Q3, minimum value, maximum value) are summarized using boxplots. The period with no data (1940 to 1972 CE) is denoted by a vertical grey bar. Two different vertical scales are used (right: reconstructed values; left: instrumental data).

Figure 10 .

 10 Figure10. Temperature and precipitation trends on Mallorca during the pre-industrial period compared with the long-term solar activity[START_REF] Solanski | Unusual activity of the Sun during recent decades compared to the previous 11,000 years[END_REF][START_REF] Steinhilber | 9,400 years of cosmic radiation and solar activity from ice cores and tree rings[END_REF]. The

Acknowledgments

This manuscript is dedicated to the memory of our friend Valérie Trichon. She was involved in the field campaign and in the design of this study. This work constitutes part of the activities of the Equip de Recerca Arqueològica i Arqueomètrica de la Universitat de Barcelona (ERAAUB), Consolidated Group (2017 SGR 1043), supported by the Comissionat per a Universitats i Recerca del DIUE de la Generalitat de Catalunya. We are also indebted to the Conselleria de Medi Ambient del Govern de les Illes Balears, Consell de Mallorca, and Albufera Natural Park for granting permission for core POL12. The Agencia Estatal de Meteorologia kindly provided instrumental climate data. The coring campaign was supported by the Consorci de la Ciutat Romana de Pollèntia within the framework of archaeological research on the Roman city of Pollentia (Alcúdia, Mallorca). We are indebted to Esther Chávez, co-director of the project and the rest of the team for their support. Financial support was provided by the Ecosystèmes urbains CNRS 2019, URBMED program. This work is a contribution to Labex OT-Med (n° ANR-11-LABX-0061) and has received funding from the Excellence Initiative of Aix-Marseille University -A*MIDEX, a French "Investissements d'Avenir" project. This research is also part of the project ARCHREMOTELANDS (HAR2017-83335-P) funded by the Ministerio de Ciencia, Innovación y Universidades, with a contribution from the European Regional Development Fund of the European Commission.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.