# ZERO-PRESSURE BALLOONS TRAJECTORY PREDICTION: DUSTER FLIGHT SIMULATIONS.

Vito Liuzzi¹, ¹, ², Vincenzo Della Corte², Alessandra Rotundi¹, ², StavroL. Ivanovski³, Zelia Dionnet¹, ², Rosario Brunetto⁴, Laura Inno¹.

1) University of Naples “Parthenope”, Department of Sciences and Technologies, CDN IC4, 80143, Naples, Italy

2) INAF – Institute for Space Astrophysics and Planetology, Area Ricerca Tor Vergata, Via Fosso del Cavaliere 100, 00133, Rome, Italy

3) Osservatorio Astronomico di Trieste, INAF, Via Tiepolo 11, 34143, Trieste, Italy

4) Institut d'Astrophysique Spatiale, CNRS, UMR-8617, Université Paris-Sud, Batiment 121, F-991405 Orsay Cedex, France

*E – mail: vito.liuzzi@uniparthenope.it

## Abstract

DUSTER (Dust from the Upper Stratosphere Tracking Experiment and Retrieval) is a balloon-born instrument designed to collect dust particles with sizes less than 30 microns in the upper stratosphere (30 – 40 km), to be analysed in laboratory. In support to the DUSTER flight campaigns, it is crucial to predict the balloon trajectories according to atmospheric circulation and to estimate the landing. For a feasibility study of an Antarctic flight campaign we predict the DUSTER trajectory with a statistical approach, considering safety issues and the pivotal instrument recovery, to retrieve the collected samples. To this aim, we used the NOAA (National Oceanic and Atmospheric Administration) meteorological data and the HYSPLIT software (HYbrid Single Particle Lagrangian Integrated Trajectory). The Antarctic Polar Vortex, every year in the period December – January, generates a counter clockwise circulation (averaged wind speeds around 36 m/s from 29 to 35 km of altitudes). We evaluated the best launch opportunity windows for DUSTER to take advantage of the Polar Vortex starting from historical data. The NOAA meteorological data (Global reanalysis of atmospheric data – GFS Analysis and Forecasts products of the National Centers for Environmental Predictions - NCEP) from 1994 to 2015 allowed us to select the period after the 10th of January as the launch window maximizing the DUSTER time of flight. We used HYSPLIT program to project the balloon simulated trajectories on ground, providing critical information for the recovery operations. Successively, with the aim of becoming autonomous within the DUSTER team, i.e. taking care of the complete procedure to realize a DUSTER flight campaign, we developed a dedicated software, the DUSTER Flight Simulation Predictor (DFSP) that considers the ascent, floating and descent phases and the meteorological conditions. DFSP considers the balloon dynamics, the heat flow transfer between the balloon and the atmosphere, and the lifting variable gas mass to simulate the complete flight. We are now able to obtain balloon trajectory simulations with a flight time up to 240 hours, starting and ending in any Earth geographical position. In the following paper, a comparison between real DUSTER flight data in 2009 and simulated DFSP data will be shown.

**Keyword:** Zero-pressure balloons; Numerical Simulations; DUSTER Flight Simulation Predictor.

## 1 Introduction

Stratospheric balloons are useful platforms dedicated to multiple aims:1) to study atmospheric phenomena, dynamics and composition;2) to perform astronomical observations;3) to test new technologies in extreme environments. Depending on the payload characteristics and aims it is critical to be able to predict the balloon trajectory, especially if its recovery is pivotal. The precise prediction of a stratospheric zero-pressure
balloon trajectory is a complex task as shown by several studies(e.g. Franco and Cathey 2004; Farley R.E., 2005; Dai et al. 2012) considering the physics and the balloon characteristics. Investigations of the uncertainties tied to the environment forecasts have been performed (Jewtoukoff et al., 2016): the uncertainty originating from the wind forecasts constitutes a major difficulty for balloon trajectory prediction, regardless the challenging physics of the flight. Jewtoukoff et al., 2016, found that errors mainly originate at the floating phase, even with limited altitude variations, because of the poor representation of mesoscale motions in numerical weather prediction models, e.g. inertia-gravity waves, with amplitudes of several meters per second. Sóbester et al., (2014) examine as sources of prediction errors the uncertainties in the wind profile, but point out as the main error source the incomplete understanding of the atmospheric drag on the balloon. They propose a new stochastic empirically derived sphere drag model based on radiosonde flights data providing a simple analytical approximation to be used for the drag computation. Forecasting balloon positions, altitudes, and landing site location as a function of the launch time is of paramount importance for DUSTER (Della Corte et al., 2012). DUSTER collects particles with size less than 30 microns in the upper stratosphere (30 – 40 km) for subsequent laboratory analysis (Della Corte et al., 2013; 2014; Rietmeijer et al., 2016). Since it is mandatory to recover DUSTER to retrieve the collected dust samples to be analysed in the laboratory, we put specific effort to simulate the balloon trajectory and the landing site. In this work, we present the DUSTER Flight Simulation Predictor (DFSP), which is specifically designed for DUSTER flights. It predicts the balloon trajectories from the take-off to the landing site taking into account also the meteorological conditions provided by the forecast services, in particular we used the Global Forecast System – GFS of the NOAA (National Oceanic and Atmospheric Administration). In the first part of the paper, we describe how we used HYSPLIT to determine the possible trajectories followed by zero-pressure balloons over the polar region following the Polar Vortex circulation. In section 3 we present the dynamical and thermal equations written for the DFSP code to introduce in section 4 the entire development of the code and a simple test case where we assume standard atmospheric profile. DFSP predictions in the test case agree with what is expected within the framework of standard assumptions. Thus, we proceed to include a more complex treatment of the atmosphere by using real meteorological data. We did that by using the meteorological data and flight records of DUSTER mission of 2009. The results are discussed in Section 5.

## 2 Modelling the balloon floating phase

By means of the NOAA Air Resources Laboratory HYSPLIT (Hybrid Single-Particle Lagrangian Integrated Trajectory) coupled with historical NOAA meteorological data (Reanalysis data on global scale, with a scale of 1.0 deg by 1.0 deg and an output time step of 6 hours) from 1994 to 2015 (Stein et al., 2015), obtained by the NOAA website<fnref n="1" />, we simulated DUSTER trajectories (Liuzzi et al., 2017) for a possible balloon launch campaign from the Antarctic Italian Station “Mario Zucchelli” (74° 41′ S – 164° 07′ E). HYSPLIT requires only U, V (the horizontal wind components) and W (vertical wind component) (Draxler et al., 2014) to produce the simulation outputs. We retrieved these information using the NOAA meteorological data. The main aim was to determine the DUSTER launch window maximizing the balloon residence time in the Polar Vortex Circulation to optimize stratospheric dust particles collection. This period is related to the wind speeds and the wind directions. During the intensification of the vortex, the circulation becomes circular over the polar region and the wind speeds reach a mean value of 36 m/s. These conditions let the balloon fly more time in the vortex. The simulation provided the balloon latitude, longitude and altitudes for two specific periods of the year: from 1 – 20 December to 10 - 20 January, when the Polar Vortex is at its maximum stability. The simulations starting time is 00:00 (UTC) with a flight time of 240 hours (10 days) considering three different

<footnote n="1">¹https://www ready.noaa.gov/ready2-bin/extract/extracta.pl</footnote>
floating levels above the mean sea level (29 km, 30 km and 31 km) and five different locations around the launching site (a circle with a radius of 2 km is considered). The mean sea level has been chosen as reference level instead of the ground level (default of HYSPLIT). HYSPLIT simulates the balloon possible trajectories around Antarctica. The simulations are repeated every 5, 48 and 150 hours and then the averaged values are displayed. To analyze data dispersion, we obtained the mean and the standard deviation of the balloon coordinates from the simulation outputs. The time intervals have been selected to understand if major trajectory variations happen on different intervals with respect to fixed data. Considering all trajectories, the first point position during the floating phase is coincident with the position of the launching site at the floating altitude. Then, the dispersion is calculated only considering the last point of the simulations. The trajectories displayed in **Figure 1**, are computed by assuming a launch date corresponding to the 1st of December, and they have a large dispersion, with landing sites that are mostly located on the ocean, where retrieving the instrument would be impossible. On the contrary, the trajectories shown in **Figure 2**, which were computed for a launch on the 10th of January are more uniform and lead to landing site mostly located on the continent. Therefore, we selected this period as the best period to launch DUSTER from Antarctica: the balloon permanence in the Polar Vortex is maximized and the stratospheric particle collection optimized.

<figure><img src="image_1.png" /><figcaption>**Figure 1** –Simulated trajectories obtained for a launch on the 1st of December by HYSPLIT (Stein et al., 2015) from 1994 – 2015; Image obtained by Google Earth – Image Landsat – Image U.S. Geological Survey – Data SIO, NOAA, U.S. Navy, NGA, GEBCO – Image PGC/NASA.</figcaption></figure>
<figure><img src="image_2.png" /><figcaption>**Figure 2** – Simulated trajectories obtained by Google Earth – Image Landsat – Image U.S. Geological Survey – Data SIO, NOAA, U.S. Navy, NGA, GEBCO – Image PGC/NASA.</figcaption></figure>

## 3Modelling the entire balloon flight

The main problem with simulating balloon trajectories using HYSPLIT model are related to the ascent and descent phases of the flights. The HYSPLIT software could be used to define the trajectories of the passive tracer in the atmosphere, but not for the stratospheric balloon flight simulation (the balloon filled with gas cannot be considered a passive tracer). By aiming at becoming independent in the numerous campaign's preparation, we decided to develop a DUSTER Flight Simulation Predictor (DFSP). In fact, DUSTER is a dust collection facility designed to fly on a regular basis, in the wake of the NASA (National Aeronautics and Space Administration) stratospheric aircraft dust collections (e.g. Rotundi et al. 2007), and up to date, since 2006, faced already 8 flight campaigns. In particular, in 2006 the qualification flight, from Kiruna, Sweden with the Centre National d'Études Spatiales - CNES, (Palumbo et al. 2007); in 2008 and 2009, two collection flights from Longyearbyen, Svalbard with the Italian Space Agency - ASI, (Della Corte et al. 2012, 2013, Rietmeijer et al., 2016; 2017, De Angelis et al., 2011); in 2011a collection flight, from Kiruna, Sweden (Della Corte et al., 2014; Rietmeijer et al., 2016); in 2014a campaign from Longyearbyen, Svalbard, but no flight was performed due to meteorology conditions; in 2017a collection campaign from M. Zucchelli base, Antarctica, with ASI (Della Corte et al. 2019; Liuzzi et al., 2017); in 2018a collection flight from Longyearbyen, Svalbard, with ASI (dust particles collection under analysis); and in September 2019 a flight campaign was held from Kiruna, Sweden, in the frame of the HEMERA project (the European Union's Horizon 2020 research and Innovation programme under grant agreement No 730970). DFSP is a model that we developed to obtain information on the complete flight profile, which considers the forces and the heat fluxes acting on the balloon. However, in this paper we focus on simulating flight conditions similar to the ones occurred during DUSTER's flight in 2009, in order to demonstrate the reliability of DFSP. The balloon is considered as a free body, a dynamical and thermal system immersed in the real environment conditions, i.e. wind, turbulence, solar radiation. The lifting/buoyancy force is considered by taking into account the air density, the density of the lifting gas, the balloon mass and volume, and the free lift. The balloon thermal behaviour is also taken into account, since the solar radiation induces temperature variations of the air and of the lifting gas, which are crucial to determine the balloon vertical speed.
In fact, temperature variations control internal and external convective fluxes during the day and night periods, thus the balloon flight profile. In addition, for the thermal modelling we consider the albedo that is the energy scattered by the ground. Note that we only considered a constant diurnal cycle in the DFSP software, because every DUSTER mission was hold during the summer period, both in the northern hemisphere and in the southern hemisphere.

To study the physics of a stratospheric balloon and the difference forces acting on it, we need to know the different balloon parameters (Table 1).

<table><caption>**Table 1** – Balloon natural shape, i.e. the balloon shape when all the acting forces are balanced, and geometrical parameters, critical values for dynamic and thermal modelling (Farley R.E., 2005).</caption><tr><td>Balloon diameter</td><td>Diameter = 1.383 · Volume<sup>1/3</sup></td></tr><tr><td>Length of gore</td><td>Lengthgore = 1.914 · Volume<sup>1/3</sup></td></tr><tr><td>Height of the balloon</td><td>Height = 0.748 · Diameter</td></tr><tr><td>Surface of the balloon</td><td>A<sub>surf</sub> = 4.94 · Volume<sup>2/3</sup></td></tr><tr><td>Top Projected Area</td><td>A<sub>top</sub> = π/4 · Diameter<sup>2</sup></td></tr><tr><td>Projected Area</td><td>A<sub>p</sub> = A<sub>top</sub> · [0.9125 + 0.0875 · cos(π)]</td></tr></table>

Where A is the Area of the balloon [m²], A<sub>p</sub> is the Projected area of the balloon [m²] and Volume is the balloon maximum volume without having a bursting event [m³].

## 3.1 Modelling balloon dynamics

To model the balloon dynamics, we make the following assumptions:

- the balloon has three degrees of freedom;

- the atmospheric air and the lifting gas follow the ideal gas law;

- the temperature of the lifting gas inside the balloon is uniform;

- the temperature of the balloon is uniform along its entire surface.

To model the dynamical behaviour of a stratospheric balloon is necessary to consider all the acting forces: the net buoyancy force, the weight and the aerodynamic drag.

The *Buoyancy force* can be written as

$$F_{Arch} = g \rho_{air} Volume \quad (\text{Eq. 01})$$

where ρ<sub>air</sub> is the air density [kg/m³], g is the gravitational acceleration [m/s²] and the Volume refers to the volume of the lifting gas:

$$Volume = M_{gas} R_{gas} \frac{T_{gas}}{p_{gas}} \quad (\text{Eq. 02})$$
Where $M_{gas}$ is the mass of balloon gas [kg], $R_{gas}$ is the Specific gas constant of the lifting gas (for Helium: 2078.5 J/kg K), $p_{gas}$ is the lifting as (Helium) pressure [N/m²] and $T_{gas}$ is the lifting gas (Helium) temperature [K].

The gravitational force is

$$W = M_{gross}g \quad (\text{Eq. 03})$$

where the gross mass, $M_{gross}$, is

$$M_{gross} = M_{film} + M_{load} + M_{gas} + M_{ballast} \quad (\text{Eq. 04})$$

where $M_{film}$ is the Mass of balloon film [kg], $M_{load}$ is the Mass of payload [kg] and $M_{ballast}$ is the Total ballast [kg].

The Drag Force, Drag, is given by:

$$\mathbf{Drag} = \frac{1}{2} \rho_{air} \cdot \mathbf{V}_{rel}^2 \cdot C_d \cdot A_{top} \quad (\text{Eq. 05})$$

where $C_d$ is the aerodynamic coefficient, $A_{top}$ is the balloon top projected area, $V_{rel}$, is the relative velocity that can be expressed as a function of the x, y and z directions:

$$V_{relx} = V_x - V_{windex} \quad (\text{Eq. 06}),$$

$$V_{rely} = V_y - V_{windy} \quad (\text{Eq. 07})$$

$$V_{relz} = V_z - V_{windz} \quad (\text{Eq. 08})$$

$$V_{rel} = \sqrt{V_{relx}^2 + V_{rely}^2 + V_{relz}^2} \quad (\text{Eq. 09})$$

Note that Drag is particularly important during the ascent of the balloon, when the relative velocity is higher and the shape of the balloon – therefore its inertia – can significantly change.

The differential equations governing the free balloon motion are:

$$a_x = \frac{Drag_x}{M_{total}} \quad (\text{Eq. 10})$$

$$a_y = \frac{Drag_y}{M_{total}} \quad (\text{Eq. 11})$$

$$a_z = \frac{F_{Arch} - W - Drag_z}{M_{total}} \quad (\text{Eq. 12})$$

where the total mass, $M_{total}$, is

$$M_{total} = M_{gross} + M_{added} \quad (\text{Eq. 13})$$

and the added mass, $M_{added}$, defined as

$$M_{added} = (\rho_{air} \cdot Volume) \cdot C_{added} \quad (\text{Eq. 14})$$

takes into account the dragged mass of air during the ascent phase. In fact, it considers the mass of the fluid displaced by a moving solid body. The parameter $C_{added}$ we used is 0.5 in accordance with Palumbo et al. (2007).

During the landing phase the gondola uses a parachute system installed since the initial departure phase.
The simplest way to evaluate the velocity of the parachute is considering the mass of payload mp, including the parachute system mass, the descent velocity of the system can be calculated using the relative velocity between the wind velocity and the parachute velocity along the vertical direction. The parachute velocity is calculated solving the following partial differential equation

$$ \frac{dv_{par}}{dt} = g - \frac{1}{2} \frac{C_{Dpar} A_{par}}{m_p} \rho_a v_{par}^2 \quad (\text{Eq. 15}) $$

Where $v_{par}$ is the parachute velocity, $A_{par}$ is the parachute section and $C_{Dpar}$ is the drag coefficient of the parachute used in the mission.

For the horizontal motion, the evaluation of the parachute velocities could be described using the same scheme reported in the Eq. 06-09.

## 3.2 Thermal Model: heat fluxes on the balloon

The thermodynamics processes inside the balloon affect the balloon shape, thus its dynamic and trajectory (Farley R.E., 2005).

In our model we consider these heat sources:

- Direct Sun Radiation, or Direct Solar Flux, $q_{Sun}$ is equal (at the balloon altitude Z), to the product of its value at the top of the atmosphere $I_{Sun}$ and the atmosphere transmittance, $\tau_{atm}$:

$$ q_{Sun} = I_{Sun} \cdot \tau_{atm} \left[ \frac{W}{m^2} \right] \quad (\text{Eq. 16}) $$

- Thermal Infrared Radiation ($q_{IRground}$), indicates the infrared flux coming from the ground. It depends on the ground emissivity, $\epsilon_{ground}$, and on the temperature of the ground, $T_{ground}$:

$$ q_{IRground} = \epsilon_{ground} \cdot \sigma \cdot T_{ground}^4 \left[ \frac{W}{m^2} \right] \quad (\text{Eq. 17}) $$

where $\sigma$ is the Stephan-Boltzmann constant, $5.67e-8$ [W/m²/K⁴]

The ground emissivity value used in the simulation is

$$ \epsilon_{ground} = 0.95 \quad (\text{Eq. 18}) $$

The infrared flux from ground at the balloon altitude Z, considering the atmosphere transmittance $\tau_{atmIR}$, is

$$ q_{IREarth} = q_{IRground} \cdot \tau_{atmIR} \left[ \frac{W}{m^2} \right] \quad (\text{Eq. 19}) $$

where $q_{IREarth}$ is the diffuse planetary longwave flux (Earth) [W/m²].

- Thermal Flux Induced by the Albedo ($q_{Albedo}$):

$$ q_{Albedo} = Albedo \cdot I_{Sun} \cdot \sin(ELV) \left[ \frac{W}{m^2} \right] \quad (\text{Eq. 20}) $$
where ELV is the Sun's elevation. The Albedo, i.e. coefficient accounting for the fraction of solar radiation reflected by the Earth surface and the atmosphere, depends on the overflown region properties and on the presence of possible cloud. For example, if the balloon flies over a zone with a high amount of ice the albedo value is 1.0; in a generic simulation it is possible to use 0.5 as mean value.

- *Convective Heat Transfer between the atmosphere and the lift gas:*

The convective heat transfer coefficients depend on the atmospheric and the lifting gas conditions.

These conditions are based on the values of the following parameters (Table 2).

<table><caption>**Table 2** – Physical parameters of the air and the gas that influence the interactions between these ones (Farley R.E., 2005).</caption><tr><td>Air dynamic viscosity</td><td>μ<sub>air</sub> = <sup>1.458 · 10<sup>-6</sup> · T<sub>air</sub><sup>1.5</sup></sup>&frasl;<sub>T<sub>air</sub> + 110.4</sub></td></tr><tr><td>Helium dynamic viscosity</td><td>μ<sub>gas</sub> = 1.895 · 10<sup>-5</sup> · (<sup>T<sub>gas</sub></sup>&frasl;<sub>273.15</sub>)<sup>0.647</sup></td></tr><tr><td>Air conductivity</td><td>k<sub>air</sub> = 0.0241 · (<sup>T<sub>air</sub></sup>&frasl;<sub>273.15</sub>)<sup>0.9</sup></td></tr><tr><td>Helium conductivity</td><td>k<sub>gas</sub> = 0.144 · (<sup>T<sub>gas</sub></sup>&frasl;<sub>273.15</sub>)<sup>0.7</sup></td></tr><tr><td>Air Prandtl Number</td><td>Pr<sub>air</sub> = 0.804 – 3.25 · 10<sup>-4</sup> · T<sub>air</sub></td></tr><tr><td>Helium Prandtl Number</td><td>Pr<sub>gas</sub> = 0.729 – 1.6 · 10<sup>-4</sup> · T<sub>gas</sub></td></tr></table>

Following the radiation law, the reflectivity, $r_\lambda$, of a body can be expressed by the absorptivity ($\alpha$) and the transmissivity ($\tau$) coefficients

$$r_\lambda = 1 - \alpha - \tau \quad (\text{Eq. 21}).$$

The temperature variations of the balloon film in the model can be expressed by this formula (Palumbo et al., 2007)

$$\frac{dT_{film}}{dt} = \frac{Q_{Sun} + Q_{Albedo} + Q_{IREarth} + Q_{IRfilm} - Q_{IROut}}{c_f \cdot M_{film}} \left[ \frac{K}{s} \right] \quad (\text{Eq. 22})$$

where $c_f$ is the Specific heat of the film material [J/Kg/K] and the heat powers acting on the balloon are reported in (Table 3).
<table><caption>**Table 3** – List of the heat powers acting on the balloon.</caption><tr><td>1.</td><td>Total Absorbed Direct Sun radiation</td><td>Q<sub>Sun</sub> = α · A<sub>p</sub> · q<sub>Sun</sub> · (1 + <sup>τ</sup>&frasl;<sub>1 &minus; r</sub>) [W]</td></tr><tr><td>2.</td><td>Total Absorbed Earth Infrared Radiation</td><td>Q<sub>IREarth</sub> = α<sub>IR</sub> · A<sub>surf</sub> · q<sub>IREarth</sub> · ViewFactor · (<sup>τ<sub>IR</sub></sup>&frasl;<sub>1 &minus; r<sub>IR</sub></sub>) [W]</td></tr><tr><td>3.</td><td>Total Absorbed Albedo Heat</td><td>Q<sub>Albedo</sub> = α · A<sub>surf</sub> · q<sub>Albedo</sub> · ViewFactor · (1 + <sup>τ</sup>&frasl;<sub>1 &minus; r</sub>) [W]</td></tr><tr><td rowspan="2">4.</td><td>Absorbed Infrared(IR) Film Radiation</td><td>Q<sub>IRAds</sub> = σ · ε · α<sub>IR</sub> · A<sub>surf</sub> · T<sup>4</sup><sub>film</sub> · (<sup>1</sup>&frasl;<sub>1 &minus; r<sub>IR</sub></sub>) [W]</td></tr><tr><td>Emitted Infrared(IR) Film Radiation</td><td>Q<sub>IREm</sub> = σ · ε · A<sub>surf</sub> · T<sup>4</sup><sub>film</sub> · (1 + <sup>τ<sub>IR</sub></sup>&frasl;<sub>1 &minus; r<sub>IR</sub></sub>) [W]</td></tr></table>

## 4 Simulating complete DUSTER flights

By means of the DUSTER Flight Simulation Predictor (DFSP) we are able to predict its balloon flight paths from any launching site using a zero-pressure. DFSP calculating chart flow is reported in **Figure 3**.

<figure><img src="image_3.png" /><figcaption>**Figure 3**– The DUSTER Flight Simulation Predictor represented by a block diagram.</figcaption></figure>

Provided the input of the launch site coordinates, i.e. latitude, longitude and altitude, DFSP starts to retrieve the atmospheric parameters from the meteorological data, available on the National Oceanic and Atmospheric Administration website<fnref n="2" />, and we can use either historical or forecast data. These atmospheric data are formatted as a matrix where the information about temperature, pressure and the wind components, are reported on a spatial grid of one by one degree in latitude and longitude.

<footnote n="2">²https://www.ncdc.noaa.gov/data-access/model-data/model-datasets/global-forcast-system-gfs</footnote>
For the vertical grid step, the terrestrial atmosphere is divided in 64 layers starting from the mean sea level. We implemented an interpolation method to evaluate the atmospheric data for each simulation point and to overcome the limitation of the original grid spacing.

To interpolate the original downloaded atmospheric data we used a linear interpolation between the nearest grid points. DFSP, for each balloon position, evaluates the different forces acting on the system using as inputs: the atmospheric data, the balloon and lifting gas temperatures, the thermal radiation balance (the radiation emitted and absorbed by the balloon). With the forces evaluated for the time tᵢ DFSP is able to compute the position at the time tᵢ+₁.

This procedure is applied iteratively obtaining the simulated position of the balloon entire flight. We performed a first check of the DFSP computational module (Figure 4) using, instead of the atmospheric data, the following standard equations of the air temperature, T<sub>air</sub>, and the air pressure, p<sub>air</sub>, and the derived pressure expression, p, as function of the altitudes (Kayhan and Hastaoglu, 2014) (Xin-Lin Xia et al., 2010) between 0 meters and 35000 meters

$$T_{\text{air}} = \begin{cases} 288.15 - 0.065 \cdot Z & \text{if } 0 \le Z \le 11000 \text{ m} \\ 216.15 & \text{if } 11000 < Z \le 20000 \text{ m} \quad (\text{Eq. 23}) \\ 216.15 + 0.0012 \cdot (Z - 20000) & \text{if } 20000 < Z \le 35000 \text{ m} \end{cases}$$

$$p_{\text{air}} = \begin{cases} 101325 \cdot \left(1 - \frac{Z}{44330}\right)^{5.26} & \text{if } 0 \le Z \le 11000 \text{ m} \\ 22605 \cdot e^{\frac{11000-Z}{6340}} & \text{if } 11000 < Z \le 20000 \text{ m} \quad (\text{Eq. 24}) \\ 2447 \cdot \left(\frac{141.89+0.003 \cdot Z}{216.65}\right)^{-11.388} & \text{if } 20000 < Z \le 35000 \text{ m} \end{cases}$$

$$p = p_0 \cdot \left(1 - \frac{\text{LapseRate} \cdot Z}{T_0}\right)^{\frac{g \cdot M_B}{R \cdot L}} \quad (\text{Eq. 25})$$

where the LapseRate indicates the mean temperature variation with the altitude and it's equal to 0.0065 [°C/m] for the troposphere, p₀ is the pressure at mean sea level [Pa], R is the average radius of the Earth, i.e. 6378 [km], L is the Temperature lapse rate [K/m] and M<sub>B</sub> is the Molecular weight of the gas (air) [kg/kmol]. But this model considers an ideal case and it doesn't take into account the atmospheric circulation and conditions. When real meteorological data are used, the simple linear differential equations don't give accurate results and more complex models have to be used.
<figure><img src="image_4.png" /><figcaption>Figure 4 – Balloon trajectory for an entire flight simulated by the DUSTER Flight Simulation Predictor using standard values of temperature (T<sub>air</sub>) and pressure (p<sub>air</sub>) profiles. In this case the motion of the balloon is only along the vertical direction and the wind components (longitudinal and meridional components) are zero in this case.</figcaption></figure>

Thus, we implemented two different mathematical methods: 1) the Fourth order Runge-Kutta method for the integration of the differential equations; 2) the Monte Carlo method to take into account the errors on the atmospheric parameters obtained from the atmospheric forecast model outputs.

For the Fourth order Runge-Kutta method all equations can be expressed in the vector form:

$$y' = f(t, y) \quad (\text{Eq. 26})$$

i.e. the time-dependent function that we are using in the model to calculate the variable y. This method works very well with a nonlinear initial value system of ordinary differential equations, having an error proportional to the fourth power of the time step set in DFSP. In our simulations it gets us the possibility to reduce the error values for the simulated positions along the horizontal and the vertical directions as described in Section 5.

Following Qiumin Dai et al. (2012) we describe the initial problem using four coefficients ($k_1, k_2, k_3$ and $k_4$):

$$y_{i+1} = y_i + \frac{\Delta t}{6} * (k_1 + 2k_2 + 2k_3 + k_4) \quad (\text{Eq. 27})$$

where

$$k_1 = f(t_i, y_i)$$

$$k_2 = f(t_i + \frac{\Delta t}{2}, y_i + \frac{\Delta t}{2} * k_1)$$

$$k_3 = f(t_i + \frac{\Delta t}{2}, y_i + \frac{\Delta t}{2} * k_2)$$

$$k_4 = f(t_i + \Delta t, y_i + \Delta t * k_3)$$

This method allows, starting from the balloon flight information about the buoyancy force and the thermal behaviour at certain coordinates, $y_i$, to derive the next step information, $y_{i+1}$. DFSP uses the Monte Carlo
probabilistic method to forecast how the balloon trajectory is influenced by the error on the atmospheric parameters. With this mathematical method we randomly change the flight scenario selecting different values for parameters (e.g. wind speeds, air temperatures) to get the highest probable balloon trajectory. This approach is critical to estimate the landing site to plan a feasible recovery of the DUSTER instrument. In the simulation the variables that are changed randomly are the air temperature, the horizontal component, the meridional component and the vertical component of the wind. For the temperature the uncertainty on the original value is less than 10% and the uncertainty on the wind components is less than 25%. Using these randomly variation, the maximum dispersion in the flight simulation discussed in **Section 5** for the total balloon position is in a range from 2.55% to 8.30%. Considering the last simulation position, the maximum dispersion value is 3.3%. Naturally, all simulation results depend on the data used. More precise input data will result in more precise output.

## 4.1 Prediction of the balloon trajectory during different phases of the flight

Heat, mass and momentum transfer models were developed and coded with different subroutines in MATLAB (see **Section 3.1** and **Section 3.2**). Three different high-level subroutines were developed to study the three balloon main flight phases:

- **Ascent Phase:** The first DFSP routine developed for the balloon ascent foreseen the following conditions as described by Yajima et al., (2009):1) through the troposphere considers a roughly constant ascent rate, 2) in the stratosphere the air temperature increases and the lifting gas temperature decreases (due to adiabatic expansion), thus their difference increases. This phenomenon increases the ascent rate of the balloon due to the lower friction of the less dense air. During the entire ascent phase, the DFSP routine considers the acting vertical forces and the ascent rate increases as the 1/6-power of the balloon volume. Combining the described effects, the balloon, after slowing down in the tropopause, starts to gradually accelerate in the stratosphere.

- **Floating Phase:** In the second DFSP routine, approaching to the floating altitude, the balloon continues to ascend due to inertial motion, then it descends for few meters before setting at the floating altitude (Figure 5). During the ascent and floating phases, the absorbed radiation heat increases with altitude inducing adiabatic expansion. During the floating phase, the balloon absorbs energy, from the Sun and from ground, and diffuses it to space in similar amount, thus the altitude remains constant (Yajima et al., 2009).
<figure><img src="image_5.png" /><figcaption>**Figure 5** – Zooming effect of **Figure 4** to better visualize the balloon altitude variation before reaching the floating phase, where the altitude remains almost constant.</figcaption></figure>

- **Descent Phase:** For the descent phase from the floating altitude the third DFSP routine implies the sudden release of the balloon and the use of the parachute (see **Figure 4**). To define the landing site, the DFSP considers the aerodynamical parachute model and all forces acting on it. The knowledge/forecast of the landing position is crucial to the payload recovery and subsequent analyses of collected samples.

# 5 Results and Discussion

To verify the accuracy of DFSP we compared simulated trajectories with areal DUSTER flight path. In particular, we considered the DUSTER campaign supported by the Italian Space Agency: launch on the 1st July 2009, at 00:57, from coordinates $\phi = 78^\circ 13' N, \lambda = 15^\circ 38' E$, Longyearbyen, Svalbard. DUSTER was launched as a piggyback payload using a zero-pressure balloon filled with helium. The weight of the main payload was 800 kg. The total floating time was around 3 days and the real maximum floating altitude of the balloon during the flight was 39801 meters. The positions (vector position along the x, y, z directions) to reconstruct the simulated trajectory were calculated using the *Fourth order Runge-Kutta*, to get a good stability in the equations solution and with the *Monte Carlo* simulation, to determine the confidence of results. The downloaded atmospheric data (downloaded by the NOAA Air Resources Laboratory archive website; see **Section 4**) are the results of atmospheric forecast models corrected by the actual measurements obtained at the time of the flight. By means of the interpolation method previously described, we retrieve for each point of the flight path simulation the atmospheric parameters needed as input for the thermal and dynamical models used. All atmospheric information is contained in this matrix. All other inputs (solar radiation, balloon shape, lifting gas temperature, emitted/absorbed energy) are evaluated at each integration step. Wind speed components were used to evaluate the balloon dynamics (speeds and positions).

In **Figure 6** and **Figure 7** we report the simulated path compared with the real flight trajectory (Latitude, Longitude and Altitude). As can be seen from **Figure 6**, the floating altitude reached is the same of the value predicted by DFSP and the balloon ascent rate is consistent with the real values experienced by DUSTER.
<figure><img src="image_6.png" /><figcaption>**Figure 6**– Comparison between the simulated ascension and floating phase (blue line) and the real corresponding phase experienced (green circles) by DUSTER during the flight campaign performed on the 1st July 2009 from Longyearbyen (Svalbard, Norway).</figcaption></figure>

In Figure 7 we compare the DFSP simulated latitude and longitude with those effectively flown by DUSTER during the July 2009 flight, showing a good correlation with a max error value on the horizontal position of 28.57 meters. For the vertical position the max error value is 14.99 meters. These results were obtained using the fourth order Runge-Kutta method in the simulation scheme. In fact, using only the simple differencing scheme, the max error value for the horizontal position is 1110.9 meters and the max error value for the vertical position is 499.08 meters, with an improvement on the simulations of the 44.92%. Then, repeating the simulations using *Monte Carlo* method we define the maximum position errors: for the horizontal position the standard deviation value is 1.16 km.

<figure><img src="image_7.png" /><figcaption>**Figure 7**– Comparison between the predicted latitudes and longitudes (blue line) and the real ones (green circles) flown by DUSTER during the flight campaign performed on the 1st July 2009 from Longyearbyen (Svalbard, Norway).</figcaption></figure>
In **Figure 8** we report a 3D view of the simulated (green) and the real (red) DUSTER flight. In addition, we report the longitudes and latitudes projected on ground, simulated (magenta) and the flown ones (blue). The evident overlapping of simulation to real data indicates that DFSP is a reliable model.

<figure><img src="image_8.png" /><figcaption>**Figure 8** – The complete 3D trajectory simulated (red line) and flown by DUSTER (green circles) during its flight on 1st July 2009 from Longyearbyen (Svalbard, Norway). In addition, the projection on ground of the latitudes and longitudes simulated (blue line) and the actually flown ones (magenta circles) are reported.</figcaption></figure>

# 6 Conclusion

HYSPLIT is a software designed by the NOAA Air Resources Laboratory to simulate the trajectories of air parcels. We used this software to perform a preliminary analysis of the behaviour of a possible DUSTER flight over Antarctica. This software is useful for studying the flight floating phase. However, since our goal was to get all information about the entire flight, to plan autonomously DUSTER flight campaigns, we developed the "DUSTER Flight Simulation Predictor - DFSP". The DFSP simulates the DUSTER flights considering the dynamical and the thermal behaviour of zero-pressure balloons up to 40 km of altitude using Helium as lifting gas, using atmospheric parameters as input. Since the recorded atmospheric parameters are not continuous in time and space, for the DFSP simulations we interpolated the meteorological data to cover the entire balloon flight region. The balloon trajectories are simulated combining both dynamical and thermal processes, considering the variations of balloon temperature, shape, lifting gas temperature and pressure, emitted and absorbed energies.

To simulate more accurate trajectories and to retrieve information on the reliability of the solution by means of DFSP, we applied both the *Fourth order Runge-Kutta* method (predicting the balloon position starting from the previous one) and the *Monte Carlo* method(to find the most probable trajectory that the balloon would follow).

As a testbed, DFSP was used to simulate the trajectory flown by DUSTER during the flight campaign that took place in July 2009, from Longyearbyen (Svalbard, Norway). We reproduced DUSTER flight trajectory
during the three different flight phases: ascent, floating and descent. The results showed an optimum
agreement between the real DUSTER flight path and the DFSP simulated trajectory: the maximum difference
between the simulated and the real trajectory in altitude is of 14.99 meters and in the position is of
28.57meters, during the entire flight. The standard deviation obtained by the Monte Carlo simulations is 1.16
km.

Acknowledgements

This work was supported by the Italian Ministero dell'Istruzione, Università e Ricerca (MIUR) through the Programma Nazionale delle Ricerche in Antartide (PNRA) "2013/AZ3.01 DUSTER (Dust in the Upper Stratosphere Tracking Experiment and Retrieval)" and the Progetti di Rilevenza Nazionale "Cosmic Dust" project (grant ID:PRIN201520158W4JZ7). DUSTER project has been funded by ASI (Italian Space Agency), MAE (Ministero degli Affari Esteri) and Regione Campania. The authors gratefully acknowledge the NOAA Air Resources Laboratory (ARL) for the provision of the HYSPLIT transport and dispersion model and READY website (http://www ready.noaa.gov) used in this publication. Thank you to Prof. Ezio Bussoletti for his amazing support for this project and our flight campaigns.
Bibliography

1. Ahrens, C.D., 2000, Essential of Meteorology. An Invitation to the Atmosphere, III Edition.

2. Cardillo, A., Memmo, A., Musso, I., Ibba, R., Spoto, D., 2008, Analysis and prediction of stratospheric balloons trajectories, Vol. 79, 915.

3. Dai, Q., X. Fang, X. Li, and L. Tian, 2012, Performance simulation of high altitude scientific balloons. Adv. Space Res., 49, 1045-1052, doi:10.1016/j.asr.2011.12.026.

4. Das, T. & Mukherjee, R., 2003, Optimal Trajectory Planning for Hot - air Balloons in Linear Wind Fields, Journal of Guidance, Control and Dynamics, Vol. 26, No. 3.

5. De Angelis, S., Della Corte, V., Baratta, G.A., Rietmeijer, F.J.M., Brunetto, R., Palumbo, P., Ciucci, A., and Rotundi, A., 2011, Raman Microspectroscopy Performed on Extraterrestrial Particles, Spectroscopy Letters, 44:7-8, 549-553.

6. Della Corte, V., Palumbo, P., De Angelis, S., Ciucci, A., Brunetto, R., Rotundi, A., Rietmeijer, F.J.M., Zona, E., Bussoletti, E., Colangeli, L., Esposito, F., Mazzotta Epifani, E., Mennella, V., Peterzen, S., Masi, S., Ibba, R., 2011, DUSTER (Dust in the Upper Stratosphere Tracking Experiment and Return): a balloon - borne dust particle collector, Vol. 16, 14.

7. Della Corte, V., Palumbo, P., Rotundi, A. De Angelis, S., Rietmeijer, F. J. M., Bussoletti, E., Ciucci, A., Ferrari, M., Galluzzi, V., Zona, E., 2012, In Situ Collection of Refractory Dust in the Upper Stratosphere: The DUSTER Facility, Space Science Reviews, Volume 169, Issue 1-4, pp. 159-180, DOI: 10.1007/s11214-012-9918-9.

8. Della Corte, V., Rietmeijer, F.J.M, Rotundi, A. and Ferrari, M., 2014, Introducing a new stratospheric dust collecting system with potential use for upper atmospheric microbiology investigations, Astrobiology, Volume 14, Number 8, DOI: 10.1089/ast.2014.1167.

9. Della Corte, V., Rietmeijer, F.J.M, Rotundi, A., Ferrari, M. and Palumbo, P., 2013, Meteoric CaO and carbon smoke particles collected in the upper stratosphere from an unanticipated source, Tellus B, 65, 20174, http://dx.doi.org/10.3402/tellusb.v65i0.20174.

10. Draxler, R.R., 1996, Trajectory Optimization for Balloon Flight Planning, Weather and Forecasting, American Meteorological Society, Vol. 11, 1.

11. Draxler, R.R.& Hess, G.D., 1998, An Overview of The HYSPLIT_4 Modeling System Of Trajectories, Dispersion and Deposition, Australian Meteorological Magazine, 47, 295-308.

12. Draxler, R.R.&Hess, G.D., 2014, Description of the HYSPLIT_4 Modeling System, NOAA Technical Memorandum ERL ARL - 224.

13. Farley, R.E., 2005, Balloon Ascent: 3-D Simulation Tool for the Ascent and Float of High - Altitude Balloons, Aerospace Engineer, Mechanical Systems Center, NASA - Goddard Space Flight Center, AIAA - 2005 - 7412.

14. Franco, H. & Cathey, H.M. Jr., 2003, Thermal performance modelling of NASA's scientific balloons, Advances In Space Research.

15. Gallice, A., Wienhold, F.G., Hoyle, C.R., Immler, F., Peter, T., 2011, Modeling the ascent of sounding balloons: derivation of the vertical air motion, Atmospheric Measurement Techniques, 4, 2235 - 2253.

16. Jewtoukoff, V., Plougonven, R., Hertzog, A., Snyder, C. and Romine, G., 2016, On the Prediction of Stratospheric Balloon Trajectories: Improving Winds with Mesoscale Simulations, Journal of atmospheric and oceanic technology, Vol. 33, pp. 1629-1647, DOI: 10.1175/JTECH-D-15-0110.1.

17. Kayhan, O. and Hastaoglu, M.A., 2014, Modeling of Stratospheric Balloon Using Transport Phenomena and Gas Compress - Release System, Journal Of Thermophysics and Heat Transfer, Vol. 28, No. 3.

18. Kreith, F. & Kreider, J.F., 1974, Numerical Prediction of the Performance of High Altitude Balloons, National Center For Atmospheric Research.
19. Liuzzi, V., Della Corte, V., Ivanovski, S.L., Rotundi, A., Bussoletti, E., 2017, Stratospheric Balloon Trajectory Simulations in the Antarctic Polar Vortex for DUSTER Flights, 23° European Space Agency on European ROCKET & BALLOON Programmes and related Research, Visby (Sweden).

20. Musso, I., Cardillo, A., Cosentino, O., Memmo, A., 2003, A balloon trajectory prediction system, Advances In Space Research.

21. Palumbo, P., Della Corte, V., Colangeli, L., Esposito, F., Rotundi, A., Ciucci, A., Aronica, A., Brucato, J.R., MazzottaEpifani, E., Mennella, V., Battaglia, R., Ferrini, G., Rietmeijer, F.J.M, Flynn, G.J., Renard, J.B., Stephens, J.R., Zona, E., and Inarta, S., 2008, DUSTER – Aerosol collection in the stratosphere, Mem. S.A.It. Vol. 79, 853.

22. Palumbo, R., Mercogliano, P., Corraro, F., De Matteis, P. P., Sabatano, R., 2008, Meteorological conditions forecast and balloon trajectory estimations, Memorie della Società Astronomica Italiana, Vol. 79, 841.

23. Palumbo, R., Russo, M., Filippone, E., Corraro, F., 2007, ACHAB: Analysis Code for High Altitude Balloons, AIAA.

24. Payne, R.E., 1972, Albedo of the Sea Surface, Journal Of The Atmospheric Sciences, Vol. 29.

25. Qiumin Dai, Xiande Fang, Xiaojian Li, Lili Tian, 2012, Performance simulation of high altitude scientific balloons, Advanced In Space Research.

26. Rietmeijer, F.J.M, Della Corte, V., Ferrari, M., Rotundi, A. and Brunetto, R., 2016, Laboratory analyses of meteoric debris in the upper stratosphere from settling bolide dust clouds, Icarus, 266, 217-234, DOI: 10.1016/j.icarus.2015.11.003.

27. Rietmeijer, F.J.M, Ferrari, M., Della Corte, V., Rotundi, A., Palumbo, P., De Angelis, S. and Galluzzi, V., 2017, A fiery birth of alumino-silica analogs of refractory dust in the upper Stratosphere, Advances in Space Research, 60 (2017) 2091-2098, https://doi.org/10.1016/j.asr.2017.07.030.

28. Rolph, G.D., 2017, Real-time Environmental Applications and Display sYstem (READY) Website (http://www ready.noaa.gov), NOAA Air Resources Laboratory, College Park, MD.

29. Rolph, G., Stein, A., Stunder, B., 2017, Real-time Environmental Applications and Display sYstem: READY. Environmental Modelling and Software, 95, 210 – 228, https://doi.org/10.1016/j.envsoft.2017.06.025.

30. Rotundi, A., Ferrini, G., Baratta, G.A., Palumbo M.E., Palomba, E., Colangeli, L., 2007, Combined micro-infrared (IR) and micro-raman measurements on stratospheric interplanetary dust particles, ESA - Workshop on Dust in Planetary Systems, pp. 149-153.

31. Sóbester, A., Czerski, H., Zapponi, N. and Castro, I., 2014, High-Altitude Gas Balloon Trajectory Prediction: A Monte Carlo Model, AIAA Journal, Vol. 52, No. 4, pp. 832-842, https://doi.org/10.2514/1.J052900.

32. Schumann, A.P., 1980, Carrier Balloon Trajectories in the Stratosphere, Journal of Applied Meteorology, American Meteorological Society.

33. Stein, A.F., Draxler, R.R., Rolph, G.D., Stunder, B.J.B, Cohen, M.D., Ngan, F., 2015, NOAA's HYSPLIT atmospheric transport and dispersion modeling system, Bull. American Meteorology Society, 96, 2059-2077, http://dx.doi.org/10.1175/BAMS-D-14-00110.1.

34. Xin - Lin Xia, De - Fu Li, Chuang Sun, Li - Ming Ruan, 2010, Transient thermal behaviour of stratospheric balloons at float conditions, Advances In Space Research.

35. Yajima, N., Izutsu, N., Imamura, T., Abe, T., 2009, Scientific Ballooning, Technology and Applications of Exploration Balloons Floating in the Stratosphere and the Atmosphere of Other Planets.