© The Authors 2025

1 Introduction

Comets are active, small Solar System bodies. When approaching the Sun, their nuclei get heated up and volatile ices within them sublimate, dragging dust particles alongside to form detectable coma and tails (Whipple 1950). The activity of comets is an intriguing topic, as it provides insights into the formation and evolution of comets, such as the cyclical variations in the relative abundance of water ice keeping the comet “alive” (De Sanctis et al. 2015; Shi et al. 2018), while heat transport and ice sublimation help reflect the properties of the uppermost layer of cometary nucleus (Gulkis et al. 2015) and seasonal effects of ice sublimation induce various topographical changes (El-Maarry et al. 2015; Tang et al. 2025). Their activity can also help us understand the possibility that comets could have delivered water to Earth (Mandt et al. 2024).

Within ∼3 au, water vapor is the predominant species in outgassing (Hässig et al. 2015; Läuter et al. 2020). The water production rate and the dust production rate are primary indicators for characterizing cometary activity and its variations. The modeling of water and dust production requires us to simulate the thermophysical processes, taking into account the orbit and the shape of the nucleus, as well as its various physical properties. There are variations in the thermophysical models for different aims (Watson et al. 1963; Cowan & A'Hearn 1979; Weissman & Kieffer 1981). They have been applied in quantitative analysis of the impact of parameters on water ice activity (Marshall et al. 2019), as well as investigating mechanisms behind the disruption of a cometary nucleus (Steckloff et al. 2015). For interpreting observational data, more realistic thermophysical models are required, such as the inversion of 67P's physical properties from non-gravitational force caused by outgassing (Davidsson & Gutiérrez 2005), the evolution of morphology of the nucleus influenced by the insolation-induced activity (Keller et al. 2015a; Groussin et al. 2025), interpretation of the measured water production (Hu et al. 2017; Blum et al. 2017; Attree et al. 2019), and the explanation of distinct activity phenomena, such as shortlived outbursts from fractured terrains (Skorov et al. 2016). To further improve the model fidelity, it is necessary to consider additional complex physical mechanisms, such as the behavior of volatile transport in the nucleus given rise by its porosity, the deposition of gas on dust particles, and so on (De Sanctis et al. 1999; Davidsson & Skorov 2002; Gortsas et al. 2011; González et al. 2014; Hu & Shi 2021).

Traditionally, thermophysical models are realized by numerical treatment based on the finite-difference technique, such as the Crank-Nicolson method utilized by Hu et al. (2017). However, it usually requires relatively high computation cost due to the iterative calculations performed on a large time-space grids, leading to limited efficiency in related studies and posing challenges for practical implementation. Firstly, the inverse extraction of physical parameters automatically from observational data is a difficult problem, typically, the forward trial-and-error method is utilized (Keller et al. 2015a; Hu et al. 2017; Attree et al. 2019). In addition, Combi et al. (2019, 2021b, 2023) have been continuously updating the water production rate curves of numerous comets, although directly solving the thermal models through numerical simulation to analyze these data is clearly a highly challenging task. The orbital evolution of the comet influenced by non-gravitational acceleration from sublimation (Davidsson & Gutiérrez 2005; Attree et al. 2019; Jewitt et al. 2020), as well as the changes in its spin state due to the non-gravitational torque (Keller et al. 2015b; Attree et al. 2024), are both complex topics in the context of forward evolution and fitting the observational data as well. If the thermophysical model is embedded within a complex physical process as a sub-mechanism, such as the shape evolution (Zhao et al. 2021) and coma formation (Shi et al. 2018), conducting large-scale or high-resolution studies is also challenging.

In our previous work, we addressed similar issues in the thermophysical modeling for asteroids by introducing deep learning techniques (Zhao et al. 2024). Here, we have applied a similar concept to develop a deep learning-based approach for the comets that can efficiently solve the thermophysical process and produce the water production rate.

This paper covers: (i) how the deep learning model can be applied to model the thermal physics of comets; (ii) whether the deep learning model can serve as a unified model for simulating the water activity of the comets in large parameter spaces; and (iii) how the deep learning model can be utilized in the analysis of real comet data.

2 Thermophysical model and water activity

In this work, we adopted the so called “dust mantle model,” in which a dry dust layer is present at the surface of the nucleus over the dust-ice mixture and the water ice sublimation occurs at the boundary between the mantle and the ice front (Hu et al. 2017, and references therein).

2.1 Dust mantle model

The nucleus is regarded as a polyhedron consisting of triangular facets. For each facet, neglecting the changes in thermal properties due to the gas from ice sublimation filling the porosity of the nucleus, the 1D heat conduction equation can be expressed as follows: $$\rho C\frac{\mathrm{\partial }}{\mathrm{\partial }t}T=\kappa \frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{x}^2}T$$(1) where ρ, C, and κ correspond to the density, specific heat capacity, and thermal conductivity, respectively. They are assumed to be different constants in dry dust mantle and dust-ice mixture, denoted as ρ_d, C_d, κ_d and ρ_m, C_m, κ_m. In fact, ρ_mC_m, κ_m can be expressed in terms of ρ_d, C_d, κ_d and ρ_i, C_i, κ_i of the pure water ice (Davidsson & Skorov 2002; Davidsson & Gutiérrez 2005): $$\begin{array}{l}{\rho }_m{C}_m=(1-f){\rho }_d{C}_d+f{\rho }_i{C}_i,\\ {\kappa }_m=h((1-f){\kappa }_d+f{\kappa }_i),\end{array}$$(2) where f is the icy area fraction and h is the Hertz factor, dependent on the morphological structure of the medium and is generally regarded as a free parameter (Davidsson & Skorov 2002).

The surface boundary condition of the dust mantle derived from conservation of energy is expressed as (Kührt & Keller 1994; Lagerros 1996): $$ϵ\sigma T_{x=0}^4-{\kappa }_d\frac{\mathrm{\partial }}{\mathrm{\partial }x}{T}_{x=0}=(1-A_B)s\psi F_{\odot }$$(3) where ∊ is the emissivity, σ is the Stefan-Boltzmann constant. (1 - A_B) SψF_θ describes the radiation flux received by the facet, A_B is the Bond albedo, s indicates whether the facet is illuminated by the Sun, considering the horizon shadows and projected shadows, ψ represents the cosine function of the solar elevation angle, and F_θ is the incident solar radiation flux determined by F_θ = S_θ/r². Here, r is the heliocentric distance of the asteroid and S_θ is the solar constant.

It is necessary to introduce an additional boundary condition at the interface between the dust mantle and dust-ice mixture, due to the energy consumption caused by water ice sublimation. Assuming the thickness of the dust mantle is X, the condition can be expressed as (Kührt & Keller 1994; Hu et al. 2017): $$-{\kappa }_d\frac{\mathrm{\partial }}{\mathrm{\partial }x}{T}_{x=X^{-}}=-{\kappa }_m\frac{\mathrm{\partial }}{\mathrm{\partial }x}{T}_{x=X^{+}}+l{Z}_{(X)}$$(4) then, l is the latent heat of water ice, Z_(X) is the mass flux of sublimation determined by (Hu et al. 2017): $$Z_{(X)}=\mathrm{\Psi }f{Z}^{\mathrm{H}-\mathrm{K}}(T_{(X)})$$(5) where $$Z^{\mathrm{H}-\mathrm{K}}(T)=\alpha (T)P_V(T)\sqrt{\frac{\hat{m}}{2\pi k_{B}T}}$$(6)

This equation comes from the Hertz-Knudsen formula, describing the sublimation flux from the surface of pure solid ice into vacuum (Kossacki et al. 1999; Gundlach et al. 2011); m̄ is the mass of a water molecules and k_B is the Boltzmann constant; P_V(T) is the saturation vapor pressure; α(T) is the sublimation coefficient, representing the reduction in sublimation flux due to the reattachment of gas molecules to the ice surface after impacting with each other. They are respectively evaluated as (Panale & Salvail 1984; Mauersberger & Krankowsky 2003; Gundlach et al. 2011): $$P_V(T)=a\exp (-\frac{b}{T}),\alpha (T)=c_0+\frac{c_1}{1+\exp (c_2-c_3/T)}$$(7) with constants a, b, c₀, c₁, c₂, and c₃ taking the values of 3.23×10¹² Pa, 6134.6 K, 0.146, 0.854, 57.78, and 11 580 K.

The parameter Ψ in Eq. (5) represents the reduction factor in the sublimation flux from the bare icy surface, accounting for the permeability of the mantle to gas flow. It is given by the expression 1/(1 + pH) (Gundlach et al. 2011; Hu et al. 2017), where p is a constant and H = X/d, with d denoting the diameter of the spherical dust aggregates in the dust mantle (Skorov & Blum 2012). We take p = 0.14 and d = 1 mm in this work, similarly to the approach taken by Hu et al. (2017). Furthermore, f in Eq. (5) is used to approximate the reduction in sublimation flux from the pure ice caused by the presence of mixed dust, as the Hertz-Knudsen formula is derived for pure ice (Crifo 1997).

The last boundary condition describes the invariability of temperature at infinitely large depth, $$\frac{\mathrm{\partial }}{\mathrm{\partial }x}{T}_{x\to \mathrm{\infty }}=0$$(8)

2.2 Normalization of thermophysical model

A direct application of deep learning methods to the model described above presents several challenges: (1) the parameters have a broad dynamic range, leading to the fact that the resulting temperature distribution could span over a wide numerical range; (2) The model is influenced by the rotational velocity of comet, which is not explicitly represented in the established equations. Instead, it is implicitly embedded in the time-dependent function Sψ and the time steps used in the subsequent numerical simulations. Therefore, normalization becomes particularly crucial in this context.

We denote the rotational angular velocity as ω and introduce the following normalized variables (Lagerros 1996), $${\overline{x}}_d=\frac{x}{l_{sd}},{\overline{x}}_m=\frac{x}{l_{sm}},\overline{t}=\omega t,\overline{T}=\frac{T}{T_e},$$(9) where l_sd and l_sm are the skin depth of the dust mantle and the dust-ice mixture, T_e is the theoretical maximum surface temperature at the heliocentric distance, r, called the characteristic temperature, given by: $$l_{sd}=\sqrt{\frac{{\kappa }_d}{{\rho }_d{C}_d\omega }},l_{sm}=\sqrt{\frac{{\kappa }_m}{{\rho }_m{C}_m\omega }},T_e=\sqrt[4]{\frac{(1-A_B)F_{\odot }}{ϵ\sigma }},$$(10) the heat conduction equation Eq. (1) and three boundary conditions Eq. (3), Eq. (4), and Eq. (8) can then be normalized as: $$\begin{array}{l}\frac{\partial \overline{T}}{\partial \overline{t}}=\frac{{\partial }^2\overline{T}}{\partial {\overline{x}}^2},\\ {\overline{T}}_{{\overline{x}}_d=0}^4-{\text{Φ}}_d\frac{\partial }{\partial {\overline{x}}_d}{\overline{T}}_{{\overline{x}}_d=0}=E(\overline{t}),\\ -{\text{Φ}}_d\frac{\partial }{\partial {\overline{x}}_d}{\overline{T}}_{{\overline{x}}_d={\overline{X}}^{-}}=-{\text{Φ}}_m\frac{\partial }{\partial {\overline{x}}_m}{\overline{T}}_{{\overline{x}}_m={\overline{X}}^{+}}+l{\overline{Z}}_{(\overline{X})},\\ \frac{\partial }{\partial {\overline{x}}_m}{\overline{T}}_{{\overline{x}}_m\to \infty }=0,\end{array}$$(11) where Φ_d and Φ_m are the thermal parameters (Spencer et al. 1989) of the dust mantle and dust-ice mixture, defined by $${\mathrm{\Phi }}_d=\frac{\sqrt{{\rho }_d{C}_d\omega {\kappa }_d}}{ϵ\sigma T_e^3},{\mathrm{\Phi }}_m=\frac{\sqrt{{\rho }_m{C}_m\omega {\kappa }_m}}{ϵ\sigma T_e^3}.$$(12)

Here, E(t̄) is the radiation flux rescaled in the range [0, 2π] for variables t̄ and [0, 1] for E. Then, Z̄_(X̄) = Z(X)/(1 - A_B)F_Θ is the sublimation flux after variation.

It is evident that the radiation flux and the derived temperature can be successfully rescaled in the range of [0, 1], with all the physical parameters explicitly represented in the mathematical expression of the equations. This transformation can significantly enhance the efficiency and ease of dataset creation for deep learning. For instance, when studying different comets, the original model for the radiation flux functions leads to varying domains of definition due to differences in rotational velocities. Additionally, the range of values fluctuates with variable parameters such as the Bond albedo, A_B, and heliocentric distance, r. This variability complicates the construction of a comprehensive set of radiation flux functions. However, the normalized model can effectively address this issue through the parameter extraction.

Excluding the radiation flux, we have identified and summarized all the independent parameters that influence the equations: $$\begin{array}{l}{p}_1=(1-A_B)F_{\odot },p_2=ϵ,p_3=X,\\ p_4={\rho }_d{C}_d,p_5=\omega ,\\ p_6={\kappa }_d,p_7={\kappa }_m,p_8=f.\end{array}$$(13)

In the following discussion, we treat all these parameters as variable within a physically reasonable range in order to obtain a generalized model.

2.3 Numerical treatment

The radiation flux needs to be calculated before solving the equations. It is a periodic function about time, determined by the angle between the normal vector of the surface element and the solar direction, incorporating both the horizon and projected shadows. The calculation method is consistent with our previous work for asteroids (secondary radiation and self-heating are not considered in this work), as described in Sect. 2.1.2 of Zhao et al. (2024). It is worth noting that the radiation flux used to train the network should be randomly generated, which will be detailed in Sect. 3, while the illumination and shadow mentioned here are solely specific to real-world situations.

The backward finite difference method is utilized to solve the partial differential equations, considering its computational efficiency. First, we divided the continuous variables into a set of space-time grids, where the i-th grid of depth is denoted by iΔx̄ and the j-th grid of time is denoted as jΔt̄. As a result, the heat conduction equation is discretized into $${\overline{T}}_{i,j+1}=(1-2\frac{\mathrm{\Delta }\overline{t}}{\mathrm{\Delta }{\overline{x}}^2}){\overline{T}}_{i,j}+\frac{\mathrm{\Delta }\overline{t}}{\mathrm{\Delta }{\overline{x}}^2}({\overline{T}}_{i+1,j}+{\overline{T}}_{i-1,j}).$$(14)

The thickness of the dust mantle is much smaller than that of the dust-ice mixture. To improve computational efficiency while ensuring accuracy, different spatial step sizes are applied in two layers, with a smaller step size Δx_d = 1 mm for the dust mantle and a larger one Δx_m = 1 cm for the dust-ice mixture. Then, Δx̄_d and Δx̄_m can be obtained by dividing Δx_d and Δx_m by the skin depth, l_sd, l_sm. The time step is taken as Δt̄ = 2π/1200. We denote the temperature at the interface between the two layers as T_m,j, then the three boundary conditions are expressed in the finite difference scheme as follows: $$\begin{array}{l}{\overline{T}}_{0,j}^4-{\text{Φ}}_d\frac{{\overline{T}}_{1,j}-{\overline{T}}_{0,j}}{\text{Δ}{\overline{x}}_d}=E(j\text{Δ}\overline{t}),\\ -{\text{Φ}}_d\frac{{\overline{T}}_{m,j}-{\overline{T}}_{m-1,j}}{\text{Δ}{\overline{x}}_d}+{\text{Φ}}_m\frac{{\overline{T}}_{m+1,j}-{\overline{T}}_{m,j}}{\text{Δ}{\overline{x}}_m}=l{\overline{Z}}_{(\overline{X})},\\ {\overline{T}}_{n,j}={\overline{T}}_{n-1,j}.\end{array}$$(15)

The surface temperature, T_{0, j}, and the interface temperature, T_m,j, can be determined from the known T_1,j, T_m-1, j, T_m+1,j by taking advantage of Newton-Raphson iterative algorithm to solve the boundary conditions.

The overall temperature needed to exhibit a stable and static distribution due to the rotation of the comet. After one complete spin of the nucleus, the temperature should be invariant at the subsolar point, which holds at all depths. This static condition is used as a convergence criterion for the iteration. In this work, the criterion is for the mean absolute error of the normalized temperature between iterations to be below 1 × 10⁻⁴.

2.4 Water activity

The water activity of the comet is effectively represented by the sublimation flux in the model. The global total outgassing flux of water from the comet nucleus is the integral of the sublimation flux over the entire surface. In the polyhedral shape model, this can be expressed as follows (Hu et al. 2017): $$Z_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}=\sum _k{Z}_k{S}_k$$(16) where Z_k indicates the sublimation flux of k-th facet and its area is denoted as S_k. Furthermore, we calculated the global water production rate by averaging the global outgassing flux of water over one rotational period, which quantifies the number of water molecules produced per second, given by (Hu et al. 2017) where P is the rotational period.

3 ThermoONet

In Zhao et al. (2024), we showed that many factors in a thermophysical model makes the equation-driven deep learning models such as physics-informed neural networks (PINNs; Raissi et al. 2019) unsuitable, despite their successful application in solving ordinary differential equations (ODEs) and partial differential equations (PDEs; Cai et al. 2020; Zobeiry & Humfeld 2021; He et al. 2021; Martin & Schaub 2022a,b; Mathews & Thompson 2023; Laghi et al. 2023). Therefore, for the asteroid thermophysical modeling, we adopted deep operator neural network (DeepONet) with strong generalization abilities (Lu et al. 2021; Garg et al. 2022; He et al. 2023; Branca & Pallottini 2024). Such a network is well suited for handling problems with variable boundary conditions. Following the same consideration, we continue to utilize DeepONet for the development of ThermoONet in this work.

The surface of an asteroid is treated as a layer of dust, whose temperature is influenced primarily by the radiation flux and the thermal parameter (a total of two independent parameters). However, for comets, the distinct physical properties of the dust-ice mixture layer in contrast to the dust mantle, along with the phase transitions of water ice, necessitate a more complex selection of parameters to describe the thermophysical process. Consequently, we identified a total of eight independent parameters (Eq. (13)). In this work, we modified the structure of the neural network to address the high dimensionality of the parameters. Furthermore, we manually optimized the distribution of the training dataset, based on the relative impact of different parameters on the results. The specific methods are introduced in detail below.

3.1 Architecture of the neural network

The principles and structure of DeepONet are explained in detail in Lu et al. (2021). Here, we just provide a brief summary. The core idea is to utilize a neural network to map the different variables from the partial differential equations into a common mathematical space for processing. The network that maps the parameters or functions of the equations is referred to as the branch network, while the network that maps the independent variables is called the trunk network. If we denote the parameters or functions as independent variables, such as as u, y, these two mapping relationships can be expressed as u → G(u) and y → G(u)(y), where G is called the operator.

In the thermophysical modeling of cometary systems, the operator mapping process involves the radiation flux function, E(t̄), and eight characteristic parameters from Eq. (13), with the depth and time serving as independent variables. A DeepONet architecture could theoretically be constructed with branch networks encoding radiation flux function and parameters, such as trunk networks handling independent variables. The large number of parameters renders the training of an accurate neural network particularly challenging. This necessitates a simplification of the network architecture through problem-specific adaptations. Theoretically the sublimation flux fundamentally depends on the temperature at the dust mantle-and-dust-ice mixture interface. This physical understanding enables significant architectural simplification by reducing the network output to this subsurface temperature. Consequently, the trunk network becomes redundant and can be eliminated, with the branch network outputs directly providing the temperature predictions. For the comprehensive modeling of non-gravitational effects induced by sublimation processes, additional consideration must be given to the surface temperature of the dust mantle, which should be incorporated as another output. This dual-temperature approach enables the simultaneous prediction of both the subsurface temperature governing sublimation dynamics and the surface temperature, which is essential for calculating the velocity of the gas (detailed analyses of these effects are provided in Attree et al. 2019 and Jewitt et al. 2020).

In the asteroid network, we model the thermal behavior of the asteroid using two separate branch networks: one for encoding the radiation flux function and the other for encoding the thermal parameter (Zhao et al. 2024). The treatment of generating random radiation flux functions can be directly adapted to the comet model (see Sect. 3.2 for more about the dataset). However, the comet model introduces a larger number of independent parameters exhibiting heterogeneous sensitivity to the solutions of the model, so direct adoption of the parameter encoding scheme for asteroids proves ineffective. A single branch network collectively encoding all parameters can reduce sensitivity to each individual parameter, especially for those parameters significantly influencing the results. Conversely, constructing separate branch networks for individual parameters would result in excessive computational memory allocation and time cost. To resolve this challenge, we propose a parameter grouping method based on sensitivity analysis. Parameters governing the thermal processes are partitioned into two distinct categories: high-sensitivity and moderate-sensitivity parameters. Each category is encoded through dedicated branch network, achieving an optimal balance between parameter sensitivity and computational efficiency.

The schematic architecture of the neural network adapted for this work is shown in Fig. 1, where the branch networks are presented in Fig. 2. Parameter group₁ consists of the high-sensitivity parameters, p₁, which exhibit a dominant influence over system behavior (more details are shown in Appendix A), while parameter group₂ consists of the remaining parameters. The neural network encodes these grouped parameters through dedicated branch networks (branch net₁ and branch net₂), followed by a Hadamard product operation between their respective outputs. This combined tensor undergoes transformation through an attention mechanism layer, designed to amplify the features through adaptive weight allocation (Vaswani et al. 2017; Woo et al. 2018; Wang et al. 2020). Notably, each branch network incorporates internal attention mechanisms to perform hierarchical feature enhancement. The target subsurface temperature is derived from a averaging operation to the Hadamard product between the amplified features and branch net₃ output.

Fig. 1 Architecture of ThermoONet. The network consists of three branch networks with distinct inputs, parameter group, p1, p2-8, and the radiation flux function, E, where E is presented as some discrete locations with m elements sampled on [t̄1,t̄2,...,t̄m]. A layer incorporating the attention mechanism processes the Hadamard product of the outputs from branch net1 and branch net2. The resulting output from this layer [$[b_1^1,b_2^1,...,b_p^1$] is then element-wise multiplied (Hadamard product), with the output from branch net3 [$[b_1^2,b_2^2,...,b_p^2$]. The average of the resulting values is taken as the final output.

Fig. 2 Detailed architectures of branch netn and attention mechanism, where FC is the full connected layer. The input of branch netn experiences the initial feature extraction through a FC layer and then carries out the feature amplification by attention mechanism, where i represents the layer number. 3, 6, and 6 layers are respectively used in branch networks from net1 to net3.

3.2 Dataset and training

The neural network developed in this work focuses on the diurnal water ice activity. Seasonal variations in thermal conditions are reflected through variations in insolation.

As mentioned above, the neural network will be trained with a nine-dimensional parameter space, comprising the scalar parameters, p₁-p₈, and one functional input, E, representing the normalized diurnal insolation curve. To ensure that ThermoONet is able to efficiently generalize across different cometary environments, we should design a training dataset comprised of a broad, diverse selection of diurnal radiation flux curves.

For training the asteroid network, we produced a dataset with random pairs of the radiation flux function, E, and the thermal parameter, Φd (Zhao et al. 2024). The corresponding temperatures were calculated by the numerical simulation. Random functions are generated by the mean-zero Gaussian random field (GRF; Lu et al. 2021), then rescaled in the range [0, 1] through the sigmoid function 1/(1 + e^−(x+d)), where d represents the translation distance along the y axis in order to ensure the completeness of the function space (Zhao et al. 2024). For the network used in this study, the training process employs a modified stochastic sampling framework, adapted for cometary thermophysical model requirements. The only difference is that all eight parameters undergo a simultaneous joint randomization to maintain proper parameter covariance structures. A special treatment can be implemented for high-sensitivity parameter p₁ through a stratified sampling strategy; p₁ is randomly generated within different value ranges, with a similar number in each range. This balanced multi-range generation mechanism, grounded in prior knowledge of the nonlinear dominance of p₁, ensures that the network effectively captures the substantial influence of various ranges on the resulting temperature. In the practical training process, we generated a dataset with 36 000 sets of data.

The loss function is the mean square error (MSE) between the temperature calculated with numerical simulation and the output from the neural network and the adaptive moment estimation method (Adam) is used as the optimization algorithm. Just as in the training process of the asteroid network in Zhao et al. (2024), the learning rate decay approach was utilized to mitigate issues related to wide oscillations of the loss function and slow convergence.

We called the trained neural network ThermoONet, designed on the basis of the core principles of DeepONet and specifically tailored for cometary thermophysical modeling.

4 Results

The cost of computation time and the level of accuracy are used to assess the performance of ThermoONet.

4.1 Computational cost

We compared the computational time consumed by three approaches: ThermoONet with GPU, ThermoONet with CPU, and a traditional numerical simulation, for calculating the global ice front temperature distribution of a spherical nucleus with varying number of facets. It is worth mentioning that dataset construction required approximately 6 hours, while the neural network training process was completed within 20 minutes. Once trained, the network is ready for immediate application and does not require further adjustment, so we only need to analyze the inference time of the network. The ThermoONet calculations were carried out on both GPU and CPU for over 20 times to get an estimate on the uncertainty. The time consumption of the numerical simulation is evaluated by multiplying the average time cost per facet (calculated based on the temperature of a large number of facets) by the total number of facets. All of the computations were performed on a personal computer¹ and under the environment of Python 3.8. Our results indicate, as shown in Fig. 3, that:

• Calculations using ThermoONet are approximately six orders of magnitude faster than the numerical simulation.

• Computational time costs by numerical and CPU-based ThermoONet increase almost linearly with number of facets, while the computation cost by GPU-based ThermoONet stays almost constant against an increasing number of facets. When we need to calculate a large number of global temperature distributions, it is worth noting that the vectorized operations can be used next to ensure that the computation time does not increase linearly with the number of facets.

Fig. 3 Comparison of computational cost between ThermoONet and traditional numerical simulation for the shape models with varying number of facets. The error bars represent the maximum and minimum time costs in the repeated calculation for more than 20 times.

4.2 Accuracy assessment

First, we used a spherical comet to assess the accuracy of temperature predictions with varying obliquity, β, (the angle between the spin axis and the normal vector of the orbital plane) and heliocentric distance, r. We also investigated how the errors change with different key model parameters such as the thickness of the dust mantle, X, and icy area fraction, f , where fixed parameters are listed in Table 1. We measure the mean absolute percentage error (MAPE) of temperature between the neural network and numerical simulation, denoted as MAPE (T).

Figure 4 shows the results. The overall error in temperature fluctuates as the radiation flux function or parameters change. With the inclusion of more parameters, the error inevitably increases compared to that of the deep learning model for asteroids (Zhao et al. 2024). However, it is still maintained at a reasonably low level, with a maximum below 4% and an average around 2%. We notice that the error does not show significant correlation with varying heliocentric distances, which confirms the effectiveness of the specially designed training dataset.

As shown in Eqs. (5)-(7), the sublimation flux of water ice exhibits an exponential relationship with the temperature, meaning that any error in temperature has an even more significant impact on the calculated sublimation flux. Moreover, given that many cometary nucleus show a concave shape, the projected shadows between facets complicate the radiation flux function. This presents a great challenge for the generalization of the neural network, which is even more severe than in the case of asteroid (Zhao et al. 2024). To validate the generalization capability of ThermoONet, it is applied to a real comet with irregular shape. In practice, we take comet 67P as the example and simulate its temperature with a shape model of 2868 facets (Jorda et al. 2016). Ephemerides of 67P is retrieved using Rosetta SPICE kernels (Acton 1996; ESA SPICE Service 2017). We evaluated both the temperature and water production rate at different orbital phases by comparing them with the results from numerical solutions (see Fig. 5). We also examined the temperature distributions at three specific orbital locations before, near, and after the perihelion (Fig. 6). The other parameters are the same as in Table 1.

The average temperature errors are about 2.5% and the global water production rates from two methods exhibit the same trend. The average relative errors are respectively 5.04% for (X = 5mm, f = 1%), 19.25% for (X = 5mm, f = 10%), and 10.26% for (X = 10 mm, f = 1%). The patterns of the error distribution shown in Fig. 5 suggest that the temperature error is the smallest in the high-temperature region, typically below 3 K. The primary source of error lies in the mid-temperature regions. However, only a few facets exhibit errors exceeding 10 K.

5 Application to the interpretation of water production rate

Interpreting the water production rate as a function of the heliocentric distance using thermophysical models has been an intriguing topic in cometary science as it helps us understand the nature of cometary activity (Keller et al. 2015a; Hu et al. 2017; Attree et al. 2019; Marshall et al. 2019; Skorov et al. 2020). It is a time-consuming process to calculate the water production rate curve (hereafter, the 'water curve') using traditional numerical simulation, which requires computations of the global temperature distribution at different orbital phases. Thus, it needs to be recalculated whenever the parameters are adjusted. Due to the large number of physical parameters involved, previous studies often fix certain parameters to some typical values, leaving only afew that are adjustable. When fitting the observational data, itis difficult to take advantage of the automated algorithms to extract parameters from the observational data; thus, the parameters are typically determined through empirical trial-and-error testing.

Generating water curves with ThermoONet is efficient, thereby enabling the application of global optimization algorithms, such as genetic algorithm (GA), particle swarm optimization (PSO), simulated annealing (SA), and so on, to find the best-fitting parameters. To do so, we implemented ThermoONet to calculate the water curve in three steps:

• A lookup table was generated by calculating illumination and shadows under varying solar positions, incorporating the shape and rotational state of the comet. This is achieved through Möller-Trumbore algorithms and interpolation algorithms.

• Eight scalar parameters are computed based on the physical properties of the comet and orbital ephemerides, the radiation flux for each facet then is retrieved from the precomputed lookup table.

• The derived nine-dimensional parameters serve as input to ThermoONet, which outputs the sub-surface temperature at the ice front of each facet. The results subsequently drive the calculation of water ice sublimation rates and associated activity levels.

Fig. 4 Differences in the output subsurface temperature between ThermoONet and numerical simulation as functions of obliquity, β, and heliocentric distance, r. The left panels present the error functions for different thicknesses of dust mantle, with the blue, black, and red lines, representing 5 mm, 10 mm, and 15 mm respectively. The right panels illustrate the error functions for varying icy area fractions, with the blue, black, and red lines respectively corresponding to 1%, 10%, and 20%. The panels are fixed to r = 2 au for the top and β = 0° for the bottom, along with f = 10% for the left panels and X = 10 mm for the right panels. For the settings of other parameters, we refer to Table 1.

Fig. 5 Differences in global water production rate derived from Ther-moONet and the numerical simulation at different orbital phases of comet 67P. The top panel is the temperature calculation errors, where the black line was derived with X = 5 mm, f = 1%, and the red line with X = 10 mm, f = 1%, the blue line with X = 5 mm, and f = 10%. The asterisks denote specific examples of temperature distributions shown in Fig. 6. The bottom panel is the water production rate, where the real lines stand for the results from ThermoONet, while the circles refer to the numerical solutions, with the same parameters as in the top one.

5.1 Fitting the water curve of 67P

We combine ThermoONet with SA to fit the observed water curve of comet 67P from Läuter et al. (2020), the shape model with 23 788 facets is used (Jorda et al. 2016). We chose five parameters to be fitted, including the thickness of the dust mantle, X, the icy area fraction, f , the product between the density of the dust mantle, ρ_d, and the specific heat capacity, C_d, the heat conductivity of the dust mantle, κ_d, and the dust-ice mixture, κ_m. The initial values of the five parameters are selected as in Table 2 and the corresponding water curve is generated by applying ThermoONet at every observation time. The loss function is MAPE between the predicted results and the average observed data and SA is applied to search for the preliminary optimal parameter combination that minimizes MAPE. Finally, a batch gradient descent (BGD) was utilized to refine the parameters, converging to a solution that further reduces the loss function. The best-fitting water curve is shown in Fig. 7 and the optimized parameter combination is shown in Table 2.

The overall fitting performance is satisfactory with MAPE = 37.08%, capturing the asymmetrical trend before and after perihelion. A slight underestimation is observed near the perihelion, which is consistent with the finding in Hu et al. (2017). This phenomenon may be explained by several factors, including the varying levels of activity across different regions of 67P (Attree et al. 2019), the gradual exposure of water ice as the comet approaches perihelion (Fornasier et al. 2016; Filacchione et al. 2020), and the temperature-dependent changes in thermal conductivity (Skorov et al. 2020). It is worth mentioning that the advantage of ThermoONet in computation efficiency could be of help in assessing aforementioned mechanisms.

Fig. 6 Sub-surface temperature distribution of comet 67P modeled by ThermoONet and the numerical simulation at three specific positions marked by asterisks in the top panel of Fig. 5. The top row shows the normalized radiation flux and the second to fourth rows show the results calculated with ThermoONet, the numerical simulation, and the differences between them.

5.2 Fitting the water curve of 21P

There is a growing database of water curves of different types of comets retrieved using observations by solar observatories Combi et al. (2019, 2021b, 2023). Most of these comets still lack well-defined parameters such as size, shape, rotation characteristics, and thermophysical parameters. It could prove difficult for traditional numerical simulations to fit these data by searching a large parameter space. However, by combining trial-and-error testing with ThermoONet, exploring a large parameter space becomes feasible. We take comet 21P as an example to demonstrate the application of the trial-and-error testing. The water curve of 21P is significantly different from that of 67P, remaining essentially flat before perihelion and exhibiting an obvious decline after perihelion (see more in Combi et al. 2021a). A preliminary consideration is that the solar illumination area of 21P gradually decreases during the observation period, which could be the combined effect of its shape and the orientation of its spin axis. This was described in Marshall et al. (2019). We can iteratively adjust the parameters based on trial-and-error results by taking advantage of ThermoONet.

Figure 8 shows preliminary fitting results. The fixed parameters are ∊ = 1, ω = 1.74 × 10⁻⁴ rad s⁻¹ (the rotation period is uncertain with a wide range from 9.5 to 19 h (Leibowitz & Brosch 1986; Moulane et al. 2020). The parameters included in the fitting are shown in Table 3. The test shapes are respectively sphere with a_x = a_y = a_z = 3.5 km, biaxial ellipsoid with a_x = 2a_y = 2a_z = 5.2 km, and triaxial ellipsoid a_x = 2a_y = 4a_z = 5.2 km. As shown in Fig. 8, the evolutionary trend of the triaxial ellipsoid better matches the observational data, confirming our previous hypothesis regarding the solar illumination area, meaning that the farther the comet is before the perihelion, the larger the region with the higher temperature.

Fig. 7 Fitting the global water production rate of 67P as a function of time. The blue rectangles represent the observation data from Läuter et al. (2020), and the red line corresponds to the best fit curve achieved by combining ThermoONet with the optimization algorithm. The black dashed line indicates the epoch of the perihelion. A detailed description of the observation data and sources thereof can be found in Läuter et al. (2020).

6 Conclusion

In this work, we have developed ThermoONet, a DeepONet-based neural network for modeling cometary activity. Due to the high-dimensional parameter space in cometary thermophysical model, we tailored the architecture of DeepONet by removing the trunk network and trained it specifically for the target temperature. With the solar flux as input, ThermoONet predicts the subsurface temperatures and water production rate. By taking advantage of the generalization ability of DeepONet, ThermoONet could be used in modeling comets with varying physical parameters and shapes.

We conducted a comprehensive evaluation of the network's performance, including its computation time, accuracy, and generalization ability. The results show that ThermoONet deployed on GPU has the highest efficiency and generates the global sublimation flux approximately six orders of magnitude faster than using numerical simulation. The accuracy of ThermoONet varies with parameters, with an average of around 2%. Ther-moONet can be successfully applied to irregular comets as exemplified by comet 67P, with the water production rate curve closely matching the numerical simulation, demonstrating a strong generalization ability and practicability.

We tested the capability of ThermoONet in computationally demanding tasks by fitting global water production rate curves of comets using global optimization algorithms. By applying a simulated annealing, we successfully fit the observed water production rate curve of comet 67P derived with the Rosetta mission data and obtained a set of optimized physical parameters of 67P. Additionally, we fit the water curve of comet 21P to demonstrate the capability of combining ThermoONet with trial-and-error testing. Although limited knowledge is available for 21P, we managed to reconstruct the unique trend observed in its water production rate by adjusting the shape, rotation, and thermal parameters.

We anticipate a broad range of application for ThermoONet in various research topics concerning small body activities. In addition to its application to large-scale and high-resolution data analysis, it can also serve as an efficient model that can be embedded into more complex physical systems.

Fig. 8 Fitting the global water production rate of 21P as a function of days from perihelion by applying ThermoONet to different shapes. The black line corresponds to a sphere model, the blue line is a ellipsoid model with ax = 2ay = 2az and the red line is a ellipsoid with ax = 2ay = 4az. Other parameters are as stated in Table 3. The scatter points with different colors and shapes represent observational data at different apparitions (Combi et al. 2019), where purple circles from 1998, blue triangles from 2005, green squares from 2012, and black asterisks from 2018. The detailed description about the observation data and sources thereof can be found in Combi et al. (2019).

Acknowledgements

The authors thank Prof. Groussin for the insightful review comments that helped to significantly improve the manuscript. The authors thank Dr. Xuanyu Hu for extensive discussions on the thermophysical modeling for comets. XS thanks members of the ISSI International Team “Timing and Processes of Planetesimal Formation and Evolution” for inspirational discussions. This work is financially supported by the National Natural Science Foundation of China (No. 12233003 and 12073011) and the China Manned Space Program with grant no. CMS-CSST-2025-A16.

Appendix A Parameter sensitivity

We conducted the parameter sensitivity analysis using the controlled variable approach with a fixed radiation flux as an example. When modifying a single parameter, the remaining parameters maintain their reference values, as specified in Table A.1. Here, the rotational angular velocity p₅ is represented as the rotation period. Figure A.1 shows the resultant variations in temperature curves induced by alterations of different parameters; p₁ primarily influences the value range of temperature distributions, whereas the remaining parameters affect the shape of temperature curves within this established range.

Fig. A.1 Comparative analysis of parameter sensitivity. The first panel depicts the selected example radiation flux, followed by eight panels illustrating different temperature curves under varying values of the eight distinct parameters. The titles of these panels show the variation ranges with units as the same as Table A.1 and the curves with darker hues correspond to larger values of parameters.

References

All Tables

All Figures

Fig. 1 Architecture of ThermoONet. The network consists of three branch networks with distinct inputs, parameter group, p1, p2-8, and the radiation flux function, E, where E is presented as some discrete locations with m elements sampled on [t̄1,t̄2,...,t̄m]. A layer incorporating the attention mechanism processes the Hadamard product of the outputs from branch net1 and branch net2. The resulting output from this layer [$[b_1^1,b_2^1,...,b_p^1$] is then element-wise multiplied (Hadamard product), with the output from branch net3 [$[b_1^2,b_2^2,...,b_p^2$]. The average of the resulting values is taken as the final output.

In the text

	Fig. 2 Detailed architectures of branch netn and attention mechanism, where FC is the full connected layer. The input of branch netn experiences the initial feature extraction through a FC layer and then carries out the feature amplification by attention mechanism, where i represents the layer number. 3, 6, and 6 layers are respectively used in branch networks from net1 to net3.
In the text

	Fig. 3 Comparison of computational cost between ThermoONet and traditional numerical simulation for the shape models with varying number of facets. The error bars represent the maximum and minimum time costs in the repeated calculation for more than 20 times.
In the text

Fig. 4 Differences in the output subsurface temperature between ThermoONet and numerical simulation as functions of obliquity, β, and heliocentric distance, r. The left panels present the error functions for different thicknesses of dust mantle, with the blue, black, and red lines, representing 5 mm, 10 mm, and 15 mm respectively. The right panels illustrate the error functions for varying icy area fractions, with the blue, black, and red lines respectively corresponding to 1%, 10%, and 20%. The panels are fixed to r = 2 au for the top and β = 0° for the bottom, along with f = 10% for the left panels and X = 10 mm for the right panels. For the settings of other parameters, we refer to Table 1.

In the text

Fig. 5 Differences in global water production rate derived from Ther-moONet and the numerical simulation at different orbital phases of comet 67P. The top panel is the temperature calculation errors, where the black line was derived with X = 5 mm, f = 1%, and the red line with X = 10 mm, f = 1%, the blue line with X = 5 mm, and f = 10%. The asterisks denote specific examples of temperature distributions shown in Fig. 6. The bottom panel is the water production rate, where the real lines stand for the results from ThermoONet, while the circles refer to the numerical solutions, with the same parameters as in the top one.

In the text

	Fig. 6 Sub-surface temperature distribution of comet 67P modeled by ThermoONet and the numerical simulation at three specific positions marked by asterisks in the top panel of Fig. 5. The top row shows the normalized radiation flux and the second to fourth rows show the results calculated with ThermoONet, the numerical simulation, and the differences between them.
In the text

Fig. 7 Fitting the global water production rate of 67P as a function of time. The blue rectangles represent the observation data from Läuter et al. (2020), and the red line corresponds to the best fit curve achieved by combining ThermoONet with the optimization algorithm. The black dashed line indicates the epoch of the perihelion. A detailed description of the observation data and sources thereof can be found in Läuter et al. (2020).

In the text

Fig. 8 Fitting the global water production rate of 21P as a function of days from perihelion by applying ThermoONet to different shapes. The black line corresponds to a sphere model, the blue line is a ellipsoid model with ax = 2ay = 2az and the red line is a ellipsoid with ax = 2ay = 4az. Other parameters are as stated in Table 3. The scatter points with different colors and shapes represent observational data at different apparitions (Combi et al. 2019), where purple circles from 1998, blue triangles from 2005, green squares from 2012, and black asterisks from 2018. The detailed description about the observation data and sources thereof can be found in Combi et al. (2019).

In the text

Fig. A.1 Comparative analysis of parameter sensitivity. The first panel depicts the selected example radiation flux, followed by eight panels illustrating different temperature curves under varying values of the eight distinct parameters. The titles of these panels show the variation ranges with units as the same as Table A.1 and the curves with darker hues correspond to larger values of parameters.

In the text