1Boundary conditions¶

1.1DONE: Soil properties¶

ORCHIDEE reads soil property maps prescribing soil texture, pH, and bulk density. The soil texture is used primarily to determine the hydraulic and thermal properties of the soil, and the soil pH is used in the soil nitrogen calculations. Bulk density is currently not used, but is required by upcoming developments.

The default soil texture map is based on the 5-arc-min (1/12°) map by Reynolds et al. (2000), while the simplified three-class scheme uses the 1° map of Zobler (1986). Both of these maps in turn are based on the 1:5,000,000 FAO/UNESCO Soil Map of the World. In addition to these, any soil texture map using the USDA classification, such as SoilGrids Hengl et al., 2014, can be used directly by setting a runtime parameter.

The soil texture map is regridded onto the model mesh by assigning to each grid cell the soil texture class that covers the largest fraction of the grid cell. In addition to this, granulometric sand, silt, and clay fractions are calculated for each grid cell as the area-weighted means of the values typical for the textural classes present in the grid cell. In both cases, soil horizons are not distinguished and all tiles present in a grid cell share the same textural class and granulometric composition.

Soil pH and bulk density are both read from 1° maps. The pH map was generated using the SoilData program from the IGBP-DIS CD-ROM International Geosphere-Biosphere Program Data and Information System, 2000. The bulk density map is derived from the reference topsoil bulk density of the Harmonized World Soil Database version 1.1 Nachtergaele et al., 2009 or 1.2 Nachtergaele et al., 2012 [which?]; The soil PH and BulkDens description should be checked !.

1.2DONE: Vegetation distribution¶

Various historical reconstructions of the vegetation distribution ($f^{veg,max}$) can be used with ORCHIDEE. We describe below the main one that is used for global applications. It combines information from the land use harmonisation database (LUH v3.1.1; Update ref LUH3 ?) Hurtt et al., 2011 with land cover information derived from satellite observations, the Medium Resolution (300 m) Land Cover product (MRLC 2.8) from the Climate Change Initiative (CCI) of the European Space Agency (ESA) (referred as CCI-MRLC, Bontemps et al. (2015)). Note that in the latest CCI-MRLC product that we use (described in Harper et al. (2023)), the original 38 land cover classes Di Gregorio, 2005 are directly re-mapped onto a set of 16 generic PFTs at 300 m resolution (see https://orchidas.lsce.ipsl.fr/dev/lccci/generic_pfts.php), covering the period 1992 to 2020. In addition, in this product the grassland PFTs are divided into separate PFTs for grassland with the C3 or C4 photosynthetic pathway. Combining LUH v3.1.1 with satellite-derived generic PFTs consists of several processing steps that are described briefly below and in details in Olivera et al., (in preparation):

1. The bioclimatic zones (Köppen–Geiger climate classification map) were used to split wide-spread generic PFTs into separate generic PFTs for the tropical, temperate and boreal zones.

2. The generic PFTs are mapped into the ORCHIDEE specific PFTs. In the default PFTs classification (Table 6) shrubs were classified as 60% tree PFT and 40% grass PFT (80% tree and 20% grass for boreal zone) and lakes were classified as bare soil, unless the Flake model is activated (see section 1.9).

3. LUH v3.1.1 contains a land cover reconstruction from 850 to 2024 and several land cover scenarios from present-day to 2100 Hurtt et al., 2011 (update ref). These maps had to be merged with the ORCHIDEE PFTs defined in the previous steps. Adjustments were as follow for each year of the CCI-MRLC period (1992 to 2020): a) for each grid cell the crop fraction was taken from LUH v3.1.1 using the proposed C3/C4 split of that product. The pasture fraction from LUH v3.1.1 was pre-assigned to ORCHIDEE grass PFTs (with further split into C3/C4 pathway) b) the remaining fraction of the grid cell ($f^{remain}$) corresponds to natural vegetation. Present day fractions of natural PFTs (ORCHIDEE PFTs from step 2) were rescaled to match $f^{remain}$. For the historical reconstruction (years 850-1991) we use the average over 1992-1996 (first five years of the satellite derived PFTs) of the present day fraction of natural PFTs. For the recent present-day maps (years 2021-2024) and future projection (years 2025-2100) we used the last available five years of satellite derived data (2016-2020). More details and illustrations are provided in the website https://orchidas.lsce.ipsl.fr/dev/lccci/tools.php.

For model simulations that only require vegetation reconstruction over the satellite-era period (i.e., 1992 - onwards) we directly use the CCI-MRLC product, without step-3 described above (except for the crops fractioning into C3/C4 pathway that is always derived from LUH v3.1.1). Note that the ORCHIDEE team follows the regular update of the CCI-MRLC product; the latest one being described in Harper et al. (2023).

1.3DONE: CO₂ forcing¶

ORCHIDEE uses the atmospheric CO₂ concentration (ppm) to compute the stomatal conductance and CO₂ assimilation (see section 1.2). By default, the annual time-series from the TRENDY intercomparison protocol is used Sitch et al. (2024). It is derived from ice core CO₂ data (starting in 1700) merged with NOAA observations from 1958 onwards. The construction of the dataset, updated each year, is detailed in Sitch et al. (2024).

1.4DONE: Climate forcing¶

ORCHIDEE requires eight meteorological variables at a default, but configurable, reference height of 2 m or 10 m above the plane of zero displacement when the stand-alone configuration is used or the lowest atmospheric layer when the land-atmosphere configuration is used. The eight meteorological variables are: air temperature (K), incoming direct shortwave and diffuse longwave radiations ($\text{W\,m}^{-2}$), liquid and solid precipitation (mm), air specific humidity ($\text{g\,g}^{-1}$), surface wind ($\text{m\,s}^{-1}$, with the possibility to provide latitudinal and longitudinal components), and surface pressure (Pa).

The geographical domain of the forcing data has to be provided on a standard latitude and longitude coordinate system. It can be any spatial domain ranging from a single point up to the entire globe. In order to reduce the disk space used by the forcing data and as ORCHIDEE only runs over land points, the two-dimensional grid can be reduced to a subset of land points. This means that only the forcing data for land points are stored in the forcing files, together with an indexing table that allows scattering the grid points back onto the regular latitude longitude grid when needed.

Although the stand-alone configuration, which makes use of a climate forcing, runs at the half-hourly time step, different climate forcings could, depending on their source, come at a half-hourly to a daily time step. If the frequency of the forcing is less than 6-hourly, ORCHIDEE interpolates the climate forcing between two data time steps. Instantaneous fields and long-wave radiation are linearly interpolated for each model time step $t$, typically half-hour, based on the values available at $t_0$ and $t_0+dt$ where $t_0$ < $t$ < $t_0+dt$. Precipitation fields are spread over the interval defined by a user defined parameter, so that this interval is the one over which the precipitation lasts when the forcing interval has rain or snow. Shortwave radiation is interpolated using a function distribution that corresponds to the solar angle distribution over a forcing time period while conserving the average flux for the interval [$t_0$, $t_0+dt$].

If the frequency of the climate forcing is greater than 6-hours, ORCHIDEE needs to reconstruct the diurnal cycle. The solar angle is calculated at each model time step—typically $t$ along with the times of sunrise, sunset, and solar noon. Based on these solar parameters and the concurrent time information, solar radiation is interpolated from its daily mean, and the corresponding air temperature is estimated using the daily max and min values. In oRCHIDEE v4.2 the diurnal variation of precipitation is reconstructed using a numerical weather generator. On any rainy or snowy day, the duration of precipitation is initially set to 2 hours and then adjusted based on air temperature. If the temperature is below 20 ° C, the duration increases up to a maximum of 8 hours. when the temperature exceeds 20 ° C, the duration remains fixed at 2 hours. Subsequently, the timing of precipitation within the day is random where the randomisation uses air temperature as the seed to ensure reproducibility. Due to the randomisation, the frequency and intensity of precipitation do not necessarily reproduce the local precipitation characteristics, e.g., drizzle at the global west coasts. Subsequently, wind speed, air pressure, and humidity are linearly interpolated between the days.

ORCHIDEE v4.2 uses CRU-JRA as its default climate forcing. It is based on the Japanese Reanalysis data Kobayashi et al., 2015 aligned with the CRU TS data from meteorological stations Harris et al., 2020Harris et al., 2014. The realignment preserves the monthly means of the CRU TS dataset and concerns temperature (Tmin, Tmax, Tmean), vapor pressure and precipitation, as detailed in Sitch et al. (2024). The forcing is defined on a 0.5 ° regular grid at 6-hourly time steps that cover the time period from 1901 to 2024. Regridded CRU-JRA forcing to 2 ° resolution is used for systematic ORCHIDEE reference simulations.

Other forcing datasets prepared to be used with ORCHIDEE include:

• GSWP has global coverage with a spatial resolution of 0.5 ° on a regular grid, and a temporal resolution of 3 hours from 1901 to 2010 Dirmeyer et al., 2006;

• WFDE5-CRU-GPCC has global coverage with a spatial resolution of 0.5 ° on a regular grid and a temporal resolution of 1 hour from 1980 to 2018 Cucchi et al., 2020;

• ERA5 has global coverage with a spatial resolution of 0.25 ° on a regular grid and a temporal resolution of 1 hour from 19XX to 20XX Bell et al., 2021;

• ERA5Land has global coverage with a spatial resolution of 0.1 ° on a regular grid and a temporal resolution of 1 hour from 19XX to 20XX Muñoz-Sabater et al., 2021;

• SAFRAN covers the metropolitan area of France with a spatial resolution of 8 km on a Lambert II extended grid with a temporal resolution of 1 hour from 1959 to 2023 Vidal et al., 2010;

• CRU-JRA55 has a global coverage with a spatial resolution of 0.5 degree on a regular grid and a temporal resolution of 6 hours from 1901 to 2022 Harris et al., 2020.

The atmospheric forcing is required only when running the ORCHIDEE model in the so-called offline mode. When coupled to the atmospheric model (mainly LMDz in case of ORCHIDEE), the atmospheric forcing is not needed.

1.5DONE: Lake properties¶

The lake fraction per grid cell is based on the HydroLakes database Messager et al., 2016. ORCHIDEE reads in an additional file containing spatialized values of effective lake depth, wind fetch, water albedo and extinction coefficients for shallow (< 5 m), medium (between 5 and 25 m), and deep lakes (> 25 m). Both files are available at a 0.25 ° and a 0.5 ° resolution for ORCHIDEE v4.2.

The HydroLakes database maps 1.4 million lakes of size larger than 0.1 km$^2$ on the global scale and documents their main properties, such as surface area and average depth. To generate the ORCHIDEE lake parameters files, all lakes available in the HydrolLake database were clustered into the three depth classes, i.e., shallow, medium and deep. For each grid cell and each depth-class, the total surface area was calculated as well as a weighted average mean of their respective depth and wind fetch which was then used to calculate the effective lake area, depth and fetch of the each depth-class separately. The albedo and extinction coefficient of freshwater are currently prescribed as 0.07 and 1 $m^{-1}$, respectively.

1.6DONE: Routing graphs¶

For the river routing scheme, ORCHIDEE reads information from a hydrological digital elevation model (HDEM). The minimal information required is elevation, flow direction, flow accumulation and distance to the ocean for each pixel. Ideally, the elevation in the input file is hydrologically consistent, in the sense that no flow direction should lead water to gain elevation. Several datasets that fulfil these criteria are commonly used with ORCHIDEE:

• STN-30p at 0.5° resolution, from Vörösmarty et al. (2000). It has been enhanced by flow directions in Antarctica in order to carry meltwater to the ocean and close the global water cycle. It is thus commonly used for global simulations.

• MERIT at 60 arcsec resolution, from Yamazaki et al. (2019). It is a global dataset, excluding Antarctica, but its high resolution makes it more suitable for regional simulations.

• HydroSHEDS at 30 arcsec resolution, from Lehner & Grill (2013). It only covers regions below 60° N and without Antarctica, and is also used for regional simulations.

These HDEM maps are then processed to provide the two minimal elements needed at each of the ORCHIDEE routing cells to create a routing graph:

1. A single water flow direction among 11 possibilities: 8 directions towards another routing cell (N, NE, E, SE, S, SW, W, NW), one direction towards the endorheic lakes - local inland convergences of the routing graph (also called lake inflow in ORCHIDEE), one direction towards the ocean from the main stream (riverflow) and one direction accounting for the small disperse flows into the ocean (coastalflow).

2. A value of the topographical water retention index $k^{Wat-TopoIndex}$ (in km), computed as:

   $$k^{Wat-TopoIndex}= \sqrt{\frac{d^3}{\Delta z}}\cdot 10^{-3}$$

   (3)

   with stream length $d$ (in m) assumed as the distance to the downstream routing cell, and $\Delta z$ (in m) the elevation drop between the two cells.

Beside this minimal information for river routing, the routing graph files can be completed by information characterizing natural hydrological elements (flood plains, swamps, ponds) or anthropogenic infrastructures (dams, reservoirs, adduction channels, gauge stations) as described below Schrapffer et al., 2023Polcher et al., 2023.

1.7Reservoirs and irrigation¶

for irrigation, floodplains, and ponds (we heave a flag doponds) refer to https://forge.ipsl.fr/orchidee/wiki/Documentation/Ancillary as well as Pedro and Patricia’s articles

For the irrigation scheme introduced in Arboleda-Obando et al. (2024), two input files are required, which the irrigation scheme can update the map for every simulated year.

The first one is a map of irrigated fractions, which can be obtained from the Historical Irrigation Dataset HID, Siebert et al., 2015, as in Arboleda-Obando et al. (2024)Tiengou et al. (2025). HID provides a map every 10 years before 1980 and every 5 years after, at 5 arcmin resolution. To run future climate simulations, Arboleda-Obando et al. (2025) used irrigated fractions from the Land Use Harmonization 2 dataset LUHv2, Hurtt et al., 2020. The dataset was used in the CMIP6 framework with historical and SSPs scenarios, and provides data at 0.25° resolution.

The second input map describes the available equipments for irrigation withdrawals, to preferentially withdraw water from the surface reservoirs (overland and rivers) or groundwater reservoir of the routing scheme. This map is derived from the inventory of areas equipped for irrigation of Siebert et al. (2010).

1.8Slope¶

Slope is required as input to the hydrology module to constrain the re-infiltration of surface runoff. A map at 15 arc-min resolution (1/4°) from the US Geological Survey is used. It contains... COMPLETE.

1.9DONE: Nitrogen inputs¶

With the implementation of the nitrogen cycle in ORCHIDEE, the nitrogen deposition rate of mineral nitrogen and the application of fertiliser onto the land surface are required as input to the model. Nitrogen fertiliser input datasets are taken from the NMIP2 project Hanqin et al., 2022. In the standard configuration we take the most recent update of gridded N application rates from the TRENDY project Sitch et al., 2024.

It includes inorganic nitrogen fertiliser application, which only started after the Haber-Bosch process was developed in the early 20th century, and manure application. Since manure is organic nitrogen that comes from waste from animals that ate vegetation, it should in principle be taken from other organic nitrogen sources in the model, which is not done yet in ORCHIDEE v4.2. For nitrogen deposition, the historical deposition rates are supplied as a time-varying spatially-varying deposition rate provided by the TRENDY project Sitch et al., 2024. The time series starts in 1850 and we thus use the 1850 gridded values for any year prior to 1850. The nitrogen deposition maps combines deposition of NO_y and NH_x species.

BNF is still missing

1.10DONE: Forest management¶

Historical simulations use a spatially explicit global reconstruction of dominant forest management strategies. The reconstruction is performed on a regular 0.25 ° x 0.25 ° grid between the years 1700 and 2022. The reconstruction distinguishes two forest management strategies: (1) unmanaged, and (2) managed. Although scientific forest management was rare to none existent before 1750 Perlin, 2005, management was thought to be an acceptable proxy for historical forest use where stand density was also reduced.

As the forest management maps have to be consistent with the vegetation distribution in ORCHIDEE, the annual vegetation maps are the basis of the management maps. All forests located in Greenland, Canada, the United States, and Russia were assigned to the unmanaged class unless the LUH2v2 reconstruction considered the pixel a secondary forest Hurtt et al., 2011. All forests outside of Greenland, Canada, the United States, and Russia were assigned to the managed class unless the LUH2v2 reconstruction considered the pixel to be primary forest Hurtt et al., 2011. All forest PFTs within a grid cell are assigned to the same forest management strategy.

1.11DONE: Albedo background¶

A background albedo map for the visible (VIS) and near infrared (NIR) parts of the solar spectrum is used to compute the overall grid-cell albedo (See 1.2). The background maps are derived from an optimization process (See 1.3.3) that uses MODIS VIS and NIR surface albedos and the background albedo maps of the Joint Research Centre Two-stream Inversion Package (JRC-TIP; Pinty et al. (2011)) as prior information.

1.12DONE: Litter raking¶

European maps of litter demand were based on historical livestock estimates Krausmann et al., 2013, taken to be equal to 0.6, 0.5, and 0.3 head of livestock person$^{-1}$ for northern, central, and southern Europe, respectively. The dividing parallels between northern, central and southern Europe were taken to be respectively 55 and 45° N latitude. It has been reported that 200 to 480 kg of dry litter were collected per livestock unit per year Bürgi, 1999. It was assumed that 480 kg litter per livestock unit per year corresponded to the peak demand in the mid 1800s and faded out afterwards.

Litter maps were generated from the livestock density maps using peak demand for all years from 1600 to 2010. Next, these initial litter maps were multiplied by a correction factor to account for the temporal evolution in litter demand McGrath et al., 2015. The correction factor was tuned to give the desired behaviour, based on historical information from Switzerland Bürgi, 1999.