# Schumann Resonances as a tool to constrain the depth of Titan's buried water ocean: Re-assessment of Huygens observations and preparation of the EFIELD/Dragonfly experiment

Paul Lagouanelle, Alice Le Gall

► To cite this version:

Paul Lagouanelle, Alice Le Gall. Schumann Resonances as a tool to constrain the depth of Titan's buried water ocean: Re-assessment of Huygens observations and preparation of the EFIELD/Dragonfly experiment. Icarus, 2025, 428 (March), pp.116372. 10.1016/j.icarus.2024.116372. insu-04683340v3

HAL Id: insu-04683340

https://insu.hal.science/insu-04683340v3

Submitted on 15 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Research Paper

# Schumann Resonances as a tool to constrain the depth of Titan's buried water ocean: Re-assessment of Huygens observations and preparation of the EFIELD/Dragonfly experiment

Paul Lagouanelle *, Alice Le Gall

LATMOS/IPSL, Université Paris-Saclay, UVSQ, Sorbonne Université, CNRS, 78280, Guyancourt, France

<figure><img src="image_3.png" /></figure>

## ARTICLE INFO

**Keywords:**

Titan, interior
Instrumentation
Ionospheres
Radio observations
Titan, atmosphere

## ABSTRACT

Among the lines of evidence for a buried ocean on Titan is the possible detection of Schumann-like Resonances (SR), in 2005, by the Permittivity, Wave and Altimetry (PWA) analyzer on board the ESA Huygens probe. SR are Extremely Low Frequency electromagnetic waves resonating between two electrically conductive layers. On Titan, it has been proposed that they propagate between the moon's ionosphere and a salty subsurface water ocean. Their characterization by electric field sensors can provide constraints on Titan's cavity characteristics and in particular on the depth of Titan's ocean which is key to better assess Titan's habitability. For this work we have developed a numerical model of Titan's electromagnetic cavity as well as a surrogate model (i.e., an approximate mathematical model) able to accurately approximate the behavior of the cavity. This surrogate model can be used to conduct simulations and sensitivity analyses at a low computational cost. It is used both to re-assess PWA/Huygens measurements and to predict the future performance of the EFIELD experiment on board the NASA Dragonfly mission. We demonstrate that the PWA/Huygens measurements, in particular due to their low spectral resolution, do not bring any meaningful constraint on Titan's ocean depth. On the other hand, the finer resolution of the EFIELD experiment and its ability to capture several harmonics of SR should provide more robust constraints on Titan's internal structure, especially if the electrical properties of the ice crust and the atmosphere can be better constrained.

## 1. Introduction

Several lines of evidence point to the presence of a global water ice ocean in Titan's interior. The strongest evidence arises from the investigation of the tidal variations of Titan's gravity fields inferred from Cassini flybys of the satellite (Iess et al., 2012; Durante et al., 2019). Indeed, the tidal Love number $k_2$ of 0.62 derived by Durante et al. (2019) is compatible with a high-density ocean while a recent re-assessment of $k_2$ (0.375) points to a low-density water or ammonia ocean (Goossens et al., 2024). Titan's measured obliquity of ~0.3° (Stiles et al., 2008; Meriggiola et al., 2016) is also significantly larger than the value expected for an entirely solid object and therefore suggests a decoupling between the outer ice shell and the interior of Titan (Baland et al., 2011, 2014; Bills and Nimmo, 2008, 2011). Based on both the values of Durante et al. (2019)'s $k_2$ and of the obliquity, Baland et al. (2014) estimate that the outer ice shell of Titan is at least 40 km and at most 170 km thick consistent with the results published in Kronrod et al. (2020) which test a wide range of internal structure models for Titan including thermal considerations.

Another, possible evidence for an internal ocean on Titan is the detection of ELF (Extremely Low Frequency) waves by the PWA/HASI (Permittivity, Waves and Altimetry analyzer, part of the Huygens Atmospheric Structure Instrument) experiment on board the Huygens interpreted as Schumann Resonances (Béghin et al., 2012). Schumann resonances (SR) are a set of ELF (from a few Hz up to 100 Hz) electromagnetic propagation modes that can develop in a planetary cavity excited by a broadband electromagnetic source (Schumann, 1952). On Earth, these modes are generated by lightning discharges and propagate between the ionosphere and the surface. Theoretically, SR could be observed on other planets and serve as a tool to obtain information on the planetary cavities, in particular on their dimensions (Simoes, 2007; Simões et al., 2008a,b).

If SR exist on Titan, they are probably not triggered by lightning as such activity is not expected to be common on Titan (Lorenz, 1997) and was never observed during the Cassini mission, nor by the Radio and Plasma Wave Science (RPWS) experiment (Fischer et al., 2007, 2020; Fischer and Gurnett, 2011) nor by PWA/HASI (Hamelin et al., 2009).

<footnote>* Corresponding author.
E-mail addresses: paul.lagouanelle@latmos.ipsl.fr (P. Lagouanelle), alice.legall@latmos.ipsl.fr (A. Le Gall).</footnote>
<figure><img src="image_4.png" /><figcaption>Fig. 1. Structure and parameters of Titan's cavity.</figcaption></figure>

In addition, the surface of Titan being very poorly conductive (Grard et al., 2006; Hamelin et al., 2016), it cannot act as the lower boundary of the resonant cavity which instead must be an internal electrically conductive layer. Béghin et al. (2012) propose that SR on Titan could be excited by interactions with Saturn's magnetosphere and that the signal detected at ~ 36 Hz by PWA/HASI is the second harmonic of a SR propagating between Titan's fully ionized atmospheric layer (at ~ 60–70 km altitude) and a buried salty ocean lying at a depth encompassed between 40 and 80 km. However, this interpretation is still debated as the 36 Hz line may actually be due to mechanical oscillations of the booms on which the PWA/HASI electrodes were installed or of other parts of the Huygens probe (Lorenz and Le Gall, 2020).

Nevertheless, if SR occur on Titan, their detection and characterization would place new and more robust constraints on the buried ocean. In particular, knowing more accurately its depth is key to estimate the likeliness of exchange between the ocean and the surface and therefore to assess Titan's habitability and astrobiological potential. That is the reason why the forthcoming mission to Titan, Dragonfly (NASA), will embark sensors to measure the time-varying electrical field, namely the EFIELD experiment which is part of the DraGMet (Dragonfly Geophysics and Meteorology) package (Barnes et al., 2021). Using two spherical electrodes mounted at different locations on the body of the Dragonfly drone, EFIELD will passively record the AC electrical field between ~ 5 and 100 Hz with a much finer spectral resolution than PWA/HASI.

In this paper, we describe the numerical model we have developed to simulate Titan's electromagnetic cavity and predict its resonant frequencies and associated quality factors (Section 2). This model is used to build a much less computationally expensive surrogate model which allows to perform an accurate sensitivity analysis of the SR characteristics to the cavity parameters. The surrogate model is then used to re-examine the PWA/HASI measurements leading to results very different from the ones published in Béghin et al. (2012) (Section 3). In Section 4, it is used to investigate the expected performance of the EFIELD/DraGMet experiment. Lastly, we conclude and discuss the implications of this work in Section 5.

## 2. Modeling Titan's resonant cavity

Following Simões et al. (2007, 2008a,b) who developed cavity models for Titan, Venus and other planetary environments, we used the COMSOL Multiphysics© tool to build a numerical model of electromagnetic wave propagation in the cavity of Titan. We then used

the numerical model to construct a surrogate model of the propagation of SR on Titan in order to conduct eigenfrequency analysis at a low computing and memory cost.

### 2.1. Numerical model

#### 2.1.1. Cavity description and parameters

The numerical code solves the Maxwell's equations in a spherical structure made of discrete slabs. Fig. 1 displays the simplified structure of Titan's cavity we considered; it consists of three concentric layers:

* the atmosphere/ionosphere layer for which an analytic conductivity profile is given to the model, namely the one proposed by Béghin et al. (2012) or, more recently, by Lorenz (2021) displayed in Fig. 2. Both conductivity models include a fully ionized layer at an altitude of about 60–70 km on which ELF waves are reflected. They rely on Huygens measurements of the electron-density performed from an altitude of 140 km down to the surface (Grard et al., 2006). In Lorenz (2021), the conductivity profile is interpolated from 140 to 750 km and better respect the upper limit of the near-surface conductivity imposed by Huygens Relaxation Probe measurements.

* the ice crust layer in which ELF waves are refracted following Fresnel's laws. This layer is assumed uniform in terms of electrical properties with a very small conductivity that allows ELF waves to propagate over a very long path (i.e., the skin depth is > 1000 km).

* the salty ocean layer which is assumed to be a perfect electric conductor and therefore on which ELF waves are fully reflected.

The parameters of the cavity model considered for parametric analysis are the followings:

* the thickness of the ice crust $z_c$: Based on Cassini and Huygens observations as well as on gravity and thermal modeling (see Section 1), we consider that $z_c$ can vary over a wide range of values from 5 to 200 km, and most likely between 40 and 170 km.

* the real part of the ice crust relative permittivity $\epsilon_c$: The relative permittivity of water ice at Titan's temperatures and ELF frequencies is ~ 3 (e.g., Mattei et al. (2014)). However, the crust permittivity also depends on its porosity and on the presence of impurities or contaminants such as ammonia. We therefore test values in the range 2–4 which encompasses the value measured at the Huygens landing site by the permittivity probe PWA-MIP/HASI, namely 2.5 ± 0.3 (Hamelin et al., 2016).

* the electrical conductivity of the ice crust $\sigma_c$: Through expected to be small, the conductivity of the ice crust of Titan remains uncertain. It is especially sensitive to the possible presence of ionic contaminants. Béghin et al. (2012) investigates the 1–4 nS m⁻¹ range while Hamelin et al. (2016) found a conductivity of 1.2 ± 0.6 nS m⁻¹ at the Huygens landing site. To account for this measurement, we consider values in the range 0.6–4 nS m⁻¹.

As outputs, the numerical model computes the eigenfrequencies of the first three modes of the resonant cavity ($f_1, f_2, f_3$) along with their corresponding Q-factors (quality factors: $Q_1, Q_2, Q_3$) which describe wave attenuation in the cavity. More specifically, the model provides the complex frequencies of the different eigenmodes from which the Q-factor is computed as follows:

$$Q_n = \frac{\mathrm{Re}(f_n)}{2\mathrm{Im}(f_n)} \approx \frac{f_n}{\Delta f_n} \qquad (1)$$

where Re and Im are respectively the real and imaginary parts of the complex eigenfrequency, $f_n$ is the peak power frequency of mode n, and $\Delta f_n$ is the width at half-power.
<caption>**Table 1**</caption>

Comparison of results with the 2D axi-symmetric approximation model with the 3D model for two different cavities (and associate relative error).

<table><thead><tr><th>Model</th><th>Q<sub>1</sub></th><th>ΔQ<sub>1</sub>/Q<sub>1</sub></th><th>f<sub>1</sub></th><th>Δf<sub>1</sub>/f<sub>1</sub></th></tr></thead><tbody><tr><td>2D axi-symmetric</td><td>3.20</td><td>4.69%</td><td>32.61 Hz</td><td>0.797%</td></tr><tr><td>3D</td><td>3.05</td><td>4.92%</td><td>32.35 Hz</td><td>0.803%</td></tr></tbody></table>

(a) First case: $z_c = 20 \text{ km}, \epsilon_c = 3.94, \sigma_c = 9.69 \times 10^{-9} \text{ S m}^{-1}$

<table><thead><tr><th>Q<sub>1</sub></th><th>ΔQ<sub>1</sub>/Q<sub>1</sub></th><th>f<sub>1</sub></th><th>Δf<sub>1</sub>/f<sub>1</sub></th></tr></thead><tbody><tr><td>2.735</td><td>2.93%</td><td>18.28 Hz</td><td>1.70%</td></tr><tr><td>2.655</td><td>3.01%</td><td>17.97 Hz</td><td>1.73%</td></tr></tbody></table>

(b) Second case: $z_c = 15 \text{ km}, \epsilon_c = 2.5, \sigma_c = 1 \times 10^{-9} \text{ S m}^{-1}$

<caption>**Table 2**</caption>

Comparison of results with the 2D axi-symmetric approximation model with the model used in Simoes (2007) (relative error $\epsilon_c$).

<table><thead><tr><th>Model</th><th>f<sub>1</sub></th><th>Δf<sub>1</sub>/f<sub>1</sub></th><th>f<sub>2</sub></th><th>Δf<sub>2</sub>/f<sub>2</sub></th><th>f<sub>3</sub></th><th>Δf<sub>3</sub>/f<sub>3</sub></th></tr></thead><tbody><tr><td>2D axi-symmetric</td><td>22.54 Hz</td><td>1.06%</td><td>39.08 Hz</td><td>1.20%</td><td>55.27 Hz</td><td>1.23%</td></tr><tr><td>Simoes (2007)</td><td>22.30 Hz</td><td>1.08%</td><td>38.61 Hz</td><td>1.22%</td><td>54.59 Hz</td><td>1.25%</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Input conductivity profiles of the atmosphere/ionosphere layer used in the numerical model of Titan's cavity from Béghin et al. (2012) (red) and Lorenz (2021) (blue).</figcaption></figure>

For a given uncertainty $\delta f_n$ on $f_n$ and $\Delta f_n$, which in practice is mainly dictated by the measurement spectral resolution, the uncertainty $\delta Q_n$ on the corresponding quality factor $Q_n$ can be derived by logarithmic differentiation of Eq. (1):

$$ \delta Q_n \approx Q_n(1+Q_n) \frac{\delta f_n}{f_n} \qquad (2) $$

### 2.1.2. Numerical approach

The numerical model uses the Finite Element Method (Zimmerman, 2006) for solving Maxwell's equations with the boundary conditions and layers properties as described above. Since layers properties are only functions of the radial distance, the resonant cavity problem can be solved in a 2D axi-symmetric configuration. We nevertheless validate our 2D model with comparison to a 3D model and results from the numerical model from Simoes (2007).

An example of a 3D model of Titan's cavity is displayed on Fig. 3 as well as a 2D cut of the mesh. Due to the level of discretization needed to accurately reproduce the behavior of the electric field in Titan's atmosphere, the complete 3D mesh consists of 774,258 domain elements, 181,100 boundary elements, and 2588 edge elements. A single resolution of Maxwell's equations using such a mesh takes ~2h (on an Intel Core i5-12500H, 2.5 GHz, 32 GB of RAM). In contrast, the design of a 2D axi-symmetric model with a mesh composed of 202,993 domain elements and 3,652 boundary elements requires ~30s of computation time which is much more reasonable for the purpose of using the model to perform an accurate parametric inversion.

Table 1 reports the results from the 3D model and the 2D approximation for two different cavities. The 2D approximation is accurate

enough so that the results from the 3D model are reproduced with relative errors smaller than 5% for both the Q-factor and the resonant frequency. Table 2 reports the results from the 2D approximation and the model from Simoes (2007) showing that the 2D approximation is able to accurately reproduce the cavity behavior in the case of study case with a relative error on the first three resonant frequencies smaller than 2%. Thus, for the remainder of the paper, we only consider the 2D axi-symmetrical model of Titan's planetary cavity.

### 2.2. Surrogate model

Even with the 2D axi-symmetrical approximation, a complete sensitivity analysis or data inversion using Titan's cavity numerical model would require to compute several thousand of datapoints. This would represent a significant amount of computation time that is not suitable even for a supercomputer. To avoid this, we built a surrogate model (also called “metamodel”) which provides an approximate description of the behavior of Titan's resonant cavity by analytical functions. The resulting metamodel can then be used, instead of the numerical model, to perform various analysis (e.g. optimization, sensitivity analysis) at a low computational cost (Van Steenkiste et al., 2016).

#### 2.2.1. Description

The surrogate model used here is an exact interpolator; the Polynomial-Chaos based Kriging (PCK) (Schobi et al., 2015) metamodeling combines both Polynomial Chaos Expansion (PCE) and Kriging to predict the variations of a given model $M(X)$. Kriging is used to interpolate the local variations of the output model while PCE is useful for the global approximation. A PCK metamodel is defined by:

$$ M(X) = \sum_{\alpha \in \mathcal{A}} y_{\alpha} \psi_{\alpha}(X) + \sigma^2 Z(X, \omega) \qquad (3) $$

where $\sum_{\alpha \in \mathcal{A}} y_{\alpha} \psi_{\alpha}(X)$ is a weighted sum of orthonormal polynomials describing the trend of the PCK model and $\sigma^2 Z(X, \omega)$ is a zero-mean stationary Gaussian process with a variance of $\sigma^2$. The computation of the metamodel parameters is performed by the UQLab framework available on Matlab (Marelli and Sudret, 2014).

The costly part of the metamodeling process is the training time which can be greatly reduced by using sequential sampling instead of classical space-filling approaches. In this work an adaptive sampling algorithm combined with PCK has been used to build the metamodel. This algorithm has already been proven useful for various electromagnetic problems (Lagouanelle et al., 2023).

#### 2.2.2. Metamodel accuracy

Once built, the surrogate model is used to predict the behavior of the cavity outside of the training data. For such a purpose, a proper metric is crucial to quantify the accuracy of these predictions. A classical approach consists in using a validation dataset outside of the training dataset and computing the mean squared error (MSE) of the metamodel prediction compared to the real input values. However, the resulting MSE is biased by the use of only one dataset and could vary greatly from one validation dataset to another. Moreover, this approach
<figure><img src="image_6.png" /><figcaption>Fig. 3. (a) 3D model of Titan's planetary cavity on COMSOL Multiphysics, (b) 2D cut of the 3D mesh on (a).</figcaption></figure>

requires additional calls of the expensive computational numerical model for building the training dataset, which ultimately increases the computation time. Thus, a better metric was chosen : the Leave-one-out cross-validation error (LOO), which does not require additional computations.

Let us consider a set $\{(X_1, Y_1), \dots, (X_N, Y_N)\}$ of $N$ input samples. Using this set, one can build a PCK metamodel $\mathcal{M}$ and evaluate the LOO as follows:

$$ LOO = \frac{1}{N} \sum_{i=1}^{N} \left( \frac{\|M_{i}(X_i) - Y_i\|}{\|Y_i\|} \right)^2 \qquad (4) $$

where $M_{i}$ is the mean predictor that was trained using all $(X, Y)$ except $(X_i, Y_i)$. For a given datapoint $(X_i, Y_i)$, a metamodel is built with all datapoints except datapoint $i$, which gives $N-1$ training datapoints. This metamodel is then used to predict the value $Y_i$ at the remaining datapoint $i$, where the difference is classified with a MSE. The process is repeated for every datapoint which, after average, provides the LOO. The use of $N$ different validation sets of one datapoint guarantees that the LOO is much less biased than a classical MSE and reduces the probability of overestimating the validation error (Elisseff et al., 2003). In this work, we therefore consider the LOO as our accuracy metric. A LOO close to 1 (or 100%) implies that the surrogate model does not provide a good approximation of the system. On the other hand, the smaller the LOO, the more accurate the surrogate model.

### 2.2.3. Sensitivity indices

Since the resulting metamodel consists of an analytical function, calling the metamodel is extremely cheap in terms of computation time. Thus, sensitivity analyses, which are usually performed by Monte-Carlo analyses over the parameter spaces, are now feasible at a low computation cost.

The sensitivity analysis we conduct relies on Sobol' indices which are scalars between 0 and 1 describing the influence of a set of inputs on a model output (Sobol, 1993). The most commonly used Sobol' indices is the first-order Sobol' index defined, for a given parameter $P_i$, as:

$$ S_i = \frac{\mathrm{Var}_{P_i}(\mathbb{E}_{X_i}(Y|X_i))}{\mathrm{Var}(Y)} \qquad (5) $$

$S_i$ is a measure of the fraction of the output variance caused by the variance of a given input parameter. In other words, it describes the impact of a parameter $P_i$ alone on the output model compared to other parameters. The closer to 1, the bigger impact $P_i$ has on the model output.

However, parameters are usually not independent and their relative effects cannot be separated from each other. This leads to the definition of higher-order Sobol' indices as, for a subset of parameters $(P_{i_1}, \dots, P_{i_s})$:

$$ S_{i_1, \dots, i_s} = \frac{\mathrm{Var}_{(P_{i_1}, \dots, P_{i_s})}(\mathbb{E}_{X_{i_1}, \dots, X_{i_s}}(Y | X_{i_1}, \dots, X_{i_s}))}{\mathrm{Var}(Y)} \qquad (6) $$

which describes the sensitivity of the model to the variations of several input parameters simultaneously.

For high dimensional output models, the interpretation of all orders Sobol' indices can be difficult due to the high number of possible combinations. Therefore, for an input parameter $P_i$, a total-effect index (or total Sobol' index) $S_i^T$ is defined by summing all the Sobol' indices as follows:

$$ S_i^T = \sum_{\{u, u \subseteq \{1, \dots, d\} \text{ and } i \in u\}} S_u \qquad (7) $$

$S_i^T$ is the most suited sensitivity tracker for our study and will be refereed to as $S_i$ in the remainder of the paper. When using PCK metamodels, the computation of the various Sobol' indices of the surrogate model can be easily extracted from the polynomial decomposition (see Eq. (3)). Therefore, no additional computation of the surrogate model are required to perform the sensitivity analysis based on the total Sobol' indices.

## 3. Re-assessment of PWA/HASI/Huygens observations

Fulchignoni et al. (2005) first reported the detection of a narrow spectral feature at ~36 Hz in the ELF spectrum measured by PWA/HASI during Huygens' descent in Titan's atmosphere in January 2005, from an altitude of 140 km down to the surface. The magnitude of this signal is especially enhanced just after the deployment of the stabilizer parachute, at an altitude of ~110 km. Béghin et al. (2007) proposed different scenarios, both natural and artificial, to explain the 36 Hz signal. In Béghin et al. (2012), a natural scenario is preferred: the signal and associated Q-factor of about ~6 would be the second harmonic of a SR propagating between Titan's ionosphere and ocean and triggered by interactions with Saturn's magnetosphere. Using an approximate analytical model of Titan's cavity, Béghin et al. (2012) further derive constraints on the physical parameters of the cavity from PWA/HASI measurements. More specifically, they conclude that the measured $f_2 = 36 \pm 3$ Hz and $Q_2 \sim 6$ are indicative of a water-ammonia ocean lying at a depth of 40–80 km.

In this section, we re-assess the PWA/HASI data using the surrogate model we have developed (see Section 2) to investigate, in a more accurate fashion, the constraints Huygens measurements bring on the thickness of the ice crust (i.e., the depth of the ocean) $z_c$. As a starting point, we adopt exactly the same hypotheses as in Béghin et al. (2012) that is the same conductivity profile in the atmosphere and ranges of variation for $\epsilon_c$ ($2-4$) and $\sigma_c$ ($1-4$ nS m$^{-1}$).

The metamodeling process estimates consistently (LOO ≈ 3.9%, $n_{samples} = 1584$) the second harmonic of SR ($f_2, Q_2$). Using the resulting surrogate model, two regular 3D grids of $f_2$ and $Q_2$ values along all three input parameters $z_c, \epsilon_c$ and $\sigma_c$ can be computed at a low computation cost: $50 \times 50 \times 50 = 125,000$ values, which would have
<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Inversion of Huygens PWA/HASI data taking the same hypotheses as in Béghin et al. (2012) (i.e., conductivity profile, $ε_c ∈ [2–4]$ and $σ_c ∈ [1–4]$ nS m⁻¹) and using the surrogate model developed in this work (a). All values of $z_c$ can be in a combination with $ε_c$ and $σ_c$ that reproduces Huygens data which were measured with a 3 Hz spectral resolution. A 1 Hz spectral resolution would drastically reduce the parameter space for $z_c$. (b) displays the same inversion exercise but considering that the detected line at 36 Hz is the fundamental of the SR instead of the second harmonic (no inversion solution for a 1 Hz spectral resolution in this case).</figcaption></figure>

taken 43 days of computation time using directly the numerical model instead of only 13 hours for training the metamodel.

By analyzing the 3D grids, every combination of $z_c$, $ε_c$ and $σ_c$ which gives $f_2 ∈ [33–39]$ Hz and $Q_2 ∈ [3–9]$ can be classified as a potential solution explaining Huygens measurements. Fig. 4(a) displays all the potential inversion of Huygens measurements in the plane ($z_c$, $ε_c$) where all the grids in the direction $σ_c$ have been stacked. When accounting for 3 Hz spectral resolution, the range of possible values for $z_c$ covers all the parameter space (i.e. 5–200 km) which means that no constraint can be deduced on $z_c$ from the PWA/HASI dataset (blue zone). On the other hand, if the signal characteristics were known with a 1 Hz resolution, only a narrow range of $z_c$ values could explain the observations (red outlined area). This hypothetical inversion would restrain the thickness of the ice crust to $z_c ∈ [165–200]$ km.

Our results are in contradiction with Béghin et al. (2012) which can be explained by the various analytical approximations they consider to solve the wave propagation equation. Furthermore, since the metamodel is sufficiently accurate (LOO ≈ 3.9%), it can be used to conduct an accurate Sobol'-based sensitivity analysis. The following total Sobol' indices are found: $S_{z_c} = 0.73$, $S_{ε_c} = 0.28$ and $S_{σ_c} = 0.74$. These indexes indicate that the ice crust thickness is the parameter that has the most significant impact on the SR characteristics. However, since all indexes are of the same order of magnitude, no parameter can be regarded as having a negligible impact. This is a further guarantee that the metamodel accuracy is correctly estimated with the LOO and that our metamodel is highly accurate. Therefore the differences from Béghin et al. (2012)'s results and the present work cannot be ascribed to an incorrect estimate of the errors from the metamodel.

The surrogate model also demonstrates that solutions (in terms of combinations of $z_c$, $ε_c$ and $σ_c$ values) in which the detected signal is not the second harmonic but the fundamental are possible (Fig. 4(b)). If so, the constraints on $z_c$ derived from PWA/HASI would be different (namely >80 km, see Fig. 4(b)).

Lastly, Fig. 5 shows how the constraints on $z_c$ are modified if the numerical model from which the surrogate model was built rather uses (Lorenz, 2021)'s conductivity profile and values of $ε_c$ and $σ_c$ from Hamelin et al. (2016). We will keep these hypotheses for the remainder of the paper and, in particular, to assess EFIELD/DraGMet future performance in the frame of the Dragonfly mission. The associated sensitivity analysis provides the following total Sobol' indices: $S_{z_c} = 0.78$, $S_{ε_c} = 0.34$ and $S_{σ_c} = 0.51$. Again and notably, $z_c$ is the parameter that have the most impact on SR characteristics. Nevertheless, the effects of $ε_c$ and $σ_c$ cannot be neglected as their respective total Sobol' index are of the same order of magnitude. This further implies

<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Inversion of Huygens PWA/HASI using the surrogate model developed in this work with up-to-date assumptions, namely the atmosphere conductivity profile from Lorenz (2021) and $ε_c ∈ [2.2–2.8]$ and $σ_c ∈ [0.6–1.8]$ nS m⁻¹ as estimated by Hamelin et al. (2016).</figcaption></figure>

that their accurate knowledge would greatly reduce the uncertainty on the inversion of Huygens measurements as well as be very valuable for the analysis of future EFIELD data. In particular, if the real part of the permittivity of the ice crust were measured as $ε_c = 2.5 ± 0.1$ instead of $2.5 ± 0.3$, it would help discriminating between the two distinctive domains in Fig. 5: $z_c ∈ [5–40]$ km and $z_c ∈ [150–200]$ km. This is further discussed below.

# 4. Anticipated EFIELD/DraGMet/Dragonfly performance

This section investigates the performance and possible outcome from the forthcoming electric-field experiment on Titan.

## 4.1. The EFIELD experiment on board Dragonfly

In June 2019, NASA selected the Dragonfly mission project for its New Frontiers program (Turtle et al., 2018; Lorenz et al., 2018; Barnes et al., 2021). The primary goal of the Dragonfly mission is to investigate
the chemistry and habitability of Titan. Starting operations in mid-2030s, the Dragonfly quadcopter drone will visit a variety of sites, from a dune field to the rim of a young impact crater, and sample materials in different geologic settings. Dragonfly includes a Geophysical and Meteorological package (DraGMet) which is a suite of sensors designed to measure e.g. the temperature, pressure, methane humidity, wind speed and direction, ground dielectric constant, thermal properties and level of seismic activity at each Dragonfly landing site.

Among these sensors is the EFIELD experiment which consists of two independent spherical electrodes (~5 cm in diameter) accommodated at the end of ~ 25 cm long stalks pointing away from the drone body. From two different locations on the drone, these electrodes will passively record the time-varying electrical field at low frequencies (~5–100 Hz) with the main goal of detecting SR, if any. As a secondary objective, the EFIELD probes will detect and characterize near-surface wind-blown charged grains flying in their vicinity (Chatain et al., 2023).

The EFIELD experiment offer many advantages over the PWA/HASI one. It will operate during an extended period of time (several times a Titan day for the 3.3 years of the nominal mission), from a stable and much mechanically-quieter platform than the Huygens probe. Further, the EFIELD design should guarantee a spectral resolution of 1 Hz (against, at best, 3 Hz for PWA/HASI) and the capture of the first three harmonics of the SR. Fig. 4(a) demonstrates how valuable a finer spectral resolution would be to bring more robust constraints on the depth of the buried ocean $z_c$. The benefit of detecting more than one SR harmonics is investigated below.

## 4.2. Multi-modal analysis

The re-assessment of Huygens data presented in Section 3 relies on the measurement of only one mode of the SR; it concludes that a wide range of values are possible for the thickness of the ice crust ($z_c \in [5-200]$ km). The multi-modal analysis of SR enabled by EFIELD should drastically reduce this range because, in addition to a better spectral resolution, the three first modes of the Schumann resonances will be associated with three different domains of possible $z_c$ values whose intersection may be narrow.

As an illustration, we numerically simulate Titan's cavity for the following parameter values: $z_c = 60$ km, $\epsilon_c = 2.5$, $\sigma_c = 1.2$ nS m⁻¹ and the up-to-date conductivity profile from Lorenz (2021). The numerical model provides the following outputs: $f_1 = 28.4$ Hz, $f_2 = 44.9$ Hz, $f_3 = 62.4$ Hz along with their corresponding quality factors: $Q_1 = 3.25$, $Q_2 = 3.58$ and $Q_3 = 3.81$. Using the metamodel, we further compute the variations of the resonant frequencies ($f_1, f_2, f_3$) and of the quality factors ($Q_1, Q_2, Q_3$) as a function of the thickness of the ice crust $z_c$, assuming fixed values for $\epsilon_c = 2.5$ and $\sigma_c = 1.2$ nS m⁻¹. These variations are displayed on Fig. 6.

Assuming a 1 Hz resolution, each harmonic can be inverted separately (using the method described in Section 3), resulting in three different possible domains for $z_c$:

$$ \bullet f_1 = 28.4 \pm 1 \text{ Hz} \Rightarrow z_c \in [19.7-21.3] \cup [56.0-73.0] \cup [81.4-85.8] \text{ km} $$

$$ \bullet f_2 = 44.9 \pm 1 \text{ Hz} \Rightarrow z_c \in [20.9-22.4] \cup [54.8-72.3] \cup [83.4-87.8] \text{ km} $$

$$ \bullet f_3 = 62.4 \pm 1 \text{ Hz} \Rightarrow z_c \in [21.7-23.5] \cup [33.2-46.3] \cup [56.7-63.2] \text{ km} $$

The intersection of these domains is: $z_c \in [56.7 - 63.2]$ km which corresponds to an uncertainty of 6% with respect to the input value of 60 km. A similar inversion is performed on the quality factors separately:

$$ \bullet Q_1 = 3.25 \pm 0.49 \Rightarrow z_c \in [13.8 - 89.4] \text{ km} $$

$$ \bullet Q_2 = 3.58 \pm 0.36 \Rightarrow z_c > 11.9 \text{ km} $$

$$ \bullet Q_3 = 3.81 \pm 0.29 \Rightarrow z_c > 14.9 \text{ km} $$

In the case considered here, no restriction can be further obtained on $z_c$ from the Q-factor values (see Fig. 6b). Nevertheless, this example well illustrates the values of measuring several modes of the SR.

## 4.3. Inversion uncertainty

The example above assumes fixed values for both $\epsilon_c$ and $\sigma_c$. To take into account our imperfect knowledge of the electrical properties of Titan's ice crust, the uncertainties on these parameters have to be propagated through the inversion to estimate their effects on the derivation of $z_c$. As another illustration, we numerically simulate Titan's cavity with a thickness of the ice crust fixed at $z_c = 60$ km. A multi-modal analysis is then performed for every point ($z_c, \epsilon_c, \sigma_c$) in a regular parameter grid with $\epsilon_c$ and $\sigma_c$ respectively varying in $2.5 \pm 0.3$ and $1.2 \pm 0.6$ nS m⁻¹ as found at the Huygens landing site by Hamelin et al. (2016). The parameter grid consists of $17 \times 25 \times 25 = 15,300$ datapoints in the 3D parameter space $\{z_c\} \times \{\epsilon_c\} \times \{\sigma_c\}$. Every datapoint on the grid which returns the desired values for the first three modes considering their measurement uncertainty of 1 Hz is a possible inversion of EFIELD.

This allows us to compute the relative frequencies with 3 km intervals of the returned values of $z_c$ assuming that the three first harmonic of the SR are detected with a resolution of 1 Hz, as Fig. 7(a) displays. Interestingly, 90% of the inversion cases fall in the range $z_c \in [51-69]$ km while $\approx 10\%$ of the cases return $z_c$ values in the ranges $[20-40]$ km or $[70-80]$ km (blue bars). Assuming a normal distribution of $z_c$ centered in 60 km, this corresponds to a standard deviation (STD) of 6.8 km. Fig. 7(b) further shows that reducing the uncertainty on $\sigma_c$ from $1.2 \pm 0.6$ nS m⁻¹ to $1.2 \pm 0.2$ nS m⁻¹ reduces the STD to the value of 5.7 km.

Such exercise was repeated for a set of $z_c$ values in the $[10-85]$ km range to produce Fig. 8 which displays the STD (in km and %) of the distribution of the inferred $z_c$ as a function of the $z_c$ value for two assumptions on the range of variation of $\sigma_c$. In most cases, the STD is smaller than 10 km. Moreover, it significantly decreases as $z_c$ increases, especially after 60 km where it becomes smaller than 10%. In contrast, for small values of $z_c$ (< 25 km), the theoretical relative error on $z_c$ can reach almost 100%. This can be partially explained by the relative coarse mesh of the numerical model; indeed the dimensions of the cells have been imposed greater than 5 km in order to reduce computation time. Further developments are required to build a surrogate model more appropriate to small thicknesses (i.e., with a finer mesh) but this is out of the scope of this paper.

The same analysis was conducted assuming a more constrained knowledge of $\sigma_c$ (namely $\sigma_c = 1.2 \pm 0.2$ nS m⁻¹) resulting in smaller STD for almost all cases (red lines in Fig. 8). In particular, the uncertainty at estimating small thicknesses drops from $\approx 80\%$ to $\approx 40\%$. This further shows the need for a more accurate knowledge of the electrical properties of the crust. This point and others are discussed in the following section.

# 5. Discussion and conclusion

For this work we have developed a numerical model of Titan's cavity to then build a less computationally expensive surrogate model able to describe how the cavity characteristics (i.e., eigenfrequency and Q-factors) vary with the main cavity parameters (i.e., Titan's ice crust thickness and electrical properties). This model (and its use for data inversion) is a powerful tool for the analysis of electric field measurements on Titan (and elsewhere). It was used to re-assess Huygens observations leading to the conclusion that the 2005 detection of a line at ~ 36 Hz, if indeed due to SR, does not provide any specific constraint on the depth of Titan's ocean in the range 5–200 km contrary to what is advanced in Béghin et al. (2012).

The surrogate model was also used to estimate the possible outcomes from the EFIELD/DraGMet/Dragonfly experiment. EFIELD is designed to detect several modes of SR with a fine spectral resolution; we have demonstrated that it has the ability to put a meaningful constraint on the thickness of the ice crust. Considering the electromagnetic properties varying in the ranges specified by Hamelin et al. (2016), the various sensitivity analysis presented throughout this work, reach the same conclusion: although the thickness of the ice crust
<figure><img src="image_9.png" /><figcaption>**Fig. 6.** Variations of the resonant frequencies ($f_1, f_2, f_3$) (a) and of the quality factors ($Q_1, Q_2, Q_3$) (b) as a function of the thickness of the ice crust $z_c$ assuming $\epsilon_c = 2.5$, $\sigma_c = 1.2 \text{ nS m}^{-1}$. The red lines embody the 1 Hz resolution and its corresponding quality factor resolution for the EFIELD measurement. In a case where the input value is $z_c = 60 \text{ km}$, the multi-modal analysis provides a range of values for $z_c$ that is between 56.7 km and 63.2 km (dashed red lines).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 7.** Relative frequencies histogram of the EFIELD inversion with 3 km intervals for an ice crust $z_c = 60 \text{ km}$ using the metamodelling process for two different uncertainties on $\sigma_c$: $\sigma_c = 1.2 \pm 0.6 \text{ nS m}^{-1}$ (a) and $\sigma_c = 1.2 \pm 0.2 \text{ nS m}^{-1}$ (b).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Standard deviation of the EFIELD inversion against the thickness of the ice crust $z_c$ for two different uncertainties on $\sigma_c$ in km (a) and percentage of $z_c$ (b)</figcaption></figure>
is the most influential parameter on SR characteristics, the resonant frequencies are dependent in the same order of magnitude on the electrical properties of Titan's crust. Therefore, in order to reduce the uncertainty at estimating Titan's crust thickness, it is crucial to reduce the uncertainty on the crust electrical properties.

Another experiment on board the Dragonfly drone will contribute to this task: the DIEL/DraGMet experiment. Acting as a mutual-impedance probe with a pair of electrodes mounted on each landing skid of the drone, DIEL will measure the complex permittivity (which includes the electrical conductivity) of the ground at several low frequencies (<10kHz) thus providing insights into the composition, moisture and porosity of the near-subsurface of Titan as well as on the spatial and temporal variations of such properties. Though all measured permittivity values may not be representative of the crust, values measured on the ejecta blanket of the geologically-young Selk crater (the final destination of Dragonfly) may be. In addition, variations of the measured complex permittivity along the Dragonfly journey to Selk and its possible correlation with otherwise inferred vicinity of the water ice bedrock in the near subsurface will provide further constraint on the ice crust electrical properties.

Nevertheless, as highlighted in Lorenz and Le Gall (2020), one of the main sources of uncertainty is and will remain our limited knowledge of the lower atmosphere conductivity structure. Unfortunately, no improvement is to be expected from forthcoming observations as Dragonfly will not perform measurements during its descent in Titan's atmosphere. As a consequence, only theoretical developments can provide further insights on the atmosphere conductivity profile and its expected variations with the local hour, solar activity and the position of Titan in Saturn's magnetosphere.

Future investigations will include the simulation of the actual EFIELD electrodes accommodated on the Dragonfly (conductive) body as well as the study of the effect on measurements of the location and polarization of the possible sources of SR. In the mid-2030s, when the Dragonfly drone will be on Titan, methane-storms are expected at the South Pole. Though still speculative, such storms might be associated with some atmospheric electric phenomena which could generate SRs that the EFIELD experiment (located at much lower latitudes) will try to detect. Given the different orientation and altitude on the drone of the two EFIELD electrodes, at least two components of the electrical field will be measured. The third component could be captured by rotating the drone and future study will also explore the value of measuring the full electrical field vector.

Lastly, the approach we have developed in this paper, based on the development and use of a surrogate model, can be readily applied to the study of SR in other bodies in the Solar System such as Venus. In particular in may help better characterize the electrical environment of Venus by comparison to Venera 11 and 12 past observations or DaVinci future electric measurements, if any.

### CRediT authorship contribution statement

**Paul Lagouanelle:** Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Methodology, Investigation, Formal analysis, Data curation, Conceptualization.
**Alice Le Gall:** Writing – review & editing, Writing – original draft, Visualization, Validation, Supervision, Software, Methodology, Investigation, Funding acquisition, Formal analysis, Data curation, Conceptualization.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgments

The authors wish to thank Ralph Lorenz (APL, PI of DraGMet/Dragonfly) and Jean-Jacques Berthelier, Audrey Chatain and Franck Montmessin from LATMOS for useful discussions on this project. The authors are also grateful to Région Ile-de-France (DIM-ACAV+ program) for funding this research.

### Data availability

Data will be made available on request.

### References

* Baland, R.-M., Tobie, G., Lefèvre, A., Van Hoolst, T., 2014. Titan's internal structure inferred from its gravity field, shape, and rotation state. Icarus 237, 29–41.

* Baland, R.-M., Van Hoolst, T., Yseboodt, M., Karatekin, Ö., 2011. Titan's obliquity as evidence of a subsurface ocean? Astron. Astrophys. 530, A141.

* Barnes, J.W., Turtle, E.P., Trainer, M.G., Lorenz, R.D., MacKenzie, S.M., Brinckerhoff, W.B., Cable, M.L., Ernst, C.M., Freissinet, C., Hand, K.P., et al., 2021. Science goals and objectives for the Dragonfly Titan rotorcraft relocatable lander. Planet. Sci. J. 2 (4), 130.

* Béghin, C., Randriamboarison, O., Hamelin, M., Karkoschka, E., Sotin, C., Whitten, R.C., Berthelier, J.-J., Grard, R., Simões, F., 2012. Analytic theory of Titan's Schumann resonance: Constraints on ionospheric conductivity and buried water ocean. Icarus 218 (2), 1028–1042.

* Béghin, C., Simões, F., Krasnoselskikh, V., Schwingenschuh, K., Berthelier, J.-J., Besser, B., Bettanini, C., Grard, R., Hamelin, M., López-Moreno, J., et al., 2007. A Schumann-like resonance on Titan driven by Saturn's magnetosphere possibly revealed by the Huygens Probe. Icarus 191 (1), 251–266.

* Bills, B.G., Nimmo, F., 2008. Forced obliquity and moments of inertia of Titan. Icarus 196 (1), 293–297.

* Bills, B.G., Nimmo, F., 2011. Rotational dynamics and internal structure of Titan. Icarus 214 (1), 351–355.

* Chatain, A., Le Gall, A., Berthelier, J.-J., Lorenz, R.D., Hassen-Khodja, R., Lebreton, J.-P., Joly-Jehenne, T., Déprez, G., 2023. Detection and characterization of wind-blown charged sand grains on titan with the DraGMet/EFIELD experiment on Dragonfly. Icarus 391, 115345.

* Durante, D., Hemingway, D., Racioppa, P., Iess, L., Stevenson, D., 2019. Titan's gravity field and interior structure after Cassini. Icarus 326, 123–132.

* Elisseff, A., Pontil, M., et al., 2003. Leave-one-out error and stability of learning algorithms with applications. In: NATO Science Series Sub Series III Computer and Systems Sciences, vol. 190. IOS Press, pp. 111–130.

* Fischer, G., Farrell, W., Gurnett, D., Kurth, W., 2020. Nondetection of radio emissions from Titan lightning by Cassini RPWS. J. Geophys. Res. Planets 125 (9), e2020JE006496.

* Fischer, G., Gurnett, D., 2011. The search for Titan lightning radio emissions. Geophys. Res. Lett. 38 (8).

* Fischer, G., Gurnett, D.A., Kurth, W.S., Farrell, W.M., Kaiser, M.L., Zarka, P., 2007. Nondetection of Titan lightning radio emissions with Cassini/RPWS after 35 close Titan flybys. Geophys. Res. Lett. 34 (22).

* Fulchignoni, M., Ferri, F., Angrilli, F., Ball, A., Bar-Nun, A., Barucci, M., Bettanini, C., Bianchini, G., Borucki, W., Colombatti, G., et al., 2005. In situ measurements of the physical characteristics of Titan's environment. Nature 438 (7069), 785–791.

* Goossens, S., van Noort, B., Mateo, A., Mazarico, E., van der Wal, W., 2024. A low-density ocean inside Titan inferred from Cassini data. Nat. Astron. 1–10.

* Grard, R., Hamelin, M., López-Moreno, J., Schwingenschuh, K., Jernej, I., Molina-Cuberos, G., Simões, F., Trautner, R., Falkner, P., Ferri, F., et al., 2006. Electric properties and related physical characteristics of the atmosphere and surface of Titan. Planet. Space Sci. 54 (12), 1124–1136.

* Hamelin, M., Grard, R., López-Moreno, J., Schwingenschuh, K., Béghin, C., Berthelier, J., Simões, F., 2009. Comment on “Evidence of electrical activity on Titan drawn from the Schumann resonances sent by Huygens probe” by J.A. Morente, J.A. Portí, A. Salinas, E.A. Navarro [2008, Icarus, 195, 802–811]. Icarus 204 (1), 349–351. http://dx.doi.org/10.1016/j.icarus.2009.01.031, URL https://www.sciencedirect.com/science/article/pii/S0019103509000748.

* Hamelin, M., Lethuillier, A., Le Gall, A., Grard, R., Béghin, C., Schwingenschuh, K., Jernej, I., López-Moreno, J.-J., Brown, V., Lorenz, R.D., et al., 2016. The electrical properties of Titan's surface at the Huygens landing site measured with the PWA-HASI mutual impedance probe. New approach and new findings. Icarus 270, 272–290.

* Iess, L., Jacobson, R.A., Ducci, M., Stevenson, D.J., Lunine, J.I., Armstrong, J.W., Asmar, S.W., Racioppa, P., Rappaport, N.J., Tortora, P., 2012. The tides of Titan. Science 337 (6093), 457–459.

* Kronrod, V., Dunaeva, A., Gudkova, T., Kuskov, O., 2020. Matching of models of the internal structure and thermal regime of partially differentiated Titan with gravity field. Sol. Syst. Res. 54 (5), 405–419.
Lagouanelle, P., Freschi, F., Pichon, L., 2023. Adaptive sampling for fast and accurate metamodel-based sensitivity analysis of complex electromagnetic problems. IEEE Trans. Electromagn. Compat..

Lorenz, R., 1997. Lightning and triboelectric charging hazard assessment for the huygens probe. In: Wilson, A. (Ed.), Huygens: Science, Payload and Mission, Proceedings of an ESA Conference. vol. 1177, p. 265.

Lorenz, R.D., 2021. The low electrical conductivity of Titan's lower atmosphere. Icarus 354, 114092.

Lorenz, R.D., Le Gall, A., 2020. Schumann resonance on Titan: A critical re-assessment. Icarus 351, 113942.

Lorenz, R.D., Turtle, E.P., Barnes, J.W., Trainer, M.G., Adams, D.S., Hibbard, K.E., Sheldon, C.Z., Zacny, K., Peplowski, P.N., Lawrence, D.J., et al., 2018. Dragonfly: A rotorcraft lander concept for scientific exploration at Titan. Johns Hopkins APL Tech. Dig. 34 (3), 14.

Marelli, S., Sudret, B., 2014. UQLab: A Framework for Uncertainty Quantification in Matlab. American Society of Civil Engineers.

Mattei, E., Lauro, S., Vannaroni, G., Cosciotti, B., Bella, F., Pettinelli, E., 2014. Dielectric measurements and radar attenuation estimation of ice/basalt sand mixtures as martian Polar Caps analogues. Icarus 229, 428–433.

Meriggiola, R., Iess, L., Stiles, B.W., Lunine, J.I., Mitri, G., 2016. The rotational dynamics of Titan from Cassini RADAR images. Icarus 275, 183–192.

Schobi, R., Sudret, B., Wiart, J., 2015. Polynomial-chaos-based Kriging. Int. J. Uncertain. Quantif. 5 (2).

Schumann, W.O., 1952. "Uber die strahlungslosen Eigenschwingungen einer leitenden Kugel, die von einer Luftschicht und einer Ionosphärenhülle umgeben ist. Z. Naturforsch. A 7 (2), 149–154.

Simoes, F., 2007. Theoretical and Experimental Studies of Electromagnetic Resonances in the Ionospheric Cavities of Planets and Satellites; Instrument and Mission Perspectives (Ph.D. thesis). Université Pierre et Marie Curie-Paris VI.

Simões, F., Grard, R., Hamelin, M., López-Moreno, J., Schwingenschuh, K., Béghin, C. Berthelier, J.-J., Besser, B., Brown, V., Chabassière, M., et al., 2007. A new numerical model for the simulation of ELF wave propagation and the computation of eigenmodes in the atmosphere of Titan: Did Huygens observe any Schumann resonance? Planet. Space Sci. 55 (13), 1978–1989.

Simões, F., Grard, R., Hamelin, M., López-Moreno, J., Schwingenschuh, K., Béghin, C., Berthelier, J.-J., Lebreton, J.-P., Molina-Cuberos, G., Tokano, T., 2008a. The Schumann resonance: A tool for exploring the atmospheric environment and the subsurface of the planets and their satellites. Icarus 194 (1), 30–41.

Simões, F., Hamelin, M., Grard, R., Aplin, K., Béghin, C., Berthelier, J.-J., Besser, B., Lebreton, J.-P., López-Moreno, J., Molina-Cuberos, G., et al., 2008b. Electromagnetic wave propagation in the surface-ionosphere cavity of Venus. J. Geophys. Res. Planets 113 (E7).

Sobol, I.M., 1993. Sensitivity analysis for non-linear mathematical models. Math. Model. Comput. Exp. 1, 407–414.

Stiles, B.W., Kirk, R.L., Lorenz, R.D., Hensley, S., Lee, E., Ostro, S.J., Allison, M.D., Callahan, P.S., Gim, Y., Iess, L., et al., 2008. Determining Titan's spin state from Cassini RADAR images. Astron. J. 135 (5), 1669.

Turtle, E., Barnes, J., Trainer, M., Lorenz, R., Hibbard, K., Adams, D., Bedini, P., Brinckerhoff, W., Cable, M., Ernst, C., et al., 2018. Dragonfly: In situ exploration of Titan's organic chemistry and habitability. In: 49th Lunar and Planetary Science Conference. vol. 49, LPI-Contribution.

Van Steenkiste, T., van der Herten, J., Couckuyt, I., Dhaene, T., 2016. Sensitivity analysis of expensive black-box systems using metamodeling. In: 2016 Winter Simulation Conference. WSC, IEEE, pp. 578–589.

Zimmerman, W.B., 2006. Multiphysics Modeling with Finite Element Methods. vol. 18, World Scientific Publishing Company.