K R Maiyelvaganan

S Kamalakannan

S Shanmugan

M Prakash
email: prakashspm@gmail.com

François-Xavier Coudert
email: fx.coudert@chimieparistech.psl.eu

M Hochlaf
email: hochlaf@univ-mlv.fr

$ S Kamalakannan

$ S Shanmugan

$ M Prakash
email: prakashm4@srmist.edu

F.-X Coudert

Identification of a Grotthuss Proton Hopping Mechanism at Protonated Polyhedral Oligomeric Silsesquioxane

Keywords: DFT, AIM, AIMD, Grotthuss, POSS, Protonated Water Clusters, IR Signatures, Eigen and Zundel

HAL is

Introduction

Proton transport and conduction of proton through water chains occur in a variety of atmospheric, biological and chemical processes. [START_REF] Schmitt | Multistate empirical valence bond model for proton transport in water[END_REF][START_REF] Voth | Computer simulation of proton solvation and transport in aqueous and biomolecular systems[END_REF][START_REF] Prakash | Stepwise hydration of protonated carbonic acid: a theoretical study[END_REF][START_REF] Achtyl | Aqueous proton transfer across single-layer graphene[END_REF][START_REF] Boussouf | Hochlaf Structure, spectroscopy and bonding within the Zn q+ -imidazolen (q = 0, 1, 2; n = 1 -4) clusters and implications for ZIFs and Zn-enzymes[END_REF][START_REF] Weichselbaum | Origin of proton affinity to membrane/water interfaces[END_REF][START_REF] Wang | Stability of hydronium ions in zeolite channels[END_REF] The formation of hydronium ion and the mobility of the proton is on the origins of these phenomena. For explanation, two scenarios were suggested: (i) Back in 1806, von Grotthuss 8 proposed a mechanism, which makes it possible to account for the high ionic molar conductivity of hydronium and hydroxide ions and thus of water liquid conductivity. Within this mechanism, the ion does not have to move entirely through the liquid to carry its charge, but only passes the excess proton. Indeed, each oxygen atom passes and at the same time receives a hydrogen atom. Obviously, Grotthuss-type of proton transfer is rapid and efficient, where water facilitates an efficient charge separation when the proton leaves the acid. If it occurs at the interfaces of liquid water with biological entities, or catalysts or material interfaces, it changes drastically the kinetics and thus the outcomes and applications of such compounds. (ii) Alternatively, Kreuer, Rabenau, and Weppner [START_REF] Kreuer | Vehicle mechanism, a new model for the interpretation of the conductivity of fast proton conductors[END_REF] established in 1982 a mechanism (denoted as "vehicle mechanism") where the protons attached to the protonated species and the neutral surrounding uncharged molecules diffuse together. The excess proton travels on top of the host molecule through the solvent. The vehicle mechanism corresponds hence mostly to the proton conduction in aqueous systems.

The solid-liquid interface proton transfers are important for fuel cell applications such as the proton conduction in electrolyte membranes. The commonly used Nafion-based membranes have however some drawbacks limiting their applications. [START_REF] Mauritz | State of understanding of Nafion[END_REF] Exploiting alternative composite fuel cell membranes without compromising proton conductivity is an active field. For that purposes, polyhedral oligosilsesquioxanes (POSS) as nano-additives are expected to enhance the properties of such cells. This is also associated with good thermal, chemical and mechanical stabilities because of the rigidity of the inner silicon and oxygen frameworks. [START_REF] Li | Polyhedral oligomeric silsesquioxane (POSS) polymers and copolymers: A Review[END_REF]12 Indeed, earlier experimental studies highlighted the benefits of POSS framework such as its non-volatility and its extensively stable nature at high temperature, pressure, odorless, extreme acid-base condition. It is also ecologically friendly. [START_REF] Shea | Bridged polysilsesquioxanes Molecular-engineered hybrid organic-inorganic materials[END_REF][START_REF] Kuo | POSS related polymer nanocomposites[END_REF][START_REF] Cordes | Recent developments in the chemistry of cubic polyhedral oligosilsesquioxanes[END_REF] It is very challenging to probe the mobility of protons at the solid-liquid interface using experimental techniques. Instead, insights into the proton transfer mechanism, pathway and conduction ability of materials can be understood by a variety of computational and dynamical techniques. [START_REF] Cuma | A multi-state empirical valence bond model for weak acid dissociation in aqueous solution[END_REF][START_REF] Prakash | Structure, stability and spectral signatures of monoprotic carborane acid-water clusters (CBWn, where n= 1-6)[END_REF][START_REF] Quaranta | Proton-transfer mechanisms at the water-ZnO interface: The role of presolvation[END_REF][START_REF] Leung | Concerted metal cation desorption and proton transfer on deprotonated silica surfaces[END_REF][START_REF] Wang | A robust zirconium amino acid metal-organic framework for proton conduction[END_REF][START_REF] Mileo | Achieving superprotonic conduction with a 2D fluorinated metal-organic framework[END_REF][START_REF] Lowe | Acid-base dissociation mechanisms and energetics at the silica-water interface: An activationless process[END_REF] For instance, the proton transfer in aqueous media and proton dissociation ability from acidic materials have been studied using the density functional theory (DFT) and ab initio electronic structure approaches. [START_REF] Wang | An ab initio study of the primary hydration and proton transfer of CF3SO3H and CF3O (CF2)2 SO3H: Effects of the hybrid functional and inclusion of diffuse functions[END_REF][START_REF] Losch | Proton mobility, intrinsic acid strength, and acid site location in zeolites revealed by varying temperature infrared spectroscopy and density functional theory studies[END_REF] There are however few studies that have been reported on POSS-based materials for proton conduction ability of sulfonic acid or phosphonic acid derivatives. [START_REF] Xing | Improved performance of sulfonated polyarylene ethers for proton exchange membrane fuel cells[END_REF][START_REF] Hartmann-Thompson | Proton-conducting polyhedral oligosilsesquioxane nanoadditives for sulfonated polyphenylsulfone hydrogen fuel cell proton exchange membranes[END_REF] At present, quantum chemical calculations using ab initio and DFT methods are performed for understanding the acidic nature of POSS-H + cage as well as the proton transfer, the proton conduction pathway and proton affinity in this solid-liquid interface. This is done after structural changes and symmetric, asymmetric, proton oscillation and strong hydrogen bond (H-bond) stretch vibrational analyses of the protonated form of polyhedral oligomeric silsesquioxane with water (W) clusters (POSS-H + @Wn, where n = 1 -6). In addition, we show that the protonated species participate into a proton hopping mechanism at this solid-liquid interface via Grotthuss-type mechanism. The excess proton dissociation from POSS-H + molecule to water molecule is associated with the successive formation of Eigen/Zundel/Eigen cations. Their vital role in proton transfer process at microscopic level is elucidated. This process was already observed at silica-water interface, [START_REF] Lowe | Acid-base dissociation mechanisms and energetics at the silica-water interface: An activationless process[END_REF] but it is found here for the first time for an interface composed by a nano-cluster and water molecules.

Methodologies

Characterization of POSS-H + @Wn (n = 1 -6) clusters

The initial structure of POSS was taken from crystal database of POSS framework. For the initial structure of POSS-H + , we used the procedure described in Figure S1. All geometry optimizations were done at the DFT level, using the B3LYP functional (Becke's three-parameter hybrid exchange functional and Lee-Yang-Parr correlation functional), in conjunction with the 6-311++G** basis set. [START_REF] Becke | Density-functional exchange-energy approximation with correct asymptotic behavior[END_REF][START_REF] Becke | Density-functional thermochemistry. III. The role of exact exchange[END_REF] Earlier reports revealed that this level of basis set is suitable to study accurately the structure, stability and spectral signatures of protonated water clusters and proton transfer mechanism. [START_REF] Kuo | Comprehensive analysis of the hydrogen bond network morphology and OH stretching vibrations in protonated methanol-water mixed clusters, H + (MeOH)1 (H2O)n (n = 1-8)[END_REF][START_REF] Jiang | Hydrogen Bond Rearrangements and Interconversions of H + (CH3OH)4 H2O Cluster Isomers[END_REF][START_REF] Wu | Protonated clathrate cages enclosing neutral water molecules: H + (H2O)21 and H + (H2O)28[END_REF][START_REF] Lin | Vibrational predissociation spectra and hydrogen-bond topologies of H + (H2O)[END_REF] Geometry optimizations were systematically followed by frequency calculations carried out at the same level of theory, to confirm the minimal nature of the stationary points on the corresponding potential energy surface. All electronic structure calculations were performed using GAUSSIAN 16. 33 GaussView 6.0 and Chemcraft softwares were used for visualization and geometrical analysis. 34,[START_REF] Zhurko | Chemcraft Program[END_REF] We also performed atoms in molecules (AIM) analysis using the AIM2000 package [START_REF] Biegler-Konig | [END_REF] to quantify the various types of noncovalent interactions in POSS-H + @Water interfaces (see Supporting Information for full details).

Ab initio molecular dynamics (AIMD) methodology

The hybrid Gaussian and plane wave method GPW was used as implemented in the CP2K software. [START_REF] Vandevondele | QUICKSTEP : Fast and Accurate Density Functional Calculations Using a Mixed Gaussian and Plane Waves Approach[END_REF] The production simulations were performed in the constant-volume (N, V, T) ensemble with fixed size and shape of the unit cell. A time step of 0.5 fs was used in the MD runs, and the temperature was controlled by using a canonical sampling through velocity rescaling thermostat with a time constant of 1 ps. The exchange-correlation energy was evaluated in the Perdew-Burke-Ernzerhof (PBE) approximation [START_REF] Perdew | Generalized Gradient Approximation Made Simple[END_REF] and the dispersion interactions were treated by incorporating Grimme corrections [START_REF] Grimme | Semiempirical Hybrid Density Functional with Perturbative Second-Order Correlation Semiempirical Hybrid Density Functional with Perturbative Second-Order Correlation[END_REF][START_REF] Grimme | A Consistent and Accurate Ab Initio Parametrization of Density Functional Dispersion Correction ( DFT-D ) for the 94 Elements H-Pu[END_REF] at the DFT-D3 level. The Quickstep module uses a multigrid system to map the basis functions onto. The default number of 4 different grids was used, along with a plane-wave cutoff for the electronic density of 800 Ry (suited for the presence of O atoms), and a relative cutoff of 40 Ry. Valence electrons were described by double-ζ valence polarized basis sets and norm-conserving Goedecker-Teter-Hutter pseudopotentials all adapted for PBE (DZVP-GTH-PBE) for H, C, O and Si. In addition, to identify the proton transfer reaction pathways and to get more insights on the ionic features of protonated water clusters at this interface, we have performed Nudged elastic band (NEB) calculations with PBE+D3 method, where the atoms were described with the triple-ζ valence polarized basis sets (TZVP-GTH-PBE) for H, O and Si atoms.

The present methodology relies on the GGA-level exchange-correlation functional PBE, with dispersion corrections, to describe the water and water-POSS interactions. While this level of theory is known to slightly over structure water systems (including the liquid phase), this is in balance with a lower computational cost that allows for better statistical sampling. Given the system sizes and MD length involved in this study, hybrid functionals could not be used. Furthermore, we note that GGA-level exchange-correlation functionals are widely used for the description of water near interfaces, including in the description of Grothuss proton transfer, and are qualitatively well validated. This includes a wide range of systems and phenomena, such as chemically reactive graphene oxide in water [START_REF] Mouhat | Structure and chemistry of graphene oxide in liquid water from first principles[END_REF] , inorganic or mineral surface/water interface [START_REF] Motta | Ab Initio Molecular Dynamics Study of the AlOOH Boehmite/Water Interface: Role of Steps in Interfacial Grotthus Proton Transfers[END_REF][START_REF] Ngouana-Wakou | An Atomistic Description of the γ-Alumina/Water Interface Revealed by Ab Initio Molecular Dynamics[END_REF] , boron nitride and graphene [START_REF] Grosjean | Spontaneous liquid water dissociation on hybridised boron nitride and graphene atomic layers from ab initio molecular dynamics simulations[END_REF] , etc. … Figure 1 displays the optimized geometries of POSS-H + @Wn (n=1-6) clusters. This figure gives also the main distances between POSS-H + and water molecules. For POSS-H + and POSS-H + @W1 and POSS-H + @W2 a unique form is found. For POSS-H + @W3, POSS-H + @W4 and POSS-H + @W5, four stable configurations are identified. For POSS-H + @W6, we obtain six stable structures. Within all species, the POSS moiety adopts a similar cube-like geometry, which is analogous to the cubane framework: Si atoms are placed at the corners of the cube and bridging oxygen atoms are in-between two Si atoms. For the optimized structure of POSS, the angles are 149° and 109° for the Si-O-Si and O-Si-O linkages, respectively (Figure S1). These parameters are in very close agreement with the experimental crystal structure of POSS. [START_REF]Applications of Polyhedral Oligomeric Silsesquioxanes[END_REF] Upon protonation, the additional proton links to an oxygen atom and the corresponding O-H bond distance is 0.970 Å. The corners of the cuboid structure in the vicinity of this protonated oxygen are affected by the presence of the H + , and we compute Si-O-Si ~132° and O-Si-O ~105°.

Results

Structures of POSS-H + @Wn (n = 1 -6) clusters

Upon addition of water molecules to POSS-H + , we found that the proton strongly interacts with the primary shell formed by water. Depending on the geometry and number of water molecules, we observe either the lengthening of the distance between the proton and the oxygen of POSS, or the proton transfer from POSS-H + to the water cluster. For instance, this distance increases from 0.970 Å in POSS-H + to 1.019 Å in POSS-H + @W1 and to 1.074 Å in POSS-H + @W2. At the same time, the H + •••Ow distance (Ow is the O of the water molecule) decreases and the H-bond becomes shorter and stronger. Additional water molecules attached to POSS-H + @W2 lead to proton dissociation from POSS-H + and attachment to the neighboring water molecule. This is associated with the formation of either an Eigen cation core or a Zundel cation core at the water chain. The latter requires minimum of four water molecules interacting with POSS-H + . This is the case for all POSS-H + @Wn (where, n = 3-6) structures except the one denoted as POSS-H + @W3d (Figure 1). These cations are stabilized by linear, branched or cyclic H-bond patterns. Specifically, we notice an excess proton on water clusters. When possible, we assist to the complete solvation of the Eigen cation to further stabilize this ion. For the Zundel ion, the excess proton is localized in the middle of the two water molecules with the formation of Zundel core (highlighted by the blue dashed rectangle in Figure 1) and the proton donor and acceptor sites in the rigid POSS materials. This is associated with a planar water six-membered ring. The excess proton is dangling between two water molecules (H2O•••H + •••OH2) with an Ow -Ow distance of ~2.4 Å, which is typical for Zundel ion. [START_REF] Park | Eigen and Zundel forms of small protonated water clusters: Structures and infrared spectra[END_REF]

Binding energies

To get insights on the stability and the mechanism of the proton transfer on the POSS-H + @Wn (n = 1-6) clusters, we computed their binding energies (BE). The BEs were calculated using the supermolecular approach and corrected for basis set superposition error (BSSE) employing the counterpoise (CP) procedure suggested by Boys and Bernardi 47

      - =  = n i i Cluster E E BE 1 (1)
where, ECluster is the total energy of the POSS-H + @Wn cluster, Ei is the energy of the monomers (i.e. POSS, POSS-H + , H2O or H3O + /H2O5 + ). To quantify the dispersion energy in these clusters, we have computed energetics both with (DFT-D3) and without (DFT) the inclusion of Grimme's dispersion correction on the DFT optimized structures. [START_REF] Grimme | A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu[END_REF] The BSSE-corrected BEs of various POSS-H + @Wn (n = 1-6) clusters are given in Table 1. As can be seen there, the BEs including the D3 corrections are larger (in absolute value), and they will be considered for the discussion as being in principle more accurate. For POSS-H + @Wn (n ≥ 3) clusters, the Eigen core cyclic clusters are more stable than the other isomers. 1 shows that the |BE|'s increase from 20 to 120 kcal/mol when increasing the number of water molecules. This is due to the improved stabilization of the larger complexes by H-bonds, yet typical H-bond energies are in the range of ~ 5-10 kcal/mol, i.e., much less than the bond energies of our interface material. The existence of Eigen core with water clusters is significantly more stable than the non-Eigen containing complexes (e.g. POSS-H + @W3d; Figure 1 and Table 1). Also, the Eigen core clusters are more stable than the corresponding Zundel type having the same H-bonds; this is due to the positive charge delocalization (i.e. equally shared) between two water molecules. Thus, the existence of Eigen core with water clusters significantly stabilizes this solid@water interface model. In contrast, the formation of Zundel cation relatively destabilizes this interface. optimized at the MP2 level of theory. 49

IR spectra

For better highlighting the effects of clustering of POSS-H + and to help their identification in the laboratory, we generated the IR spectra of POSS-H + @Wn. The vibrational frequencies are reported after scaling the harmonic frequencies by a factor of 0.9577, for ease of comparison with experimental data. [START_REF] Nibu | IR spectroscopy of 2-fluoropyridine-water clusters in the electronic excited states[END_REF] The full analysis of the data is given in the Supporting Information (Table S2). For illustration, Figure 2 presents the simulated spectrum of POSS-H + @Wn (n = 1-4) clusters. The calculated characteristic intense peaks of POSS framework materials for Si-O-Si bending and Si-H stretching are observed at ~1104, 2260 cm -1 , respectively in good agreement with the corresponding experimental stretching frequencies (1100 cm - 1 and 2293 cm -1 ). [START_REF]Applications of Polyhedral Oligomeric Silsesquioxanes[END_REF] Similarly, the calculated IR spectrum of O-Si-O linkage of νas (380 cm -1 ), νss (533 cm -1 )

and Si-H oscillations peak (862 cm -1 ) values are also in very good agreement with experimental values which are 389, 557 and 860 cm -1 , respectively. This confirms the good accuracy of the present predictions.

Note that the vibrational results analyses are in good agreement with the geometrical parameters and AIM electron density topography analyses. (See Supporting Information for details). The calculated IR spectra of POSS-H + @Wn (n = 1-4) clusters are shown in Figure 2. For POSS-H + @W1, the peaks appear at 3640 and 3729 cm -1 due to the symmetric and asymmetric O-H stretching of water molecule in two free O-H groups. Compared to a free water molecule, these values are slightly red shifted by ~ 16 and 28 cm -1 respectively. The intense band at 2701 cm -1 corresponds to the H-bonded acceptor water (A) molecule with the POSS-H + moiety. The corresponding red shifted value is 897 cm -1 .

Such large redshift is a signature of the strength of the H-bonds in these clusters. In POSS-H + @W2 cluster, the H-bonded AD water molecule gives intense band at 1925 cm -1 and resulting in a red shift of 1674 cm -1 .

The large red shift is due to the presence of charge movement towards H-bonded water chain and the subsequent O-H bond order decrease. Earlier reports revealed that the calculated red shift value of acceptor (A) is significantly higher than the less acceptor (A) or donor (D) types of H-bond. [START_REF] Prakash | Stepwise hydration of protonated carbonic acid: a theoretical study[END_REF][START_REF] Prakash | Structure, stability and spectral signatures of monoprotic carborane acid-water clusters (CBWn, where n= 1-6)[END_REF][START_REF] Prakash | Ab initio and density functional theory (DFT) studies on triflic acid with water and protonated water clusters[END_REF] For POSS-H + @W3 clusters, there is an Eigen core stretching band around 2347 cm -1 , which is slightly red shifted with real Eigen core in aqueous medium (i.e. 2420 cm -1 ). Both AD and A types H-bond stretching frequencies are observed for water trimer.

The IR spectra of POSS-H + @W4 isomers exhibit both Eigen and Zundel features. The bands at 1982 and 2373 cm -1 for POSS-H + @W4a and POSS-H + @W4b clusters are due to Eigen core stretching A and AD water molecules. The larger red shift corresponds to linear type clusters, which is more favorable for the proton transport through water chain. The excess proton always stabilizes in the middle of the water chain, irrespective the n value. Also, we observe Zundel features in POSS-H + @W4c and POSS-H + @W4d clusters. The respective O-H + -O oscillation vibrational modes appear at 1688 and 1685 cm -1 . In addition, low frequency modes are computed at 904 and 882 cm -1 for the Zundel characteristic peaks. Since the linear chain Eigen water interactions are stronger than the cyclic and branched water clusters, a large red shift is computed.

Dynamics of proton transfer at the POSS-Water interface

To probe the dynamics of the POSS-H + @water proton transfer, we then performed ab initio molecular dynamics (AIMD) simulations, also known as first-principles MD. The protocol followed was the following: (i) generate an initial configuration by random packing the POSS-H + with 53 water molecules in a 12 Å cubic box; (ii) perform a short geometry optimization of the water molecules, to relax possible unphysical configurations due to random packing; (iii) run a constant-pressure dynamics (NPT MD) with the POSS-H + frozen, so that the density can relax to its equilibrium value, and the water can reorient and solvate the POSS-H + ; (iv) run a "production" constant-temperature dynamics (NVT MD) at the average volume, allowing all atoms to move freely. This computational protocol ensures that the POSS-H + @water dynamics start from a physically meaningful configuration, and we can focus observation on the dynamics of the proton. We studied the system at two temperatures, 100 K and 300 K, and performed 3 simulations at each temperature. In the majority of cases, we observe a spontaneous transfer for H + from POSS-H + to the water. This is clearly manifested in the plot of O*-H + distance over time (Figure 3), if we label O* the POSS oxygen that initially bears the proton: in this 100 K dynamics, we see that the initial structure (solvated POSS-H + ) is stable in water for the first 2.8 ps of the trajectory, after which we observe a proton transfer to form POSS + H3O + . We also see that the hydronium cation remains weakly hydrogen-bonded to the neutral POSS, even in solution, as the O*-H + distance fluctuates around 1.75 Å (characteristic of a H bond). Further simulation shows that the excess proton can then diffuse by a Grotthuss-like mechanism, which AIMD is known to be well-suited to reproduce, both in bulk [START_REF] Corradini | Carbon dioxide transport in moltencalcium carbonate occurs through an oxo-Grotthuss mechanism via a pyrocarbonate anion[END_REF] and at interfaces 52 . We have observed several events of proton transfer between water molecules (4 times during the dynamics at 300 K), and the simulation was not continued after that: as the excess proton goes further away from the POSS, its dynamics will be that of the excess proton in liquid water.

In some cases, we have also observed reactive events involving the POSS itself, i.e. strong deformation leading to bond breaking inside the Si-O framework, in a relatively short time (100 to 300 fs). This was observed in one case at each temperature, and while the details of the dynamics leading to this chemical event vary, the even itself is the same in both cases: if we name O* the POSS oxygen that initially bears the proton, then one of the Si-O* bonds breaks (see Figure 4). On one side it forms -Si-OH, which reorients to some extent but does not react further. On the other side, the nowundercoordinated Si will 13 bond to one of the water molecules, forming -Si-OH2 + ; this will then very rapidly transfer its excess proton to a neighboring water molecule, ultimately leading to the formation of -Si-OH and H3O + (as depicted on Figure 4). The excess proton then diffuses by a Grotthuss-like mechanism. Therefore, the AIMD study confirms the conclusions of a Grotthuss-type mechanism for transfer of proton in this POSS-H + @water system. It also suggests that, in liquid water (or in presence of more than a few water molecules), there are other possible evolutions of the POSS-H + system, including reactivity of the POSS framework itself. This will need to be further investigated by systematic exploration of the potential energy surface of the POSS-H + @Wn configurations taking into account such opening, or by more computational expensive free energy methods based on the AIMD description. To get information about the conduction of proton, we performed nudged elastic band analyses for POSS-H + @W3a, POSS-H + @W4a and POSS-H + @W4c clusters. Computations suggest that POSS-H + @W3a and POSS-H + @W4a isomers are more stable than the other clusters. Thus, we selected POSS-H + @W4c model for cyclic pathway. The detailed reaction pathway with energy barrier analysis are shown in Figure 5 and Figure S6 for POSS-H + @W3, POSS-H + @W4a and POSS-H + @W4c clusters. Our study reveals that the activation energy required for Eigen cation is ~240 kcal/mol, whereas Zundel ion formation requires less energy (i.e. ~184 kcal/mol). This documents the destabilization of Zundel ion at this interface favoring the formation of an Eigen core. This is also in line with the BE analysis and IR shifting values discussed in previous section. The present work identifies a proton transfer pathway through water wire, which corresponds to the well-known Grotthuss mechanism as illustrated in Figure 6. This figure shows the important steps with formation of Eigen  Zundel  Eigen features during hydration process. The present geometrical analysis reveals that the formation of Eigen feature requires a minimum of three water molecules whereas Zundel feature requires a minimum of four water molecules in the POSS-H + @Wn clusters. Similar observation also found for the experimental and theoretical studies on proton exchange membrane fuel cell (PEMFC) materials (i.e. Nafion), [START_REF] Li | theoretical studies of proton transfer in water and model polymer electrolyte systems[END_REF] a hydrated form of protonated carbonic acid (PCA), [START_REF] Prakash | Stepwise hydration of protonated carbonic acid: a theoretical study[END_REF] carborane acid (CA) [START_REF] Prakash | Structure, stability and spectral signatures of monoprotic carborane acid-water clusters (CBWn, where n= 1-6)[END_REF] and trifluoromethanesulfonic acid (TA) [START_REF] Prakash | Ab initio and density functional theory (DFT) studies on triflic acid with water and protonated water clusters[END_REF] density analysis of POSS-H + @Wn clusters reveals that Zundel-type feature is observed when an even number of water molecules are arranged in cyclic clusters in conjunction with the formation of the symmetric configuration. However, the unequal/odd amount of hydration patterns are more favorable for the transfer of proton from Zundel to Eigen. Then it forms a Zundel core again for even number of water in the hydration process. This is in line with earlier reports, which revealed that proton migration in linear type clusters is more favorable than in the branched and cyclic motifs. [START_REF] Prakash | Stepwise hydration of protonated carbonic acid: a theoretical study[END_REF][START_REF] Prakash | Structure, stability and spectral signatures of monoprotic carborane acid-water clusters (CBWn, where n= 1-6)[END_REF][START_REF] Prakash | Ab initio and density functional theory (DFT) studies on triflic acid with water and protonated water clusters[END_REF] This is because of the formation of cyclic H-bonded water chain, of a large number of H-bonds and of the perfect localization of the Eigen core in the middle of the framework. A similar mechanism was recently found by Lowe et al. [START_REF] Lowe | Acid-base dissociation mechanisms and energetics at the silica-water interface: An activationless process[END_REF] at the silicawater interface.

Discussion

The BEs of the POSS-H + @Wn clusters agree well with the earlier reports on proton transfer models such as PCA, CA and TA with water clusters (where n = 1-6). [START_REF] Prakash | Stepwise hydration of protonated carbonic acid: a theoretical study[END_REF][START_REF] Prakash | Structure, stability and spectral signatures of monoprotic carborane acid-water clusters (CBWn, where n= 1-6)[END_REF][START_REF] Prakash | Ab initio and density functional theory (DFT) studies on triflic acid with water and protonated water clusters[END_REF] The BEs depend on the proton leaving ability of the parent molecule. They will reflect the number of required water molecules to dissociate the proton from the acid. At the B3LYP/6-311++G** level, the stabilization of the Eigen cation occurs at hexameric water clusters in linear model is computed ~ -119, -124, -196 and -118 (in kcal/mol) for POSS, PCA, CA and TA, respectively. For the branched types lower energies are derived: POSS (~ -117 kcal/mol), PCA (~ -112 kcal/mol) and TA (~ -117 kcal/mol). Note that the BEs of POSS-H + @W6 clusters are in close agreement with the Nafion based triflic acid-water cluster model. [START_REF] Prakash | Ab initio and density functional theory (DFT) studies on triflic acid with water and protonated water clusters[END_REF] This confirms that the proton transfer ability of POSS moiety at solid-liquid interface is online with well-known PEMFCs such as Nafion. The main advantage of the POSS framework is that it is insoluble in water, chemicals and it presents a thermal stability even at elevated temperatures. The POSS moiety acts as H + source/carrier and at the same time water molecules act as transporter for the PEMFC materials.

Conclusion

This work presents the gas-phase geometries, binding energies, and infra-red spectral signatures of protonated POSS -water (POSS-H + @Wn, where n = 1-6) clusters as calculated using DFT based B3LYP/6-311++G** method. It is the first report on the prediction of Zundel features in POSS-H + @Wn clusters interface and provide valuable information to detect the Zundel ion in such clusters. We reveal also the existence of Zundel  Eigen  Zundel features that are trapped at the POSS-water interface. Moreover, close scrutiny of linear chain isomers shows that the proton transfer occurs from one water molecule to another molecule via Grotthuss-type mechanism and these clusters are stabilized by the Eigen core with the Wn clusters at the middle of the chain. It is interesting to note that the formation of the Eigen cation requires a minimum of three water molecules, whereas that of the Zundel ion needs four water molecules at the POSS-water interface. At the same time, both POSS moiety and water molecules act as an H + source and better transporter for the PEMFC materials, respectively. These features can enhance proton transport through the water chain at the microscopic hydration level. Such water molecules can be qualified as active water molecule, since it can accept a proton from POSS with simultaneous transfer of a proton to a neighboring water molecules. Since such mechanism was already found at the silica-water interface [START_REF] Lowe | Acid-base dissociation mechanisms and energetics at the silica-water interface: An activationless process[END_REF] and since the present work hints at its occurrence at nano-water clusters, it may be more general than usually thought, and we suggest that it should be investigated further, both by computational and experimental means, to improve our understanding of proton transport in water near interfaces.

Figure 1 :

 1 Figure 1: B3LYP/6-311++G** optimized geometries of protonated polyhedral oligomeric silsesquioxane water clusters (POSS-H + @Wn; where n = 1 -6). Distances are in Å. The dashed blue circle and rectangle highlights the formation of H3O + and H5O2 + cations, respectively.

Figure 2 :

 2 Figure 2: B3LYP/6-311++G** calculated gas phase IR spectrum of the POSS-H + @Wn (n=1-4) clusters. The three component bands (A) (black), (B) (blue) and (C) (brown) correspond to free water O-H stretching, to Eigen core to water stretching and to Zundel core proton oscillation, respectively.

Figure 3 :

 3 Figure 3: O*-H + distance as a function of time during a 100 K AIMD simulation, where O* is the oxygen atom of the POSS originally bonded to the excess proton in POSS-H + .

Figure 4 :

 4 Figure 4: Final structure after reactivity of the POSS-H + system with water, as observed in AIMD (see text for details); 52 water molecules not directly involved in the chemical event are not displayed, for the sake of clarity.

Figure 5 :

 5 Figure 5: Reaction pathway and respective energy barrier (in kcal/mol) for the different conformation of the Eigen cation at POSS-W3 interface (at PBE+D3/TZVP method). (a), (b) and (c) correspond to snapshots at the beginning, at the maximum of the energy barrier and at the end of the reaction, respectively. The blue dotted circle indicates the proton localized sites.

Figure 6 :

 6 Figure 6: Proposed mechanism for the proton transfer from the POSS framework in to water clusters. The excess proton in the POSS-H + @Wn clusters is shown in blue color.

  with water clusters (Wn, where n = 1-6). Moreover, the formation of Zundel ion requires a symmetrical arrangement of water molecules, whereas asymmetric hydration patterns dislocate into Eigen ion. The changes in bond distances upon attaching water molecules to POSS-H + suggest that the excess proton can jump from one molecule to other through H-bonded water chain during the hydration. Since the proton localized in the center of the water chain, the charge dissociation / association mechanism of water molecule plays a vital role in the formation of hydronium ion and the same interacts with the neighbor. This clearly reveals that the primary shell interaction of water molecule strength reduces during the migration of a proton.The energetics of the clusters allow determining the minimal energy pathway to conduct proton through water chain via the identification of the most stable isomers, important for the proton hopping mechanism. The origin of the proton hopping mechanism is the existence of Eigen and Zundel cationic defects. The excess proton localized in the middle with the formation of Zundel feature illustrates the origin of proton conduction via the Zundel to Eigen transformation. When slight distortion happens in the molecular arrangement, H-bonds can favor the formation of Eigen or other short-lived species from symmetrical Zundel core. The transformation of Zundel-Eigen-Zundel is favored by increasing the number of water molecules within these clusters. Close scrutiny of the geometries, IR spectra and AIM electron

Table 1 :

 1 BSSE corrected binding energies (BE in kcal/mol) with and without considering the dispersion correction (+D3) for the POSS-H + @Wn (where n=1-6) clusters using the B3LYP/6-311++G** method.

	NHB corresponds to the number of Hydrogen bonds.
	POSS-H + @Wn Ionic defect in water clusters	BE BSSE BSSE+D3	NHB	BE/NHB BSSE+D3
	POSS-H + @W1	H-Bond with W1	-21.4	-23.1		-23.1
	POSS-H + @W2	H-Bond in W2 chain	-40.8	-43.7		-21.8
	POSS-H + @W3a	Eigen in W3 chain	-84.7	-89.1		-29.7
	POSS-H + @W3b Eigen in POSS-W interface	-81.9	-86.6		-28.9
	POSS-H + @W3c	Eigen in cyclic W3	-86.0	-93.0		-23.3
	POSS-H + @W3d	H-Bonded Wn Cluster	-44.8	-51.3		-17.1
	POSS-H + @W4a	Eigen in W4 chain	-95.3	-100.2		-25.1
	POSS-H + @W4b Eigen in POSS-W interface	-94.7	-100.3		-25.1
	POSS-H + @W4c	Formation of Zundel	-58.7	-65.3		-16.3
	POSS-H + @W4d	Formation of Zundel	-58.0	-64.5		-16.1
	POSS-H + @W5a	Eigen in W5 chain	-111.6	-117.3		-23.5
	POSS-H + @W5b	Eigen cyclic W5	-109.7	-117.8		-19.6
	POSS-H + @W5c POSS-H + @W5d	Eigen in POSS-W interface	-107.2 -106.5	-119.2 -112.9		-23.8 -22.6
	POSS-H + @W6a	Eigen in W6 chain	-119.1	-125.5		-20.9
	POSS-H + @W6b POSS-H + @W6c	Eigen in POSS-W interface	-119.2 -117.5	126.3 -124.8		-21.1 -20.8
	POSS-H + @W6d Zundel formation in W6 Ring -81.3	-84.0		-14.0
	POSS-H + @W6e	Eigen in W6 Book	-119.0	-129.3		-16.2
	POSS-H + @W6f	Eigen in W6 Ring	-120.8	-129.1		-18.4
	Table

  Table1shows indeed that there is a reduction of |BE| for the corresponding clusters. For illustration, Table1gives the BE/NHB quantity which corresponds to the BE per H-bond within the complex. They are in the range -15 to -30 kcal/mol with minimal values for clusters having Zundel ion structure or H-bonded Wn cluster. We further note that both clusters having Zundel motifs possess BE/NHB of ~ -16.0 kcal/mol which is close to the value for the Zundel ion (H + (H2O)2)

Acknowledgment K. R. M. and S. K. thank SRM Institute of Science and Technology (SRM-IST) Research Fellowship for his research work. M. P. thanks the Department of Science and Technology-Science and Engineering Research Board (DST-SERB) of India for the financial support (Grant number: 17 ECR/2017/000891). The authors also thank SRM Supercomputer Centre (HPCC), SRM Institute of Science and Technology for providing the computational facility and financial support, and acknowledge access to HPC platforms provided by a GENCI grant (A0090807069).

Appendix A. Supplementary data

Supporting Information: This file contains the AIM Analysis and the vibrational Analysis of POSS-H + @Wn (where n=1 to 6) and the corresponding tables and figures. Reaction pathway and energy barrier analysis of POSS-H + @W4 clusters using NEB method.


Global climate models like any in sillico numerical experiments are affected by different types of bias. Uncertainties quantification remains a challenge in any climate attribution analysis. A fundamental methodological question is to determine which statistical summaries, while bringing relevant signals, can be robust with respect to multi-model errors. In this paper, we propose a simple statistical framework that significantly improves signal detection in climate attribution studies. Under a well-defined and verifiable assumption, the complex bias correction step can be entirely bypassed. To illustrate our approach, we infer emergence times in precipitation from the CMIP5 and CMIP6 archives. The anthropogenic signal in yearly maxima of daily precipitation clearly emerges at the beginning of the 21st century. In addition, no CMIP model seems to outperform the others and a convex combination of all improves the estimation of emergence times.

Significance statement. We show that the bias in multi-model global climate simulations can be efficiently handled when the appropriate metric is chosen. This metric leads to an easy-toimplement statistical procedure based on a checkable assumption. This allows us to demonstrate that optimal convex combinations of CMIP outputs can improve the signal strength in finding emergence times. Our data analysis procedure is applied to yearly maximum of precipitation from CMIP5 and CMIP6 databases. The attribution of the anthropogenic forcing clearly emerges in extreme precipitation at the beginning of the 21st century.

Introduction

Global climate model outputs like any numerical simulations correspond to an approximation of the true system under study, here the climate system. In the realm of Detection and Attribution (D&A), either in a transient setup or in the context of extreme event attribution (EEA), numerous review studies (see, e.g. [START_REF] Chen | Recent progress and emerging topics on weather and climate extremes since the fifth assessment report of the intergovernmental panel on climate change[END_REF][START_REF] Stott | Attribution of extreme weather and climate-related events[END_REF][START_REF] Shepherd | A common framework for approaches to extreme event attribution[END_REF]) listed different sources of variability, uncertainties and errors. In particular, these reviews highlighted that model error in numerical experiments like the Coupled Model Intercomparison Project (CMIP) can be large and has to be taken into account in any D&A statistical analysis (see, e.g. [START_REF] Knutti | Uncertainty quantification using multiple models -prospects and challenges[END_REF]National Academies of Sciences and Medicine 2016).

To address the issue of multi-model error in attribution studies, we need to go back to the origin of D&A. This research field aims at answering questions related to relative changes between two worlds. In EEA, a factual scenario of conditions that occurred around the time of a specific event is compared to the probability of the same event but under a counterfactual scenario in which anthropogenic emissions had never occurred (see, e.g. [START_REF] Angélil | An independent assessment of anthropogenic attribution statements for recent extreme temperature and rainfall events[END_REF]). In D&A with transient runs, the two worlds correspond to global coupled climate runs with all forcings (ALL) and with only natural forcings (NAT), respectively (see, e.g. [START_REF] Hegerl | Use of models in detection and attribution of climate change[END_REF]. To combine model error uncertainties, various authors (see, e.g. [START_REF] Lorenz | Prospects and caveats of weighting climate models for summer maximum temperature projections over north america[END_REF] noticed that giving equal weight to each available model projection may be suboptimal. In addition, model inter-dependencies have been identified as an important issue in uncertainties analysis (see, e.g. [START_REF] Abramowitz | Esd reviews: Model dependence in multi-model climate ensembles: weighting, sub-selection and out-of-sample testing[END_REF]. To integrate multi-model error into the EEA or transient D&A context, we leverage a hypothesis from the bias-correction community to propose an easy-to-implement strategy that, under well identified conditions, has the main advantage of bypassing multi-model error. In addition, the problem of model inter-dependencies is handled by fixing a robust referential invariant for this issue. Our main application is the inference of emergence times in yearly maxima of daily precipitation and temperatures from the CMIP database.

Methods

a. Notations and assumptions

A common hypothesis in most D&A studies is that, although numerical models may not be able to exactly reproduce the true world, one can expect that any appropriate bias correction technique for a given numerical model should be applied in the same way to correct the factual and counterfactual worlds. For example, if factual temperature runs from a given model are too warm with respect to recorded measurements during a specific period and have to be corrected, say by one Kelvin in the factual world, then this bias of one Kelvin has to be corrected in the counterfactual world of this model during the same specific epoch. In practice, such an hypothesis can be challenged and needs to be assessed with caution (see, e.g. [START_REF] Maraun | Towards process-informed bias correction of climate change simulations[END_REF]. Due to the lack of records in the counterfactual world, the assessment cannot be done without carefully designing the appropriate metric. To explain this point, we need a few mathematical notations. Let Z t and X t be the same real-valued continuous random variable of interest for the year t but from the hypothetical true factual and counterfactual worlds with cumulative distribution functions (cdf) F t (z) = P(Z t ≤ z) and G t (x) = P(X t ≤ x), respectively. In simulation studies, we never have accessed to these two random variables because perfect factual and counterfactual distributions cannot be reproduced exactly. Instead, imperfect ensemble outputs, say from M different numerical model experiments are available, and we denote this by Z (m) t and X (m) t the factual and counterfactual versions from model m and cdfs

F (m) t (z) = P(Z (m) t ≤ z) and G (m) t (x) = P(X (m) t ≤ x), respectively.
A popular approach in the bias correction literature for univariate continuous random variables is the quantile mapping transform that matches two random variables with different distributions (see, e.g. [START_REF] Maraun | Towards process-informed bias correction of climate change simulations[END_REF][START_REF] Cannon | Multivariate quantile mapping bias correction: an n-dimensional probability density function transform for climate model simulations of multiple variables[END_REF]. A positive aspect of quantile mapping is its theoretical basis. It can be viewed as the solution of an optimal transport problem (see, e.g. [START_REF] Robin | Multivariate stochastic bias corrections with optimal transport[END_REF] and it has been tested in various settings (see, e.g. for precipitation downscaling, [START_REF] Kallache | Nonstationary probabilistic downscaling of extreme precipitation[END_REF]. Quantile mapping can be adapted to our D&A framework in the following way. As any real-valued random variable, the climate model output X (m) t can always be transformed, in a distributional sense, into the unobserved true counterfactual random variable X as follows

X t d = G ← t • G (m) t X (m) t
, where d = corresponds to a equality in distribution and G ← t (.) represents the inverse of G t , i.e. its quantile function. The same type of operator can be implemented in the factual world, i.e.

Z t d = F ← t • F (m) t Z (m) t .
Then, the main assumption that we will test can be expressed in the following way. We say that model m verifies assumption (A) if

A : F ← t • F (m) t =G ← t • G (m) t , for t ∈ 1, . . .,T. (1)
(1) means that the bias between the simulated and unobserved counterfactual worlds is the same as between the simulated and unobserved factual worlds. Note that (1) allows for non-linear bias correction. For example, heavy rainfall can have different tail behaviors in the observed and model runs (see, e.g. [START_REF] Coles | An Introduction to Statistical Modeling of Extreme Values[END_REF], for an introduction to upper tails modeling). It is also important to notice that, although none of the available climate models may exactly satisfy assumption (A), it makes sense to weigh models according to their capability to satisfy (A). Another key point is that assumption (A) is related to a specific variable. For the same climate model, (A) can be valid for mean hemispheric temperatures, but incorrect for heavy rainfall over a specific region.

Another feature of (A) is the temporal indexing. The subscript t makes the notation complex, but it allows for having different bias corrections for different years, and brings flexibility. The cooling effect of volcanic forcing, like Pinatubo in the early nineties, can be included in our bias correction approach. The same could be said for slow changes due to solar forcing. So, the hypothesis of temporal stationarity is not needed in our framework.

b. Quantities of interest

In D&A studies, most researchers (see, the bibliography in the review articles [START_REF] Stott | Attribution of extreme weather and climate-related events[END_REF][START_REF] Naveau | Statistical methods for extreme event attribution in climate science[END_REF] aim to contrast differences between two worlds, e.g. the factual and counterfactual worlds in EEA, or NAT and ALL worlds in long transient climate simulations analysis. In particular, making the distinction between the two survival functions

P(X t > u) and P(Z t > u) (2)
has been a recurrent theme in this field. To explain our statistical approach, we can start by asking a typical hydrological D&A question. For the current year t, are precipitation intensities in the factual world heavier than the ones produced in the counterfactual one? If two random precipitation intensities would have the same distribution in the factual and counterfactual world, then the probability to observe the event {Z t > X t } would be as likely as {Z t < X t } for any given year, and consequently P(Z t > X t ) = 0.5 in this case1. If instead, one answers positively to the question, this will imply that the chance of Z t being greater than X t is greater than 0.5. To make the link with the classical (2), we can look at the special case where the threshold u in (2) is chosen to be equal to a random draw from G t . This choice leads to our definition of two simple probabilities q 0 = 1 2 and q t = P(Z t > X t ).

(3)

These two probabilities have many advantages. They are invariant with respect to non-decreasing changes. For example, if both X and Z are simultaneously multiplied by two, then q t remains the same. This is a critical feature when climate models are bias corrected. More precisely, under assumption (A), we always have, see SI for a proof,

q t = P Z (m) t > X (m) t
, for all t = 1, . . .,T .

The fact that the left-hand part of this equation does not depend on m anymore is fundamental in this work. The remaining part of this article is to explain its consequences, its applicability and its validity within the CMIP database.

Under (A), we do not need to have access to X t and Z t to compute q t . This probability can be directly obtained from numerical outputs encapsulated by X (m) t and Z (m) t . The most important consequence of (A) is that biased models outputs like (X (m) t , Z (m) t ) T do not need to be corrected.

Practically, this also implies that we do not need to know the cdf's, F t and G t , to estimate q t . In addition, q t in (3) is always equal to .5 whenever the true factual and counterfactual worlds are exchangeable. Consequently, no inference is required in this case. Last but not least, the value of .5 can be used a reference point when the factual and counterfactual start to differ.

1The independence between X t and Z t is not necessary to get P(Z t > X t ) = 0.5. As the sum P(Z t > X t ) + P(Z t < X t ) is always equal to one, the only requirement is that P(Z t > X t ) = P(Z t < X t ). This is always true under exchangeability, i.e. whenever(Z t , X t ) d = (Z t , X t ). For example, a standardized and correlated bivariate Gaussian vector is exchangeable.

c. Inference in the transient case

The temporal indexing t in the probability q t defined by (2) can be interpreted in two different ways. One can freeze the time to the current year and this leads to the so-called event attribution realm. Large ensembles of this given year are classically drawn from both factual and counterfactual worlds (see, e.g. [START_REF] Stott | Attribution of extreme weather and climate-related events[END_REF] In this work, the 16 model runs used in this work are listed in Table 12. Their main drawbacks are the model uncertainties, their small ensemble sample sizes and their spatial resolution, which can be too coarse for some applications. Statistically, a subtle point is the transient nature of these simulations. This implies that factual runs in CMIP contain some non-linear trends that should be taken into account in the statistical analysis (see, e.g. [START_REF] Kharin | Changes in the extremes in an ensemble of transient climate simulations with a coupled atmosphere-ocean GCM[END_REF]. As the ALL run Z (m) t distribution may change over time, we estimate the time-varying q (m)

t = P Z (m) t > X (m)
t in a non-parametric regression manner. A classical kernel regression approach [START_REF] Nadaraya | On estimating regression[END_REF][START_REF] Watson | Smooth regression analysis[END_REF]) leads to the following estimator q(m)

t = 1 K h (t -t j ) J j=1 K h (t -t j ) G (m) I Z (m) t j , (5)
2After the year 2006, we analyze CMIP precipitation under their respective worst-case scenarios (RCP8.5 and SSP5-8.5). Our approach can be easily implemented under different scenarios. Concerning emergence times, they occur during the first two decades of the 21st century.

Consequently, this result remains valid under other scenarios.

where the positive function K h (.) corresponds to a weighting Kernel with bandwidth h, and G (m) I (.)

represents any estimator of G (m) t (.). In the SI, the choices of the kernel and G (m) I (.) are discussed and the statistical arguments to build asymptotic confidence intervals for q t are given.

As an example, the inferred probability q(m) t for yearly maxima of CMIP daily precipitation for the grid point near Oxford in Great Britain is displayed in Figure 1. Each panel represents one CMIP climate model from Table 1. In each panel, the x-axis spans the year t = 1850 to t = 2100 and the y-axis corresponds to the probability q(m) t obtained from ( 5). The departure from the horizontal line at q t = .5 indicates a change between factual runs and counterfactual trajectories without anthropogenic forcing. All panels show a smooth increase in q(m) At this stage, assumption (A) has not been used yet. If (A) was satisfied, it would be straightforward to move, via (4), from q(m) t to a common estimator of q t . To combine climate model outputs, the capability to satisfy (A) for a given model m will be key.

Merging climate model simulations

The main roadblock to asses the quality of simulated runs is that we will never observe draws from X t , only measurements with observational errors at spatial scales different from climate gridded data. This lack of data is even worse from Z t , the hypothetical world of an unperturbed and never observed climate. To bypass this difficulty, we assume that there exists a time period, say T , during which the ALL and NAT worlds were identical (in distribution). In practice, this corresponds to the pre-industrial period for which we assume

F t = G t , for all years t ∈ T . (6)
Again, we do not need to assume that F t = F t+1 as internal climate variability and natural forcings may change in time, even at the annual scale, e.g. like after the 1883 Krakatoa eruption. In our Application section, we define the pre-industrial epoch T as T = {1850, . . ., 1900}. Then, we can always write

q t = P(Z t > X t ) = 1 2 , for any year t ∈ T . ( 7
)
This equality is parameter-free and it does not depend on model m. There is nothing to estimate, so no inferential error needs to be taken into account. As already highlighted in Section b, the bivariate vector (X t , Z t ) does not have to be stationnary in t, only exchangeability between X t and Z t . This last condition is always satisfied as X t and Z t are not computer simulated. They just represent conceptual independent draws from a thought experiment of two possible climate trajectories in pre-industrial times. For simulated runs from model m, we also expect to have

P(Z (m) t > X (m) t ) = .5 for any year t ∈ T if Z (m)
t and X (m) t are exchangeable (label free) during the pre-industrial period. Exchangeability is a weak hypothesis with respect to the issue of model interdependencies studied by [START_REF] Abramowitz | Esd reviews: Model dependence in multi-model climate ensembles: weighting, sub-selection and out-of-sample testing[END_REF]. Concerning the validity of

P(Z (m) t > X (m) t ) = .5
in regard to a given model, we leverage the following fact. During the time period T , we have F t = G t , and consequently, the following equivalence is always true ranges strongly vary among models, e.g. compare the precipitation spread between IPSL-CM5A-LR and CERFACS-CNRM-CM6. By construction, this difference is not an issue because we always look at relative changes within a CMIP run. According to (7), we also expect that q (m) t to be close to half for t ∈ T . This can be used to weight and merge our 16 models.

F (m) t =G (m) t , ∀t ∈ T ⇐⇒ (A) holds for any t ∈ T . Hence, checking that F (m) t is equal to G (m) t for
Binary events like {Z (m) t > X (m) t } during the preindustrial period are the building blocks of

q (m) t = P Z (m) t > X (m) t .
In our setup, the reference density is the Bernoulli sequence under the success probability q t = .5 and the competitor corresponds to a Bernoulli sequence with the success probability q (m) t . A statistical tool is needed to differentiate these two Bernoulli distributions. The Kullback Leibler divergence (see, e.g. [START_REF] Burnham | Model Selection and Inference: a Practical Information-theoretical Approach[END_REF][START_REF] Naveau | A non-parametric entropy based approach to detect changes in climate extremes[END_REF]) compares two distributions by calculating the expectation of the logarithmic difference between a target probability and a competitor, where the expectation is obtained with respect to the target density; see the SI method section for our Bernoulli and (e.g. [START_REF] Haughton | Weighting climate model ensembles for mean and variance estimates[END_REF][START_REF] Knutti | Challenges in Combining Projections from Multiple Climate Models[END_REF]).

As any Kullback Leibler divergence is convex, it is natural to merge our estimates as a convex combination defined in the following way

qt = M m=1 w m × q(m) t , with w i ≥ 0 and w 1 + • • • + w M = 1. ( 8
)
To compute the weights (w 1 , . . ., w M ) T , we implement a two-step procedure. Small p-values of the two-sample Anderson-Darling test in the previous section, see yellow-orange diamonds in Figure 2, highlighted a poor fit. For this reason, our first step is to remove all models that have Anderson-Darling p-values smaller than 20%, a very conservative rejection rate. This first step allows to treat the rare but possible case when all models are wrong and to give give them null weights. Our second step is to simply find the weights of selected models that minimize the Kullback Leibler divergence between c(.5, . . ., .5) and qt∈T under the constraint

w 1 + • • • + w M = 1.
Concerning CMIP yearly maxima of precipitation around Oxford, the 90% green confidence band in Figure 3 represents the estimate of the convex weighted combination, i.e. of qt with the weights shown in Figure 1. In Figure 3 , the 90% red confidence band corresponds to the model with the highest Anderson-Darling p-value. As expected, combining estimates of q t reduces the confidence bandwidth. By construction, qt follows well the referential horizontal line centered at .5 during the pre-industrial period. The departure from this horizontal grey line is significant around the year 2000. Overall, the estimate qt has a smooth trajectory over time and the detected signal is strong in 2100. A clear indication that anthropogenic forcing can cause changes in precipitation intensities at this location.

Emergence times in yearly precipitation maxima

Our methodology can be applied to any type of real-valued continuous atmospheric variables:

temperatures, wind speeds, precipitation intensities and others. As many D&A studies have already focused on temperatures, we chose to study annual maxima of daily precipitation (results about temperatures are available upon request). Statistically, precipitation intensities are skewed and sometimes heavy tailed. So, this type of non-gaussian random variable represents a challenging testbed for our approach. In addition, multi-error analysis with annual precipitation maxima can reveal key information for impact studies. Details about our selected CMIP models can be found in the SI data section.

To illustrate our statistical method in the previous section, we focused on the Oxford grid point in figures 1, 2 and 3. The same type of analysis can be done for each individual grid point. The nine panels of Figure 4 provide global snapshots of the estimates qt for nine different years: 1850, 1900, 1940, 1970, 2000, 2020, 2030, 2050 and 2100. As already pointed out in Figure 3 Figure A1 shows that no specific model appears to outperform, in terms of weights, the other models. This complements our understanding at the Oxford grid point in Figure 2 where four among the 16 models had strong weights. For other grid point locations, other models are chosen.

Overall, if all models were equiprobable, then, in average, each model should have a weight around 1/16 = 0.0625. This value basically corresponds to the the fourth column of Table 1 indicates that the global estimate average weights for each model. Concerning the first step of our weight procedure during which we only kept models that have a Anderson-Darling p-value above 0.2 at a given grid point, the last column of Table 1 confirms that, at the global scale, the average number of grid point rejected by each model is, as expected, around 20%. This points toward the fact that no model appears to be superior (less rejected) at the global scale.

Emergence times definition

As already mentioned, a consequence of our definition q t is that, when there is no difference between the factual (ALL) and counterfactual (NAT) worlds, it has to be always equal to .5. This robust referential allows us to define a emergence time in the following way τ p = argmin{t for all t ≥ t, we have P ( qt > 0.5) > p}.

This means that, at the confidence level p, all years after the emergence time τ p have a q t significantly higher than 0.5. So, an increase in precipitation is detected at year τ p and this signal remains present after this specific year. From Figure 5, one can deduce that, in most regions, the anthropogenic signal becomes detectable in precipitation around the year 2000 at the confidence level of 90%. In northern latitudes (below 50N), the detection starts to emerge even as early 1950's. The grey areas do not mean that the signal is not attributable, but, instead of an increase, these regions correspond to a precipitation decrease, see the green areas in the bottom three panels of Figure 4 and (e.g., [START_REF] Pfahl | Understanding the regional pattern of projected future changes in extreme precipitation[END_REF][START_REF] Tandon | Understanding the dynamics of future changes in extreme precipitation intensity[END_REF][START_REF] Dong | Attribution of extreme precipitation with updated observations and cmip6 simulations[END_REF].

Conclusions and discussions

From a climatological point of view, our analysis clearly indicates changes in precipitation intensities. Emergence times and precipitation patterns are spatially coherent, see figures 4 and 5.

In 2100, precipitation in most regions of the globe appears, under scenario RCP8.5, to be either strongly or unequivocally affected by the prescribed anthropogenic forcing. We also found that a convex combination of all models perform better than any single one.

The optimization of model weights is tuned during the pre-industrial period. This calibration is assumed to be valid over other time periods and one can wonder if our results are robust with respect to weights changes. To explore this possibility, we also studied two other ways of setting weights. We implemented the classical form of weights in variable selection problems (see, e.g. [START_REF] Burnham | Model Selection and Inference: a Practical Information-theoretical Approach[END_REF] and used in the climate literature (see, e.g. [START_REF] Lorenz | Prospects and caveats of weighting climate models for summer maximum temperature projections over north america[END_REF].

Another approach is a so-called expert aggregation approach (e.g., [START_REF] Gaillard | Forecasting electricity consumption by aggregating experts; how to design a good set of experts[END_REF].

Overall, all three approaches gave similar emergence times, with well structured spatial patterns with the expert aggregation technique and our approach. For these two methods, estimated of qt are robust to weight changes during the pre-industrial period. For other time periods, the absence of counterfactual perfect observations makes it impossible to remove of the assumption of stationary weights in a D&A context. Other research avenues could be explored to integrate observational data in our statistical approach. A delicate issue, especially for precipitation, is to define the spatial scale of interest. In most CMIP based D&A studies like in [START_REF] Ribes | Making climate projections conditional on historical observations[END_REF], the analysis is done at global or regional scales. This facilitates the integration of observations and improve the signal-noise ratio. But, gridded spatial features like the ones in Figure 4 are lost. This also leads to the open question of how to find optimal regions that maximize the attribution power in a multivariate context (see, e.g. [START_REF] Kiriliouk | Climate extreme event attribution using multivariate peaksover-thresholds modeling and counterfactual theory[END_REF]. A related question is how to adapt our approach to a multivariate framework to attribute compounds events (see, e.g. [START_REF] Zscheischler | The effect of univariate bias adjustment on multivariate hazard estimates[END_REF]. Coupling the field of counterfactual theory and multivariate analysis could help in the direction (see, e.g. [START_REF] Hannart | Counterfactual causality theory for the attribution of weather and climate-related events[END_REF]. A last point is our use of the CMIP worst-case scenarios (RCP8.5 and SSP5-8.5). Our proposed technique can be easily applied to other scenarios. Still, most inferred emergence times of yearly maxima of precipitation span the period 2000-2020, see the blue to yellow regions in Figure 5. This period is prior to any strong differences among scenarios. Hence, our estimated emergence times for these regions will remain valid with other scenarios. From the CMIP5 archive, we select 12 models for which we found a complete set of precipitation simulations for the historical (1850-2005), historicalNat (1850-2012) and RCP8.5 (2006[START_REF] Zscheischler | The effect of univariate bias adjustment on multivariate hazard estimates[END_REF] experiments. Additionally, we also treat four CMIP6 models (see, e.g. [START_REF] Eyring | Overview of the coupled model intercomparison project phase 6 (CMIP6) experimental design and organization[END_REF]) for which we found historical simulations in the CMIP6 deck, hist-nat simulations in DAMIP and SSP585 projection. The historical simulations combined with the RCP8.5 and SSP585 simulations represent the factual world whereas the historicalNat simulations correspond to the counter-factual world. All runs have been remapped to a common 5 o × 5 o HadCRUT grid (cdo rmapcon operator).

Acknowledgments

Table 1 below provides the list of the 16 CMIP models used in this study.

b. Methods

Proof of Eq. ( 4): To simplify the proofs, we have drop the temporal indexing t whenever it was possible. We have by definition of Z t and X t ,

q t = P (Z t > X t ), = P F ← t • F (m) t (Z (m) t ) > G ← t • G (m) t (X (m) t ) , = P G ← t • G (m) t (Z (m) t ) > G ← t • G (m) t (X (m) t ) , from (A), = P G (m) t (Z (m) t ) > G (m) t (X (m) t ) , as G t (.) non-decreasing, = P Z (m) t > X (m) t , as (G (m) t ) ← (.) non-decreasing.

Computation of confidence intervals of q(m)

t

For each climate model m, the year-varying probability q (m)

t = P Z (m) t > X (m) t = E G (m) t (Z (m) t )
can be inferred from the Nadaraya-Watson estimate q(m) t defined by Eq. ( 5). To simplify the inference process, we consider the cdf G t of the NAT as a stationary process3. In this case, G t can be estimated by the classical empirical estimator G (m) I defined as

G (m) I (x) = 1 I I I=1 I(X (m) t i ≤ x),
where I(A) represents the indicator function equal to one if A is true and zero otherwise and (X t 1 , . . ., X t I ) T corresponds to a NAT run trajectory of I time steps. To derive confidence intervals around this estimate, we need to introduce q(m)

t = 1 K h (t -t j ) J j=1 K h (t -t j ) G (m) I Z (m) t j
that corresponds to a simpler version of q(m) t in which the cdf G (m) is supposed to be known. In this case, the asymptotic behavior of the difference A = q(m) tq (m) t is well known (see, e.g. [START_REF] Härdle | Smoothing Techniques: With Implementation in S. Springer Series in Statistics[END_REF]. It converges towards a Gaussian limit law with known asymptotic mean and variances, (nh) 1/2 q(t j )q(t j )

(σ 2 (t j ) K 2 2 / f (t j )) 1/2 k j=1 converges into toward N [q (t j ) + 2q (t j ) f (t j )/ f (t j )] ∫ s 2 K(s)ds k j=1 , I
where σ 2 (t) is the conditional variance of {G t (Z t )} and f is the density of the variable T. Hence, we can write q(m) tq (m)

t = A + B, with B = 1 K h (t -t j ) J j=1 K h (t -t j ) G (m) I Z (m) t j -G (m) Z (m) t j .
3This assumption could easily remove by adding a kernel estimator for G (m) t . But, this was not necessary for our precipitation application Following [START_REF] Naveau | Revising return periods for record events in a climate event attribution context[END_REF], we consider that

√ I G(m) (Z (m) t j ) -G (m) (Z (m) t j ) behaves (almost) like a Brownian Bridge B = 1 √ I K h (t -t j ) J j=1 K h (t -t j ) √ I EG (m) I Z (m) t j -G (m) Z (m) t j ≈ 1 √ I K h (t -t j ) J j=1 K h (t -t j ) √ I B(G (m) (Z (m) t j )),
where B(u) represents a classical Brownian Bridge on [0, 1] with E[B(u)] = 0 and covariance:

Cov(B(u), B(v)) = min(u, v) -uv. Thus, E[B] = E[ q(m) t -q(m) t ] ≈ 0E and Var[B] ≈ 1 I( K h (t -t j )) 2 J j=1 J i=1 K h (t -t j )K h (t -t i ) × E min G (m) (Z (m) t i ), G (m) (Z (m) t j ) -G (m) (Z (m) t i )G (m) (Z (m) t j ) .
Assuming that the terms A and B are independent and that their mean is negligible, we build confidence intervals for q (m) t assuming that asymptotically:

q(m) t ∼ N q (m) t , Var(A) + Var(B) .
To go one step further, we need to take into account that q (m) t may not be centered around the real quantity of interest q t . It is centered if assumption (A) is satisfied, otherwise a bias exists between q (m) t and q t . To deal with this issue, we additionally assume that such a bias is randomly distributed between all our climate models in the following way

q (m) t ∼ N q t , σ 2 t q(m) t |q (m) t ∼ N q (m) t , (σ (m) t ) 2
where q (m) t and (σ (m) t ) 2 = Var(A) + Var(B) from the previous paragraph. The variance σ 2 t captures the inter-model variability between the different climate runs. By marginalising over q (m) t , we then obtain by Bayesian conjugation for normal laws that q(m)

t ∼ N q t , σ 2 t + (σ (m) t ) 2 . (A1)
This variance decomposition simply means that the variance of q(m) t can be divided into the intra-model variance (σ (m) t ) 2 and the inter-model variance σ 2 t which can be estimated as:

σ2 t = 1 M M m=1 q(m) t -qt 2 , where qt = M m=1 w m × q(m) t .
The choice of the weights w m depends on the climate model capability to satisfy during preindustrial period. Given the weights, the weighted estimator of q t follows qt

∼ N q t , M m=1 w 2 m σ 2 t + (σ (m) t ) 2 . ( A2)
(A1) and (A2) were used to obtain the confidence intervals displayed in our figures. The kernel used in q(m) t , see Eq. ( 5), is the classical Epanechnikov kernel [START_REF] Epanechnikov | Non-Parametric Estimation of a Multivariate Probability Density[END_REF]) with a bandwidth of 60.5 years. The bandwidth has been determined on using a leave-one-out cross-validation scheme to find out the bandwidth that minimises the root mean square error (RMSE) between the estimated q1 (t i ) and the GI (Z t i ). More precisely, the cross-validation has been performed for each model individually and then we select the median of bandwidths optimised for each model.

Note that the derivation of the confidence bands relies on the assumption of independence between the estimates of q t . Although various studies (see e.g., [START_REF] Knutti | Challenges in Combining Projections from Multiple Climate Models[END_REF][START_REF] Haughton | Weighting climate model ensembles for mean and variance estimates[END_REF] pointed out that climate models are not necessarily independent, this issue may not be too prevalent in our case for the following reasons. We do not require independence between raw atmospheric variables, but between events like

{Z (i) t > X (i) t } and {Z ( j) t > X ( j)
t } for climate models i and j. These events are based on increments like

Z (i) t -X (i) t and Z ( j) t -X ( j)
t and consequently, by removing additive error, increments are more likely independent than raw data. Concerning the later, our main focus are gridded annual maxima of precipitation. Such variables have high variability and inter-model dependence of increments is secondary compared to the signal/noise ratio issue. In addition, Figure 1 indicates that the differences within a research center, see e.g.

IPSL, appear to be as important than the ones between different research laboratories.

c. Appendix Figure

The weights used in the aggregated estimator of q t defined by Eq. ( 8) are displayed in Figure A1.

The x-axis corresponds to the labels of the 16 CMIP models used in the aggregation, see Table 1.

The y-axis represents grid points locations of each model. The black (yellow) color corresponds a weight near one (zero) in Eq. ( 8).

d. SI Kullback-Leibler divergence computation

The log-likelihood of independent4 Bernoulli sequences, B (m) t = {Z (m) t > X (m) t }, over t ∈ T can be written as

t∈T log 1 -q (m) t 1-B (m) t q (m) t B (m) t .
For two binomial distributions with respective success rates p and q, the Kullback Leibler divergence is equal to D (q; p) = q log q p + (1q) log 1q 1p .

4As we analyze yearly maxima of daily values in our application, the hypothesis of year-to-year independence is reasonable. If it was not the case, then the full times series B (m) t with t ∈ T will have to be modeled as a multivariate binary random vector with a memory component (see, e.g. [START_REF] Tuel | Skillful prediction of multidecadal variations in volcanic forcing[END_REF][START_REF] Dai | Multivariate bernoulli distribution[END_REF]). This will lead to a more complex likelihood. Albeit added complexity, the principles based on the Kullback Leibler divergence exposed in this section will remain valid.

This leads to the following Kullback Leibler divergence that compares the T-dimensional vector q (m) t∈T with q t∈T = c(.5, . . ., .5) D c(.5, . . ., .5); q (m) t∈T = -T × log 2 -

1 2 t∈T log q (m) t 1 -q (m) t .
For each model m, we can measure the departure from the term T × log 2 with T = 50 years. A KL divergence is always non-negative and equals to zero if q (m) t = c(.5, . . ., .5). The estimate of q (m) t from Eq. ( 5) can be plugged in the expression D(.; .), and consequently each model m can be evaluated with respect to condition Eq. ( 7). The optimization to find the weights in Eq. ( 8) is numerically done using the solnp function in the R package Rsolnp.

In the model worlds, we also have that

P Z (m) t > X (m) t = 1 2
for all years t during which the two variables Z (m) t and X (m) t are exchangeable. The exchangeability assumption simply means that the labelling of the numerical model type (X and Z) are uninformative, i.e. P Z (m) t > X (m)

t = P X (m) t > Z (m) t .
This implies that, although climate runs from the same laboratory may be dependent and share the same code, see e.g. [START_REF] Knutti | Uncertainty quantification using multiple models -prospects and challenges[END_REF]), they are likely to be exchangeable (labelling free) during time periods with similar forcing. Another aspect is that the bivariate vector Z (m) t , X (m) t does not have to be stationary in time to have P Z

(m) t > X (m) t = 1 2 .
In particular, during the early period where the anthropogenic forcing was weak (pre-industrial period), natural forcings were still a source of non-stationarity in our climate system, but q t = .5 for t within this pre-industrial period.

I K

One can notice that, if q (m) t is stationary over T for model m, then the divergence D q t∈T ; q (m) t∈T can be expressed as a function of a variance ratio

D q t∈T ; q (m) t∈T = - T 2 × log VB (m) t VB t ,
where VB (m) t = q (m) t (1q (m) t ) is always smaller that VB t = 1/4. Equivalently, we can write the variance ratio as a function of the divergence

VB (m) t VB t = exp -2 × D q t∈T ; q (m) t∈T /T ≤ 1. If this
variance ratio is close to one (zero), then q (m) t is close (far) from q t = .5.

LIST OF FIGURES

Fig. 1. Analysis of yearly maxima of daily precipitation from 16 CMIP climate runs at the grid point near Oxford, Great Britain. In each panel, the x-axis corresponds to the years from t = 1850 to t = 2100 (RCP8.5 or SSP5-8.5 after 2006) and the y-axis represents the estimate of q t defined by (5) for the closest grid point location near Oxford. Each panel corresponds to a specific CMIP model setup, see Table 1. The shaded area denotes 90% confidence intervals around the mean estimate of q t . The horizontal black line centered on . Oxford, Great Britain. In each panel, the x-axis corresponds to the years from t = 1850 to t = 2100 (RCP8.5 or SSP5-8.5 after 2006) and the y-axis represents the estimate of q t defined by (5) for the closest grid point location near Oxford. Each panel corresponds to a specific CMIP model setup, see Table 1. The shaded area denotes 90% confidence intervals around the mean estimate of q t . The horizontal black line centered on .5 corresponds to the null-hypothesis where the factual and counterfactual worlds are indistinguishable.The diamond colors and sizes are proportional to the weight values.

q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q qq q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqq q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q qq q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q qq q MRI MRI-ESM2-0(CM6)

t

  over time, but each model appears to give a different speed and amplitude of change. Inter-model variability appears to be large and combining model errors is necessary.

  each model m during the pre-industrial epoch appears as the appropriate step. A simple scatter-plot between the ranked X (m) t and the ranked Z (m) t should be close to a straight line. Figure 2 indicates that, overall, most CMIP climate model runs behave appropriately with respect to simulated precipitation at the Oxford grid point. P-values of the two-sample Anderson-Darling test (Pettitt 1976; Anderson and Darling 1952) are represented by diamonds. Their respective size and color indicate the quality of the fit. As expected, precipitation

  , a signal starts to emerge in 2000 and, for some areas, becomes very clear in 2020. In 2050, qt departs from the referential .5 in vast regions (red and orange), especially southern and northern, and there annual precipitation maxima can be attributed to change in anthropogenic forcing. This statement integrates larger zones in 2100. It is noteworthy that a few patches of green colors indicate that anthropogenic forcing has a reverse impact, i.e. a decrease in precipitation intensities. These spatial precipitation patterns are consistent with the results of previous studies which showed that, under continued greenhouse gas emissions, heavy precipitation magnitude is expected to increase over much of the world, except in the subtropics where robust declines are projected (e.g.,[START_REF] Pfahl | Understanding the regional pattern of projected future changes in extreme precipitation[END_REF][START_REF] Tandon | Understanding the dynamics of future changes in extreme precipitation intensity[END_REF][START_REF] Dong | Attribution of extreme precipitation with updated observations and cmip6 simulations[END_REF]).

.

  Part of this work was supported by the DAMOCLES-COST-ACTION on compound events and the French national program: FRAISE-LEFE/INSU, ANR-Melody,ANR-TRex, and european H2020 XAIDA. We also thank Dr. Tuel for fruitful discussions on global precipitation patterns and other related topics. APPENDIX a. Datasets

Fig. 2 .Fig. 3 .Fig. 4 .Fig. 5 .F . 1 .

 23451 Fig. 2. Visual check to know if, within each of the 16 models, NAT and ALL runs have similar distributions over the pre-industrial period defined as years in T = {1850, . . ., 1900} for the grid point near Oxford. The x-axis corresponds to the ranked yearly maxima of daily precipitation for the ALL run, the same quantity is displayed on the x-axis but for the NAT run. Each panel represents a CMIP model. If the points are aligned along the black diagonal, then NAT and ALL distributions can be considered similar. The diamond colors and sizes are proportional to the p-values of the two-sample Anderson-Darling test (Pettitt 1976; Anderson and Darling 1952). . . . . . . . . . . . . . . . . . . . . . 32 Fig. 3. Oxford grid point: comparison between the best individual model (red), with respect to the Anderson-Darling statistic, and the best convex combination (blue/green) qt defined by (8) that merges the 16 estimates. The probability q t should be equal to .5, see the horizontal black line, over the pre-industrial period when t ∈ T = {1850, . . ., 1900}. . . . . . . . 33 Fig. 4. Values of the multi-model estimates defined by (8) at each grid point for yearly maxima of daily precipitation from CMIP runs in Table1. A red (green) color indicates a significant increase (decrease) in precipitation intensities due to anthropogenic forcing. . . . . . . 34 Fig. 5. Emergence times: years after which all q t 's are significantly higher the 0.5 at the confidence level of 90%. The grey regions correspond to a detected decrease in precipitation, see green patches in Figure 4. . . . . . . . . . . . . . . . . . . . . . 35 Fig. A1. Weights of each model in Eq. (8) to obtain Figure 4. The x-axis represents the labels of our 16 selected CMIP runs and the y-axis corresponds to grid point locations over the full globe. . 36

  qq-plot(Xm, Zm) F . 2. Visual check to know if, within each of the 16 models, NAT and ALL runs have similar distributions over the pre-industrial period defined as years in T = {1850, . . ., 1900} for the grid point near Oxford. The x-axis corresponds to the ranked yearly maxima of daily precipitation for the ALL run, the same quantity is displayed on the x-axis but for the NAT run. Each panel represents a CMIP model. If the points are aligned along the black diagonal, then NAT and ALL distributions can be considered similar. The diamond colors and sizes are proportional to the p-values of the two-sample Anderson-Darling test[START_REF] Pettitt | A two-sample anderson-darling rank statistic[END_REF][START_REF] Anderson | Asymptotic theory of certain "goodness of fit" criteria based on stochastic processes[END_REF].

LIST OF TABLES

pr, yearmax, q(t) multimodel F . 4. Values of the multi-model estimates defined by (8) at each grid point for yearly maxima of daily precipitation from CMIP runs in Table 1. A red (green) color indicates a significant increase (decrease) in precipitation intensities due to anthropogenic forcing.