# Modelling methanol and hydride formation in the JWST Ice Age era

Izaskun Jiménez-Serra, Andrés Megías, Joseph Salaris, Herma Cuppen, Angèle Taillard, Miwha Jin, Valentine Wakelam, Anton Vasyunin, Paola Caselli, Yvonne Pendleton, et al.

► To cite this version:

Izaskun Jiménez-Serra, Andrés Megías, Joseph Salaris, Herma Cuppen, Angèle Taillard, et al. Modelling methanol and hydride formation in the JWST Ice Age era. Astronomy & Astrophysics - A&A, 2025, 695, pp.A247. 10.1051/0004-6361/202452389. hal-05345499

HAL Id: hal-05345499

https://hal.science/hal-05345499v1

Submitted on 6 Nov 2025

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
# Modelling methanol and hydride formation in the JWST Ice Age era

Izaskun Jiménez-Serra¹,*, Andrés Megías¹, Joseph Salaris², Herma Cuppen², Angèle Taillard¹,
Miwha Jin³,⁴, Valentine Wakelam⁵, Anton I. Vasyunin⁶, Paola Caselli⁷, Yvonne J. Pendleton⁸,
Emmanuel Dartois⁹, Jennifer A. Noble¹⁰, Serena Viti¹¹,¹²,¹³, Katerina Borshcheva⁶,¹⁴, Robin T. Garrod¹⁵,
Thanja Lamberts¹⁶,¹¹, Helen Fraser¹⁷, Gary Melnick¹⁸, Melissa McClure¹¹, Will Rocha¹¹,
Maria N. Drozdovskaya¹⁹, and Dariusz C. Lis²⁰

*(Affiliations can be found after the references)*

Received 27 September 2024 / Accepted 6 February 2025

## ABSTRACT

**Context.** Recent JWST observations have measured the ice chemical composition towards two highly extinguished background stars, NIR38 and J110621, in the Chamaeleon I molecular cloud. The observed excess of extinction on the long-wavelength side of the H₂O ice band at 3 µm has been attributed to a mixture of CH₃OH with ammonia hydrates (NH₃·H₂O), which suggests that CH₃OH ice in this cloud could have formed in a water-rich environment with little CO depletion. Laboratory experiments and quantum chemical calculations suggest that CH₃OH could form via the grain surface reactions CH₃ + OH and/or C + H₂O in water-rich ices. However, no dedicated chemical modelling has been carried out thus far to test their efficiency. In addition, it remains unexplored how the efficiencies of the proposed mechanisms depend on the astrochemical code employed.

**Aims.** We modelled the ice chemistry in the Chamaeleon I cloud to establish the dominant formation processes of CH₃OH, CO, CO₂, and of the hydrides CH₄ and NH₃ (in addition to H₂O). By using a set of state-of-the-art astrochemical codes (MAGICKAL, MONACO, Nautilus, UCLCHEM, and KMC simulations), we can test the effects of the different code architectures (rate equation vs. stochastic codes) and of the assumed ice chemistry (diffusive vs. non-diffusive).

**Methods.** We consider a grid of models with different gas densities, dust temperatures, visual extinctions, and cloud-collapse length scales. In addition to the successive hydrogenation of CO, the codes' chemical networks have been augmented to include the alternative processes for CH₃OH ice formation in water-rich environments (i.e. the reactions CH₃ + OH → CH₃OH and C + H₂O → H₂CO).

**Results.** Our models show that the JWST ice observations are better reproduced for gas densities ≥10⁵ cm⁻³ and collapse timescales ≥10⁵ yr. CH₃OH ice formation occurs predominantly (>99%) via CO hydrogenation. The contribution of reactions CH₃ + OH and C + H₂O is negligible. The CO₂ ice may form either via CO + OH or CO + O depending on the code. However, KMC simulations reveal that both mechanisms are efficient despite the low rate of the CO + O surface reaction. CH₄ is largely underproduced for all codes except for UCLCHEM, for which a higher amount of atomic C is available during the translucent cloud phase of the models. Large differences in the predicted abundances are found at very low dust temperatures (T_dust<12 K) between diffusive and non-diffusive chemistry codes. This is due to the fact that non-diffusive chemistry takes over diffusive chemistry at such low T_dust. This could explain the rather constant ice chemical composition found in Chamaeleon I and other dense cores despite the different visual extinctions probed.

**Key words.** ISM: clouds – dust, extinction – ISM: molecules

## 1. Introduction

Interstellar ices represent the main reservoir of organic and volatile components in molecular clouds and star-forming regions (see e.g. the review by Boogert et al. 2015, and references therein). Ices grow during the collapse of molecular clouds, largely influencing the chemical composition of protoplanetary disks and, ultimately, of planetesimals and cometesimals (see e.g. Quénard et al. 2018a; Coutens et al. 2020). A fraction of the chemical compounds measured in comet 67P/Churyumov-Gerasimenko by the ESA Rosetta mission could be of pre-stellar origin, which suggests that they could have formed during the dense molecular cloud phase (Altwegg et al. 2016; Drozdovskaya et al. 2019, 2021).

The composition of ices in molecular clouds and star-forming regions is studied through absorption spectroscopy of

interstellar ices in the near- and mid-IR against bright background sources (see e.g. Whittet et al. 1983; Pendleton et al. 1999; Boogert et al. 2013; Goto et al. 2018). However, these observations have been scarce because they are challenging to perform from the ground, due to telluric features and high sky-background noise in the IR. In addition, dense molecular clouds and dense cloud cores can have very high levels of visual extinction (A_V ≥ 20 mag), for which detailed ice absorption spectroscopic observations in the near- and mid-IR require IR-bright background field stars and very high sensitivity.

From the ground, H₂O, CO, and methanol (CH₃OH) ices have been detected towards a small sample of dense molecular clouds such as Taurus, Lupus, Serpens, and the Pipe Nebula (Whittet et al. 1983; Boogert et al. 2013; Goto et al. 2018; Perotti et al. 2020) and towards dense cloud cores such as the star-forming core L483, and the L1544 and L694-2 pre-stellar cores (Chu et al. 2020; Goto et al. 2021). From space, the Infrared

<footnote>* Corresponding author; ijimenez@cab.inta-csic.es</footnote>
Space Observatory (ISO), the Spitzer Space Telescope, and the AKARI Infrared Satellite provided the first look at the chemical composition of ices in the mid-IR in dense molecular clouds and cores (Whittet et al. 1998; Dartois 2005; Knez et al. 2005; Whittet et al. 2009; Boogert et al. 2011; Noble et al. 2013), which included information on additional molecular species such as CO₂ or OCN⁻. However, the sensitivity or spectral resolving power of these observations was rather limited, which allowed only the detection of simple and abundant species towards the lines of sight in molecular clouds with moderate levels of visual extinction (A_V ≤ 30 mag; see e.g. the summary by Boogert et al. 2015).

The James Webb Space Telescope (JWST), with its unprecedented sensitivity and high spectral resolving power (up to ~3000), is providing the most comprehensive view of the chemical composition of interstellar ices during the initial conditions of star and planet formation. The Ice Age Director's Discretionary Early Release Science programme has obtained ice absorption spectra towards hundreds of field stars located behind the Chameleon I dense molecular cloud (McClure et al. 2023). The goal of the project is to establish the chemical composition of interstellar ices concurrently with the process of solar-type system formation, from the molecular cloud stage to pre-stellar cores, young stellar objects (Class 0 and I sources), and, ultimately, protoplanetary disks (Class II-type sources; McClure et al. 2023; Sturm et al. 2023, 2024; Dartois et al. 2024; Noble et al. 2024). The first results have reported the ice chemical composition towards the most highly extinguished background stars observed to date in two molecular clouds, NIR38 and J110621, with pre-flight estimated visual extinctions of A_V = 60 mag and A_V = 95 mag, respectively (McClure et al. 2023), although these values are subject to large uncertainties (see Dartois et al. 2024; Noble et al. 2024, and Sect. 4.7). These observations reveal that the ice composition is very similar towards the two lines of sight: CO ice is ~30% and ~45%, respectively, relative to water; and CO is more abundant than CO₂ (between 13% and 20% with respect to water) and methanol (between 4 and 9%). Hydrides such as CH₄ and NH₃ comprise only up to 5% of ice with respect to water.

Based on the expected total CO abundance in molecular clouds for these levels of extinction (Lacy et al. 2017), it was proposed that the carbon-bearing ices observed with the JWST accounted for ≤50% of the total molecular CO in gas or ice (McClure et al. 2023). This scenario was consistent with the prominent absorption feature in the red wing of the water 3µm band, which was attributed to the presence of hydrides suggesting a water-rich ice environment (see Extended Data Fig. 9 in McClure et al. 2023). Under these conditions, the formation of methanol at 10–20 K is thought to be driven by the grain surface reaction chain CH₄ + OH → CH₃ + H₂O and CH₃ + OH → CH₃OH, as proposed by laboratory experiments and quantum chemical calculations of ice formation (Qasim et al. 2018; Molpeceres et al. 2021). Additionally, the CO₂ ice feature profile was better matched with laboratory spectra of a CO₂–H₂O ice blend (Extended Data Fig. 4; see McClure et al. 2023), and weak features between 7 and 7.5 µm were consistent with ethanol blended with H₂O suggesting a similar formation environment for these species. The formation of hydrides such as CH₄ and NH₃ could also proceed in water-rich ice environments under the same cold interstellar conditions as shown by laboratory experiments (Hiraoka et al. 1995; Qasim et al. 2020; Lamberts et al. 2022), which supports the scenario of efficient CH₃OH ice formation in a water-ice rich environment towards the lines of sight of NIR38 and J110621.

This idea, however, has been re-evaluated recently, and new radiative transfer calculations show that much of the red wing of the 3 µm water band observed with JWST towards these two field stars, could alternatively be due to grain growth and not hydrides, suggesting that ices in the Chamaeleon I cloud may not be as water-rich as initially thought (Dartois et al. 2024). In that case, the production of methanol via the reaction CH₃ + OH would be less efficient, and thus other chemical mechanisms for methanol formation (as e.g. the successive hydrogenation of CO; see Watanabe & Kouchi 2002) are expected to take over. It remains unclear, however, to what amount these two mechanisms contribute to the production of methanol in ices, and hence dedicated chemical modelling is needed to investigate the possible and relative importance of different chemical routes for the production of methanol ice.

In a recent work, Jin et al. (2022) modelled the ice chemistry towards the Class 0 source Cha-MMS1, also observed within the Ice Age programme, where the gas densities reach values ≥10¹⁰ cm⁻³ and gas and dust temperatures can be as high as ≥100 K. However, this modelling did not focus on the formation of ices towards the quiescent NIR38 and J110621 sightlines in the Chamaeleon I molecular cloud, with typical gas densities of ~10⁴–10⁶ cm⁻³ and temperatures of ~10–15 K, and therefore the chemical modelling of the JWST data towards these two background stars is still missing. In addition, while different chemical codes do exist at present, they are built following different approaches and techniques (rate-equation approach vs. stochastic Monte Carlo techniques) or they consider different chemical processes on the surface of dust grains (diffusive vs. non-diffusive reaction mechanisms and thermal and non-thermal desorption processes; see e.g. the review by Cuppen et al. 2017). The Ice Age results therefore provide a unique opportunity to compare different astrochemical codes so that the dominant chemical routes yielding the ice chemical inventory detected with JWST towards the Chamaeleon I molecular cloud can be constrained.

For this work we modelled the chemical composition of the ices observed with the JWST towards the background stars NIR38 and J110621 within the Ice Age ERS programme. For this modelling effort, we used five different astrochemical codes to identify the features (or chemical processes) within these models that best reproduce the observed ice chemical compositions. We note that the last community effort to compare astrochemical models goes back to the review work of Cuppen et al. (2017).

The paper is organised as follows. In Sect. 2 we describe the astrochemical models used and the parameter space of the models used in this work. In Sect. 3 we report the results of the models and compare the ice molecular abundances predicted by the models with those measured with the JWST. In Sect. 3.3, we identify the models from each astrochemical code that best reproduce the ice chemical compositions observed within Ice Age. Section 3.4 discusses how our results change if lower dust temperatures are considered in our best models. In Sect. 4, we present the general trends obtained from the models and discuss the uncertainties and caveats. Finally, our conclusions are summarised in Sect. 5.

## 2. Description of the chemical models

Our aim is to model the ice chemical formation and evolution in the Chameleon I dense molecular cloud using different astrochemical codes. These codes were built following different philosophies. We have rate equation codes, such as MAGICKAL (Garrod 2019), MONACO (Vasyunin & Herbst 2013; Vasyunin et al. 2017), Nautilus (Ruaud et al. 2016), and UCLCHEM
<table><caption>**Table 1.** Main features of the astrochemical codes used in this work.</caption><thead><tr><th></th><th>MAGICKAL</th><th>MONACO</th><th>Nautilus</th><th>UCLCHEM</th><th>KMC</th></tr><tr><th>Code type</th><th>Rate Equation</th><th>Rate Equation</th><th>Rate Equation</th><th>Rate Equation</th><th>Stochastic</th></tr></thead><tbody><tr><td><strong>Thermal desorption</strong></td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td></tr><tr><td><strong>Non-thermal desorption</strong></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>UV photodesorption</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td></tr><tr><td>CR-induced photodesorption</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td></tr><tr><td>CR-induced heating</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td></tr><tr><td>Chemical (reactive) desorption</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td></tr><tr><td>H<sub>2</sub>-formation-induced desorption</td><td>–</td><td>–</td><td>–</td><td>✓</td><td>–</td></tr><tr><td>CR sputtering</td><td>–</td><td>–</td><td>✓</td><td>–</td><td>–</td></tr><tr><td><strong>Diffusive chemistry</strong></td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td></tr><tr><td><strong>Non-diffusive chemistry</strong></td><td>✓</td><td>✓</td><td>–</td><td>–</td><td>✓</td></tr></tbody></table>

<table><caption>**Table 2.** Physical parameters assumed for the grid of models.</caption><thead><tr><th rowspan="2">Model</th><th colspan="4">NIR38</th><th colspan="4">J110621</th></tr><tr><th>n<sub>H</sub><br>(cm<sup>-3</sup>)</th><th>T<sub>dust</sub><br>(K)</th><th>A<sub>V</sub><br>(mag)</th><th>Collapse time<br>(yr)</th><th>n<sub>H</sub><br>(cm<sup>-3</sup>)</th><th>T<sub>dust</sub><br>(K)</th><th>A<sub>V</sub><br>(mag)</th><th>Collapse time<br>(yr)</th></tr></thead><tbody><tr><td>1</td><td>2 × 10<sup>4</sup></td><td>14</td><td>60</td><td>10<sup>4</sup></td><td>2 × 10<sup>4</sup></td><td>12</td><td>95</td><td>10<sup>4</sup></td></tr><tr><td>2</td><td>2 × 10<sup>4</sup></td><td>14</td><td>60</td><td>10<sup>5</sup></td><td>2 × 10<sup>4</sup></td><td>12</td><td>95</td><td>10<sup>5</sup></td></tr><tr><td>3</td><td>2 × 10<sup>4</sup></td><td>14</td><td>60</td><td>10<sup>6</sup></td><td>2 × 10<sup>4</sup></td><td>12</td><td>95</td><td>10<sup>6</sup></td></tr><tr><td>4</td><td>2 × 10<sup>5</sup></td><td>14</td><td>60</td><td>10<sup>4</sup></td><td>2 × 10<sup>5</sup></td><td>12</td><td>95</td><td>10<sup>4</sup></td></tr><tr><td>5</td><td>2 × 10<sup>5</sup></td><td>14</td><td>60</td><td>10<sup>5</sup></td><td>2 × 10<sup>5</sup></td><td>12</td><td>95</td><td>10<sup>5</sup></td></tr><tr><td>6</td><td>2 × 10<sup>5</sup></td><td>14</td><td>60</td><td>10<sup>6</sup></td><td>2 × 10<sup>5</sup></td><td>12</td><td>95</td><td>10<sup>6</sup></td></tr><tr><td>7</td><td>2 × 10<sup>6</sup></td><td>14</td><td>60</td><td>10<sup>4</sup></td><td>2 × 10<sup>6</sup></td><td>12</td><td>95</td><td>10<sup>4</sup></td></tr><tr><td>8</td><td>2 × 10<sup>6</sup></td><td>14</td><td>60</td><td>10<sup>5</sup></td><td>2 × 10<sup>6</sup></td><td>12</td><td>95</td><td>10<sup>5</sup></td></tr><tr><td>9</td><td>2 × 10<sup>6</sup></td><td>14</td><td>60</td><td>10<sup>6</sup></td><td>2 × 10<sup>6</sup></td><td>12</td><td>95</td><td>10<sup>6</sup></td></tr></tbody></table>

**Notes.** $T_{\text{gas}}$ is assumed to be the same as $T_{\text{dust}}$ in all models.

(Holdship et al. 2017). In addition, we have the stochastic Kinetic Monte Carlo (KMC) model of Cuppen & Herbst (2007). These codes also include different levels of complexity for the formation and chemistry of the ices, and thus the idea of this work is to test the features from these astrochemical codes that best reproduce the JWST observations obtained towards NIR38 and J110621, which are field stars with the highest extinction observed thus far (McClure et al. 2023). For instance, while MAGICKAL, MONACO and KMC include diffusive and non-diffusive chemistry, Nautilus and UCLCHEM only include diffusive chemistry. On the other hand, Nautilus and UCLCHEM include some non-thermal desorption mechanisms that are not considered in the other codes. See Table 1 for the main characteristics and features of each of these codes. In this section we briefly describe the astrochemical codes employed in this work (listed in alphabetical order) and we report the physical parameter space used to run the models.

## 2.1. Kinetic Monte Carlo (KMC) simulations

The KMC method denotes a class of algorithms aimed at modelling the time evolution of a given system by generating a stochastic trajectory using random sampling. This is the only code of the stochastic type used in this work. This method provides detailed information about the microscopic processes taking place on the surface of dust grains. For an overview of the theoretical background and its application to interstellar ice layers, see Cuppen et al. (2013).

In the present study, a microscopic type of KMC was used called “continuous-time random-walk kinetic Monte Carlo” (CTRW-KMC). This particular simulation program is an adaptation from the one used in Cuppen & Herbst (2007), Cuppen et al. (2009), Lamberts et al. (2014), and Simons et al. (2020). The program exclusively simulates grain surface chemistry, but the results can be coupled to the gas by submitting gas-phase abundances as input, from which the adsorption fluxes are calculated. For the present paper, the code was adapted to be able to treat densities and gas abundances that vary with time. These were simulated using the UCLCHEM code, where collapse was assumed to follow free-fall and the input parameters were varied in accordance with the physical parameter grid of Table 2 and described in Sect. 2.6. The resulting time-dependent densities and a selection of gas abundances were then passed on to the KMC grain code. We note that there is no feedback from the thermal desorption of grain species onto the gas abundances, as we assume the effect to be negligible at these low temperatures. The selection of gas abundances transferred to the KMC model was limited to a small set of species: H, CO, O, C, and N, i.e. the atomic and molecular species that contribute to the formation of the main ice species studied here (CO<sub>2</sub>, CH<sub>3</sub>OH, and the hydrides). The selection was kept minimal, because we do not expect the results from the UCLCHEM grain surface chemistry to conflate with those from the KMC simulations. This is due to the low dust temperatures and UV radiation fields considered in our models, which induce little desorption from the grain surface.
The surface of the grain is represented by a 50 × 50 lattice with periodic boundary conditions upon which species can accrete on top and layers can form. Apart from thermally activated processes such as diffusion (hopping) and desorption, the lattice KMC model also includes the non-thermal processes of direct photodissociation and cosmic-ray dissociation. Diffusion and desorption are made site-specific processes where the binding energy is neighbour-dependent. The sequence of processes and corresponding time advances are determined in a stochastic manner from the rates. See Cuppen & Herbst (2007) for a more detailed account of the calculation of the rate coefficients for these various processes. The chemical reaction networks used were confirmed experimentally or by quantum mechanical calculations in Lamberts et al. (2014) (water network), Simons et al. (2020) (CO network) and Ioppolo et al. (2020) (glycine network).

## 2.2. MAGICKAL code

The MAGICKAL code is a three-phase astrochemical model that solves a set of rate equations describing chemical kinetics in the gas, on the grain surface, and within the bulk of the ice mantle. Grain-surface chemistry is coupled with the gas through accretion and desorption of chemical species. The latter can proceed both thermally and non-thermally. Non-thermal desorption mechanisms include chemical (reactive) desorption, and photo-desorption by both external and CR-induced UV photons.

In this model, grain-surface chemistry is controlled by both diffusive and non-diffusive mechanisms. The former represents the chemical kinetics described by the thermal diffusion of surface species by which surface species may meet and react. MAGICKAL uses the modified-rate method of Garrod (2008) for diffusive surface reactions (see also Garrod et al. 2009); this method adjusts reaction rates to replicate the stochastic behaviour of the system when particular reactions fall into the so-called “small-grain” limit, i.e. when reactant populations are small and the diffusion rate is large. In these cases, the reaction rate should be limited by the accretion (or production) rate of the reactants on the surface.

The diffusion rate is governed by the temperature of the dust and the diffusion barrier ($E_{diff}$) of the relevant species. The diffusion barriers of the surface species are assumed to be some constant fraction of their binding energies to the surface ($E_{des}$), initially based on work on metal surfaces where often a ratio between the diffusion barrier and the binding energy of 0.3 was observed. However, astrochemically relevant systems appear to show a larger scatter for these data based on the few systems for which both the diffusion barrier and the binding energy are known (Karssemeijer & Cuppen 2014; Yang et al. 2022). Because of this, astrochemical models usually assume a fixed value in the approximate range of 0.3–0.8. In MAGICKAL, we assumed 0.35 for the molecules, and 0.55 for other atomic species, as these values best reproduce the production of interstellar CO$_2$ ice (Garrod & Pauly 2011). For atomic H, average binding energy and diffusion barrier values of 661 K and 243 K, respectively, are used, based on the calculations of Senevirathne et al. (2017). The bulk diffusion barriers for H and H$_2$ are twice the values used for surface diffusion, as per Garrod (2013).

The kinetics of non-diffusive mechanisms consider the probability by which a newly formed chemical species reacts with other potential reactants in the proximity without diffusion. This mechanism is important for accurately treating the production of complex organic molecules (COMs) at low temperatures (Garrod 2019; Jin & Garrod 2020). For the details of the four

non-diffusive mechanisms included in MAGICKAL, we refer to Garrod et al. (2022).

The surface layer is coupled to the bulk ice phase through the net gain/loss rate of material to/from the surface, which either acts to cover up the existing surface or to expose the underlying bulk ice material (Hasegawa & Herbst 1993; Garrod & Pauly 2011). The chemical kinetics within the bulk ice is controlled by both tunnelling-mediated bulk diffusion for H and H$_2$ (see Garrod et al. 2022) and non-diffusive chemistry. The chemical networks used in these models are based on that presented by Garrod et al. (2022).

## 2.3. MONACO code

The gas-grain astrochemical model MONACO (Vasyunin et al. 2017; Jiménez-Serra et al. 2021, Borscheva et al., in prep.) is designed to simulate chemical evolution under interstellar conditions using a rate equations approach. Similar to MAGICKAL, the MONACO astrochemical code is also a three-phase model that considers the gas phase, the surface of the icy mantle and the bulk of the ice. We note that the surface of the icy mantle corresponds to the four uppermost surface layers of the ice.

Chemistry in the gas phase and on grains is coupled via accretion and desorption. Four types of desorption processes are included in the MONACO code: thermal sublimation, photodesorption, cosmic ray-induced desorption and reactive (chemical) desorption. Although several treatments of reactive desorption are included in the code (following Garrod et al. 2007; Vasyunin & Herbst 2013; Minissale et al. 2016), here we employ the approach of Garrod et al. (2007).

In the code, diffusive and non-diffusive chemical processes on grains are both considered. Diffusive chemistry is implemented following the classical approach by Hasegawa et al. (1992) with several modifications. First, the treatment of reaction/diffusion competition is added following Chang et al. (2007). Next, time-dependent adjustment of binding energies of species on the surface due to the presence of H$_2$ molecules is added following Garrod & Pauly (2011). We note, however, that while Garrod & Pauly (2011) adjusted binding energies for all species, we do it only for H and H$_2$. Heavier species are considered to penetrate through the top monolayer of H$_2$ molecules on ice surface during accretion, and bind to the regular ice surface beneath. For diffusive chemistry, only thermal hopping is considered as a source of mobility of species. Diffusion-to-binding energy ratios are assumed different for atomic and molecular species. For atoms, $E_{diff}/E_{des} = 0.5$, while for molecules $E_{diff}/E_{des} = 0.3$. In the bulk, diffusion energies are twice higher than on the surface for all species except for H and H$_2$. For H and H$_2$, their bulk diffusion energies are 1.5 times higher than values in the ice surface. Tunnelling through reaction activation barriers is allowed.

Non-diffusive chemistry is considered following the approach proposed by Jin & Garrod (2020). Excited three-body reactions are excluded. Ice photochemistry and chemistry induced by cosmic rays is treated as in Vasyunin et al. (2017), although rates of some processes are updated according to recent estimates (see e.g. Jiménez-Serra et al. 2018, 2021).

## 2.4. Nautilus code

Nautilus is an astrochemical gas-grain code based on the rate equation method and available through a git repository<fnref n="1" /> (Ruaud

<footnote n="1">¹ https://forge.oasu-u-bordeaux.fr/LAB/astrochem-tools/pnautilus</footnote>
et al. 2016; Wakelam et al. 2024). It is also a three phase model meaning that it takes into account the gas-phase chemistry and it makes a distinction between the few upper layers (four monolayers) of molecules and the rest of the mantles (bulk) on top of interstellar refractory grains. The diffusion of the species is less efficient in the bulk than on the surface. The ratio between the diffusion energies and binding energies is assumed to be 0.8 in the bulk and 0.4 on the surface. Diffusion is assumed to occur via thermal hopping but also through tunnelling depending on the mass of the species. Grains are assumed to be of one size and composition (silicate grains of 0.1 µm in radius). Non-diffusive chemistry is not considered in Nautilus.

In addition to thermal desorption, a number of non-thermal desorption processes are included: whole grain heating induced by cosmic-rays (Hasegawa & Herbst 1993), photodesorption (Ruaud et al. 2016), chemical reactive desorption (Minissale et al. 2016), and sputtering by cosmic-rays (Wakelam et al. 2021). We note that this is the only code in this work that considers the effects of cosmic-rays sputtering as a desorption mechanism following the experimental results of Dartois et al. (2020) and Dartois et al. (2021). The chemical networks and parameters are the same as in Clément et al. (2023).

## 2.5. UCLCHEM code

UCLCHEM (Holdship et al. 2017) is a time-dependent gas-grain chemical model based on the rate equation approach. UCLCHEM is a public open source code<fnref n="2" /> that considers the gas phase as well as the bulk and the surface of dust grains (i.e. it is also a three-phase model). Gas phase reactions are two-body reactions, including interactions with cosmic rays and UV photons. Gas species can accrete onto grain surfaces. Diffusive chemistry is considered as described in Quénard et al. (2018b), using a constant ratio between diffusion and binding energies of 0.5. However, UCLCHEM does not include non-diffusive processes in the ices. Both thermal and non-thermal desorption mechanisms are considered in this code. For thermal desorption, UCLCHEM follows the formalism of Viti et al. (2004), while for the non-thermal desorption processes UCLCHEM includes desorption produced by the energy released in H₂ formation, incident cosmic rays, UV photo-desorption, cosmic-ray-induced UV photo-desorption (as described in Holdship et al. 2017) and chemical reactive desorption following the formalism of Minissale et al. (2016) (see Quénard et al. 2018b). The default grain chemical network is based on the one developed by Quénard et al. (2018b), and the gas-phase chemical network is taken from the UMIST 2013 database (McElroy et al. 2013). The default grain chemical network was updated by adding the new reactions for CH₃OH ice formation proposed in Sect. 2.6 and by changing the activation energy of reaction H₂CO + O → CO₂ + H₂ from 0 to 335 K (Minissale et al. 2015). We also added a set of reactions involving HOCO and CO₂, which are relevant for the production of CO₂ in the ice (see Appendix A and references therein). We note that these reactions were already included in the other codes.

## 2.6. Physical parameters and chemical reactions explored by the models

The chemistry of the Chamaeleon I molecular cloud is modelled considering its evolution from its translucent phase to the dense cloud stage. To do this, we first consider a “phase 0” stage in which the chemistry of a translucent cloud is modelled for a

duration of 10⁶ yr assuming that the cloud is static and adopting a number density of total H nuclei of $n_H = 2 \times 10^3$ cm⁻<fnref n="3" />, a dust and gas temperature $T_d = T_{gas} = 15$ K, visual extinction $A_V = 2$ mag, UV radiation field $G_0=1$ hab, and a cosmic-ray ionisation rate $\zeta = 1.3 \times 10^{-17}$ s⁻¹.

The chemical evolution of the cloud collapse is modelled in a subsequent “phase 1” stage, in which the output gas and ice (both surface and bulk) abundances from “phase 0” are used as input abundances for the “phase 1” simulations. For the cloud collapse, we assume a free-fall collapse treatment with the gas density ($n_H$) evolving with time ($t$) as (Brown et al. 1988; Rawlings et al. 1992)

$$ \frac{d n_H}{d t} = b \left( \frac{n_H^4}{n_{H,0}} \right)^{\frac{1}{3}} \left[ 24 \pi G m_H n_{H,0} \left( \left( \frac{n_H}{n_{H,0}} \right)^{\frac{1}{3}} - 1 \right) \right]^{\frac{1}{2}} \quad (1) $$

where $n_{H,0}$ is the initial hydrogen nuclei central density, $m_H$ is the mass of a hydrogen nucleus, and $G$ is the gravitational constant. We also include a factor $b$ in the equation so that we can test the effects of an accelerated collapse in the chemistry of the ices (Holdship et al. 2017). The collapse finishes when $n_H$ reaches the final density values at the final collapse times shown in Table 2. As the density $n_H$ evolves during the collapse, the visual extinction $A_V$ increases as a function of density (and therefore time), following Eq. (1) of Jiménez-Serra et al. (2018).

Table 2 reports the grid of physical parameters considered for the models of the two field stars NIR38 and J110621. The assumed visual extinctions are $A_V = 60$ mag and $A_V = 95$ mag, respectively; however, these values may be lower (see Dartois et al. 2024, and Sect. 4.7). The considered gas/dust temperatures are 14 K for NIR38 and 12 K for J110621, as inferred from the Herschel-based dust temperatures maps of Väisälä et al. (2014)<fnref n="3" />. For the MAGICKAL and MONACO codes, the dust temperature is assumed to decrease from 15 K during the core collapse with increasing visual extinction as described in Garrod & Pauly (2011). If the temperature reaches the minimum value of the line of sight (12 K or 14 K depending on the field star; Table 2), the temperature is then fixed. For the rest of the codes, $T_{dust}$ is kept constant at either 12 K or 14 K (see Table 2). The gas temperature is set to be the same as the dust temperature in all runs and for all codes.

The final $n_H$ densities at the end of the “phase 1” collapse stage are $2 \times 10^4$, $2 \times 10^5$ and $2 \times 10^6$ cm⁻<fnref n="3" />. These densities have been selected from the H₂ volume gas densities estimated by Belloche et al. (2011) towards the Cha1-C1 and Cha1-C2 cores (see their Table 6), which are consistent with the results of Väisälä et al. (2014). These densities lie between some $10^5$ and $10^6$ cm⁻<fnref n="3" />. For completeness, we also explore the lower-density case with $n_H = 2 \times 10^4$ cm⁻<fnref n="3" />.

For the final collapse times, we assume timescales of $10^4$, $10^5$ and $10^6$ yr. Free-fall collapse times<fnref n="4" /> correspond to ~$10^6$ yr for an initial gas cloud density of $n_{H,0} = 2 \times 10^3$ cm⁻<fnref n="3" /> (i.e. the one used in the “phase 0” stage of the simulations). The collapse times of $10^4$ and $10^5$ yr are tested by modifying the $b$ parameter in Eq. (1), which adopts values from 1 (free-fall, ~$10^6$ yr) to 120 (accelerated collapse, ~$10^4$ yr). The three different collapse times assumed in our grid of models are included to test the effects of a fast or slow collapse on the chemical composition

<footnote n="3">³ The uncertainties associated with these $T_{dust}$ values are as high as ±0.6 K, after considering a $\beta$ range between 1.8 and 2.2 and the uncertainties due to photometric noise (Väisälä et al. 2014).</footnote>

<footnote n="4">⁴ Free-fall collapse times are calculated as $t_{free-fall} \approx 4 \times 10^7 / \sqrt{n_{H,0}}$ yr, with $n_{H,0}$ the initial gas cloud density.</footnote>

<footnote n="2">² https://uclchem.github.io/</footnote>
<table><caption>**Table 3.** Initial elemental abundances with respect to $n_H$.</caption><thead><tr><th>Species</th><th>Abundance</th></tr></thead><tbody><tr><td>He</td><td>1.40 × 10<sup>-1</sup></td></tr><tr><td>N</td><td>2.14 × 10<sup>-5</sup></td></tr><tr><td>O</td><td>1.76 × 10<sup>-4</sup></td></tr><tr><td>C</td><td>7.30 × 10<sup>-5</sup></td></tr></tbody></table>

**Notes.** Abundances correspond to the EA1 low-metal abundance case of Wakelam & Herbst (2008). These values consider the partial depletion of these elements into dust grains.

of the ices. We note that short timescales of the order of $10^4$ yr are typically not observed for pre-stellar cores (Keto & Caselli 2010). However, accelerated collapse is considered in this study to investigate the possible effects of gravitational collapse driven by stellar feedback from nearby star formation<fnref n="5" />. We note that we deliberately do not explore a larger grid of models because the KMC simulations are computationally very expensive. The cosmic-ray ionisation rate and UV radiation field assumed for the “phase 1” stage in all models are, respectively, $\zeta=1.3 \times 10^{-17}$ s⁻¹ and $G_0=1$ hab.

As initial elemental abundances, we consider the low-metal abundance case (the “EA1” case) from Wakelam & Herbst (2008, see Table 3). All models consider a constant grain size of 0.1 µm. The output fractional abundances are given with respect to the number density of total H nuclei. The conversion from ice molecular abundances to ice molecular column densities is performed using the column densities of total H nuclei towards NIR38 and J110621 calculated<fnref n="6" /> as $N_H = 1.6 \times 10^{21} \times A_V \text{ cm}^{-2}$.

For the formation of methanol, in addition to the classical route of multiple CO hydrogenation reactions in the ice, we also include the following chemical processes: 1) the chain of neutral-neutral grain surface reactions $CH_4 + OH \rightarrow CH_3 + H_2O$ and $CH_3 + OH \rightarrow CH_3OH$ proposed by Qasim et al. (2018); and 2) the reaction $C + H_2O \rightarrow H_2CO$ and subsequent hydrogenation reactions proposed by Molpeceres et al. (2021) and Potapov et al. (2021).

### 3. Results

In this section, we analyse the results obtained with each code, presenting the evolution of the chemical composition of the ices and the main processes that lead to their formation for the grid of models described in Sect. 2.6. We focus on some of the molecular species observed with JWST towards NIR38 and J110621, namely $H_2O$, CO, $CO_2$, $CH_3OH$, NH₃, and CH₄. The model results are thus compared to the observational data reported in McClure et al. (2023) and presented in Table 4. In Sect. 3.1, we compare the output from all models for their phase 0 runs to identify any possible difference arising from this stage of the simulations. In Sects. 3.2.1, 3.2.2, 3.2.3, and 3.2.4 we present the results for the phase 1 runs from the rate equation codes MAGICKAL, MONACO, Nautilus, and UCLCHEM, respectively, while in Sect. 3.2.5, we report the results from the stochastic KMC model.

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Time evolution of the gas and ice molecular abundances for most relevant species predicted by the different astrochemical codes for the phase 0 runs. The solid lines indicate gas-phase abundances, while the dashed lines show ice abundances. The colours indicate the different molecular species (see legend at top).</figcaption></figure>

#### 3.1. Comparison of the Phase 0 results from all models

Figure 1 shows the abundances with respect to $H_2$ resulting from the phase 0 stage of the simulations, where we report the time evolution of the relevant species modelled by the different astrochemical codes. For the gas-phase, we only plot the relative abundances of atomic C and CO since these two species are essential to understanding the amount of C that will be locked into CO, $CO_2$, and $CH_4$ ices during the phase 1 of the simulations. We note that we do not show the phase 0 results for the KMC model, because they are similar to the ones from the UCLCHEM code. The largest differences found between the UCLCHEM and KMC phase 0 results are the abundances of $CO_2$ for which several reactions have been added within UCLCHEM (see Appendix A).

From Fig. 1, it is clear that for all models atomic C has been efficiently converted into CO by the end of the phase 0 stage with a CO gas-phase abundance that goes up to ~$7 \times 10^{-5}$ with respect to the total H nuclei at the end of the simulations. As shown in Sects. 3.2.1, 3.2.2, 3.2.3, 3.2.4, this will largely affect the subsequent production of $CO_2$ and $CH_4$ during the phase 1 stage of the models. For the Nautilus and MAGICKAL codes, the temporal evolution of the gas and ice molecular abundances are similar, although some differences can be seen for $CH_3OH$ and $CH_4$. In these models, $H_2O$ and $CO_2$ are the most abundant ice species with relative abundances lying between $10^{-7}$ and $10^{-6}$. Due to the dust temperature of $T_{dust}= 15$ K considered for the phase 0 runs, any CO molecule sticking on the grains is converted into $CO_2$, which explains the high abundance of $CO_2$ ice (even higher than that of CO) for these models. For the MONACO code, the final abundances of the phase 0 runs are somewhat higher than for the other codes (see Fig. 1), with an ice thickness of ~10 monolayers, versus ~1 monolayer or less.

For the UCLCHEM code, Fig. 1 shows that the ices grow at a slower effective rate compared to the other codes. This is produced by the higher efficiency of non-thermal desorption in UCLCHEM compared to other codes. In particular, UV-photodesorption is the mechanism responsible for the slow net ice growth at these low visual extinctions ($A_V=2$ mag). We note, however, that if the assumed UV yield in UCLCHEM (i.e. the

<footnote n="5">⁵ A 10 km s⁻¹ shock propagating through a pre-shock gas with density of ~$10^4$ cm⁻³ would induce a density compression in the postshock gas by only a factor ~3 for timescales ~$10^4$ yr (see Jiménez-Serra et al. 2008; Gusdorf et al. 2008).</footnote>

<footnote n="6">⁶ This expression differs only by a factor of ~1.1 with respect to the one reported by Lacy et al. (2017) of $N_H = 1.8 \times 10^{21} A_V \text{ cm}^{-2}$.</footnote>
<table><caption>**Table 4.** Observed column densities and column density ratios with respect to H₂O towards NIR38 and J110621 (McClure et al. 2023).</caption><thead><tr><th rowspan="3">Molecule</th><th colspan="2">NIR38<br>(A<sub>V</sub> = 60 mag)</th><th colspan="2">J110621<br>(A<sub>V</sub> = 95 mag)</th></tr><tr><th>N(molecule)<br>(10<sup>18</sup> cm<sup>-2</sup>)</th><th>N(molecule)/N(H<sub>2</sub>O)<br>(in %)</th><th>N(molecule)<br>(10<sup>18</sup> cm<sup>-2</sup>)</th><th>N(molecule)/N(H<sub>2</sub>O)<br>(in %)</th></tr></thead><tbody><tr><td>H<sub>2</sub>O</td><td>6.9<sup>+1.9</sup><sub>-1.1</sub></td><td>100</td><td>13.4<sup>+1.3</sup><sub>-1.8</sub></td><td>100</td></tr><tr><td>CO</td><td>3.0<sup>+0.6</sup><sub>-0.4</sub></td><td>43<sup>+12</sup><sub>-10</sub></td><td>3.7<sup>+0.6</sup><sub>-0.4</sub></td><td>28<sup>+6</sup><sub>-4</sub></td></tr><tr><td>CO<sub>2</sub></td><td>1.38<sup>+0.20</sup><sub>-0.20</sub></td><td>20<sup>+5</sup><sub>-5</sub></td><td>1.74<sup>+0.21</sup><sub>-0.22</sub></td><td>13.1<sup>+2.5</sup><sub>-2.0</sub></td></tr><tr><td>CH<sub>3</sub>OH</td><td>0.61<sup>+0.11</sup><sub>-0.11</sub></td><td>8.7<sup>+2.5</sup><sub>-2.1</sub></td><td>0.51<sup>+0.19</sup><sub>-0.09</sub></td><td>3.9<sup>+1.5</sup><sub>-0.8</sub></td></tr><tr><td>NH<sub>3</sub></td><td>0.30<sup>+0.22</sup><sub>-0.03</sub></td><td>4.6<sup>+3.0</sup><sub>-1.1</sub></td><td>0.66<sup>+0.15</sup><sub>-0.06</sub></td><td>5.0<sup>+1.3</sup><sub>-0.7</sub></td></tr><tr><td>CH<sub>4</sub></td><td>0.180<sup>+0.017</sup><sub>-0.013</sub></td><td>2.6<sup>+0.5</sup><sub>-0.6</sub></td><td>0.250<sup>+0.010</sup><sub>-0.030</sub></td><td>1.84<sup>+0.30</sup><sub>-0.24</sub></td></tr></tbody></table>

number of molecules desorbed per UV photon) is reduced from 0.1 molecules photon⁻¹ to 0.03 molecules photon⁻¹, the predicted ice abundances become similar to those given by the other codes. For the remaining of this work, we use the default value of 0.1 molecules photon⁻¹ given in UCLCHEM, to evaluate the effects of a less efficient ice growth during the static phase 0 stage in a translucent cloud (see Sect. 3.2.4).

In the following sections, we present the results from the phase 1 calculations: first, for all rate equation codes (presented in alphabetical order); and second, for the KMC code.

## 3.2. Collapse Phase 1 results

### 3.2.1. Results from the MAGICKAL code

In Fig. 2, we plot the time evolution of the gas and ice relative abundances of H₂O, CO, CO₂, CH₃OH, CH₄, and atomic C during the phase 1 stage of the collapse, for four representative models obtained with the MAGICKAL code. As examples, we report the most extreme cases with the lowest-density and fastest-collapse models (models #1; left panels) and the highest-density and slowest-collapse models (models #9; right panels) for NIR38 and J110621, to provide information about the maximum changes expected within the grid of models. In the low-density and fast-collapse models (left panels, Fig. 2), the chemical evolution is nearly identical for NIR38 and J110621 regardless of temperature, except for CH₃OH ice. The gas phase abundances of CO remain largely unchanged over time due to the collapse timescale being much shorter than the accretion timescale of CO. Consequently, this leads to almost constant abundances of CO ice, and thus CO₂, as CO₂ primarily forms from the association of CO and OH on the surface. However, in the model of J110621 with T<sub>dust</sub>=12 K, hydrogenation becomes more efficient, resulting in the formation of more CH₃OH via the hydrogenation of CO.

In the models with high density and slow collapse (right panels; Fig. 2), gas-phase CO and C efficiently deplete on grains as density evolves, significantly impacting surface chemistry. As can be seen from this figure, the abundance of CO ice remains relatively constant (or even decreases) until the density of the gas becomes high enough for the catastrophic depletion of CO to take place. Surface CO then rapidly reacts forming primarily CO₂ through the grain surface reaction CO + OH → CO₂ + H (see Garrod & Pauly 2011; Clément et al. 2023), where OH forms via the hydrogenation of adsorbed atomic oxygen. Once

CO₂ is formed, it accumulates in the ice due to its chemically inert characteristics. Consequently, the ice abundance of CO₂ in these models gradually increases over time, while CO does not, resulting in a CO₂ / CO ratio > 1. The abundance of CH₃OH ice exhibits an evolutionary trend similar to that of CO since it primarily forms from the hydrogenation of CO on the surface. In the model with T<sub>dust</sub>=12 K, methanol is more abundant due to the higher abundance of CO on the grain surface and the higher efficiency of CO hydrogenation. The contribution of the reactions CH₃ + OH and C + H₂O to the production of methanol in these models is minor compared to grain-surface CO hydrogenation. This is due to the fact that atomic C is mainly converted into gas-phase CO, which quickly freezes out onto dust grains to form either CH₃OH or CO₂. This leaves little atomic C in the gas phase to adsorb on grains and subsequently hydrogenate forming CH₄, or to react with H₂O in the ice forming H₂CO (see also Sects. 4.2 and 4.4).

The final ice chemical compositions towards NIR38 and J110621 for the whole grid of models can be found in Figs. 3 and 4. In these figures, we provide the absolute column densities of H₂O predicted in the ice at the end of the phase 1 simulations, together with the ice abundance ratios of the rest of species analysed in this study (CO, CO₂, CH₃OH, NH₃, and CH₄) with respect to H₂O in %. From Figs. 3 and 4, we find, regardless of the model used, that both lines of sight consistently underpredict water ice column densities and the ice abundances of CH₄. However, the water ice column density of the model with the longest collapse timescale (t<sub>col</sub> = 10⁶ yr) and the highest density (n<sub>H</sub> = 10⁶ cm⁻³) aligns with observations within the margin of error. This is due to the fact that models with longer collapse timescales enable the accumulation of larger amounts of H₂O in the ice mantle. CO, CO₂, and CH₃OH are the most sensitive species to variations depending on the model parameters. Ice abundances of CO₂ surpass observational values in all models for both lines of sight. CO₂ becomes more abundant in the ice than even CO yielding CO₂ / CO ratios ≥ 1, which is not observed towards either the field star NIR38 or J110621. For NH₃ and CH₄, we find similar ice abundances across all models given that both species are formed through the consecutive hydrogenation of atomic N and C.

To summarise, despite the observed differences among models, the ones that better match the observations are those with higher densities and longer collapse timescales. In Sect. 3.3, we evaluate the model that best matches the ice chemical compositions observed within Ice Age.
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Time evolution of gas (solid) and ice (dashed) abundances in *MAGICKAL* during the collapse phase 1 stage. Models with $T_{\text{dust}}=14$ K and $A_v=60$ mag correspond to NIR38, while models with $T_{\text{dust}}=12$ K and $A_v=95$ mag correspond to J110621.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Ice compositions predicted by MAGICKAL for NIR38 for the models listed in Table 2. Left panel: predicted water ice column densities. Wide right panel: ice column density ratios with respect to water for the remaining species. The histograms are coloured and drawn as indicated in the upper right part of the panel. The visual extinction and dust temperatures assumed in the models are reported in the upper part of the panel. The black wide histograms report the observed Ice Age column density of water and ice column density ratios obtained by McClure et al. (2023). The $3\sigma$ error bars are also shown.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Fig. 4. Predicted ice compositions by MAGICKAL for J110621. The legend details are as in Fig. 3.</figcaption></figure>

### 3.2.2. Results from the MONACO code

The evolution of the ice chemical composition during the collapse (phase 1) stage obtained with the MONACO code is presented in Fig. B.1. The evolution of the ice molecular abundances along the collapse is similar to that predicted by the MAGICKAL code. The final ice compositions for all considered models with the MONACO code are presented in Figs. C.1 and C.2.

As can be seen from these figures, the abundances of major ice species ($H_2O$, CO, $CO_2$) behave differently in the models with $T_{dust} = 14$ K from those with $T_{dust} = 12$ K. At $T_{dust} = 14$ K, solid CO is converted into $CO_2$ ice more efficiently than at $T_{dust} = 12$ K. The major route for $CO_2$ ice formation is the reaction $CO + OH \rightarrow CO_2 + H$. At $T_{dust} = 14$ K, this reaction proceeds more efficiently both due to the higher probability to overcome the reaction activation barrier and also due to the non-zero mobility of CO molecules at this dust temperature. Interestingly, the contribution of the surface reaction $O + CO \rightarrow CO_2$ to the formation of $CO_2$ is minor in the MONACO model both at $T_{dust} = 14$ K and $T_{dust} = 12$ K. Solid $H_2O$, on the contrary, forms more efficiently on colder grains, which results in twice higher final abundance of $H_2O$ ice at 12 K. The abundance of solid $CO_2$ surpasses that of solid $H_2O$ at 14 K, despite the continuous $H_2O$ ice accumulation. At 12 K, however, the result is the opposite: $H_2O$ is the most abundant ice species, with a predicted solid abundance that is a factor of ~3 higher than those of CO and $CO_2$ ices.

From Figs. C.1 and C.2, we find that the fast-collapse models ($t_{col}=10^4$ years) exhibit the worst agreement with the Ice Age observations both in terms of absolute $H_2O$ ice column densities and ice chemical composition. The models of NIR38 significantly overproduce solid $CO_2$, yielding abundances that clearly exceed those of $H_2O$ ice, which is not observed (McClure et al. 2023). The models of J110621, with lower dust temperatures ($T_{dust} = 12$ K), produce in general more reasonable results. For the J110621 runs, CO and $CO_2$ are produced in similar amounts (i.e. 20–40% the abundance of $H_2O$ ice) in models with intermediate or slow collapse times ($t_{col}=10^5$ and $10^6$ yr) and intermediate to high final densities ($2 \times 10^5$ and $2 \times 10^6$ cm$^{-3}$).

However, for the slow-collapse models ($t_{col}=10^6$ yr), the CO/$CO_2$ ice abundance ratio reaches values ~2, which is similar to the one measured with the Ice Age data. We note that the rest of models results in CO/$CO_2$ abundance ratios <1.

Methanol ice is severely underproduced in the models with $T_{dust} = 14$ K, which can be related to the overproduction of $CO_2$ ice (Figs. C.1 and C.2). The formation of $CH_3OH$ in the MONACO model runs occurs via successive hydrogenation of CO, which is driven by diffusive H-addition reactions. No significant impact on the final $CH_3OH$ ice abundance is observed when including the chain of grain-surface reactions $CH_4 + OH$ and $C + H_2O$ proposed by Qasim et al. (2018), Molpeceres et al. (2021) and Potapov et al. (2021). The final abundance of $CH_3OH$ ice is higher in the model with $T_{dust} = 12$ K as expected from the low binding energy of atomic hydrogen. As for $CH_3OH$, $CH_4$ ice is also underproduced in all models, although the slow-collapse models tend to provide better agreement with the Ice Age observations.

In summary, models with higher final densities and collapse timescales of $10^5–10^6$ yr provide a closer match to the ice chemical composition obtained by the Ice Age observations. This conclusion is similar to those obtained for the MAGICKAL chemical code (see Sects. 3.2.1 and 3.2.3).

### 3.2.3. Results from the Nautilus code

The evolution of the ice and gas phase abundances obtained by the Nautilus code for the phase 1 stage is reported in Fig. B.2. The evolution of the ice molecular abundances is similar to that predicted by the MAGICKAL and MONACO codes as described in Appendix B.

The predicted column densities of water ice and the ice abundance ratios obtained by the Nautilus code for phase 1 for the different collapse models are shown in Figs. C.3 and C.4. As for the MAGICKAL and MONACO codes, CO is formed in the gas-phase and subsequently depletes onto the surface of dust grains. The CO ice molecular abundances increase with density and time, although not linearly. Because of the dust temperatures $T_{dust} \ge 12$ K, a large fraction of the CO that sticks on the grains is converted into either $CH_3OH$ or $CO_2$. For $CH_3OH$, the main
formation process is hydrogenation of CO. However, for CO₂, its conversion from CO does not occur via the reaction CO + OH → CO₂ + H as found for the MAGICKAL and MONACO codes, but through the grain surface reaction CO + O → CO₂. This conversion is very efficient, and the fraction of CO ice is clearly underestimated while the CO₂ and CH₃OH fractions tend to be overestimated with respect to the Ice Age observations, independently of the n_H gas density and the collapse timescales considered in the models. However, note that the ice abundances of CO increase for the highest-density, slowest collapse models, which induce a decrease in the predicted CO₂ ice abundances. As found for the MAGICKAL and MONACO codes, the predicted ice abundances get closer to the observed values for the latter models, and hence they are the ones that best reproduce the Ice Age observations.

For NH₃ and CH₄, these species are mostly formed by the successive hydrogenation of N and C on the grain surfaces. While our models do a good job at reproducing the measured ice abundances of NH₃, CH₄ is systematically underproduced in all cases except for the slow-collapse models with t_col=10⁶ yr (see Figs. C.3 and C.4).

### 3.2.4. Results from the UCLCHEM code

In Fig. B.3 we present the evolution of the relevant species both in the ice and in the gas phase for the phase 1 runs (collapse stage). The growth of the ices mainly occurs at the end of the collapse phase with the rapid increase of the density. The rate of ice growth in the UCLCHEM models is, however, higher than for the previous codes (MAGICKAL, MONACO, and Nautilus) because the initial abundances for the phase 1 runs are lower as a result of the efficient non-thermal desorption during the static (phase 0) translucent cloud stage (see Sect. 4.6). This is likely the reason for the underproduction of water ice in the UCLCHEM runs, although the predicted values for some of the models still lie within a factor of 10 with respect to the ice abundances observed with the JWST (see Figs. C.5 and C.6).

For CO, CO₂, and CH₃OH, the predicted ice abundances are either underestimated or overestimated depending on the model. Like for the other codes, the overproduction of CO₂ is detrimental for CH₃OH formation as the majority of atomic C is locked into the inert CO₂ species. According to our models, the main formation route of CH₃OH on grains at the high visual extinctions of the field stars NIR38 and J110621 is the hydrogenation of CO on dust grain surfaces. Regarding CO₂, its formation mainly occurs on the surface of dust grains via the reaction CO + OH → H + CO₂. We find no difference in the final ice abundance of CO, CO₂, and CH₃OH between the case of adding the reactions C + H₂O → H₂CO, and CH₄ + OH → CH₃ + H₂O and CH₃ + OH → CH₃OH to the chemical network and the case of not adding them.

For NH₃ and CH₄, UCLCHEM tends to overpredict the ice abundance of these molecules by factors ~2–8. It is interesting to note that UCLCHEM provides the best match to the observed ice CH₄ abundance. This is due to the fact that gas-phase atomic C is roughly one order of magnitude more abundant in the UCLCHEM runs than in the runs for the other codes (see Fig. B.3), which yields the depletion of larger amounts of atomic C onto dust grains, enhancing the production of CH₄ ice.

Overall, the models that best reproduce the observed ice compositions are those with intermediate collapse times ($t_{col} = 10^5$ yr) and intermediate or high final densities ($n_H = 10^5, 10^6$ cm⁻³) (see Sect. 3.3), although there are larger mismatches for the field star J110621 (see Figs. C.5 and C.6).

### 3.2.5. Results from the KMC model

Figure B.4 shows the evolution of the ice chemical composition during the phase 1 runs for the low-density fast-collapse models and for the high-density slow-collapse models. The time-dependent gas abundances of CO and C – modelled by UCLCHEM gas-phase simulations – are also included in Fig. B.4. The high initial CO abundance compared to C (10⁻⁴ vs. 10⁻⁶ with respect to total H nuclei) drives the CO-based chemistry on the grain. With fast collapse and low density, no significant freeze-out of CO occurs, while with the slow-collapse high-density models, depletion of CO starts at later times (>10⁵ yr) induced by the rapid increase of density. A steady increase of ice (Bulk) abundances is generally observed throughout the collapse. We note that the increase rate seems unaffected by the accelerated increase of density, except for CO in the slow-collapse models, when the above-mentioned freeze-out occurs.

The heavy fluctuations observed for CH₄ is a finite-size effect, due to the very low concentration numbers generated by the KMC on the grain. That is, the CH₄ abundance is found to fluctuate between 1 and 0 on our 50 × 50 lattice. Hence, not much value should be given to the variation of CH₄ abundance shown in this figure. We note that the 50 × 50 size surface model is chosen as a compromise between simulation time and size effect. The algorithm scales at least with N² log N, and hence much larger sizes become prohibitively expensive. Fortunately, because of the large number of species on the grain the finite size effect is rather limited and a few short test simulations show that the 100 × 100 size grains lead to the same results as the chosen 50 × 50 grain. This is for the initial stage where the number density is low and the largest effect would be expected.

The final ice abundances predicted by the KMC simulations are shown in Figs. C.7 and C.8. We note that due to the long computation times not all models of NIR38 reached the final collapse times and the results shown correspond to slightly shorter collapse times (at least 90% of the final collapse time; see caption in Fig. C.7). See Sect. 4.10 for a discussion of the limited applicability of the KMC models.

As shown in Figs. C.7 and C.8, the absolute column densities of H₂O ice tend to be underproduced for all models and towards both sources. While for the shortest collapse time models (10⁴ yr) there is a clear mismatch, for the intermediate and longer collapse times (≥10⁵ yr), the H₂O abundance is still within a factor 10 of the observed column density. In terms of ice composition, an increasing trend in abundance is observed for CO, CO₂, and CH₃OH with increasing density and collapse time, where the models with collapse times ≥10⁵ yr match better the observations with the exception of CH₄. The ice fractional abundances with respect to H₂O, however, remain fairly similar for most species for a given collapse time independently of the density considered in the model. Notable differences are seen for CO and CO₂, where we observe an increasing trend with density, although with the largest collapse time ($t_{col}=10^6$ yr) the fractional abundances level off. Hence, if enough time is given to the ices to grow, these species converge to some final fractional value with respect to H₂O, which is reached faster for higher densities.

Comparing CO₂ to CO, we find that their final column density ratio is relatively close to 1 (i.e. CO₂/CO~0.3–3.7). Both species have slightly lower fractional abundances than observed for shorter collapse times and lower densities, while the opposite is true for longer collapse times and higher densities. For methanol, some differences are found between the simulations for J110621 (Fig. C.8) and for NIR38 (Fig. C.7). For the higher dust temperature models ($T_{dust} = 14$ K), the methanol fraction
<table><caption>**Table 5.** MAE results for the models for NIR38 and J110621.</caption><thead><tr><th rowspan="2">Chemical code</th><th colspan="3">NIR38<br>(A<sub>V</sub> = 60 mag, T<sub>dust</sub>=14 K)</th><th colspan="3">J110621<br>(A<sub>V</sub> = 95 mag, T<sub>dust</sub>=12 K)</th></tr><tr><th>n<sub>H</sub> (cm<sup>-3</sup>)</th><th>t<sub>col</sub> (yr)</th><th>MAE</th><th>n<sub>H</sub> (cm<sup>-3</sup>)</th><th>t<sub>col</sub> (yr)</th><th>MAE</th></tr></thead><tbody><tr><td>MAGICKAL</td><td>2×10<sup>6</sup></td><td>10<sup>6</sup></td><td>0.43<sup>+0.11</sup><sub>-0.09</sub></td><td>2 × 10<sup>6</sup></td><td>10<sup>6</sup></td><td>0.56<sup>+0.10</sup><sub>-0.09</sub></td></tr><tr><td>MONACO</td><td>2 × 10<sup>6</sup></td><td>10<sup>6</sup></td><td>0.63<sup>+0.15</sup><sub>-0.10</sub></td><td>2 × 10<sup>6</sup></td><td>10<sup>6</sup></td><td>0.45<sup>+0.10</sup><sub>-0.08</sub></td></tr><tr><td>Nautilus</td><td>2 × 10<sup>4</sup></td><td>10<sup>6</sup></td><td>0.54<sup>+0.12</sup><sub>-0.09</sub></td><td>2 × 10<sup>6</sup></td><td>10<sup>6</sup></td><td>0.46<sup>+0.09</sup><sub>-0.09</sub></td></tr><tr><td>UCLCHEM</td><td>2 × 10<sup>5</sup></td><td>10<sup>5</sup></td><td>0.48<sup>+0.17</sup><sub>-0.14</sub></td><td>2 × 10<sup>5</sup></td><td>10<sup>5</sup></td><td>0.62<sup>+0.11</sup><sub>-0.12</sub></td></tr><tr><td>KMC</td><td>2 × 10<sup>5</sup></td><td>10<sup>5</sup></td><td>0.80<sup>+0.14</sup><sub>-0.11</sub></td><td>2 × 10<sup>5</sup></td><td>10<sup>5</sup></td><td>0.76<sup>+0.09</sup><sub>-0.08</sub></td></tr></tbody></table>

**Notes.** The errors associated with the MAE results correspond to the 3σ level.

is well reproduced, except for the longest collapse times for which it is overproduced. For the lower dust temperature simulations ($T_{dust}$ = 12 K), methanol is always slightly overproduced, regardless of collapse times and densities. This might be due to the limited network that misses some reactions leading to more complex molecules such as ethanol. We find that our main formation route for methanol is through successive hydrogenation of CO since there is not enough atomic carbon or CH$_4$ to lead to efficient formation of CH$_3$OH from these reactants (Sect. 4.2). Therefore, the C + H$_2$O → H$_2$CO pathway towards CH$_3$OH plays only a minor role.

Unlike the other main ice species, the relative abundances of NH$_3$ show rather homogeneous results for the different physical conditions with only a minor dependence on dust temperature, collapse time, and density. For CH$_4$, however, this species shows larger variations. It is heavily underproduced when compared to the observed value for both sources and it does not obey the upward trend in fractional abundance with collapse time. We find that the formation of CH$_4$ is hindered by the lack of available atomic carbon on the grain, since elemental carbon is mainly locked up in CO in the gas phase (see Sect. 3.1). As reported in Fig. B.4, the gas-phase abundances of CO and atomic C show CO/C ≫ 1 for all models during the entire duration of collapse. We note that this ratio is already rather high at the start of collapse (~50), as CO is efficiently produced in the gas during the phase 0 translucent cloud stage.

### 3.3. Best-fit models

To identify the model that best matches the Ice Age observations, we calculate the median absolute error (MAE) between the model predictions and the observed column densities in two steps as described in Appendix D. This method is also called mean distance of disagreement, as introduced in Wakelam et al. (2006). In the initial step, we identify a sub-sample of best models characterised by H$_2$O column densities consistent with the Ice Age observations within a factor of a few. The models that are favoured are those with longer collapse timescales and higher densities because the efficiency of H$_2$O ice accumulation is enhanced and photo-dissociation, a primary mechanism for H$_2$O ice destruction in low-density environments, is mitigated. In a second step, we obtain the MAE between the predicted and the observed ice abundance ratios for the rest of species with respect to H$_2$O and on logarithmic scale, in order to select the model (or models) with the lowest MAE value (see Appendix D). Table 5 summarises the results of this analysis, and Fig. 5 presents the comparison of the predicted ice column densities between all codes and with the Ice Age observations.

For MAGICKAL, the code generally shows best agreement with observations for higher densities and longer collapse timescales. For NIR38 ($A_V$=60 mag, $T_{dust}$=14 K), the most suitable model conditions are $n_H$=2 × 10$^6$ cm$^{-3}$, $t_{col}$ = 10$^6$ yr, with logarithmic MAE value for the ice abundances with respect to H$_2$O of 0.43$^{+0.11}_{-0.09}$. For J110621 ($A_V$ = 95 mag, $T_{dust}$ = 12K), the observed ice abundance ratios are best explained by the same model conditions, with a corresponding MAE value of 0.56$^{+0.10}_{-0.09}$.

Similarly, for MONACO, the best-fitting models for both NIR38 and J110621 are the ones with $n_H$ = 2 × 10$^6$ cm$^{-3}$ and $t_{col}$ = 10$^6$ yr. The MAE values are 0.63$^{+0.15}_{-0.10}$ and 0.45$^{+0.10}_{-0.08}$, respectively.

For Nautilus, the best-fitting models for NIR38 and J110621 are respectively ($n_H$ = 2 × 10$^4$ cm$^{-3}$, $t_{col}$=10$^6$ yr) and ($n_H$ = 2 × 10$^6$ cm$^{-3}$, $t_{col}$=10$^6$ yr). For NIR38, the model underproduces the H$_2$O column density by a factor of ~1.5–2, with an associated MAE of 0.54$^{+0.12}_{-0.09}$. The second best model is the one with the highest density and longest collapse phase ($n_H$ = 2 × 10$^6$ cm$^{-3}$, $t_{col}$=10$^6$ yr) with a MAE of 0.95$^{+0.12}_{-0.10}$. Although the latter gives in general closer values to the observed H$_2$O column density and ice abundance ratios of CO, CO$_2$, CH$_3$OH, and NH$_3$, it fails to produce any CH$_4$ and this is the reason for obtaining a worse MAE for this model. The observations towards J110621 are better reproduced by the models with an overproduction of H$_2$O only by a factor of ~1.2–1.4 between the modelled and the observed value. The associated MAE is 0.46$^{+0.09}_{-0.09}$. Slower collapse overall reproduces more closely the observed column densities.

For UCLCHEM, the two best models are the ones with intermediate collapse times ($t_{col}$ = 10$^5$ yr) and intermediate or high final densities ($n_H$ = 2 × 10$^5$, 2 × 10$^6$ cm$^{-3}$). Between these two models, we cannot favour one over the other, as both the errors in column density and the logarithmic MAE for the abundance ratios with respect to H$_2$O show similar values and range between 0.48$^{+0.17}_{-0.14}$ and 0.68$^{+0.12}_{-0.12}$.

Similarly to UCLCHEM, the intermediate physical conditions yield the best results for the KMC code. That is, models with collapse time $t_{col}$=10$^5$ yr and final density $n_H$ = 2 × 10$^5$ cm$^{-3}$ are the ones that best fit the observed ice abundances. These models underproduce the H$_2$O column density by a factor of ~5 and 8 and show a MAE of 0.8$^{+0.14}_{-0.11}$ and 0.76$^{+0.09}_{-0.08}$ for the fractional abundances of NIR38 and J110621, respectively. Other models with lower MAE (~0.6–0.7) are found for shorter collapse times, but they are off in terms of absolute H$_2$O column densities with an underproduction by factors ~20–70.

From Fig. 5, we find that the largest discrepancies between the best-fit models with the observed ice column densities are found for CH$_3$OH and CH$_4$. These species tend to be
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Comparison between the ice column densities predicted by all codes for the best-fit models, and with the Ice Age observed values. The different colours indicate different molecular species and the different symbols label the distinct astrochemical codes used in this work (see legend at top of figure). The horizontal lines show the observed ice column densities reported in Table 4 and the shaded areas indicate their associated errors.</figcaption></figure>

underproduced in the ice, especially for the models of NIR38. This is related to the majority of atomic C getting locked into CO in the gas-phase early in the simulations (Sect. 3.1), which is subsequently converted into CO₂ ice upon accretion, preventing the formation of CH₃OH and CH₄ in the ice.

To summarise, the models that better reproduce the ice chemical composition observed towards NIR38 and J110621 are those with intermediate–high densities and intermediate–slow collapse times; models with fast collapse and low densities are clearly discarded. These results nicely match the densities and timescales required by grain growth models ($n_H \ge 10^5$ cm⁻³ and $t_{col} \ge 10^6$ yr) to obtain molecular cloud grains of 0.9 µm in size, which is the size inferred from the JWST Ice Age observations (see Dartois et al. 2024).

### 3.4. Models with $T_{dust}=10$ K

In the model runs explored above, we have considered $T_{dust}=12$ K and 14 K as estimated from Herschel-based dust temperature maps (see Sect. 2.6). However, it is well known that these maps are representative of the outer and slightly warmer layers of molecular clouds (see e.g. Hocuk et al. 2017; Clément et al. 2023). To test the influence of the dust temperature on the ice chemical composition towards the field stars NIR38 and J110621, we have repeated the simulations for the best models identified in Table 5, but assuming a dust temperature of $T_{dust} = 10$ K. The results differ significantly for the codes that include only diffusive chemistry (Nautilus and UCLCHEM) from those that include both diffusive and non-diffusive chemistry (MAGICKAL, MONACO, and KMC).

For Nautilus, the main impact on the results comes from the severe decrease in the production of CO₂ by factors of 55 and 65 for NIR38 and J110621, respectively. At the same time, CO is more abundant in the 10 K models (by factors of 3.6 and 1.8, respectively) as a consequence of a slower diffusion of both CO

and O on the surface of dust grains<fnref n="7" />. CH₃OH production is thus increased in both models by factors of 2 and 25 for NIR38 and J110621, while CH₄ is also increased in the $T = 10$ K model, for which CH₄ production is in competition with methanol. Similarly to Nautilus, the decrease of $T_{dust}$ in the UCLCHEM simulations induces a drastic reduction in the formation of CO₂ as well as an increase in the formation of CH₃OH. This behaviour was hinted from the results for J110621 ($A_V = 95$, $T = 12$ K), for which we already observed an underproduction of CO₂ due to the decrease in the efficiency of the reaction CO + OH → CO₂ + H.

In contrast to UCLCHEM and Nautilus, we do not observe any significant differences when the temperature decreases down to 10 K for the MAGICKAL, MONACO, and KMC codes. For the MAGICKAL code, the chemistry hardly changes for dust temperatures <12 K because non-diffusive chemistry takes over diffusive chemistry at these low dust temperatures. As an example, CO₂ predominantly forms at $T_{dust}=14$ K via diffusive association of CO and OH. However, for lower $T_{dust}$, the reaction rate of diffusive association of CO and OH declines, and non-diffusive mechanisms begin to dominate the formation of this molecular species. The formation rate of both mechanisms becomes equivalent at ~13 K. Thus, the decrease in CO₂ (by a factor of 2) is not as drastic as observed for UCLCHEM or Nautilus, because CO₂ can still form through non-diffusive mechanisms at $T_{dust} < 12$ K.

In the case of CH₃OH, a moderate increase in CH₃OH ice abundance is predicted (by a factor of 2), which is related to the increase of CO ice abundance as a result of its slower conversion into CO₂ at the lower dust temperature. As shown by the KMC simulations, CH₃OH is predominantly formed through CH₃O + H → CH₃OH, but CH₃O + H₂CO → CH₃OH + HCO also plays a role (Simons et al. 2020; Santos et al. 2022). Therefore, the dominant mechanisms responsible for ice formation at $T_{dust}=10$ K for the MAGICKAL, MONACO, and KMC codes are non-diffusive, which can operate with the same efficiency independently of the dust temperature considered in the model.

## 4. Discussion

### 4.1. Depletion of CO onto ice grains

From the Ice Age observations of McClure et al. (2023), CO is the second most abundant species in the ices along the lines of sight of NIR38 and J110621. Despite the significant abundance of CO in solid state (~30–40% with respect to H₂O; see Table 4), it has been proposed that the cloud local gas density towards these field stars may not exceed the limit of $n_H \sim 10^5$ cm⁻³ required for the catastrophic depletion of CO (McClure et al. 2023). This idea is supported by the similar ice chemical composition of the NIR38/J110621 field stars to Elias 16 (McClure et al. 2023), whose visual extinction is $A_V \sim 19$ mag.

By using the results from the best models reported in Table 5, we can provide information about the level of CO depletion predicted for the dense cloud regions towards field stars NIR38 and

<footnote n="7">⁷ The diffusion rate of each species is proportional to $\exp(-E_{bind}/kT_{dust})$, where $E_{bind}$ is the binding energy and $k$ is the Boltzmann constant. This implies that the diffusion rate of atomic O, for instance, decreases by 11 orders of magnitude at 10 K with respect to the case at 12 K.</footnote>
<table><caption>**Table 6.** Predicted CO depletion factors for the best models.</caption><thead><tr><th rowspan="2">Code</th><th colspan="2">NIR38 (A<sub>V</sub> = 60 mag)</th><th colspan="2">J110621 (A<sub>V</sub> = 95 mag)</th></tr><tr><th>χ(CO)<sub>gas</sub> (10<sup>-4</sup>)</th><th>f<sub>D</sub><sup>(a)</sup></th><th>χ(CO)<sub>gas</sub> (10<sup>-4</sup>)</th><th>f<sub>D</sub><sup>(a)</sup></th></tr></thead><tbody><tr><td>MAGICKAL</td><td>0.04</td><td>25</td><td>0.03</td><td>30</td></tr><tr><td>MONACO</td><td>0.0007</td><td>1360</td><td>0.001</td><td>790</td></tr><tr><td>Nautilus</td><td>0.3</td><td>~3</td><td>0.0007</td><td>≥1400</td></tr><tr><td>UCLCHEM</td><td>0.52</td><td>~2</td><td>0.53</td><td>~2</td></tr><tr><td>KMC</td><td>0.5</td><td>~2</td><td>0.6</td><td>1.7</td></tr></tbody></table>

**Notes.** (a) The depletion factor, $f_D$, is calculated as $10^{-4}/\chi(\text{CO})_{\text{gas}}$, with $\chi(\text{CO})_{\text{gas}}$ the abundance of CO in the gas phase predicted by the best model from each code (see Table 5).

J110621. The depletion factor of CO, $f_D$, is calculated as

$$f_D = \frac{10^{-4}}{\chi(\text{CO})_{\text{gas}}}, \qquad (2)$$

where $\chi(\text{CO})_{\text{gas}}$ is the gas-phase CO abundance predicted by our best models, and where $10^{-4}$ is the assumed canonical value for the gas-phase abundance of CO, as measured towards nearby star-forming regions (see e.g. Frerking et al. 1982; Fuente et al. 2019), and that should be considered as a lower limit.

In Table 6, we report the values of $\chi(\text{CO})_{\text{gas}}$ at the end of the collapse and the CO depletion factors predicted by each best model. From Table 6, we find that for UCLCHEM and KMC codes, the CO depletion factors in our best models range between ~2–3. This is produced by the intermediate timescales for the collapse of their best models, and it would be consistent with the idea of CO not being severely depleted towards the lines of sight of NIR38 and J110621 (McClure et al. 2023). Something similar is found for the Nautilus results of NIR38. However, for J110621, Nautilus gives a very large depletion factor ($f_D \ge 1400$), which is due to the catastrophic depletion of CO in its highest-density, slowest-collapse model. The same applies to the MAGICKAL and MONACO codes, for which their predicted depletion factors are clearly $\ge 10$. Maps of CO ice recently obtained with the JWST towards the Chamaeleon I molecular cloud estimate that between 45%–85% of CO is depleted onto dust grains, although large uncertainties do exist (Smith et al. 2025). Large-scale maps of gas-phase CO and of its $^{13}$C and $^{18}$O isotopologues towards the Chamaeleon I cloud, will be analysed in the future to accurately measure the level of CO depletion towards the lines of sight of NIR38 and J110621.

## 4.2. Formation of CH<sub>3</sub>OH ice

In all the codes used in this modelling effort, we have considered the same formation routes of methanol on the surface of dust grains. These formation routes are: 1) successive CO hydrogenation in the reaction chain CO $\rightarrow$ HCO $\rightarrow$ H<sub>2</sub>CO $\rightarrow$ CH<sub>3</sub>O $\rightarrow$ CH<sub>3</sub>OH (Watanabe & Kouchi 2002); 2) the reaction chain CH<sub>4</sub> + OH $\rightarrow$ CH<sub>3</sub> + H<sub>2</sub>O and CH<sub>3</sub> + OH $\rightarrow$ CH<sub>3</sub>OH (Lamberts et al. 2017; Qasim et al. 2018; Lamberts et al. 2022); and 3) the formation of formaldehyde through C + H<sub>2</sub>O and subsequent hydrogenation processes (Potapov et al. 2021; Molpeceres et al. 2021).

For MAGICKAL, CH<sub>3</sub>OH is predominantly (> 99.99%) formed via the successive hydrogenation of CO. The formation of methanol via CH<sub>3</sub> + OH is negligible, given that the hopping rate of heavy species such as radicals is too slow at the low $T_{\text{dust}}$ of the simulations. The diffusion of heavy radicals is very slow and the chemistry is dominated by atomic hydrogen at this temperature regime (<20 K; Garrod et al. 2008). For $T_{\text{dust}}=10–16$ K, all hydrogen atoms landing on the grain will quickly be

involved in a hydrogenation or abstraction reaction and will not desorb, because the surface concentration of reactants is large and the hydrogen atom will only need to diffuse a short distance to react. The same applies to the Nautilus and MONACO codes, for which CH<sub>3</sub>OH is only (>99.99%) formed by hydrogenation of CO on grains. The methanol formation from CH<sub>3</sub> + OH $\rightarrow$ CH<sub>3</sub>OH is insignificant, regardless the temperature, density and final collapse time.

For the KMC simulations, we also find that methanol is predominantly formed through the successive hydrogenation of CO. Around 30% of the total methanol present in the ice has been formed through hydrogenation of CH<sub>3</sub>O whereas the radical-neutral reaction CH<sub>3</sub>O + H<sub>2</sub>CO $\rightarrow$ CH<sub>3</sub>OH + HCO accounts for ~10% only. The remainder 60% comes predominantly from hydrogenation of CH<sub>2</sub>OH, which should not be considered as a formation reaction since CH<sub>2</sub>OH is exclusively formed through H-abstraction of methanol. The reaction CH<sub>3</sub> + OH $\rightarrow$ CH<sub>3</sub>OH hardly occurs in the simulations of both sources (only once in $10^5$ yr) and the contribution of the C + H<sub>2</sub>O reaction to the formation of formaldehyde is very marginal: only ~0.2% of the reaction flux towards H<sub>2</sub>CO is found to be of this type.

One may think that the low efficiency in the formation of CH<sub>3</sub>OH ice through the reaction chain starting from CH<sub>4</sub> + OH, is due to the underproduction of CH<sub>4</sub> in the models (see Sect. 3). However, we note that even when CH<sub>4</sub> is efficiently formed, as in the case of the UCLCHEM code, the CH<sub>3</sub>OH ice formation reaction CH<sub>3</sub> + OH $\rightarrow$ CH<sub>3</sub>OH contributes $\le 1.4\%$ to the total CH<sub>3</sub>OH ice produced for the best models of NIR38 and J110621.

Therefore, all codes show that the most efficient formation mechanism of methanol for the high extinction conditions probed towards the lines of sight of NIR38 and J110621, is the hydrogenation of CO on grain surfaces. We note, however, that we do not exclude that the other two mechanisms of methanol formation (namely, the chain of reactions starting from CH<sub>4</sub> + OH and C + H<sub>2</sub>O) are efficient at lower extinctions and densities, but that the dominant mechanism towards the densest regions in the Chamaeleon I cloud is CO hydrogenation. The JWST observations reveal that a CH<sub>3</sub>OH:CO mixture alone is not sufficient to reproduce the observed 9.8 µm band (McClure et al. 2023), which suggests that some mixing with water may be present towards the outer shells of the cloud (see also Sect. 4.8).

## 4.3. Conversion of CO into CO<sub>2</sub>

Multiple laboratory and theoretical works have been carried out to understand the conversion of CO into CO<sub>2</sub> in ices. The most studied surface reactions are: (i) CO + O $\rightarrow$ CO<sub>2</sub> (Goumans & Andersson 2010; Raut & Baragiola 2011; Ioppolo et al. 2013; Minissale et al. 2013); (ii) HCO + O $\rightarrow$ CO<sub>2</sub> + H (Ruffle & Herbst 2001; Ioppolo et al. 2013); and (iii) CO + OH $\rightarrow$ CO<sub>2</sub> + H (Oba et al. 2010; Ioppolo et al. 2011; Noble et al. 2011; Molpeceres et al. 2023). While reactions (i) and (iii) have been found to form CO<sub>2</sub> in different quantities, reaction (ii) has been proven to be inefficient (Ioppolo et al. 2013), and hence it is not included in the astrochemical codes. Therefore, the main routes for CO<sub>2</sub> ice formation are reactions (i) and (iii).

A general result for all models that use directly the values of $T_{\text{dust}}$ inferred from Herschel data (14 K and 12 K for NIR38 and J110621, respectively; see Sect. 3), is that CO<sub>2</sub> is overproduced compared to the CO<sub>2</sub> ice abundances measured with the JWST. We note that the abundances in McClure et al. (2023) could be a factor of ~1.45 higher, due to uncertainties in the determination
of the CO₂ 4.7 µm band strength<fnref n="8" />. The overproduction of CO₂ is a well-known effect that was pointed out by e.g. Garrod & Pauly (2011) and, more recently, by Clément et al. (2023). The formation of CO₂ ice in the ices is strongly dependent on T_dust and it is extremely efficient once T_dust≥12 K. However, note that the reaction mechanisms responsible for this CO₂ ice overproduction is different in different codes.

For the MAGICKAL code, the dominant mechanism for CO₂ ice formation is the reaction CO + OH → CO₂ + H when T_dust > 12 K. At higher T_dust, CO becomes mobile enough to reach and react with OH radicals before H or H₂, efficiently converting CO into CO₂ (see Sect. 3.2.1). This reaction is also the dominant CO₂ ice formation process in the MONACO and UCLCHEM chemical codes. In contrast, for the Nautilus code, CO₂ ice is produced via the grain-surface reaction CO + O → CO₂, after CO and atomic O being adsorbed onto the icy mantles from the gas phase. It has been proposed that this difference could be due to distinct binding energies of atomic oxygen used in the codes (Clément et al. 2023). However, in our case, the binding energy used in all codes is ~1600–1660 K, which is consistent with that recently revisited by He et al. (2015) and Minissale et al. (2022), and with the lack of detection of O₂ in the Orion Bar (Melnick et al. 2012). Alternatively, the model discrepancies could be produced by differences in the activation energy barrier assumed for the CO + O → CO₂ reaction. For example, while MAGICKAL uses 1000 K and 1.27 Å for the height and width of the reaction barrier (chosen to reproduce the tunnelling rate in the quantum chemical calculations of Goumans & Andersson 2010), Nautilus uses 627 K and 1 Å, respectively, following the experimental results of Minissale et al. (2013). This would explain why this reaction is more efficient in Nautilus than in MAGICKAL. Given all uncertainties, it thus remains unclear which of these two reactions (i.e. CO + OH → CO₂ + H or CO + O → CO₂) is responsible for the main production of CO₂ ice in the dense molecular regions towards the NIR38 and J110621 field stars.

To investigate this problem in detail, we used the KMC model where we have included all these different reactions for CO₂ ice formation. First, we considered the reactions where i) CO₂ can form directly from CO + OH — CO₂ + H; and ii) CO₂ can form through a minor channel in a two-step fashion CO + OH — HOCO and HOCO + H — CO₂ + H₂. The model adopts a branching ratio of 0.8 and 0.2 between the two routes. Quantum chemical calculations and laboratory experiments (Arasa et al. 2013; Molpeceres et al. 2023; Ishibashi et al. 2024) suggest that the latter route should be more important, since HOCO is found to stabilise in the ice rather than to fall apart into CO₂ and H. This latter process is endothermic.

Additional simulations were performed with an updated network that included CO + O — CO₂ (k = 10⁻¹⁰ s⁻¹; Goumans & Andersson 2010), HCO + O — CO₂ + H (no barrier; Westenberg & DeHaas 1972; Campbell & Handy 1978) and only the indirect HOCO route with a rate constant of 10¹⁰ s⁻¹ (Molpeceres et al. 2023). Interestingly, our results show that all routes contributed roughly equally to the formation of CO₂, while the rate constants dictate that CO reacting with O should be less favourable by 20 orders of magnitude. The most probable explanation for this is the higher mobility of atomic O on the grain. The species diffuse rapidly in the shallow sites until they reach high binding sites that can act as a sink where species

reside for a long time, long enough to react even for low reaction rate constants. This mechanism is more effective for O than for OH, as a result of its lower diffusion barrier distribution. Therefore, the chemical modelling of CO₂ on grain surfaces should consider the two formation mechanisms, CO + O — CO₂ and the indirect HOCO route.

### 4.4. Formation of hydrides in the ice: H₂O, NH₃ and CH₄

As shown in Sect. 3, the models tend to underproduce water for all codes. The models that generate the largest amounts of water are the ones with a slow-collapse phase and with high densities, reproducing the observed water ice abundances within a factor of 10.

Similarly, CH₄ also tends to present lower abundances in all model runs compared to those measured by the Ice Age observations. An exception is the UCLCHEM code, for which the models reproduce within factors of 2–4 the CH₄ ice abundances observed towards NIR38 and J110621 (see Fig. 5). The formation of CH₄ ice is intimately linked to the formation of CO₂ ice since both molecular species compete for the conversion of atomic C depleted onto dust grains. This conversion, however, is expected to take place at lower extinctions during the translucent cloud (phase 0) stage (see e.g. Garrod & Pauly 2011). As explained in Sect. 3.1, the amount of gas-phase atomic C available for depletion in the translucent cloud (phase 0) stage is roughly one order of magnitude higher in the UCLCHEM models than for the other codes. For NH₃, its formation in the ice occurs predominantly via the hydrogenation of atomic N. The models reproduce fairly well the ice abundances observed with the JWST, which differ by typically factors of a few.

We finally note that lowering the temperature from 12 K/14 K to 10 K in our simulations has, in general, little to no effect on the resulting ice abundances of H₂O, CH₄, and NH₃. This is particularly true for the codes that include non-diffusive chemistry such as MAGICKAL, MONACO, and the KMC codes.

### 4.5. Oxygen budget in the models

For our modelling, we have considered an initial abundance of 1.76×10⁻⁴ for atomic oxygen, as proposed in the EA1 low-metal abundance case of Wakelam & Herbst (2008, see Table 3). We can check whether this abundance is consistent with the oxygen budget locked into H₂O ice, COice, and CO₂ice in the models. Taking into account that our best models provide values close to those measured with the JWST, we use the abundances shown in Table 4 to make this estimate. For NIR38, the oxygen elemental abundance locked into the ice can be calculated from the sum N(H₂O ice) + N(CO ice) + 2×N(CO₂ice), which gives ~1.3×10¹⁹ cm⁻². Since the column density of hydrogen nuclei towards NIR38 is ~10²³ cm⁻², this yields an oxygen abundance of 1.34×10⁻⁴, which corresponds to 76% of the initial oxygen assumed in the models. However, the EA1 abundance of 1.76×10⁻⁴ considers depletion onto dust grains, and hence when compared to the cosmic abundance of oxygen (~4.9×10⁻⁴; see Asplund et al. 2009), the oxygen abundance locked into ices represents 27% of the total atomic oxygen. On the other hand, the derived oxygen abundance in the ice towards the J110621 line of sight is ~2.1×10⁻⁴, i.e. higher than the initial atomic oxygen abundance assumed in our models. Therefore, the EA2 and EA3 cases proposed by Wakelam & Herbst (2008), provide a more realistic initial abundance for elemental oxygen for astrochemical modelling of dense molecular clouds.

<footnote n="8">⁸ The band strength used in McClure et al. (2023) was 1.1×10⁻¹⁶ cm molec⁻¹ as revisited by Bouilloud et al. (2015), but a value of 7.6×10⁻¹⁷ cm molec⁻¹ can be found in earlier works (Gerakines et al. 1995).</footnote>
## 4.6. Importance of the length of the collapse

In our models, we have considered different collapse time lengths to establish the impact of this parameter on the ice chemical composition. As shown in Sect. 3, the abundances of major ice species such as H₂O, CO, or CO₂ change little in the fast-collapse models because the collapse time length ($t_{col}=10^4$ yr) is orders of magnitude shorter than the depletion timescale ($t_{freeze-out} \approx 10^9/n_H$ yr). Only for the highest-density models with $n(H)=10^6$ cm⁻³, for which the depletion timescales are similar to the collapse time length, the predicted ice abundances reach values closer to those measured in the Ice Age observations. As discussed in Sect. 3.3, the best agreement with the observations is achieved for models with longer collapse time lengths (i.e. with timescales $t_{col}=10^5$ yr and $t_{col}=10^6$ yr; see Table 5), although note that shorter collapse times favour the formation of species such as CH₄ (see Figs. C.8 and C.7). The high densities and long timescales for ice formation are broadly consistent with the $n_H=10^5$ cm⁻³ and $t=10^6$ years required by grain growth models for molecular clouds (Ormel et al. 2009) in order to reproduce the 0.9 µm grain sizes found for NIR38 and J110621 by Dartois et al. (2024). The convergence of models for two different physical and chemical phenomena adds weight to these results.

## 4.7. Uncertainties in the determination of the extinction towards NIR38 and J110621

The visual extinction values first reported in McClure et al. (2023) towards the lines of sight of NIR38 and J110621 ($A_V=60$ mag and $A_V=95$ mag, respectively), were pre-flight estimates. A more recent analysis of the NIRSpec spectra obtained with the JWST towards these two field stars, has however yielded $A_K$ values of 4.1 for NIR38 and 5.6 mag for J110621 (Dartois et al. 2024). If we assume $A_V/A_K = 8.4$, these values correspond to visual extinctions of $A_V=34$ mag and $A_V=47$ mag, respectively. These values are different from those given by McClure et al. (2023), because these authors used the pre-launch extinctions determined photometrically under the assumption that both background stars were both G or K giants. Instead, Dartois et al. (2024) fitted the photospheric absorption features to determine that both stars were K7 dwarfs.

The new $A_V$ values are therefore factors of ~2 lower than those reported in McClure et al. (2023) and than those used in our chemical modelling simulations (see Sect. 2 and Table 2), which explains the similar ice compositions of NIR38/J110621 and Elias 16. However, we note that we do not expect significant differences for these lower visual extinctions because UV photo-chemistry (the only one that could be affected by the lower extinction values) is inefficient at $A_V>10$ mag. The only impact is expected for the predicted column densities of water ice reported in Figs. 3–C.8, whose values would be lower by roughly a factor of ~2. The ice column density ratios would not suffer any change because they are calculated directly from the ice abundances predicted by the models. The predicted CO depletion factors would not be affected either by the change in visual extinction for NIR38 and J110621, because $f_D$ are inferred from the CO gas-phase abundance given by the models (see Sect. 4.1).

## 4.8. Other possible effects not considered in the models: Grain growth, cloud dynamics and line-of-sight effects

The aggregation and clustering of icy dust grains during the collapse of a cloud will modify their size distribution, as explored numerically in e.g. Weingartner & Draine (2001);

Ormel et al. (2011); Silsbee et al. (2020); Marchand et al. (2023); Lebreuilly et al. (2023). For the two background sources behind Chamaeleon I studied here, NIR38 and J110621, it has been shown by simulating the observed absorption spectra that the largest grains in the distribution have increased in radius by a factor of about four (Dartois et al. 2024) with respect to the largest grains in the expected diffuse ISM distribution (Mathis et al. 1977), giving rise to clear scattering features. Enhanced scattering observed more widely across this cloud with Spitzer photometric filters has also been attributed to grain growth (Lefèvre et al. 2014). Grain growth is accompanied by a lowering of the visual extinction per dust mass (see Sect. 4.7) and is expected to lower the peak dust temperatures and modify the radiative and sublimation-induced cooling efficiencies of dust grains (Herbst & Cuppen 2006). In addition, the surface to volume ratio might be affected, and thus could influence, possibly diminishing, the available total surface for chemistry and the surface available for depletion from the gas phase. Clues as to the local ice structure and morphology were given by the observation of the “dangling OH” absorption features of H₂O ice along both lines of sight (Noble et al. 2024). These features, while not providing a direct measure of porosity/compaction of the ice, trace H₂O molecules not fully bound to surrounding H₂O molecules in the ice, whether due to interaction with other molecules such as CO, CO₂, or due to a non-bound OH in the bulk or at the ice surface. The surface to volume parameter will depend on the degree of compaction of the ice during the aggregation process. Modification of the grain surface morphology due to energy dissipation following grain collision or cosmic ray hits can also impact the chemistry that occurs in and on the grain, as the availability and nature of binding sites changes. All of these physical modifications to the grain will, in turn, modify binding energy distributions at the surface and in the bulk, and thus activation barriers to diffusion and adsorption/desorption (Noble et al. 2012; Cuppen et al. 2013; Karssemeijer & Cuppen 2014; Nguyen et al. 2018; Minissale et al. 2022), impacting the availability of reactants for, and the relative rates of, the various chemical processes discussed above (Herbst & Cuppen 2006; Penteado et al. 2017; Iqbal & Wakelam 2018).

Another parameter that could affect the results of our simulations is the history of the evolution of the molecular cloud. In the simulations presented in this paper, we have considered a simple evolution of the physical properties of the cloud (temperature and density) during its collapse phase (we have assumed free-fall collapse; see Sect. 2.6 and Brown et al. 1988; Rawlings et al. 1992). However, by doing so, we have only considered a small portion of the possible evolutionary pathways that a parcel of gas in the Chamaeleon I cloud could undergo. By using hydrodynamical simulations, Clément et al. (2023) have shown that the passage from the diffuse to dense cloud phase in these simulations, was key to producing CO₂ very early in the cloud’s evolution when the dust temperature is larger than 12 K, but also to generate enough hydrides before atoms are fully converted into other molecules in the gas-phase. In addition, we note that, for the sake of simplicity, we have kept the dust and gas temperature equal in our simulations, and have varied them only between 14 and 10 K. As explored by Clément et al. (2023), the decoupling between the gas and dust temperatures may lead to different ice chemical compositions.

Finally, the model results compared to the JWST observations are those at the end of the phase 1 collapse stage, which is characterised by constant physical conditions. In reality, JWST observations towards the two field stars NIR38 and J110621 are probing all material along the two lines of sight from the outer
edges to the central zone of the molecular cloud. Therefore, line-of-sight effects due to the inhomogeneous structure of the molecular cloud are not taken into account by the modelling. This could be one of the reasons behind the general difficulty of the models to reproduce the observed CH₃OH and CH₄ ice composition simultaneously, as these two observed ice components may probe different locations along the line of sight. For example, CH₄-rich ice could be present in the outer shell of the cloud, where the C to CO conversion is not complete (so that C atoms can be efficiently hydrogenated once accreted onto dust grains). In contrast, CH₃OH ice should copiously form towards the central colder and denser regions of the cloud, where most atomic C is locked into CO and CO can efficiently freeze out onto the surface of dust grains to be transformed into CH₃OH via successive hydrogenation. However, note that visual extinctions A<sub>V</sub>≤2 mag are required for carbon to remain in its atomic form in the gas phase. As shown in Fig. 1, even in the translucent cloud regime with visual extinctions A<sub>V</sub>=2 mag, atomic C is almost completely converted into gas-phase CO. This implies that the region where atomic C could avoid conversion into CO extends at most A<sub>V</sub>=2 mag into the cloud from the surface. Considering that every cloud has two surfaces – the one on the observer's side and the second one on the far side – the surface layer would extend for 4 mag in total, which represents ≤12% of the NIR38 line of sight, assuming that its extinction is at least A<sub>v</sub>=34 mag (as mentioned in Sect. 4.7). Detailed physical and chemical modelling of the cloud is needed to reconstruct such line-of-sight effects, and this requires a dedicated observational effort, including dust continuum and spectral line emission mapping, which is undergoing.

## 4.9. Differences between diffusive versus non-diffusive chemistry

In this work, we have used astrochemical codes that consider diffusive chemistry only (e.g. Nautilus, UCLCHEM) and codes that include both diffusive and non-diffusive chemistry (MAGICKAL, MONACO, KMC). The largest differences between these codes are mainly found for models with very low dust temperatures of T<sub>dust</sub>=10 K. In codes with diffusive-chemistry-only, diffusion of grain-surface species becomes activated at higher T<sub>dust</sub>. This is the reason why the production of species such as CO<sub>2</sub> is largely decreased, in favour of CH<sub>3</sub>OH and CH<sub>4</sub>. In contrast, for codes that include diffusive and non-diffusive chemistry, the formation of new molecular species at T<sub>dust</sub>=10 K occurs exclusively by non-diffusive processes (i.e. in-situ on the grain surface upon the formation of the reactants), and therefore it does not require the species to diffuse across the surface. This implies that non-diffusive chemistry is efficient at T<sub>dust</sub>=10 K. As shown in Sect. 3.4, no differences are found between models with T<sub>dust</sub>=10, 12 and 14 K for codes with non-diffusive chemistry, which indicates that the formation of H<sub>2</sub>O, CO, CO<sub>2</sub>, CH<sub>3</sub>OH, CH<sub>4</sub> and NH<sub>3</sub> in these models occurs mainly non-diffusively at these T<sub>dust</sub>. This could explain the rather constant ice chemical composition observed not only in Chamaeleon I but also in other dense molecular clouds (as e.g. towards the Elias 16 background star; Knez et al. 2005), as found by McClure et al. (2023).

## 4.10. Differences between rate equation codes versus stochastic KMC models

Rate equation codes are macroscopic models in which every different chemical reaction results in a term in the rate equation that

combines the meeting of the reactants and the subsequent reaction. Whether the reactants meet by diffusion or through some non-diffusive mechanism should be explicitly included in the rate equations, thereby requiring an estimation of the coupling between the different processes and the average outcome.

In contrast, the lattice KMC model is microscopic in the sense that it simulates every move of each individual species on top of the grain surface, monitoring their exact location. In this model, the effect of diffusion on reactions is observed as an outcome of the simulation of each individual move of a species. This can pose an advantage, especially when the coupling between processes is difficult to model or when assumptions for the rate equations do not hold throughout the entirety of evolution of the system. Another advantage is that the surface structure of the grain can be taken into account. The rates of various processes (e.g. diffusion, reaction and desorption) is made site-specific and so the effect of surface morphology on the surface chemistry can be modelled in detail.

On the other hand, the microscopic approach of the lattice KMC grain model comes with a cost in terms of practical applicability. Certainly, when compared to macroscopic rate equation models, the computation time is very long. Even so much so, that for the highest final densities and longest collapse times of the grid studied in this work, the computation time becomes inapplicable (estimated time is on the order of months). Therefore, for some combinations of parameters, we have not been able to include results for the entire collapse time duration (instead at least 94% of collapse time). Moreover, as a single KMC run simulates one possible (stochastic) evolution of the system (grain), ideally one would like to perform multiple runs and average over the results. Evidently, this is currently unfeasible for the densities and collapse times studied here.

Long computation times also place restrictions on the resolution that we can use with our KMC model. As is obvious from the CH₄ abundance, which we found to fluctuate between 1 and 0 on the grain, we would like to use larger lattice size N. However the simulation timescales with N²log(N), making higher resolution simulations unattainable for the present effort. Conversely, if one would like to reduce the computation time, the discrete fluctuations in abundance for some species shows that we are already at the limit of being able to resolve one particle in our system.

# 5. Conclusions

In this work, we have used the state-of-the-art astrochemical codes (MAGICKAL, MONACO, Nautilus, UCLCHEM, and KMC) to model the ice abundances measured with the JWST within the Ice Age ERS programme towards the highly extinguished background stars NIR38 and J110621 in the line of sight of the Chamaeleon I dense molecular cloud (see McClure et al. 2023). To do this, we have considered a small grid of physical parameters covering gas densities n<sub>H</sub>=2×10<sup>4</sup>, 2×10<sup>5</sup>, and 2×10<sup>6</sup> cm<sup>-3</sup>; dust temperatures of T<sub>dust</sub>=10, 12, and 14 K; and cloud collapse timescales of t<sub>col</sub>=10<sup>4</sup>, 10<sup>5</sup>, and 10<sup>6</sup> yr. The comparison of the model predictions with the observed JWST ice abundances towards NIR38 and J110621 allow us to understand the chemical processes behind the ice formation of CH<sub>3</sub>OH (and hence of CO and CO<sub>2</sub>, which are closely related to the ice chemistry of CH<sub>3</sub>OH) and of hydride species such as H<sub>2</sub>O, CH<sub>4</sub>, and NH<sub>3</sub>. The main conclusions of our work are as follows:

1. The best match between the predicted and the observed ice abundances is found for the higher-density and slower collapse models with n<sub>H</sub>≥2×10<sup>5</sup> cm<sup>-3</sup> and t<sub>col</sub>≥10<sup>5</sup> yr.
Low-density ($n_H=2 \times 10^4 \text{ cm}^{-3}$) fast models ($t_{\text{col}}=10^4$ yr) are ruled out because their collapse timescales are shorter than the depletion timescales preventing the efficient growth of the ices. These results nicely match the densities and timescales required by grain growth models to grow molecular cloud dust grains to sizes of 0.9 µm as measured towards the Chamaeleon I cloud;

2. The models provide mixed values for the CO depletion factor, which range from $f_D \sim 2-3$ to $f_D \ge 25$. Large-scale maps of gas-phase CO are therefore needed to constrain the level of CO depletion onto dust grains;

3. All our models find that CH$_3$OH ice formation towards the highly extinguished lines of sight of NIR38 and J110621 in the Chamaeleon I cloud, predominantly occurs via the successive hydrogenation of CO. The contribution of other proposed CH$_3$OH formation mechanisms, such as the reaction chains starting from CH$_4$ + OH or from C + H$_2$CO, is minor for the high-extinction conditions probed towards the NIR38 and J110621 field stars;

4. Our modelling effort shows that CO$_2$ ice abundance tends to be overpredicted in all models, but significant discrepancies are found across codes regarding the dominant reaction for CO$_2$ ice formation (i.e. CO + OH reaction forming an HOCO intermediate vs. CO + O → CO$_2$ reaction). Our KMC simulations, however, reveal that both reactions should be considered in any ice chemical modelling given that they contribute equally to the formation of CO$_2$ ice as a result of the higher mobility of atomic O on the grain compared to OH;

5. The formation of H$_2$O ice tends to be underpredicted in all models, although the predicted values lie within a factor of 10 with respect to the JWST observations. CH$_4$ is largely underproduced for all codes except for UCLCHEM, which shows a higher gas-phase abundance of atomic C during the translucent cloud phase due to more efficient non-thermal desorption. The models reproduce fairly well the observed ice abundance of NH$_3$;

6. Models with dust temperatures $T_{\text{dust}}=10$ K show significant changes in the results for codes with diffusive chemistry only compared to the results with $T_{\text{dust}}=12/14$ K (i.e. Nautilus and UCLCHEM). However, no differences are found for codes that consider non-diffusive chemistry. This is due to the fact that non-diffusive mechanisms take over at lower $T_{\text{dust}}$, compensating for the less efficient diffusiveness of atomic species and radicals at these low dust temperatures.

As discussed in Sect. 4.8, our chemical modelling does not consider the effects of grain growth and/or ice morphology (Dartois et al. 2024; Noble et al. 2024), cloud dynamics (Clément et al. 2023), or the presence of multiple ice components along the line of sight. Future chemical modelling will have to be carried out to assess the importance of all these effects in detail. Furthermore, spatial information on the gas-phase abundances of these major ice compounds will be obtained in the near future to complement the JWST Ice Age observations, which will help us constrain key parameters in the modelling of the chemistry of ices (such as the depletion of CO, $f_D$) towards the Chamaeleon I molecular cloud.

**Acknowledgements.** I.J.-S. would like to thank all the modellers involved in this project, who made it an enjoyable and stimulating experience. We also would like to thank Sergio Ioppolo and Giulia Perotti for insightful discussions and useful comments. I.J.-S. acknowledges funding from the European Research Council (ERC) under European Union's Horizon Europe research and innovation programme OPENS (101125858). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European

Union nor the granting authority can be held responsible for them. I.J.-S. and A.M. also acknowledge funding from grant PID2022-136814NB-I00, and grant PRE2019-091471, funded by the Spanish Ministry of Science, Innovation and Universities / State Agency of Research, MCIN/AEI/10.13039/501100011033, and by “ERDF/EU, A way of making Europe” and “ESF, Investing in your future”. V.W., A.T., J.A.N. and E.D acknowledge support from the French program “Physique et Chimie du Milieu Interstellaire” (PCMI) of the CNRS/INSU with the INC/INP cofunded by the CEA and CNES. MJ was supported by the NASA Planetary Science Division Internal Scientist Funding Program through the Fundamental Laboratory Research work package (FLaRe). R.T.G. thanks the National Science Foundation for funding through the Astronomy & Astrophysics program (grant number 2206516). The work by A.I.V. and K.B. has been supported via the project number FEUZ-2025-0003. M.N.D. acknowledges the Holcim Foundation Stipend. S.V. acknowledges funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program MOPPEX 833460. Part of this research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (80NM0018D0004). D.C.L. acknowledges financial support from the National Aeronautics and Space Administration (NASA) Astrophysics Data Analysis Program (ADAP). Astrochemistry at the Open University is currently supported by STFC under grant agreement numbers ST/X001164/1, ST/T005424/1 and ST/Z510087/1 and H.J.F. also thanks The Open University for supporting the OU Astrochemistry Group. Astrochemistry in The Netherlands is supported by the NWO Dutch Astrochemistry Network (grant no. ASTRO.JWST.001).

## References

Altwegg, K., Balsiger, H., Bar-Nun, A., et al. 2016, *Sci. Adv.*, 2, e1600285

Arasa, C., van Hemert, M. C., van Dishoeck, E. F., & Kroes, G. J. 2013, *J. Phys. Chem. A*, 117, 7064

Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, *ARA&A*, 47, 481

Belloche, A., Schuller, F., Parise, B., et al. 2011, *A&A*, 527, A145

Boogert, A. C. A., Huard, T. L., Cook, A. M., et al. 2011, *ApJ*, 729, 92

Boogert, A. C. A., Chiar, J. E., Knez, C., et al. 2013, *ApJ*, 777, 73

Boogert, A. C. A., Gerakines, P. A., & Whittet, D. C. B. 2015, *ARA&A*, 53, 541

Bouilloud, M., Fray, N., Bénilan, Y., et al. 2015, *MNRAS*, 451, 2145

Brown, P. D., Charnley, S. B., & Millar, T. J. 1988, *MNRAS*, 231, 409

Campbell, I. M., & Handy, B. J. 1978, *J. Chem. Soc., Faraday Trans.*, 1, 74, 316

Chang, Q., Cuppen, H. M., & Herbst, E. 2007, *A&A*, 469, 973

Chu, L. E. U., Hodapp, K., & Boogert, A. 2020, *ApJ*, 904, 86

Clément, A., Taillard, A., Wakelam, V., et al. 2023, *A&A*, 675, A165

Coutens, A., Commerçon, B., & Wakelam, V. 2020, *A&A*, 643, A108

Cuppen & Herbst 2007, *ApJ*, 668, 294

Cuppen, van Dishoeck, E. F., Herbst, E., & Tielens, A. G. G. M. 2009, *A&A*, 508, 275

Cuppen, Karssemeijer, L. J., & Lamberts, T. 2013, *Chem. Rev.*, 113, 8840

Cuppen, H. M., Walsh, C., Lamberts, T., et al. 2017, *Space Sci. Rev.*, 212, 1

Dartois, E. 2005, *Space Sci. Rev.*, 119, 293

Dartois, E., Chabot, M., Bacmann, A., et al. 2020, *A&A*, 634, A103

Dartois, E., Chabot, M., Id Barkach, T., et al. 2021, *A&A*, 647, A177

Dartois, E., Noble, J. A., Caselli, P., et al. 2024, *Nat. Astron.*, 8, 359

Drozdovskaya, M. N., van Dishoeck, E. F., Rubin, M., Jørgensen, J. K., & Altwegg, K. 2019, *MNRAS*, 490, 50

Drozdovskaya, M. N., Schroeder I, I. R. H. G., Rubin, M., et al. 2021, *MNRAS*, 500, 4901

Frerkin, M. A., Langer, W. D., & Wilson, R. W. 1982, *ApJ*, 262, 590

Fuente, A., Navarro, D. G., Caselli, P., et al. 2019, *A&A*, 624, A105

Garrod, R. T. 2008, *A&A*, 491, 239

Garrod, R. T. 2013, *ApJ*, 765, 60

Garrod, R. T. 2019, *ApJ*, 884, 69

Garrod, R. T., & Pauly, T. 2011, *ApJ*, 735, 15

Garrod, R. T., Wakelam, V., & Herbst, E. 2007, *A&A*, 467, 1103

Garrod, R. T., Widicus Weaver, S. L., & Herbst, E. 2008, *ApJ*, 682, 283

Garrod, R. T., Vasyunin, A. I., Semenov, D. A., Wiebe, D. S., & Henning, T. 2009, *ApJ*, 700, L43

Garrod, R. T., Jin, M., Matis, K. A., et al. 2022, *ApJS*, 259, 1

Gerakines, P. A., Schutte, W. A., Greenberg, J. M., & van Dishoeck, E. F. 1995, *A&A*, 296, 810

Goto, M., Bailey, J. D., Hocuk, S., et al. 2018, *A&A*, 610, A9

Goto, M., Vasyunin, A. I., Giuliano, B. M., et al. 2021, *A&A*, 651, A53

Goumans, T. P. M., & Andersson, S. 2010, *MNRAS*, 406, 2213

Goumans, T. P. M., Uppal, M. A., & Brown, W. A. 2008, *MNRAS*, 384, 1158

Gusdorf, A., Cabrit, S., Flower, D. R., & Pineau Des Forêts, G. 2008, *A&A*, 482, 809

Hasegawa, T. I., & Herbst, E. 1993, *MNRAS*, 261, 83

Hasegawa, T. I., Herbst, E., & Leung, C. M. 1992, *ApJS*, 82, 167
He, J., Shi, J., Hopkins, T., Vidali, G., & Kaufman, M. J. 2015, ApJ, 801, 120

Herbst, E., & Cuppen, H. M. 2006, PNAS, 103, 12257

Hiraoka, K., Yamashita, A., Yachi, Y., et al. 1995, ApJ, 443, 363

Hocuk, S., Szűcs, L., Caselli, P., et al. 2017, A&A, 604, A58

Holdship, J., Viti, S., Jiménez-Serra, I., Makrymallis, A., & Priestley, F. 2017, ApJ, 154, 38

Ioppolo, S., van Boheemen, Y., Cuppen, H. M., van Dishoeck, E. F., & Linnartz, H. 2011, MNRAS, 413, 2281

Ioppolo, S., Fedoseev, G., Lamberts, T., Romanzin, C., & Linnartz, H. 2013, Rev. Sci. Instrum., 84, 073112

Ioppolo, S., Fedoseev, G., Chuang, K.-J., et al. 2020, Nat. Astron., 5, 197

Iqbal, W., & Wakelam, V. 2018, A&A, 615, A20

Ishibashi, A., Molpeceres, G., Hidaka, H., et al. 2024, ApJ, 976, 162

Jiménez-Serra, I., Caselli, P., Martín-Pintado, J., & Hartquist, T. W. 2008, A&A, 482, 549

Jiménez-Serra, I., Viti, S., Quénard, D., & Holdship, J. 2018, ApJ, 862, 128

Jiménez-Serra, I., Vasyunin, A. I., Spezzano, S., et al. 2021, ApJ, 917, 44

Jin, M., & Garrod, R. T. 2020, ApJS, 249, 26

Jin, M., Lam, K. H., McClure, M. K., et al. 2022, ApJ, 935, 133

Karssemeijer, L. J., & Cuppen, H. M. 2014, A&A, 569, A107

Keto, E., & Caselli, P. 2010, MNRAS, 402, 1625

Knez, C., Boogert, A. C. A., Pontoppidan, K. M., et al. 2005, ApJ, 635, L145

Lacy, J. H., Sneden, C., Kim, H., & Jaffe, D. T. 2017, ApJ, 838, 66

Lamberts, T., de Vries, X., & Cuppen, H. M. 2014, Faraday Discuss., 168, 327

Lamberts, T., Fedoseev, G., Kästner, J., Ioppolo, S., & Linnartz, H. 2017, A&A, 599, A132

Lamberts, T., Fedoseev, G., van Hemert, M. C., et al. 2022, ApJ, 928, 48

Lebreuilly, U., Vallucci-Goy, V., Guillet, V., Lombart, M., & Marchand, P. 2023, MNRAS, 518, 3326

Lefèvre, C., Pagani, L., Juvela, M., et al. 2014, A&A, 572, A20

Marchand, P., Lebreuilly, U., Mac Low, M. M., & Guillet, V. 2023, A&A, 670, A61

Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425

McClure, M. K., Rocha, W. R. M., Pontoppidan, K. M., et al. 2023, Nat. Astron., 7, 431

McElroy, Walsh, C., Markwick, A. J., et al. 2013, A&A, 550, A36

Melnick, G. J., Tolls, V., Goldsmith, P. F., et al. 2012, ApJ, 752, 26

Minissale, M., Congiu, E., Manicò, G., Pirronello, V., & Dulieu, F. 2013, A&A, 559, A49

Minissale, M., Loison, J.-C., Baouche, S., et al. 2015, A&A, 577, A2

Minissale, M., Dulieu, F., Cazaux, S., & Hocuk, S. 2016, A&A, 585, A24

Minissale, M., Aikawa, Y., Bergin, E., et al. 2022, ACS Earth Space Chem., 6, 597

Molpeceres, G., Kästner, J., Fedoseev, G., et al. 2021, J. Phys. Chem. Lett., 12, 10854

Molpeceres, G., Enrique-Romero, J., & Aikawa, Y. 2023, A&A, 677, A39

Nguyen, T., Baouche, S., Congiu, E., et al. 2018, A&A, 619, A111

Noble, J. A., Dulieu, F., Congiu, E., & Fraser, H. J. 2011, ApJ, 735, 121

Noble, J. A., Congiu, E., Dulieu, F., & Fraser, H. J. 2012, MNRAS, 421, 768

Noble, J. A., Fraser, H. J., Aikawa, Y., Pontoppidan, K. M., & Sakon, I. 2013, ApJ, 775, 85

Noble, J. A., Fraser, H. J., Smith, Z. L., et al. 2024, Nat. Astron., 8, 1169

Oba, Y., Watanabe, N., Kouchi, A., Hama, T., & Pirronello, V. 2010, ApJ, 712, L174

Ormel, C. W., Paszun, D., Dominik, C., & Tielens, A. G. G. M. 2009, A&A, 502, 845

Ormel, C., Min, M., Tielens, A., Dominik, C., & Paszun, D. 2011, A&A, 532, A43

Pendleton, Y. J., Tielens, A. G. G. M., Tokunaga, A. T., & Bernstein, M. P. 1999, ApJ, 513, 294

Penteado, E. M., Walsh, C., & Cuppen, H. M. 2017, ApJ, 844, 71

Perotti, G., Rocha, W. R. M., Jørgensen, J. K., et al. 2020, A&A, 643, A48

Potapov, A., Krasnokutski, S. A., Jäger, C., & Henning, T. 2021, ApJ, 920, 111

Qasim, D., Chuang, K. J., Fedoseev, G., et al. 2018, A&A, 612, A83

Qasim, D., Lamberts, T., He, J., et al. 2019, A&A, 626, A118

Qasim, D., Fedoseev, G., Chuang, K. J., et al. 2020, Nat. Astron., 4, 781

Quénard, D., Ilee, J. D., Jiménez-Serra, I., et al. 2018a, ApJ, 868, 9

Quénard, D., Jiménez-Serra, I., Viti, S., Holdship, J., & Coutens, A. 2018b, MNRAS, 474, 2796

Raut, U., & Baragiola, R. A. 2011, ApJ, 737, L14

Rawlings, J. M. C., Hartquist, T. W., Menten, K. M., & Williams, D. A. 1992, MNRAS, 255, 471

Ruaud, M., Loison, J. C., Hickson, K. M., et al. 2015, MNRAS, 447, 4004

Ruaud, M., Wakelam, V., & Hersant, F. 2016, MNRAS, 459, 3756

Ruffle, D. P., & Herbst, E. 2001, MNRAS, 324, 1054

Santos, J. C., Chuang, K.-J., Lamberts, T., et al. 2022, ApJ, 931, L33

Senevirathne, B., Andersson, S., Dulieu, F., & Nyman, G. 2017, Mol. Astrophys., 6, 59

Silsbee, K., Ivlev, A. V., Sipilä, O., Caselli, P., & Zhao, B. 2020, A&A, 641, A39

Simons, Lamberts, T., & Cuppen, H. M. 2020, A&A, 634, A52

Smith, Z. L., Dickinson, H. J., & Fraser, H. 2025, Nat. Astron., submitted

Sturm, J. A., McClure, M. K., Beck, T. L., et al. 2023, A&A, 679, A138

Sturm, J. A., McClure, M. K., Harsono, D., et al. 2024, A&A, 689, A92

Väisälä, M. S., Harju, J., Mantere, M. J., et al. 2014, A&A, 564, A99

Vasyunin, A. I., & Herbst, E. 2013, ApJ, 769, 34

Vasyunin, A. I., Caselli, P., Dulieu, F., & Jiménez-Serra, I. 2017, ApJ, 842, 33

Viti, S., Collings, M. P., Dever, J. W., McCoustra, M. R. S., & Williams, D. A. 2004, MNRAS, 354, L141

Wakelam, V., & Herbst, E. 2008, ApJ, 680, 371

Wakelam, V., Herbst, E., & Selsis, F. 2006, A&A, 451, 551

Wakelam, V., Dartois, E., Chabot, M., et al. 2021, A&A, 652, A63

Wakelam, V., Gratier, P., Loison, J. C., et al. 2024, A&A, 689, A63

Watanabe, N., & Kouchi, A. 2002, ApJ, 571, L173

Weingartner, J. C., & Draine, B. T. 2001, ApJ, 548, 296

Westenberg, A. A., & DeHaas, N. 1972, J. Phys. Chem., 76, 2215

Whittet, D. C. B., Bode, M. F., Longmore, A. J., Baines, D. W. T., & Evans, A. 1983, Nature, 303, 218

Whittet, D. C. B., Gerakines, P. A., Tielens, A. G. G. M., et al. 1998, ApJ, 498, L159

Whittet, D. C. B., Cook, A. M., Chiar, J. E., et al. 2009, ApJ, 695, 94

Yang, Y.-L., Green, J. D., Pontoppidan, K. M., et al. 2022, ApJ, 941, L13

<footnote n="1">¹ Centro de Astrobiología (CAB), CSIC-INTA, Ctra. de Ajalvir km 4, 28850 Torrejón de Ardoz, Spain</footnote>

<footnote n="2">² Institute for Molecules and Materials, Radboud University, Heydaaladsweeg 135, 6525 AJ Nijmegen, The Netherlands</footnote>

<footnote n="3">³ Astrochemistry Laboratory, Code 691, NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA</footnote>

<footnote n="4">⁴ Department of Physics, Catholic University of America, Washington, DC 20064, USA</footnote>

<footnote n="5">⁵ Laboratoire d'astrophysique de Bordeaux, Univ. Bordeaux CNRS, B18N, allée Geoffroy Saint-Hilaire, 33615 Pessac, France</footnote>

<footnote n="6">⁶ Institute of Natural Sciences and Mathematics, Ural Federal University, 19 Mira Str., 620075 Ekaterinburg, Russia</footnote>

<footnote n="7">⁷ Max-Planck-Institut für extraterrestrische Physik, Gießenbacherstrasse 1, 85748 Garching bei München, Germany</footnote>

<footnote n="8">⁸ Department of Physics, University of Central Florida, Orlando, FL 32816, USA</footnote>

<footnote n="9">⁹ Institut des Sciences Moléculaires d'Orsay, CNRS, Univ. Paris-Saclay, 91405 Orsay, France</footnote>

<footnote n="10">¹⁰ Physique des Interactions Ioniques et Moléculaires, CNRS, Aix Marseille Univ., 13397 Marseille, France</footnote>

<footnote n="11">¹¹ Leiden Observatory, Leiden University, PO Box 9513, 2300 RA Leiden, The Netherlands</footnote>

<footnote n="12">¹² Transdisciplinary Research Area (TRA) ‘Matter’/Argelander-Institut für Astronomie, University of Bonn, Germany</footnote>

<footnote n="13">¹³ Physics and Astronomy, University College London, Gower Street, London WC1E 6BT, UK</footnote>

<footnote n="14">¹⁴ Institute of Astronomy of the RAS, Pyatnitskaya st. 48, 119017, Moscow</footnote>

<footnote n="15">¹⁵ Department of Astronomy and Chemistry, University of Virginia, Charlottesville, VA 22904, USA</footnote>

<footnote n="16">¹⁶ Leiden Institute of Chemistry, Gorlaeus Laboratories, Leiden University, PO Box 9502, 2300 RA Leiden, The Netherlands</footnote>

<footnote n="17">¹⁷ School of Physical Sciences, The Open University, Kents Hill, Milton Keynes MK7 6AA, UK</footnote>

<footnote n="18">¹⁸ Center for Astrophysics | Harvard & Smithsonian, 160 Concord Avenue, Cambridge 02138 MA, USA</footnote>

<footnote n="19">¹⁹ Physikalish-Meteorologisches Observatorium Davos und Weltstrahlungszentrum (PMOD/WRC), Dorfstrasse 33, 7260 Davos Dorf, Switzerland</footnote>

<footnote n="20">²⁰ Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91109, USA</footnote>
## Appendix A: Additional grain-surface chemical reactions included in UCLCHEM for HOCO and CO₂

For the UCLCHEM calculations, we have added several chemical reactions involved in the grain-surface chemical network of the hydrocarboxyl radical (HOCO) and CO₂. These reactions are listed in Table A.1, and they have been extracted from the works by Goumans et al. (2008), Goumans & Andersson (2010), Arasa et al. (2013), Ruaud et al. (2015), Qasim et al. (2019), and Molpeceres et al. (2023).

<table><caption>**Table A.1.** Grain-surface reactions included in UCLCHEM for HOCO and CO₂.</caption><thead><tr><th colspan="2">Reaction</th><th>Energy barrier (K)</th></tr></thead><tbody><tr><td>CO + O →</td><td>CO₂</td><td>2500</td></tr><tr><td rowspan="2">CO + OH →</td><td>HOCO (50%)</td><td>150</td></tr><tr><td>H + CO₂ (50%)</td><td>150</td></tr><tr><td rowspan="2">N + HOCO →</td><td>NH + CO₂ (50%)</td><td>0</td></tr><tr><td>OH + OCN (50%)</td><td>0</td></tr><tr><td>O + HOCO →</td><td>OH + CO₂</td><td>0</td></tr><tr><td rowspan="3">H + HOCO →</td><td>H₂ + CO₂ (70%)</td><td>0</td></tr><tr><td>CO + H₂O (20%)</td><td>0</td></tr><tr><td>HCOOH (10%)</td><td>0</td></tr></tbody></table>

**Notes.** The third reaction, forming H and CO₂, was already present in the default chemical network of UCLCHEM, but we updated the energy barrier from 1000 K to 150 K (Ruaud et al. 2015).

## Appendix B: Evolution of the ice and gas phase abundances for each code

We present here the evolution of the ice abundances obtained during the collapse (phase 1) stage by the MONACO, Nautilus, UCLCHEM and KMC models for H₂O, CO, CO₂, CH₃OH, CH₄, and atomic C. Four representative models are shown corresponding to the cases of a low-density fast-collapse model and of a high-density slow-collapse model for both field stars NIR38 and J110621. In these figures we also show the evolution of the gas phase abundances of CO and atomic C.

The evolution is similar for all codes: the ices undergo a fast growth by the end of the collapse for the high-density slow-collapse models, when the density rapidly increases and the catastrophic depletion of CO takes place (see right panels in these figures). We note that the evolution of the density follows the free fall equations, meaning that the increase of density is very rapid towards the end of the calculations. As a consequence, for the models with final $n_H$ densities of $2 \times 10^4$ cm⁻³, the increase from about $2 \times 10^3$ cm⁻³ to $2 \times 10^4$ cm⁻³ is done very rapidly (regardless the collapse timescales), and thus the ice chemical composition does not have time to evolve and it reflects the chemistry at very low densities inherited from the beginning of the collapse. Therefore, the ices do not experience significant growth.

## Appendix C: Predicted ice composition and comparison with Ice Age observations

Figures C.1-C.8 show a comparison between the ice composition predicted by the MONACO, Nautilus, UCLCHEM and KMC codes, and the ice abundances reported by McClure et al. (2023) towards the NIR38 and J110621 field stars.

## Appendix D: Evaluating the error of the models

To evaluate the performance of the models in a systematic way, we used the same algorithm for all models. First, we select the models that predict H₂O column densities closer to the Ice Age observations, that is, the models with the lowest difference or error. Then, we calculate the mean absolute error (MAE) between the predicted and observed abundances with respect to H₂O on a logarithmic scale. We note that we have used the H₂O ice column densities for the first step of the MAE calculations, because the H₂O ice feature is the one with the highest signal-to-noise ratio in the JWST observations, and therefore it presents the lowest associated uncertainties.

The equation used for the MAE is

$$ \text{MAE} = \frac{1}{N} \sum_{i=1}^{N} \left| \log_{10}(y_i) - \log_{10}(\hat{y}_i) \right|, \quad (\text{D.1}) $$

where $y_i$ represents the observed value of species i (in abundance with respect to H₂O), $\hat{y}_i$ is the value predicted by the astrochemical model, and N is the number of species observed with the JWST. We propagate the uncertainty in the values reported by McClure et al. (2023) to derive the uncertainty in the MAE. The best model corresponds to the one that presents the lowest MAE value within the estimated uncertainty coming from the observations.

As an example, Fig. D.1 shows the H₂O error and the MAE for all the NIR38 models using UCLCHEM, while Fig. D.2 presents the same but for J110621. The error bars correspond to the 3σ level. For NIR38, we can see that there are four models with the lowest MAE value (within the uncertainties), but the one with the lowest final density ($2 \times 10^4$ cm⁻³) presents an error for H₂O column density which is too high and clearly larger than the rest, so we discard it. Then, the best models are the ones with intermediate collapse times ($t_{col} = 10^5$ yr) and intermediate to high final densities ($n_H = 2 \times 10^5$, $2 \times 10^6$ cm⁻³), and also the one with high final density and fast collapse. However, if we now repeat the same procedure but with J110621 (Fig. D.2), we obtain the same models except for the last one with fast collapse. Therefore, we conclude that the two best models for UCLCHEM are, for both sources, the ones with intermediate collapse times and intermediate to high final densities (second and fifth models from the right).

We obtained analogous figures for the remaining codes used in this paper. The results of the MAE analysis are reported in Table 5 of Sect. 3.3. We also note that we performed some tests using the root mean squared error (RMSE) and found similar results. The main difference is that the RMSE penalises the models that have bigger discrepancies between observations and predictions for some molecules. Therefore, it would yield a big error if one species were very wrongly predicted, while it would favour models with a more uniform error between the different species predictions. For the RMSE, the equation used is

$$ \text{RMSE} = \left( \frac{1}{N} \sum_{i=1}^{N} (\log_{10}(y_i) - \log_{10}(\hat{y}_i))^2 \right)^{\frac{1}{2}}. \quad (\text{D.2}) $$
<figure><img src="image_7.png" /><figcaption>**Fig. B.1.** Time evolution for the species closely related to the formation of CH₃OH and hydrides in the ice predicted by some of our MONACO models (the extreme cases) for the collapse phase.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. B.2.** Time evolution for the species closely related to the formation of CH₃OH and hydrides in the ice predicted by some of our NAUTILUS models (the extreme cases) for the phase 1 collapse phase.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Fig. B.3.** Time evolution for the species closely related to the formation of CH₃OH and hydrides in the ice predicted by some of our UCLCHEM models (the extreme cases) for the collapse phase.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. B.4.** Time evolution for the species closely related to the formation of CH₃OH and hydrides in the ice predicted by some of our KMC models (extremes) for the collapse phase. We note that the heavy fluctuations observed for CH₄ is a finite-size effect, due to the very low concentration numbers generated by the KMC on the grain (CH₄ abundance is found to fluctuate between 1 and 0 on our 50 × 50 lattice). For the high-density slow-collapse model (model 9) towards NIR38 (*T* = 14 K, *A*<sub>V</sub> = 60 mag), the results are shown for 95.9 % of the intended collapse time.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Fig. C.1. Predicted ice compositions by MONACO for NIR38. Details are as for Fig. 3.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. C.2. Predicted ice compositions by MONACO for J110621. Details are as for Fig. 3.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Fig. C.3.** Predicted ice compositions by Nautilus for NIR38. Details are as for Fig. 3. Missing bars indicate ice abundance ratios ≤0.03% relative to H₂O.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. C.4.** Predicted ice compositions by Nautilus for J110621. Details are as for Fig. 3. Missing bars indicate ice abundance ratios ≤0.03% relative to H₂O.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Fig. C.5. Ice compositions predicted by UCLCHEM for NIR38. Details are as in Fig. 3.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. C.6. Ice compositions predicted by UCLCHEM for J110621. Details are as in Fig. 3. Missing bars indicate ice abundance ratios ≤0.03% relative to H₂O.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Fig. C.7.** Predicted ice compositions by the KMC simulations for NIR38. For models 6 ($n_H = 2 \times 10^5$ cm$^{-3}$, $t_{col}=10^6$ yr), 8 ($n_H = 2 \times 10^6$ cm$^{-3}$, $t_{col}=10^5$ yr) and 9 ($n_H = 2 \times 10^6$ cm$^{-3}$, $t_{col}=10^6$ yr), the final abundances are taken from simulations that completed 96.8, 99.5 and 94.7 % of their intended collapse times, respectively. Missing bars indicate ice abundance ratios $\le 0.03\%$ relative to H$_2$O.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Fig. C.8.** Predicted ice compositions by the KMC simulations for J110621. Missing bars indicate ice abundance ratios $\le 0.03\%$ relative to H$_2$O.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**Fig. D.1.** Logarithmic errors between the predicted ice compositions by the different models of UCLCHEM for $A_V$ = 60 mag with $T_{\text{dust}}$ = 14 K. Errors for only the H$_2$O column densities are indicated with blue rhombus, while black squares mark the errors calculated with the abundances with respect to water, both superimposed to the bars that represent the models in previous figures. Error bars correspond to $3\sigma$.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Fig. D.2.** Same as Fig. D.1, but for J110621 ($A_V$ = 95 mag). The two best models correspond to intermediate collapse times and intermediate/high final densities (second and fifth models from the right).</figcaption></figure>