Francois Sanson
email: francois.sanson@inria.fr

Charles Bertorello

Jean-Marc Bouilly

Pietro Marco Congedo

Breakup prediction under uncertainty: application to Upper Stage controlled reentries from GTO orbit

Keywords: Space Object reentry, Uncertainty Quantification, Breakup Prediction

come

are extremely challenging to predict and subject to uncertainties. In this work, we present an original model for robustly predicting the breakup of a reentering space object. This model is composed of a set of individual solvers that are coupled together such as each solver resolves a specific aspect of this multiphysics problem. This paper deals with two levels of uncertainties. The first level is the stochastic modelling of the breakup while the second level is the statistical characterization of the model input uncertainties. The framework provides robust estimates of the quantities of interest and quantitative sensitivity analysis. The objective is twofold: first to compute a robust estimate of the breakup distribution and secondly to identify the main uncertainties in the quantities of interest.

Due to the significant computational cost, we use an efficient framework par-

Introduction

Since the beginning of the conquest of space, the number of space rocket launches never stopped increasing. As a result, crucial orbits such as the geostationary orbit, are saturated with end-of-life objects usually referred to as space debris. Upper Stages (US) orbiting in GTO are good examples of threat for active satellites as they cross several essential orbits. According to the French legislation, such debris must be deorbited back on Earth with minimal risk for human assets [START_REF] Dias | Derivation of the french space operation act requirements in the specifications of the future european launcher ariane 6[END_REF]. Assessing the risk implies predicting US reentries robustly.

Unfortunately, such predictions are incredibly challenging. When they renter the atmosphere, US go through intense heating and stress due to aerodynamic effects leading to one or several breakups into several fragments with potentially unknown shapes and material characteristics. Even the breakup time is complex to predict and requires the resolution of hypersonic flows on complex geometries coupled with finite element simulations to estimate the thermal and mechanic response of the US. Such accurate modelling is rarely available, and space object reentry predictors use further simplifying assumptions permitting to obtain predictions within the specific time constraint.

There are two major classes of space object reentry predictors developed in the literature. i) The Object-oriented models use predefined fragments collected inside a parent object such that, at a given altitude chosen by the user, the child fragments are released, and their trajectories are computed. For each debris, an ablation module calculates the mass loss rate. Examples of such pieces of software are for instance DEBRISK developed at Centre d'Etude Spa-tiale (CNES) [START_REF] Omaly | a Tool for Re-Entry Risk Analysis[END_REF] coupled with a risk estimation tool [START_REF] Lazare | Electra launch and re-entry safety analysis tool[END_REF]. Similarly, ORSAT [START_REF] Dobarco-Otero | The object reentry survival analysis tool (orsat)-version 6.0 and its application to spacecraft entry[END_REF] is an object-oriented spacecraft that has been used to predict the reentry of UARS spacecraft [START_REF] Rochelle | Analysis of reentry survivability of uars spacecraft[END_REF]. Similar approaches are implemented in DRAMA [START_REF] Martin | A debris risk assessment tool supporting mitigation guidelines[END_REF] from the European Space Agency (ESA) and DRAPS [START_REF] Ziniu | Space debris reentry analysis methods and tools[END_REF] that uses a Monte Carlo approach to account for the uncertainties in the breakup prediction and the debris survivability. In all those solvers, the breakup time has to be chosen a priori by the user and is necessarily deterministic. The object oriented predictors have a limited computational cost but having a fixed, a priori breakup altitude introduces a significant bias in the ground footprint estimate. ii) More complex models are implemented in vehicle oriented software such as SCARAB [START_REF] Koppenwallner | Scarab -a Multi-Disciplinary Code for Destruction Analysis of Space-Craft during Re-Entry[END_REF], assuming a thermally driven breakup. Such approaches do not model fragmentation as a single catastrophic event but as a sequence of fragment releases where the shape of the fragments are estimated based on the flight conditions and the object structure. Additionally, SCARAB uses an extensive material library allowing to model many materials used for aerospace applications [START_REF] Koppenwallner | Scarab -a Multi-Disciplinary Code for Destruction Analysis of Space-Craft during Re-Entry[END_REF].

The software has been successfully applied to predict ATV and Ariane 5 EPC reentries [START_REF] Koppenwallner | Scarab -a Multi-Disciplinary Code for Destruction Analysis of Space-Craft during Re-Entry[END_REF][START_REF] Lips | Spacecraft destruction during re-entry latest results and development of the scarab software system[END_REF] or the demisabibility of propellant tanks [START_REF] Lips | About the demisability of propellant tanks during atmospheric re-entry from leo[END_REF]. Recently, Onera developed a spacecraft oriented software [START_REF] Prvereaud | Contribution la modlisation de la rentre atmosphrique des dbris spatiaux[END_REF], called MUSIC/FAST, that include the modelling of debris aerodynamic interaction [START_REF] Prevereaud | Debris Aerodynamic Interaction and its Effect on Reentry Risk Assessment[END_REF]. Despite significant efforts from the community to improve breakup models, it is acknowledged that those simplifying models include mismodelling errors [START_REF] Annaloro | Space debris atmospheric entry prediction with spacecraft-oriented tools[END_REF] and should include a quantification of the uncertainties to increase the robustness of the predictions. Considerable efforts have been made to include uncertainty quantification tools in the reentry predictors and orbital propagation [START_REF] Morand | CNES activities on polynomial chaos expansion for uncertainty propagation[END_REF][START_REF] Lips | A comparison of commonly used re-entry analysis tools[END_REF][START_REF] Tardioli | Comparison of non-intrusive approaches to uncertainty propagation in orbital mechanics[END_REF][START_REF] Vasile | Analysis of the de-orbiting and re-entry of space objects with high area to mass ratio[END_REF][START_REF] Mehta | Sensitivity analysis and probabilistic re-entry modeling for debris using high dimensional model representation based uncertainty treatment[END_REF][START_REF] Reyhanoglu | Estimation of debris dispersion due to a space vehicle breakup during reentry[END_REF]. In [START_REF] Frank | A probabilistic paradigm for spacecraft random reentry disassembly[END_REF], a probabilistic paradigm for modelling breakup is proposed. Instead of using a single breakup time, the authors propose to use a predefined breakup distribution. Monte Carlo simulations are then run to compute the ground impact zone. The main drawback of Monte Carlo methods is their computational cost and limited precision when not considering enough samples. For instance, in [START_REF] Frank | A probabilistic paradigm for spacecraft random reentry disassembly[END_REF], only 2000 points were used. In [START_REF] Morand | CNES activities on polynomial chaos expansion for uncertainty propagation[END_REF][START_REF] Lips | A comparison of commonly used re-entry analysis tools[END_REF][START_REF] Fritsche | Uncertainty quantification for re-entry survivability prediction[END_REF], a Polynomial Chaos approach is used to propagate model parameters uncertain-ties through the solver. In their work, the use of a polynomial surrogate model allows to cut down the computational cost of propagating uncertainties, compared to Monte Carlo approaches. However, they modelled a limited number of uncertainties, since getting a good polynomial approximation often requires a large number of runs to obtain reliable results especially in high dimensions [START_REF] Le Maître | Spectral Methods for Uncertainty Quantification, Scientific Computation[END_REF].

In [START_REF] Morio | Estimation of a launch vehicle stage fallout zone with parametric and non-parametric importance sampling algorithms in presence of uncertain input distributions[END_REF][START_REF] Morio | Missile target accuracy estimation with importance splitting[END_REF], efficient statistical methods are proposed to include uncertainties in the evaluation of fallout regions.

In this work, we propose a new breakup predictor for space object reentry with two original features: i) a probabilistic modelling of breakup event to account for the modelling complexity of breakup; ii) the statistical modelling of a large number of uncertain parameters. The main advantages of this model are the robustness of the prediction as it includes the breakup modelling uncertainties and parameters uncertainties. In particular, we consider two level of uncertainty modelling. The first one comes from the probabilistic breakup model: even at nominal input values, the output (flight conditions at breakup) of our model is stochastic. Additionally, we also include a second uncertainty modelling coming from the uncertainties in the model input parameters. In contrast to earlier work, we take into account a large number of uncertain parameters propagated for a low computational cost using surrogate models based on the framework developed in [START_REF] Sanson | Uncertainty Propagation Framework for Systems of Solvers[END_REF].

The breakup predictor proposed in this work is an assembly of interconnected solvers as represented in Fig. 1 and presented in Section 2. Due to a large number of uncertain variables, classical uncertainty propagation methods such as Polynomial Chaos [START_REF] Xiu | The Wiener-Askey Polynomial Chaos for Stochastic Differential Equations[END_REF] methods or Gaussian Processes [START_REF] Chen | Uncertainty propagation using infinite mixture of gaussian processes and variational bayesian inference[END_REF][START_REF] Rasmussen | Gaussian processes for machine learning[END_REF][START_REF] Bilionis | Multi-output separable gaussian process: towards an efficient, fully bayesian paradigm for uncertainty quantification[END_REF] are not optimal. This difficulty is mitigated here by using a framework proposed in [START_REF] Sanson | Uncertainty Propagation Framework for Systems of Solvers[END_REF] and briefly introduced in Section 3. Finally, Section 4 provides qualitative and quantitative analysis about the importance of input uncertainties through global sensitivity analysis. We study the evolution of the uncertainties at different points of the trajectory before breakup. Finally, we derive a breakup distribution of an Upper Stage that can be used by object-oriented software to assess the survivability of object fragments.

Numerical tools

In the case of a controlled reentry, the end-of-life object receives a final deorbiting boost so that it quickly falls back on Earth. While falling back, the object goes through intense heating and stress from aerodynamic effects such that at some point, the object breaks up. The phenomena leading to its breakup are modelled here by four solvers coupled as shown in Fig. 1 The second class is constituted by the intermediate variables, which are at the same time outputs of an upstream solver and inputs of a downstream solver.

As observed in Figure 1, the computation starts with the deorbiting solver.

The deorbiting solver takes as inputs the initial orbit and the deorbiting manoeuvre characteristics. It computes the trajectory of the US until it reaches an arbitrary altitude of 120 km. The deorbiting solver returns the position and the velocity of the object at 120km of altitude denoted as X 120km and V 120km , which are intermediate variables. Meanwhile, by using the aerodynamic solver, the aerodynamic coefficient table of the object (which are intermediate variables) are generated for all possible object orientations, altitudes and Mach observed during the reentry from 120km to breakup.

Using the deorbiting solver outputs at an arbitrary altitude of 120 km (X 120km and V 120km ), the aerodynamic coefficient table, the breakup model parameters, the atmosphere model parameters and the material characteristics, the trajectory/thermal solver is launched to predict the breakup time, position and velocity of the object at breakup. Finally, using the trajectory, i.e. X(t) and V (t), the probabilistic breakup solver predict the breakup probabilistic distribution.

In the following, we provide a detailed description of each solver.

Deorbiting solver

The deorbiting solver is an orbit propagator that computes trajectories after a deorbiting boost is released. It takes as input the orbital elements of the initial orbit and the deorbiting characteristics previously defined. At t boost , a retro boost is fired with amplitude D boost and an orientation characterized by two angles, i.e. α boost and β boost . The new orbit characteristics are evaluated, and the orbit propagated until the object reaches 120 km of altitude.

Aerodynamics Solver

The module computes aerodynamic forces and moments for a given object at given flow conditions and orientation. The aerodynamic forces F aero F aero F aero and moments M aero M aero M aero are defined as:

F aero F aero F aero = S (p p p + τ τ τ )dS (1)
M aero M aero M aero = S (r r r × p p p + r r r × τ τ τ )dS (2)
where p p p is the local pressure, normal to the surface and τ τ τ is the shear stress, tangent to the surface, integrated over the object surface S. The solver uses a local panel method: the object is represented by a set of elementary surface panels for which p p p and τ τ τ are computed. The force and moment are obtained by summing all the contributions from each panel.

The formulation of the pressure and shear stress depends on the flow regime.

Three flow regimes are observed during a reentry: the free molecular flow, the rarefied flow and the continuous flow.

The free molecular flow exists at high altitudes and low densities, where the molecular mean free path is comparable to the object size. In this regime, analytical expressions of the pressure and shear stress hold using [START_REF] Gombosi | Gaskinetic theory[END_REF] (eq. 7.88-7.89) and injected in Eqs. 1 and 2.

In the continuous flow regime, the mean free path is much smaller than the object and the flow can be seen as continuous such as Navier Stokes equations are valid. In this case, the modified Newton law holds [START_REF] Anderson | Hypersonic and High Temperature Gas Dynamics[END_REF] for the pressure forces while the shear stress is set to zero. The modified Newton law is valid for hypersonic flow (Mach> 4) as usually encountered during the early phase of reentry (before breakup).

The rarefied flow is the transition between the free molecular flow and the continuous flow. Computing the resulting forces usually requires to run expensive Direct Simulation Monte Carlo (DSMC). In this work, we use bridging functions that define the pressure and shearing stress as a weighted sum of their expression in the free molecular regime and the continuous regime.

The aerodynamic contributions (force and moment) are not directly used by the trajectory/thermal solver but used as the aerodynamic coefficients. In a six degree-of-freedom setting, there are several aerodynamic coefficients to consider:

the drag coefficient C D , the lift coefficient C L , the side force coefficient C Y , the rolling moment coefficient C l , the pitching moment coefficient C m and the yawing moment C n . The relationship between the aerodynamic coefficients and the corresponding aerodynamic force is given by:

F * = p dyn C * S ref (3)
with C * the aerodynamic coefficient, p dyn the dynamic pressure and S ref the object reference area chosen by the user.

The relationship between the aerodynamic coefficients and the corresponding aerodynamic moment is given by:

M * = p dyn C * S ref L ref (4)
with L ref the object reference length.

The aerodynamic solver is used to build a complete aerodynamic table containing the aerodynamic coefficients for any flight condition encountered during the reentry. In particular, the constructed table depends on the object orientation, flow Mach number and Knudsen.

Trajectory/thermal solver

The trajectory/thermal solver is composed of two distinct modules: a trajectory module and a thermal module. It computes the following quantities along the trajectory of the object: the position X X X(t), the velocity V V V (t) = Ẋ Ẋ Ẋ(t), the temperature T (t) and mass m(t). These quantities are a function of time from 120 km until the object reaches an altitude of 30 km. The 30 km altitude is chosen low enough so that the object breaks up before reaching this altitude. In our model, we consider two parts of the object structure: the internal structure and the shell. The internal structure is composed of various components (with a global mass of m int ): internal tanks, pipes and many others, whereas the shell (of mass m shell ) is the external structure of the object. The internal components are not included in the thermal model since we assume that the shell of the object absorbs all the heat flux and that the conductive flux between the internal components of the shell is neglected. This approach is classical and used for instance in object-oriented models where the shell would be a parent object, and the components would be the child fragments.

Trajectory module

The equations governing the trajectory are the following:

d 2 X X X dt 2 = A e + A c + g + F aero F aero F aero m shell + m int (5)
I • dΩ Ω Ω dt = M aero M aero M aero -Ω Ω Ω ∧ (I • Ω Ω Ω) ( 6
)
dT dt c p m shell = S ref q aero (X X X, V V V ) - dm shell dt H f us ( 7
)
where X is the position of the object in the Earth centered inertial frame, g is the gravity constant, A e the centrifugal acceleration and A c the Coriolis acceleration. Note that the term F aero F aero F aero represents the aerodynamic forces precomputed by the aerodynamic solver (see Eq. 1) computed using Eq. 3 and the aerodynamic table computed with the aerodynamic solver. Moreover, Ω Ω Ω is the instantaneous rotation vector expressed in the engine inertial frame; I is the inertia matrix associated with the object in the engine inertial frame;

the term M aero M aero M aero represents the moment induced by the aerodynamic forces and precomputed by the aerodynamic solver (see Eq. 4).

The last equation is the energy conservation equation for the object shell with T the temperature of the object shell, c p its heat capacity, q aero the heat flux received from the atmosphere, S ref is the reference area and H f us the fusion enthalpy. Note that Eq. 7 is coupled with Eqs. 6 and 5 through q aero and the mass of the object shell. Equation 7requires the computation of the heat flux q aero , which is done employing the software ADRYANS developed by ArianeGroup and presented in the next Section.

Thermal module

ADRYANS [START_REF] Bertorello | Adryans v5.0 a 1-dimension object oriented survivability tool[END_REF] is a thermal model used along a trajectory to compute the heat flux q aero (used in Eq. 7). The chosen model assumes uniform object temperature and neglects chemical reactions at the surface of the object such as metal oxidation or molecular recombination. Hence, q aero is the sum of two contributions:

q aero = q conv + q rad , (8)
where q conv is the convective heat flux and q rad is the radiated heat flux. Note that the conductive heat flux from the shell to the internal structure is neglected.

To recover the averaged convective heat flux over the US, we simplify the geometry of the object into a cylinder in random tumbling. Then ADRYANS uses correlation formulae close to [START_REF] Klett | Drag coefficients and heating ratios for right circular cylinders in free-molecular and continuum flow from mach 10 to 30[END_REF] to relate the convective heat flux received by the US with the heat flux at the stagnation point of a flat plate for which analytical formulations exists. For instance, in the free molecular flow regime, the convective heat flux at the stagnation point of a perpendicular flat plate is [START_REF] Klett | Drag coefficients and heating ratios for right circular cylinders in free-molecular and continuum flow from mach 10 to 30[END_REF] qF

M sp = αρ ∞ V 3 ∞ 2 ( 9
)
where α is the accommodation coefficient, ρ ∞ and V ∞ are the free stream flow density and velocity.

In continuous regime, at the stagnation point of a perpendicular flat plate, the Detra-Hidalgo correlation gives [START_REF] Detra | Generalized heat transfer formulas and graphs for nose cone reentry into the atmosphere[END_REF]:

qcont sp = C √ R N V ∞ V c 3.15 ρ ∞ ρ s 0.5 (10)
with ρ s being the atmosphere density at sea level,V c the circular orbit velocity and R N the object nose radius.

For the transitional regime the Matting bridging function is used [START_REF] Matting | Approximate bridging relations in the transitional regime between continuum and free-molecule flows[END_REF]:

qrar = qcont 1 -exp - qF M qcont . (11)
The radiative heat flux is computed using Stefan-Boltzmann law:

qrad = σ (T 4 -T 4 env ) ( 12
)
where σ is the Stefan Boltzmann constant, is the material emissivity, T is the temperature and T env the environment temperature set to 4K.

Breakup module

The trajectory/thermal solvers computes the object position, velocity, mass and its temperature (X(t), V (t), m shell (t), T (t)). The breakup module directly uses those quantities to compute the object flight conditions (position X X X breakup and velocity V V V breakup ) at the time when the breakup occurs. As in [START_REF] Frank | A probabilistic paradigm for spacecraft random reentry disassembly[END_REF], the prediction is probabilistic. For a given deterministic object trajectory (X(t),

V (t), m shell (t), T (t)),
we derive a distribution of breakup conditions. This assumption constitutes the first level of uncertainty modelling. Note that later in this work, we also consider a stochastic representation of the solver inputs, yielding a second level of uncertainty modelling. Here, the probabilistic breakup time is sampled from an underlying distribution to represent the uncertainties in the breakup process.

We define the breakup time distribution T f rag defined as:

T f rag = U(t init , t end ), (13)
where U denotes a uniform distribution and t init and t end are the lower and upper bounds of the time interval between which the breakup is bound to happen. In this work, thermal criteria are used to predict t init and t end . As the US flies through denser and denser layers of the atmosphere, its temperature increases and its strength decreases until it reaches the fusion temperature. Once weakened, the shell protecting the internal structure can be broken long before reaching fusion temperature due to the aerodynamic forces and the degraded structural performance due to the temperature increase [START_REF] Finzi | Simulation of the ariane 5 epc reentry with the fragmentation tool suite[END_REF]. Therefore, we assume that breakup can begin as early as when the shell temperature reaches

T f rag with T f rag < T f us .
Mathematically, t init is defined as:

T (t init ) = T f rag , (14)
where T f rag is a model parameter considered as uncertain (see uncertainty definitions in Section 4.2). It depends on the material and also on the atmospheric conditions during the reentry.

Once most of the shell has been ablated, it is safe to say that the breakup has already occurred. Hence, we define t end , the latest time for which breakup can occur, as the time where a certain mass proportion of the shell p abl has been ablated. As for T f rag , p abl is an uncertain model parameter. For a given p abl , t end is obtained by solving the following equation:

m shell (t end ) = m shell (0)p abl , (15)
where m shell (0) is the initial mass of the US shell.

Recall that X X X(t) and V V V (t) are the position and velocity of the US at time t, the breakup position, defined as X X X breakup , is given by:

X X X breakup = X X X(T f rag ) (16)
and the breakup velocity, defined as V breakup , by

V V V breakup = V V V (T f rag ). ( 17
)
Note that after the breakup, the user-defined fragments should be released with the initial conditions (X X X breakup , V V V breakup ). Since T f rag is a random variable, X X X breakup and V V V breakup are random variables. To make a prediction, one needs to generate a sample denoted t f rag from T f rag and the predicted breakup conditions are X X X(t f rag ) and V V V (t f rag ).

Discussion of the breakup distribution. The breakup distribution presented is voluntarily simple, but it can be modified to include more a priori information about the breakup or even to model several breakups. For instance, instead of using a uniform distribution that represents our complete lack of knowledge about the breakup time between t init and t end , one could use more informative distributions such as the beta distribution with suited a, b parameters based on a priori structure analysis. Moreover, for simplicity, the model assumes that there is only one breakup time but, as numerous in-flight experiments have shown [START_REF] Beck | Cleansat building block 12-Demisable joints[END_REF], there are several breakups and fragment releases along the trajectory.

This assumption however can easily be removed by considering distinct release times for each fragments, that is to define T i f rag = U(t init , t end ) associated with fragment i. The joint distribution of (T i f rag ) can be chosen independent (each fragment is released independently from the others) but if a chronology of the fragment release is available such that we know that

T 1 f rag < T 2 f rag < . . . T n f rag , (18)
then a more informative breakup distribution can be sampled. For instance, assuming that the T i f rag are identically distributed, then the sampling process consists in sampling n times from T 0 f rag to constitute t 1 f rag , . . . , t n f rag labelled such as

t 1 f rag < • • • < t n f rag .
For the application presented in this work, we use the simplest formulation as in Eq. 13.

Uncertainty Quantification tools

In addition to modelling breakup as a random event, one objective of this paper is to characterize the uncertainties and to propagate them through the system of solvers presented in the previous section. Uncertainty Propagation (UP) is often a computationally challenging task because it requires a significant number of runs of the SoS. Several methods exist in the literature to accurately propagate uncertainties. Sampling methods, for instance, are the most straightforward strategies but require a large number of runs and are often out of reach due to their prohibitive computational load. Another classical approach in UP consists in building a cheap surrogate model of the system that accurately emulates the system outputs. Sampling methods can then be used with the cheap surrogate in place of the actual SoS at negligible computational cost. For this method to be efficient, the surrogate needs to accurately emulate the real system output. In general, the surrogate learns the mapping between the inputs and the output using an (ideally) limited number of training points. Several methods for building surrogates exist, such for example Polynomial Chaos Expansions (PCE) or Gaussian Processes (GP). For most techniques, the increase in the number of input dimension usually degrades their predictive ability. This effect is known as the curse of dimensionality [START_REF] Bellman | Applied dynamic programming[END_REF]: more training points are necessary to maintain a certain level of accuracy when the input dimension increases. Our case is exceptionally challenging for state-of-the-art UP methods as it features a large number of dimensions.

Moreover, our problem features several solvers while most surrogate model building strategies are designed for individual solvers. In this work, two solvers (the deorbiting solver and the breakup solver) are composed together to predict the breakup conditions. In the case of two coupled solvers, they generally yield a highly non-linear mapping between the inputs and the outputs that can be difficult to learn.

In this work, we use the framework developed in [START_REF] Sanson | Uncertainty Propagation Framework for Systems of Solvers[END_REF] that tackles both the problem of high input dimension and multiple solvers. It efficiently builds a surrogate model of each constituting solver and uses the composition of the surrogate models to make predictions. The set of surrogate models built upon the system solvers is called a System of Gaussian Processes (SoGP). The interested reader can find details on how to build and train such surrogate models in [START_REF] Sanson | Uncertainty Propagation Framework for Systems of Solvers[END_REF]. The rest of this section illustrates the central concepts for training making predictions with SoGPs.

Predictions with SoGP

Predictions with a single GP

The elementary components of a SoGP are single Gaussian Processes, denoted here as GPs. GPs are stochastic processes used as a probabilistic approximation of functions f : R n -→ R m . GPs have been extensively used for uncertainty propagation [START_REF] Chen | Uncertainty propagation using infinite mixture of gaussian processes and variational bayesian inference[END_REF] or sensitivity analysis [START_REF] Oakley | Probabilistic sensitivity analysis of complex models: a bayesian approach[END_REF]. A complete introduction to GPS can be found in [START_REF] Rasmussen | Gaussian processes for machine learning[END_REF][START_REF] Cressie | The origins of kriging[END_REF]. GPs are entirely characterized by their mean and covariance functions. A classical covariance function is the squared exponential covariance function:

k(x, x ) = σ 2 exp - (x -x ) 2 2l 2 (19)
in the case where n = m = 1. Parameters σ and l are hyper-parameters referred as variance and length-scale. We consider the case where the mean function is 0. Conditioned on a noiseless set A of observations x i , f (x i ), the predictive distribution of the GP model of f , at input location x is Gaussian distributed with mean and variance :

µ(x) = k A (x)K -1 A y A , (20)
σ 2 (x) = k(x, x) -k A (x)K -1 A k A (x) T . ( 21
) with k A (x) = (k(x, x 1 ) • • • k(x, x p )) T and K Ai,j = k(x i , x j ).
The choice of covariance functions is critical and is extensively discussed in [START_REF] Stein | Interpolation of spatial data: some theory for kriging[END_REF]. In our problem, the covariance functions are sums of uni-variate squared exponential covariance functions and a joint non-isotropic squared exponential covariance function.

The covariance function hyper-parameters are a priori unknown and are fitted by maximizing the GP likelihood with respect to A [START_REF] Rasmussen | Gaussian processes for machine learning[END_REF]. Practical implementation is based on the python software GPy [START_REF] Gpy | A gaussian process framework in python[END_REF] 3.1.2. Predictions with a multiple GPs SoGPs are obtained by substituting the actual solvers in the SoS by GP models. The composition of Gaussian Processes have been used for machine learning application in [START_REF] Damianou | Deep gaussian processes[END_REF]. In general the resulting predictive distribution is not Gaussian distributed and it has no simple analytical expression except for very particular kernels [START_REF] Girard | Gaussian Processes: Prediction at a Noisy Input and Application to Iterative Multiple-Step Ahead Forecasting of Time-Series[END_REF]. In this work we use the composition of averages of each GP to make predictions. [25]

Training strategies

The main difficulty in training SoGP comes from the generation of training sets for intermediate solvers. While for single solver problems, the input variables are often independent and their distributions well defined, for solvers using other solvers outputs as inputs, their input distributions are often unknown and dependent. As a consequence, the use of low discrepancy methods or Latin Hyper Cube sampling is not straightforward, since they are designed for independent probability distributions [START_REF] Mckay | Comparison of three methods for selecting values of input variables in the analysis of output from a computer code[END_REF]. Instead, a more sensitive approach is to generate a LHS training set on the global input space, propagate it through the SoS and to use propagated samples as a training set for the subsequent solver (see [START_REF] Sanson | Uncertainty Propagation Framework for Systems of Solvers[END_REF] for details). An alternative training sample generation techniques based on clustering can also be used [START_REF] Sanson | Efficient clustering based training set generation for systems of solvers[END_REF]. The two strategies were tested on the SoGP under investigation and yield very similar performance in terms of accuracy, hence we arbitrarily present the result obtained with the first one.

Sensitivity analysis

In order to assess the sensitivity of the quantities of interest with respect to the input uncertainties, we use the Sobol indices. The Sobol indices are classical tools [START_REF]Theorems and examples on high dimensional model representation[END_REF][START_REF] Saltelli | Variance based sensitivity analysis of model output. design and estimator for the total sensitivity index[END_REF] based on the variance decomposition of the output variance based on each input contribution. The Sobol indices are global indices. They are defined as :

V i = V ar(E X∼i [Y |X i ] V ar(Y ) (22)
V T,i = E X∼i [V ar(Y |X ∼i ] V ar(Y ) ( 23
)
where X are the inputs of the system and Y its output. The vector X ∼i corresponds to the vector of all entries of X except X i . The first order Sobol indices V i accounts for the contribution of input X i to the variance of the output neglecting its interactions with other variables. The total Sobol indices V T,i are the contributions of X i to the variance of the output including the interaction with other variables. It holds that V i ≤ V T,i . In this work, only the first order indices are used as, in our case, they were found to be close to the total Sobol indices. Several methods are available to estimate the Sobol indices in particular with Polynomial Chaos [START_REF] Le Maître | Spectral Methods for Uncertainty Quantification, Scientific Computation[END_REF]. In this work, since a surrogate model is available, they are estimated with a Monte Carlo estimator (Ref. from [START_REF] Saltelli | Variance based sensitivity analysis of model output. design and estimator for the total sensitivity index[END_REF][START_REF] Saltelli | Making best use of model evaluations to compute sensitivity indices[END_REF]) at minimal computational cost.

Application to a US controlled reentry from a GTO orbit

The objective of this work is to provide a qualitative and quantitative analysis of the uncertainties before any breakup and to provide flight condition distributions at breakup in the case of a controlled reentry from a GTO orbit.

Those distributions are crucial inputs of survivability tools such as DEBRISK [START_REF] Omaly | a Tool for Re-Entry Risk Analysis[END_REF]. We focus on the altitude, longitude, latitude and velocity at several critical moments in the reentry such as t init t end and breakup time. Additionally, the object velocity, position and temperature uncertainties are studied at different altitudes to understand the evolution of the uncertainties along the reentry trajectory. This section starts with the characterization of the uncertainties (4.1 and 4.2). The use of a surrogate model is then validated in Section 4.3. We present several studies: the evolution of the uncertainties along the trajectory

Definition of the case study

In a GTO controlled reentry, the US is deorbited from its equatorial GTO orbit. The deorbiting manoeuvre consists in a final backward boost fired to slow down the US and make it fall back on Earth. Fig. 2 illustrates the US. It is mostly composed of aluminium and weighs a total of 7000 kg (see Table 1).

The shell susceptible to ablation weighs 4000 km made of aluminium.

Uncertainty definitions

In this study, we consider four primary sources of uncertainties: the material characteristics, the initial flight conditions, the atmosphere model parameters and the breakup model parameters, for a global number of uncertainties of twenty-two [START_REF] Le Maître | Spectral Methods for Uncertainty Quantification, Scientific Computation[END_REF]. This list does not include the intrinsic uncertainty in the breakup model. The list of uncertainties is summarized in table 2.

The material characteristics. The thermal module uses the fusion temperature (T f us ), the fusion enthalpy (H f us ) and the emissivity ( ε) and density (ρ M ).

The material characteristics such as the density, the fusion temperature and fusion enthalpy are directly measured and therefore are modelled as Gaussian variables and reported in table 2. The emissivity uncertainty is more complex to model as it varies along the trajectory. Studies have shown that the reacting flow alters the surface of the aluminium shell during a reentry [START_REF] Prvereaud | Contribution la modlisation de la rentre atmosphrique des dbris spatiaux[END_REF][START_REF] Barka | Influence of oxidation and emissivity for metallic alloys space debris during their atmospheric entry[END_REF]. While the aluminium emissivity usually varies between 0.3 and 0.5 depending on the surface polishing, due to oxidation, it can reach 0.8. [START_REF] Barka | Influence of oxidation and emissivity for metallic alloys space debris during their atmospheric entry[END_REF]. For this reason, we assume a uniform distribution ranging from 0.3 to 0.8 for the emissivity uncertainty.

Initial conditions. The deorbiting manoeuvre conditions are not perfectly known for the following reasons:

• The orbit is not perfectly known, these uncertainties are reflected through independent uniform uncertainties in the orbital elements.

• The deboost amplitude and direction is only partially known (uncertain quantity of propellant, orientation of the object before deorbiting etc.) hence they are modelled as uniform distributions.

• The deorbiting solver does not compute the orientation of the object at 120 km and therefore they are modelled as uncertain inputs of the trajectory/thermal solver. We model the object orientation uncertainties with a uniform distribution in the space of the rotations. Since the object possesses one symmetry axis (Ox) in Fig. 2 it is equivalent to sampling the euler angles (in degrees) according to :

θ ∼ U(-180, 180) (24)
ψ ∼ arccos(2 × u -1) × 180 π where u ∼ U(0, 1) (25)
where θ is the Euler angle with respect to (Oz) and φ the Euler angle with respect to (Oy). Note that considering a uniformly distributed object orientation does not correspond to taking uniformly distributed Euler angles.

Atmosphere model parameters. The atmosphere conditions (temperature, density and pressure) are fluctuating quantities depending on the time of the day, the day of the year, the solar activity or the Earth magnetic activity that are not a priori unknown in our case. The chosen atmosphere model (MSIS00 [START_REF] Picone | NRLMSISE-00 empirical model of the atmosphere: Statistical comparisons and scientific issues[END_REF])

takes as inputs the time of the year in days, the time of the day an hour, the solar activity or the Earth geomagnetic index. We consider all those parameters as random variables. As the time of reentry is a priori unknown, we assume time distributions as uniform distributions over the year. We also consider the solar flux intensity and the Earth magnetic activity index as uniformly distributed random variables which ranges are chosen based on historical measurements [START_REF] Delorme | Caracteristation de l' environement naturel pour un lanceur spatial decollant du centre spatial guyanais (csg)[END_REF] and reported in table 2.

Breakup model parameters.

As mentioned in the previous section, our probabilistic breakup model uses two uncalibrated (i.e. not estimated from experimental data) parameters T f rag and p abl , defining t init and t end . T f rag is the limiting temperature above which the object strength is too weak to withstand aerodynamic forces. Necessarily, T f rag is smaller than T f us = 800K. Besides, experimental studies found that the strength of aluminium greatly reduced as early as 400K [START_REF] Summers | Overview of aluminum alloy mechanical properties during and after fires[END_REF][START_REF] Finzi | Simulation of the ariane 5 epc reentry with the fragmentation tool suite[END_REF]. In this work, T f rag is arbitrary assumed to be a uniform distribution ranging from 400 K to 700 K. p abl is the minimum percentage of remaining mass such as the object is still intact. In this work, the p abl uniformly varies between 50 % and 70 %.

Surrogate model convergence analysis

An accurate surrogate model should be constructed and validated to obtain the distribution of the QoI. In this section, we perform an analysis of the convergence of the SoGP strategy illustrated in this work. In particular, we compare the performance of the SoGP and a single classical GP used as surrogate models. We illustrate the convergence of the US altitude computed at t init . The L2 error is calculated using 500 validation samples and normalised by the output standard deviation. Fig. 3 shows the evolution of the L2 error of the two surrogates models trained with the same LHS training plans.

The results show that the SoGP framework outperforms the single GP by almost an order of magnitude. Using 1000 training samples, the SoGP reaches an error of 0.1 % of the output standard deviation. We observe a similar behaviour for the predictions of other quantities of interest presented in this work.

This low level of error justifies the use of the SoGP as a surrogate model of the system. Monte Carlo sampling technique can then be used with the SoGP to compute the QoI statistics and distribution and carry out sensitivity analysis, at a negligible computational cost.

To highlight the gain of building a surrogate model compared to standard Monte Carlo approaches, Table 3 summarizes the results obtained when computing the mean and the variance of the object altitude at t init , using a fixed computational budget of 600 solver evaluations representing a cost of around one hour of computation on 28 cores IntelXeon at 2.4 Ghz. We compare three different methods. i) The standard Monte Carlo approach that directly uses the 600 samples to estimate the mean and variance. ii) The single GP and iii)

the SoGP use the 600 samples to build a surrogate model and then run the surrogate model with 500 000 samples at no additional computational cost. The 500 000 samples evaluated using the surrogate models are then used to compute the mean and variance of the altitude. Additionally, a Bayesian statistical error estimate [START_REF] Oliphant | A bayesian perspective on estimating mean, variance, and standard-deviation from data[END_REF] for each estimation is included in Table 3. Note that the surrogate model-based approaches are orders of magnitudes more precise than the standard Monte Carlo approach since many more samples (500 000 vs 600) could be used with the surrogate assisted strategies for the same computational budget.

When comparing two surrogate model-based approaches, the best indicator remains the L2 error computed on a validation set (Fig. 3). In the rest of this section, the results presented are obtained using a SoGP trained with 600 Table 3: Estimation of the altitude first two statistical moments at t init using 600 SoS evaluations using Monte Carlo, a single GP trained on 600 samples using evaluated 500 000 times or a SoGP trained on 600 samples using evaluated 500 000 times.

evaluations.

Breakup range prediction

As stated in Section 2.3 (see Eq. 13), t init and t end are crucial times in the breakup prediction.

t init represents the earliest moment at which breakup can occur and defines the lower bound of the time interval where breakup occurs. Physically, at t init , the US reaches a critical temperature for which the US shell is weakened and may not resist aerodynamic forces.

t end is the other bound delimiting the time interval. At t end a significant portion of the US envelope has been ablated, and the chances that the US is still intact are meagre. As a consequence, we assume that breakup has occurred by that time.

t init prediction

In this section, the objective is to compute the uncertainties in altitude, longitude, latitude and velocity at t init .

In Table 4, the principal statistics of the distributions of altitude, longitude, latitude and velocity at t init are reported, namely, the mean, the standard deviation and different quantiles at 0.5 % 2.5% 50 % 97.5 % and 99.5%. Quantiles are cut points in a distribution delimiting intervals of same probability. They provide additional information about the shape of the distribution compared to the mean and standard deviation, in particular, low and high quantiles provide information about the tails of the distribution. The distribution of altitudes is widespread since 99 % of the distribution (between the 0.5 % and 99.5 %) lies between 63km and 81km with a mean and a median around 71km of altitude.

The results show that breakup cannot occur before 81km. This fact is coherent with most reentry models that assume breakup altitudes between 70 km and 80 km [START_REF] Omaly | a Tool for Re-Entry Risk Analysis[END_REF][START_REF] Morand | CNES activities on polynomial chaos expansion for uncertainty propagation[END_REF].

Note however, that due to the nature of GTO orbit, the slope at which the upper stage reenters is high compared to typical reentry of LEO satellites. As a consequence, the US quickly reaches low altitudes before reaching T f rag . Fig. 4 represents the probability density function of the latitude (4a), of the altitude (4b), of the longitude (4c) and of the velocity (4d) at t init . In Fig. 4b, the altitude distribution is skewed toward low altitudes. Due to the exponential increase in atmosphere density when altitude decreases, the received heat flux given in Eq. ( 10), and the object temperature increases more rapidly at low altitudes than high altitudes. As a result, the object temperature increases more rapidly at a lower altitude and is more likely to reach T f rag . Looking at Eq. ( 10), we see that velocity is a direct factor leading to heat production and temperature increase. Hence high-velocity cases lead to high heat flux and high altitudes at t init . Consequently, there is a positive correlation between the velocity and the altitude of t init . The velocity distribution is also highly t end is the latest time at which the object can breakup, therefore the object is very likely to have fragmented before reaching t end . As for t init , table 5 summarizes the distribution characteristics of altitude, longitude, latitude and velocity at t end . The altitude distribution ranges from 51780 to 65750 meters such that the support of the distribution is a bit smaller than for t init . As for the t init case, the altitude pdf represented in Fig. 6b is slightly skewed toward low altitudes. Concerning the velocity, the uncertainties have significantly increased since the standard deviation is 81 m/s at t init and 173 m/s at t end . The shape of the distribution plotted in Fig. 6d is also very different for the two times. At t end , the PDF is much less skewed but still positively correlated to the altitude.

The uncertainties in longitude and latitude are also larger at t end than at t init as shown in Fig. 6a and6c. Contrary to t init , at t end , the longitude and latitude are highly correlated with the velocity and the altitude. High velocity cases are likely to lead to high altitude breakup than slower one. Finally note that the object velocity remains extremely high even at altitudes as low as 50 km.

Figure 7 shows the first order Sobol sensitivity indices for the four quantities of interest at t end , namely altitude, latitude, longitude and velocity. The deorbiting manoeuvring uncertainties are the predominant sources of uncertainties for all quantities of interest. In the case of velocity (and altitude to a lesser extent), the uncertainty in the breakup model parameter p f rag explains 15 % and 2 % of the output variances. As for t init , the uncertainties in the apogee position of the initial orbit influences the longitude at t end . While the range of variation of p f rag is voluntarily chosen very large (between 50 % and 70% of the shell mass) for robustness, its influence is still secondary compared to the deorbiting uncertainties. This result contrasts with the sensitivity results at t init where most of the uncertainties come from the uncalibrated T f rag parameter model. Note that the uncertainties in the material characteristics have a minor influence in the QoIs, even at t end and despite the vast uncertainty range used to model the emissivity changes due to surface oxidation. This result shows As expected, the range breakup conditions covered by our breakup model is widespread from altitudes as high as 81km down to 52km. The sensitivity analysis showed that most of the uncertainties in the estimation of the flight conditions for t init and t end are mostly driven by the poorly calibrated breakup model parameters p abl and T f rag and the deorbiting manoeuvre uncertainties.

With the aim to improve the predictive ability of the model, the model parameters should be calibrated either with experimental results or high fidelity simulations of the breakup. In this section, we consider the evolution of the uncertainties in the position, velocity and thermal properties of the US along the trajectory where the breakup has a low probability of occurring (i.e. before 70km). The objective is to understand better how the uncertainties in the object position evolve with time. We consider several altitudes: 100km, 90km, 85km, 80km and 70km. In each case, we observe the longitude, latitude, velocity, convective heat flux and radiative heat flux.

Longitude. In Fig. 8, the pdf for longitude at the altitudes of 100km, 90km, 85km, 80km and 70 km are represented. While at 100km, the pdf is fairly symmetrical, it becomes larger and more skewed as the object falls Earth. The sensitivity analysis presented in Fig. 9 shows that the longitude exclusively depends on the deorbiting uncertainties (mostly the boost orientation and the boost amplitude). While it is not surprising that the material characteristics do not influence the position of the object before the breakup, it is interesting to observe that the atmosphere uncertainties also have no influence on the longitude at any altitude. The sensitivity indices vary with altitude significantly. The uncertainties in the initial orbit (apogee and perigee positions) have a significant influence at high altitudes, but their influence reduces at lower altitudes. Meanwhile, the impact of the deorbiting boost amplitude increases as the altitude decreases.

Latitude. Fig. 10 illustrates the pdfs of the object latitude for different altitudes.

As for the longitude, the distributions become wider and more skewed as the altitude decreases. Fig. 11 illustrates the sensitivity indices for the different altitudes. As for the longitude, most of the uncertainties come from the unknown deorbiting conditions. On the other hand, the initial orbit uncertainties have a much smaller influence on the latitude than the longitude. Note that further analysis showed that the first order Sobol indices decrease with the altitude, while higher order interactions between the deorbiting boost orientation angles increase.

Velocity. The reentry considered here is a high-speed reentry compared other human-made object reentries. For instance reentries from LEO would be around 6.5 km/s to 7km/s while this one is around 9.5 km/s. Fig. 12 represents the The residual variance comes from higher-order interactions between the boost amplitude and the boost orientation.

Thermal analysis of the uncertainties

The uncertainties in the object temperature have a significant impact on our ability to accurately predict breakup since the breakup model presented The second heat flux contribution of our model (see eq. 8), is the radiated heat flux. The heat radiated by the object is directly dependent on the temperature and the emissivity of the object through Eq. 12. We report the pdf of the heat flux radiated at different altitudes on a log scale in Fig. 17. As for the convective heat flux, it increases as the altitude decreases. Note that it remains one order of magnitude smaller than the convective heat flux at this state of the reentry. At 100km, the radiative heat flux distribution is almost uniform, but it shifts toward a more bell-shaped curve as the altitude decreases.

The sensitivity analysis results summarised in Fig. 16 show that most of the variance is due to the emissivity variance. This influence is strong at 100km, and it diminishes with the altitude to reach less than 40 % of the variance at 70 km. This observation is coherent with the uniform-like shape of the distribution above 80 km since the emissivity, modelled as a uniform distribution, is directly proportional to the heat flux (eq. 12). The influence of the deorbiting manoeu- total variance. The rest is due to interactions higher-order interactions between the deorbiting boost orientation uncertainties α boost and β boost .

Figure 19 shows the evolution of the object temperature distribution at different altitudes. The uncertainty increases as the object reaches lower altitudes.

In fact, at 100 km of altitude, 95 % of the probability density lies between 307 and 315 K whereas at 80 km it ranges between 518 K and 677 K. Note that at this altitude, the distribution is skewed toward high temperatures.

The sensitivity analysis presented in Fig. 18 also changes with altitude. For high altitudes, the time of reentry is the largest contributor to variance with more than 85 % of the variance explained at 100km but its influence vanishes to only a few percents of the variance at 70 km. On the other hand, the uncertainties in the deorbiting boost become more and more important as the object moves along the trajectory. Note the minimal influence of the initial orbit uncertainties and in particular the perigee position. The uncertainty in the object emissivity does not influence temperature for this portion of the tra- jectory although it dramatically affects the radiated heat flux. This observation is because most of the temperature increase is due to the convective heat flux that is one order of magnitude higher than the radiated heat flux.

by i) is denoted as (X(t m , ξ i )) m∈[0,n] for the position and (V (t m , ξ i )) m∈[0,n] for the velocity. The computational and memory cost of generating that many trajectories is prohibitive since (X(t m , ξ i )) m∈[0,n] and (V (t m , ξ i )) m∈[0,n] can be very large matrices, depending on the time sampling rate. In our case several millions trajectories are necessary to compute the Sobol indices. To overcome this difficulty, the trajectories are represented on a reduced basis obtained with a Karhunen Loeve Expension (KLE) [START_REF] Betz | Numerical methods for the discretization of random fields by means of the karhunenlove expansion[END_REF]. The decomposition expresses the large matrix (X(t m , ξ i )) i,m in a reduced basis such that :

(X(t m , ξ i )) i,m N k=1 α k (ξ i )φ k (t m ) (27)
where the φ k are the reduced basis functions that depend on time only and α k the coefficients that only depend on the uncertain inputs ξ that generated the trajectory. Similar expressions are obtained with (V (t m , ξ i )) m∈[0,n] . In our case, this approach allows to represent the full trajectory of dozens of thousand points with only a couple of basis functions with an average relative error lower than 1 × 10 -4 .

In order to generate a large number of trajectories at low computational cost, a surrogate model is built to create a mapping between uncertain input parameters ξ and the coefficients α k of the reduced basis. Then using Eq. 27, a large number of trajectories can be reconstructed from the reduced basis. The overall procedure yields a complete trajectory (X(t m , ξ)) m∈[0,n] , (V (t m , ξ)) m∈[0,n] from any ξ with a very satisfying level of error. The overall absolute error estimated on a validation set of 500 additional trajectories are reported in table 6 and are much smaller than the observed uncertainties. Once the trajectory (X(t m , ξ)) m∈[0,n] , (V (t m , ξ)) m∈[0,n] is computed, t init (ξ) and t end (ξ) are derived from Eqs. 14 and 15 and t f rag (ξ) is derived from eq. 13 or equivalently eq. 26.

At breakup, the quantities of interest are the altitude, longitude, latitude and velocity (their distributions are represented in Figs. 20a 20b, 20d and20c).

According to our model, the breakup occurs at random between t init and t end according to a uniform distribution. The range of possible breakup goes from 50 Certain studies assume a constant breakup at 80-78 km of altitude [START_REF] Morand | CNES activities on polynomial chaos expansion for uncertainty propagation[END_REF]. While this study confirms that, in the case of a GTO re-entry of a US, a breakup can take place at 80km, it is more likely to happen at lower altitudes. Note that the breakup velocity is extremely high ranging from 8250 m/s to 9800 m/s, in comparison with the average reentry velocity at 120 km of 9800 m/s. High-velocity breakups usually occur at high altitude for reentries with relatively low slope while the low velocity scenarii lead to breakups at lower altitudes.

Figure 21 summarizes the results of the sensitivity analysis. The main contributors to the breakup conditions can be divided into two groups: the deorbiting conditions (A boost ,α boost ,β boost ) and the breakup model parameters (T f rag , p f rag and u f rag ). Recall that u f rag represents the uncertainty due to the random nature of the breakup model (see eq. 26). For the altitude and velocity, the main contributors are u f rag and to a lesser extent T f rag which triggers t init .

Hence, the stochastic nature of the breakup model is the most significant contributor to the velocity and altitude variance compared to the other uncertain inputs. Altogether, u f rag and t init explain almost 80 % of the variance in the altitude and more than 90 % of the velocity variance. For those two variables, it is clear that the breakup model is responsible for most of the uncertainties and that higher fidelity breakup models could reduce the uncertainties significantly. Regarding, the longitude and the latitude of the breakup, the analysis of variance shows that the breakup model parameters account for 30% of variance The computed distributions can be highly valuable when assessing solvers for human risk such as ELECTRA ( [START_REF] Lazare | Electra launch and re-entry safety analysis tool[END_REF]) that requires a breakup altitude to compute the risk for the human population. Finally, the proposed breakup model is flexible in the sense that it could easily exploit a more precise characterisation of the breakup distribution if expert judgements or a priori analyses are available. This additional information can be implemented by either changing the shape of the distribution or introducing multiple breakup times. Future works could also investigate more informative distributions to improve the breakup predictions by using mechanical breakup models and a better characterisation of the material composing the reentering space object.

Figure 1 :

 1 Figure 1: Graphical representation of the breakup predictor with its composing solvers

  : a deorbiting solver, an aerodynamic solver , a trajectory/thermal solver and a probabilistic break up solver. We consider two classes of input variables. The first class is constituted by the ensemble of parameters denoted as global, i.e. which are not the outcomes of an upstream solver. The global inputs are: • Deorbiting solver -The initial orbit elements: apogee altitude (apo ), perigee altitude (per ), inclination (inc), ascending node longitude (anl), perigee argument (pa), true anomaly (tan). -The deorbiting boost characteristics: boost time (t boost ), boost amplitude (A boost ), boost inclination (α boost ), boost bearing (β boost ) • Aerodynamic solver -The object geometry represented by a mesh -The object orientation at 120km represented by the euler angles (denoted θ and φ) • Trajectory/thermal solver -The object material characteristics : material density (ρ M ), emissivity (ε), fusion temperature (T f us ) and fusion enthalpy (H f us ) -The atmosphere model parameters: the solar flux (solarf), the magnetic activity index (ag), the day of reentry (atmoday), the hour of reentry (atmohour) • Probabilistic breakup solver -The breakup model parameters: the fragmentation temperature T f rag and the ablation percentage (p abl )

Figure 2 :

 2 Figure 2: Upper Stage mesh

Figure 3 :

 3 Figure 3: L2 error evolution for a SoGP and single GP (state of the art) surrogate model

Fig. 5

 5 Fig.5represents the first order Sobol indices for altitude, longitude, latitude and velocity. For all four variables, the same inputs account for most of the variance while the majority of the inputs have no significant influence on the output variance. The predominant variables are T f rag which is the temperature at t init and the deorbiting manoeuvre uncertainties such as the boost orientation and amplitude. The variance in altitude and velocity are mostly due to T f rag (respectively 89 % and 95%). The rest is due to the boost amplitude and orientation uncertainties. For latitude and longitude, T f rag explains only 10 % of the output variance while most of the variance is due to the uncertainties associated with the deorbiting boost. Note that the uncertainty in the apogee position has a significant influence in the longitude variance only. Although relatively small, the reentry time affects the altitude at t init . The atmosphere density is dependant on the reentry time in particular at high altitude, and in turn, it changes the convective heat flux.

Figure 4 :

 4 Figure 4: Distributions at t init

Figure 5 :

 5 Figure 5: Sensitivity indices at t init

Figure 6 :

 6 Figure 6: Flight conditions distributions at t end

Figure 7 :

 7 Figure 7: Sensitivity indices at t end

Figure 8 :Figure 9 :

 89 Figure 8: Longitude distribution at 100km, 90km, 85km, 80km, 70km

Figure 10 :Figure 11 :

 1011 Figure 10: Latitude distribution at 100km, 90km, 85km, 80km, 70km

Figure 12 :Figure 13 :

 1213 Figure 12: Velocity distribution at 100km, 90km, 85km, 80km, 70km

Figure 14 :

 14 Figure 14: Convective heatflux sensitivity indices at 100km, 90km, 85km, 80km, 70km

Fig. 15

 15 Fig. 15 reports the pdf of the convective heat flux distributions on a log scale.Although the distribution supports seem to shrink as the altitude decreases due to the log scale, the standard deviation doubles from a 364 kW at 100km to 670 kW at 70 km.

Figure 15 :

 15 Figure 15: Convective heatflux distributions at 100km, 90km, 85km, 80km, 70km

  vre uncertainties become more and more important as the altitude decreases to reach a total of 45 % at 70 km. The atmosphere uncertainties have limited influence in the radiated heat flux around 90 km and they become completely negligible at 70km. At 80 km the first order Sobol indices represent 80 % of the

Figure 16 :

 16 Figure 16: Radiated heatflux sensitivity index at 100km, 90km, 85km, 80km, 70km

Figure 17 :Figure 18 :

 1718 Figure 17: Radiated heatflux distributions at 100km, 90km, 85km, 80km, 70km

Figure 20 :

 20 Figure 20: Flight conditions distributions at breakup

Figure 21 :

 21 Figure 21: Sensitivity indices at breakup

Table 1 :

 1 Case description

	Variable	Value
	Internal mass	4000 kg
	Shell mass	3000 kg
	Material composing the shell	Aluminium
	Dimensions	diameter 4.5, length (with nozzle) 7 m

Table 2 :

 2 List of uncertainties

Table 4 :

 4 Main statistics of the flight conditions at t init

		Altitude [m] Longitude [deg] Latitude [deg] Velocity [m/s]
	mean	71412	-171.97	8.627	9722.8
	std	4327	1.3	0.129	81.9
	0.5 %	63352	-174.84	8.26	9524.
	2.5 %	64367	-174.29	8.35	9553.
	50 %	70962	-172.05	8.64	9739
	97.5 %	79804	-169.22	8.86	9828
	99.5 %	81223	-168.37	8.91	9831

skewed toward high speeds that correspond to cases where t init is reached at high altitudes. Longitude and latitude distributions are relatively symmetrical although slightly skewed toward high latitudes and low longitudes.

Table 5 :

 5 Main statistics of the flight conditions at t end

		Altitude [m] Longitude [deg] Latitude [deg] Velocity [m/s]
	mean	56714.	-169.35	8.37	8644.
	std	2778.	1.45	0.16	173
	0.5%	51782.	-172.2	7.9	8247
	2.5%	52466.	-171.8	8.0	8323
	50%	56289.	-169.5	8.4	8637
	97.5%	63260.	-166.2	8.6	9001
	99.5%	65750.	-165.1	8.7	9097
	that the material characteristics are known well enough compared to the other
	sources of uncertainties.

Table 6 :

 6 L2 absolute error averaged over the whole trajectory

	Variable	absolute error
	Altitude[m]	20
	longitude [deg]	0.02
	latitude [deg]	0.002
	Velocity [m/s]	0.9
	km to more than 80 km. The most probable breakup altitude is around 62km.

Acknowledgement

We thankfully acknowledge the financial support of ArianeGroup and Region Nouvelle Aquitaine. We also thank Celia Finzi and Gregory Pinaud for their fruitful advice on breakup modelling and Francois Lemaire for providing reentry uncertainties.

Flight conditions at breakup

An additional objective of this work includes deriving the flight conditions when breakup occurs. Those quantities are of primary interest for survivability tools in order to compute the ground impact area. In this work, the breakup model bypasses complex and expensive computational models and assumes it occurs at a random time t f rag between t init and t end according to a uniform distribution. There are therefore two sources of randomness in this model. The first one comes from the input uncertainties ξ that characterize the distribution of t init (ξ) and t end (ξ) presented in Section 4.4.1 and 4.4.2. The second level of randomness comes into play when t f rag is sampled from U(t init (ξ), t end (ξ)) (see eq. 13). In other words

where u f rag is a uniform random variable between 0 and 1 that represents the second level of uncertainty. In order to derive the breakup position and velocity distribution, one has to sample a large number of ξ i and for each ξ i to compute the trajectory associated. Then, for each trajectory, a time breakup time is sampled between t init (ξ i ) and t end (ξ i ). Each trajectory is a sequence of position and velocity computed at different times, from t 0 to t n . Each trajectory (indexed