domain: chim phys

text:
Timothée Devergne
email: timothee.devergne@sorbonne-universite.fr

Théo Magrino

Fabio Pietrucci

A Marco Saitta

Combining Machine Learning Approaches and Accurate Ab Initio Enhanced Sampling Methods for Prebiotic Chemical Reactions in Solution

HAL is

Introduction Prebiotic chemistry and ab initio molecular dynamics

The study of chemical reactions in solution by means of ab initio molecular dynamics (AIMD) is a full branch of computational chemistry, which is of particular interest for prebiotic chemistry, the chemistry of the origins of life. [START_REF] Pérez-Villa | Prebiotic chemistry and origins of life research with atomistic computer simulations[END_REF] Many hypotheses on how bio molecules formed from elementary bricks exist. One of them is that life building blocks were brought to earth from the interstellar medium: they were either formed in comets [START_REF] Napier | The origin of life in comets[END_REF][START_REF] Bar-Nun | Shock Synthesis of Amino Acids in Simulated Primitive Environments[END_REF] or in interstellar ice via UV radiation. [START_REF] Burton | Understanding prebiotic chemistry through the analysis of extraterrestrial amino acids and nucleobases in meteorites[END_REF][START_REF] Marcellus | Aldehydes and sugars from evolved precometary ice analogs: Importance of ices in astrochemical and prebiotic evolution[END_REF] Other works show that elementary biomolecules were formed on earth in a primordial soup, hydrothermal submarine vents [START_REF] Martin | Hydrothermal vents and the origin of life[END_REF] or mineral surfaces. [START_REF] Wächtershäuser | Groundworks for an evolutionary biochemistry: The iron-sulphur world[END_REF][START_REF] Brack | Developments in Clay Science[END_REF][START_REF] Ferris | Mineral Catalysis and Prebiotic Synthesis: Montmorillonite-Catalyzed Formation of RNA[END_REF] All these possible prebiotic synthesis environments imply very different conditions to investigate. This is where AIMD comes into play: it allows one to calculate interatomic forces from the quantum-based potential energy surface (PES) of a given systems, in order to simulate the motion of atoms in organic molecules and biomolecules, while explicitly including the non-trivial presence of the chemical environment, i.e. the solvent and/or the mineral surfaces. In order to be predictive and quantitative, it needs to be used in conjunction with enhanced sampling methods such as metadynamics, [START_REF] Laio | Escaping free-energy minima[END_REF][START_REF] Laio | Metadynamics: a method to simulate rare events and reconstruct the free energy in biophysics, chemistry and material science[END_REF] umbrella sampling (US) [START_REF] Torrie | Nonphysical sampling distributions in Monte Carlo freeenergy estimation: Umbrella sampling[END_REF] or transition path sampling. [START_REF] Bolhuis | PATH SAM-PLING: Throwing Ropes Over Rough Mountain Passes, in the Dark[END_REF] This is necessary to overcome the energy barriers of reactions and thus the time limitations of computer simulations in order to explore the possible pathways in a given environment.

In our group, some successful results were obtained in the past few years using these methods in different contexts. The famous Miller-Urey experiment [START_REF] Miller | A Production of Amino Acids Under Possible Primitive Earth Conditions[END_REF] was computationally reproduced [START_REF] Saitta | Miller experiments in atomistic computer simulations[END_REF] and the mechanism of formation of glycine explained through a novel pathway, with formamide as a crucial intermediate. The decomposition pathways of amino acids in the supposed hydrothermal conditions of meteoritic parent bodies was studied, [START_REF] Pietrucci | Hydrothermal Decomposition of Amino Acids and Origins of Prebiotic Meteoritic Organic Compounds[END_REF] also to obtain clues about a possible origin of their chirality. Recently, we have thoroughly studied the classic Strecker-cyanohydrin synthesis of glycine, [START_REF] Magrino | Step by Step Strecker Amino Acid Synthesis from Ab Initio Prebiotic Chemistry[END_REF] generally considered the main synthesis pathway of glycine, the simplest amino acid, in the primordial Earth and universe. seven elementary chemical steps, and to complete the sparse existing experimental data on thermodynamics by accurate ab initio results.

To this end, as in our previous studies, we calculated the free energy surfaces and barriers for each intermediate chemical step using an AIMD-based enhanced-sampling protocol the exploration of the chemical pathways, the identification of the transition states (TS), the definition of precise RCs, and the accurate calculation of free-energy landscapes via umbrella sampling (US), which in turn requires the calculation of a large number of ab initio forces. [START_REF] Magrino | Step by Step Strecker Amino Acid Synthesis from Ab Initio Prebiotic Chemistry[END_REF][START_REF] Pietrucci | Formamide reaction network in gas phase and solution via a unified theoretical approach: Toward a reconciliation of different prebiotic scenarios[END_REF] This latter part, although necessary to ensure the statistical convergence of the FES to within 1-2 kcal/mol, is the computationally most expensive one, and thus the bottleneck significantly limiting the extent of potential applications, and the possibility to upscale the protocol to more realistic prebiotic environments, possibly including mineral surfaces. [START_REF] Brigiano | Mechanism and Free-Energy Landscape of Peptide Bond Formation at the Silica-Water Interface[END_REF] In the past few years, some attempts to overcome this obstacle in computational prebiotic chemistry were carried out exploiting density functional tight-binding approaches. [START_REF] Kroonblawd | Generating Converged Accurate Free Energy Surfaces for Chemical Reactions with a Force-Matched Semiempirical Model[END_REF] More recently, machine learning approaches to describe the local PES along chemical processes in solution have emerged as a very appealing opportunity. [START_REF] Yang | Using metadynamics to build neural network potentials for reactive events: the case of urea decomposition in water[END_REF]

Machine learning potentials with ab initio accuracy

In the past few years, a number of machine learning methods and frameworks have been developed to tackle the problem of obtaining ab initio-level PES at a reasonable computational cost. [START_REF] Bartók | Gaussian Approximation Potentials: The Accuracy of Quantum Mechanics, without the Electrons[END_REF][START_REF] Bartók | On representing chemical environments[END_REF][START_REF] Behler | Generalized Neural-Network Representation of High-Dimensional Potential-Energy Surfaces[END_REF][START_REF] Lindsey | ChIMES: A Force Matched Potential with Explicit Three-Body Interactions for Molten Carbon[END_REF][START_REF] Wang | DeePMD-kit: A deep learning package for manybody potential energy representation and molecular dynamics[END_REF][START_REF] Zhang | End-to-end Symmetry Preserving Inter-atomic Potential Energy Model for Finite and Extended Systems[END_REF][START_REF] Bogojeski | Quantum chemical accuracy from density functional approximations via machine learning[END_REF] Although these methods have significant differences, the basic principle is common: using AIMD trajectories, a model is trained and then used to infer ab initio quality PES, and thus perform ab initio quality molecular dynamics simulations for a much lesser computational cost. This approach can be used to access larger system sizes, [START_REF] Lu | PFLOPS Deep Potential Molecular Dynamics simulation of 100 million atoms with ab initio accuracy[END_REF][START_REF] Morawietz | How van der Waals interactions determine the unique properties of water[END_REF] and/or to perform simulations longer than the time reachable using traditional AIMD simulations, in order to observe interesting physical transformations and/or improve statistical sampling. [START_REF] Gerrits | Accurate Probabilities for Highly Activated Reaction of Polyatomic Molecules on Surfaces Using a High-Dimensional Neural Network Potential: CHD3 + Cu(111)[END_REF] dynamics: how to build optimal training sets, and how to critically assess the quality of machine-learning potentials in chemical reactions, a difficult setting where the system is led to explore very energetic configurations far from the geometries of the metastable minima.

An interesting tool in this respect is the "neural-network committee" method. [START_REF] Schran | Committee neural network potentials control generalization errors and enable active learning[END_REF] It consists in training several NNPs on the same training set but with different random seeds. In this way, for the same configuration, the NNPs will give different results. The standard deviation of the predictions of the different members of the committee on some observable is used to assess the reliability of the average prediction. Indeed a good agreement between the different NNP means the configuration is close-enough to the training set for the NNP to be accurate, while a higher value means that the prediction cannot be trusted. This technique can be used to build a training set using an iterative training [START_REF] Behler | Constructing high-dimensional neural network potentials: A tutorial review[END_REF][START_REF] Zhang | A concurrent learning platform for the generation of reliable deep learning based potential energy models[END_REF] procedure. First, a small set of NNPs is trained with existing AIMD data; then, a MD simulation is run using one of the NNP of this set. After this step, configurations that display a standard deviation over the set of NNP on the prediction of some observable above a certain threshold, are recomputed using single point ab initio calculations and added to the training set. This is repeated until no more configurations are evaluated as mis-predicted.

The iterative learning framework has led to thorough studies of systems with ab initio quality at reduced computational cost. [START_REF] Xu | A Deep-Learning Potential for Crystalline and Amorphous Li-Si Alloys[END_REF][START_REF] Zhang | Phase Diagram of a Deep Potential Water Model[END_REF] More recently, the committee method has been also used to quantify the error on an observable computed using NNP-based molecular dynamics simulations, [START_REF] Imbalzano | Uncertainty estimation for molecular dynamics and sampling[END_REF] as well as a way to iteratively select configurations from an AIMD trajectory to build optimal training sets. [START_REF] Schran | Machine learning potentials for complex aqueous systems made simple[END_REF] A recent in-depth study of committee methods [START_REF] Kahle | On the Quality of Uncertainty Estimates from Neural Network Potential Ensembles[END_REF] has shown that in iterative schemes randomly selecting additional configurations to be evaluated at DFT level to improve the training set (random sampling) is equivalent to a selection based on committee disagreement beyond a threshold, but that the latter has to be carefully calibrated.

In this work we adopt NNP-driven enhanced-sampling molecular dynamics simulations to study chemical reactions in solution. This is a challenging goal as it requires the NNP to explore high-energy configurations far from equilibrium, with highly distorted chemical bond geometries. In particular, for a given A → B reaction, we aim at extensively sampling along the RC connecting the two basins, with the aim of reconstructing the accurate free energy landscape through US simulations. To this end, it is crucial to train NNPs capable to yield locally-accurate and well-behaved PES throughout the relevant reaction space.

Although enhanced sampling has been combined with machine-learning potentials in a few recent studies, [START_REF] Yang | Using metadynamics to build neural network potentials for reactive events: the case of urea decomposition in water[END_REF][START_REF] Imbalzano | Uncertainty estimation for molecular dynamics and sampling[END_REF][START_REF] Xu | Accelerating Metadynamics-Based Free-Energy Calculations with Adaptive Machine Learning Potentials[END_REF] including a combination of US with NNP, [START_REF] Pan | Machine-Learning-Assisted Free Energy Simulation of Solution-Phase and Enzyme Reactions[END_REF] the critical assessment and systematic use of NNP for chemical reactions in solution are still lacking. In the same spirit as this work, but using a different angle of approach, the ∆ -M L idea has already been used along with enhanced sampling in a previous study to obtain high accuracy free energy for a more affordable computational cost. [START_REF] Bučko | Ab Initio Calculations of Free Energy of Activation at Multiple Electronic Structure Levels Made Affordable: An Effective Combination of Perturbation Theory and Machine Learning[END_REF] In the present study we present benchmarks and construction principles for training sets. We carefully assess the total computational cost of the training and data production trajectories, with the goal of limiting the total amount of ab initio calculations without losing accuracy. We also introduce a simple approach to ensure long stable trajectories with high NNP-committee agreement: at variance with a previous method where the error is evaluated from deviations between DFT and NNP predictions, [START_REF] Young | A transferable activelearning strategy for reactive molecular force fields[END_REF] our simple scheme avoids the burden of additional ab initio calculations.

We apply the new scheme to two reaction steps of the Strecker-cyanohydrin synthesis of glycine in water, previously studied at the DFT level. [START_REF] Magrino | Step by Step Strecker Amino Acid Synthesis from Ab Initio Prebiotic Chemistry[END_REF] Our results indicate that a surprisingly reduced amount of suitably-chosen ab initio samples is sufficient to train a NNP potential able to sample the full reaction coordinate space between minima, leading to accurate free-energy profiles and barriers at a significantly reduced cost compared to purely DFT simulations.

Methods

Neural network potentials

In this paper, we use the Behler-Parrinello neural networks approach, [START_REF] Behler | Generalized Neural-Network Representation of High-Dimensional Potential-Energy Surfaces[END_REF] in which the potential energy of the system is decomposed as the sum of individual atomic contributions as in equation 1, where N atoms is the number of atoms and E i is the atomic contribution:

E = Natoms i=1 E i (1)
Each atomic contribution is computed using a neural network, and to ensure permutation invariance, the contributions corresponding to atoms of the same species are computed using the same NN.

Here, we use the deepmd/deepPot-Se code [START_REF] Wang | DeePMD-kit: A deep learning package for manybody potential energy representation and molecular dynamics[END_REF][START_REF] Zhang | End-to-end Symmetry Preserving Inter-atomic Potential Energy Model for Finite and Extended Systems[END_REF] based on this architecture. In this framework, a first NN is used to compute symmetry-invariant descriptors for each atom, based on the coordinates of its neighboring atoms within a cutoff sphere. Since this is not the scope of this study, the hyperparameters of the models were not optimized and were taken from other studies. [START_REF] Yang | Using metadynamics to build neural network potentials for reactive events: the case of urea decomposition in water[END_REF] Here, to compute the descriptors we chose a cutoff radius of 6.0 Å and a smoothing cutoff of 1.5 Å. Then, these symmetry-preserving descriptors are given as input to the energy neural networks. The optimization process of the neural networks weights is performed using the Adam optimizer 44 with a learning rate going from 10 -3 to 5.10 -8 during 10 6 iterations on a defined loss decreasing by 5% every 5000 iterations.

Our target system is quite heterogeneous, including 251 atoms in a 13. training. To overcome this potential difficulty, we defined a loss function (equation 2) capable to take into account this heterogeneity:

L(w) = 1 |B| l∈B p E |E l -E w l | 2 + p f 1 N elem N elem i=1 N atoms n i |F l -F w l | 2 (2)
, where n i is the number of atoms of type i in the system, N elem is the number of different elements in the system, N atoms is the total number of atoms in the system, E i and F i denote the DFT energies and forces of the training set, while E w l and F w l are the forces computed by the NNP, and B is the batch size (i.e., the number of trajectory frames).

This choice of the loss function allows weighting equally each element type via the error on forces. The weight p f progressively increase from 1 to 10 while the opposite happens for p e during the training, according to the protocol implemented in the deepmd-kit package. [START_REF] Wang | DeePMD-kit: A deep learning package for manybody potential energy representation and molecular dynamics[END_REF] The neural network used to compute the descriptors was made with 3 layers of 25, 50 and 100 nodes while the energies networks were made with three layers of 240 nodes each. All the AIMD configurations come from a previous study of the group. [START_REF] Magrino | Step by Step Strecker Amino Acid Synthesis from Ab Initio Prebiotic Chemistry[END_REF]

Free energy calculations

Enhanced sampling approaches are necessary to overcome the timescale limitations of molecular dynamics in the simulation of high-barrier rare events like chemical reactions. [START_REF] Pietrucci | Strategies for the exploration of free energy landscapes: Unity in diversity and challenges ahead[END_REF] The technique employed in this study, umbrella sampling, focuses the sampling in the desired regions of configuration space exploiting an external bias potential as a function of a collective variable [START_REF] Torrie | Nonphysical sampling distributions in Monte Carlo freeenergy estimation: Umbrella sampling[END_REF][START_REF] Roux | The calculation of the potential of mean force using computer simulations[END_REF] For the latter, we employ the path collective variables 50 (path CVs) defined as

         s(t) = 1 N -1 N α=1 α exp(-λD[x(t), X α ]) N α=1 exp(-λD[x(t), X α ]) -1 z(t) = -1 λ log N α=1 exp(-λD[x(t), X α ]) . (3)
We introduced a scaling factor 17 so that s ∈ [0, 1]. The first variable s measures the progress of a reaction on a path given by N configurations X α . The second variable z measures the deviation of the configuration at time t from the given path. The metric D used in this study has been extensively used in previous works [START_REF] Pietrucci | Hydrothermal Decomposition of Amino Acids and Origins of Prebiotic Meteoritic Organic Compounds[END_REF][START_REF] Pietrucci | Formamide reaction network in gas phase and solution via a unified theoretical approach: Toward a reconciliation of different prebiotic scenarios[END_REF][START_REF] Pérez-Villa | Synthesis of RNA Nucleotides in Plausible Prebiotic Conditions from ab Initio Computer Simulations[END_REF] and is based on the comparison of coordination patterns:

D[X 1 , X 2 ] = i σ (C iσ (X 1 ) -C iσ (X 2 )) 2 (4)
C iσ is defined in equation 5 for atom i of atom type α with respect to atom type σ as

C iσ = j∈σ 1 - r ij r ασ 0 m 1 - r ij r ασ 0 n (5)
where the term within the sum is a smooth function going from 1 to 0 when the interatomic distance becomes larger than a typical distance r ασ 0 . This distance was set to 1.8 Å for heavy atom-heavy atom pairs and to 1.4 Å for heavy atom-hydrogen pairs. We chose m = 8

and n = 14. The exact definition for the reference structures, explained in our previous work, 17 exploits committor-analysis trajectories performed after a preliminary metadynamics simulation 10 (using a pathCVs defined with only the reactants and the products). 10 regularly-spaced reference structures are selected along the latter trajectories using an inhouse algorithm [START_REF] Magrino | Step by Step Strecker Amino Acid Synthesis from Ab Initio Prebiotic Chemistry[END_REF] based on the nudge elastic band approach. [START_REF] Jónsson | Nudged elastic band method for finding minimum energy paths of transitions[END_REF] Finally, two more references are added at the start and at the end of the path by linearly extrapolating the coordination numbers of the last and first segments to avoid metastable states to appear as spikes in the free-energy landscape: the result are pathCVs s, z based on 12 reference structures that we use throughout this work.

We adopt a set of one-dimensional, quadratic umbrella sampling potentials applied on the s variable and centred at positions s j expressed as:

V bias,j (s) = k 2 (s -s j ) 2 (6)
Windows are equally spaced by s j+1 -s j = ∆s based on k = k B T /(∆s/2.5) 2 in order to have sufficient overlap between two windows (see Ref. [START_REF] Magrino | Step by Step Strecker Amino Acid Synthesis from Ab Initio Prebiotic Chemistry[END_REF] ). The s path-CV measures the progress along a given reaction pathway, the sampling is therefore performed on that coordinate. The z path-CV measures the deviation from the pathway and helps to detect possible anomalies in the sampling. Once a specific reaction pathway is determined and targeted, US is performed on it, in order to determine the corresponding free-energy profile. Due to the intrinsically high-energy, unstable character of configurations explored close to the barrier top, in US simulations of chemical reactions it is occasionally observed that the system can deviate from the targeted pathway to explore a different one. To focus the sampling on the reaction mechanism under study, a restraining potential is sometimes applied along z, as in the case of the step (1) → (2 ′ ) of reference. [START_REF] Magrino | Step by Step Strecker Amino Acid Synthesis from Ab Initio Prebiotic Chemistry[END_REF] We compute the free energy surface, according to:

F (s) = -k B T ln ⟨ρ(s)⟩ (7)
where ⟨ρ(s)⟩ is the unbiased distribution function of s, globally estimated over the whole set of windows using the weighted histogram analysis method 53 implemented in Grossfield's code. [START_REF] Grossfield | WHAM: the weighted histogram analysis method version 2[END_REF] We use a convergence criterion of 10 -7 kcal/mol and 150 bins in s space. We estimate the statistical uncertainty as the deviation between the free energy profiles computed using the third and fourth quarters of each trajectory.

The obtained free energy profile is then used to compute the activation barrier by taking the free energy difference between the highest point in the profile and the reactants free energy; the free energy difference between reactants and products can also be computed.

The value of the activation barriers can then be compared to experimental results as it was done in reference 17. A correction term can also be added to this free energy difference as it was done in reference. [START_REF] Bučko | Dissociative iodomethane adsorption on Ag-MOR and the formation of AgI clusters: an ab initio molecular dynamics study[END_REF] However, since the scope of the present work is to compare NNP results to the AIMD ones of the original study, 17 this additional term was not in this work.

From a ML point of view, the splitting of the RC space into independent US windows simplifies the selection of different data sets to form the training set, and the evaluation of the performance of the NNP.

Results and discussion

We first consider the initial step of the Strecker-cyanohydrin synthesis of glycine, particularly emblematic in prebiotic chemistry studies: [START_REF] Miller | A Production of Amino Acids Under Possible Primitive Earth Conditions[END_REF][START_REF] Saitta | Miller experiments in atomistic computer simulations[END_REF] HCN + COH 2 + NH 3 → CH 2 OH(NH 3 ) + + CN - as illustrated in Fig. 1. The CVs regarding this reaction will be noted with a subscript

(1) → (2 ′ ). [START_REF] Magrino | Step by Step Strecker Amino Acid Synthesis from Ab Initio Prebiotic Chemistry[END_REF] The aim is to assess the behavior of a NNP-system, trained on a minimal amount of AIMD US trajectories, carefully selected along an optimized reaction pathway.

In the last section we will cross-validate the protocol on another, more complex, chemical reaction step of the Strecker mechanism.

Detecting the frontiers of accurate predictions in a neural network potential

The performance of a NNP is customarily measured via the error with respect to DFT on a training set and a testing set. However, this is not sufficient with respect to our goal, i.e., accurately computing the full free energy landscape. In this respect, it is important Figure 1: Free energy diagram of the first step of the Strecker-cyanohydrin synthesis of glycine obtained in our previous work. [START_REF] Magrino | Step by Step Strecker Amino Acid Synthesis from Ab Initio Prebiotic Chemistry[END_REF] The error is of the order of the kcal/mol and was estimated using block averages. The reactants, products, and transition state configurations are also reported on the graph.

to assess the capacity to generate long, stable and accurate MD trajectories. We base our assessment on the committee approach, 32 including a time-dependent metric similarly to a recent study. [START_REF] Young | A transferable activelearning strategy for reactive molecular force fields[END_REF] As shown in figure 2, when, during a simulation, the NNP-generated trajectory exits the "safe" region of configuration space, where forces are accurately predicted, the value of the z pathCV significantly increases, indicating a large deviation from the reference transition pathway (see Methods section) and suggesting the likelihood of unphysical configurations. This is confirmed by the analysis of the corresponding structural properties of the system.

An inspection of the C-O pair correlation functions g(r) (figure 2 c)), before and after this jump, reveals that an unphysical short-distance peak has appeared, not present in the ab initio data. This confirms that the NNP sampled an unphysical region and is trapped in it, as the predicted energy decreases (figure 2 b), a spurious stable configuration. On the other hand, this transition corresponds to a sudden jump in the standard deviation in energy predictions among committee members as well as in the maximum standard deviation on the predicted atomic forces, defined in the following way: where F j i is the force on atom j computed by NNP i, and F j avg is the average over all the neural networks. [START_REF] Yang | Using metadynamics to build neural network potentials for reactive events: the case of urea decomposition in water[END_REF][START_REF] Xu | A Deep-Learning Potential for Crystalline and Amorphous Li-Si Alloys[END_REF][START_REF] Zhang | Phase Diagram of a Deep Potential Water Model[END_REF] We remark that σ max is a more general and accurate indicator of the loss of predictive power than the g(r) or the pathCV z.

σ max = max j∈[1,Natoms] 4 i=1 ∥F j i -F j avg ∥ 2 (8)
This approach allows detecting the frontiers of the configuration space region where the NNP gives accurate prediction. We can thus define a simulation lifetime τ as the time at which σ max surpasses a reliability threshold, simply corresponding to the full simulation time if no pathological behaviors occur. All the trajectory before τ is physically sound, and can be employed to collect statistics. The next logical step is to devise a procedure to maximise τ all along the RC space as a function of the training set composition and extension, in order to carry out reliable US simulations.

Generating training sets suited for free energy calculations

In this section, we critically assess the effect of the composition and size of the training set on the error of the free energy profile recontruced with the NNP potential. The aim is to retain ab initio accuracy while minimizing the amount of DFT calculations necessary to train the potential. Following the approach employed in our previous full-ab initio work [START_REF] Magrino | Step by Step Strecker Amino Acid Synthesis from Ab Initio Prebiotic Chemistry[END_REF] on the reaction considered here, we identify the following algorithm: built upon the reactants and products as the only references: 18 this allows a prejudicefree exploration of a reactive pathway (quick, without the need to converge a free-energy estimate), to be refined through committor analysis and leading to the definition of improved pathCVs based on multiple reference structures (see Methods and Ref. • Plot the maximum lifetime τ max over each US window (see figure 3): if in some RC region τ max is smaller than the autocorrelation time of the pathCV, i.e., if it is impossible to generate uncorrelated samples, it is necessary to generate an additional AIMD US trajectory in that region (since the NNP is locally unreliable) to be added to the training set.

• Repeat the two last steps until a satisfactory NNP is obtained, i.e. with an acceptable lifetime τ across the full RC space, hence capable of producing a dense sampling and a converged US free energy landscape.

The procedure is illustrated in figure 3, showing the training configurations for each iteration (about 600 structures saved every 20 fs) and the resulting NNP samples in the (s, z) plane and lifetime for every US window. Clearly, adding training points according to our scheme allows to progressively increase the simulation lifetime and the overlap between NNP US simulations. In the next two sections, we discuss the error of the NNP on a testing set, and the accuracy of the NNP free energy landscape with respect to the ab initio one.

Error estimate on the benchmark reaction Although this error-evaluation step gives us an idea on the quality of the NNP produced using our new procedure, it cannot be applied on a system that has not been extensively studied with AIMD.

Calculation of the NNP free energy surface for the benchmark reaction

We now assess the accuracy of the free energy landscape reconstructed from the NNP. Figure 5 shows that it is possible to obtain a first-principle quality free energy surface using, in the training set, only one US AIMD simulations and the two end-point AIMD equilibration simulations of the reactants and the products, and to recover the full FES mostly via NNP-US windows. This is a major gain of computational time, which could allow to study, at the same level of ab initio accuracy, larger and more complex systems.

The full original ab initio FES in Ref. [START_REF] Magrino | Step by Step Strecker Amino Acid Synthesis from Ab Initio Prebiotic Chemistry[END_REF] was obtained using 55 US AIMD windows; our procedure significantly reduces, by more than an order of magnitude, the number of ab initio MD simulations to be carried out. From the computational point of view, the training of 1 neural network takes approximately 24 GPUh; while one 5-ps simulation takes about 0.05 GPUh. Hence, for the full study of this test reaction, we needed a total of 161 GPUh, to which one needs to add the 40k CPUh used to build the ab initio training set. The fact that all the calculations were run on GPU makes non-trivial the comparison with respect to the CPU simulations. However, the whole ab initio study of this reaction in Ref. [START_REF] Magrino | Step by Step Strecker Amino Acid Synthesis from Ab Initio Prebiotic Chemistry[END_REF] took around 700k CPUh.

In figure 5, the free energy was computed using short simulations of 5 ps per window. In the next section, we propose a method to generate longer, but still stable, NNP trajectories.

Figure 5: Panel a): US-obtained free energy profiles using full AIMD (black line), and using a NNP trained on AIMD equilibration of reactants and products (green line). The shaded zones correspond to the estimated statistical errors. Panel b): US-obtained free energy profiles using full AIMD (black line), and using a NNP trained on AIMD equilibration of reactants, products, and the transition state US windows (green line). The red line represents the free energy which would be obtained using the three AIMD US windows (reactants, products, transition state). The reactants and products endpoint-configurations are also reported on the graphs.

Generating stable NNP trajectories

After properly training our NNP, our aim is to be able to generate long and stable MD simulations. However, it is known that, at some point, the NNP-system will eventually explore untrained regions of the configuration space, where energies and forces will necessarily be extrapolated, likely failing to describe the correct system anymore. To overcome this common drawback, one possibility to restrain the sampling to low σ max regions is to add a quadratic potential on σ max , as it was proposed in reference. [START_REF] Schran | Committee neural network potentials control generalization errors and enable active learning[END_REF] In this work, instead, we find it more effective to completely avoid high-σ max values, which would correspond to an infinite spring constant with respect to the above-mentioned reference. Therefore, we can define such a "mirror reflection operation" as follows:

We identify two regions in the configuration space sampled by each US window: the region where he NNP behaves properly, i.e. where the maximum standard deviation on the prediction of forces σ max is low, and the region where significant deviations, hence unphysical behavior, are observed.

Whenever σ max exceeds a given threshold, we define for each atom α the following vector G α , the value of this vector in the k direction is:

G αk = 1 n comm ncomm i=1 (F i αk -F avg αk ) 2 (9)
With n com the number of neural networks in the committee (n comm = 4 in this study). This quantity provides the normal vector of a frontier between the reliable, "interpolation region" of the NNP, and a NNP-unknown "extrapolation region". We can therefore mirror-reflect the trajectory of the atoms on this hypersurface when σ max displays a pathological jump to high values.

In order to do so, and force the NNP-system to remain in the interpolation region, the simulation is stopped a t step number of timesteps before the instability, and is restarted after imposing the following transformation on the velocities:

v new α = v old α - 2v old α .G α ∥G α ∥ 2 G α (10)
where v old α is the velocity of atom α at the moment in which we mirror-reflect the trajectory, This operation preserves the kinetic energy of the system, without pathological effects on the US procedure. We consider that it is a useful and legitimate complement to our NNPbased free-energy calculation protocol whenever the typical lifetime τ is significantly longer than the auto-correlation time of the CV, so that uncorrelated samples can be collected between successive reflections. In the opposite case, it is advisable to improve the NNP via a larger training set before employing it to perform statistical sampling.

To illustrate the choice of the parameters of the protocol, the time evolution of σ max during a US simulation is reported in figure 6 along with the NNP-autocorrelation function of the CV. The threshold of σ max is chosen as 0.6eV/A, while choosing the number of timesteps at which the reflection is performed before the instability between 200-500 gives similar results. Clearly, the typical time lapse between two reflections (ranging here from 0.3 to 5 ps) is much larger than the autocorrelation time (∼ 0.1 ps). Reflections therefore appear not harmful from a statistical viewpoint, and allow to carry out long, stable US simulations localized in the region of configuration space where the NNP is reliable, leading to free-energy estimations retaining ab initio accuracy (Figure 7 and table 3). The profile obtained in figure 7 lies within the typical 2 kcal/mol uncertainty of the PBE functional for this kind of studies, as well as the predicted free energies, shown in table 3.

Table 1: Activation barrier (∆F † ) and free energy difference between reactants and products(∆F Application of the protocol to a more complex reaction [START_REF] Magrino | Step by Step Strecker Amino Acid Synthesis from Ab Initio Prebiotic Chemistry[END_REF] The error is of the order of the kcal/mol and was estimated using block average. [START_REF] Magrino | Step by Step Strecker Amino Acid Synthesis from Ab Initio Prebiotic Chemistry[END_REF] The relevant configurations are also reported on the graph.

The aim of this part of our work is to challenge the method previously devised to reproduce a complex FES. To this end, we chose another intermediate step of the Streckercyanohydrin synthesis, presenting a more complex free energy landscape, consisting of a two-step process, as shown in figure 8: a hydrogen bond breaking between the imine nitrogen, a hydrogen atom and the cyanide, followed by the addition of the cyanide to the imine.

This process displays a small barrier of 5 kcal/mol, a small drop of 3 kcal/mol towards a metastable state, another small barrier of 3 kcal/mol with respect to the latter step, followed by a large 23 kcal/mol drop to the products, as obtained from our original AIMD FES in Ref. [START_REF] Magrino | Step by Step Strecker Amino Acid Synthesis from Ab Initio Prebiotic Chemistry[END_REF] by using 44 US windows. Such small free-energy barriers are likely more difficult to be quantitatively and qualitatively described by a NNP, than in the previous reaction.

We underline that, although the stoichiometry of the chemical species is the same as in the previous case, the reactants and the products are different, and thus the NNPs are generated from scratch, independently from the previous case.

Results of US simulations employing a series of three NNPs with different training sets are presented in figure 9. We progressively enlarged the training set as previously described, by using at first 1, then 3 and finally 6 AIMD US windows, besides the reactants and products ones (respectively 3, 5 and 8 training points) .

The NNP FES obtained starting from the different training sets are shown in figure 10.

Judging from the accuracy of the FES, NNP training is optimal when 6 AIMD US windows are used, which is reasonable considering the fact that this reaction consists of two successive steps. Instead, unsurprisingly, the FES obtained with only 3 training points compares very poorly with the ab initio one.

Adding AIMD US windows to the training set leads to a progressive increase of the accuracy of the NNP FES, until a satisfactory agreement with the benchmark AIMD one is achieved. As in the case of the previous reaction, the cumulative duration of the AIMD trajectories necessary for training the NNP is one order of magnitude shorter than the total duration required for computing an accurate AIMD FES, despite the fact that the present reaction is clearly more challenging to be reproduced in its fine details. The free energy profile obtained using the velocity transformation technique is shown in figure 11. We report in Table 2 the CPU time comparison between a full ab initio FES calculation and a NNPbased one. As we have no guarantee that our NNPs are transferable to another reaction, we chose to add the CPU time needed to perform the training simulations in the comparison with the ab initio study 4) ) and free energy difference between reactants and products(∆F We also provide an approach to exploit the disagreement between predictions of equallytrained NNPs not only as diagnostics but also to render robust and stable the long-time dynamics, hence to achieve satisfactory statistical sampling, a crucial advantage for freeenergy reconstruction.

Our investigation suggests a new protocol for the accurate, ab initio quality, NNP calculation of FES of chemical reactions in solution: the computational load of pure and costly AIMD simulations would be limited to preliminary metadynamics exploration with mechanism-agnostic general-purpose CVs, [START_REF] Pietrucci | Formamide reaction network in gas phase and solution via a unified theoretical approach: Toward a reconciliation of different prebiotic scenarios[END_REF] followed by committor analysis to define mechanism- We expect our approach to be easily generalizable to a range chemical reactions in solution, allowing a limited and incremental use of costly AIMD calculations only if and when needed. In perspective, this controlled and efficient scheme will help to exploit NNPs to overcome a significant computational bottleneck in the accurate calculation of free-energy profiles in solution chemistry.

  4 Å-side cubic box under periodic boundary conditions, with the following composition: 2 C, 2 N, 6 H and 1 O atom in the reactive molecules, solvated with 81 water explicit water molecules, whose H and O atoms are indistinct from the ones in the reactive part of the system, and treated exactly on the same footing by the NN. As a consequence, the water molecules would comparatively acquire an overwhelming weight, with respect to the important reactive subsystem, during

  We used LAMMPS 45 for all the NNP-based molecular dynamics simulations, employing hydrogen atoms instead of deuterium in the original study to have a quicker, more realistic dynamics. All the enhanced sampling simulations performed with neural networks were carried out using the open-source, community-developed PLUMED library, 46 version 2.5.0 47 .

Figure 2 :

 2 Figure 2: Panel a): left (red line): time evolution of z path CV; right: time evolution of the maximum standard deviation on the prediction of forces(σ max , blue line) and standard deviation on the prediction of energies(σ energy , green line). Panel b): time evolution of the potential energy predicted by the four neural networks. Panel c): Carbon-Oxygen radial distribution function for the AIMD simulation (black line) and for the NNP simulation, before (green line) and after (red line) the unphysical jump in the predictions of energy and forces.

Figure 3 :

 3 Figure 3: Panel a): maximum lifetime τ max (red dots) in each US window throughout the path CV for a NNP trained only on the AIMD equilibrations of the reactants and the products; the AIMD training points are represented in blue in the (s,z) plane. Panel b): corresponding instantaneous location of the NNP configurations in the (s,z) path CV coordinate plane; panel c) maximum lifetime τ max (red dots) in each US window throughout the path CV for a NNP trained only on the AIMD equilibrations of the reactants, the products, and a transition state US window; the AIMD training points are represented in blue in the (s,z) plane. Panel d): corresponding instantaneous location of the NNP configurations in the (s,z) path CV coordinate plane

17

 17

••

  Generate one AIMD US trajectory (of about 15 ps) in the transition state window along the optimized pathCV (defined at the previous step), and include it in the training set along with AIMD equilibration trajectories (of about 15 ps) of the reactants and the products. These three trajectories represent ab initio "milestones" of the RC to train and test the NNPs. Train four models on the same training set, with different random seeds. Define a range of US windows densely spanning the full RC range, and for each window perform 50 short (about 5 ps) NNP-based US simulations, with the same starting point but different random initial velocities taken from the Boltzmann distribution.

Figure 4 :

 4 Figure 4: Panel a): Root mean squared error (RMSE) of the NNP-estimated energy along the chemical path for a NNP trained only on the AIMD equilibrations of the reactants and the products (blue line) and for a NNP trained on those configurations and on the US transition state window (orange line). The test set is built only with configurations that are not present in any training set. Panel b): same as panel a) but the RMSE was computed on forces instead of energies

Figure 6 :

 6 Figure 6: Panel a): Time evolution of the maximum standard deviation on the prediction of forces (σ max ) during a NNP umbrella sampling simulation. We report with black crosses its instantaneous values before the the unphysical behavior sets in, and with red crosses its behavior after that. Panel b): NNP-autocorrelation function of the s reaction coordinate computed in the same NNP-umbrella sampling window.

Figure 7 :

 7 Figure7: Panel a): US-obtained free energy profiles using full AIMD (black line), and using NNP trained on AIMD equilibration of reactants, products and transition state US windows when the velocity transformation is applied (green line). The red line represents the free energy which would be obtained using only the three AIMD US windows (reactants, products, transition state). The reactants and products endpoint-configurations are also reported on the graph. Panel b): corresponding instantaneous location of the NNP configurations in the (s,z) path CV coordinate plane. As in the original study 17 z path CV was confined to within 0.12 via a semiparabolic wall.

  (1)→(2 ′ ) ) obtained in the ab initio study along with the ones obtained in this work Ab initio study machine learning study ∆F † (1)→(2 ′ ) (kcal/mol) 10 8 ∆F (1)→(2 ′ ) (kcal/mol) -14 -14.5 ± 1

Figure 8 :

 8 Figure 8: Free energy diagram of the reaction (3) → (4) of the Strecker-cyanohydrin mechanism.[START_REF] Magrino | Step by Step Strecker Amino Acid Synthesis from Ab Initio Prebiotic Chemistry[END_REF] The error is of the order of the kcal/mol and was estimated using block average.[START_REF] Magrino | Step by Step Strecker Amino Acid Synthesis from Ab Initio Prebiotic Chemistry[END_REF] The relevant configurations are also reported on the graph.

Figure 9 :

 9 Figure 9: Illustration of the iterative procedure followed to build a converged NNP for reaction (3) → (4). Panel a): maximum lifetime τ max (red dots) in each US window throughout the path CV for a NNP trained only on the AIMD equilibrations of the reactants, the products, and the transition state; the AIMD training points are represented with blue crosses in the (s,z) plane. Panel b): same as a) but adding two additional AIMD US windows in the NNP training set. Panel c): same as b) but adding three additional AIMD US windows in the NNP training set.

20 Conclusions

 20 (3)→(4) obtained in the ab initio study along with the ones obtained in this work Ab initio study machine learning study ∆F † (3)→(4) (kcal/mol) 3 3.5 ± 1 ∆F (3)→(4) (kcal/mol) -20 -In this work, we perform a critical assessment of different procedures to build a NNP training set starting from AIMD US trajectories of chemical reactions in solution, exploiting two reactions along the Strecker pathway 17 as test cases. The systematic comparison of AIMD training sets increasingly more dense along the reaction coordinate clearly indicates a threshold for achieving NNP free energy profiles of ab initio accuracy.

Figure 10 :

 10 Figure 10: Panel a): US-obtained free energy profiles using full AIMD (black line), and NNP trained on AIMD equilibration of reactants, products and transition state windows (green line). The red line represents the free energy which would be obtained using only the three AIMD US windows (reactants, products, transition state). The shaded zones correspond to the estimated statistical errors. Panel b): same as a) but using the five training points in panel b) of figure 9. Panel c): same as b) but using the eight training points in panel c) of figure 9. The reactants and products endpoint-configurations are also reported on the graphs.

Figure 11 :

 11 Figure 11: Panel a) free energy profile obtained as in panel c) of figure 10 but applying our velocity transformation (equation (10)). Panel b) corresponding instantaneous location of the NNP configurations in the (s,z) path CV coordinate plane.

Table 2 :

 2 Summary and comparison of the computational times for the (3) → (4) reaction between the pure AIMD protocol and the combined AIMD-ML one.

		Ab initio study	machine learning study
	training time	0	24h(GP U ) × 4 × 3 iterations
	US ab initio time	15 × 44 = 660ps	15 × 6 = 90ps
	CPU/GPU simulation time	540k CPU.h	100k CPU.h + 200 GPU.h
	Total CPU/GPU time	540kCPU.h	488 GPU.h + 100k CPU.h

Table 3 :

 3 Activation barrier (∆F † (3)→(

Acknowledgement

This work was granted access to the HPC ressources of IDRIS (Institut du Développement et des Ressources en Informatique Scientifique) under the allocations 2020-A0080901387 and 2021-A0100910143 attributed by GENCI (Grand Equipement National de Calcul Intensif).


domain: spi phys

text:
Luc Frison
email: luc.frison@onera.fr

Julien Queyrel
email: julien.queyrel@onera.fr

Laurent Feral
email: laurent.feral@aero.obs-mip.fr

Characterization and modeling of tropospheric propagation at high latitudes for satcom systems

teaching and research institutions in France or abroad, or from public or private research centers.

Characterization and modeling of tropospheric propagation at high latitudes for satcom systems

Introduction

Due to global warming, the polar regions are increasingly accessible and offer new opportunities that require the deployment of terrestrial or space telecommunications systems. For example, Earth observation satellites in polar orbit could take advantage of their increased daily visibility to transfer data to high latitudes ground stations more often; ships using the new maritime routes also need access to Satcom services provided by satellites in geostationary orbit. However, at these latitudes, the radio links are established at very low elevation angles (a few degrees at most). Consequently, the associated radioelectric paths through the troposphere are considerably increased, which intensifies the propagation effects (attenuation, refraction, scintillation) classically observed at lower latitudes. On the other hand, in order to reach high data rates, the satcom systems considered must operate beyond 20 GHz. However, at these frequency bands, the total tropospheric attenuation due to oxygen, water vapor, clouds and hydrometeors in liquid phase (rain) or mixed phase (wet snow) can become prohibitive [START_REF] Queyrel | Preliminary Results of the THOR7 Propagation Experiment in the North Pole Region[END_REF]. To these attenuations are added small scale turbulent fluctuations of the tropospheric refractive index. They are responsible for scintillation effects that become very significant for such low elevation links. Finally, at larger scales, specific refraction effects (surface or altitude duct) are likely to appear along these extended slant paths.

In order to quantify these propagation effects for the optimal design of polar Telecom systems, empirical propagation models are usually derived from experimental campaigns. Unfortunately, such propagation experiments are costly and difficult to maintain in Polar Regions. Following [START_REF] Queyrel | Preliminary Results of the THOR7 Propagation Experiment in the North Pole Region[END_REF], [START_REF] Mire | these: Modelling of Earth-space propagation in Ka-Band in tropical and equatorial regions[END_REF] or [START_REF] Fayon | Modelisation statistique de la diversité multi site aux fréquences comprises entre 20 et 50 GHz[END_REF], our approach consists in modeling the various tropospheric propagation effects by coupling a highresolution meteorological model, PWRF-ARW (Polar WRF [START_REF] Polar | Polar WRF, Model Item[END_REF], polar version of the Advanced Research Weather Research and Forecasting WRF-ARW model [START_REF] Skamarock | A Description of the Advanced Research WRF Version 3[END_REF]) to an ElectroMagnetic Module (EMM). The set PWRF-EMM defines the Atmospheric Numerical Simulator ANS which is presented in Section 2. Its performances will be evaluated on a single day by comparing with meteorological and propagation data, presented in Section 3. Finally, the results are given in Section 4.

The Atmospheric Numerical Simulator (ANS)

The PWRF-ARW model is used to produce a 3D description of the atmosphere every 5 min over an area of about 200×200 km² (horizontal resolution 2×2km²) centered halfway through the ground projection of the SvalSat (78.25°N, 15.50°E) -THOR7 tropospheric slant path. The outputs are sampled over 38 vertical model levels (ranging from ground level to ~20km). The ERA5 re-analysis data ( §3) are used for the initial and boundary conditions of the simulation. In order to mimic realistic past states of atmosphere, the simulations length is limited to one day with an initialization period of 12 hours (the results of which are later discarded) allowing the model to provide stable outputs. In compliance with Figure 1, to reach the final resolution of 2x2km², two nested domains are used: a coarser domain d01 (horizontal resolution of 10x10km² constructed with a polar stereographic projection and static fields initialized directly by interpolation from the ERA5 database), finer domain d02 nested at the center of the former, at the desired resolution (2x2km²). Each output 3D voxel of a datacube provides information on several meteorological variables (temperature T, pressure P, humidity H, concentration of precipitating hydrometeors, etc.) which depend on an exhaustive parametrisation, part of which is the microphysical scheme (interaction models between the various water phases). In our case, a cumulus scheme is activated for the large domain d01 only (Tiedtke diagram [START_REF] Zhang | Projected Future Changes of Tropical Cyclone Activity over the Western North and South Pacific in a 20-km-Mesh Regional Climate Model[END_REF]). To reduce the computational time, d01 is hydrostatic and d02 non-hydrostatic (defining d01 as a non-hydrostatic domain is timeconsuming but does not change the simulation results significantly). The spectral nudging model provided within PWRF is also used in order to generate realistic simulations when compared to the ERA5 reanalysis data.

An ElectroMagnetic Module (EMM) written in Python then converts the atmospheric parameters of interest (pressure P, temperature T, water vapor content 𝜌 𝑣𝑎𝑝 , cloud water content 𝜌 𝑐𝑙𝑜𝑢𝑑 , rain water content 𝜌 𝑟𝑎𝑖𝑛 ) within the 3D meteorological voxels simulated every 5 min by PWRF into specific attenuations (in dB/km) for each tropospheric contributor to the total attenuation: oxygen, water vapor, cloud, rain and snow. The integration of the specific attenuation along the line of sight towards THOR7 allows to recover the attenuation suffered by the signal through the troposphere. The physical models used by the EMM are: the line-by-line spectral model of recommendation ITU-R P676-11 for oxygen and water vapor and the model based on Rayleigh scattering of ITU-R P840 -7 for cloud droplets.

The granulometry of the precipitating hydrometeors (rain, dry and wet snow) is estimated from the content (g/m3) and the number of particules given by PWRF in each voxel. Here, the Thompson microphysical scheme [START_REF] Morrison | Impact of Cloud Microphysics on the Development of Trailing Stratiform Precipitation in a Simulated Squall Line: Comparison of One-and Two-Moment Schemes[END_REF] is used because it is well suited for high latitudes simulations, providing information on several types of polar hydrometeors such as snow, ice and graupel in addition to rain. The granulometry of each water phase component 𝑥 is then given by:

𝑵 𝒙 (𝑫 𝒙 ) = 𝑵 𝑻𝒙 𝟏 𝚪(𝟏 + 𝜶 𝒙 ) [ 𝑵 𝑻𝒙 𝑴 𝒔𝒙 𝝅 𝟔𝝆 𝒙 𝚪(𝟒 + 𝜶 𝒙 ) 𝚪(𝟏 + 𝜶 𝒙 ) ] 𝟏+𝜶 𝒙 𝟑 𝒆 -[ 𝑵 𝑻𝒙 𝑴 𝒔𝒙 𝝅 𝟔𝝆 𝒙 𝚪(𝟒+𝜶 𝒙 ) 𝚪(𝟏+𝜶 𝒙 ) ] 𝟏 𝟑 𝑫 𝒙 (𝟏)
With x ∈ [snow, rain, ice, graupel], N x the number of hydrometeors with diameter D x per unit volume, α x = 0 [START_REF] Polar | Polar WRF, Model Item[END_REF], M sx the density of the hydrometeor x, Γ the gamma function, ρ x the mixing ratio and N Tx the concentration of component x. Importantly, note that ρ x and N Tx are PWRF outputs.

The phase of the hydrometeors (rain, ice, dry snow, wet snow) is also taken into account by using the Sihvola equivalent permittivity model [START_REF] Sihvola | Mixing Rules with Complex Dielectric Coefficients[END_REF]. It makes it possible to determine the equivalent permittivity of a mixed hydrometeor as a function of the mass fractions composing it. The single-phase permittivity is computed from the double Debye model [START_REF] Ray | Broadband complex refractive indices of ice and water[END_REF]. Then, the associated specific attenuation 𝛾 𝑥 is calculated using the theory of scattering by spherical (Mie formulation for rain) or spheroidal (T-matrix formulation for dry and wet snow) scatterers for which the axis ratio is a function of the size and type of hydrometeor.

𝜸 𝒙 = 𝟒. 𝟑𝟒𝟑 × 𝟏𝟎 𝟑 ∫ 𝝈 𝒆𝒙𝒕 (𝑫)𝑵(𝑫)𝒅𝑫 (𝟐) 𝑫 𝒎𝒂𝒙

𝟎

With 𝛔 𝐞𝐱𝐭 (𝐃) the extinction cross section derived from Mie theory or from the Tmatrix for respectively x=rain or x=snow, 𝐃 𝐦𝐚𝐱 is the maximum diameter of a hydrometeor of category x. 𝐃 𝐦𝐚𝐱 = 6mm for rain and 𝐃 𝐦𝐚𝐱 = 8mm for snow [START_REF] Oguchi | Electromagnetic wave propagation and scattering in rain and other hydrometeors[END_REF].

Following the approach proposed for rain in [START_REF] Pruppacher | A wind tunnel investigation of the internal circulation and shape of water drops falling at terminal velocity in air Quarterly[END_REF] and considering the limit values of axis ratio found by Garret [START_REF] Garrett | Orientations and aspect ratios of falling snow[END_REF] in measurement -𝒂 = 𝟎. 𝟒𝟓 for 𝑫 𝒆,𝒎𝒂𝒙 = 8mm [START_REF] Oguchi | Electromagnetic wave propagation and scattering in rain and other hydrometeors[END_REF], the following relation for the axis ratio for snow has been adopted: 𝒂(𝑫 𝒆 ) = 𝟎. 𝟗𝟔𝟒𝟑 -𝟎. 𝟎𝟔𝟒𝟑 * 𝑫 𝒆 𝟎 < 𝑫 𝒆 < 𝟖𝒎𝒎 With 𝐃 𝐞 , the equivalent diameter of water content In order to calculate the total attenuation suffered by the signal along its path towards the beacon receiver at low elevation angles, large-scale geometric refraction effects are also taken into account. First, for each time step, a vertical slice of atmosphere, directed towards the satellite is extracted from the PWRF data cube. Next, a dichotomy algorithm is used to estimate the true elevation angle of the geometric path between the beacon receiver and the satellite by propagating a geometric ray through the different layers of atmosphere within the slice. The refractive index (computed from P, T and 𝛒 𝐯𝐚𝐩 with Rec. ITU-R P.453-13) is used to apply the geometric laws of Snell/Descartes between two successive layers as suggested in ITU-R P.676-11. In order to propagate the signal to THOR7, starting at SvalSat the process is repeated for different elevation angles until the beacon receiver and satellite are joined within an acceptable margin of error. The specific attenuation is then integrated along the newly found geometric path to recover the total attenuation (dB) at the current time step.

Experimental and re-analysis database for validation.

In order to assess the performance of the ANS, various datasets are used as a basis for comparison.

ERA5 [13]

Although used as initial and boundary conditions in PWRF, the ERA5 reanalysis data produced by the European Center for Weather Prediction (ECMWF) are also considered to evaluate the weather outputs produced by our model. ERA5 data provide an hourly estimate of various meteorological quantities (temperature, pressure, cumulative precipitation, etc.) at any point on the globe with a horizontal resolution of 0.25x0.25deg² on 38 vertical pressure levels as well as single level and surface quantities. [START_REF] Seity | The AROME-France Convective-Scale Operational Model[END_REF] PWRF simulations will be compared with Arome-Arctic reanalysis provided by the Norwegian Meteorological Institute. These reanalysis data come from the Arome meteorological model developed by Meteo-France centered on the Arctic zone (horizontal resolution 2.5×2.5km², 65 vertical pressure levels).

Arome Arcitc

CARRA [15]

The simulations will also be compared with the Copernicus Arctic Regional Reanalysis (CARRA) dataset. It contains 3-hour analyses and short-term hourly forecasts of atmospheric and surface weather variables such as temperature and humidity (but no rainfall data) at 2.5x2.5km² resolution over 23 vertical pressure levels.

3.4 THOR7 experiment (WX) [START_REF] Houts | Design of a Ka-band propagation terminal for atmospheric measurements in polar regions[END_REF] In order to study the specificities of the Ka-band propagation channel at polar latitudes, ESA financed a two-year campaign (April 2016 to March 2018) of propagation at 20 GHz between the geostationary satellite THOR-7 and a NASA beacon receiver installed at the Svalbard Satellite Teleport (SvalSat) operated by KSAT on the island of Spitsbergen in the Svalbard archipelago of Norway. With the help of NASA and additional funding from CNES, ONERA was in charge of analyzing the propagation data [START_REF] Queyrel | Preliminary Results of the THOR7 Propagation Experiment in the North Pole Region[END_REF]. The experimental ground station consists in a beacon receiver designed by the NASA Glenn Research Centre [START_REF] Houts | Design of a Ka-band propagation terminal for atmospheric measurements in polar regions[END_REF] that operates in the Ka-band at the 20.198 GHz beacon frequency of THOR7which is seen with a 2.65° elevation angle from SvalSat. A Ka-band radiometer (Radiometric PR-2230) that can observe multiple channels between 22 and 30 GHz is also available on site. Other meteorological instruments such as tipping bucket rain gauge and a weather station equipped with an IR cloud sensor complement the equipment deployed for this experiment (see Figure 2). During the measurement campaign, it was noted that meteorological conditions (precipitation and low temperature) could interfere with the proper functioning of the tipping bucket (freezing problems) and generate faulty measurements.

ANS simulator validations

Due to the large computation time of the PWRF simulations, we focus in this study on a single daily simulation. The 2017-02-21 was selected by searching for a snowy day with attenuation events within the precipitation and temperature data from Arome-Arctic and in the in-excess attenuation time series obtained by plotting THOR7 data completed with the concurrent WX meteorological measurements collected at the SvalSat station (temperature during the day always below 273K and precipitation measurement).

Meteorological outputs

In order to verify the ability of our model to reproduce the states of the polar atmosphere, we will first focus on meteorological variables from PWRF simulations. In particular, the evolution of temperature and humidity during the 2017-02-21 at the coordinates of the SvalSat station (latitude: 78.22, longitude: 15.42, height: 465m) are shown in 3 and 4 respectively. Note that the reference data and the simulated data do not have the same temporal and spatial sampling: Table 1 summarizes the differences.

Table 1 : Temporal and spatial sampling for ERA, CARRA, Arome Arctic, WX and PWRF

In Figure 3 and4, the results of the PWRF simulations are plotted in black, those from the ERA5 reanalysis in red, in yellow those related to Arome-Arctic, and in blue the CARRA data. The meteorological data (WX temperature in Figure 3, WX humidity in Figure 4) collected at the SvalSat station during the THOR7 propagation campaign are reported in green. All simulations and reanalysis data (PWRF, CARRA, ERA5 and Arome Arctic) are interpolated at the height of the THOR7 experimental station (465m) in order to obtain comparable curves. Figure 3 shows a good agreement between the PWRF simulations and the reanalysis or the WX weather station data: PWRF seems to produce temperature outputs consistently with the four sets of references. However, Figure 4 shows that the simulated humidity time series at the position and height of the SvalSat station are less comparable with the reanalysis or WX time series. Nevertheless, despite a specific moisty event in the first part of the day, the model seems to reproduce the experimental and reanalysis data quite faithfully.

The next step towards the PWRF weather model validation effort focuses on its ability to reproduce accumulated precipitation contents during the whole day. The microphysical scheme selected to run PWRF in this study makes it possible to isolate four distinct types of precipitations: rain, snow, ice and graupel. We consider the amount of total precipitation as the addition of the number of millimeters of equivalent water of each kind. Figure 5 shows the time series of cumulative precipitations within the pixel centered on the SvalSat station coordinates for the 2017-02-21. The data from ERA5 are also plotted for information only (the 0.25°x0.25° horizontal resolution of ERA5 is very different from the 2x2km² horizontal resolution of Arome-Arctic and PWRF. This result highlights the difficulties for PWRF to accurately reproduce precipitating events. It is however observed that the simulations (PWRF, AROME and ERA5) are in agreement with the measured data (THOR7) on the presence of a snow (precipitations et temperature below à 0) event occurring during the first 600 minutes of the day (00h-10h UTC) and followed by a calm afternoon. The occurrence of mainly stratiform precipitation events (uniform over larger scales) in the region could explain the fairly good agreement between the ERA5 data and the other sources. The two simulations (PWRF and Arome-Arctic) seem to overestimate the total daily accumulated amount of precipitation compared to the WX rain gauge measurements, though freezing issues (see §3.4) may have cause some problems.

Attenuations outputs

The ability of the PWRF meteorological model coupled with the EMM module to reproduce the propagation channel behavior is presented in this section. In this study, the in-excess attenuation is considered as it can be easily obtained from the propagation experiment, without having to resort to a radiometer for the clear sky contribution. On the simulation side, the attenuations due to the hydrometeors only (cloud, rain, snow etc…) are considered in order to mimic a synthetic in-excess attenuation.

In the previous study [START_REF] Queyrel | Preliminary Results of the THOR7 Propagation Experiment in the North Pole Region[END_REF], the results showed good correlations between the simulated attenuation time series and the experimental ones. Without a snow model, the ANS seemed capable of reproducing rainfall events, but proved unable to reproduce the CCDF (Complementary Cumulative Distribution Function) of in-excess attenuation. Although snow is not known to trigger strong attenuation events, we intend to see if the addition of a snow attenuation model makes it possible to fix this underestimation.

Figure 6.a shows the experimental (green) and synthetic (black) time series of in-excess attenuation. Because the experimental time series at 10 Hz is very noisy, it is first dowsample to 1Hz with a moving average. Then, a Butterworth low-pass filter (cuttof frequency=0.01 and sampling frequency fs=1) is used to remove the scintillation. To clarify the results a bit, we only display one point out of ten in the figure, but all the points are taken into account in figure 6.b. We clearly observe a concordance between the synthetic time series and the experimental one: the precipitation event is accounted for at the same time and with similar intensities, though the internal dynamics do not exactly match. Some of the differences could be explained by the different temporal sampling: a point every 5 min for the PWRF simulations and a point every second for the experimental data. Figure 6.b represents the CCDF of in-excess attenuation derived from the PWRF synthetic time series and from the experimental. The statistics for the day are comparable even though the PWRF generates more attenuation. The probability of occurrence is higher with the experimental data because of the residual noise and scintillation encountered in the time series: noisy values below a few tenths of dB cannot be associated with actual in-excess attenuation. These statistics will have to be verified over longer study periods in order to really conclude on the capability of the model to produce reliable attenuation statistics.

Conclusions

In order to understand and quantify the effects of the polar meteorological conditions on microwave tropospheric propagation, the approach consisting in coupling a meteorological model with an EM module ([1], [START_REF] Mire | these: Modelling of Earth-space propagation in Ka-Band in tropical and equatorial regions[END_REF], [START_REF] Fayon | Modelisation statistique de la diversité multi site aux fréquences comprises entre 20 et 50 GHz[END_REF]) is picked-up again and developed further with the introduction of a new multiphasic hydrometeor attenuation model. The WRF-ARW model previously used in [START_REF] Queyrel | Preliminary Results of the THOR7 Propagation Experiment in the North Pole Region[END_REF] is traded for PWRF, a model specifically developed for high latitudes. First, the meteorological simulations are compared to reference reanalysis dataset as well as meteorological measurements at the location of the THOR7 beacon receiver. The results for temperature and humidity are quite satisfactory. The simulations also seem to fairly well reproduce precipitation events despite some inconsistencies in precipitation intensity (accumulated precipitations at the end of the day). Second, the simulated in-excess attenuation time series are compared to the experimental ones obtained from the THOR7 beacon. The ANS simulations present a good temporal/intensity agreement with the data as far as in-excess attenuation events are concerned. Although one day is a short period, comparing the attenuation statistics (CCDF) show good promise for the of the study. The next steps consists in selecting several representative days of polar climatology for further simulations, allowing to improve the parametrization of PWRF and of the snow attenuation model . In the end, the objective is to simulate the 2 years of the THOR7 propagation measurements with the best possible PWRF and EMM parametrization in order to obtain a good match between synthetic and measured statistics.

Figure 1 :

 1 Figure 1: The two nested d01 end d02 domains of the PWRF-ARW simulations (a) and position of the high-resolution domains d02 with the horizontal projection of the line of sight towards the THOR7 satellite (b)

Figure 2 :

 2 Figure 2 : SvalSat station (78.25°N, 15.50°E)

Figure 3 :

 3 Figure 3: Temperature temporal evolution, SvalSat, 02-21-2017 from PWRF simulations and from the various reference meteorological databases.

Figure 5 :

 5 Figure 5 : Cumulative precipitation amount (water equivalent) on the SvalSat station on 02/21/2017, from PWRF simulations and from the various reference meteorological databases.

Figure 4 :

 4 Figure 4 : Humidity temporal evolution, SvalSat, 02-21-2017 from PWRF simulations and from the various reference meteorological databases.

Figure 6 :

 6 Figure 6 : Figure 6 : (a): In-excess attenuation time series (synthetic from PWRF simulations and experimental from THOR7 propagation measurements), SvalSat, 02-21-2017; (b) associated in-excess attenuation CCDF


text:
Étienne Suquet
email: etienne.suquet@onera.fr

Laurent Castanet

Julien Queyrel

Laurent Féral

Xavier Boulanger

Carlo Riva

FREQUENCY SCALING OF RAIN ATTENUATION FOR VERY HIGH THROUGHPUT SATELLITE SYSTEMS

In this paper, we present a methodology to compare frequency scaling techniques of rain attenuation. This methodology consists in obtaining rain attenuation from beacon measurements, selecting a subset of events, and applying metrics on the obtained time series. Are also discussed methodologies to separate the contributions of clouds and rain to attenuation, as well as a method to identify events affected by noise. The relevance of the approaches is evaluated on a database of satellite beacon attenuation: Italsat beacon measurements from Spino d'Adda (7 years).

Introduction

Some Very High Throughput Satellite (VHTS) systems in preparation plan to use Q/V band for their feeder links and Ka band for their user links in order to maximise global system capacity. At such frequencies, according to the local meteorological conditions, the tropospheric attenuation can become prohibitive. Ever since, in order to provide an acceptable link availability, UpLink Power Control (ULPC) and Adaptive Coding and Modulation (ACM) techniques are not sufficient and gateway diversity needs to be implemented. Switches of gateway to counteract propagation at Q/V band may be activated from estimation of the attenuation at 50 GHz [START_REF] Jeannin | Smart gateways for terabit/s satellite[END_REF]. This estimation can be carried out from measurements of attenuation from a downlink beacon at 40 GHz and/or at 20 GHz scaled at 50 GHz using frequency scaling techniques.

Frequency scaling is a technique to infer the attenuation at a given frequency from the knowledge of the attenuation at another frequency. It gives more accurate results than prediction methods using solely rainfall rate measurements [START_REF]Propagation data and prediction methods required for the design of Earth-space telecommunication systems[END_REF]. It can be used either to convert Cumulative Distribution Functions (CDF) of attenuation from the lower frequency toward the higher or to synthetize attenuation time series for system design, or to perform real time frequency scaling to trigger ULPC or to switch from a gateway to another one using N+P diversity. At microwave frequencies, among all the tropospheric effects (rain, clouds and gases), rain attenuation is the most difficult to scale because its frequency scaling ratio can exhibit strong variations, in particular during a rain event [START_REF] Sweeney | Hysteresis effects in instantaneous frequency scaling of attenuation on 20 and 30 GHz satellite links[END_REF].

A major part of the frequency scaling models presented in the open literature are empirically derived. This is the case of the model developed by L. Boithias, which is one of the model recommended by Rec. ITU-R P.618-13 [START_REF]Propagation data and prediction methods required for the design of Earth-space telecommunication systems[END_REF], section 2.2.1.3.2, for scaling equiprobable attenuations (CCDF). A statistical model is proposed in section 2.2.1.3.1. Recently, two rain frequency scaling techniques have been proposed, both relying on a physical approach:  Frequency scaling relying on the synthetic storm technique (SST) [START_REF] Luini | Enhancement of the Synthetic Storm Technique for the Prediction of Rain Attenuation Time Series at EHF[END_REF].

 Frequency scaling relying on dual beacon measurements with an inversion of the rain Drop

Size Distribution (DSD) [START_REF] Queyrel | Frequency Scaling Relying on the Estimation of the Drop Size Distribution from Dual Frequency Beacon Measurements[END_REF].

Because empirical models have been derived from measured propagation data, they have a limited frequency range (7 to 50 GHz for Boithias' model). On the contrary, physical models of rain frequency scaling rely on hypotheses on the rain drop size distribution, and should therefore have a wider frequency validity range (k and α coefficients of Rec. ITU-R P.838 are obtained from a Laws-Parsons rain DSD and are accurate up to 1000 GHz [START_REF]Specific attenuation model for rain for use in prediction methods[END_REF]). As a consequence, these physical models are expected to be more accurate for the design of future systems, such as W band for which no propagation data are available so far.

This paper aims to introduce a methodology to objectively compare frequency scaling models using error metrics. The models are briefly recalled in Section 2. Then, in Section 3, the full process to obtain a clean database of rain attenuation is explained. The error metrics used to evaluate the frequency scaling models are presented in Section 4. Their performances are finally compared in Section 5 by comparisons with experimental attenuation time series collected at Spino d'Adda (Italy, 7 years, frequencies 18.7 GHz, 39.6 GHz and 49.5 GHz).

Frequency scaling models

In this paper, three frequency scaling models are compared:

 The model proposed by L. Boithias, currently used by Rec. ITU-R P.618-13 (section 2.2.1.3.2) for frequency scaling of equiprobable attenuations. It was derived empirically using a fitting procedure on measurements at various frequencies.

 A frequency scaling technique based on the Synthetic Storm Technique (SST) for slant paths [START_REF] Matricciani | Physical-mathematical model of the dynamics of rain attenuation based on rain rate time series and a two-layer vertical structure of precipitation[END_REF], or rather on an updated version, the so-called E-SST [START_REF] Luini | Enhancement of the Synthetic Storm Technique for the Prediction of Rain Attenuation Time Series at EHF[END_REF]. The SST is a physical model to reconstruct rainfall rate fields along the direction of the wind. Using a power relationship 𝜸 𝑹 = 𝒌. 𝑹 𝜶 to convert rain rate 𝑹 [mm/h] into specific rain attenuation 𝜸 𝑹 [dB/km], it is possible to calculate a frequency scaling ratio to scale the rain attenuation measurements from one frequency to the other.

 A model using measurements at two frequencies, based on an inversion of the drop size distribution of the rain along the propagation path [START_REF] Queyrel | Frequency Scaling Relying on the Estimation of the Drop Size Distribution from Dual Frequency Beacon Measurements[END_REF].

Data and processing

Sources of Data

The data used for this analysis are:

 Propagation data: a database of beacon attenuation, measured in Spino d'Adda, from 1994 to 2000 (7 years), for the three frequencies of the Italsat beacons (18.7, 39.6 and 49.5 GHz), collected and processed by Politecnico di Milano. The time series are sampled at 1 Hz. Beacon attenuation is obtained by measuring the power level of the signal received by a satellite beacon, from which is subtracted a so-called 0 dB level corresponding to the level of power received in clear sky conditions. It is therefore the combination of clouds, rain and tropospheric scintillation.

 Rainfall rate time series from a tipping bucket rain gauge (60s integration time), colocated with the beacon receiver in Spino d'Adda.

 Reanalysis data (ERA5) provided by the European Centre for Medium-Range Weather Forecasts (ECMWF): wind speed at 700 hPa and 0°C isotherm altitude as an input for the SST, and columnar cloud liquid water content to be used to compute the CDF of cloud attenuation. The conversion between columnar cloud liquid water content and cloud attenuation is done using Rec. ITU-R p.840-8, section 3.2 [START_REF]Attenuation due to clouds and fog[END_REF]. These data are interpolated at the location of the beacon receiver from a 0.25° grid in longitude and latitude. The time resolution is 1 hour.

All data are used with a strict sample-by-sample concurrency.

Selection of high quality periods for beacon attenuation

Beacon attenuation measurements are used, after a procedure described in the following section, as an input for rain frequency scaling models (at 18.7 and 39.6 GHz), but also as a reference (at 49.5 GHz). Consequently, these measurements must not be too affected by unwanted signals in order to obtain significant results. The amplitude of tropospheric scintillation can be quite high in beacon measurements [START_REF] Matricciani | 18.7 GHz tropospheric scintillation and simultaneous rain attenuation measured at Spino d'Adda and Darmstadt with Italsat[END_REF]. Issues with the beacon receiver, for instance due to aging of the components, might cause errors in the data as well. This section describes an automated procedure to identify the parts of the database with the highest quality.

It is generally accepted that the power spectral density (PSD) of beacon attenuation is composed of several parts. The PSD for frequencies below 0.02 to 0.1 Hz are the result of rain and cloud attenuation, and display a typical -20 dB/decade slope [START_REF] Matricciani | Physical-mathematical model of the dynamics of rain attenuation with application to power spectrum[END_REF]. The PSD for higher frequencies is ordinarily caused by scintillation, with a flat part up to a cut-off frequency, followed by a -80/3 dB/dec slope [START_REF] Karasawa | Characteristics of fading on low-elevation angle Earth-space paths with concurrent rain attenuation and scintillation[END_REF]. Figure 1 shows examples of PSD of attenuation derived from daily experimental time series. The blue curve corresponds to a rain event with very low noise levels, whereas the red curve corresponds to a day affected by strong scintillation, which happens more often during the summer period [START_REF] Karasawa | Characteristics of fading on low-elevation angle Earth-space paths with concurrent rain attenuation and scintillation[END_REF]. For the classical modelling of scintillation, with a flat spectrum below its cut-off frequency, a measurement of the power spectral density in a narrow band around 0.1 Hz is expected to be representative of the relative level of noise in the measurements. Eventual issues of the receiver causing noise with high frequency spectral components (such as white noise) could be detected as well.

3.3.

Processing of experimental beacon measurements

The conventional method [START_REF] Boulanger | Four Years of Total Attenuation Statistics of Earth-Space Propagation Experiments at Ka-Band in Toulouse[END_REF] to obtain validated rain attenuation (in-excess attenuation) time series consists in several steps. First, the high frequency components of the signal, dominated by tropospheric scintillation or noise, are removed with a low-pass filter (with a typical cut-off frequency of 0.025 Hz). Then, the rain events are marked (often with a visual inspection of beacon attenuation and rain gauge time series), to remove a value linearly interpolated between the attenuations at the beginning and at the end of the rain event. Nevertheless, the identification of rain events is not an easy task, especially on a slant path for which the presence of rain on the beacon receiver (measured punctually by a rain gauge) may not be representative of the presence of rain along the propagation path. In addition, cloud attenuation is likely to increase with rain intensity, which would cause rain attenuation to be overestimated.

To overcome this difficulty, this study uses the statistical method proposed in [START_REF] Luini | Enhancement of the Synthetic Storm Technique for the Prediction of Rain Attenuation Time Series at EHF[END_REF] to remove the contribution of clouds from the beacon attenuation. This method consists in defining the attenuation due to clouds as a function of the beacon attenuation. The function is defined so that the distribution of cloud attenuation matches as well as possible a reference distribution. In [START_REF] Luini | Enhancement of the Synthetic Storm Technique for the Prediction of Rain Attenuation Time Series at EHF[END_REF], the reference cloud distribution is computed from radiosonde measurements (RAOBS). Particularly, the cloud attenuation 𝑨 𝑪 is calculated from the beacon rain+cloud attenuation 𝑨 𝑹𝑪 from: 𝒂, 𝒃, 𝒄, 𝑨 𝑪 𝒎𝒂𝒙 and 𝑨 𝑹𝑪 𝒎𝒂𝒙 : where 𝒂, 𝒃, 𝒄, 𝑨 𝑪 𝒎 and 𝑨 𝑹𝑪 𝒎 are parameters determined by an optimization procedure to fit the equiprobable values of 𝑨 𝑪 and the reference cloud attenuation derived from RAOBS. The rain attenuation 𝑨 𝑹 = 𝑨 𝑹𝑪 -𝑨 𝑪 is then obtained by subtracting the cloud attenuation 𝑨 𝑪 to the beacon attenuation 𝑨 𝑹𝑪 .

𝑨 𝑪 = { 𝑹 𝑹𝑪 ⋅
In the present paper, as ERA5 data are used instead of RAOBS measurements, a new formulation is proposed for 𝑨 𝑪 :

𝑨 𝑪 = 𝐥𝐨𝐠 𝟏𝟎 (𝜷 ⋅ 𝑨 𝑹𝑪 + 𝟏) 𝜷 (𝟑) 𝜷 = 𝒂 ⋅ (𝟏 -𝒆 -𝒃⋅𝑨 𝑹𝑪 ), (𝟒)
where 𝒂 and 𝒃 are obtained by fitting the values of 𝑨 𝑪 exceeded for various time percentages to the reference values of cloud attenuation, exceeded for the same time percentages, derived from ERA5 data. Equations ( 3) and ( 4) have been derived empirically, and gave good results for the site of Spino d'Adda (at 18.7, 39.6 and 49.5 GHz) as can be observed on Figure 2 where the cloud attenuation Complementary CDF (CCDF) is represented from ERA5 data (the reference), from equations ( 1) and ( 2) (methodology of [START_REF] Luini | Enhancement of the Synthetic Storm Technique for the Prediction of Rain Attenuation Time Series at EHF[END_REF]) and from equations ( 3) and (4) (our approach). Similar results have been observed using beacon attenuation collected in Toulouse. Nevertheless, it is not obvious that this particular function will provide a good fit for all datasets, especially if the climatology is different. Using tabulated values of 𝑨 𝑪 and 𝑨 𝑹𝑪 could be a more practical approach for a larger number of sites. It is worthwhile to note that above probabilities of exceedance of 10%, the cloud attenuation computed from ERA5 data is higher than the beacon attenuation, presumably because of the spatial resolution of ERA5 (the columnar cloud liquid water content is averaged over 0.25°x0.25° pixels). Figure 3 shows the rain attenuation CCDF (in-excess attenuation) derived from the experimental measurements once the cloud contribution is removed with equations (3) and (4). [1994][1995][1996][1997][1998][1999][2000], at Spino d'Adda, 39.6 GHz, for clouds from ERA5 data (reference CCDF), from [START_REF] Luini | Enhancement of the Synthetic Storm Technique for the Prediction of Rain Attenuation Time Series at EHF[END_REF] with ( 1) and ( 2), and from ( 3) and ( 4).

Error metrics

Usual metrics used to compare synthetic attenuation time series (i.e. derived from frequency scaling algorithms) with experimental ones are the following: the root mean square error (RMSE), the standard deviation, and the mean error (or bias). One common approach is to compute these metrics on the full set of samples of the database, and to give the result as a single number [START_REF] Luini | Enhancement of the Synthetic Storm Technique for the Prediction of Rain Attenuation Time Series at EHF[END_REF], [START_REF] Paiement | Instantaneous Frequency Scaling for Efficient Mitigation of Uplink Atmospheric Attenuation at EHF[END_REF].

While this approach has the advantage of being simple, it is problematic for several reasons:

 Due to the lognormal distribution of rain attenuation [START_REF] Boulanger | A Rain Attenuation Time-Series Synthesizer Based on a Dirac and Lognormal Distribution[END_REF], most samples refer to low attenuation values: typically, half of the identified rain events in the database of Spino d'Adda have peak attenuations below 5 dB at 49.5 GHz. Consequently, usual single-figure error metrics will be mainly determined by the performance of frequency scaling models at low attenuations, while higher attenuations are of grand interest for fade mitigation techniques.

 It is expected that the mean error of the frequency scaling models depends on the attenuation (the frequency scaling model currently proposed by the ITU predicts that the frequency scaling ratio decreases with increasing attenuation [START_REF]Propagation data and prediction methods required for the design of Earth-space telecommunication systems[END_REF]). Averaging the mean error on the entire database could cancel out attenuation dependent biases of opposite signs.

Giving the RMSE and bias for bins of attenuation solves these issues. In this paper, error figures use equal frequency binning (all bins contain the same number of samples), so that all points have the same statistical significance.

Results

The three frequency scaling models previously mentioned have been applied on the database of Italsat beacon attenuation, recorded in Spino d'Adda between 1994 and 2000 by Politecnico di Milano. The procedure to remove noisy events described in section 3.2 significantly reduced the RMSE of all frequency scaling models (specifically, Figure 4 shows improvements in terms of RMSE for the model of the ITU). All later results use a database with noisy events excluded. A rain event for the three beacon frequencies is shown in Figure 5, with the three frequency scaling models previously mentioned.

The prediction frequency is 49.5 GHz, the models from ITU and SST use measurements at 39.6 GHz (frequency scaling techniques give better results when the frequencies are close). The bifrequency scaling technique uses measurements at 18.7 GHz and 39.6 GHz. Figure 6 gives the performance of the three frequency scaling models on the full database. In terms of RMSE (Figure 6, left), the three models give similar results. The bi-frequency model presents a slightly lower RMSE than the model of Rec. ITU-R P.618-13 for attenuations above 5 dB, but it is biased: it overestimates attenuations lower than 5 dB, and underestimates attenuations above 10 dB (Figure 6, right). The SST is the less biased, but its RMSE is slightly higher. An analysis showed that high values of error came from convective events, for which the rain rate measured by the rain gauge was low when compared to the level of attenuation.

Conclusion

Three frequency scaling models of rain attenuation have been compared on the propagation database of Italsat, in Spino d'Adda (7 years). To this aim, a new formulation of the statistical method to derive the attenuation due to rain from satellite beacon measurements has been proposed. A procedure to automatically detect noisy events has also been implemented, which results in a cleaner database with reduced errors for the frequency scaling models. The two physical frequency scaling models (one using the Synthetic Storm Technique, the other based on an inversion of the rain drop size distribution) show promising results when compared with the current empirical frequency scaling model of the Recommendation ITU-R P.618-13. Their physical basis place them as better candidates for higher frequency bands for which no measurements are currently available, though their validity to other climatological areas has still to be demonstrated.

Figure 1 :

 1 Figure 1: Power spectral density of beacon attenuation of 2 distinct days, at Spino d'Adda, 39.6 GHz.

Figure 2 :

 2 Figure 2: CCDF of attenuation over 7 years(1994)(1995)(1996)(1997)(1998)(1999)(2000), at Spino d'Adda, 39.6 GHz, for clouds from ERA5 data (reference CCDF), from[START_REF] Luini | Enhancement of the Synthetic Storm Technique for the Prediction of Rain Attenuation Time Series at EHF[END_REF] with (1) and (2), and from (3) and (4).

Figure 3 :

 3 Figure 3: CCDF of rain attenuation at 39.6 GHz derived from beacon measurements once (3) and (4) are applied.

Figure 4 :

 4 Figure 4: Binned root mean square error of the frequency scaling technique of Rec. ITU-R P.618-13 for (blue curve) the full database of rain events and (orange curve) a reduced set of events selected using an analysis based on the power spectral density of the rain events.

Figure 5 :

 5 Figure 5: Example of a rain event in Spino d'Adda (11th of June, 1995), with the rain attenuation at the three frequencies, and attenuation from frequency scaling models.

Figure 6 :

 6 Figure 6: Binned error metrics for frequency scaling models at 49.5 GHz, Italsat, Spino d'Adda. Left: RMSE. Right: mean error.

  𝑨 𝑹𝑪 , 𝑨 𝑹𝑪 ≤ 𝑨 𝑹𝑪

		𝒎𝒂𝒙
	𝑨 𝑪 𝒎 ,	𝑨 𝑹𝑪 > 𝑨 𝑹𝑪 𝒎𝒂𝒙	(𝟏)
	𝑹 𝑹𝑪 = 𝒂 𝒆𝒙𝒑(-𝒃𝑨 𝑹𝑪 ) + (𝟏 -𝒂)𝒆𝒙𝒑(-𝒄𝑨 𝑹𝑪 ),	(𝟐)

Acknowledgements

The authors would like to acknowledge the Agenzia Spaziale Italiana (A.S.I.) that supported the ITALSAT propagation experiment and CSTS-CNR for data pre-processing.