Europlanet Science Congress 2022
Palacio de Congresos de Granada, Spain
18 – 23 September 2022
Europlanet Science Congress 2022
Palacio de Congresos de Granada, Spain
18 September – 23 September 2022
Computational astrophysics and numerical models of small bodies and planets


The goal of this session is to cover numerical simulations and
relevant laboratory investigations related to the Small Bodies
(comets, KBOs, rings, asteroids, meteorites, dust), their formation
and evolution, and the instruments of their exploration.
This session is specially focused on the interdisciplinary approach
in the development of models (formal descriptions of physical phenomena),
experiments (on ground and in micro-gravity), and mathematical simulations
(computational methods and algorithms of solution) of various astrophysical
phenomena: (i) dusty gas cometary atmospheres; (ii) volcanic activity on icy
satellites (e.g. Enceladus and Io); (iii) planetary body formation
(e.g. via pebbles growth), and planetesimal dynamics.

This session will include an introduction and discussion of new and/or existing
laboratory studies in simulated space-like environments and models, experimental
techniques, computational methods that can address the results of analytical,
experimental and numerical analysis (with respect to computational methods and
algorithms of solution) on the above described studies.
Abstracts on thermophysical evolution models of small bodies interiors as well
as on the modeling of atmosphere and exosphere are welcome.

Conveners: Vladimir Zakharov, Stavro Lambrov Ivanovski, Raphael Marschall | Co-conveners: Luis Diego Pinto, Michelangelo Formisano, Diego Turrini
| Thu, 22 Sep, 10:00–11:30 (CEST)|Room Albéniz+Machuca
| Attendance Thu, 22 Sep, 18:45–20:15 (CEST) | Display Wed, 21 Sep, 14:00–Fri, 23 Sep, 16:00|Poster area Level 2

Session assets

Discussion on Slack

Orals: Thu, 22 Sep | Room Albéniz+Machuca

Chairpersons: Stavro Lambrov Ivanovski, Raphael Marschall
Christian Schuckart, Dorothea Bischoff, Bastian Gundlach, Johanna Bürger, Nicholas Attree, and Jürgen Blum

The activity evolution of comets along their orbits around the Sun has not yet been fully understood. Especially around perihelion, comets show an abnormally high dust and water activity, growing far more rapidly than the increase in energy from the incoming sunlight (e.g. Skorov et al., 2020; Combi et al., 2020). To better understand this complex phenomenon, we developed a thermophysical model that encompasses mass-transfer of volatile species and the resulting build-up of pressure. We based our model on the assumption of a comet nucleus that was formed by gravitational collapse of dust-ice pebbles (Blum et al., 2017).

With these conditions in mind, this model was designed for numerically solving the one-dimensional heat transfer equation (HTE) for porous media. It solves the HTE through an explicit finite-difference scheme to more easily incorporate the temperature-dependent thermal conductivity. The thermal conductivity for a pebble system was implemented in accordance with Gundlach and Blum, 2012. The model calculates sublimating and resublimating particles via the Hertz-Knudsen equation and uses the resulting values to calculate the latent heat transfer. Sub-surface gas pressure levels become available for analysis through this microphysical approach. For ice-dust layers, we calculate the reduced gas flux through dust layers via a reduction function, as used in Gundlach, Fulle and Blum, 2020. Furthermore, dust ejection through pressure build-up for ultra-low tensile strength layers is an included effect, which can be used for calculating the ejected dust mass of a comet.

The current model is designed for applications with dust-ice pebbles, where H2O and CO2 are the main volatile species, but the implementation is adaptable for different scenarios and allows applications for non-porous media and addition of multiple volatile species. 

Figure 1: An example of the resulting temperature profile of our model. It shows the temperature and the volatile species content of the layers as a function of depth for a simulation of comet 67P/Churyumov-Gerasimenko close to perihelion. The two dotted lines mark the current position of the sublimation fronts of the volatile species.

Figure 2: An example of an ejection event of the top three layers during a simulation of comet 67P. The first plot shows the temperature profile and volatile species content on the time step before the ejection occurred, the second plot shows these parameters immediately after the ejection occurred

Figure 1 shows the effect of the latent-heat transfer on the temperature profile. The profile shows a  dent at the diurnal skin depth and another dent where the sublimation front of the water ice is located. Figure 2 shows the effect of a dust-ejection event on the temperature profile. We will discuss these effects on the temperature profile as well as different applications of the model. Furthermore, we will discuss the way our model includes the microphysical properties of the pebbles and the latent-heat transfer numerically. 



Skorov Y. et al., 2020, MNRAS 494.3, pp. 3310-3316.

Combi M. R. et al., 2020, The Planetary Science Journal 1.3, p. 72.

Blum J. et al., 2017, MNRAS 469, pp. 755-773.

Gundlach B. and Blum J., 2012, Icarus 219.2, pp. 618-629.

Gundlach B., Fulle M. and Blum J., 2020, MNRAS 493.3, pp. 3690-3715.

How to cite: Schuckart, C., Bischoff, D., Gundlach, B., Bürger, J., Attree, N., and Blum, J.: Thermophysical modelling of comets: one dimensional microphysical modelling of cometary activity, Europlanet Science Congress 2022, Granada, Spain, 18–23 Sep 2022, EPSC2022-148,, 2022.

Sunny Laddha, Wolfgang Macher, Günter Kargl, Stephan Zivithal, Antoine Pommerol, Jürgen Blum, and Bastian Gundlach and the CoPhyLab Team


After the visit of comet 1P/Halley by the Giotto mission, “KOSI” (comet simulation) – the first large scale laboratory campaign for cometary research – investigated various physical phenomena occurring on comets by using analogue materials. The experiments were also accompanied by intensive modeling activity, which yielded several thermophysical models for the description of comets [1]. Following the revolutionary findings of the Rosetta mission, however, a new laboratory campaign became necessary to adapt existing models and to incorporate the new knowledge acquired. Thus, the “CoPhyLab – Comet Physics Laboratory” project was initiated in 2018 by an international consortium (

In the framework of this project so far, we have used several individual experiments to characterize various materials as potential ingredients for a new cometary analogue. In parallel, the construction of a large cryogenic vacuum chamber (“L-Chamber”) was completed, enabling the measurement of refractory-volatile mixtures by multiple instruments simultaneously. While cometary activity involves various physical processes, such as sublimation and gas diffusion in porous materials, understanding the comet’s thermal behavior is crucial. A rigorous thermophysical model (TPM) that considers the porous structure of the cometary material, as well as the volatile phases, and processes connected to them (sublimation, phase changes etc.), would significantly improve the understanding of the evolution of comets.

Here, we present a strategy for the development of a new TPM that in the first phase is able to describe the laboratory experiments and ultimately can be scaled to cometary environments to assist the interpretation of data collected by observations and space missions.



Previous modelling activity related to the “KOSI”-campaign was successful in the sense that the experiments could be well described by the models [2, 3]. However, a large variation of parameters between individual experiments impeded a proper understanding of phenomena that occur differently on comets than in the laboratory. Therefore, we pay special attention in the CoPhyLab campaign to repeating the most important baseline experiments. Moreover, we evolve the experiments from a simple form, where well characterized albeit idealized samples are used (e.g. glass beads), to more advanced iterations with more realistic analogue materials. On the one hand, this is to ensure that we understand the fundamental physical processes before considering complex interactions. On the other hand, the simple experiments provide us with the functional dependencies of thermophysical parameters (e.g. thermal conductivity) on the material properties, such as porosity, grain size distribution, grain shape and temperature.


Numerical Model

We approach the TPM analogously to the experiments, by starting with a very simplified macrophysical model using the finite element method (FEM). Thereby, the material-specific input parameters, which contain the microphysical relations, are taken from the characterization experiments. While most TPMs for comets are one-dimensional [2, 4], the consideration of special spatial features (e.g. sample boundaries, sensor hardware etc.) may require a 3D model [5]. To this end, we use the commercial FEM simulation software “COMSOL Multiphysics”, as it is ideally suited for our application. After the verification of our basic TPM with COMSOL, we will translate and evolve the model as a 1D code in python. This will provide a complementary, resource efficient method to simulate processes such as material ejection.


First Phase

In the first phase, we will establish the mathematical model that is required to describe the temperature evolution in a dry, well-defined sample (e.g. glass beads), illuminated in vacuum. A schematic of the corresponding experiment is shown in Figure 1. After satisfactory agreement between the model and experiment, more complex iterations will follow, for example by adding volatile phases and modelling their sublimation.

These investigations lead to a better understanding of effects in the cometary analogue material during measurements, and in addition, of effects related to the measurement setup and its interaction with the sample. As an extra benefit, a higher accuracy of future laboratory experiments in this context can be achieved.


This work is carried out in the framework of the CoPhyLab project funded by the D-A-CH programme (DFG GU 1620/3-1 and BL 298/26-1 / SNF 200021E 177964 / FWF I 3730-N36).



[1] Sears D. W. G. et al.: Laboratory simulation of the physical processes occurring on and near the
surface of comet nuclei
, Meteoritics & Planet. Sci., 34, pp. 497-525, 1999.

[2] Spohn, T. and Benkhoff, J.: Thermal history models for KOSI sublimation experiments, Icarus Vol. 87, pp. 358-371.

[3] Benkhoff, J. and Spohn, T.: Thermal histories of the KOSI samples, Geophys. Res. Letters Vol. 18, pp. 261-264, 1991.

[4] Steiner, G. and Kömle, N. I.: A model of the thermal conductivity of porous water ice at low gas pressures, Planet. Space Sci. Vol. 39, pp. 507-513, 1991.

[5] Macher, W et al.: 3D thermal modeling of two selected regions on comet 67P and comparison with Rosetta/MIRO measurements, Astron. Astrophys., Vol. 630, id. A12, 2019.

How to cite: Laddha, S., Macher, W., Kargl, G., Zivithal, S., Pommerol, A., Blum, J., and Gundlach, B. and the CoPhyLab Team: Thermophysical Modelling of the CoPhyLab Experiments, Europlanet Science Congress 2022, Granada, Spain, 18–23 Sep 2022, EPSC2022-528,, 2022.

Matthias Läuter, Tobias Kramer, Martin Rubin, and Kathrin Altwegg

During the apparition of comet 67P/Churyumov-Gerasimenko (67P/C-G) solar irradiation causes varying rates for sublimation of volatile species from the cometary nucleus. Because sublimation processes take place close to the cometary surface, the relative abundance of volatiles in the coma and the ice composition are related to each other. To quantify this relation we assume a model for the expansion of a collisionless gas from the surface into the surrounding space. We use an inverse model approach to relate the in situ measurements of gas densities from the two Rosetta instruments COPS (COmet Pressure Sensor) and DFMS (Double Focusing Mass Spectrometer) at the positions of the spacecraft to the locations of surface gas emissions during the Rosetta mission 2014-2016. We assume the temporally integrated gas emissions to be representative for the ice composition close to the surface. Our analysis shows characteristic differences in the ice compositions between both hemispheres of 67P/C-G. In particular CO2 ice has a reduced abundance on the northern hemisphere. In contrast to the hemispherical differences, the two lobes do not show significant differences in terms of their ice composition.

How to cite: Läuter, M., Kramer, T., Rubin, M., and Altwegg, K.: Determination of the ice composition near the surface of comet 67P/Churyumov-Gerasimenko, Europlanet Science Congress 2022, Granada, Spain, 18–23 Sep 2022, EPSC2022-826,, 2022.

Investigations on the formation of organic species in the gas phase cometary coma
Sana Ahmed and Kinsuk Acharyya
Fabiola Antonietta Gerosa, Héloïse Méheut, and Jérémie Bec


Planetary systems form from sub-micron sized dust embedded in a mostly gaseous protoplanetary disk. Dust has to grow into planetesimals (∼ km sized objects) considered to be the building blocks of planets. We see reminiscences of these objects in the solar system small bodies, such as asteroids and comets, that have been left over from the planets genesis process. However, planetesimal formation is still one of the major open questions in planet formation theory. Solids can’t grow up to asteroid size relying only on sticking after pairwise collisions, due to the fragmentation barrier and the drift barrier [1]. The favoured solution is to form dense particle clumps, with low velocity-dispersion, that can then collapse under self-gravity. However the way in which the solids can be concentrated, e.g. streaming instability [2] or Rossby wave instability [3], is highly debated and can be modelled as a turbulent dynamic in the disk. To understand planetesimal formation it is then fundamental to study jointly the turbulent dynamics. In this context, we model the dynamics of particles in turbulent flows with Keplerian rotation and shear. To address this astrophysical problem we use fluid-dynamics methods to provide innovative perspectives on this challenging question.


We model the dynamical evolution of the protoplanetary disk’s dust coupled to the gas phase. We perform 2D direct numerical simulations using a pseudo-spectral solver and the shearing box approach, in which the incompressible Navier-Stokes equation is solved locally with shear periodic boundary conditions. The gas is maintained turbulent through an external stochastic forcing. We use an Eulerian approach for the fluid while we implement the dust as Lagrangian particles. We explore various values of the rotation frequency Ω and the solid response time τ , a parameter related to the particle size. We analyse the results using tools borrowed from the study of dynamical systems.

Fractal and strong clustering of solids

In order to characterize the dust dynamics in the flow and to identify in which case there is a strong concentration of the solids that could lead to small body formation, the Lyapunov dimension dL is calculated for each run (Fig. 1a). This quantity gives an estimation of the fractal attractor dimension in the phase space [4]. We find three different regimes. At low rotation rate and for large particle size we obtain dL > 2, therefore the inertial particles fill the whole space. Focusing instead on intermediate values of the response time, for small rotation rates the particles are expelled from the eddies and form fractal structures (Fig. 1b), while they tend to concentrate inside the anticyclones for larger Ω (Fig. 1c). Particles eventually form a pointwise cluster for dL = 0 (Fig. 1d).

We then inspected the mass distribution of solids in the domain. Analysing a case that falls into the strong clustering regime (Ω=64/3, τ=0.05) we find that the different clusters of particles tend to merge with time to form a unique pointwise cluster at the end of the simulation (Fig. 2a). In order to estimate the interaction strength between sets of n particles, we then calculated the fractal dimensions Dn for different values of the rotation frequency and the solid response time (Fig. 2b, Fig. 2c). This dimensions characterize the anomalous scaling typical of multifractal attractors [5]. We notice that for the strong clustering case (Fig.2b, Ω=64/3) the attractor is homogenous with its fractal dimension Dn = 0 for all n. We instead find that when particles cluster fractally the attractors are multifractals that show decreasing values of Dn with increasing n. This behaviour at large n is particularly interesting for the gravitational collapse of dust: it means that particles in numerous groups will interact strongly, therefore possibly collapsing under self-gravity to form a planetesimal.

Conclusions and future perspectives

We have identified promising tools for the understanding of planetesimal formation. We also showed how the behaviour of dust particles in a turbulent flow drastically change with increasing rotation rate, leading to the formation of strong dust clumps. In the future the back-reaction from dust on gas will be added, particularly important for triggering the streaming instability. Self-interaction between solids particles (e.g. collisions, gravity) will also be considered to study how dust in vortices can accrete and form large aggregates.

Figure 1: a) Phase diagram of the rotation frequency Ω vs the particle response time τ . The colour-code is for the Lyapunov dimension dL. b) Snapshot of dust particles position and fluid vorticity at t=200 for Ω=4/3 and τ=0.1 c) for Ω=32/3 and τ=0.1. d) for Ω=80/3 and τ=0.1.

Figure 2: a) Path of different clusters of particles along time for Ω=64/3 and τ=0.05. b) Fractal dimension Dn as a function of the order of the fractal dimension n for different  values of rotation frequency Ω and τ=0.025 c) and τ=0.5.


[1] L. Testi, T. Birnstiel et al., Dust Evolution in Protoplanetary Disks, Protostars and Planets, IV, 339-361 (2014).
[2] A. N. Youdin and J. Goodman, Streaming Instabilities in Protoplanetary Disks, The Astrophysical Journal, 620, 459–469 (2005).
[3] R. V. E. Lovelace, S. A. Colgate, A. F. Nelson, Rossby Wave Instability of Keplerian Accretion Disks, The Astrophysical Journal, 513, 805-810 (1999).
[4] J.L. Kaplan and J.A. Yorke, Chaotic behavior of multidimensional difference equations, Proceedings on Functional Differential Equations and Approximation of Fixed Points, 730, 204-227 (1979).
[5] G. Paladin and A. Vulpiani, Anomalous scaling laws in multifractal objects, Physics Reports, 156, 147-225 (1987).

How to cite: Gerosa, F. A., Méheut, H., and Bec, J.: New approach to planetesimal formation: clusters of heavy particles in two-dimensional Keplerian turbulence, Europlanet Science Congress 2022, Granada, Spain, 18–23 Sep 2022, EPSC2022-462,, 2022.

Marc Brouwers and Rico Visser

Many objects that form via a gravitational collapse or contraction appear to rotate around their own axis (spin) in a manner that aligns with their orbit around larger parent structures. Systematic prograde rotation is found in the planets and larger asteroids of the Solar System [1], in the stars of several open clusters [2], as well as in the molecular clouds of different galaxies [3,4]. Remarkably, the fact that such a spin-orbit alignment exists across different scales has received little scientific attention in the last decades, while a satisfactory answer remains to be found.

Indeed, when the rotational direction of an object is assumed to be equal to that of the cloud from which it formed, only specific birth environments with rising orbital velocity curves lead to systematic prograde spin [5]. We present a new mechanism that describes how collections of particles or clouds gain a prograde rotational component when they collapse or contract while subject to an external, central force. We show that because particles that orbit in any non-rigid cloud shear away from one another over time – and do so on curved paths – their combined center-of-mass moves inward. The orbital angular momentum that is thus liberated, adds a prograde component to the spin of the object that forms. 

We first visualize the mechanism of prograde spin-up in Fig. 1a, where we plot the motion of a disk of particles on circular orbits around a central mass. For illustrative purposes, mutual interactions like self-gravity are not included here. Over time, the shear over curved orbits morphs the particles into an increasingly extended arc, whose combined center-of-mass moves inward to an orbit with a reduced semi-major axis. Hence, if the cloud of particles were to gravitationally collapse, the object forms at an orbit interior to the cloud, with a lower orbital angular momentum. The deficit goes into rotation, and the object attains a prograde spin. Being a geometric effect, the mechanism of prograde spin-up persists around any central potential that triggers shear, even those where the shear is strongly retrograde.

Fig. 1: Visual representation of the prograde spin-up mechanism

To verify this mechanism of prograde spin-up, we perform N-body simulations of gravitationally bound, collapsing clouds with REBOUND [6]. As seen in Fig. 1b, the rotational angular momentum of the cloud increases quadratically during the collapse (δLrot/LH ∝ t2), before slowing down when the collapse completes around t ∼0.8 tff in this example. We provide further visualization in Fig. 1c, where shear is seen to initially deform the cloud into a curved, bar-like structure whose center moves inward, spinning up the material. As the cloud contracts, collisions prevent further relative motion and eventually a stable, prograde binary forms.