Max Planck Institute for Meteorology

Abstract. Computational models play an increasingly vital role in scientific research by enabling the numerical simulation of complex processes. Such models are also fundamental in geosciences. For instance, they offer critical insights into the impacts of global change on the Earth system today and in the future. Beyond their value as research tools, models are also software products and should therefore adhere to certain established software engineering standards. However, scientists are rarely trained as software developers, which can lead to potential deficiencies in software quality like unreadable, inefficient, or erroneous code. The complexity of models, coupled with their integration into broader workflows, also often makes it challenging to reproduce results, evaluate processes, and build upon them. In this paper, we review the state and current practices of the development processes of the state-of-the-art land surface models used by the Global Carbon Budget. We combine the experience of modelers from the respective research groups with the expertise of software engineers from tech companies to outline key principles and tools for improving software quality in research. We explore four main areas: (1) model testing and validation, (2) scientific, technical, and user documentation, (3) version control, continuous integration, and code review, and (4) the portability and reproducibility of workflows. Our review reveals that while modeling communities are incorporating many best practices, significant room for improvement remains in areas such as automated testing, automated documentation, and reproducibility. Therefore, we here identify and promote essential software engineering practices, including numerous examples of practices from within the community that can serve as guidelines for other models and could help streamline processes across the entire community. We conclude with an open-source example implementation of these principles, demonstrating portable and reproducible data flows, a continuous integration setup, and web-based visualizations. This example may serve as a practical resource for model developers, users, and all scientists engaged in scientific programming.

Abstract. Although shrubs employ distinct water- and carbon-use strategies compared to trees and are increasingly expanding across warming tundra and grassland, they remain insufficiently represented in global land surface models. Here, we incorporated two shrub types, deciduous and evergreen, into the nutrient-enabled terrestrial biosphere model QUINCY, which features a state-of-the-art treatment of soil nutrient dynamics and carbon exchange. We investigate the change in carbon fluxes and storage due to shrub cover, it's response to climate and CO2 fertilization effect and the role of nitrogen availability. With this new implementation, shrubs showed reasonable seasonal cycle of gross primary production (GPP) at 50 % of the Arctic study sites. The model achieved mean R2 values of 0.5 and 0.6, when compared with in situ measurements and remote sensing products for modeled shrubs. However, at 50 % of the study sites the model underestimated observed GPP due to too strong simulated nitrogen limitation. Compared to needle leaved evergreen forest the modeled gross primary production of shrubs is similarly distributed with a non-significant difference in the median. Compared to graminoids the carbon fluxes of shrubs are 40 % higher. Shrubs produce a substantial, though lower, above-ground biomass than needle leaved trees and show phenological patterns that are distinct from those of trees. Although CO2 fertilization generally benefits all plant types, shrubs appear to maintain a particularly strong growth response under elevated CO2 concentrations. We also demonstrated that the modeled deciduous shrubs reduce their nitrogen sources substantially more than evergreen shrubs, generally resulting in a 50 % decrease in gross primary production. Providing the plants with unlimited nitrogen and thus doubling gross primary production at most sites improved the model-measurement agreement by 15 %. A similar effect occurred when initializing nitrogen and carbon contents best on permafrost profiles, resulting in partly alleviating nitrogen limitation in the model. These finding underlines the importance of including evergreen and deciduous shrub PFTs in global land surface models to accurately predict ongoing changes in the Arctic carbon cycle. However, the strong nitrogen limitation of Arctic shrub productivity when using the standard model parametrizations suggests that the Arctic contribution to global land carbon is underestimated by global models.

Accurate estimates of land-use change CO2 fluxes (FLUC) are essential for understanding the terrestrial carbon cycle and for informing national and global climate reporting. FLUC can be quantified using a range of approaches, including bookkeeping models, dynamic global vegetation models (DGVMs), national greenhouse gas inventories, Earth observation-based products, and atmospheric inversions. However, systematic comparison of FLUC estimates across approaches remains challenging due to fundamental methodological differences. Additionally, individual approaches are often associated with substantial uncertainties.Here, we provide an overview of the main sources of uncertainty and methodological discrepancies across approaches, while highlighting recent advances and identifying promising directions to further harmonize and improve FLUC estimates. Key sources of uncertainty and inconsistency include differing objectives and definitions across approaches, differences in the separation of natural and anthropogenic drivers (leading to FLUC differences of up to 2.8 PgC yr-1), incomplete process representation (contributing up to 30% uncertainty in FLUC), uncertainties in land-use data (up to 30% uncertainty), and constraints related to spatial resolution. At the same time, several important methodological improvements have been achieved recently, including the consideration of environmental changes in bookkeeping-based FLUC estimates and the correction for replaced sinks and sources (RSS) in DGVM-based FLUC estimates. Ongoing research projects are addressing several remaining challenges, including the improvement of regrowth estimates, the quantification of disturbance impacts, and the consolidation and harmonization of different land-use datasets. These developments build primarily on combining high-resolution Earth observation-based databases with the flexibility and traceability of semi-empirical modelling. Together, these advances provide an important contribution towards more consistent and robust estimates of land-use change CO2 fluxes within and across approaches.

The land–ocean warming contrast—where global land mean warms more than the global ocean mean —is one of the most prominent features of global warming. In idealized CO₂-increase experiments, although there is inter-model spread, it is known that land warming is typically about 1.6 times larger than ocean warming. However, it remains unclear what determines the magnitude of this land–ocean warming contrast. For the 1980–2014 trend, observed land warming exceeds ocean warming by more than a factor of 2.3, and such a large observed land–ocean warming contrast cannot be reproduced by CMIP6 historical simulations. It is also well known that during this period there is a mismatch in the SST trend pattern between observations and historical simulations. In this study, we investigated how differences in sea surface temperature (SST) patterns affect the magnitude of land–ocean warming contrast. Using the climate model MPI-ESM, we conducted AGCM experiments forced by (i) a globally uniform SST warming (+2 K, +4 K, and +6 K), (ii) the same global-mean SST warming superimposed with the observed 1980–2014 SST trend pattern, and (iii) the same global-mean SST warming with the sign of the 1980–2014 SST trend pattern reversed. The results show that, despite having the same global-mean SST warming, land warming differs significantly among the SST patterns, and that the observed SST trend pattern tends to enhance the global-mean land warming. Under the observed SST pattern experiments, warming tends to be amplified across the entire Eurasian continent, which is a major contributor to the enhanced global-mean land warming. The strong warming over the mid-to-high-latitude Eurasian continent is explained primarily by pronounced Atlantic warming and warming in the northwestern Pacific, whereas the cooling tendency in the eastern equatorial Pacific affects the land-warming pattern over the North America through teleconnections. This study demonstrates that SST patterns exert a substantial influence on the factors controlling the magnitude of the land–ocean warming contrast, and suggests that the coupling between ocean and land temperature changes varies markedly depending on the future SST pattern change.

Some major volcanic eruptions, such as the one of Mt. Pinatubo in 1991, can inject large amounts of sulfur dioxide (SO2) into the stratosphere, leading to the formation of a volcanic aerosol cloud. This dense aerosol cloud induces radiative heating of the stratosphere, causing ozone and water vapour changes, thereby altering middle atmospheric dynamics and chemistry. The scale of these impacts on stratospheric temperature anomalies is still highly uncertain.In this study we analyse data from the Historical Eruptions SO2 Emission Assessment Protocol (HErSEA) under the Interactive Stratospheric Aerosol Model Intercomparison Project (ISA-MIP). The results from eight global interactive-aerosol models confirm our general understanding of the stratospheric aerosol forcing due to SO2 injection following a volcanic eruption. As direct observations are sparse we compare models to three widely used reanalyses (ERA5, MERRA2, and JRA55). This analysis shows that while the multi-model mean temperature anomalies agree well with reanalyses, differences among individual models can be large. Our study shows that agreement in the median occurs through error compensation when averaging across models. The analysis and the sensitivity tests for model selection presented here highlight that by far the most important factor driving both magnitude and spread of the multi-model distribution in temperature response to volcanic aerosol forcing is model choice. Differences in transport, radiative transfer, and microphysics as well as the characterization of aerosol size distributions play a crucial role for the simulated spread in the temperature response.Another candidate to explain the spread in the ISA-MIP models, is the use of interactive aerosol schemes. To test this hypothesis, we compared the ISA-MIP multi-model distribution with those obtained from CCMI-2022 and CMIP6-AMIP model intercomparisons, which use prescribed SADs. If indeed interactive aerosol treatment would be a key contributor, one would expect smaller multi-model temperature anomaly distributions from CCMI-2022 and CMIP6-AMIP. Interestingly, this hypothesis has to be rejected, as no reduction in the multi-model spread is found. Hence, we argue for caution in attribution studies and the interpretation of stratospheric aerosol injection experiments relying on individual or few models.

The Arctic permafrost is warming widely, and the circumpolar region is heating up four times faster than the global average. Common features of permafrost landscapes are drained lake basins (DLBs) which play an important role for the geomorphological, hydrological and the ecological development of those landscapes. In addition, rivers transport significant amounts of freshwater, dissolved organic carbon, and other materials into the Arctic Ocean. The impact of close-by streams for DLB drainage and refilling events needs to be quantified.Here, we focus on DLBs and their connections to river floodplains on the Yamal Peninsula in northern Siberia, a region underlain by both discontinuous and continuous permafrost and covered with tundra vegetation, thaw lakes, and wetlands. For this study, we investigated manually selected DLBs using DEM derivatives. The lakes represent a North-South climatic gradient and different drained lake basin development stages. DEM derivatives are also used to represent related small-scale landscape features in climate models. The utility of metrics based on the Copernicus DEM, which were recently implemented in the ICON-Land model, was assessed. Potential floodplain linkage could be identified in many cases, leading also to diverting patterns of landcover evolution after drainage.The results provide a foundation for further analysis of lakes and their connections to streams, which play a role in the wetting and drying processes and subsequent impacts on the carbon cycle across the Arctic.

One of the main sources of uncertainty in climate modeling is microphysics. While Lagrangian approaches are promising, their computational cost is not yet suited to global climate simulations, leaving global simulations to bulk schemes. Starting from simple one-moment microphysics formulations that prognose only the specific mass of each hydrometeor 𝑞x, physically richer approaches like double- and higher-moment bulk schemes were introduced to narrow model–observation gaps, but they deliver mixed, metric-dependent gains over simpler one-moment schemes [1,2]. This mixed record suggests that the link between microphysical processes and other atmospheric processes remains poorly understood. Few studies have probed this link directly [3], and the intricacy of bulk schemes has hindered clear attribution of changes in climate statistics to specific microphysical processes. Here we aim to make microphysics more transparent by simplifying it to a Kessler-like three-category scheme. We thus collapse Lin’s six-category scheme [4] to three prognosed variables: vapor 𝑞v, condensates 𝑞ci (cloud water 𝑞c + ice 𝑞i ), and precipitates 𝑞rgs (rain 𝑞r + graupel 𝑞g + snow 𝑞s), where their precise category is defined by local thermodynamics (temperature 𝑇 and relative humidity RH) and updraft velocity 𝑤 (Fig. 1). Using a global 5-km atmosphere-only ICON simulation with a Lin-like microphysics scheme, we ask whether accurate mappings 𝑓x , 𝑓y exist such that Condensates: 𝑞x ≈ 𝑞x* = 𝑓x (𝑞ci , 𝑇, RH, 𝑤), for x ∈ {c, i} , Precipitates: 𝑞y ≈ 𝑞y* = 𝑓y (𝑞rgs , 𝑇, RH, 𝑤), for y ∈ {r, g, s} ,where 𝑞x,y are the specific masses from our ICON simulation, and 𝑞*x,y are the predicted specific masses from the partitioning functions 𝑓x,y , given 𝑞ci and 𝑞rgs (see Fig. 1). We construct histograms in the (𝑇, 𝑤, RH)-space and fit simple partition functions, like sigmoids, to build our partitioning functions. We present here results from this mapping, as well as an evaluation of its performance. We measure the R², and we compare global instantaneous outputs from our ICON simulation to the predictions provided by our partitioning functions, such as LWP and IWP comparisons. Fig.1: Collapse state of microphysics and its link to the predicted partitioning. [1] Seiki, Kodama, Noda, Satoh, J. Climate, 28, 2405–2419 (2015).[2] Song, Sunny Lim, Weather and Climate Extremes, 37, 2212-0947 (2022).[3] Proske, Ferrachat, Neubauer, Staab, Lohmann, Atmos. Chem. Phys., 22, 4737–4762 (2022).[4] Lin, Farley, Orville, J. Appl. Meteorol. Climatol., 22, 1065–1092 (1983).

In recent years, German forests have faced numerous stress episodes, such as droughts, heatwaves and insect outbreaks, leading to a significant rise in tree mortality. Consequently, deadwood production increased, partly outpacing extraction capacities. This trend coincides with a shift towards more selective forest management methods that partially retain deadwood on logged sites. As a result, the average amount of deadwood per hectare in Germany has substantially increased, climbing from 19.9 m3 ha-1 yr-1 in 2012 to 29.4 m3 ha-1 yr-1 in 2022.Besides providing valuable species habitats and sequestering carbon, deadwood plays a significant role in the water cycle, by reducing surface runoff and acting as a water reservoir.Despite its growing abundance and ecological relevance, deadwood is largely neglected in current land surface models such as in JSBACH (ICON-LAND main component), particularly with respect to its impact on hydrology and microclimates.This omission is rooted in the historical removal of deadwood under conventional management and in limited understanding of deadwood's hydrological properties, which vary with factors such as soil characteristics, canopy closure, deadwood position and species.Here, we present an approach to quantify and model the influence of deadwood on hydrology and microclimates, by combining experimental field measurements from Bavarian Nature Reserves and the data from the LabForest project in the University forest of the LMU in southern Germany with JSBACH outputs. We compare different methodological pathways for deriving observation-based parameterizations suitable to model integration. Preliminary results indicate that deadwood exerts measurable effects on near-surface microclimate and soil moisture dynamics, highlighting the need to explicitly represent deadwood in land surface modeling frameworks.

The land subsurface stored around 6% of the Earth’s energy imbalance in the last five decades (of around 0.5 Wm-2, equivalent to 380 ZJ), being the second contributor to the energy partitioning after the ocean (90%). Previous studies have shown that state-of-the-art Earth System Models (ESMs) remarkably underestimate the observational land heat uptake values. This underestimation stems from Land Surface Models (LSMs) within ESMs imposing too shallow zero-flux bottom boundary conditions to correctly represent the conductive propagation and land heat uptake with depth. When realistically deep boundary conditions are prescribed, land heat uptake increases by a factor of five. However, changes in ground surface temperature are negligible. The reasons for this lack of impact of the LSM depth on surface temperatures are assessed herein.An ensemble of eight historical and RCP8.5 land-only simulations with different subsurface depths was conducted with the LSM of the Max Planck Institute for Meteorology ESM (MPI-ESM), JSBACH. Simulation-derived latent (LHF), sensible (SHF), and ground heat fluxes (GHF) were compared across simulations, and GHF was additionally evaluated against estimates from a one-dimensional heat conduction forward model. Results show that, for a global warming of 1.5 ºC with respect to 1850-1900, GHF increases from 0.04 to 0.07 Wm-2 when deepening the LSM from 10 to 22 m, saturating at around 0.12 Wm-2 when the boundary condition is placed at approximately 100 m. The increase in the incoming GHF is mainly compensated by a global decrease in the outgoing SHF, a small decrease of the LHF in wet regions, and a decrease in the surface net radiation in arid and semi-arid regions. These quantities, yet small, evidence that an insufficient LSM depth induces to an inaccurate resolution of the long-term surface energy balance, which may have implications for land-atmosphere interaction. Their accumulation over time also produces biases in the terrestrial energy partitioning.

Arctic permafrost has been identified as a critical element in the global climate system, since it stores a vast amount of carbon that is at high risk of being released under climate change. The feedbacks between permafrost carbon and climate change are moderated by complex interactions between physical, hydrological, biogeochemical, and ecological processes. Many of these are not well understood to date, and therefore also only rudimentarily represented in current Earth System Models (ESMs).The Q-ARCTIC project funded by the European Research Council (ERC) follows a synergetic approach by combining remote sensing and local-scale observations with modeling on scales from a few meters to hundreds of kilometers. The primary objective of Q-ARCTIC is to close the gap between process scales and the much coarser grid resolution of Earth System Models (ESMs), with a particular focus on the net effect of disturbance processes and associated changes in hydrology on the pan-Arctic scale. To close this gap, we developed new ESM modules representing subgrid through stochastic parameterizations, trained and evaluated with high-resolution remote sensing data and site-level observations.We will present novel results from in-situ experiments that quantify carbon fluxes and environmental response functions at patch-level within heterogeneous Arctic ecosystems, resulting in optimized strategies to integrate data streams for upscaling. Satellite remote sensing products investigate fine scale (few meters) patterns in Arctic landscapes that are undergoing modifications linked to climate change, including e.g. InSAR data to constrain ground ice content, and related subsidence patterns. Targets investigated include for example sinking surfaces, wetness gradients in heterogeneous landscapes, and drained lake basins. Assimilation of these new datasets supported the development of new ESM model components that successfully capture the statistics of small-scale features, including e.g. lateral connections between hydrologic landscape elements across scales, or thermokarst lake dynamics. Our results demonstrate that the ability to project the response of the high-latitude water, energy and carbon cycles to rising global temperatures may strongly depend on the ability to adequately represent the soil hydrology and disturbance effects in permafrost affected regions.