CM2Mc-LPJmL v1.0: Biophysical coupling of a process-based dynamic vegetation model with managed land to a general circulation model

The terrestrial biosphere is exposed to land-use and climate change, which not only affects vegetation dynamics, but also changes land-atmosphere feedbacks. Specifically, changes in land-cover affect biophysical feedbacks of water and energy, therefore contributing to climate change. In this study, we couple the well established and comprehensively validated Dynamic Global Vegetation Model LPJmL5 to the coupled climate model CM2Mc, which is based on the atmosphere model AM2 and the ocean model MOM5 (CM2Mc-LPJmL5). In CM2Mc, we replace the simple land surface model LaD (where vegetation 5 is static and prescribed) with LPJmL5 and fully couple the water and energy cycles using the Geophysical Fluid Dynamics Laboratory (GFDL) Flexible Modeling System (FMS). Several improvements to LPJmL5 were implemented to allow a fully functional biophysical coupling. These include a sub-daily cycle for calculating energy and water fluxes, a conductance of the soil evaporation and plant interception, a canopy-layer humidity, and the surface energy balance in order to calculate the surface and canopy layer temperature within LPJmL5. Exchanging LaD by LPJmL5, and therefore switching from a 10 static and prescribed vegetation to a dynamic vegetation, allows us to model important biosphere processes, including fire, mortality, permafrost, hydrological cycling, and the impacts of managed land (crop growth and irrigation). Our results show that CM2Mc-LPJmL has similar temperature and precipitation biases as the original CM2Mc model with LaD. Performance of LPJmL5 in the coupled system compared to Earth observation data and to LPJmL offline simulation results is within acceptable error margins. The historic global mean temperature evolution of our model setup is within the range of CMIP5 models. The 15 comparison of model runs with and without land-use change shows a partially warmer and drier climate state across the global land surface. CM2Mc-LPJmL opens new opportunities to investigate important biophysical vegetation-climate feedbacks with a state-of-the-art and process-based dynamic vegetation model. 1 https://doi.org/10.5194/gmd-2020-436 Preprint. Discussion started: 18 February 2021 c © Author(s) 2021. CC BY 4.0 License.

use the code version 5.1.0 from the MOM5 project's git repository 1 . The model configuration is based on the accompanying test case named CM2M_coarse_BLING.

The Flexible Modeling System (FMS)
The Flexible Modeling System (FMS) is the coupler between the different model components of CM2Mc and has been developed at GFDL (Balaji, 2002). 2 FMS is a software framework for supporting the efficient development, construction, execution 90 and scientific interpretation of atmospheric, oceanic and coupled climate model systems. The infrastructure is prepared to handle the data interpolation between various model grids in a parallel computing infrastructure. It standardizes the interfaces between various model components and handles the fluxes between them. The flexibility of FMS allows for the relatively simple exchange of model components. All model components are simulated on different spatial and temporal scales and the coupler is the interface directly connected to the different parts. It interpolates the different scales to a common grid and adapts 95 the respective fluxes to the grid of the receiving model component. Usually the variables are not directly exchanged between model components. For instance, the land model calculates the humidity of the canopy layer, and the atmosphere the humidity of the lowest atmospheric layer. The coupler calculates the moisture flux between both layers and provides them to the different models on their respective spatial and temporal scale, while the different humidity variables are not exchanged. By tracking these explicit fluxes of energy and water, the coupler ensures the conversation of these quantities. CM2Mc employs GFDL's Modular Ocean Model (MOM) version 5 in a nominally 3x3 • lateral grid, with 28 vertical levels (Galbraith et al., 2011). Meridional grid resolution increases to a maximum of 0.6 • at the equator to allow the explicit simulation of some equatorial currents. The model uses re-scaled pressure vertical coordinates (p*), with the uppermost 8 layers having a thickness of 10 dbar, which increases with depth to a maximum layer thickness of 506 dbar (Galbraith et al., 2011). MOM5 105 utilises the tri-polar model grid of Murray (1996) to avoid a singularity at the North Pole and the use of partial bottom cells for a more accurate representation of bottom topography. Where the grid fails to resolve important exchanges of water between ocean basins, the cross-land mixing scheme of Griffies et al. (2005) is employed. MOM5 in CM2Mc is coupled to the GFDL thermodynamic-dynamic sea ice model (SIS, Delworth et al. 2006). We refer to Galbraith et al. (2011) for a more complete description of the model setup. 110 Enclosed in the ocean component MOM5, the Biogeochemistry with Light, Nutrients and Gases (BLING) model is run. It was developed at Princeton/GFDL as an intermediate-complexity tool to approximate marine biogeochemical cycling of key elements and their isotopes. More details can be found in Galbraith et al. (2011). 1 https://mom-ocean.github.io/ 2 https://www.gfdl.noaa.gov/fms/

Atmospheric Model 2
The atmospheric module in CM2Mc is GFDL's Atmospheric Model version 2.1 (AM2, Anderson et al., 2004). It uses the 115 finite volume dynamical core as in Lin (2004), as implemented in CM2.1 (Delworth et al., 2006) with a latitudinal resolution of 3°and a longitudinal resolution of 3.75°and 24 vertical levels, the lowest being at 30 m and the top at about 40 km above the surface. For the coupled setup, we use a general atmospheric time step of 1 hr at which variables are exchanged with the coupler. Dynamic motion and the thermodynamic state of the atmosphere are calculated on a 9 min time step, while the radiation scheme has a time step of 3 hrs. The coupled model includes an explicit representation of the diurnal cycle of solar 120 radiation. For a more detailed description of the model and its configuration, see Galbraith et al. (2011) and Delworth et al. (2006).

LPJmL5
The LPJmL5 (Lund-Potsdam-Jena managed Land) DGVM simulates the surface energy balance, water fluxes, carbon fluxes and stocks in natural and managed ecosystems globally and has been intensively evaluated (Von Bloh et al., 2018;Schaphoff 125 et al., 2018a, b). The model is driven by climate, atmospheric CO 2 concentration and soil texture data. Since its original implementation by Sitch et al. (2003), LPJmL has been improved by a better representation of the water balance (Gerten et al., 2004), the introduction of agriculture , and new modules for fire (Thonicke et al., 2010), permafrost (Schaphoff et al., 2013) and phenology (Forkel et al., 2014). In this study, we use the updated version of the fire model SPITFIRE as described in Drüke et al. (2019). All LPJmL (sub-)versions that build on the LPJmL5 version published by 130 Von Bloh et al. (2018), include the nitrogen and nutrient cycle. Because further adaptations would be necessary to include the nitrogen cycle in the coupled model, we concluded that it is beyond the scope of this study and deactivated it in this study.
LPJmL5 simulates global vegetation distribution as the fractional coverage (foliage projective cover or FPC) of plant functional types (PFTs, Appendix B) which changes depending on climate constraints and plant performance (establishment, growth, mortality). Plants establish according to their bioclimatic limits (adaptation to local climate) and survive depending on their 135 productivity and growth, their sensitivity to heat damage, light and water limitation as well as fire-related mortality. The interaction of these processes describes the simulated vegetation dynamics in natural vegetation. The model also simulates land use, i.e. the sowing, growth and harvest of 14 crop functional types and managed grassland (Rolinski et al., 2018). The proportion of potential natural vegetation and land-use within one grid cell is determined by the prescribed land-use input.
Each type of land cover, i.e. natural vegetation, managed grassland or crops, have their own respective stand. While receiving 140 the same climate information, soil and water properties as well as carbon-related processes are simulated separately.
In standard settings the model operates on a global grid with a spatial resolution of 0.5 • × 0.5 • . However, the actual resolution can be changed according to the spatial resolution of the model input.
To bring vegetation and soil carbon pools into equilibrium with climate, the model is run for a uncoupled spin up time of 5000 years, where the first 30 years of the given climate data set are repeated.

Adapting LPJmL5 to implement it into the FMS coupling framework
While Section 2.2 described the standard LPJmL5 model as previously published we introduce in Section 2.3 our adaptations to LPJmL5 in order to be coupled with the FMS coupling framework. An overview of our coupling approach between LPJmL5 and the CM2Mc model is provided in Fig. 1. The coupling software FMS, and hence the atmosphere model, expects a certain set of variables for full dynamic coupling. We consider canopy humidity, soil and canopy temperature, roughness length and albedo 150 as essential variables to allow dynamic vegetation to fully interact with the atmosphere, and describe their implementation in this Section. All these variables are exchanged with the atmosphere on the so-called "fast time step", for which we currently set one hour. Because the offline-version of LPJmL5 simulates carbon and water fluxes only at a daily time step, we introduced a sub-daily time step of the same duration as the fast time step and ensured a diurnal cycle for temperature and humidity which is important to stabilise the atmosphere and the coupled model system (Randall et al., 1991;Kim et al., 2019). These processes 155 included calculations of the water and energy cycles, i.e. surface temperature, evapotranspiration and water stress. Albedo and roughness lengths are expected to be less dynamic and are thus independent of the diurnal cycle. Hence, they are calculated in the original daily time step within LPJmL5, but still exchanged every hour. For ecosystems that are temporarily covered by snow, sublimation is implemented building on the simple snow model in LPJmL5, which also operates at the fast time step. In the fast time step, the coupling variables are sent from LPJmL5 to the FMS coupler. The coupler then provides the synoptic 160 climate variables (temperature, precipitation, radiation) as the input for LPJmL5 in the next (fast) time step.
In this Section we describe our coupling approach at the interface between the land model (LPJmL5) and the FMS coupler.
FMS calculates the fluxes between the different model components and provides these information to the sub-components.
The tasks of the coupler also include the calculation of air stability and surface drag, hence it has some functionality of a land surface model. Because it is beyond the scope of this paper to explain the processes within FMS in detail, we refer to Milly 165 and Shmakin (2002) and Anderson et al. (2004) for further details.

Interface between FMS and LPJmL5
The C main function of LPJmL5 used in the offline version is replaced by a coupler function providing the interface between the internal C functions of LPJmL5 and the Fortran functions of the CM2Mc model. The coupler function is called by FMS on an hourly time step and calls itself the specific update functions of LPJmL5 at the end of each hour, day, month or year, 170 respectively. Ingoing and outgoing data are transferred as array arguments of this function. The mapping of the coarse resolution of the CM2Mc model to the 0.5 • × 0.5 • resolution of LPJmL5 is done by the FMS coupler. We found that the FMS land model component must be run at LPJmL5 resolution, which is 0.5 • , so that all model components and the FMS coupler agree on which cells belong to land which to the ocean. This yields slight changes of the land-sea-mask from the original CM2M_coarse_BLING setup.

175
CM2Mc as well as LPJmL5 can use the Message Passing Interface (MPI) to run the simulation in parallel on a compute cluster.
CM2Mc uses FMS to set up a 2-dimensional domain decomposition, i.e. it splits the global grid into rectangular domains which are mapped to concurrent MPI tasks. In contrast, the LPJmL5 grid is represented by an unsorted 1-dimensional array of land cells, which is evenly distributed onto the MPI tasks. Since this LPJmL5 grid is not compatible with the FMS grid exchange framework, a small wrapper library for the data exchange between LPJmL5 and FMS domains was developed. The wrapper 180 library is called for the ingoing and outgoing data and the time overhead for this data exchange is negligible. The coupler function as well as the wrapper library are part of the LPJmL5 distribution.

New canopy module
The stand-alone version of LPJmL5 does not calculate the essential coupling variables canopy temperature and humidity, which is remedied in the coupled configuration via the addition of a new canopy module. In this new module, the canopy humidity 185 and canopy temperature and some further quantities linked to those variables are calculated (Fig. 2). In this setting, the canopy layer corresponds to the lower boundary for the temperature in the atmosphere. The atmospheric diurnal cycle as well as the seasonal changes depend on the surface energy balance. The canopy humidity, on the other hand, is the lower boundary for the atmospheric humidity and hence, sets the moisture content and the amount of precipitation in the atmosphere, as well as the potential for evapotranspiration on the surface. A schematic overview over the different calculation steps is provided in Fig. 3. 190 In the stand-alone version of LPJmL5 climatic input is prescribed, and therefore calculations of processes and fluxes, such as evapotranspiration, do not feed back to the atmosphere. In the coupled version, however, a small perturbation in a positive feedback loop can influence the climate and push the process towards an even larger perturbation. Therefore, special attention has to be given to ensure the stability of the model by either ignoring the feedback and implementing a simple, empirical and stabilizing relationship or by increasing the complexity of the implementation, in order to get a more realistic representation of 195 the vegetation embedded in the Earth system. The latter was done in CM2Mc-LPJmL by replacing the former simple Priestley-Taylor approach for calculating potential evapotranspiration ET 0 with the more complex and process-based Penman-Monteith evapotranspiration (Monteith, 1965). The Penman-Monteith approximation also accounts for additional parameters, such as humidity, that were previously not available in stand-alone LPJmL5 (Fig. 2):

200
where L v is the volumetric latent heat of vaporization of 2453 MJ m −3 , ET 0 is the evapotranspiration in m day −1 , dqsat dT the slope of the vapor pressure curve in kPa°C −1 , R n the net radiation at the surface in MJ m −2 day −1 , G the soil heat-flux density in MJ m −2 day −1 , 86400 the conversion factor from seconds to daily values, ρ a the air density in kg m −3 , C p the specific heat of dry air (1.013 · 10 −3 MJ kg −1°C−1 ), e 0 s the saturated water vapor pressure in kPa, e a the actual water vapor pressure in kPa, τ aν the bulk surface aerodynamic resistance for water vapor in s m −1 and τ s the canopy surface resistance in s m −1 . γ is the 205 psychrometric constant in kPa°C −1 and is calculated as: where P is the atmospheric pressure at the surface in kPa, λ the latent heat of vaporization of 2.45 MJ kg −1 and µ the ratio of molecular weight of water vapor to dry air, which is 0.622. ET 0 is presented here in the general daily form, but applied to the model on the subdaily timescale, therefore divided by the number of time steps per day (in the current version 24).

210
Eq. 1 uses the canopy surface resistance τ s , which is the reciprocal of the non-waterstressed canopy conductance g p in mm s −1 .
g p was also slightly changed, compared to Schaphoff et al. (2018a)   vapor pressure deficit (D).

215
where g 0 (mm s −1 ) is a PFT-specific minimum canopy conductance scaled by FPC, occurring due to other processes than photosynthesis. p a is the ambient partial pressure of CO 2 in Pa, A dt denotes the daily net daytime photosynthesis and 1000 is the unit conversion factor from mm to m. D (in Pa) can be obtained by the canopy humidity q ca and the saturation humidity q sat :

220
While the new potential evapotranspiration is calculated in the subdaily time step, the non-water-stressed canopy conductance is calculated in a daily time step, due to the daily calculation of the photosynthesis in LPJmL5. Since climate data from FMS is available on a subdaily basis, the photosynthesis routine uses a diurnal average of air temperature and photosynthetic active radiation. The newly calculated potential evapotranspiration, accounting for g p , is then also used in several LPJmL5 routines (e.g. bare soil evaporation or interception) instead of the equilibrium evapotranspiration (E q ), which was based on the Priestley- 225 Taylor formula (Schaphoff et al., 2018a).
As a next step, we calculate the water-stressed transpiration E tr , by using the supply-demand functions of LPJmL5 as follows: The demand is calculated by the newly implemented potential evapotranspiration (Eq. 1, corrected by the fraction used for interception) and the supply is driven by vertical root distribution and phenology (as in Schaphoff et al., 2018a). The initial transpiration is then a function of the minimum of supply and demand for water. The transpiration is then subtracted from the 230 various soil layers, depending on water availability. If the available water is not sufficient, transpiration decreases. The adjusted transpiration is consequently used in an inverse version of the Penman-Monteith formula in order to calculate the actual canopy conductance, linked to transpiration g tr .
The total canopy conductance is additionally influenced by the conductance of soil evaporation (g e ) and plant interception (g i ). Therefore, we use a simple approach taking into account the maximum rainfall interception conductance (GI MAX = 10 235 mm s −1 ) and by considering the fraction of rainfall i stored in the canopy of a biome-dependent rainfall regime (Gerten et al., 2004): where f v is the vegetated grid cell fraction and P r the daily precipitation. P r and ET 0 are here applied in mm s −1 . The soilevaporation conductance is calculated for the non-vegetated area of a grid cell and depends on the maximum soil conductance 240 (GE MAX = 10 mm s −1 , Huntingford and Monteith 1998), and an empirical scaling factor for the dependency of soil-evaporation conductance on soil-water status (α 0 = 10, Zhou et al. 2006): where w evap is the soil water content relative to the water holding capacity available for evaporation defined for a certain soil depth (Schaphoff et al., 2018a). Both conductances are calculated in the daily timestep. 245 We then calculate the total canopy conductance g c by adding g tr , g i , g e and using τ aν following Milly and Shmakin (2002).
where β ph is the water available for photosynthesis: with W r as the actual soil water and W * r as the maximum available soil water. The increment of the canopy humidity q ca per 250 time step is then calculated as following, using g c : where q flux is the water flux from the canopy layer to the atmosphere, provided by the FMS coupler, dT dt the gradient of the surface temperature over time and ET the final evapotranspiration, consisting of transpiration, evaporation, interception and sublimation from surface or vegetation into the canopy layer. For the calculation of ET we used the Penman-Monteith equation
The total canopy conductance and the final increment of the canopy humidity, which is important for the FMS coupler, are calculated in the subdaily time step. Eq. 9 is based on Milly and Shmakin (2002) and derived in Appendix C.
It was further necessary to implement the calculation of surface/canopy temperature within LPJmL5, therefore requesting major adaptions to the energy cycle in LPJmL5. Stand-alone LPJmL5 calculates the temperature of different soil layers by 260 employing a temperature transport scheme and taking into account air temperature as climatic input. In CM2Mc, however, the energy balance is calculated on the surface and then passed to the coupler and the atmosphere. Therefore, we had to implement this energy balance analogously in the coupled version of LPJmL5. While this surface temperature depends on several inputs from the coupler, as for instance radiation, it also uses several variables connected to the water cycle in LPJmL5 (evaporation, sublimation and melted water). Since our approach does not account for a height dependent canopy temperature, we used here 265 the surface temperature as an approximation for the canopy temperature, which is needed to calculate canopy humidity and evapotranspiration. Hence, surface temperature and canopy temperature are assumed the same, following the approach in the LaD model (Milly and Shmakin, 2002).
The soil temperature is still important for internal processes in LPJmL such as permafrost but not needed in the coupler to calculate fluxes from the land to the atmosphere. The calculation of heat transfer in the soil layers uses the heat-convection 270 scheme as in stand-alone LPJmL5 (Schaphoff et al., 2018a) by taking into account the air temperature, which highly depends on the canopy temperature. Both temperature calculations, for the surface/canopy temperature and for the soil temperature, operate on the fast time step.
In order to calculate the surface/canopy temperature within LPJmL5, we employed a simple energy-balance formulation for the incremental change of temperature ∆T for each time step (adapted from Milly and Shmakin, 2002): where m is the melted ice transformed to water in kg m −2 s −1 , LE f the latent heat of the conversion of ice into water in J kg −1 , LE v the latent heat of the conversion of water into vapor in J kg −1 , Q sn the released energy by snow in W m −2 , H the sensible heat provided by FMS in W m −2 , C s the heat capacity of the soil in J kg −1 and ∆ t the fast time step duration in seconds. R n is used here in the unit W m −2 . While the temperature is calculated individually for each stand, a weighted average over all 280 stands within one grid cell is used in the humidity calculation and passed to the coupler. The heat balance of snow is calculated as performed for the soil layers (see Schaphoff et al., 2018a) where snow temperature changes (∆T snow ) depend on the thermal conductivity (λ snow = 0.2 W m −2 K −1 ) and heat capacity (C snow = 630000 J m −3 K −1 ) of snow as follows: and heat flux from snow (Q snow ) is calculated: where z snow is the snow depth, T air is the air temperature and T soil [0] is the soil temperature of the first layer.

Albedo and roughness length
Albedo (β), the average reflectivity of the grid cell, is calculated as in Schaphoff et al. (2018a), based on a first implementation by Strengers et al. (2010) and later improved by considering several drivers of phenology by Forkel et al. (2014): where the albedo for bare soil β soil is defined as 0.3 and for snow β snow as 0.7. β PFT is calculated for each PFT depending on the foliage projective cover (FPC) and the stem, litter and leaf albedo of the respective PFT. The value for each parameter is as in Schaphoff et al. (2018a). F snow and F bare are the snow coverage and the fraction of bare soil, respectively. Water bodies as lakes and rivers have a constant albedo value of 0.1.

295
Roughness length z 0m is calculated according to Strengers et al. (2010): and where z b is the height of the boundary layer in stable conditions, set to 100m (Ronda et al., 2003), z i 0m is the PFT-specific 300 roughness length, and FPC i the foliage projective cover of each PFT, respectively. The coupler uses the roughness length to calculate aerodynamic resistance and surface drag and provides these variables to the different sub-models of the ESM.

Further changes in the coupled LPJmL5
For a global model we also need to consider Antarctica, which has not been part of the standard grid of the stand-alone LPJmL5 modelling configuration. It was implemented in a simplified approach, and will be replaced with the Parallel Ice Sheet Model 305 (PISM, Winkelmann et al. 2011) in the future. For now Antarctica is assigned the soil type ice and a constant albedo of 0.7.
The temperature balance is calculated as on the other continents.
In stand-alone LPJmL5, sublimation is subsumed by a constant global value of 0.1 mm per day, likely underestimating the sublimation at high latitudes. Especially in winter times, we do not expect much evapotranspiration, and hence the sublimation changes with meteorological conditions and becomes an important process. For this reason, we implemented the calculation of 310 sublimation E s by using the formula from Gelfan et al. (2004): where u is the wind speed in m s −1 from the coupler, e s the saturated vapor pressure in mbar and e a the air vapor pressure in mbar.
the effective rooting depths of the tropical-tree PFTs to 2.3 m in order to counter a negative AM2 precipitation bias in northern South America. Therefore, we increased the beta-value of each tropical tree PFT describing their vertical fine root distribution in the soil column from 0.96 as in Schaphoff et al. (2018a) to 0.99 in this study.

Model setup and forcing
In the stand-alone version, as well as in the coupled version, LPJmL5 is forced with gridded soil texture data (Nachtergaele code can employ hybrid MPI+OpenMP parallelism, but computational costs for LPJmL5 remain unaffected.

Modelling protocol
Soil carbon and vegetation biomass need timescales of hundreds to several thousand years to reach an equilibrium with climate, which would require extremely long spin-up simulations in the coupled model. Hence we produce a first spin-up for 5000 years with the more computational efficient stand alone LPJmL5, using climate input from GLDAS and an earlier CM2Mc-LaD run.

340
To bring vegetation, soil and climate into a consistent equilibrium (stand-alone LPJmL5 spin-up and the restart files from CM2Mc using LaD), we perform afterwards a fully coupled run of 500 simulation years under pre-industrial conditions with land use deactivated. The climate of this run is then used as forcing for another stand-alone LPJmL5 spin-up run of 5000 years, producing restart conditions much closer to the state of the coupled model. This multi-step spin-up approach minimizes the time for the computationally expensive coupled model to reach a stable state.

345
To account for changed dynamics in the coupled system, the LPJmL5 spin-up is then followed by a coupled spin-up, which runs for 500 years at pre-industrial and potential natural vegetation (PNV, i.e. without land use) conditions in a fully coupled  inter-comparison CMIP5 data (Taylor et al., 2012) and LPJmL5-offline (Section S4).
As evaluation metrics we used the normalized mean error (NME Kelley et al., 2013): where y i is the simulated and x i the observed value in grid cell i. x is the mean observed value. The NME is 1 if the model is as good as using the data mean as a predictor, larger than 1 for worse performance and zero for perfect agreement. We use this 390 metric for the evaluation of the performance of temperature, precipitation and above ground biomass.

Results
The evaluation of the model performance is provided in Section 3.1, while the impact of land-use change on the results of the coupled CM2Mc-LPJmL model is analyzed in 3.2.

Model performance 395
Here, we evaluate the performance of CM2Mc-LPJmL against climate and biosphere observations, by first looking into the long-term stability of global mean surface temperature (referred to as temperature, hereafter) and precipitation (Section 3.1.1) from the pi-Control experiment, before evaluating the historic temperature increase of the coupled model, using the TR exper-

Model stability
The analysis of the model stability was based on the pi-Control experiment, which ran over 750 years in total (see Section 2.5 for details). Here, we evaluate temperature and precipitation in terms of absolute values as well as rate of change over time and the variability.
After the initial 300 years, the global temperature remains relatively stable at ca. 14.7°C over the remaining simulation period 405 of 400 years with a slight drift of less than of 0.05°C per 100 years (Fig. 4a). The interannual variability in this period is ca. 0.1-0.2°C. The decreasing temperature over most of the 750-year simulation period can be explained by the energy uptake of the ocean, since deep ocean layers are not yet in equilibrium. The average precipitation follows a similar trend as temperature and reaches a relatively stable state at around 2.88 mm/day after ca. 400 years, changing less than 0.01 mm/day over the remaining period (Fig. 4b). The interannual variability is 0.01-0.02 mm/day.

Temperature evolution over the historical period
The temperature evolution over the historical period, hence the climate sensitivity to changes in atmospheric forcing, is evalu- In the PNV experiment, climate change is also well captured, but weaker as compared to having included land use in the model (Fig. S5).

Surface temperature evaluation 425
Basic climate patterns are well captured in the annual mean surface temperature (Fig. 6a), as temperatures are increasing from polar temperatures of below −10°C towards the equator with a maximum of ca. 25-30°C in the tropics. Desert regions are usually warmer, while mountainous regions are colder than the surrounding area. In the high latitudes ocean cells are usually a bit warmer than land cells, due to the ocean's ability to store heat.
Between 1994 and 2003 the average global temperature is 15.6°C compared to 14.3°C in the ERA5 data set with a NME 430 of 0.16. While the temperatures in the tropics and temperate zone are slightly overestimated (by ca. 1°C), the poles and the boreal zone show a large negative temperature bias (up to −10°C) (Fig. 6b). The Southern Ocean has a significant positive temperature bias (ca. 3°C on average). Large differences between CM2Mc-LPJmL and ERA5 are also visible for mountainous areas, where the temperature bias is partly due to the coarse resolution of the model, not adequately capturing the orographic influence of most mountain ranges on climate (e.g. Andes or Himalaya).

435
While the seasonal cycle is usually well captured in CM2Mc-LPJmL, especially in Antarctica a strong seasonal temperature bias is partly balanced out in the annual mean temperature. Temperature over Antarctica is largely overestimated during the southern-hemisphere summer, while being underestimated during the southern-hemisphere winter (Figs. S1 and S2).
The latitudinal distribution of modeled mean temperature between 1994 and 2003 (Fig. 6c) shows similar values compared to ERA5 data from high to mid-latitudes in the northern hemisphere, but a slight overestimation in parts of the temperate zone 440 and the tropics (between 70°S and 40°N). Specifically, the cold bias in the boreal zone leads to a slight underestimation of temperature between 60°N and 90°N.
The comparison of CM2Mc-LPJmL (pi-Control) and pi-CM2Mc-LaD (as in Galbraith et al., 2011) shows that similar biases in relation to ERA5 are present in both model versions. For example, both model versions slightly overestimate global temperature (Fig. S6). The strong regional biases as compared to ERA5 data are also present in both model setups (Fig. S6), hence not due 445 to the implementation of LPJmL5.

Precipitation evaluation
The spatio-temporal pattern of global precipitation is well simulated with a global average of 2.86 mm/day and a maximum of up to 10 mm/d in the tropics close to the Inter-Tropical Convergence Zone (ITCZ, Fig. 7a). Regions with little to no vegetation, such as deserts and polar areas, receive very little precipitation throughout the year.
Precipitation biases with respect to ERA5 data are, however, stronger than temperature biases with an NME of 0.50 compared to 0.16 for temperature (Fig. 7b). The biases are strongest at the equator with an apparent shift of the ITCZ. While precipitation in the Pacific is underestimated directly at the equator, it is overestimated north and south of the equator (Fig. 7b). Also northern South America shows a large negative precipitation bias.
The seasonal patterns (Figs. S3 and S4) confirm the imprecise modeling of the ITCZ, which remains for a large part of the 455 year north and south of the equator, while passing the equator region relatively swift. While precipitation south of the equator is overestimated, it is underestimated north of it.
The latitudinal annual mean precipitation between 1994 and 2003 (Fig. 7c)  Both models also show a large dry bias in northern South America (Fig. S6). Simulated AGB shows overall a good pattern, with largest values in the tropics, decreasing biomass in the subtropics and a local maximum in the temperate and boreal zone (Fig. 8d). In vegetation-free areas such as deserts or polar regions, simulated AGB 470 is zero or very close to zero (less than 200 gC/m 2 ). When comparing AGB against GlobBiomass (Fig. 8a), spatial differences emerge (Fig. 8c). While simulated AGB is slightly overestimated in boreal North America and Asia, it is underestimated in the European temperate zone and in Scandinavia, extending into eastern Europe and West-Siberia. In most of the other temperate, Mediterranean-type and subtropical regions, AGB matches the observed values. In the tropics, AGB is overestimated in semi-arid regions, whereas wet-tropical rainforests are mostly underestimated, especially the eastern Amazon. AGB shows 475 good agreement in the seasonal-dry Cerrado region in South America, but appears overestimated in the Caatinga in northeastern Brazil. In central Australia, AGB matches observations, but being overestimated in the north, and underestimated in the southeastern part of the continent (Fig. 8c).  The geographic distribution of dominant PFT cover in CM2Mc-LPJmL follows the spatial pattern of the biomass distribution 485 (Fig. 9a). The tropics are mostly dominated by the evergreen tree PFT. In the tropical savanna areas the tropical deciduous tree PFT dominates, along with the C 4 -grass PFT. The temperate zone is dominated by land-use with some summergreen trees most common in, e.g., Europe. The boreal zone is correctly covered by boreal needle-leaved and boreal summergreen trees and the tundra zone with polar grasses. To better visualize the model error for the PFT distribution, we produced an error map, which consists of the sum of the square error for each PFT per cell (Fig. 9b). In tropical rainforests, the error with respect to 490 the evaluation data is relatively small. Drier savanna areas show a much larger error, as well as parts of the temperate and the boreal zone. Areas with a small FPC fraction show a small error, because the error metric takes absolute errors into account.
This applies to desert regions in Africa, the Arabian peninsula and central Australia.
3.2 Impact of land-use changes on the coupled system 495 In order to isolate the impact of land-use change, we kept the climate constant and allowed land-use to change (LU-only, see Tab. 1). We compared precipitation, temperature and AGB for the years 2006-2015 of the LU-only experiment against the last 10 years of pi-Control to evaluate the absolute impact of changing land use. Most regions with a decreasing biomass and an increasing temperature show decreasing precipitation, e.g. the Brazilian Cerrado or southern Africa. This is due to reduced evapotranspiration of agriculture and pasture compared to natural vegetation 500 (Fig. 10a). Precipitation increases in regions where natural vegetation benefits from increased temperatures, for instance in mountainous regions, in India and in parts of southeast Asia (Fig. 10a).
Due to the replacement of natural vegetation by crops and managed grass, the total biomass is decreasing compared to the pi-Control experiment in regions with large land-use areas, e.g. Europe or the USA (Fig. 10c). As a consequence, surface temperature increases in these areas (Fig. 10b), leading to a global increase of ca. 0.5°C of average land-surface air tempera-505 ture. In the LU-only experiment, temperature additionally increases in regions where little to no land-use change occurred, e.g.
over northern Australia and Siberia (Fig. 10b). Over several sparsely vegetated areas, as in the Sahara, northeastern Canada and Greenland, temperature decreases. Temperature in tropical regions, e.g. in the Amazon basin and central Africa, are unaffected, as well as most desert and polar regions. For these regions, the amount of biomass remains the same as for the pi-Control experiment (Fig. 10 c).

Discussion
In this study we show the successful biophysical coupling of the whole-ecosystem DGVM LPJmL5 into the coarse-resolution version of GFDL's CM2 coupled climate model (CM2Mc), replacing the simple land-surface model of CM2Mc with LPJmL5.
In order to couple the stand-alone LPJmL5 to CM2Mc, some well-functioning model elements and structures had to be revised 515 and modified to work in a fully coupled climate model and to meet the essential coupling variables required by the coupler and the atmosphere modules. Even though LPJmL was developed as a stand-alone DGVM, its coupling to CM2Mc does not significantly change the temperature and precipitation patterns, but enables us to explore biophysical climate-vegetation feedbacks. The resulting model is furthermore in the range of CMIP5 models as stated in the Assessment Report 5 (Kattsov et al., 2013, Fig. S7).

520
In Section 4.1 we discuss the challenges of coupling LPJmL5 to CM2Mc and the evaluation of the coupled system, in Section 4.2 we examine the model application to simulate historic climate and land-use change, and in Section 4.3 we present an outlook on how the advantages of our modeling approach can be used best in future work.

Challenges of coupling LPJmL5 into CM2Mc
The results shown in Section 3 demonstrate that we achieved a stable model performance with respect to climate-biosphere 525 interactions after a potential natural vegetation spin-up period of 500 years. By achieving a stable climate in terms of surface temperature and precipitation, other variables in the model, as for instance carbon stocks of the biosphere (see Fig. S8 in the Supplement) and ocean carbon stocks, are also assumed to stabilize (even though possibly on a different time scale).
The climate variables temperature and precipitation show very similar biases as CM2Mc with LaD (see Figs. 6, 7 and S6). In other words, the relatively large bias in CM2Mc in certain regions occurs also when using the prescribed and idealized vegeta-530 tion cover from LaD, and is therefore not introduced by the coupling to LPJmL. The distribution of plant functional types and above-ground biomass are well simulated in most regions (Figs. 8 and 9).
The performance of the coupled LPJmL5 is directly sensitive to biases in the climate input produced by the AM2 atmosphere model. These biases can lead to a different vegetation state, which affects vegetation feedbacks to the atmosphere with possible increasing biases in AM2. This feedback loop is responsible for the deviations in our LPJmL vegetation results compared to 535 stand-alone simulation experiments without such feedbacks to the atmosphere. In the latter case, an error propagation from the climate input is avoided by forcing the model with bias-corrected climate data (Frieler et al., 2017). In our model approach we abstained from bias or flux corrections within the coupled model to maintain more realistic feedbacks, and allow its application to future as well as paleo-climate conditions. Furthermore, small problems in the parameterization of important processes can lead to larger problems in the whole state of the modeled Earth system. For instance, the temperature and water cycle calcu-540 lations have a strong interconnection and hence, a small error in the calculation of the water or energy cycle could lead to a runaway temperature and cause vegetation dieback for the wrong reasons. By adapting, e.g., the calculation of evapotranspiration and sublimation (see 2.3.2 and 2.3.4) we managed to keep the model relatively stable.
CM2Mc, when coupled either with LaD or LPJmL5, has a positive temperature bias of 1.3°C, which is within the range of published Earth system models (Kattsov et al., 2013). The temperature biases in CM2Mc are especially large in the polar and 545 in other at least partially snow-covered regions. In the northern latitudes a negative temperature bias led to a large mortality of vegetation in, e.g., Scandinavia in a previous model version (not shown). By adapting the simple snow model within LPJmL we obtained a stable vegetation of polar grasses and boreal trees in boreal Eurasia (see Section 2.3.4 for methods and Fig. 9 for results). A completely revised snow model or even a parallel ice sheet model could improve the modeling performance further.
Globally, the biomass cover is captured well by CM2Mc-LPJmL (Fig. 8). However, in an early development version of 550 CM2Mc-LPJmL a dry bias in northern South America led to a strong underestimation in the biomass productivity. The modeling was improved by using the above described Penman-Monteith parameterization for evapotranspiration (Section 2.3.2) and by increasing the tropical rooting depths and hence, the soil water access of the trees (Sakschewski et al., 2020). Global biomass patterns are now also comparable with the stand-alone LPJmL5 version (Fig. 8d).
Additionally, the coarse resolution of AM2 contributes to the simulated climate and vegetation anomalies, which can be usually 555 expected, when running fully coupled ESMs (Galbraith et al., 2011). While LPJmL runs in the native resolution of 0.5 • × 0.5 • , the atmosphere and hence the climatic input to LPJmL, has a resolution of 3 • × 3.75 • . While this resolution is necessary for a low computational cost, it can decrease the model accuracy over, e.g., mountain ranges such as over the Andes. The model smooths the height of the Andes to the coarse grid cell size, which leads to warmer temperatures on the high mountain areas and to a colder temperature on the low areas. Small biomes, such as the Caatinga in Brazil, have the size of a few grid cells or 560 are even smaller than one grid cell and hence, their unique climate can not be sufficiently captured by the coarse resolution of the atmosphere model. This could be improved by using a smaller grid size, but at the drawback of larger computational costs.
Since LPJmL accounts for large carbon stores, such as soil carbon, a long spin-up of several thousand years is necessary to get the carbon pools into equilibrium (Schaphoff et al., 2018a). To save computation time, this spin-up has been calculated with stand-alone LPJmL. Due to differences in the forcing of the stand-alone LPJmL version and the fully coupled model,

565
there is still a small offset in the beginning of the fully coupled spin-up run. After ca. 300 years, temperature and precipitation have reached a state close to an equilibrium (Fig. 4), and the model can be used for further scenarios and possible applications.
Without using the multi-step spin-up, as described in the methods (Section 2.5), the time to reach a stable state would be several times larger.

570
In addition to regional temperature patterns, the global temperature trends in historic climate and land-use change simulations are often used as another important evaluation metric, closely related to the climate sensitivity of Earth system models (Kattsov et al., 2013). Compared to GISTEMP evaluation data (Lenssen et al., 2019), the global temperature evolution over the historic period from 1860 until 2018 is well captured in CM2Mc-LPJmL (Fig. 5). The temperature increase in this period is also comparable to Kattsov et al. (2013). Therefore the model is able to model the response of the climate system and, hence, the 575 response of the biosphere to historic climate change.
To realistically model regional responses to climate change, the spatial temperature biases have to be taken into account. Temperature biases on land, which are sometimes up to 2 degrees Celsius, are larger than temperature increases during historic climate change. These biases have to be considered, when interpreting results from future model runs. Furthermore, the model does not account for climate modes and extreme events (e.g. El Niño Southern Oscillation), hence the interannual variabil-580 ity is smaller than expected. The interpretability of future runs is also hampered by the uncertain effect of CO 2 fertilization (Clark et al., 2013;Körner, 1993). This effect is relatively strong in LPJmL, leading to an increase in vegetation productivity at increasing CO 2 and temperature. The CO 2 fertilization effect under current climate has a stronger impact in LPJmL5 than heat stress in a warming climate. Activating the nitrogen cycle in LPJmL5, could reduce this strong effect by taking nitrogen limitation on vegetation productivity into account (Von Bloh et al., 2018). Historic biomass increase resulting from the CO 2 585 fertilization effect agrees, however, with previous studies (e.g. Zhu et al., 2016). A decrease in biomass in the historic period occurs almost exclusively in regions with land-use expansion.
Land use and land use management are often neglected in Earth system models, which leads to a inaccurate modeled temperature impact through land-use changes (Luyssaert et al., 2014). Since only ca. 30% of the land surface remains untouched by humans, a correct representation of land-use practises is important for modeling climate change of the 21st century (Levis, harvest and irrigation) for 12 different crop types.
By including land-use change in CM2Mc-LPJmL, natural vegetation is partially replaced by pasture and crops over time. This decreases biomass which affects the climate in three different aspects: 1) Less vegetation transpires less water, which decreases the water flux to the atmosphere, cooling by latent heat, humidity and precipitation (Gkatsopoulos, 2017), 2) the albedo of crops 595 is larger than that of closed forest, hence leading to a lower temperature (Unger, 2014), 3) the roughness lengths decreases, which increases temperature (Hoffmann and Jackson, 2000). While these effects mostly consist of a cooling through larger albedo and a warming through a smaller flux of latent and sensible heat, the net effect in CM2Mc-LPJmL is a warming climate in most areas. Especially in the tropics the latent and sensible heat fluxes outweigh a potential cooling by albedo increases.
The biophysical effect of land-use changes is furthermore highly sensitive to changes in roughness lengths and albedo for the 600 different PFTs and crop functional types, as well as different management options as, for instance, a different irrigation scheme (Kueppers et al., 2007).

Outlook
Using the advanced land use scheme of LPJmL5 and the capability of CM2Mc to accurately model climate change, the combined model CM2Mc-LPJmL is a powerful tool to model future trajectories of the Earth system. It allows to calculate various land-use change scenarios or management practises under changing climate in a computational efficient way. It is further 615 possible to separately investigate different biophysical processes and feedbacks, while forcing the model with representative concentration pathways (RCPs). Given the speed and relatively low computational cost of the model, even long term equilibrium experiment of several hundred years can be completed within days to a few weeks.
While CM2Mc-LPJmL is fully biophysically coupled, the biogeochemical coupling is not yet included. Each submodel accounts for a local carbon cycle and balance, but the carbon cycle is not yet closed for the whole model. For this study we 620 prescribed the atmospheric CO 2 concentration in all model runs and therefore a closed carbon cycle was not necessary. A fully closed carbon cycle is in the scope of future studies.
The key advantages of CM2Mc-LPJmL are the relatively fast and computational inexpensive atmosphere-ocean general circulation model (due to its relatively low spatial resolution) and the ability to investigate detailed feedbacks of the biosphere using the state-of-the-art DGVM LPJmL5. While LPJmL5 is constantly improved, recent new features such as a process-based nitrogen cycle (Von Bloh et al., 2018), a tillage system for land use (Lutz et al., 2019) or variable root growth (Sakschewski et al., 2020) can be integrated in the modelling framework consecutively and tested in the Earth system model. The coupled model also remains flexible for new model compartments such as a new atmosphere or a new ocean model, which are compatible with FMS. GFDL has already released the newest AM4 atmospheric model (Zhao et al., 2018), as well as MOM6 such as a state-of-the-art ocean model (Adcroft et al., 2019). Both could be integrated in the already existing modeling framework, and 630 are expected to further reduce model bias.

Conclusions
In this study we demonstrate the successful biophysical coupling of the state-of-the-art DGVM LPJmL5 into the coupled cli- The fully coupled energy and water cycle allows investigating the impact of biophysical atmosphere-biosphere feedbacks on global climate trajectories and quantifying impacts of deforestation or afforestation scenarios. CM2Mc-LPJmL might further help in identifying tipping points and planetary boundaries especially in the biosphere. By using LPJmL5 we can make, e.g., use 645 of its advanced land use scheme, the sophisticated process-based fire model SPITFIRE (Thonicke et al., 2010), a representation of permafrost and a state-of-the-art water cycling (Schaphoff et al., 2018a) and incorporate future model developments.

Appendix C: Derivation of humidity increment
Assuming equilibrium conditions the flux entering the canopy layer from soil and vegetation through evapotranspiration ET or E in equals the flux leaving the canopy layer into the atmosphere q f lux or E out .
The water fluxes for the next time step t+1 yield: using Using (Milly and Shmakin, 2002) and Eq. 7 from this paper yields for E: where ρ is the air density, r a the aerodynamic resistance, g c the canopy conductance, q sat the saturation humidity and q a the actual humidity. The derivation of Eq. C4 can be used for dEin dt . Eq. C2 then yields: Expanding dEout dt with q a yields: Rearranging Eq. C7 yields: Expanding dqs dt with dT for the temperature change yields: which is the final form for the change of actual humidity over a timestep. By using ET for E in , q f lux for E out and de dq for dEout dqa the final form yields: