| 04 Aug 2025
biogeodyn-MITgcmIS (v1): a biogeodynamical tool for exploratory climate modelling
Abstract. Modelling the climate system is challenging when slow-response components, such as the deep ocean, vegetation and ice sheets, must be evolved alongside fast-response ones. This is crucial for investigating, for example, climate tipping elements and their interactions on the global spatial scale over multimillennia timescales. While Earth system models, such as those used in the Coupled Model Intercomparison Project (CMIP), are too computational expensive for simulations spanning thousands of years, simplified parameterizations and coarse resolutions in Earth Models of Intermediate Complexity (EMICs) can significantly affect the nonlinear interactions among climate components. Here, we describe a new tool, biogeodyn-MITgcmIS, which has a complexity level intermediate between EMICs and CMIP-class models. The core of biogeodyn-MITgcmIS is a coupled MITgcm setup that includes atmosphere, ocean, thermodynamic sea ice, and land modules. To this, we have added offline couplings with a vegetation model (BIOME4), a hydrological model (pysheds), and a new global-scale ice sheet model (MITgcmIS). The latter is implemented on the same cubed-sphere grid as MITgcm, using the shallow-ice approximation, as well as MITgcm outputs and a modified Positive Degree Day method to estimate the ice-sheet surface mass balance. Here, we describe in detail the new ice sheet model and the coupling procedure. We evaluate biogeodyn-MITgcmIS using simulations for the pre-industrial period and the 1979–2009 period. These two experiments allow us to assess the model's performance against CMIP-class models, as well as a combination of reanalyses and observations. biogeodyn-MITgcmIS successfully reproduces the large-scale climate and its major components, with results comparable to those of two CMIP models with dynamical vegetation. We discuss its potential applications and future developments.
Summary
Moinat et al. present a modelling tool that consists of the MITgcm augmented with additional offline components for vegetation, runoff and ice sheets. The model is run at intermediate resolution and mainly intended for deep time steady state simulations. Simulations under pre-industrial and present-day forcing are compared to other models and observations.
General comments
I have struggled quite some time to understand what this modelling tool is supposed to be used for. It is not until page 6 that this is clearly spelled out ("Since our purpose is to investigate climatic steady states in which ice sheets are in balance with the ocean, atmosphere, and biosphere"). If this is true, the statement should already be present in the abstract and the introduction should be geared towards that purpose. However, it is important to understand that in the real world, an ice sheet is never really in balance with the prevailing climate. While there is also mention of transient simulations, I believe in practice the model can only be used for time slice simulations, and only for those where assuming ice sheets in steady state with that climate is an acceptable assumption. The premise of the presented modelling (operating with climate and ice sheet steady states) excludes the application of the tool to interactions between climate and ice sheets that happen on timescales shorter than the relaxation time of the slowest component (the ice sheets), which is 40-100 kyr (l289). This may e.g. be illustrated by the Greenland ice sheet losing large amounts of its volume in run2 as it is relaxed to a steady state in the present climate, rather than attaining a transient state.
I can see two options for the manuscript to continue, both involving considerable work:
1. There may be value in this tool being applied to simulating deep time climate states (as hinted at e.g. in the title and model name) with steady state ice sheet configurations. In that case, the manuscript has to be strongly refocussed, aspects of transient simulations removed, comparison focussed on the preindustrial and a deep time use case added and discussed.
2. It may be possible to modify the model setup to allow for transient ice sheet simulations. It is common practice to run asynchronously coupled simulations between climate and ice sheet models where the climate forcing is accelerated. In that case the ice sheet component is the 'pacemaker' of the experiment and determines the physical time evolution. In that case more work would have to go into improving the modelling. Adding another example for a climate state different than the PI/PD would also be needed here.
The following comments are mostly independent of this choice.
What I miss in the introduction is a clarifying discussion of limiting factors in ESM modelling: spatial resolution, complexity and timescale. With limited HPC resources compromises are typically made on all three of them with different weight depending on the specific application of interest. What is the benefit of the specific choices made here? There is a relatively high resolution climate component compared to traditional EMICs but very low resolution of the ice sheet components. I don't understand the latter choice, since ice sheet models are typically the least expensive component and it would be feasible to run at a reasonable resolution without slowing down the system considerably. If the ice sheet resolution is not improved, it should at least be made clear that the current setup has a spatial resolution of the ice sheet models that is well below of what is considered appropriate to resolve ice sheet dynamics.
I think it should be demonstrated that the model can be set up for a considerably different climate. For a tool that supposedly can run multi-millennial simulations with ice sheet components it would be reasonable to at least be able to show an LGM configuration. Ideally I would also like to see a warmer climate configuration in addition like the last interglacial or the Pliocene. The setup of the tool may be in principle the same for any other time slice, but as modellers we know that there are always more complications once it is actually done. As is, the pre-industrial and present-day climate are too similar to serve as good validation experiments alone.
The analysis is heavily focussed on comparison to other CMIP models and observations, (where the latter is questionable due to the steady state nature of the modelling). Since I imagine that the uncoupled MITgcm has been validated extensively before, I believe there should be more focus in 3.2.1 and 3.2.2 on comparison to a pre-industrial state of the reference MITgcm without the additional components. This would show how the presented developments change the model behaviour.
Specific comments
Title and model name "a biogeodynamical tool"
I am a bit confused about the term "geodynamical". According to the Wikipedia definition (https://en.wikipedia.org/wiki/Geodynamics), this should involve mantle convection and plate tectonics. If a deep time application with different continent distribution and climate could be shown, this would make a bit more sense.
l2 "climate tipping elements"
The mention of climate tipping elements early in the abstract and in the introduction is confusing. Studying climate tipping points most people are concerned about requires full physical coupling between different climate components, which involves specific time scales, not steady states.
l7 "a complexity level intermediate between EMICs and CMIP-class models"
Consider another description of the complexity level of your model.
EMIC already contains the term intermediate. I would say your model is by definition an EMIC, albeit maybe of relatively high resolution of the climate component compared to many other EMICs (add references).
l9 "offline couplings"
Depending on where the modelling goes, I would suggest to use the term "asynchronous coupling" instead. "Offline" may be read as one-way coupling without feedback on the core climate model.
l10 "(MITgcmIS)"
The naming is a bit confusing. If MITgcm in itself is the climate model, MITgcmIS sounds like it could be the coupled climate-ice sheet model.
l10 "The latter is implemented on the same cubed-sphere grid"
This should be motivates at some point. It also has to be noted somewhere that this resolution is well below of what is considered appropriate to resolve ice sheet dynamics.
l12 "the new ice sheet model and the coupling procedure"
Focus is on the ice sheet component, which is not represented in the title. Consider modifying the title accordingly.
##l12 "We evaluate biogeodyn-MITgcmIS" and l14 "biogeodyn-MITgcmIS successfully reproduces the large-scale climate and its major components"
Offline coupling of the new components (l8) would suggest that the base climate state of biogeodyn-MITgcmIS is identical to the core MITgcm. I suppose MITgcm has been validated before. Clarify what is evaluated here specifically and why another response is expected.
l13 "pre-industrial period and the 1979-2009 period"
Difficult to imagine how to evaluate climate-ice sheet feedbacks for these rather steady periods. How does the model perform when the ice sheets are changing considerably?
l18 "CO2 concentration raises under the present-day climate crisis"
Strange combination in this sentence of future perspective and "present-day". Reformulate.
l32 "simulations cannot reach stationarity in slow deep-ocean and ice sheet dynamics"
Focus only on stationarity seems like an insufficient or even problematic goal when interest is in tipping, feedbacks and dynamics.
l35 "the nonlinear interaction among climatic components"
In your setup, this could only be addressed with high enough frequency of the asynchronous coupling and transient ice sheets.
l46 Add reference to Smith et al. 2021 (UKESM)
l56 "complexity that is in between EMICs and CMIP-like models"
See comment l7
l59 "We will provide a complete description of all the components"
A complete description of all the components of a coupled ESM does not fit in this paper. Maybe you meant a brief description of the climate components and a full description of the ice sheet component?
l74 "SPEEDY is in the intermediate complexity class"
See comment l7
l80 "planetary boundary layer"
Is this level used to provide boundary conditions for the ice sheet model? Usually PDD models ingest 2m air temperature. Discuss this as caveat of the coupling.
l90 "runoff map"
Is this runoff mapping, like a runoff routing scheme? It prescribed the pathway but not the amplitude? Clarify.
l90 "salinity and sea temperature for all the ocean levels"
Sounds like salinity and temperature are prescribed rather than dynamically calculated. Is that the case? Or is this just for the initial state? Clarify.
l90 "Orbital parameters can be set by specifying the obliquity, the duration of the day and the radiative influx from the sun"
Orbital parameters are typically obliquity, precession and eccentricity. I don't think "duration of the day" and "the radiative influx from the sun" are strictly speaking orbital parameters. Is it maybe possible to say that the effect of the orbital configuration can be set by specifying ...?
l91 "200 years per CPU day using 25 cores."
Does this mean you can run 200 years in 24h or 5000 years? Please present two numbers: 1) how fast is the model practically, i.e. how many years per day do you run with your typical configuration, 2) what is the efficiency of the model, e.g. how many CPU hours does a 100 year run consume.
l93 "as well as deep-time climates"
Since the focus of the section is not different climates but different palaeogeographies, I would reformulate here.
l107 BIOME4 is described here as a purely diagnostic component forced by output of MITgcm. Is there any aspect of the vegetation model that would feed back on climate? I am mainly thinking about albedo, but also runoff, since you mention water holding capacity.
l120 "28 including land ice"
What does that imply for the interaction between BIOME4 and the land ice component? You do not mention ice sheet extent as input to the vegetation model. What happens when the ice sheet retreats or expands? Can vegetation grow on land that has newly become ice free? What if vegetation is overrun by expanding ice sheets?
l130 "interested in horizontal resolutions of the order of 2◦ or coarser"
A bit strange to state that you are "interested" in such low resolutions. I think you should aim for a setup that resolves the processes you are interested to represent.
l130 "we can neglect basal melting and other fine-scale processes, as calving and ice streams"
I would rather say that you neglect those processes because they cannot be represented at the coarse resolution you are modelling. I still do not understand why you chose such a low resolution for the ice sheet components.
l132 "STREAMICE"
Could you explain why you are not using this module?
l135 "improved Positive Degree Day"
Is your method improved compared to Braithwaite (1977) or compared to other more recent PDD implementations?
l135 "as we also want to apply our coupled framework to paleoclimates."
Not sure to understand the reasoning. Is this an argument for using PDD instead of prescribing a given SMB?
Also here, I would like to see how the PDD results would look like for a considerably different climate.
l137 "We choose an improved PDD approach"
Same question as l135
l139 "This is an important limiting factor to consider that makes it difficult to apply SEB to paleoclimate simulations"
I am not convinced about the reasoning here. Similar limitations (too coarse resolution, missing important elements of the surface energy budget) certainly also apply to the PDD method, except that it is not as obvious/explicit as for an SEB method.
l163 "where τd = ρigHα, α= −∇zS"
I don't see the purpose of defining alpha here. It is also easily confused with "a" defined in l160
l173 "Since our purpose is to investigate climatic steady states in which ice sheets are in balance with the ocean, atmosphere, and biosphere"
This should be declared much earlier in the manuscript!
I am trying to make sense of the implications of this statement: an ice sheet fully adjusts to new climate conditions on a 100 kyr time scale. In the real world, an ice sheet is therefore never really in balance with the prevailing climate. I suppose this means that you are not interested in the coupled interaction between ice sheets and climate, which involves specific time scales?
l174 "we choose not to represent these higher-order processes"
In ice sheet terms, "higher-order processes" has a specific meaning. Please use another description.
l192 "ODE"
Introduce abbreviation.
l198 "The required input is the air surface temperature"
See comment l80
l199 "significant improvement in capturing early-season melting"
After reading l173 and in the light of all the other simplifications in the model this seems like a quite unimportant refinement to me. But I guess it is good to use recent developments. I would be concerned if the tuning of PDD model parameters presumably done at much higher resolution by Tsai and Ruan and e.g. with a geometry that resolves the ablation areas is translating to your setup. Have you considered retuning the PDD model?
Maybe more importantly for the deep-time application you envision, the assumption that the PDD with present day parameters translates to other orbital configurations is highly problematic
l209 "In summary,"
This summary concerns both 2.4.2 and 2.4.3 but appears in 2.4.3. I would suggest to describe ablation and accumulation in a joint section.
l212 "accumulation rates"
"accumulation" and "accumulation rate" are used interchangeably. Please be consistent in terminology.
l215 "Taking into account the isostatic correction"
Why "correction"? Maybe "isostatic adjustment"?
l216 "where a time delay is included"
Why not use instantaneous adjustment in line with your interest in steady states? This does not seem consistent with the approach in other components.
l232
Define Tnew and Tpickup
l233 "lapse rate is computed at each ice-sheet grid point using the MITgcm output"
Can you explain how this is done?
l235 "freshwater flux to or from the ocean is computed and included at restart at the ocean boundary of the ice sheet"
I am confused about how that works. Does the freshwater flux enter the ocean instantaneously? What determines the time scale of release?
l236 "To guarantee the conservation of salt, a compensation is performed at the global scale"
How does that work precisely? Is the model operating with physical freshwater fluxes or with (negative) salt fluxes?
l283 "the other with daily frequency (for MITgcmIS)"
It seems counterintuitive that the ice sheet component needs a higher frequency than the vegetation model. In some PDD implementations the seasonal cycle is parameterised to reduce the input to monthly or even annual. This should be mentioned and you should consider it to improve the performance of your model.
l285 "if the sea ice thickness"
Do you mean "ice sheet thickness"?
l289 "MITgcmIS is run for an equivalent of 40-100 thousand years"
Why "equivalent"? Isn't it run for an actual 40-100 kyr?
l290 "isostatic, lapse-rate and sea level corrections (Sec. 2.4.4) are included at this stage"
Are these corrections applied after the ice sheet model has run in MITgcm? Clarify. If so, how are feedbacks between isostasy and ice sheet SMB incorporated in this way?
L291 "pysheds is applied using these files"
What are "these files"?
l299 "the pickup files"
What is that?
l303 "describes in a consistent way the evolution of all the climatic variables"
I do not I agree with your notion of "consistent" evolution. The resulting evolution is not physically consistent in that the time scales of various processes are modified or removed. In the real world, e.g. the ice sheets do not have 40-100 kyr time to adjust to a new climate condition.
l304 Figure 1 "Transient simulations: repeat the same procedure by changing the forcing each N years"
See point l303. This setup is not made to run transient simulations with time steps shorter then 40-100 kyr.
l317 "with the addition of an isostatic correction"
Could you give a short summary of how that works? I am mostly interested if this is in line with the LLRA model you are using.
l321 "The resulting isostatic map"
What does "isostatic map" mean? Is this the bedrock adjusted for unloading of the ice sheets? If so, it would be interesting to see the difference to the reference bedrock in addition.
l324 "PANALESIS or other reconstructions directly provide bedrock elevations"
Not sure about the ice sheet loading in that case. Is that always included?
l325 "pre-industrial conditions at 280 ppm"
I am confused about combining present-day ice sheet reconstructions with pre-industrial climate. Is that the idea here?
l331 "set to a= 1.2·10−15 Pa−3 s−1"
What was the tuning procedure/target to find this value?
l333 Figure 2. I find this projection unusual and not very intuitive to interpret. Could you give additional figures with traditional lat-lon projection for global views and polar stereographic projections for the ice sheets?
l335 "The second simulation (run2)"
Can you give more details on how this run is set up? What is exactly the difference to run1 other than another CO2 concentration? Is this also a steady state run.
l336 "This run will be evaluated against the reanalysis and observational data"
Are you assuming that the real climate is in equilibrium with a 360 ppm CO2 forcing over this period? Are the ice sheets? How is it conceptionally possible to compare the steady state model agains transient observations?
l352 "5 hours of CPU time due to the daily stepping of the ODEs"
That seems very expensive for such a low resolution ice sheet model. Consider optimizing the SMB calculations.
l352 "for all land points"
Maybe you can exclude some points from the daily calculations when no ice is present?
l351 "one week to reach equilibrium"
Could you describe further how many years have been run for the different components?
l353 Table 2.
Specify meaning of variables in first column in the caption.
l356 Figure 4
Add latitude tick labels for NorESM2-LM
l361 "Excluding the ice sheet formation"
It would be interesting to see an analysis how the inclusion of the new components (vegetation, ice sheets) changes the model behaviour compared to the reference MITgcm.
l369 "which underestimates the elevation of Antarctica and Greenland"
Maybe this should be verified with a biogeodyn-MITgcmIS run where the ice sheet elevations are prescribed. Having that option in the model would anyhow be a good development.
l371 "the model capability to correctly reproduce the Hadley cells"
In line with comment l361, I seem to understand that MITgcm has been validated extensively. If true, this suggests strongly to focus the present analysis on how including the new components changes the model behaviour compared to the uncoupled reference in addition to comparison with other models.
l398 "the energy used for ice sheet growth is not included in this diagnostics"
Could it be? Shouldn't it be? Does that mean energy is not conserved? What is the result when you include it?
l412 "due to higher surface temperatures"
Could this be related to equilibrium climate vs transient climate. That's an important shortcoming of this comparison that should be discussed.
l438 "3.2.3 Vegetation"
I seem to understand that the offline vegetation model BIOME4 has been extensively validated. In addition you state "Offline coupling between the coupled MITgcm atmosphere-ocean-sea ice-land setup and BIOME4 has been already success-fully applied in Ragon et al. (2024, 2025)". I am wondering about the added value of section 3.2.3 in that context. Instead, it would be interesting to learn if including the ice sheet and runoff components changes the model behaviour. Is there an interaction between ice sheet change and vegetation? What vegetation grows in a much warmer climates where the ice sheets retreat?
l465
Why is there no evaluation of the runoff component here? That should be added.
l466 "3.2.4 Ice sheet"
I find it very difficult to justify the presented detailed comparison with observations of ice sheet models at such low resolution. This is particularly true for Greenland, which consist of only 25 grid points in the PI simulation.
l469 "evaluating the SMB produced for Antarctica is a prerequisite to calibrate the Glen’s law parameter"
The text until line 484 is not well placed in the Results section. I think this should be described earlier in the manuscript.
l471 Figure 10
Very unusual projection to display ice sheet results. Please consult some recent publications for inspiration. The standard are polar stereographic projections EPSG:3413 for Greenland and EPSG:3031 for Antarctica.
l486 "the Greenland ice sheet decreases strongly in height and extent"
This may be explained by the ice sheet being relaxed to a steady state in your model, while it is in a transient state in the real world. This just shows once more that the model is not meant to simulate transients with time scales of interest lower than ~50 kyr. In addition, I am concerned if the remaining 18 grid cells do a great job in simulating the dynamics of the ice sheet correctly.
l489 "after smoothing to the same spatial resolution"
How does smoothing change the total volume?
l514 "further improvements are possible"
There are quite a few that should be discussed here.
l527 "In future iterations of MITgcmIS, sliding and basal heat balance could easily be implemented"
In my view, the most important improvements are to increase the spatial resolution of the ice sheet models and make the SMB model more efficient.
l528 "this would allow study of nonlinear processes"
Before any of that, your model setup has to be modified to be able to run transient ice sheet simulations and not only steady states.
References:
Smith, R. S., Mathiot, P., Siahaan, A., Lee, V., Cornford, S. L., Gregory, J. M., Payne, A. J., Jenkins, A., Holland, P. R., Ridley, J. K., and Jones, C. G.: Coupling the U.K. Earth System Model to Dynamic Models of the Greenland and Antarctic Ice Sheets, Journal of Advances in Modeling Earth Systems, 13, e2021MS002520, https://doi.org/10.1029/2021MS002520, 2021.