© The Authors 2025

1 Introduction

WASP-107b is a transiting exoplanet orbiting a K7V-type host star at a distance of 0.055 AU, resulting in an equilibrium temperature of ∼740 K (Anderson et al. 2017). With a mass of 0.12 M_Jup and a radius of 0.90 R_Jup (Piaulet et al. 2021), its extremely low gravity and large scale height make it well suited for atmospheric characterization. A series of transit observations have revealed the structure of WASP-107b’s deeper atmosphere in unprecedented detail. Observations with the Hubble Space Telescope (HST) in the near-infrared show water vapour features, compressed by high-altitude clouds (Kreidberg et al. 2018). Further observations with James Webb Space Telescope’s (JWST) MIRI instrument (LRS) were obtained by Dyrek et al. (2024), who further constrained the type of aerosols to silicate cloud particles (i.e. MgSiO₃, SiO₂, and SiO). Additionally, the detection of SO₂ is indicative of photochemistry in a metal-enriched atmosphere (Tsai et al. 2023; Dyrek et al. 2024).

Neither dataset shows proof of CH₄, despite predictions from equilibrium chemistry showing that methane should be abundant at temperatures below ∼1000 K (Moses et al. 2011). Vertical mixing drives the atmosphere out of equilibrium and can quench CH₄ in the deep atmosphere, thereby fixing the methane content in higher layers to the volume mixing ratio (VMR) at the quench point. In addition to possible carbon depletion by a sub-solar carbon-to-oxygen ratio, Dyrek et al. (2024) show that the combination of a high intrinsic temperature and strong vertical mixing reduces the CH₄ content to undetectable levels without affecting the formation of SO₂ in the upper layers (see their Extended Data – Figure 3).

More recently, Welbanks et al. (2024) analysed JWST NIRCam spectra in combination with previous datasets and managed to detect and constrain the abundances of H₂O, SO₂, CO, CO₂, NH₃, and even traces of CH₄. These results were accompanied by the publication of Sing et al. (2024), who used a JWST NIRSpec spectrum to detect H₂O, SO₂, CO, CO₂, and CH₄. Both of the above studies conclude that the depleted CH₄ is the result of a high intrinsic temperature (>350 K) and atmospheric mixing causing deep quenching. Additionally, Welbanks et al. (2024) attribute this hot interior to eccentricity-driven tidal heating, which had previously been identified as a potential explanation for the inflated radius of WASP-107b (Thorngren et al. 2019; Fortney et al. 2020; Yu & Dai 2024).

Atmospheric retrievals aim to infer basic properties from observations by navigating a large, degenerate parameter space using parametrized, easy-to-compute models and Bayesian inference techniques (Madhusudhan 2018). However, such simplifications can lead to biases for inferred parameters, which have been documented and investigated in the past (e.g. Benneke & Seager 2012; Griffith 2014; Line & Parmentier 2016; Heng & Kitzmann 2017; Fisher & Heng 2018; Welbanks & Madhusudhan 2019; MacDonald et al. 2020; Novais et al. 2025). Alternatively, one can construct more complex forward models that take into account, for example, chemical kinetics (e.g. in Dyrek et al. 2024), radiative-convective equilibrium, cloud microphysics, and fluid dynamics. A major downside of forward models is the large computational cost, which makes them unsuited for retrievals. Combining the best of both worlds, Bell et al. (2023) introduced a 1D radiative-convective photochemicalequilibrium (1D-RCPE) retrieval (Beatty et al. 2024; Fu et al. 2024; Welbanks et al. 2024). By pre-computing a grid of physics-informed forward models, and allowing for interpolation in between, 1D-RCPE retrievals can leverage the strength of Bayesian inference with on-the-fly computation of low-cost forward models. This technique was applied to WASP-107b by Welbanks et al. (2024), who deduced a high intrinsic temperature based on the CH₄ deficiency.

Given WASP-107b’s inflated radius and low density, the atmospheric elemental composition and intrinsic temperature can help to unravel properties such as the planet’s core mass, and formation history (Madhusudhan et al. 2016; Anderson et al. 2017; Piaulet et al. 2021; Welbanks et al. 2024; Sing et al. 2024). To achieve this, it is crucial to deduce the atmospheric properties from transit spectra as accurately as possible. Various methods have been adopted in the above studies, which lead to different results. More specifically for WASP-107b, reported values for atmospheric metallicity range from ≳ 6 Z_⊙ (Dyrek et al. 2024) up to 43 ± 8 Z_⊙ (Sing et al. 2024), and for carbon-to-oxygen ratios from ${0.51}_{-0.21}^{+0.27}$ (solar) (Sing et al. 2024) down to ${0.33}_{-0.05}^{+0.06}$ (Welbanks et al. 2024). Additionally, multiple studies on WASP107b point towards a high intrinsic temperature and deep quenching based on the lack of strong CH₄ features.

In this paper, we aim to assess how well 1D-RCPE retrievals can reliably infer atmospheric properties such as metallicity (Z), carbon-to-oxygen ratio (C/O), intrinsic temperature (T_int), and eddy diffusion coefficient (K_zz) for WASP-107b. More specifically, we explore the sensitivity of our results to different cloud parametrizations, and investigate the information content of the available observations. A secondary goal of this paper is to consistently analyse the relevant chemical pathways and summarize the chemical trends with respect to the four grid parameters (Z, C/O, T_int, and K_zz).

In Sect. 2, we describe the methods and numerical codes used in this study. We start by constructing a grid of 900 1D-RCPE models and analyse the results in Sect. 3. We then proceed by running 1D-RCPE retrievals with various set-ups in Sect. 4. Finally, in Sect. 5, we discuss how our results affect previously inferred properties of WASP-107b before concluding in Sect. 6.

2 Methods

We used a number of different codes to build and analyse our grid of 1D-RCPE models. First, we outline the design of our model grid (Sect. 2.1). This is followed by description of our model set-up for the pressure-temperature profiles with PICASO (Sect. 2.2) and the chemical kinetics simulations with VULCAN (Sect. 2.3). We analysed the resulting grid of 900 1D-RCPE models with Dijkstra’s algorithm (Sect. 2.4). To generate synthetic transit spectra, we used petitRADTRANS (Sect. 2.5) with various cloud parametrizations (Sect. 2.6). Finally, we fitted the 1D-RCPE model grid to different datasets (Sect. 2.7) with UltraNest (Sect. 2.8).

2.1 Grid design

We varied four parameters in our grid of forward models: Z, C/O, T_int, and K_zz. We ran atmospheric models with metallicities 1, 2, 3, 5, 7, 10, 15, 20, and 30 Z_⊙ (Asplund et al. 2009). We discarded sub-solar values based on previous findings and the detection of SO₂, which points to metal enrichment (Kreidberg et al. 2018; Dyrek et al. 2024). The upper limit was chosen based on the interior structure models that estimate the planet’s maximal bulk metallicity (Kreidberg et al. 2018; Thorngren & Fortney 2019). We explored C/O ratios of 0.05, 0.11, 0.23, and 0.46, corresponding to roughly 10%, 25%, 50%, and 100% of the solar value, respectively (Lodders 2010), in agreement the lower boundary suggested by planet formation models (Öberg et al. 2011; Mordasini et al. 2016; Espinoza et al. 2017). We considered intrinsic temperatures of 100, 200, 300, 400, and 500 K. The lower values correspond to the predictions where only deposition of stellar radiation is considered (Thorngren et al. 2019), while higher values correspond to increasing effects of tidal heating (Leconte et al. 2010; Fortney et al. 2020; Millholland et al. 2020). We explored eddy diffusion coefficients of 10⁶, 10⁷, 10⁸, 10⁹, and 10¹⁰ cm²/s, corresponding to atmospheres with and without quenching. Combined, this results in a grid of 900 forward models. Throughout the grid, we used parameters that are representative for the WASP-107(b) system. For WASP-107b, these include a surface gravity of 270 cm/s², radius of 0.94 R_Jup, mass of 0.10 M_Jup, orbital separation of 0.055 AU, and equilibrium temperature of 740 K. For its host star, we adopted ${\mathrm{T}}_{\text{eff}}^{\star }=4430\mathrm{K}$, a radius of 0.66 R_⊙, and log g(cgs) ≃ 4.6 (Piaulet et al. 2021).

2.2 Radiative-convective equilibrium: PICASO

We used the open-source 1D radiative-convective code PICASO (Mukherjee et al. 2023) to compute 180 temperature-pressure profiles of the atmosphere of WASP-107b, each with a different metallicity, carbon-to-oxygen ratio, and intrinsic temperature. We ran PICASO with 91 pressure log-spaced levels between 300 and 10⁻⁷ bar. The code uses pre-mixed correlated-k opacities that are calculated under the assumption of equilibrium chemistry (Lupu et al. 2023; Mukherjee et al. 2023). As is listed in Lupu et al. (2023), these include C₂H₂ (Wilzewski et al. 2016; Rothman et al. 2013), C₂H₄ (Rothman et al. 2013; Marley et al. 2021), C₂H₆ (Rothman et al. 2013; Marley et al. 2021), CH₄ (Yurchenko et al. 2013; Yurchenko & Tennyson 2014; Pine 1992; Wenger & Champion 1998), CO (Rothman et al. 2010; Gordon et al. 2017; Li et al. 2015), CO₂ (Huang et al. 2014), CrH (Burrows et al. 2002), Fe (Ryabchikova et al. 2015; O’Brian et al. 1991; Fuhr et al. 1988; Bard et al. 1991; Bard & Kock 1994), FeH (Dulick et al. 2003; Hargreaves et al. 2010), H₂ (Gordon et al. 2017), ${\mathrm{H}}_3^{+}$ (Mizus et al. 2017), H₂O Polyansky et al. (2018); Barton et al. (2017), H₂S (Azzam et al. 2016; Rothman et al. 2013; Marley et al. 2021), HCN (Harris et al. 2006; Barber et al. 2014; Gordon et al. 2022), LiCl (Bittner & Bernath 2018), LiF (Bittner & Bernath 2018), LiH (Coppola et al. 2011), MgH (GharibNezhad et al. 2013; Yadin et al. 2012; Gharib-Nezhad et al. 2021), N₂ (Rothman et al. 2013; Marley et al. 2021), NH₃ (Yurchenko et al. 2011; Wilzewski et al. 2016; Marley et al. 2021), OCS (Gordon et al. 2017), PH₃ (Sousa-Silva et al. 2014; Marley et al. 2021), SiO (Barton et al. 2013; Kurucz 1992; Marley et al. 2021), TiO (McKemmish et al. 2019; Gharib-Nezhad et al. 2021), and VO (McKemmish et al. 2019; Gharib-Nezhad et al. 2021; Marley et al. 2021), in addition to alkali metals (Li, Na, K, Rb, Cs) (Ryabchikova et al. 2015; Allard et al. 2019, 2016, 2007a, b). For a planet such as WASP107b, H₂O and potentially CH₄ are expected to be the dominant infrared absorbers (Mollière et al. 2015).

In order to converge to a solution with a smooth transition between the radiative and convective layers, PICASO requires an initial guess of the pressure level of the radiative-convective boundary zone (RCB). As the RCB is highly sensitive to T_int (Sarkis et al. 2021), we varied the RCB pressure level for each model with different intrinsic temperatures to ensure said convergence. The pre-computed opacity mixtures in PICASO do not include metallicities 7 and 15 Z_⊙. Therefore, we linearly interpolated these PT structures from neighbouring metallicities (5 to 10 Z_⊙, and 10 to 20 Z_⊙, respectively). For models with C/O=0.05, we used the same computation as C/O=0.11 as such low carbon-to-oxygen ratios are not available in the opacity data. Finally, we chose to compute day-side averages by setting r_st=1, which controls the contribution of stellar radiation to the net radiative flux in every atmospheric layer (see Eq. (20) in Mukherjee et al. 2023), so that the temperature in the upper, quasi-isothermal layers matches the reported values in Dyrek et al. (2024). We further discuss this assumption in Sect. 5. Fig. 1 shows a selection of the computed PT profiles with PICASO.

2.3 (Photo)-chemical modelling: VULCAN

We used the open-source chemical kinetics code VULCAN (Tsai et al. 2017, 2021) to model the molecular composition of WASP-107b. We used the default chemical network SNCHO_PHOTO_NETWORK that includes 89 species composed of H, C, O, N, S, coupled by 1030 forward and backward thermochemical reactions. Vertical mixing was included via diffusion, with accompanying eddy diffusion coefficient K_zz. We used a constant K_zz profile for each model, with the values specified in Sect. 2.1, to reduce complexity and maintain good interpretability of our results.

Photochemistry was included by way of 60 photodissociation reactions for the main UV absorbing species in H₂−dominated atmospheres. We used a spectral energy distribution (SED) of HD 85512, as a proxy for host star WASP-107. HD 85512 is a K7V star with a similar effective temperature as WASP-107, and its SED is available in the MUSCLES database (France et al. 2016; Youngblood et al. 2016; Loyd et al. 2016). Note that the SEDs of stars with similar spectral types are not necessarily comparable, as high-energy (X-ray, UV) radiation fluxes are linked to stellar age, rotation period and overall magnetic activity, rather than photospheric emission. With a rotation period of 17 ± 1 days (Anderson et al. 2017) and age of 3.4 ± 0.7 Gyr (Piaulet et al. 2021), WASP-107 is expected to be more active compared to HD 85512, which has a ∼47 day rotation period and a resulting age estimate of 5.6 ± 0.1 Gyr (Tsantaki et al. 2013). However, observations in the NUV and X-ray have shown that WASP-107 is slightly less active than HD 85512, although this difference is less than a factor of two in the UV so that HD 85512 remains a good proxy for WASP-107 (Dyrek et al. 2024).

To initialise VULCAN, we used 150 log-spaced vertical layers between 10 and 10⁻⁷ bar for models with high T_int to ensure convergence in the deep layers, while this upper pressure boundary can be extended to 300 bar for low T_int to capture the quenching points. To fix individual convergence issues, we lowered parameters ATOL and/or RTOL (absolute and relative tolerance thresholds for errors during the numerical integration), or systematically decreased the pressure boundary at the interior while making sure the quench points of key species (i.e. CH₄ and NH₃) are captured in the model. Otherwise, we used the default set-up of VULCAN. The above adjustments were especially necessary for models with low K_zz, where the quench points are located much higher than their high- K_zz counterparts.

The elemental mixture was adjusted based on the given combination of Z and C/O, before the simulation is started from chemical equilibrium abundances. We first scaled the solar elemental mixture (Asplund et al. 2009) with the Z, before lowering the carbon abundance to the requested C/O. We did not include processes such as condensation, ion chemistry, or atmospheric escape. Finally, we highlight that VULCAN assumes a fixed thermal background during the numerical integration. Studies have shown that the introduction of disequilibrium chemistry can lead to temperature differences up to 100 K in the PT profile (Drummond et al. 2016; Agúndez 2025).

Fig. 1 Selection of the 180 computed PT profiles with PICASO for C/O=0.46. The line transparency corresponds to intrinsic temperature, from a full line for Tint=100 K to almost transparent for Tint=500 K. The shaded grey region indicates layers where the radiative-convective boundaries are located in all 180 PT profiles.

2.4 Chemical pathway analysis: Dijkstra’s algorithm

Alongside of chemical abundances, we are interested in the reaction rates in the steady-state solution of VULCAN. This enables us to analyse the chemical pathways that govern the production and destruction of key species in our model. In particular we are interested in the dominant or fastest chemical pathways, which translates to finding the shortest chain of inverse reaction rates that connect two molecules. This is equivalent to finding the shortest path in graph theory, where the atmospheric constituents are nodes and reactions the connecting edges. Therefore, we implement Dijkstra’s algorithm (Dijkstra 1959), following the implementation outlined in Appendix B of Tsai et al. (2018).

2.5 Radiative transfer: petitRADTRANS

We used the open-source radiative transfer code petitRADTRANS (version 2.7.7) (Mollière et al. 2019) to compute synthetic transmission spectra from the chemical and radiative-convective models. We included correlated-k line absorption opacities of SO₂ (Underwood et al. 2016), H₂S (Azzam et al. 2016), CO₂ (Yurchenko et al. 2020), CO (Kurucz 1993; Rothman et al. 2010), NH₃ (Coles et al. 2019), CH₄ (Yurchenko et al. 2017), and H₂O (Polyansky et al. 2018) available via ExoMolOP (Chubb et al. 2021) and HITEMP (Rothman et al. 2010). To reduce computational cost, we left out line absorption opacity of SO, C₂H₂, OH, and HCN as we found no substantial contribution to the spectra of our grid. Furthermore, we included Rayleigh scattering opacity of H₂ (Dalgarno & Williams 1962) and He (Chan & Dalgarno 1965), and continuum opacity from collision induced absorption (CIA) by H₂−H₂ and H₂−He (Borysow et al. 1988, 1989; Borysow & Frommhold 1989; Borysow et al. 2001; Borysow 2002; Richard et al. 2012).

2.6 Cloud opacity parametrizations

We considered different types of clouds when post-processing our grid of models to synthetic spectra. First, we considered a grey opacity source based at a certain pressure, P_cloud, which effectively blocks the layers below from forming spectral features. Second, we considered non-grey clouds by adding an additional opacity source to the radiative transfer models. Following Dyrek et al. (2024), we added silicates consisting of crystalline MgSiO₃, SiO₂, or amorphous SiO cloud particles that are irregularly shaped, for which the opacity was computed with the DHS (distribution of hollow spheres) method (Min et al. 2005). Due to computational restrictions, we only considered one of the above species at a time. At the cloud base pressure (P_cloud), we assumed a mass fraction (X_cloud) that decreases with altitude following X_cloud(P/P_cloud)^f_sed. Further based on Ackerman & Marley (2001), petitRADTRANS computes a distribution of cloud particle sizes based on the sedimentation efficiency f_sed, width of the log-normal size distribution σ_g (=1+2 x_σ), and vertical eddy diffusion coefficient ${\mathrm{K}}_{\mathrm{z}\mathrm{z}}^{(\text{cloud})}$. Although K_zz is already a fitting parameter in our grid of 1D-RCPE models, we considered a different eddy diffusion coefficient for the clouds (Zhang 2020). The sedimentation efficiency is defined as the ratio of cloud particle sedimentation velocity to turbulent vertical mixing speed. Lower values (f_sed ≲ 1) are indicative of an extended cloud deck, composed of small sub-micron particles while higher values (f_sed ≳ 3–5) point to a compact cloud deck where larger particles (≳ 1 μm) have settled efficiently. Initial tests revealed that x_σ has minimal effect on the final fit. To reduce computation time, we set x_σ = −1. This totals an additional four free parameters (P_cloud, X_cloud, f_sed, and ${\mathrm{K}}_{\mathrm{z}\mathrm{z}}^{(\text{cloud})}$) to be retrieved.

Third, we adopted an alternative cloud treatment, presented by Welbanks et al. (2024), who parametrize a Gaussian centred around the 10 μm silicate cloud feature. Such a parametrization acknowledges the presence of an additional non-grey opacity source without further specifying its origin. The parametrization is given by $$K_{\text{cloud}}(\lambda )=2{\kappa }_{\text{opac}}\varphi \left(\frac{\lambda -{\lambda }_0}{\omega }\right)\mathrm{\Phi }\left(\xi \frac{\lambda -{\lambda }_0}{\omega }\right),$$(1)

where κ_opac is the base cloud opacity, ϕ the standard normal probability density function, Φ the standard normal cumulative density function, λ₀ the mean of the Gaussian, ω the standard deviation, and ξ a skewness parameter. After initial testing, and in agreement with Welbanks et al. (2024), we fixed ξ=0 to reduce the run time of our retrievals with this setup. This effectively reduces the parametrization of Eq. (1) to a Gaussian function, and adds three additional parameters to the fitting routine (κ_opac, λ₀, and ω). Henceforth, we refer to this parametrization as ‘Gaussian’ clouds.

Finally, we parametrized a wavelength-dependent scattering slope by $$K_{\text{scatt}}(\lambda )={\kappa }_{\text{scatt}}{\left(\frac{\lambda }{{\lambda }_{\text{ref}}}\right)}^{{\gamma }_{\text{scatt}}},$$(2)

where κ_scatt is the opacity at λ_ref=0.35 μm, and γ_scatt = −4 in the case for Rayleigh scattering.

2.7 Data source and model offsets

We fitted the synthetic transmission spectra of our chemical models to the archival transit spectra of WASP-107b. Specifically, this includes HST WFC3 (1.1–1.6 μm) (Kreidberg et al. 2018), JWST NIRCam (2.5–5.0 μm) (Welbanks et al. 2024), and JWST MIRI/LRS (5.0–12.2 μm) (Dyrek et al. 2024). For the latter dataset, we used the reduction presented in Welbanks et al. (2024).

Aside from the four grid parameters (Z, C/O, T_int, and K_zz), and cloud properties, we also fit for offsets between the HST (Δ_HST) and NIRCam (Δ_NIR) data with respect to the MIRI data. The latter is necessary to compensate for differences in instrument sensitivities and assumptions during the data reduction process (Carter et al. 2024), which can affect the retrieval outcome (Yip et al. 2021). Finally, we fitted a reference pressure (P_ref) to the planetary radius to account for an offset between the model spectra and data. We note that some studies choose to fit a planetary radius at a fixed reference pressure, but this choice has no significant effect on the retrieval outcome (Welbanks & Madhusudhan 2019).

2.8 Nested sampling: UltraNest

To perform the 1D-RCPE retrieval, we used the Python package UltraNest (Buchner 2021), which implements the nested sampling algorithm for Bayesian inference. It estimates the Bayesian evidence (marginal likelihood) and generates approximate posterior distributions of all model parameters. The algorithm maintains a set of live points that explore the parameter space, iteratively replacing the live point with the lowest likelihood by a new sample with a higher likelihood drawn from the (remaining) prior distributions. This approach enables UltraNest to navigate complex likelihood surfaces, making it particularly suitable for atmospheric retrieval problems with degenerate and multi-modal posterior distributions.

At each likelihood evaluation, a set of parameters is drawn from the prior distributions. For each combination of Z, C/O, T_int, and K_zz, we quad-linearly interpolate the PT profile, mean-molecular weights, and VMRs for each chemical species from our pre-computed grid of 1D-RCPE models on a common pressure grid (150 log-spaced layers between 10 and 10⁻⁷ bar) before computing the spectrum on-the-fly with petitRADTRANS. Our default set-up consists of 100 to 200 livepoints (depending on the size of the parameter space), and a convergence criterion DLOGZ =0.8. We use high-performance computing facilities and run our retrievals in parallel on a AMD EPYC 76433.64 GHz CPU, on 80 cores with 250 GB RAM memory or on a Intel Xeon Platinum 8360Y 2.4 GHz CPU, on 130 cores with 232 GB RAM memory. Run times vary from less than a day to up to two weeks.

3 Disequilibrium chemistry and pathway analysis

Before presenting our results of the atmospheric retrievals, we analyse the chemistry simulations of our 1D-RCPE grid of forward models. More specifically, we discuss the dominant production and destruction pathways of key molecular species that have been detected in WASP-107b. These include CH₄, SO₂, NH₃, CO, CO₂, and H₂O (Kreidberg et al. 2018; Dyrek et al. 2024; Welbanks et al. 2024; Sing et al. 2024). Additionally, we discuss trends with respect to the grid parameters (Z, C/O, T_int, and K_zz) and highlight regions where the composition deviates from chemical equilibrium (i.e. the minimal Gibbs free energy state at local pressure and temperature) by disequilibrium effects such as vertical mixing and photochemistry. The goal of this section is to develop a comprehensive understanding of the chemistry in WASP-107b, building on existing literature that has examined a broader range of planetary atmospheres beyond those that resemble WASP-107b.

Fig. 2 shows the behaviour of SO₂, CH₄, and NH₃ with respect to the four different grid parameters (Z, C/O, T_int, and K_zz), as these molecules provide information on the sulphur (Sect. 3.1), carbon (Sect. 3.2), and nitrogen chemistry (Sect. 3.3) in WASP-107b.

3.1 Sulphur chemistry: H₂S, SO, and SO₂

The presence of SO₂ in WASP-107b’s atmosphere is evidence of active photochemistry. Outlined by Tsai et al. (2023) for WASP39 b, a photochemical production pathway of SO₂ is $${\mathrm{H}}_2\mathrm{O}\stackrel{\mathrm{h}v}{⟶}\mathrm{O}\mathrm{H}+\mathrm{H}$$(3) $${\mathrm{H}}_2\mathrm{O}+\mathrm{H}\to \mathrm{O}\mathrm{H}+{\mathrm{H}}_2$$(4) $${\mathrm{H}}_2\mathrm{S}+\mathrm{H}\to \mathrm{S}\mathrm{H}+{\mathrm{H}}_2$$(5) $$\mathrm{S}\mathrm{H}+\mathrm{H}\to \mathrm{S}+{\mathrm{H}}_2$$(6) $$\mathrm{S}+\mathrm{O}\mathrm{H}\to \mathrm{S}\mathrm{O}+\mathrm{H}$$(7) $$\mathrm{S}\mathrm{O}+\mathrm{O}\mathrm{H}\to {\mathrm{S}\mathrm{O}}_2+\mathrm{H}$$(8) $$\text{Net:}\mathrm{S}+2{\mathrm{H}}_2\mathrm{O}\to {\mathrm{S}\mathrm{O}}_2+3{\mathrm{H}}_2.$$(9)

This shows that OH, produced by hydrogen abstraction and photolysis of H₂O, oxidizes atomic sulphur, which is thermochemically produced from H₂S via hydrogen abstraction (Hobbs et al. 2021; Tsai et al. 2023). However, for lower temperature atmospheres (<1000 K), Reaction (4) would produce insufficient OH to form detectable oxidized S-species (Tsai et al. 2023). The detection of SO₂ in WASP-107b, with T_eq ≈ 740 K, was enabled by its low gravity (∼270 cm/s²) and favourable stellar UV insolation (Dyrek et al. 2024). With water photolysis as the dominant source in the upper atmosphere (≲ 10⁻⁵ bar), Dyrek et al. (2024) argue that the primary source of OH between 10⁻⁵ and 1 bar stems from photolysis of various abundant molecules (including H₂O, NH₃, N₂) aided by Reaction (4) (see Suppl. Inf. – Figure 14 in Dyrek et al. 2024). A recent study by de Gruijter et al. (2025) further investigated the chemical production pathways of SO₂ and, while confirming the two distinct pressure-dependent regimes for OH production, revealed an additional complexity regarding the source of atomic sulphur. Depending on the UV flux and pressure level, sulphur is liberated by hydrogen abstraction of SH (Reaction (6)) and/or by photolysis of SH and S₂ that liberate sulphur from the H₂S reservoir.

From Fig. 2, it is clear that SO₂ is most sensitive to the atmospheric metallicity. The VMRs show a positive correlation with metallicity that primarily results from the larger sulphur reservoir (H₂S), but also from the increase in hydroxide (OH), in line with previous studies (Polman et al. 2023; Dyrek et al. 2024). We validate the formation pathway of SO₂ along this dimension of our grid, using models with C/O=0.46, T_int=400 K, and log K_zz=9 (cgs), similar parameters as presented in Dyrek et al. (2024). We find that OH radicals are produced by hydrogen abstraction of H₂O (Reaction (4)) throughout the atmosphere and direct photolysis of H₂O (Reaction (3)), where the latter dominates the upper layers. At solar metallicity (Z = 1 Z_⊙), the region where Reaction (3) dominates extends up to several millibar. For extremely metal-rich models (Z > 10 Z_⊙), this region is limited to only the very upper layers (≲ μbar), while Reaction (4) dominates OH production elsewhere. This indicates that radiation penetrates to shallower pressures in higher-metallicity atmospheres, which results from the additional metal opacity that increases the overall optical depth at any given pressure layer. Hence, for higher atmospheric metallicities, the layer at which the atmosphere becomes optically thick moves upward.

We find that the production pathway from OH to SO₂ depends on metallicity and pressure in our models. At solar metallicity (Z=1 Z_⊙), the entire atmospheric column follows the production pathway summarized by Reaction (9). Given the quadratic dependence of sulphur atoms on metallicity, we see an increased reaction rate for reactions with exclusively sulphur reactants for higher Z models (>7 Z_⊙). In deeper layers of such high-metallicity atmospheres, SO is formed predominantly by $$\mathrm{S}+\mathrm{S}\mathrm{H}\to \mathrm{H}+{\mathrm{S}}_2$$(10) $${\mathrm{S}}_2+\mathrm{O}\mathrm{H}\to \mathrm{S}\mathrm{O}+\mathrm{S}\mathrm{H},$$(11)

where S arises predominantly from Reaction (6) or SH recombination $$\mathrm{S}\mathrm{H}+\mathrm{S}\mathrm{H}\to {\mathrm{H}}_2\mathrm{S}+\mathrm{S}.$$(12)

Finally, in high-metallicity environments, the subsequent production of SO₂ results primarily from $$\mathrm{S}\mathrm{O}+\mathrm{S}\mathrm{O}\to {\mathrm{S}\mathrm{O}}_2+\mathrm{S}$$(13)

rather than Reaction (8).

One can further analyse the last step of the production pathway (SO → SO₂) by looking at the reaction rates in Fig. 3, which shows the fastest reactions (by reaction rate) involving SO₂ and SO. At solar metallicity (left panel, Fig. 3), Reaction (8) (forward) is faster than Reaction (13) (forward), making it the dominant SO₂−producing reaction. By dominant, we mean that removing this particular reaction from the network will significantly disturb the steady-state equilibrium abundances. Furthermore, both reactions run much slower in their reverse directions. This indicates that, when isolating these reactions, a net production of SO₂ is achieved. To reach a steady state, and constant VMR, direct photolysis (${\mathrm{S}\mathrm{O}}_2{\displaystyle \stackrel{\mathrm{h}v}{⟶}\mathrm{S}\mathrm{O}+\mathrm{O}}$) is the main sink of SO₂ and balances the system almost entirely.

For Z = 30 Z_⊙ (middle panel, Fig. 3), between ∼ mbar and ∼100 μbar, the opposite situation occurs where Reaction (13) is slightly faster (in both directions) than Reaction (8), making it the dominant reaction that controls the VMR of SO₂. This statement is valid for both production and destruction of SO₂, as Reaction (13) is nearly perfectly balanced in both directions as a result of little to no direct photolysis in these regions. This reasoning does not hold for a small region between 10⁻² and 10⁻³ bar, where it appears that both Reactions (8) and (13) are destroying SO₂ at a faster reaction rate than producing it. It is also evident that no other thermo- or photochemical reactions can compensate this net destruction of SO₂, as we have plotted the fastest reaction at play. An explanation for the imbalance in SO₂ production and destruction is found in the disequilibrium effect of vertical mixing. Indeed, when comparing to simulations with lower K_zz (right panel, Fig. 3) where the atmosphere is not quenched in this particular model, we see that both reactions are perfectly balanced in both directions. For quenched atmospheres, such as with log K_zz=9(cm²/s), we see the immediate effect of SO₂ replenishment due to vertical upward mixing. We conclude that the dominant SO₂ production pathway is influenced by atmospheric metallicity and that vertical mixing plays a crucial role in transporting SO₂ upward.

One can also further analyse the species whose sulphur ends up in SO (and later SO₂). Both pathways via Reaction (7) as well as Reactions (10) and (11), require atomic sulphur whose origin is debated in literature (Tsai et al. 2023; de Gruijter et al. 2025). Fig. 4 shows that we find Reaction (6) to be dominant throughout the atmosphere with low metallicity (e.g. Z=1 Z_⊙), while Reaction (12) becomes faster at higher metallicity for pressures above ∼1 mbar. Furthermore, both the forward and backward reaction rates of Reaction (6) are comparable, thus acting as the primary source and sink of SH in the regions where the VMR of SO₂ peaks (see Fig. 2). Further analysis of Reaction (6) reveals that it is even so slightly out of equilibrium (<2%) favouring the backwards direction. However, we believe this is not an argument against the ability of Reaction (6) to liberate S from SH (de Gruijter et al. 2025). Instead, this indicates that additional S-producing reactions (of secondary importance) compensate for this slight imbalance to reach a steady-state equilibrium. Finally, our models show that hydrogen abstraction of H₂S (Reaction (5)) is the fastest way to produce SH from the H₂S reservoir. The different chemical pathways that we have identified are summarized in Table 1.

This paints a complete picture of the dominant chemical pathways that lead to SO₂ in our models of WASP-107b. One final remark can be made regarding Reaction (4), which converts H₂O into OH by hydrogen abstraction. It is clear that an excess of H is available throughout the atmosphere, whose origin is photochemical. However, Dyrek et al. (2024) revealed the somewhat puzzling result that models with photochemistry of non-hydrogen bearing species can still produce significant amounts of SO₂. For example, the authors show (Suppl. Inf. – Figure 14 in Dyrek et al. 2024) that only N₂−photolysis is sufficient to produce SO₂, but do not specify the chemical pathway. This calls for further investigation on how sufficient amounts of H, and subsequently OH via Reaction (4), can be produced in the SO₂−producing layers from photolysis of non-hydrogen bearing species. In Fig. 5, we use our model set-up to run three models where we consider only N₂, H₂O, or NH₃ photochemistry. When considering only N₂ as a photon absorber, we validate that the products of N₂ photolysis can enable SO₂ production, with the help of the subsequent thermochemical reaction,

$${\mathrm{H}}_2+\mathrm{N}\to \mathrm{H}+\mathrm{N}\mathrm{H},$$(14)

which acts as the dominant H-producing reaction by interacting with the abundant H₂ reservoir. We do note, however, that there are several faster H-producing reactions when the final steady state is reached in our simulations. However, we are specifically interested in pathways that start from species that are abundant in chemical equilibrium, from which the chemistry simulation is initialized, such as H₂.

When we only include H₂O photolysis (Reaction (3)), Fig. 5 shows that both OH and H are deposited deep in the atmosphere, which activates Reaction (4) in both directions. In models with similar parameters where we include all available photolysis reactions (Sect. 3.1), stellar radiation is absorbed much higher up in the atmosphere. When we only allow for H₂O photolysis, radiation can penetrate much deeper due to the absence of other photon absorbers in this specific model.

Compared to Dyrek et al. (2024), we produce much less SO₂ when only activating NH₃ photolysis. However, our network analysis reveals that radiation penetrates deep in the atmosphere, releasing H there via both of the following reactions: $${\mathrm{N}\mathrm{H}}_3\stackrel{\mathrm{h}v}{⟶}{\mathrm{N}\mathrm{H}}_2+\mathrm{H}$$(15) $${\mathrm{N}\mathrm{H}}_3\stackrel{\mathrm{h}v}{⟶}\mathrm{N}\mathrm{H}+\mathrm{H}+\mathrm{H}.$$(16)

Once again, we note that such deep penetration depth can only be reached in absence of other photon absorbers. The lack of H deposition in the upper layers, and resulting lack of SO₂ with respect to Dyrek et al. (2024) can be explained due to the absence of sub- 100 nm photoabsorption cross-sections in our chemical network. We also note that vertical mixing cannot efficiently transport H atoms upwards, due to their short chemical timescales, but can successfully quench SO₂ in these layers. Additionally, the photochemical product NH₂ can locally produce H atoms by NH₂+H₂ → NH₃+H. In reality, a mixture of N₂, H₂O, NH₃ and other abundant photochemically active species all contribute to liberating H atoms throughout the atmosphere although Reaction (4) remains the dominant source.

Fig. 2 VMR profiles of key species CH4, SO2, and NH3 in our grid of disequilibrium chemistry models of WASP-107b. Chemical equilibrium is plotted with thin, dashed lines. The nominal model of Z=10Z⊙, C/O=0.46, Tint=400 K, and log Kzz=9(cm2/s) was chosen to be in agreement with previous studies (Dyrek et al. 2024; Welbanks et al. 2024).

Fig. 3 Dominant reactions (by fastest reaction rate) that involve SO and SO2, for three different combinations of metallicity and eddy diffusion coefficient. Forward rates are shown by solid lines; backward rates by dotted lines. The other values of these simulations are Tint=400 K and C/O=0.46.

Fig. 4 Dominant reactions (by fastest reaction rate) that produce/destroy S, for three different metallicities. Forward rates are shown by solid lines; backward rates by dotted lines. The other values of the simulations are Tint=400 K, C/O=0.46 and log Kzz=9(cm2/s).

Fig. 5 VMR of SO2 and dominant H-producing reactions for a nominal model (Z=10 Z⊙, C/O=0.46, Tint=400 K, log Kzz=9 cm2/s) that includes all available photodissociation cross-sections (upper left panel), and three limiting cases where only the photolysis of one particular source is considered (N2, H2O., or NH3), copying the experiment of Dyrek et al. (2024). Panels that show thermochemical reaction rates show both their forward (solid lines) and backward (dashed) rates.

3.2 Carbon chemistry: C₂H₂, CO, CO₂, and CH₄

From Fig. 2, it is clear that all four dimensions of the parameter space (Z, C/O, T_int, K_zz) impact the VMR profile of CH₄. Furthermore, each panel shows the three well-known regimes from the top to bottom layers: photochemical destruction, transport-induced quenching, and thermochemical equilibrium. Aside from the general carbon budget, controlled by metallicity and carbon-to-oxygen ratio, the VMR of CH₄ is mainly controlled by the VMR around the quench pressure, i.e. the point at which the chemical timescale matches the vertical mixing timescale.

The absence of abundant CH₄ in WASP-107b has led to the suggestion of a tidally heated core (high T_int), in combination with strong vertical mixing (high K_zz) that quenches CH₄ deep in the atmosphere and thus reduces the overall CH₄ content in the observable layers (Fortney et al. 2020; Dyrek et al. 2024; Welbanks et al. 2024; Sing et al. 2024). This is also apparent from Fig. 6, which shows the quenched VMR and pressures of CH₄ in our chemistry grid. At high T_int(>300 K), faster vertical mixing indeed quenches CH₄ more deeply (higher pressure), resulting in an overall lower VMR. The opposite is true for low T_int(≤ 200 K), where we see that deeper quenching results in more CH₄. This is due to a negative VMR gradient for layers in chemical equilibrium, as compared the positive gradient in models with high T_int (see also Fig. 2).

Furthermore, Fig. 6 shows an upwards shift in altitude (downwards in pressure) of the quench point at higher metallicity, together with an overall lower quenched VMR of CH₄. The negative correlation between metallicity and CH₄ is counterintuitive as one expects a positive correlation between the available carbon budget and eventual VMR of CH₄. However, this is explained by Fig. 1, which shows that the deep layers are warmer at higher Z as a result by of the increased abundance (and thus opacity) of strong absorbers (Drummond et al. 2018). Additionally, chemistry in high-metallicity environments shifts the CH₄−CO boundary to lower temperatures (Lodders & Fegley 2002). A warmer interior, for similar eddy diffusion coefficients, and high metallicity explain the upwards shift in altitude of the quench point, resulting in less CH₄.

The formation of CH₄ is coupled to CO, as both molecules can act as the dominant carbon reservoir under different circumstances. We find that in the majority of our models, CO dominates over CH₄. Only in set-ups with (i) quasi-solar metallicity and/or (ii) shallow quenching by low T_int or weak K_zz, is CH₄ more abundant than CO. For our nominal model of Z = 10 Z_⊙, C/O=0.46, T_int=400 K, and log K_zz = 9 (cgs), we analyse the conversion pathway from CH₄ to CO at the quench pressure and find $${\mathrm{C}\mathrm{H}}_4+\mathrm{H}\to {\mathrm{C}\mathrm{H}}_3+{\mathrm{H}}_2$$(17) $${\mathrm{C}\mathrm{H}}_3+\mathrm{O}\mathrm{H}\to {\mathrm{C}\mathrm{H}}_2\mathrm{O}\mathrm{H}+\mathrm{H}$$(18) $${\mathrm{C}\mathrm{H}}_2\mathrm{O}\mathrm{H}+\mathrm{M}\to \mathrm{H}+{\mathrm{H}}_2\mathrm{C}\mathrm{O}+\mathrm{M}$$(19) $${\mathrm{H}}_2\mathrm{C}\mathrm{O}+\mathrm{H}\to \mathrm{H}\mathrm{C}\mathrm{O}+{\mathrm{H}}_2$$(20) $$\mathrm{H}\mathrm{C}\mathrm{O}+\mathrm{M}\to \mathrm{H}+\mathrm{C}\mathrm{O}+\mathrm{M}$$(21) $$\text{Net:}{\mathrm{C}\mathrm{H}}_4+{\mathrm{H}}_2\mathrm{O}\to \mathrm{C}\mathrm{O}+3{\mathrm{H}}_2$$(22)

taking into account that OH is produced by Reaction (4). This pathway has been well documented in previous studies (Moses et al. 2011; Venot et al. 2012; Zahnle & Marley 2014; Tsai et al. 2017, 2018; Venot et al. 2020). Some of the above papers have reported dominant pathways that run via CH₃OH, effectively replacing Reaction (18) with $${\mathrm{C}\mathrm{H}}_3+\mathrm{O}\mathrm{H}+\mathrm{M}\to {\mathrm{C}\mathrm{H}}_3\mathrm{O}\mathrm{H}+\mathrm{M}$$(23) $${\mathrm{C}\mathrm{H}}_3\mathrm{O}\mathrm{H}+\mathrm{H}\to {\mathrm{C}\mathrm{H}}_2\mathrm{O}\mathrm{H}+{\mathrm{H}}_2$$(24)

for high-pressure (deep quenching) environments. We do not find a model in which this pathway is faster than Reaction (18). A faster pathway exists via COS, a stable sulphur molecule present in low amounts at chemical equilibrium (VMR < 10⁻⁷ in our nominal model), which produces CO via COS+CH₃ → CO+CH₃ S. However, COS itself is predominantly kept in equilibrium by reactions H+COS ↔ CO+SH and COS+SH ↔ CO+HS₂, which both have CO as a reactant, disqualifying COS from the interconversion pathway CH₄−CO. Finally, we find that multiple other pathways can come close to competing with the pathway of Reaction (22). This includes, for example, stripping CH₄ to atomic carbon before oxidization to CO. We refer to Tsai et al. (2018) and Hu (2021) for a detailed description on the various CH₄−CO pathways and regimes in which they occur.

Other carbon-species of interest are CO₂ (detected in WASP-107b) and C₂H₂ (undetected). Production and destruction of CO₂ is equilibrated by $$\mathrm{O}\mathrm{H}+\mathrm{C}\mathrm{O}\to {\mathrm{C}\mathrm{O}}_2+\mathrm{H}$$(25)

throughout the entire atmosphere across all dimensions of our model grid, and its VMR is hardly altered by disequilibrium chemistry (see Fig. A.1). However, we note photochemically produced OH has been shown to increase CO₂ in some exoplanetary atmospheres (Moses et al. 2011; Zahnle & Marley 2014; Hu 2021). Finally, we briefly discuss C₂H₂ as it another important photochemical product in exoplanetary atmospheres (Line et al. 2010; Moses 2014; Baeyens et al. 2022; Konings et al. 2022). We find the average VMR of C₂H₂ to be below 10⁻¹⁰ in all our models, making it too scarce to be observable. We do note that its primary formation pathway in our models is well known, following CH₄ → CH₃ → C₂H₆ → C₂ H₅ → C₂H₄ → C₂ H₃ → C₂H₂, as was described in Moses et al. (2011).

Fig. 6 Quenching VMR (left panel) and estimated quenching pressure (right panel) of CH4 versus intrinsic temperature Tint for various eddy diffusion coefficients in low (Z=1 Z⊙) and high (Z=10 Z⊙) metallicity atmospheres. The models with log Kzz = 6 (cgs) are not included as not all atmospheres are quenched at such weak vertical mixing.

3.3 Nitrogen chemistry: HCN, N₂ and NH₃

Dyrek et al. (2024) reported a tentative detection of NH₃ (2–3 σ) in WASP-107b at a VMR of ≲ 10^−5.47, based on the 10.5 μm feature in the MIRI data. Subsequently, Sing et al. (2024) found marginal evidence of NH₃, based on a 3 μm feature in the NIRSpec data. Most recently, Welbanks et al. (2024) could confirm the presence of NH₃ with a 6 σ detection based on a panchromatic analysis covering both the 3 and 10.5 μm feature at ≲ 10⁻⁵ VMR, while also requiring an enhancement of a factor ∼30 to match the observed VMR to their best-fit forward model. This sparks an interest in the behaviour and chemical pathway of NH₃ in our models.

Returning to Fig. 2, we observe only a slight dependency of NH₃ on metallicity despite the linear dependence with the available nitrogen budget. Similarly to the behaviour of CH₄, outlined in Sect. 3.2, one must take into account the PT profile and its vertical upward shift for higher metallicity as well as the shift in NH₃−N₂ balance (Lodders & Fegley 2002). Given the quasiconstant VMRs in the layers around the NH₃ quench point, the upward (or downward) shift of the quench point altitude (or pressure) does not substantially affect the quenched VMR. A similar explanation can be given to the behaviour with respect to K_zz as a stronger eddy diffusion coefficient, and thus deeper quenching, does not affect the quenched NH₃ VMR too much due to quasi-constant VMRs in chemical equilibrium (Moses et al. 2011; Moses 2014).

Similar to the CH₄−CO interconversion described in Sect. 3.2, the ratio between NH₃ and N₂ is frozen at the quench point and sets the dominant nitrogen reservoir in the middle atmosphere. From Fig. 7, it is clear that a higher intrinsic temperature results in a lower quenched VMR of NH₃ as the quench point moves up in altitude. Additionally, the minimal dependence on metallicity and K_zz is clearly shown here as well.

In all of our chemistry models, we find a clear overabundance of N₂ with respect to NH₃ by several orders of magnitude. Only in a limiting case with Z = 1 Z_⊙, T_int = 100 K, and high K_zz, the overabundance of N₂ is limited to a factor two. We identify the primary pathway for NH₃−N₂ conversion as $$2\times ({\mathrm{N}\mathrm{H}}_3+\mathrm{H}\to {\mathrm{N}\mathrm{H}}_2+{\mathrm{H}}_2)$$(26) $$\mathrm{H}+{\mathrm{N}\mathrm{H}}_2\to {\mathrm{H}}_2+\mathrm{N}\mathrm{H}$$(27) $$\mathrm{N}\mathrm{H}+{\mathrm{N}\mathrm{H}}_2\to {\mathrm{N}}_2{\mathrm{H}}_2+\mathrm{H}$$(28) $${\mathrm{N}}_2{\mathrm{H}}_2+\mathrm{H}\to {\mathrm{N}}_2\mathrm{H}+{\mathrm{H}}_2$$(29) $${\mathrm{N}}_2\mathrm{H}+\mathrm{M}\to {\mathrm{N}}_2+\mathrm{H}+\mathrm{M}$$(30) $$\text{Net:}2{\mathrm{N}\mathrm{H}}_3\to {\mathrm{N}}_2+3{\mathrm{H}}_2$$(31)

around the quench pressure, which is in agreement with the literature (Moses et al. 2011; Hu 2021). An alternative pathway can replace Reaction (28) with $${\mathrm{N}\mathrm{H}}_2+{\mathrm{N}\mathrm{H}}_3\to {\mathrm{N}}_2{\mathrm{H}}_3+{\mathrm{H}}_2$$(32) $${\mathrm{N}}_2{\mathrm{H}}_3+\mathrm{M}\to {\mathrm{N}}_2{\mathrm{H}}_2+\mathrm{H}+\mathrm{M},$$(33)

at higher quenching pressures but fails to be the quickest way to form N₂ from NH₃ in our models.

A final nitrogen-species of potential interest is HCN, as photochemistry can enhance its presence significantly from chemical equilibrium (Baeyens et al. 2024). However, HCN remains undetected in WASP-107b, and is not expected to be prominent in transit spectra at this temperature regime (Baeyens et al. 2022). In our models, HCN reaches an maximal VMR of ∼10⁻⁶ in some limiting cases (see Fig. A.1). It forms primarily from the net reaction CH₄+NH₃ → HCN+3 H₂, after hydrogen abstraction of both reactants (Moses et al. 2011), thus depending heavily on their abundance in the atmospheres.

Fig. 7 Same as Fig. 6, now for NH3.

4 1D-RCPE retrievals of WASP-107b

In this section, we run retrievals based on the grid of 1D-RCPE models. We start by briefly discussing the main molecular features in cloud-free transit spectra (Sect. 4.1). In Sect. 4.2, we fit with a simplified ‘grey’ cloud set-up for transit observations below 5 μm to gradually build up our understanding of the 1D-RCPE retrieval method. Afterwards, in Sect. 4.3, we combine all information of the previous sections and run retrievals with complex non-grey clouds on all available observations. For both Sects. 4.2 and 4.3, we present a nominal case in which the results serve as a benchmark.

4.1 Transmission spectra

With a thorough understanding of the chemical processes and pathways, we are equipped to transform these models into observables. As is described in Sect. 2, we computed transmission spectra of our models on the fly during the fitting routine. Fig. 8 shows cloud-free transmission spectra that probe the parameter space of our grid.

A number of molecular features show a clear variation in amplitude. The spectra are dominated by several broad H₂O bands, a distinct SO₂ bump around 4 μm and a double SO₂ feature at 7–9 μm, a CO₂ feature at 4.3 μm (and contribution at 2.8 μm), a small CO bump around 4.6 μm, and a narrow-peaked CH₄ feature at 3.4 μm. Note that the latter molecule has broad bands and can dominate the spectra in more methane-rich models. The VMRs of CO₂, CO, H₂O, and SO₂ all scale with metallicity to a certain degree and can therefore provide more opacity in the observable layers at higher metallicity. At the same time, a high-metallicity atmosphere is more compact due to the increased mean molecular weight, which tends to reduce feature amplitudes as the scale height decreases. Although both effects counteract each other, the upper panel of Fig. 8 shows a clear positive correlation between feature amplitude and metallicity. The opposite is true for CH₄(3.4 μm), for which the VMR decreases with metallicity (also see Fig. 2) but the feature strength remains the same for all considered metallicities.

At low C/O, less CH₄ is present in the atmosphere and more oxygen is stored in SO₂ and H₂O, compared to CO and CO₂. As a result, a sub-solar C/O gives a weaker CH₄ feature and more pronounced SO₂ and H₂O bands. Models with high T_int generally have less CH₄ and NH₃, which is directly visible in Fig. 8. The CH₄ feature at 3.4 μm diminishes in amplitude at higher T_int. In general, we do not see substantial contribution of NH₃ in these cloud-free spectra, with the exception around 10.5 μm in the models with the lowest T_int.

Finally, we briefly discuss the bottom panel of Fig. 8, where different values of K_zz are explored. While higher K_zz results in deeper quenching, and an overall lower amount of CH₄ in the bulk atmosphere, the effect on the spectral feature at 3.4 μm is small. The SO₂ features around 7–9 μm scale with K_zz, which is easily explained by its VMR in the upper layers. At higher K_zz, vertical mixing is more efficient in replenishing the upper layers with SO₂ as it competes with the photochemical destruction. Given the high transit depth of these SO₂ features, it is sensitive to exactly these upper layers in the chemistry models, which explains its behaviour in the spectrum. As we shall see later, the same mechanism can affect the 3.4 μm feature of CH₄, especially in the presence of a high altitude cloud deck that blocks the bulk CH₄ content below from forming a stronger spectral feature. We emphasize that Fig. 8 was computed for cloud-less spectra, whereas the next sections include clouds that substantially affect the (amplitude of) visible molecular features.

4.2 Retrievals with grey clouds on transit spectra below 5 μm

We proceed by fitting our grid of 1D-RCPE models to HST (Kreidberg et al. 2018) and JWST’s NIRCam (Welbanks et al. 2024) observations of WASP-107b. In this wavelength range, we can highly simplify the cloud treatment and avoid the discussion on more complex (silicate) clouds for now. The most simple parametrization of a cloud deck can be achieved by applying a grey opacity with a base on certain pressure (P_cloud). We ran our fitting routine (described in Sect. 2.8) on JWST’s NIRCam and HST data and allow for an offset between the two (Δ_HST). We included the full range of grid parameters for Z(1 to 30 Z_⊙), C/O(0.05 to 0.46), and T_int(100 to 500 K), with the exception of K_zz(10⁷ to 10¹⁰ cm²/s). We excluded models with K_zz=10⁶ cm²/s as these atmospheres are not always quenched in our grid, which makes a linear interpolation between such models inconsistent. Finally, we fitted a reference pressure (P_ref) at the planetary radius, which gives a total of seven free parameters (Z, C/O, T_int, K_zz, P_ref, Δ_HST, and P_cloud) with flat priors that were fitted for.

Fig. 8 Cloud-free synthetic spectra of our grid of 1D-RCPE models for WASP-107b computed with petitRADTRANS. The top row indicates the dominant opacity sources, with secondary contributors in brackets.

4.2.1 Results of the nominal retrieval: Case 1

Fig. 9 shows the best-fit spectrum of our nominal retrieval (Case 1), summarized in Table D.1. By analysing the order of convergence, we see that the metallicity $\mathrm{Z}={3.77}_{-0.39}^{+0.37}{\mathrm{Z}}_{\odot }$ and $\log {\mathrm{P}}_{\text{cloud}}(\mathrm{b}\mathrm{a}\mathrm{r})=-{3.54}_{-0.06}^{+0.07}$ are constrained the easiest, together with offset parameters P_ref and Δ_HST. Upon closer examination (Fig. B.1), a degeneracy exists between Z, P_ref, and P_cloud due to the dominant presence of H₂O, both due to the amount of spectral features and high VMR, and the strong feature of CO₂ at 4.3 μm. Depending on the metallicity, and thus overall H₂O and CO₂ content, one can vary both P_cloud and P_ref to correctly match the scale height of these features. Note that an increased metallicity affects the spectrum in two ways. First, the amplitude of the many H₂O features (and CO₂) will scale with their abundance, which can be fitted with a vertical offset (controlled by P_ref) to match the top of the feature and P_cloud, which blocks the layers below. Second, an increased metallicity also implies a higher overall mean molecular weight. This effectively decreases the pressure scale height and hence the physical extent of the atmospheric layers through which star light is absorbed. This reduces the amplitudes of all spectral features, including those of H₂O and CO₂, and has a similar effect to fitting a high-altitude cloud deck (low P_cloud). It is the 4 μm SO₂ feature that breaks this degeneracy, despite slightly underestimating the transit depth in this region, as SO₂ is highly sensitive to metallicity. In Sect. 4.3, we further discuss the ability of SO₂ features to constrain metallicity.

The constraint of P_cloud to low pressures is crucial for the retrieval of $\mathrm{C}/\mathrm{O}={0.13}_{-0.02}^{+0.02},{\mathrm{T}}_{\mathrm{i}\mathrm{n}\mathrm{t}}={163}_{-44}^{+62}\mathrm{K}$, and $\log {\mathrm{K}}_{\mathrm{z}\mathrm{z}}={8.68}_{-0.51}^{+0.56}$ (cgs). The most straightforward explanation of such low C/O would be that an extremely carbon-poor atmosphere removes excessive CH₄ that would otherwise dominate the full wavelength range. However, also oxygen-rich species are affected by such low carbon-to-oxygen ratios. Compared to a solar C/O, there is also less CO₂ and more H₂O, which are both dominant contributors in these spectra. Fitting the relative proportion between both H₂O and CO₂ features cannot be done by adjusting the general metallicity. Instead, one must lower the carbon-to-oxygen ratio to downscale the CO₂ feature with respect to the H₂O bands, which explains the preference towards such low C/O value.

So far, we have not discussed the 3.4 μm CH₄ feature that is fitted nicely by our retrieval. Other than a sub-solar C/O, we stated in Sect. 3.2 that a combination of high T_int and high K_zz causes deep quenching and a low CH₄ abundance (Fig. 6). However, we fit our grid to a very low T_int, which leaves a large amount of methane in the atmosphere (see Fig. 2). Fig. 9 also shows the very limited impact of T_int on the spectrum for this combination of Z, C/O and P_cloud. We find no substantial difference (< 20 ppm) across the full prior range of considered intrinsic temperatures, a very peculiar result considering our findings in Sect. 3.2 (Fig. 6). However, by fitting the cloud deck at P_cloud ≃ 10^−3.54 bar, one effectively shields the layers below where the bulk methane is present. One still needs to fit the 3.4 μm methane feature in some way, which is achieved via the vertical mixing parameter K_zz. Although a larger value for K_zz implies deeper quenching, and thus a lower quenched VMR, it has a different effect in the upper atmosphere. There, it competes with photochemical destruction to mix additional CH₄ in the upper layers. Given the high-altitude cloud deck that is retrieved exactly around these layers, the retrieval algorithm prefers models with high K_zz to provide the necessary amount of methane that can fit the 3.4 μm feature. This is also in agreement with a lower T_int, as this sets the bulk amount of methane beneath the cloud deck, which is mixed upwards. Note that this explanation does not contradict our findings of Sect. 3.2, as the situation is drastically affected by the cloud deck. This means that, in this particular case, parameters T_int and K_zz were fitted to maximize the CH₄ feature at 3.4 μm instead of minimize it.

Fig. 9 Best-fit spectrum of Case 1 with molecular contributions (upper panel). The bottom panels show spectra within 1σ, 2σ, and 3σ uncertainty levels of the mean posterior values. The data plotted in the upper panel consists of HST WFC3 (Kreidberg et al. 2018) and JWST NIRCam (Welbanks et al. 2024).

4.2.2 Sensitivity analysis: Cases 2, 3, 4, 5, and 6

We continue with the same retrieval set-up, but now we systematically fix a parameter beforehand and exclude it from the run. In this way, we enhance our understanding of the 1D-RCPE retrieval method and demonstrate how our results are influenced by variations in the specific model parameters. The results of these retrievals are summarized in Table D.1. When assuming Z=10 Z_⊙ (Case 2) based on previous findings (Welbanks et al. 2024; Dyrek et al. 2024), we find a different result compared to our nominal retrieval (Case 1). Due to the strong presence of the CO₂, SO₂, and H₂O features at this metallicity, the cloud deck was placed at a higher altitude (lower P_cloud) to reduce their amplitudes. Note that the reference pressure, P_ref, was adjusted accordingly, being degenerate with P_cloud and metallicity. Our retrieval constrains K_zz to a value close to the lower boundary of the prior (<10^7.5 cm²/s). This is again a peculiar result given the fact that previously we found that the high cloud deck requires a higher K_zz value to mix sufficient amounts of CH₄ into the upper layers to form a visible spectral feature. However, the same upward mixing also affects other species such as H₂O, CO₂, and in particular SO₂. After fixing P_cloud based on the high metallicity, the retrieval minimizes K_zz to further reduce the scale height of these features. We note that the carbon-to-oxygen ratio is, again, fitted based on the relative amplitude of the CO₂ and H₂O bands, which explains the extremely low value. Finally, we note that T_int remains unconstrained, as the combination of high clouds and low K_zz prevents any amount of methane from being visible in the upper layers. The overall fit is of a much lower quality than the nominal model (${\chi }_{\text{red}}^2=2.75$ compared to 2.18 for Case 1). The reason is twofold as (i) no combination of Z, P_ref and P_cloud can sufficiently reduce the 4 μm SO₂ feature and (ii) the low value of K_zz prevents CH₄ from forming a feature at 3.4 μm above the cloud deck.

Running the retrieval with a fixed solar C/O (Case 3), gives $\mathrm{Z}={2.71}_{-0.17}^{+0.22}{\mathrm{Z}}_{\odot },{\mathrm{T}}_{\text{int}}={490}_{-15}^{+7}\mathrm{K}$ and $\log {\mathrm{K}}_{\mathrm{z}\mathrm{z}}={7.18}_{-0.13}^{+0.22}(\mathrm{c}\mathrm{g}\mathrm{s})$. By fixing the C/O beforehand, and thus the relative amplitude between the H₂O and CO₂ features, a slightly lower metallicity is preferred to fit these features. This is a consequence of the nonlinear relation between the VMR of CO₂ and Z, while H₂O does scale linearly with metallicity. Towards higher metallicities, the increase in CO₂ is greater than H₂O, requiring a much lower C/O to fit the resulting feature amplitudes. Note that this results in a poor fit of the SO₂ feature, as there is none at such low metallicity in these spectra. The slightly lower value for Z results in a deeper cloud deck (higher P_cloud), as there is less of a need for these high clouds to reduce the amplitudes of the spectral features. These relative deep clouds now expose more layers containing a substantial CH₄ content, which in turn explains the extremely high T_int and low K_zz. We repeat that, although a low K_zz results in more shallow quenching (Fig. 6), there is less methane mixed in the upper (visible) layers (see Fig. 2).

Another exercise is to fix T_int=500 K (Case 4), a parameter that is typically constrained with great difficulty as it only affects the 3.4 μm CH₄ feature. The retrieval first finds a solution between the metallicity, P_ref, and P_cloud, after which it converges to a sub-solar C/O and low K_zz. Interestingly, both Z and C/O are much less tightly constrained than before, encompassing the values of Case 1 within the 1 to 2 σ level. As a high T_int does not affect the VMR of H₂O, CO₂, or SO₂, the only difference stems from starting with a fixed T_int that affects the VMR of CH₄. When we fix log K_zz=10 (cgs) in Case 5, we find results that closely resemble our nominal model with no fixed parameters (Case 1). Do note the exception that a slightly higher intrinsic temperature is fitted for, which removes additional CH₄, as more material gets mixed in the upper observable layers as a result of the higher K_zz. The BIC values for Case 4, and especially Case 5, are close to that of the nominal model (Case 1). This confirms that T_int and K_zz are less important for the retrieval outcome compared to Z and C/O, as they are more difficult to constrain.

Finally, we fixed the cloud deck to P_cloud=10⁻³ bar (Case 6), which is deeper in the atmosphere than has been previously fitted for. As usual, this directly affects the metallicity (and P_ref), which was now fitted to a relatively low value of $\mathrm{Z}={1.74}_{-0.08}^{+0.08}{\mathrm{Z}}_{\odot }$ as there was no high cloud deck to reduce the scale height of the spectral features. We find a slightly higher C/O compared to other runs, which is a result of fitting a lower Z that affects the relative VMR of CO₂ with respect to H₂O. Again, this lower Z affects the quality of fit as the SO₂ bump at 4 μm is not fitted well. The deep cloud deck fails to block the CH₄−rich layers. Therefore, a combination of high T_int and low K_zz minimizes the CH₄ opacity contribution to the spectrum, explaining the final results.

4.2.3 A wavelength-dependent opacity slope: Case 7

A consistent result among all runs is the very low C/O value. The retrievals achieve these values by trying to reduce the CO₂ feature (4.3 μm) in amplitude with respect to the H₂O band (2.5–3.0 μm) or vice versa. Rayleigh scattering can introduce a slope towards the optical part of the transit spectrum, effectively altering the fit for the HST data (and thus Δ_HST). Given the detection of larger silicate droplets in WASP-107b (Dyrek et al. 2024), a non-Rayleigh type of short-wavelength absorption could be present that follows a much more shallow wavelengthdependence. In Case 7, we fit for a free slope and find ${\gamma }_{\text{scatt}}=-{0.75}_{-0.09}^{0.09}$, which is more indicative of Mie scattering caused by mid-sized to larger aerosols. The effect on C/O is minimal, however, as we retrieve a similar value as before. We do obtain a better fit compared to the nominal model (Case 1), as well as a better BIC, indicating that this retrieval manages to explain the data better. Note that P_cloud loses its meaning as the opacity source with the shallow slope now replaces the grey cloud deck (P_cloud). Therefore, in future simulations with a fitted γ_scatt, we do not fit for P_cloud.

Fig. 10 Best-fit spectrum for Case 10. The data plotted in the bottom panel consists of HST WFC3 (Kreidberg et al. 2018), JWST NIRCam (Welbanks et al. 2024), and JWST MIRI/LRS (Dyrek et al. 2024).

4.2.4 Information content of HST and JWST NIRCam data: Cases 8 and 9

Finally, we briefly consider different combinations of datasets in these retrievals with grey clouds, excluding the JWST MIRI data of WASP-107b for now. A fit to the HST dataset alone (Case 8) only constrains a high cloud deck, as no features other than of H₂O are visible. In particular, the metallicity remains unconstrained, which provides further evidence that features from only H₂O are not sufficient to constrain Z. Furthermore, fitting to the NIRCam dataset separately (Case 9) gives similar results as to fitting it in combination with HST data (Case 1). We conclude that, as expected (e.g. Edwards et al. 2023), the retrieval needs a feature-rich spectrum such as JWST NIRCam to constrain atmospheric properties.

4.3 Retrievals with non-grey clouds on transit spectra up to 12.2 μm

Observations with the JWST MIRI instrument (LRS) provide us with additional wavelength coverage up to ∼12.2 μm. Given the additional opacity source between 8 and 10 μm, stemming from (silicate) clouds (Dyrek et al. 2024; Welbanks et al. 2024), we ran retrievals with non-grey cloud parametrizations described in Sect. 2.6. We included uniform priors that span the entire range of the four grid parameters (Z, C/O, T_int, K_zz) with the exception of log K_zz=6 (cgs) as we only want to include quenched atmospheres. Next to an offset with HST data (Δ_HST), we also fitted an offset with the NIRCam data (Δ_NIR) in reference to the MIRI data. Together with the model offset (P_ref), this makes a total of seven free parameters without the cloud treatment.

4.3.1 Results with Gaussian clouds: Cases 10, 11, 12, 13, and 14

We started by adopting the cloud parametrization of Welbanks et al. (2024), which introduce three more fitting parameters (κ_opac, λ₀, ω) as we fixed ξ = 0 after initial testing to reduce the computational cost. Additionally, we included a wavelength-dependent opacity slope as described in Eq. (2), which introduces two additional fitting parameters to the retrieval (κ_scatt, γ_scatt). Note that we do not fit for a grey cloud deck (P_cloud) if we consider γ_scatt in the retrieval.

Figs. 10 and B.2 show the best fit spectrum marginal posterior distributions of this retrieval (Case 10). Overall, we achieve a good fit for all datasets with the exception of the noisy region above 10 μm. At 4 μm, we find that our best-fit model contains mostly SO₂ opacity, although it does not reproduce the narrow peak in the data completely. As is established in Sect. 4.2, the metallicity is degenerate with the model offset (P_ref) and cloud opacity (e.g. P_cloud) when fitting to the H₂O and CO₂ features in the transit spectrum below 5 μm. Only the 4 μm feature of SO₂ helps constrain the metallicity as higher values for Z (and thus higher VMR of SO₂) would result in a SO₂ feature with too high amplitude. Interestingly, the inclusion of two broad SO₂ features around 7–9 μm does not change the retrieved metallicity ($\mathrm{Z}={3.78}_{-0.53}^{+0.53}{\mathrm{Z}}_{\odot }$ significantly compared to previous retrievals with grey clouds ($\mathrm{Z}={3.77}_{-0.39}^{+0.37}{\mathrm{Z}}_{\odot }$, see Table D.1). This could indicate that any of the SO₂ features are a robust indicator of atmospheric metallicity, or that the Gaussian cloud feature is degenerate with the SO₂ features at 7–9 μm, and thus metallicity is constrained only on the 4 μm SO₂ feature.

If SO₂ is a strong indicator of the metallicity in a photochemically active atmosphere, its features in the MIRI data should cause the retrieval to converge to a similar value for Z. Indeed, re-running the retrieval without HST and NIRCAM data (Case 14, with P_cloud instead of γ_scatt) results in $\mathrm{Z}={3.16}_{-0.43}^{+0.78}{\mathrm{Z}}_{\odot }$, which is within the 1 σ level as before (Case 10). Further investigation confirms that the quality of fit decreases with a fixed higher metallicity (see BIC of Case 11 with Z=10 Z_⊙ in Table D.2). We conclude that the SO₂ features, either at 4 μm or 7–9 μm, are robust features to constrain the metallicity. In Sect. 5, we further discuss their dependence on model assumptions such as the temperature in the upper layers.

The inclusion of MIRI data introduces features of two oxygen-bearing species (H₂O and SO₂) and (a lack of) one carbon-bearing molecule (CH₄), which can further impact the C/O. We find that the nominal retrieval (Case 10) has difficulties constraining C/O (next to T_int and K_zz) compared to other parameters such as Z, P_ref, and the cloud properties. Similar to our findings before (Sect. 4.2), we find that the easy constraints on cloud (and scattering) opacity affect the further convergence of C/O. As the high-altitude layer of aerosols shields the layers with substantial CH₄, and thus the formation of any features in the spectrum, the $\mathrm{C}/\mathrm{O}\left({0.20}_{-0.02}^{+0.03}\right)$ is constrained by the relative amplitude of the 4.3 μm CO₂ feature with respect to features of oxygen species (e.g. H₂O and SO₂ in NIRCAM). Within the MIRI data, lowering the C/O effectively reduces the amplitude of both SO₂ and H₂O feature as no carbon-species has a strong feature. Therefore, we find that the C/O is still mainly controlled by the 4.3 μm CO₂ feature. To further assess the sensitivity of C/O to the scattering treatment, we run the same retrieval without scattering and find a marginally lower C/O of ${0.15}_{-0.02}^{+0.02}$ (Case 13, Table D.2). Furthermore, a simulation with a fixed solar C/O gives a satisfactory fit (see BIC of Case 12 in Table D.2), indicating that an adjusted low-wavelength opacity slope can reduce the need to fit a low C/O. We carefully conclude that the above indicates a degeneracy between C/O and the wavelength-dependent opacity slope.

Parameters T_int and K_zz are least constrained by the retrieval. We find a relatively high intrinsic temperature of ${369}_{-70}^{+60}\mathrm{K}$ and intermediate log K_zz of ${8.25}_{-0.68}^{+0.64}(\mathrm{c}\mathrm{g}\mathrm{s})$. As was previously established, a high-altitude cloud deck shields the CH₄−rich layers and prevents it from forming features in the observable spectrum. Once again, this hypothesis is confirmed by finding no substantial difference (< 15 ppm) between model spectra across the full range of T_int, while otherwise adopting the best-fit values. The CH₄ feature at 3.4 μm is fitted by increasing K_zz so that additional CH₄ is mixed upwards and reaches above the cloud deck. While T_int mostly affects CH₄ features in the transmission spectrum, we do note that K_zz also affects the abundances of SO₂ and CO₂ (e.g. Fig. 2 and Fig. A.1), and therefore also the scale height of their spectral features.

Finally, we note that none of the above fits contain a significant contribution from NH₃, a species detected by Welbanks et al. (2024) at 3 μm (and 10.5 μm to a lesser extent) in WASP107b. Although we do not obtain a good fit overall above 10 μm with our models, the region around 3 μm is well fitted by mostly H₂O opacity instead of NH₃. As is discussed in Sect. 3.3 and shown in Fig. 2, NH₃ is most abundant in atmospheric models with low T_int. Similar to CH₄, high altitude clouds could obscure NH₃−rich layers and limit the amplitude of any potential spectral features. Meanwhile, at a high K_zz additional material could be mixed upwards above the cloud deck and become visible in the transit spectrum. However, we do not find substantial opacity contributions of NH₃ in cloud-free models (Fig. 8). To further investigate the impact of NH₃, we run a retrieval with Gaussian clouds and a wavelength-dependent opacity slope (Similar to Case 10) without NH₃ opacity in the radiative transfer computation. We find that this does not affect the retrieval outcome and conclude that NH₃ opacity does not play a role in any of the above retrievals.

4.3.2 Retrievals with silicate condensate clouds: Cases 15, 16, and 17

Given the detection of silicate condensates in the MIRI data (Dyrek et al. 2024), we run retrievals that use optical properties of SiO, SiO₂, or MgSiO₃ to explore their effect on other fitting parameters. As is detailed in Sect. 2.6, this introduces four more parameters ${\mathrm{X}}_{\text{cloud}},{\mathrm{P}}_{\text{cloud}},{\mathrm{f}}_{\text{sed}},{\mathrm{K}}_{\mathrm{z}\mathrm{z}}^{\text{(cloud)}}$) as we fix x_σ=−1. Next to the four grid dimensions (Z, C/O, T_int, K_zz) and three offsets (P_ref, Δ_HST, Δ_NIR), this totals 11 fitting parameters. Note that the implementation of particle size and mass fraction distribution removes the need for separate parametrized scattering. As expected, we retrieve metallicities and carbon-to-oxygen ratios that are in line with previous runs, confirming that both these parameters are not sensitive to the way we choose to parametrization clouds in our set-up. In contrast to our results with Gaussian clouds, we now find lower intrinsic temperature (≲ 250 K). In line with our previous findings, a low T_int is caused by the altitude at which the aerosols become optically thick and whether or not this exposes the CH₄−rich layers. Another notable difference with the Gaussian cloud parametrized retrievals is the quality of fit, which is much lower in this case. Closer inspection reveals that the JWST MIRI data is not properly fitted for (Fig. C.1), which confirms that a single silicate species cannot explain the broad 8 to 10 μm cloud feature (Dyrek et al. 2024). Although determining which (combination of) silicate cloud species or other condensates can fit the data better is out of the scope of this study, we briefly discuss the debated presence of silicate clouds in Sect. 5.3.

5 Discussion

In this section we discuss the implications of our results for WASP-107b. First, we discuss our findings on the elemental composition (Z and C/O) in Sect. 5.1 before discussing the reason why our retrievals cannot constrain T_int and K_zz properly in Sect. 5.2. Finally, we review the presence of high-altitude clouds that consist of silicates particles in Sect. 5.3.

5.1 Retrieved composition of WASP-107b

We consistently find metallicities of 3–5 Z_⊙ in retrievals that achieve a good fit. As is discussed extensively in Sect. 4, the metallicity is constrained based on the SO₂ features at 4 and 7–9 μm as SO₂ production is dependent on the available sulphur (H₂S) and oxygen reservoirs (H₂O) (see Sect. 3.1). Sulphur dioxide effectively breaks the degeneracy between metallicity (based on features of H₂O, and CO₂), model offset (fitted P_ref at R_pl), and top of the cloud deck (P_cloud). Our derived metallicities are lower than values reported in the literature, which range between 10 and 18 Z_⊙ (Welbanks et al. 2024) and 43 ± 8 Z_⊙ (Sing et al. 2024). As SO₂ is a photochemical product, its VMR profile is determined by (i) the SED emitted by the host star (in particular the UV), and (ii) the temperature structure (Tsai et al. 2023; Dyrek et al. 2024; de Gruijter et al. 2025). We note that Welbanks et al. (2024) use the same SED (HD 85512 from the MUSCLES survey) as we do in this study. Therefore, we attribute part of the difference in derived metallicity to their lower upper-layer temperature (450–550 K) with respect to ours (650–750 K). At lower temperature, the thermochemical conversion from H₂O into OH (Reaction (4)) becomes slower, which suppresses the efficiency of sulphur oxidization (e.g. Reaction (8)). Hence, a higher metallicity is required to produce sufficient SO₂ in their models.

The strong sensitivity of SO₂ to temperature inspires us to revisit an assumption made in the construction of PT profiles (Sect. 2.2). By fixing r_st = 1 (effectively a dayside average PT structure), we computed a profile with a quasi-isothermal upper atmosphere of 650–750 K while setting r_st = 0.5 (planetwide average) would lower this by ∼100 K. To test the effect of this assumption on our results, we ran an additional grid of 900 models with these slightly cooler PT profiles. We then ran retrievals using the same set-ups as Case 1 (Table D.1) and Case 10 (Table D.2); we present the results in Table D.4. We find metallicities of ${3.92}_{-0.36}^{+0.34}{\mathrm{Z}}_{\odot }$ (Case 18) and ${4.49}_{-0.27}^{+0.32}{\mathrm{Z}}_{\odot }$ (Case 19), respectively, which indicates only a marginal increase in retrieved metallicity compared to our initial retrievals with r_st=1. Fig. 11 shows the best-fit spectrum of Case 19 and Case 10, as well as the PT profiles and VMRs of key species. Despite a general decrease in SO₂ in the upper regions, a similar fit can be achieved by leveraging the R_pl–P_ref degeneracy and cloud opacity. This indicates that the predictive power of SO₂ for metallicity indeed depends on the assumed temperature in the upper layers. What further complicates the correct estimation of the upper layer temperature is the fact that observations do not probe a 1D atmosphere and thus we cannot assume a uniform temperature on the terminator region. A slight asymmetry between the evening and morning limb has been detected on WASP-107b, resulting from a ∼100 K difference (Murphy et al. 2024). Limb asymmetry not only affects the temperature and chemistry but changes the scale heights as well, which can bias the observed transit spectrum (Caldas et al. 2019; Falco et al. 2022; Pluriel et al. 2022; Lee et al. 2022; Nixon & Madhusudhan 2022).

Another consistent finding of our retrievals is the carbon depletion with respect to oxygen (C/O), typically in the range of 0.10 to 0.20. We find that this parameter is mainly determined to the relative amplitude of the 4.3 μm CO₂ feature with respect to broad features of oxygen species, in particular H₂O. If the retrieval of this parameter is reliant on a single CO₂ feature, one must question the relatability and investigate the underlying chemistry that produces CO₂. As is detailed in Sect. 3.2, we find that the VMR of CO₂ is mainly set by thermochemistry and is less sensitive to disequilibrium effects in our models of WASP-107b. Retrievals with cooler upper layers (r_st=0.5) prefer models with a slightly lower C/O (Table D.4), indicating an even higher carbon depletion is needed to fit the spectra properly.

Fig. 11 Best-fit PT profile (upper left), VMRs of key species (upper right), and model spectra (bottom) of retrievals Case 19 and Case 10 that have Gaussian clouds and a free wavelength-dependent opacity slope. The grid of models used for Case 19 are approximately 100 K cooler in the isothermal layers (rst=0.5) compared to those used for Case 10(rst = 1). Note that both retrievals converge to a different set of grid parameters (Z, C/O, Tint, Kzz), which affects the VMRs and PT profiles, but a fit of similar quality can be achieved.

5.2 High-altitude aerosols can shield CH₄−rich layers

We retrieve intrinsic temperatures that span the entire range of our priors (100 to 500 K). Together with K_zz, our retrievals struggle to constrain T_int and both parameters often have large uncertainties. We attribute the large spread in T_int values to high-altitude clouds that block the layers below from forming spectral features, like the 3.4 μm CH₄ feature. Additionally, we find that strong mixing can mix material above the cloud deck so that a higher K_zz implies a stronger feature. This is counter-intuitive from a chemical standpoint, which shows that stronger mixing results in deeper quenching of CH₄ and thus lower VMR (Fig. 6). We conclude that the feature at 3.4 μm is fitted by a delicate balance of P_cloud (or equivalent aerosol opacity) and K_zz, while T_int is less important than expected.

Previous studies have suggested that a high intrinsic temperature, in combination with deep quenching by efficient vertical mixing, is responsible for the lack of abundant methane in WASP-107b (Kreidberg et al. 2018; Fortney et al. 2020; Dyrek et al. 2024; Welbanks et al. 2024; Sing et al. 2024). In particular, Welbanks et al. (2024) constrains T_int>345 K and log K_zz = 8.4–9.0 cm²/s, which our retrievals reproduce when a similar set-up is used (Case 10, Table D.2). Although this paper shows that these results can vary with different modelling choices and retrieval set-ups, we can neither confirm nor reject the high internal temperature scenario for WASP-107b. A hot interior is consistent with the highly inflated radius and low density, likely caused by tidal inflation (Leconte et al. 2010; Fortney et al. 2020; Welbanks et al. 2024), and the suggested presence of high-altitude silicates (Dyrek et al. 2024). Instead, we argue that hybrid retrieval approaches must be interpreted with caution and their robustness verified against various assumptions.

The delicate balance between T_int, K_zz, and clouds relates directly to the underlying grid of forward models. If the cloud deck is fitted slightly deeper in the atmosphere, and more CH₄ is exposed to the observable layers, the retrieval tends to maximize T_int and minimize K_zz, while the opposite is true for a slightly higher cloud deck. This behaviour could be a consequence of the disconnect between our treatment of disequilibrium chemistry and post-processing of clouds in the radiative transfer. In particular, it can be expected that the methane homopause is shifted upwards as high clouds act as an additional UV absorber, shielding lower layers from photochemistry and increasing local heating rates (Molaverdikhani et al. 2020; Sánchez-Lavega et al. 2023). Additionally, clouds affect the gas-phase chemistry (Mahapatra et al. 2017; Helling 2019; Kiefer et al. 2024), which in turn also affects the PT structure (Marley et al. 2013; Lee et al. 2024). This calls for a general framework that allows feedback between (kinetic) cloud formation, gas-phase disequilibrium chemistry, and radiative transfer.

Finally, we note that other processes can deplete CH₄ as well. Zonal quenching can be induced by equatorial mixing (Bell et al. 2024), although this effect could be small for WASP-107b (Baeyens et al. 2021). Furthermore, vertical mixing is treated as eddy diffusion, which assumes length scales that are much smaller than the pressure scale height, thereby mimicking diffusion (Tsai et al. 2017). However, given the potentially large interior temperatures and resulting steep temperature-pressure gradients, a large part of WASP-107b’s atmosphere could be dominated by convective energy transport. Therefore, additional large-scale mixing might take place that is not captured by the constant eddy diffusion coefficients (Baeyens et al. 2021).

5.3 Silicate clouds at high altitudes

Dyrek et al. (2024) report a statistical preference for silicate species SiO, SiO₂, and MgSiO₃ with respect to a simpler cloud set-up. Particularly notable is the high altitude at which the silicate cloud deck was retrieved, while condensation of such silicates occurs at high temperatures (e.g. 1200–1500 K) (Lodders & Fegley 1999; Visscher et al. 2010) and thus much deeper in the atmosphere. Given this expected deep condensation, Welbanks et al. (2024) show that vertically extended clouds are required to become visible in the transit spectrum (f_sed < 0.5). However, the resulting small particles give an unsatisfactory fit to the observations. Dyrek et al. (2024) find a good fit with large silicate particles (f_sed>3) so that the question remains how such a compact cloud deck of large aerosols can be suspended at high altitudes, far above the expected condensation region. One possibility is that efficient vertical mixing can loft particles to high altitudes continuously before they grow large and settle down. Given the intermediate value for f_sed that Dyrek et al. (2024) find, atmospheres with high K_zz would require particles to grow large quickly and form a relative compact cloud at such high altitudes. WASP-107b’s low gravity could facilitate this slow settling and maintain high-altitude silicate clouds.

Although the focus of this study is to assess the robustness of retrieval outcomes against different assumptions, including the cloud parametrization, we briefly put our results in context of the above. In retrievals with individual silicate cloud species, we find f_sed values roughly between 1.50 and 3, moderate values for ${\mathrm{K}}_{\mathrm{z}\mathrm{z}}^{\text{(cloud)}}$, and relatively high altitudes of the cloud pressure (P_cloud). This supports the case of a high-altitude thin cloud deck composed of moderate-sized particles, although this challenges the hypothesis of strong vertical mixing that replenishes these layers and reduces the sedimentation efficiency. When fitting with Gaussian clouds and a wavelength-dependent opacity slope (Welbanks et al. 2024), we retrieve ${\gamma }_{\text{scatt}}=-{0.83}_{-0.09}^{+0.09}$ (Case 10). The absence of a steeper scattering slope (e.g. γ_scatt = −4) indicates that sub-micron particles are not dominant in the atmosphere. Instead, the fitted shallow slope resembles more Mie-type scattering, caused by particles with comparable sizes to the photon wavelength, supporting larger grains of one to several micron.

6 Summary and conclusion

In this paper, we have conducted an in-depth study on the performance and reliability of 1D-RCPE retrievals on transmission spectra of WASP-107b. We have built a grid of 900 atmospheric models that take into account radiative heating, convection, vertical mixing, and thermo- and photochemistry and explore a range of different metallicities, carbon-to-oxygen ratios, intrinsic temperatures, and eddy diffusion coefficients. We have analysed the behaviour of WASP-107b’s gas-phase chemistry against different parameters and mapped important chemical pathways for key species in the atmosphere. We have subsequently used a nested sampling algorithm to perform Bayesian inference on archival datasets of WASP-107b and assessed the robustness and reliability of our outcomes against different cloud treatments.

Below, we summarize our main findings on the chemistry in our models of WASP-107b. The main formation pathway of SO₂ in WASP-107b is through oxidization of atomic sulphur or H₂S+2H₂O → SO₂+3H₂, with variations at a high atmospheric metallicity involving oxidization of di-atomic sulphur (S₂) and/or direct combination of two SO compounds into SO₂. Hydroxide radicals are available due to hydrogen abstraction and photolysis of H₂O. Atomic sulphur is typically liberated by hydrogen abstraction of H₂S. The production of SO₂ is therefore also sensitive to atmospheric metallicity. Furthermore, we explain how photolysis of non-hydrogen-bearing species such as N₂, in the absence of other photon absorbers, can interact with the H₂ reservoir to enable additional SO₂ production. Regarding CH₄, we reproduce the well-known CH₄−CO conversion pathway via hydroxy-methyl (CH₂OH) that affects the quenching pressure and VMRs deep in the atmosphere. We find a negative correlation between the quenched CH₄ content and atmospheric metallicity due to the additional heating of the deep layers. At high T_int, we find that strong mixing (K_zz) results in deep quenching and decreases the quenched VMR of CH₄, while the opposite is true for models with low T_int. Finally, we find that NH₃ quenching occurs much deeper in the atmosphere and is less sensitive to K_zz due to the steep temperature gradient around the quenching point.

Below, we summarize our main findings of the 1D-RCPE retrievals:

• We find that SO₂ features at 4 and 7–9 μm are most informative to determine the atmospheric metallicity, which is consistently retrieved to be between 3 and 5 Z_⊙ in retrievals that achieve a good fit. This is the case when we consider JWST NIRCam data separately, for JWST MIRI/LRS data if the 10 μm silicate cloud feature is taken into account, or a combination of both. We note that the amplitude of the SO₂ features is, next to the stellar SED, dependent on the assumed temperature in the upper layers as this affects the available hydroxide that oxidizes Sulphur;

• In the absence of strong SO₂ features, we find a degeneracy between metallicity, the fitted reference pressure, P_ref, and the altitude of the clouds (P_cloud). This is because the latter two parameters can control the feature amplitude regardless of the VMRs of dominant absorbers such as H₂O and CO₂, which scale with Z;

• The carbon-to-oxygen ratio (C/O) was fitted to sub-solar values (≲ 0.20) based on the relative amplitude of the large 4.3 μm CO₂ feature with respect to broad bands of H₂O (and SO₂), and much less to the 3.4 μm CH₄ feature as the latter is often shielded by high-altitude clouds. This sensitivity highlights the need to properly model the VMR profile of CO₂, which is sensitive to thermochemistry for WASP-107b and disequilibrium chemistry in some other cases;

• In general, we find a slightly higher C/O when a wavelength-dependent opacity slope is included that mimics scattering towards the optical transit spectrum. The main reason for this is that such opacity increases the transit depth of the broad H₂O band around 2.5–3.0 μm with respect to the 4.3 μm CO₂ feature. Therefore, less carbon depletion is necessary to achieve a good fit to the NIRCam data. We stress that this dependence is weak and only changes our results by less than 1 σ with respect to a nominal model without such a sloped opacity source;

• We find that the retrieval has difficulty constraining T_int and K_zz. Depending on the fitted altitude of the aerosol opacity (e.g. by P_cloud, κ_scatt, γ_scatt), T_int can even remain unconstrained. The reason for this is that high-altitude aerosols shield the CH₄ layers, preventing models with low T_int from forming significant methane features in the spectrum other than the one at 3.4 μm. In retrievals with deeper clouds, a high intrinsic temperature is retrieved. This suggests that exoplanet atmospheres with high-altitude clouds prevent us from constraining the intrinsic temperature. The delicate balance between the methane homopause and cloud altitude in our models highlights the need for a more consistent treatment of cloud formation and disequilibrium chemistry;

• Our retrievals fit K_zz mainly based on the 3.4 μm methane feature. In the case of high-altitude clouds, a higher K_zz is preferred as it can mix additional material above the cloud deck. Therefore, we find that a high value for K_zz results in a larger amplitude feature. This is counter-intuitive as cloud-free chemistry models show that strong mixing quenches the atmosphere more deeply, resulting in overall less CH₄ in the quenched layers;

• We find no significant contribution of NH₃ at 3 or 10.5 μm in our cloudy model spectra, although a detection has been reported in WASP-107b. While high-altitude clouds could shield the NH₃−rich layers, we do not find a substantial contribution of NH₃ in cloud-free spectra either.

Our study shows that hybrid retrievals can be a useful tool for analysing high-quality data. However, the use of forward models introduces additional complexity and necessitates a thorough analysis before they can be incorporated into a fitting routine.

Acknowledgements

We thank Paul Mollière for assisting with the use of petitRADTRANS, and Shang-Min Tsai for insightful discussions related VULCAN and our results. We also thank Mats Esseldeurs for his valuable advice throughout the development of this research. Finally, we thank the anonymous referee for their insightful comments that highly improved this work. T.K. and L.D. acknowledge funding from the KU Leuven Interdisciplinary Grant ESCHER (IDN/19/028). L.H. and L.D. acknowledge funding from the European Union H2020-MSCA-ITN-2019 under Grant no. 860470 (CHAMELEON). R.B. acknowledges support from the FNWI Origins investment incentive. K.H. acknowledges funding from the BELSPO FED-tWIN research program STELLA (Prf-2021-022) and the FWO research grant G014425N. V.C. thanks the Belgian F.R.S.-FNRS, and the Belgian Federal Science Policy Office (BELSPO) for the provision of financial support in the framework of the PRODEX Programme of the European Space Agency (ESA) under contract number 4000142531.

Appendix A Additional VMR profiles

Fig. A.1 Same as Fig. 2 for HCN, CO2, C2H2, SO, and H2S.

Fig. A.2 Same as Fig. 2 for H2O, OH, H, CO, and H2.

Fig. A.3 VMR profiles of CO and CH4 in our grid of disequilibrium chemistry models of WASP-107b.

Fig. A.4 VMR profiles of N2 and NH3 in our grid of disequilibrium chemistry models of WASP-107b.

Appendix B Corner plots of Cases 1 and 10

Fig. B.1 Corner plot for Case 1.

Fig. B.2 Corner plot for Case 10.

Appendix C Best-fit spectra of Cases 15, 16, and 17

Fig. C.1 Best-fit spectra for Cases 15, 16, and 17 in which silicate condensates were fitted to the observations.

Appendix D Tabulated retrieval outcomes

References

All Tables

All Figures

	Fig. 1 Selection of the 180 computed PT profiles with PICASO for C/O=0.46. The line transparency corresponds to intrinsic temperature, from a full line for Tint=100 K to almost transparent for Tint=500 K. The shaded grey region indicates layers where the radiative-convective boundaries are located in all 180 PT profiles.
In the text

	Fig. 2 VMR profiles of key species CH4, SO2, and NH3 in our grid of disequilibrium chemistry models of WASP-107b. Chemical equilibrium is plotted with thin, dashed lines. The nominal model of Z=10Z⊙, C/O=0.46, Tint=400 K, and log Kzz=9(cm2/s) was chosen to be in agreement with previous studies (Dyrek et al. 2024; Welbanks et al. 2024).
In the text

	Fig. 3 Dominant reactions (by fastest reaction rate) that involve SO and SO2, for three different combinations of metallicity and eddy diffusion coefficient. Forward rates are shown by solid lines; backward rates by dotted lines. The other values of these simulations are Tint=400 K and C/O=0.46.
In the text

	Fig. 4 Dominant reactions (by fastest reaction rate) that produce/destroy S, for three different metallicities. Forward rates are shown by solid lines; backward rates by dotted lines. The other values of the simulations are Tint=400 K, C/O=0.46 and log Kzz=9(cm2/s).
In the text

Fig. 5 VMR of SO2 and dominant H-producing reactions for a nominal model (Z=10 Z⊙, C/O=0.46, Tint=400 K, log Kzz=9 cm2/s) that includes all available photodissociation cross-sections (upper left panel), and three limiting cases where only the photolysis of one particular source is considered (N2, H2O., or NH3), copying the experiment of Dyrek et al. (2024). Panels that show thermochemical reaction rates show both their forward (solid lines) and backward (dashed) rates.

In the text

	Fig. 6 Quenching VMR (left panel) and estimated quenching pressure (right panel) of CH4 versus intrinsic temperature Tint for various eddy diffusion coefficients in low (Z=1 Z⊙) and high (Z=10 Z⊙) metallicity atmospheres. The models with log Kzz = 6 (cgs) are not included as not all atmospheres are quenched at such weak vertical mixing.
In the text

	Fig. 7 Same as Fig. 6, now for NH3.
In the text

	Fig. 8 Cloud-free synthetic spectra of our grid of 1D-RCPE models for WASP-107b computed with petitRADTRANS. The top row indicates the dominant opacity sources, with secondary contributors in brackets.
In the text

	Fig. 9 Best-fit spectrum of Case 1 with molecular contributions (upper panel). The bottom panels show spectra within 1σ, 2σ, and 3σ uncertainty levels of the mean posterior values. The data plotted in the upper panel consists of HST WFC3 (Kreidberg et al. 2018) and JWST NIRCam (Welbanks et al. 2024).
In the text

	Fig. 10 Best-fit spectrum for Case 10. The data plotted in the bottom panel consists of HST WFC3 (Kreidberg et al. 2018), JWST NIRCam (Welbanks et al. 2024), and JWST MIRI/LRS (Dyrek et al. 2024).
In the text

Fig. 11 Best-fit PT profile (upper left), VMRs of key species (upper right), and model spectra (bottom) of retrievals Case 19 and Case 10 that have Gaussian clouds and a free wavelength-dependent opacity slope. The grid of models used for Case 19 are approximately 100 K cooler in the isothermal layers (rst=0.5) compared to those used for Case 10(rst = 1). Note that both retrievals converge to a different set of grid parameters (Z, C/O, Tint, Kzz), which affects the VMRs and PT profiles, but a fit of similar quality can be achieved.

In the text

	Fig. A.1 Same as Fig. 2 for HCN, CO2, C2H2, SO, and H2S.
In the text

	Fig. A.2 Same as Fig. 2 for H2O, OH, H, CO, and H2.
In the text

	Fig. A.3 VMR profiles of CO and CH4 in our grid of disequilibrium chemistry models of WASP-107b.
In the text

	Fig. A.4 VMR profiles of N2 and NH3 in our grid of disequilibrium chemistry models of WASP-107b.
In the text

	Fig. B.1 Corner plot for Case 1.
In the text

	Fig. B.2 Corner plot for Case 10.
In the text

	Fig. C.1 Best-fit spectra for Cases 15, 16, and 17 in which silicate condensates were fitted to the observations.
In the text