CAMEMBERT: A Mini-Neptunes GCM Intercomparison, Protocol Version 1.0. A CUISINES Model Intercomparison Project

With an increased focus on the observing and modelling of mini-Neptunes, there comes a need to better understand the tools we use to model their atmospheres. In this paper, we present the protocol for the CAMEMBERT (Comparing Atmospheric Models of Extrasolar Mini-neptunes Building and Envisioning Retrievals and Transits) project, an intercomparison of general circulation models (GCMs) used by the exoplanetary science community to simulate the atmospheres of mini-Neptunes. We focus on two targets well studied both observationally and theoretically with planned JWST Cycle 1 observations: the warm GJ~1214b and the cooler K2-18b. For each target, we consider a temperature-forced case, a clear sky dual-grey radiative transfer case, and a clear sky multi band radiative transfer case, covering a range of complexities and configurations where we know differences exist between GCMs in the literature. This paper presents all the details necessary to participate in the intercomparison, with the intention of presenting the results in future papers. Currently, there are eight GCMs participating (ExoCAM, Exo-FMS, FMS PCM, Generic PCM, MITgcm, RM-GCM, THOR, and the UM), and membership in the project remains open. Those interested in participating are invited to contact the authors.


INTRODUCTION
Super-Earths and mini-Neptunes represent a demarcation, albeit a nebulous one, between the giant planets with their thick atmospheres dominated by hydrogen and helium and the terrestrial planets with thinner secondary atmospheres (Lopez & Fortney 2014). The planets that have retained their hydrogen-dominated atmospheres are believed to have undergone runaway accretion during their formation in order to accumulate a thick atmosphere (Pollack et al. 1996;Lee et al. 2014) but have also been able to retain some or all of that atmosphere in the presence of irradiative evaporation (Owen & Jackson 2012).
Simulating these planets with general circulation models (GCMs) presents a unique set of challenges not necessarily seen in the Earth sciences community. The primitive equations, which assume hydrostatic balance, a thin atmosphere, and a constant gravitational acceleration with height, may not yield accurate results for cases where the thickness of the modelled atmosphere becomes significant relative to the radius of the planet, limiting their applicability to planets with thick atmospheres (see Tokano 2013, Tort et al. 2015, and Mayne et al. 2019 for a discussions related to Venus and Titan, Earth, and mini-Neptunes, respectively). It has also been argued that simulations of mini-Neptunes specifically may have extremely long convergence times, potentially of 50,000 Earth days of model time or more (Wang & Wordsworth 2020), which raises questions about the accuracy of simulations of only a few thousand days, as are common in the exoplanetary modelling community (e.g., Charnay et al. 2015;Mayne et al. 2017Mayne et al. , 2019. These issues motivate a better understanding of the tools we use to model these planets. While intercomparison studies have been somewhat common in the Earth sciences community (see e.g. Eyring et al. 2016;Pincus et al. 2016;Haarsma et al. 2016;Ullrich et al. 2017), it has not been until recently that intercomparisons have been done with a focus on exoplanetary targets. Although not an intercomparison of multiple GCMs, Heng et al. (2011) performed a comparison of the spectral and finite-difference dynamical cores in the Geophysical Fluid Dynamics Laboratory (GFDL) Princeton Flexible Modelling System (FMS) using the hot Jupiter HD209458 b as a test case. The first true intercomparison of GCMs used in exoplanetary science, Polichtchouk et al. (2014), looked at highly idealised configurations -a steady state jet, a baroclinic wave, and diabatic forcing -for five GCMs with the in-tention of better understanding their respective dynamical cores. Increasing the complexity, Yang et al. (2019) compared GCMs in the context of terrestrial planets, specifically focusing on the cases of an Earth-like planet orbiting a G-star and a tidally-locked planet around an M-star. More recently, the TRAPPIST-1 Habitable Atmosphere Intercomparison (THAI; Fauchez et al. 2020) compared GCM models of TRAPPIST-1e, investigating the dynamics Sergeev et al. 2021) and the synthetic observations ) resulting from the simulations. In the final paper, they propose a "GCM uncertainty error bar" of ∼ 50% when interpreting transmission spectra with the uncertainty explained mostly by the cloud differences found between GCMs. The ability to provide this form of context to synthetic observations highlights the value of projects like THAI.
Based on the success of the THAI project, we propose here the CAMEMBERT (Comparing Atmospheric Models of Extrasolar Mini-neptunes Building and Envisioning Retrievals and Transits) intercomparison of GCMs modelling mini-Neptunes under the umbrella of the CUISINES (Climates Using Interactive Suites of Intercomparisons Nested for Exoplanet Studies) framework for intercomparisons for exoplanets (Sohl et al. in prep.). The broad objectives of CUISINES are twofold. First, it provides a meta-framework to quantify, and potentially mitigate, differences between exoplanet model outputs. Second, it aims to assess how these output differences affect the synthetic observations that are used to predict the detectability of atmospheric constituents and to interpret data from ground and space telescopes. With the increased focus on mini-Neptunes with the launch of TESS (e.g., Trifonov et al. 2019;Lacedelli et al. 2021;Burt et al. 2021), CHEOPS (e.g., Bonfanti et al. 2021;Leleu et al. 2021), JWST (e.g., Greene et al. 2017;Bean et al. 2021;Hu & Damiano 2021) and in anticipation of increased efforts to model these planets, we believe the timing is appropriate for an intercomparison of GCMs modelling mini-Neptunes to provide a foundational understanding of how our models behave and how our model choices may impact the interpretation of observations.
In this first paper, we outline the protocol for the CAMEMBERT model intercomparison project (MIP), providing both the motivations for the test cases as well as sufficient details to reproduce them, in the hopes that the results of this intercomparison can be used as a calibration for future GCMs. While tests of the protocol were run using the UM and Exo-FMS to gauge the viability of the protocol, the results from these tests are not presented here, with the intention being to present the results from all participating GCMs in one or more followup papers. The outline of this paper is as follows: an overview of the GCMs currently participating in the intercomparison and a discussion of the target planets are found is Section 2. The protocol and associated test cases are described in Section 3. The outputs and diagnostics are found in Section 4 and a final summary and discussion is found in Section 5.

THE PARTICIPATING GCMS AND THE TARGET PLANETS
With the goal of the intercomparison being to understand the differences between GCMs, the progression of simulations can be divided up into two stages: first, a simple investigation of the dynamics and the dynamical cores at the heart of each of the GCMs. GCM simulations of GJ 1214b, for example, show up to three zonal jets, with differing amplitudes depending on the study (Menou 2012;Kataria et al. 2014;Charnay et al. 2015;Drummond et al. 2018;Mayne et al. 2019;Wang & Wordsworth 2020). In understanding the origins of these differences, it is essential to be able to disentangle differences between dynamical cores and differences between radiative transfer schemes. Once this baseline understanding is established, we progress to a second stage where we compare models with radiative transfer.
The GCMs currently participating, based on an expressed interest at the BUFFET workshop (https:// nexss.info/buffet-registration/) and during the protocol development process, are listed in Table 1. Although these are the participants at the time of the publication of this protocol, participation remains open and other teams are welcome to join.

Choice of Targets
As mini-Neptunes encompass a class of planets with a wide range of orbital and planetary parameters, we opt to focus on a warm, close-in case -GJ 1214b -and a cooler case -K2-18b. The specific target planets were selected based on the existence of previous GCM modelling efforts so as to reduce the barriers to participation as well as the existence of past and planned observations. GJ 1214b is the first mini-Neptune discovered (Charbonneau et al. 2009) and is the archetypal warm mini-Neptune. Due to its proximity to its host star, it is expected to be tidally locked and possess a constant warm dayside and cool nightside, with the possibility of multiple zonal jets, as has been seen in GCM models by Menou (2012); Kataria et al. (2014); Charnay et al. (2015); Drummond et al. (2018); Mayne et al. (2019); Wang & Wordsworth (2020). Questions about the applicability of the primitive equations in modelling these objects (Mayne et al. 2019) as well as questions about the exact flow structure and the convergence timescale (Wang & Wordsworth 2020) further motivate a mini-Neptune intercomparison for GJ 1214b. In addition, a JWST MIRI/LRS phase curve as well as multiple JWST NIRCam transit observations are planned (Bean et al. 2021;Greene et al. 2017).
The more recently discovered K2-18b (Montet et al. 2015) is a temperate mini-Neptune in the habitable zone of its host star. Its place in the habitable zone combined with a possible detection of water vapour in its atmosphere (Tsiaras et al. 2019;Benneke et al. 2019) makes it a tempting target for both characterisation and observation, although the interpretation of the 1.4 µm signal has been disputed (Bézard et al. 2020;Barclay et al. 2021). Thus far, it has been modelled in 3D by Charnay et al. (2021) and Innes & Pierrehumbert (2021), and there are planned JWST transit observations using NIRSpec (Hu & Damiano 2021;Madhusudhan et al. 2021), MIRI (Madhusudhan et al. 2021), and NIRISS (Madhusudhan et al. 2021). While it is unclear whether or not K2-18b is tidally locked (see e.g., Leconte et al. 2015;Charnay et al. 2021), we assume that it is for simplicity.

PROTOCOL
In this section we outline the simulations associated with the intercomparison and their motivations. In general, we seek to maintain consistency in parameters throughout the protocol to facilitate comparison as complexity is increased. While the various cases were initially envisioned as a progression, we encourage participants to join the cases their GCMs are capable of completing regardless of whether or not their GCMs are capable of completing cases earlier in the series.
For all models, we adopt a lower boundary pressure of 3 × 10 6 Pa. While previous studies have studied higher pressures - Charnay et al. (2015) and Mayne et al. (2019) used 8 × 10 6 Pa and 2 × 10 7 Pa, respectively -to simplify the test cases and limit differences between the primitive and less simplified equations of the dynamics for mini-Neptunes (Mayne et al. 2019), we limit ourselves to 3 × 10 6 Pa. We adopt an upper boundary pressure p ≤ 10 Pa for codes which employ a pressure-based vertical grid. For GCMs which use a height-based vertical grid (e.g., the UM), we require a 3.513 × 10 3 1.732 × 10 3 Specific heat capacity c P (J kg −1 K −1 ) 1.200 × 10 4 6.682 × 10 3 Mean molar mass µ (g mol −1 ) 2.367 4.801 2.17 × 10 4 1.37 × 10 3 Tint (K) 100 90 Case 3  domain height sufficient to include pressures of 10 Pa throughout the simulations; however, this may require varying the domain height on a per-case basis. This will also result in those GCMs using height-based grids having potentially significantly lower pressures on the nightside upper boundary, which may necessitate modifications to ensure stability. Following the lead of the THAI project, we do not place specific requirements on the timestep or grid spacing, instead encouraging participants to adopt parameters they would commonly use for exoplanet studies as requirements for stability may differ between GCMs. For simplicity, we assume friction-free, impermeable boundaries to avoid complicating the tests with boundary-layer friction or mass-exchange, and heat exchange is limited to a fixed internal heat flux with effective temperature T int . We do not, however, exclude forms of dissipation that may be required for numerical stability (e.g., sponge layers, artificial viscosity).
As mini-Neptunes may have enhanced atmospheric metallicities relative to solar (Fortney et al. 2013), we adopt a value of 100× solar for K2-18b and use parameters and profiles consistent with this throughout the intercomparison. While there may be similar motivation to use this value for GJ 1214b, due to the number of previous simulations using solar metallicity and the differences between simulations already shown in the literature, specifically in the number and speed of the zonal jets, we adopt the solar value as it allows us to probe a part of the parameter space where we know differences between GCMs already occur.
Simulations are run for a fixed number of Earth days instead of specifying a convergence condition. The chosen simulation lengths do not ensure convergence in simulations which include the deep atmosphere, this may not be sufficient time for the deep atmosphere to have converged to a steady state (Mayne et al. 2017); however, as we are limited to regions with p ≤ 3 × 10 6 Pa we should not be significantly impacted. While Wang & Wordsworth (2020) have found in their simulations of GJ 1214b that over integration times of 50,000 Earth days or more their model atmospheres transition from two off-equatorial jets to a single equatorial jet, including such long integration times in an intercomparison would likely prove computationally prohibitive and limit participation. We instead focus initially on comparing GCMs over shorter timescales, with hopes of extending the work in the future to look at these longer timescales. For cases 1 and 3, we run the simulations for 4,000 Earth days. For case 2, however, we run each simulation for 10,000 Earth days as the dual grey case is the case investigated by Wang & Wordsworth (2020). Although this does not approach the long integration times of Wang & Wordsworth (2020), an integration time of 10,000 Earth days may be sufficient to understand differences between GCMs as Menou (2012) observed the formation of the central equatorial jet in their 7800 Earth day simulation. It may be inevitable, however, that an understanding of a possible delayed formation of an equatorial jet may have to wait for a follow-up study with fewer participants investigating longer simulation times.
Although mini-Neptunes have the potential for clouds, and GJ 1214b in particular has been shown to have strong signs of clouds or hazes (Kreidberg et al. 2014), we do not include a cloudy benchmark as a part of the protocol as currently there are an insufficient number of GCMs capable of participating. We do hope that followup studies and intercomparisons will be able to include a cloud component as clouds will undoubtedly represent an important constituent of future mini-Neptune models.

Initial Conditions
To initialize our simulations, we use one-dimensional pressure-temperature profiles (see Figure 1) with no initial winds. For K2-18b, we use a profile generated using Exo-REM from Charnay et al. (2021) and for GJ 1214b we use a profile generated using ATMO from Drummond et al. (2018). For all cases, we do not include any initial latitudinal or longitudinal variation. Each of these profiles, along with the profiles for the chemical abundances needed for Case 3, is publicly available as an ASCII text file in the CAMEMBERT repository (see Section 4.2).

Case 1: Temperature Forcing
The first case investigated as a part of the intercomparison is similar to the temperature-forced benchmark of Held & Suarez (1994). The motivation is to compare the robustness of the dynamical cores without also comparing differing radiative transfer schemes. In this case, we opt to use a Newtonian cooling prescription where the temperature T is forced to an equilibrium temperature profile T eq on a given radiative timescale τ rad , where t is time. We adopt a radiative timescale τ rad for a hydrogen-dominated atmosphere (Zhang & Showman 2017) τ rad,H2 (p) = 10 4 s p ≤ 10 2 Pa 10 5/2 p 3/4 s 10 2 Pa < p < 10 6 Pa 10 7 s p ≥ 10 6 Pa , (2) with the specific values of R, c p and µ found in Table 2.
The equilibrium three-dimensional temperature profile T eq is generated from the initial temperature profile T 0 (p) and temperature difference ∆T eq (p) intended to mimic qualitatively the day-night temperature contrast expected from tidally-locked planet, where the longitude is λ and the latitude is φ. The temperature difference ∆T eq is taken to be ∆T eq,max at p ≤ 10 Pa and to decrease linearly with log p until p = 10 6 Pa where is becomes zero. The range of equilibrium temperatures at a given pressure is shown in Figure 1. For GJ 1214b, the contrast parameter ∆T eq,max is chosen to maintain consistency with the temperature forcing tests of Mayne et al. (2019). As published temperature forcing tests do not exist for K2-18b, we instead look to Charnay et al. (2021) which shows a more modest ∼ 50 K temperature contrast at the top of the atmosphere for their 100× solar metallicity atmosphere, motivating the chosen value of ∆T eq,max .
We note that the temperature forced case of K2-18b has presented significant difficulty in terms of numerical stability for the UM and Exo-FMS in tests of the protocol. Rather than remove it, we retain it as a part of the protocol with a note of caution to participants.

Case 2: Grey Radiative Transfer
For the initial investigation of the impact of radiative transfer, we employ a dual band approximation with the shortwave and longwave absorption coefficients given in Table 2 to compute heating rates, with the previous forcing scheme no longer included. The values of κ sw and κ lw have been calculated by fitting the initial profiles to the analytic profiles in Guillot (2010) 1 . As participating GCMs may offer different methods to attenuate the incoming stellar irradiation, simulations in Cases 2 and 3 are to be run using the plane-parallel approximation. We adopt this intermediate step before transitioning to non-grey radiative transfer as previous studies show that disagreements between GCMs may already exist at this point (e.g., Menou 2012; Wang & Wordsworth 2020).

Case 3: Non-Grey Radiative Transfer with Fixed Abundance Profiles
To model atmospheric chemistry, we limit ourselves to H 2 /He-dominated atmospheres with H 2 O, CH 4 , NH 3 , CO and CO 2 as well as H 2 and He collisionally-induced absorption (CIA) and Rayleigh scattering as opacity sources. Volume mixing ratios for each species as a function of gas pressure, taken from the same simulations that generated the initial conditions, are shown in Figure 2 and are provided as a part of the initial conditions archive (see Section 3.1). All participating GCMs are to use these abundance profiles, as this allows GCMs without coupled chemistry solvers to participate.
For the stellar spectra, we use the model stellar spectra from PHOENIX BT-Settl (Allard et al. 2012) which closest matches the target star (see Table 2). These spectra are made available in the CUISINES repository along with the other required inputs. The specific linelists, calculation method, and spectral resolution are left to the individual groups.

OUTPUTS, DIAGNOSTICS, AND ARCHIVING
In this section we outline the procedures for formatting and archiving data. The goal is to standardise the output and storage of data as much as possible to facilitate not only the initial analysis but also future analysis by third parties.

Outputting and Formatting of Data
In order to provide sufficient data for comparison, we output diagnostic fields with a 1000 Earth day frequency 10 9 10 7 10 5 10 3 Volume Mixing Ratio  without any averaging applied to track the evolution of the simulated atmospheres. For the final 1000 Earth days, outputs should be every 50 Earth days for the purposes of averaging. For the comparison of atmospheric dynamics, we require pressure and temperature as well as eastward (u), northward (v), and vertical (w) velocity fields for all cases. For cases 2 and 3, we additionally require longwave and shortwave heating rates as well as the top of atmosphere (TOA) outgoing longwave radiation (OLR), outgoing shortwave radiation (OSR), and incoming shortwave radiation (ISR). As the participating GCMs have different dissipation schemes which may lead to differing results, dissipation rates for all drag/damping/sponge schemes should be output as well. The required outputs are summarised in Table 3.
To facilitate the sharing of data and subsequent analysis, we ask that all GCM outputs be in the netCDF format with data stored in SI units. The metadata associated with each variable should include a description of the variable and the associated units, following the netCDF Climate and Forecasting Metadata Conventions 2 . In the case of vector quantities, the sign conventions should be such that northward, eastward, and upward all constitute positive directions for their respective components. Similarly, data should be stored such that increasing grid indices corresponds to a positive change in direction. Fields should be placed on a single rectangular latitude-longitude grid, with longitudes beginning at λ = 0 • and the poles located at latitudes φ = ±90 • . As the protocol considers only tidallylocked cases, the anti-stellar point should be located at (λ, φ) = (0 • , 0 • ). Data along cyclic coordinates should not appear more than once within the dataset. Along the vertical coordinate, a single pressure, potential temperature, or height grid should be adopted. As a number of the participating GCMs are capable of outputting scalar diagnostics with a higher frequency than other outputs, we ask that those GCMs with this capability provide the total axial angular momentum, kinetic energy, and maximum values of each velocity component with an output frequency they would normally adopt. These data can be provided in a text file at the time of submission. We do not make these data a requirement of the protocol as it would require additional development by a significant number of participants. Instead, for those GCMs that do not provide these data separately, these scalar diagnostics will be derived via post-processing of the provided outputs.

The CAMEMBERT Repositories
The data resulting from the simulations and analysis will be uploaded to the CAMEMBERT permanent repository at https://ckan.emac.gsfc.nasa.gov/ organization/cuisines-camembert by participating scientists. These data will be made available for public access upon the publication of the results. Prepublication access can be requested by contacting the authors. Inputs described in this protocol and scripts related to the analysis of data and production of plots for the publications will be made available on the CAMEMBERT GitHub repository at https://github. com/projectcuisines/camembert. Inputs will be available immediately while scripts to reproduce results will be made publicly available upon the publication of the results.

Simulated Observables
As discussed in Section 2.1, GJ 1214b and K2-18b were chosen in part because there are planned observations with JWST in Cycle 1 (Greene et al. 2017;Bean et al. 2021;Hu & Damiano 2021;Madhusudhan et al. 2021). The permanent repository will host the results of the analysis of the results and the post-processed synthetic observations. Consistent with the other intercomparisons that are a part of CUISINES, we will use the Planetary Spectrum Generator (PSG; Villanueva et al. 2018Villanueva et al. , 2022 to simulate JWST spectra for instruments and modes used in Cycle 1 observations of GJ 1214b and K2-18b, and subsequent cycles when available, using the atmospheric outputs provided by each GCM for each of the three cases.

Environmental Impact
While not related to the accuracy of the GCMs, we add as a part of the protocol the requirement that participating scientists include estimates of power consumption and CO 2 emissions associated with each production run included in the intercomparison. We include this requirement not as point of comparison between GCMs, as the environmental impact will primarily depend on the methods of energy generation in the local power grids. Instead, we include this to highlight the environmental impact of supercomputing and to encourage providers of supercomputing resources to transition to environmentally-sustainable energy sources. These data will be reported in the first results paper.

SUMMARY
In this paper, we have presented the protocol for the CAMEMBERT project which seeks to compare GCMs used by the exoplanetary science community, with models of mini-Neptunes being the primary focus. Two benchmarks were chosen -the warm GJ 1214b and the relatively cooler K2-18b -based on the volume of prior modelling work and observational potential, and a series of simulations of increasing complexity are described to calibrate and compare the participating GCMs, with all of the requisite parameters provided here and in the CUISINES repository. Membership in CAMEMBERT remains open, and other groups interested in participating are invited to contact the authors. Collaboration meetings will be held in 2022 and beyond, and the aim is to present the results of the intercomparison in one or more follow-up papers. It is hoped that the results from this intercomparison will provide a strong foundation for follow-up studies exploring elements not including in this initial protocol such as chemistry, clouds, and convergence timescales with any model differences due to these elements being more easily isolated and interpreted. As CUISINES brings together researchers employing a diverse range of tools and approaches, we anticipate that it will act to create new collaborations and stimulate progress in understanding exoplanets. For CAMEMBERT specifically, connections with MALBEC (Modeling Atmospheric Lines By the Exoplanet Community, an intercomparison of radiative transfer codes) and PIE (Photochemical model Intercomparison for Exoplanet science, an intercomparison of 1D photochemistry codes) will begin in 2022 with the goal of comparing results and better informing our own models with the insights from the other projects.

ACKNOWLEDGEMENTS
CAMEMBERT belongs to the CUISINES metaframework, a Nexus for Exoplanet System Science (NExSS) science working group. The testing of the protocol was done using Met Office Software. The contributions of D. A. Christie, N. J. Mayne, D. Sergeev, and M. Zamyatina (Grant agreement No. 740963). The work of B. Charnay was granted access to the HPC resources of MesoPSL financed by the Region Ile de France and the project Equip@Meso (reference ANR-10-EQPX-29-01) of the programme Investissements d'Avenir supervised by the Agence Nationale pour la Recherche.