2.1 Climate of Chile

Chile's location bordering the Pacific Ocean, its vast meridional extent of 4345 km, and its steep topographic longitudinal profile (Fig. 1b) result in highly variable climatic conditions. The large-scale subsidence of air masses over the southeast Pacific Ocean and other regional factors (e.g., relative cold coastal ocean currents) yield extremely arid desert conditions in northern Chile with as little as 2–20 mm of annual precipitation (Garreaud and Aceituno, 2007). In the south, a 1000 km long narrow band of Mediterranean-type climate exists on the western side of the Andes (Armesto et al., 2007). According to a general bioclimatic classification by Luebert and Pliscoff (2017), the tropical/Mediterranean boundary is located from 23^∘ S (coast) to 27–28^∘ S (inland) and the Mediterranean/temperate boundary is located at 36^∘ S (in both the Chilean Coastal and Andes mountain ranges) to 39^∘ S (the “Central Depression”), while the temperate/boreal boundary is found from 50.5 to 56^∘ S. From a climatologic standpoint, the Mediterranean bioclimatic region has a warm temperate climate dominated by winter rain (the mean annual precipitation varies from 300 to 1500 mm from north to south, respectively) and hot, dry summers with dry periods varying from 7 months (north) to less than 4 months (south) in duration (Uribe et al, 2012). To the south, midlatitude westerly winds and orographic uplift by the coastal mountains and the Andes lead to an annual precipitation of up to 3000 and 5000 mm, respectively (Veblen et al., 1996). El Niño occurrences generally lead to above average precipitation rates in the austral winter and spring in the Mediterranean zone and reduced precipitation at 38–41^∘ S in the following austral summer (Garreaud and Aceituno, 2007; Montecinos and Aceituno, 2003).

Table 1Major vegetation zones, associated plant functional types (PFT), and representative species (modified after Escobar Avaria, 2013; Luebert and Pliscoff, 2017; the letters a–e are given to help the reader identify the PFT and the representative species).

Download Print Version | Download XLSX

2.2 Vegetation of Chile

Considering the ecological divisions of South America (Young et al., 2007), Chile is represented by three noticeable areas: the Peruvian–Chilean desert, Mediterranean Chile, and the moist Pacific temperate area. These areas reflect the main floristic characteristics and vegetation types of the country; moreover, the distribution of vegetation within these areas is constrained by thermal and hydrological climatic factors that vary according to latitude and longitude (see Table 1 for an outline of the characteristic species, the biome classification, and modeled vegetation types). In the north, longitudinal variations in climate are the result of geomorphological and ombroclimatic changes. Between 17 and 28^∘ S, the coastal zone is exposed to the influence of fog and orographic precipitation, allowing vegetation such as columnar cacti (Eulychnia genera) and a diverse group of shrubs (e.g., Nolana, Heliotropium, Euphorbia, Tetragonia) and succulents (e.g., Deuterocohnia, Tillandsia, Puya, Neoporteria) to develop (Luebert and Pliscoff, 2017). To the west, and between the Chilean Coastal and the Andes mountain ranges, a hyper-arid desert zone exists, which is characterized by the absence of rainfall or coastal-fog water inputs; therefore, no vascular plants are usually found in this area. However, with punctual water inputs due to local substrate conditions, such as the presence of a water table, some halophytic shrubs (e.g., Prosopis sp.) can develop (Luebert and Pliscoff, 2017). Vegetation in the Andes mountain range is constrained by altitude due to the decrease in temperature and the increase in precipitation. The major development regarding vegetation is observed at intermediate altitudes where shrubs tend to dominate (Luebert and Pliscoff, 2017).

Vegetation increases to the south with increasing winter precipitation (Rundel et al., 2007), allowing Mediterranean-type shrubland and woodland ecosystems to develop. In these ecosystems various plant species, commonly denominated sclerophyllous plants, have small, rigid, xeromorphic leaves adapted to hot, dry summers and wet, cool winters (Young et al., 2007). Sclerophyllous woodlands and forests extend from 30–31 to 37.5–38^∘ S (Fig. 1a) and range from xeric thorn savanna elements to dense herbaceous cover in the Central Depression, whilst evergreen sclerophyllous trees and tall shrubs are found at higher elevations in which mesic conditions predominate. In these ecosystems, slope aspect can modify local moisture conditions, affecting the structure and composition of vegetation (Luebert and Pliscoff, 2017). To the south, these vegetation types transition into temperate deciduous Nothofagus forests (Maule or Nothofagus parklands. The dominant species in these areas being Nothofagus obliqua, N. glauca, and N. alessandrii) at the Coastal Cordillera and lower Andes ranges, before forming a broader vegetation zone (Donoso, 1982; Villagrán, 1995). The deciduous Nothofagus forests then grade into mixed deciduous–evergreen Nothofagus forests at 36^∘ S (Young et al., 2007). With increasingly hydric conditions evergreen broad-leafed species begin to dominate forest stands at approximately 40^∘ S and form the Valdivian rainforest (the northernmost rainforest type, ranging from 37^∘45^′ to 43^∘20^′ S) with high biomass and arboreal biodiversity and evergreen, deciduous, and needleleaf species (Veblen, 2007). Further south, the less diverse North Patagonian rainforest is mainly dominated by Nothofagus betuloides (Veblen, 2007) and transitions into the Magellanic rainforest at approx. 47.5^∘ S and Magellanic moorland with waterlogged soils and a poor nutrition status at the coast (Arroyo et al., 2005). Cold deciduous forests stretch from 35 to 55^∘ S along the Andes covering cooler and dryer sites as compared to the coastal rainforests. These forests occur at altitudes of approx. 1300 m and gradually descend to sea level in Tierra del Fuego (Pollmann, 2005). In Tierra del Fuego and east of the low Andes in southern Patagonia a gramineous steppe exists (Moreira-Muñoz, 2011), and a high-Andean steppe also extends north at higher altitudes.

3.1 Vegetation model

The Lund-Potsdam-Jena General Ecosystem Simulator (LPJ-GUESS; Smith et al., 2001, 2014) is a state-of-the-art dynamic vegetation model that also simulates detailed stand dynamics using a gap-model approach (Bugmann, 2001; Hickler et al., 2004). The model is developed by an international community of scientists and has been used in more than 200 peer-reviewed international publications, including model evaluations against a large variety of benchmarks such as vegetation type distribution, vegetation structure and productivity, as well as carbon and water cycling at regional to global scales (https://www.nateko.lu.se/lpj-guess, last access: 15 September 2018). Vegetation development and functioning is based on the explicit simulation of photosynthesis rates, stomatal conductance, phenology, allometric calculations, and carbon and nutrient allocation. The model simulates the growth and competition of different plant functional types (PFT, see Bonan et al., 2002) based on their competition for space, water, nutrients, and light. Population dynamics are then simulated as stochastic processes that are influenced by current resource status, life history, and demography for each PFT. To enable a representative description of average site conditions within a landscape each grid cell is simulated as a number of replicate patches in order to allow different (stochastic) disturbance histories and development (successional) stages (see Hickler et al., 2004; Wramneby et al., 2008). Fire occurrence is determined by the model using temperature, fuel load, and soil moisture levels (Thonicke et al., 2001). Soil hydrology in LPJ-GUESS is represented by a simple two-layer bucket model with percolation between layers and deep drainage (see Gerten et al., 2004).

Vegetated surface area and runoff are affected by a range of model parameters in LPJ-GUESS. In “cohort” mode, an average individual from each PFT with a given age and development status is used to characterize vegetation state. Depending on PFT-specific parameters (i.e., maximum crown area, sapling density, allometric properties, leaf-to-sapwood area), age, and competition for light, water, space, nutrients, and demographic processes (establishment and mortality) individual cohorts can develop different states.

Using Eq. (1), we approximate the fraction A of the land surface covered by vegetation using the foliar projected cover (FPC) – the vertical projection of leaf area onto the ground (see Wramneby et al., 2010). In LPJ-GUESS the FPC is derived from daily leaf area index (LAI, the leaf area to ground area ratio,  m² m⁻²) summed for all simulated PFTs (n refers to the number of PFTs) using the Lambert–Beer extinction law (originally proposed by Monsi and Saeki in 1953 for the estimation of light extinction in plant canopies, see translation in Monsi and Saeki, 2005; Prentice et al., 1993):

Thus, depending on the composition of PFTs and the disturbance regime, varying levels of ground cover can be simulated that not only reflect the environmental conditions but also vegetation diversity and development.

In addition, the hydrological cycle is also affected by PFT-specific interception and transpiration rates that are a function of PFT-specific parameters and the development stage. Thus, vegetation modulates infiltration via interception (that is a function of vegetation cover) and runoff (via plant uptake and transpiration of water) under the given environmental constraints. In LPJ-GUESS water enters the top soil layer as precipitation until this layer is fully saturated (excess water is lost as surface runoff, and evaporation removes water from a 20 cm sub-horizon of the top layer). During precipitation days, water can percolate from the top to the lower layer until the lower layer is saturated (excess water is lost as drainage). In addition, water from the lower layer can drain as baseflow with a fixed drainage rate (Gerten et al., 2004; Seiler et al., 2015). The model does not consider lateral water movement between grid cells nor does it take routing in a stream network into account (in this study we report the surface runoff component only).

3.2 Landform classification

To bridge the gap in spatial resolution between LPJ-GUESS (typical spatial resolution of $\mathrm{0.5}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$) and the LEM Landlab (Hobley et al., 2017; Schmid et al., 2018, typical spatial resolution ∼100 m) and facilitate the future coupling of these models, we introduced the concept of landform disaggregation of grid cell conditions to smaller sub-pixel entities (Figs. 2, B1). The advantages of introducing sub-pixel entities as opposed to simply performing higher-resolution simulations are twofold. Firstly, higher-resolution simulations require climate forcing data of the desired output resolution that are not available for the intended simulation periods or region (and are not generally available for past time periods). Secondly, using sub-pixel entities incurs smaller additional computation costs than higher-resolution simulations.

Figure 2Schematic procedure of simulations. Coarse resolution model driving data (1) is disaggregated using a high-resolution elevation model and topographic landform classification (2), and the ecosystem model LPJ-GUESS then simulates vegetation state and dynamics using the landform classification to simulated topographic-adjusted patch composition (3). Vegetation cover and surface runoff results can then be passed on to a coupled landscape evolution model (LEM) (not implemented in this study, see Schmid et al. (2018) for a description of the LEM).

Download

Sub-pixel entities, hereafter termed “landforms”, were derived for each grid cell using SRTM1-based elevation models (Kobrick and Crippen, 2017). Pixels from the elevation model (30 m spatial resolution) were classified based on their elevation (200 m bands) and their association with topographic features (ridges, mid-slope positions, valleys, and plains – based on slope and aspect), and similar pixels were grouped to form the landforms. The elevation of the landforms was used to modify the temperature at a given landform. Using the elevation difference between the average elevation of a landform as derived from the high-resolution elevation model and the reference elevation obtained from the $\mathrm{0.5}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$ grid, a temperature delta was calculated using the lapse rate from the International Standard Atmosphere model of −6.5 ^∘C km⁻¹ (e.g., Vaughan, 2015). However, it should be noted that this lapse rate is a global average rate that can substantially differ from local conditions and over multiple timescales (ranging from sub-daily to climatological), as it is controlled by various atmospheric thermodynamics and dynamics (i.e., radiative conditions, moisture content, large-scale circulation conditions). While a higher lapse rate would potentially be a better approximation for drier sites (e.g., Pan de Azúcar), this might not be the case for other sites or past time periods with different atmospheric conditions. Furthermore, the lack of defining environmental data and a mismatch of scales prohibits the use of more specific regional lapse rates in our simulations.

The slope and aspect were utilized to adjust the incoming radiation received by the landform (see Appendix B). The general topographic features were used to modify the depth of the lower soil layer (deeper soils in valleys and on plains, shallow soils on ridges) and to identify areas (valleys and plains) with a newly implemented time-buffered deep-water storage pool that is only accessible by tree PFTs. The classification resulted in 2 to 56 (the mean was 17) landforms for each grid cell depending on the topographic complexity.

In the proposed future coupled model, we envisage that the landform classification will be performed using elevation information from the LEM. The resulting per-landform (but non-spatially explicit) vegetation simulation results will be matched back to spatially explicit grid cells in the LEM, which will allow us to bridge the scale gap between the two models.

In the simulations presented here, all classified landforms with an area > 1 % of the total land area in a grid cell were simulated using 15 replicate patches each, and simulation results were aggregated by area-weighting the results to the grid cell level. For a summary of the implementation details see Fig. B1.

3.3 Parameterization of plant functional types and biome classification

In a previous study Escobar Avaria (2013) implemented the first regional simulation for Chilean ecosystems (also using LPJ-GUESS) for present day (PD) climate conditions using a region-specific parameterization, which the presented study adapts and builds upon. Eleven PFTs – three shrub types, seven tree types, and one herbaceous type – were defined in order to describe the major vegetation communities of Chile (Table C1). The definition of these PFTs are generally based on the proposed macro-units of Chilean vegetation (Luebert and Pliscoff, 2017) and follow the concept of representative/or dominant species for describing a physiognomic unit. Apart from growth habit and associated traits, PFTs were designed to differentiate between leaf morphology and strategy, shade-tolerant and shade-intolerant varieties, their adaption to water access (mesic, xeric type), and root distribution. An overview of major eco-zones, associated PFTs, and representative species is given in Table 1.

Using a biomization approach (see Prentice and Guiot, 1996 for the general concept), we classified the simulated PFTs into discrete vegetation types (referred to as biomes throughout the paper; see Fig. C1 for details of the classification procedure). Based on LAI thresholds and the ratio of certain key PFTs or PFT groups (i.e., boreal tree PFTs, xeric PFTs) to their peers, a cascading decision tree was implemented that led to 11 biomes resembling the general vegetation zones of Chile (Fig. 1a), but also included additional biome classes to capture the finer nuances of transitions between semi-arid and mesic vegetation communities and open woodlands. To keep the number of vegetation types reasonable we designed multiple decision paths for biomes that exist as dense forest ecosystems but also transition into lower canopy woodlands or transition into more open woodlands (i.e., Magellanic forest, cold deciduous forest; see Fig. C1). The classification was conducted at the landform level after the simulated 15 patches were averaged, and the grid cell classification was derived from picking the area-dominant class of the landforms of a grid cell.

Figure 3(a) Average annual temperature and (b) precipitation derived from the downscaled and bias-corrected TraCE-21ka transient paleoclimate data (Liu et al., 2009) for the Last Glacial Maximum (LGM), mid-Holocene (MH), and present day (PD) time slices (data is the average of 30 year monthly data; 1 represents Pan de Azúcar, 2 represents Sta. Gracia, 3 represents La Campana, and 4 represents Nahuelbuta).

Download

3.4 Environmental driving data and modeling protocol

The climate forcing data for LPJ-GUESS is derived from TraCE-21ka (Liu et al., 2009), which is a transient coupled atmosphere–ocean simulation from LGM to PD using the Community Climate System Model version 3 (CCSM3; Collins et al., 2006). In this study we present time-slice simulations for the LGM, MH, and PD (1960–1989), using perpetual climate forcing data from 30-year monthly climatologies (Fig. 3). In addition, we show results from a transient model simulation that utilizes the full time-series data from the LGM to PD. All simulations were preceded by a 500-year spin-up period with de-trended climate data until vegetation and soils reached steady state. For time-slice simulations the last 30 years of the simulation were used.

Monthly temperature, precipitation, and downward shortwave radiation from the TraCE-21ka dataset (resolution T31; $\sim \mathrm{3.7}{}^{\circ }$) were downscaled to a $\mathrm{0.5}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$ spatial resolution and bias-corrected using a monthly climatology from the ERA-Interim reanalysis (Dee et al., 2011, years 1979–2014). We used an additive bias correction for the temperature and multiplicative corrections for the precipitation and shortwave radiation (see Hempel et al., 2013; this technique was also used in, e.g., O'ishi and Abe-Ouchi, 2011). The multiplicative corrections for the precipitation and radiation are necessary because these fields can not have negative values. The resulting bias-corrected anomalies were subsequently downscaled to the ERA-Interim grid using a bilinear interpolation technique. The number of rain days within each month (used by LPJ-GUESS to distribute monthly precipitation totals to daily time steps internally) was derived from the monthly mean precipitation in the TraCE-21ka data and the day-to-day precipitation variability from ERA-Interim (see Appendix A). The [CO₂] for each simulation year was obtained from Monnin et al. (2001) and Meinshausen et al. (2017). A comparison climate dataset for the LGM (referred to as ECHAM5 in this paper) was provided by Mutz et al. (2018). Soil texture data used for bare-ground initialization of the model was obtained from the ISRIC-WISE soil dataset (Batjes, 2012) and a default soil depth of 1.5 m (0.5 m topsoil, 1 m subsoil) was assumed.

We simulated vegetation dynamics using the process-based dynamic vegetation model LPJ-GUESS (Smith et al., 2001) version 3.1 (Smith et al., 2014) with the model additions outlined above. The model runs were carried out without nitrogen limitation, using the CENTURY carbon cycle model (see Smith et al., 2014). Patch destroying disturbance and establishment intervals were defined as 100 and 5 years, respectively, and fire dynamics were enabled. Further details of the PFT specific parameterizations are given in Table B1. The transient site-scale model runs were only conducted for the four focus sites of the EarthShape SPP due to (i) computing resource constraints, (ii) better comparability with other EarthShape SPP work (see Schmid et al., 2018), and (iii) better interpretability.

4.1 Distribution of simulated plant functional types

Vegetation communities (expressed as assemblages of PFTs in LPJ-GUESS) establish spatially depending on (a) environmental controls, (b) competition, and (c) stochastic events (i.e., fire incidents and mortality). An overview of the simulated vegetation distribution expressed as the simulated LAI for each PFT under PD climate conditions is given in Fig. 4 (for key PFT properties see Table C1; the PFT distribution maps for the MH and the LGM are given in the Supplement as Figs. S1 and S2 for completeness). Temperate broad-leaved evergreen trees (TeBE_tm, TeBE_itm; t represents shade-tolerant, it represents shade-intolerant, and m represents mesic) dominate the coastal and central areas from latitudes 40 to 46^∘ S but also extend north into the Mediterranean zone. Further to the north, temperate broad-leaved summergreen PFTs (TeBS_tm, TeBS_itm) occur, with the shade-tolerant type dominating a relative small area between 37 and 40^∘ S. Northward, and in coastal areas, the sclerophyllous temperate evergreen PFT (TeBE_itscl; scl represents sclerophyllous) starts to dominate, and with dryer conditions the total LAI (LAI_tot) is dominated by evergreen and raingreen shrubs (TeE_s, TeR_s; s represents shrub). South of 40^∘ S, and in higher terrain as well as further north, boreal broad-leaved summergreen and evergreen tree PFTs (BBS_itm, BBE_itm) as well as boreal evergreen shrubs (BE_s) establish. In the Andean ranges and (to a lesser extent as secondary components) in the lowlands and coastal ranges from 35 to 50^∘ S, temperate needle-leaved evergreen trees (TeNE) are simulated. Herbaceous vegetation (C3G) dominates LAI_tot at high altitudes along the Andean ranges and is also a substantial contributor to the total LAI in the Mediterranean zone (30 to 38^∘ S). Herbaceous vegetation also contributes to a lesser extent in most other regions, except for in the hyper-arid desert areas. LAI_tot is highest in the zone from 36 to 50^∘ S. While LAI decreases substantially at sea level towards the Atacama Desert, higher values of the LAI are found further north at higher altitudes (see inset in Fig. 4).

Figure 5(a) Spatial and (b) altitudinal distribution of biomes for Last Glacial Maximum (LGM), the mid-Holocene (MH), and the present day (PD) (for the biome classification decision tree see Fig. C1).

Download

4.2 Distribution of simulated biomes at LGM, MH, and PD

The simulated biome distribution changes spatially (Fig. 5a), vertically (Fig. 5b), and over time. Under PD climate, the simulated Valdivian temperate rainforests extend from 38 to 46^∘ S at the coast and transition into the Magellanic forests/woodlands, that are dominated by the boreal PFTs. A small zone of deciduous “Maule” forest occurs to the north of the Valdivian rainforest at meso-temperate climates. With even dryer and warmer conditions the sclerophyllous woodland type establishes and (as the fraction of trees and LAI_tot is reduced even further with increasing temperatures and even lower annual rainfall) eventually give way to shrub dominated matorral and finally the arid shrubland biome type. The cold deciduous forest type is classified for parts of Tierra del Fuego, and at higher elevations in the lower Andes in Patagonia. It also forms larger zones at altitude between 30 and 40^∘ S. Mesic woodland occurs between 34 and 30^∘ S at altitude (above the sclerophyllous woodland zone dominating the lowlands), and high-Andean steppe occurs between 18 and 30^∘ S. A cold desert is present above the tree line in Patagonia and the highest Andean ranges, whereas hot desert (LAI_tot < 0.2) is simulated for areas from 20 to 26^∘ S. The model was able to simulate the general distribution of the biomes of Chile for most regions (Figs. 1a, 5), although accuracy for the occurrence of deciduous PFTs was low (deciduous “Maule” forest and cold deciduous forest, also previously reported by Escobar Avaria, 2013).

Owing to the similarity in climatic conditions (Fig. 3), the distribution simulated for the MH does not differ substantially from the PD (Table 2). The northern border of sclerophyllous woodland shifts to approx. 34^∘ S, giving way to a matorral zone. In addition, the cold deciduous forest biome covers larger areas in Patagonia at the expense of Magellanic woodland (+27.5 % and −9.1 %, respectively; Table 2). However, the spatial and vertical distribution of biomes for the LGM is markedly different (Fig. 5a, b). The substantially lower temperatures (Fig. 3) lead to an expansion of cold deserts up to 45^∘ S (coastal areas) and 40^∘ S (higher altitudes), respectively. The boreal PFT dominated Magellanic woodland biome is substantially reduced in extent (16.5 % in the LGM versus 32.6 % in the PD of simulated area, Table 2) and is shifted further north (40 to 45^∘ S at the coast, 43 to 35^∘ S at higher altitudes inland). The area covered by temperate rainforest is restricted to a small lowland area from 36 to 40^∘ S, and the larger areas of cold deciduous forest at altitude are also substantially smaller (Fig. 5b; Table 2). Lowland steppe and mesic woodland biomes are simulated instead of matorral and sclerophyllous woodlands, and desert cover larger areas of the high-Andes to the north.

Table 2Areal extent of biomes in the simulation domain (units: percent of total area).

Download Print Version | Download XLSX

Figure 6Spatial distribution of (a) foliar projected cover and (c) surface runoff simulated for the Last Glacial Maximum (LGM), the mid-Holocene (MH), and the present day (PD). Difference plots of LGM versus  PD (left) and MH versus PD (right) for foliar projected cover (b) and runoff (d), respectively (the number 1 represents Pan de Azúcar, 2 represents Sta. Gracia, 3 represents La Campana, and 4 represents Nahuelbuta).

Download

4.3 Foliar projected cover and surface runoff

The percentage of ground covered, and thus shielded from strong denudation, and surface runoff (as a major driver of erosion rates) are both influenced by the composition and state of vegetation communities. Therefore, we evaluate the regional and temporal changes of these important variables as simulated by LPJ-GUESS for the LGM, MH, and PD time slices. The simulated LAI of PFTs can be aggregated and converted to foliar projected cover (see Eq. 1); thus, this allows us to estimate the surface covered by vegetation. However, it should be noted that this is only an approximation of true ground cover, as small-scale vegetation variations are not simulated in a location-specific fashion (spatial lumping effects are not considered for instance). Nevertheless, the implementation of a sub-grid landform scale was, in part, motivated to improve the models' prediction of smaller scale differences, as it should allow for the differentiation of different sub-grid conditions for the major landforms within one simulation cell (see Sect. 5.1 and Appendix B for further information). Low FPC clearly coincides with the distribution of hyper- to semi-arid biomes (Figs. 5a, 6a). Under PD climate conditions, average FPC for the semi-arid and Mediterranean biome types, including arid shrubland, matorral, and sclerophyllous woodland are 16 %, 35 %, and 66 % (Table 3), respectively, and cover increases southward with increases in annual precipitation rates (Fig. 3b). South of 35^∘ S, FPC values > 70 % are simulated for most grid cells except for high-altitude locations, the glacier fields of North Patagonia, and parts of the Magellanic moorland at the coast (see also Table 3).

Table 3Average foliar projected cover (FPC) and annual surface runoff for simulated biomes for present day (PD) and relative change within these PD areas for the Last Glacial Maximum (Δ LGM, LGM-PD) and the mid-Holocene (Δ MH, MH–PD).

Download Print Version | Download XLSX

Simulated FPC was lower than satellite-based estimates by the MODIS Vegetation Continuous Fields product (Dimiceli et al., 2015), likely due to methodological differences in foliar projected cover and total satellite-observed vegetation cover, but the general patterns were represented (Fig. 7). The most distinct regional discrepancy can be observed in coastal areas between 30 to 36^∘ S and the Andean highlands in the north of the model domain. For the entire model domain the mean average error (MAE) of the FPC was 11.9 %. The best agreements of 6.1 % and 6.3 % MAE were achieved for the arid shrubland and Valdivian rainforest biomes, while the discrepancies were greatest for mixed forest and cold deciduous forest biomes with 22.4 % and 22.1 % MAE, respectively (see Fig. 4a for biome locations; note that LPJ-GUESS often tends to underestimate the FPC of deciduous PFTs). A strong positive correlation between annual rainfall and FPC can be observed for all semi-arid, Mediterranean, and seasonally dry temperature eco-zones (Fig. 8a) and increases in annual rainfall lead to a strong rise of ground covered until 250–300 mm rainfall per year (arid shrubland: +14.5 %–16.3 % 100 mm⁻¹; matorral: +12.9 %–16.5 % 100 mm⁻¹; steppe: +13.5 %–17.8 % 100 mm⁻¹; Fig. 8a). At these levels, Mediterranean and temperate woodland biomes start to dominate but increases in precipitation only raise FPC by +2.7 %–3.3 % 100 mm⁻¹ and +5.4 %–6.4 % 100 mm⁻¹, respectively (sclerophyllous woodland and mesic woodland).

Figure 7Comparison of satellite-derived vegetation cover and simulated average foliar projected cover of potential natural vegetation for present day (data: MODIS MOD44B VCF v6, total vegetation, 2001–2016 average; Dimiceli et al., 2017; 1 represents Pan de Azúcar, 2 represents Sta. Gracia, 3 represents La Campana, 4 represents Nahuelbuta).

Download

Figure 8(a) Foliar projected cover (FPC) and (b) surface runoff as a function of annual average rainfall for the Last Glacial Maximum (LGM), the mid-Holocene (MH), and present day (PD) (“ASh” is arid shrubland, “CDF” is cold deciduous forest, “DMF” is deciduous “Maule” forest, “MFW” is Magellanic forest/woodland, “Mat” is matorral, “MeW” is mesic woodland, “MixF” is mixed forest, “SclW” is sclerophyllous woodland, “St” is steppe, “VRF” is Valdivian rainforest; temperate and boreal forest biomes (CDF, MixF, VRF, MFW) excluded from subplot (a) as they are also strongly dependent on temperature.

Download

A similar spatial pattern of FPC can be observed for the MH (Fig. 6a) and for the relationship to precipitation rates (Fig. 8a). No significant difference in FPC to PD conditions is apparent for the Atacama Desert or for most temperate forest areas from 36 to 42^∘ S. FPC in high-Andes locations in the north and large parts of the Magellanic forest and cold deciduous forest biomes in the south is 10 %–20 % lower than under PD climate, whereas in the Mediterranean zone FPC is reduced by 5 %–10 % (Fig. 6b, Table 3). The reduction of MH FPC (47–54^∘ S) coincides with a zone of lower temperature (> −0.5 ^∘C, see Fig. 6b) and the areal extent of TeBE PFTs (Fig. S1). Lower [CO₂] at MH compared to PD might also be the reason for a general reduction of biomass productivity; however regional differences of the other climate drivers (i.e., temperature, precipitation) do alter this general trend.

Due to lower temperatures and reduced precipitation rates at high altitudes and in the southern part of the country (Fig. 3), FPC during the LGM in these areas is simulated to be strongly reduced (< −30 %, cold desert conditions for large areas south of 46^∘ S); however, the FPC was lower for all areas of Chile (Fig. 6b, Table 3) which was likely also the effect of substantially lower [CO₂] (see Sect. 4.5 for a sensitivity analysis of the effect of the atmospheric CO₂ concentration). While the general correlation of FPC to precipitation can also be observed for the LGM (Fig. 8a), the variability in vegetation cover in mesic and xeric woodlands appears to be larger – indicating the potential for greater variability in erosion rates within the same biome. The strongest correlations between annual precipitation and FPC were observed for sclerophyllous woodland (adjusted r² values of 0.84, 0.91, and 0.9 for LGM, MD, and PD, respectively; p< 0.001).

Annual average runoff varies greatly from north to south (Fig. 6c) and coincides with annual precipitation (Figs. 3, 8b). For PD and MH climates, LPJ-GUESS simulated almost no surface runoff for arid and Mediterranean areas of Chile to approx. 32^∘ S (see also Table 3). Thus, correlation between precipitation and runoff in the steppe and arid shrubland biomes was low (adjusted r²: 0.16–0.27, 0.12–0.37, 0.11–0.35 for LGM, MH, and PD, respectively; p < 0.01). For all other systems (excluding matorral and mixed forests), correlation coefficients were generally > 0.85 for all time slices; p < 0.001).

Runoff rates gradually increase southward and reach their peak (> 2000 mm yr⁻¹) in areas of hyper-humid conditions along the Pacific coast. MH runoff rates are higher for areas of northern Patagonian and Magellanic rainforest (40–46^∘ S), but lower for coastal areas of Magellanic moorlands (Fig. 6d). LGM runoff rates are higher for most areas (Table 3) and especially south of 34^∘ S, with the strongest differences occurring from 40 to 46^∘ S. However, we want to highlight that the low temporal variability of the TraCE-21ka precipitation data likely leads to a substantial underrepresentation of episodes of high hygric variability (see discussion).

Figure 9Transient simulations for the (a) Pan de Azúcar, (b) Sta. Gracia, (c) La Campana, and (d) Nahuelbuta sites. The insets show the location within the simulated $\mathrm{0.5}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$ grid cell and the area and location covered by the landform plotted in this figure (results given in other figures are an area-weighted aggregation of individual landform simulation results). Panels 1–3 (from top to bottom): “T_avg” is annual average temperature, “Precip” is annual precipitation, “Fire RI” is fire return interval. Panel 4: “PFT” is plant function type; (grouped) PFT abbreviations: “TeSh” is temperate shrubs, “Grass” is herbaceous vegetation, “TeBEX” is sclerophyllous temperate broad-leaved evergreen trees, “TeBS” is temperate broad-leaved summergreen trees (hashed represents shade-intolerant), “BSh” is boreal evergreen shrubs, “TeBE” is temperate broad-leaved evergreen trees, “TeNE” is temperate needle-leaved evergreen trees, “BBS” is boreal broad-leaved summergreen trees, “BBE” is boreal evergreen trees, “Various” refers to other PFTs (LAI <0̇.05); plain represents shade-intolerant, hatched represents shade-tolerant, cross-hatched represents raingreen. Panel 5: “FPC” is foliar projected cover; trees and shrubs > 5 m are denoted using dark gray, herbaceous vegetation and small trees and shrubs are denoted using light gray. “Runoff” is simulated surface runoff; ^* Biomization: for a legend on the applied classification see Figs. C1 and 4 for a color-coded legend. All data smoothed using 100-year averaging.

Download

4.4 Transient changes of simulated vegetation from LGM to present

In this section we present transient simulation results for grid cells that contain the EarthShape focus sites (Fig. 9). The results are given for a single landform of the simulated grid cells in order to preserve successional transitions between PFTs that might otherwise be lost through averaging. The field site location and the represented area of the chosen landform within the grid cell are marked in the insets of Fig. 9a–d. The diversity of simulated vegetation cover for all landforms of the four EarthShape focus sites is illustrated in Fig. B2. The area-averaged mean of landform enabled LPJ-GUESS simulations closely matches the default model results at the Sta. Gracia and La Campana sites, but FPC at the Pan de Azúcar and Nahuelbuta sites is approx. 3 % greater than the default simulation result (Fig. B2a). However, results from all four sites show that the inter-landform variability of FPC varies by 10 %, 15 %, 25 %, and 18 % for the Pan de Azúcar, Sta. Gracia, La Campana, and Nahuelbuta sites, respectively.

Temperatures and [CO₂] start to increase at 18 000 BP and a marked pullback in temperatures during the Antarctic cold reversal (∼14 500 BP) is present in the TraCE-21ka data for all four sites (Fig. 9a–d). Annual precipitation at the Pan de Azúcar site (21.11^∘ S 70.55^∘ W, 320 m a.s.l.) is extremely low (38–40 mm, hyper-arid) for the entire simulation period (Fig. 9a). The annual average temperature for the landform presented increases from 13.9 ^∘C (LGM) to 16.8 ^∘C (PD). As a result of these arid conditions only evergreen and raingreen shrubs and herbaceous vegetation can establish, and LAI, and consequently FPC, remains very low (LAI_tot < 0.3). For most episodes of the simulation this location is classified as desert according to the implemented biomization scheme, with only two periods switching to another state (arid shrubland, LAI_tot > 0.2) from approx. 17 000 to 15 000 BP, and more permanently, the late Holocene. Fire return intervals (expressed as the number of years between fire incidents) fluctuate greatly as fuel production is substantially limited by low vegetation growth. No surface runoff is simulated and FPC ranges from 8 % (LGM) to 13 % (PD).

Climatic conditions for the Sta. Gracia site (29.75^∘ S 71.16^∘ W, 579 m a.s.l.) show a similar temporal progression from the LGM to PD, but average temperatures are lower and range from 11.9 ^∘C (LGM) to 14.9 ^∘C (PD) (Fig. 9b). Annual precipitation does not change substantially from the LGM to PD (152 versus 122 mm). Vegetation from the LGM to approx. 18 000 BP is dominated by herbaceous vegetation, but (raingreen) shrubs increase with rising temperatures, leading to a biome shift from steppe to matorral. Fire return intervals decrease with increased (arboreal) litter production due to the encroachment of shrubs but FPC only increases by 7 % from 33 % (LGM) to 40 % (PD). Simulated surface runoff is insignificant for most simulation years.

Temperatures at the La Campana site (32.93^∘ S 71.09^∘ W, 412 m a.s.l.) increase from 11.0 ^∘C (LGM) to 14.0 ^∘C (PD), but annual precipitation amounts decrease from 446 mm (LGM) to approx. 320 mm before slightly increasing again to 355 mm at PD. The simulated LGM vegetation for the La Campana site is dominated by herbaceous plants and small fractions of temperate evergreen deciduous trees mixed with small fractions of boreal shrubs which leads to a steppe biome classification with short episodes of mesic woodlands (Fig. 9c). With the decrease of annual precipitation and increasing temperatures (approx. 17 500 BP) sclerophyllous woodlands displace the deciduous trees and evergreen and raingreen shrubs start to appear. The fire-return intervals shorten in this phase and reach values of less than 15 years for the remainder of the simulation. Raingreen shrubs expand at approx. 12 000 BP and push back on herbaceous vegetation and, in part, evergreen shrubs. The LAI of sclerophyllous broad-leaved evergreen trees and shrubs further increases during the last 5000 simulation years, which leads to a shift in our biome classification from matorral (LAI_tot > 0.5) to sclerophyllous woodland (LAI_woody > 1). Despite pronounced changes in vegetation composition, FPC only increases from approx. 51 % (LGM) to 59 % (PD), which translates to a relative stable vegetation cover for these regions over time; thus, there is likely only a low impact of biome shifts on erosive processes.

Climatic conditions at the Nahuelbuta site (37.81^∘ S 73.01^∘ W, 1234 m a.s.l.) are markedly different to the three previously mentioned sites, as this location receives substantially higher annual precipitation throughout the time series (> 1200 mm) and average temperatures at this latitude and elevation are substantially lower (5.1 ^∘C for LGM and 8.6 ^∘C for PD, Fig. 9d). However, it can be seen that the landform is only representative of a small fraction of the $\mathrm{0.5}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$ grid cell as it is located on mountainous terrain, whereas most areas in the cell are covered by coastal lowlands with higher annual temperatures; thus, the site simulation results presented here differ from the total grid cell results presented in previous sections (see marked landform cover in inset; Figs. 9d and B2a, b). LPJ-GUESS simulates a diverse composition of PFTs and a transition from boreal, Magellanic woodland conditions at the LGM, to a period of cold deciduous forest (17 500–12 000 BP), followed by 12 000 years of Valdivian rainforest and mesic woodland alternations. During the LGM, boreal broad-leaved evergreen, deciduous tree, and shrub PFTs dominate and form a forest. Annual precipitation (> 1300 mm) and surface runoff (> 440 mm) is high during this period. With rising temperatures, the boreal shrubs and evergreen tree PFTs are displaced by temperate needle-leaved evergreen trees (TeNE) and increases in herbaceous vegetation. Temperate evergreen PFTs establish approx. 17 000 BP and, after another retreat at approx. 14 500 BP (coinciding with the Antarctic cold reversal), start to dominate the forest at this location. Fire frequency is low for the first 4000 simulation years and only rises to approx. 1 fire in 100 years afterwards. FPC remains at constantly high values (> 75 %) indicating a largely closed forest for the entire simulation period and thus low dynamics of erosive processes due to constant high vegetation cover in this area.

4.5 Sensitivity of foliar projected cover to [CO₂]

The effectiveness of photosynthesis, and thus the plants ability to build up biomass, has a large dependence on [CO₂] (Farquhar et al., 1980; Hickler et al., 2015). Increases in vegetation biomass (expressed in LAI) were observed in our simulations but coincide with a simultaneous rise in temperatures and [CO₂] (Fig. 9b–d). To assess the direct effect of changes in [CO₂], we conducted a sensitivity simulation under LGM climate conditions. We compared the default LGM simulations (TraCE-21ka, [CO₂] $=\sim \mathrm{180}$ ppm) with a preindustrial [CO₂] level of 280 ppm (Fig. 10). This change led to an expansion of vegetated areas (high-Andean steppe), an increase of forest biomes at the expense of herbaceous vegetation (i.e., northward expansion of mesic woodland, cold deciduous forests, and Valdivian rainforest), and the establishment of small pockets of sclerophyllous woodland (Fig. 10a). In all vegetated areas FPC increased with higher [CO₂], most notably between 30 and 40^∘ S (+5–10 % FPC).

Figure 10(a) Effect of atmospheric CO₂, represented as [CO₂], concentrations on the simulated distribution of biomes for the Last Glacial Maximum (LGM) time-slice simulations (left panel: default LGM setup of this study with [CO₂] =180 ppm; right panel: control run with preindustrial [CO₂] =280 ppm). (b) Difference plot of resulting foliar projected cover (280–180 ppm). The number 1 represents Pan de Azúcar, 2 represents Sta. Gracia, 3 represents La Campana, and 4 represents Nahuelbuta.

Download

5.1 Effects of landform simulation mode

As previously mentioned, the novel landform approach illustrated in this study (see Appendix B) was motivated by the fact that standard DVMs do not account for sub-grid heterogeneity controlled by topography and micro-climatic conditions. A common solution for this is to run these models at the same spatial resolution as a LEM, but this is generally computationally wasteful and often violates the scale assumptions of DVMs (i.e., patch size assumption in LPJ-GUESS; see Appendix B). The most striking advantage of landform implementation regarding the realism of simulated FPC is that it allows for the simulation of sub-grid topographic units (i.e., valleys, ridges, slopes) in a computationally efficient way as similar topographic locations are grouped to landforms and only simulated once (or rather n times for the patches of a given landform; see Fig. B1). North and south facing landforms receive a modified amount of solar radiation depending on their average slope and aspect orientation (which affects photosynthesis and transpiration in addition to water availability), and previous studies have showed that this can lead to distinct differences (i.e., Sternberg and Shoshany, 2001). Variable soil depth of landforms is associated with the topographic position (ridges are assumed to feature more shallow soils while valleys have deeper soils due to alluvial material). This results in a differing water storage potential which is especially relevant in dryer locations (Kosmas et al., 2000). Local temperature differences are also considered, as higher landforms are exposed to colder temperatures depending on the elevation difference to the reference elevation of the grid cell (see Appendix B). This leads to a more diverse PFT composition within a grid cell due to the fact that altitudinal zonation of vegetation within a grid cell can be represented (see Figs. 5b and B2b).

However, it has to be noted that we do not currently account for precipitation variations in the landform approach, as a good representation of these effects would require many additional model drivers (i.e., wind speed and direction, sea surface temperature for local fog precipitation; see Gerreaud et al., 2008, 2016; Lehnert et al., 2018) that are either not available at appropriate resolutions or for the time periods considered in this study. Future work may try to incorporate sophisticated statistical downscaling schemes (i.e., Karger et al., 2017) to address this shortcoming.

As can be observed (Fig. B2a), the implemented landform approach has differing effects for the four focus sites. The area-weighted average FPC at the Sta. Gracia site closely resembles the results from the default simulation mode (apart from landform 810 – high altitudes – that only covers 1.7 % of the grid cell area, Fig. B2b). Average-landform and default results for the La Campana site also only differ marginally. However, here a set of landforms of higher altitudes has a substantially lower FPC than the average ($\sim -\mathrm{15}$ %). Variation at the hyper-arid Pan de Azúcar site is lower (as is the FPC), but it is generally higher than the default simulation. The larger vegetation cover simulated in the landform approach also better aligns with MODIS observations for the site, where the default model underestimates satellite-observed cover. The higher FPC in the new model setup is likely a result of the deeper soil profiles of flat and valley landforms that allow for longer water storage versus the default uniform 1.5m soil assumption of LPJ-GUESS. The larger variability of FPC at the Nahuelbuta site can be attributed to the relatively large altitudinal variation in this grid cell (coast to mountainous terrain) and is likely a temperature effect. From these results it can be postulated that the stronger heterogeneity of the simulated FPC using the landform approach will lead to more diverse denudation rates for the topographic units when linked to a LEM and should better resemble observed patterns of vegetation associations in the landscape (we will investigate this in detail in a future study using a coupled DVM–LEM setup).

In order to account for micro-site differences of temperature versus the average climate information per grid cell we stratified the landforms into elevation bins (Sect. 3.1 and Appendix B). Site temperature is an important environmental control for LPJ-GUESS as it impacts vegetation establishment (controlled for individual PFTs by their specific bioclimatic limits), photosynthesis rate, evapotranspiration, and soil decomposition (Smith et al., 2014). We selected 200 m as the bin width after a sensitivity analysis in order to keep the total number of landforms in complex terrain reasonable (potential temperature difference between 100 and 200 m bins < 0.3 ^∘C; data not shown). Furthermore, using more topographic units (separating slopes into lower, mid, and upper slopes) also did not affect the simulation outcome and we therefore opted for a four-unit classification (see Appendix B). These configurations are subjective and should be tested before application under very different terrain or model scales.

5.2 Comparison to vegetation change proxy data

In general, regional vegetation cover of the past is difficult to quantify, as there are limited site-scale pollen, lake level, and midden proxy datasets available in Chile (Marchant et al., 2009). Pollen data from rodent midden (Quebrada del Chaco, 25.5^∘ S, 2670–3550 m, Maldonado et al., 2005) suggest higher winter precipitation at the LGM, higher annual precipitation at 17–14 ka BP, and higher summer precipitation at 14–11 ka BP in northern Chile, but results from another study undertaken at lower elevations indicate absolute desert conditions throughout the Quaternary (Diaz et al., 2012) – which could be caused by regional fog-precipitation. TraCE-21ka precipitation for the Pan de Azúcar site (located at similar latitudes but in coastal lowlands) does not vary from PD conditions throughout the time series and as a result LPJ-GUESS simulates very little vegetation cover (Fig. 9a). Other studies also suggest that LGM conditions were extremely dry, but that the transition to modern climatic conditions in this area occurred as a nonlinear process of multiple moisture pulses over several centuries (Grosjean et al., 2001). Thus, this would imply that a coupled DVM–LEM model would underestimate episodes of high erosion potential due to poor model forcing data.

Pollen reconstructions from Laguna de Tagua Tagua 34.5^∘ S 71.16^∘ W (Heusser, 1990; Valero-Garcés et al., 2005) indicate the presence of extensive temperate woodlands at the LGM that match the simulated mesic woodland biome for this location in our simulations (Fig. 4a). Multiple lake sediment records of the region indicate dry and warm conditions in the MH and the onset of more humid, but strongly seasonal, conditions (winter rainfall) in the late Holocene (Heusser, 1990; Villa-Martínez et al., 2003; Valero-Garcés et al., 2005). While LPJ-GUESS simulates Mediterranean vegetation types (matorral and sclerophyllous woodland) for the MH and a shift to denser xeric woodlands in the late Holocene for this region, substantial variations in the precipitation regime as suggested by the lake records are not present in the TraCE-21ka data used.

Pollen records from the “Chilean Lake District” show a warming trend starting at 17 780 BP (Moreno and Videla, 2016), followed by a trend reversal with major cooling events (14 500 and 12 700 BP). In TraCE-21ka, a sharp cooling event is present at 14 500 BP (likely as an effect of the simulation setup of Meltwater Pulse 1A, Liu et al., 2009), but the second event is not. In our transient simulations this event is reflected in changes in the vegetation composition at three out of the four sites. This climatic change to colder and dryer conditions leads to a reduction of shrub PFTs at the Sta. Gracia (Fig. 9b) and La Campana sites (Fig. 9c) that established with rising temperature and [CO₂] starting at 18 000 BP. As a result, the site biome classification for Sta. Gracia briefly swings back to the initial steppe state and the reduced fuel accumulation leads to fewer fires. However, surface runoff is not affected due to the low precipitation rates. In contrast, the same event at the La Campana site does not result in changes of the fire frequency which is probably due to the presence of sufficient fuel owing to higher annual rainfall rates (> 350 mm). At the Nahuelbuta site this event briefly delays a transition from cold deciduous PFTs to a temperate evergreen/needle-leaved forest (Fig. 9d).

The early and mid-Holocene were the driest periods for central Chile (Heusser, 1990; Valero-Garcés et al., 2005; Maldonado and Villagrán, 2006). This is also visible in the TraCE-21ka dataset (La Campana and Nahuelbuta; Fig. 9c, d). However, the annual precipitation differences in the TraCE-21ka data from this period compared to the LGM and PD are relatively minor; thus, we might not fully capture the true vegetation dynamics indicated by the proxy data of the region (centennial shifts from dry xeric to humid mesic vegetation).

Pollen proxy data from Villagrán (1988) indicate that coastal cool evergreen rainforests were shifted northwards by 5^∘ during the LGM relative to their current position, and that these areas were covered by forest mosaics/parklands instead of the dense forest that is currently present. LPJ-GUESS simulations for the LGM reproduce this latitudinal biome shift (Fig. 5a). While we could not classify an open parkland vegetation type, the simulations show a 5 %–20 % lower FPC for this area at the LGM (Fig. 6b) which can be interpreted as more open conditions. For the MH, temperate Valdivian rainforests similar to PD conditions existed (Villagrán, 1988) and are also simulated by LPJ-GUESS (Fig. 5a)

LGM temperature reductions of up to 12 ^∘C have been proposed for high altitudes in the tropics (Thompson et al., 1995), and LPJ-GUESS also simulates a more expansive belt of high-elevation cold desert for the LGM at the expense of high-Andean steppe (Fig. 5a, b). Furthermore, significantly lower tree lines and vegetation zones of up to 1500 m are reported for the LGM compared to PD (Marchant et al., 2009), which in our simulations is reflected in smaller areas of the cold deciduous forest biome and a shift to lower elevations (Fig. 5b).

In summary, the main forest transition in south-central Chile from the LGM to PD was the postglacial expansion of forests (∼15 ka BP). However, strong climate seasonality and forest clearing during early land utilization (as indicated by increased fire activities from 12 to 6 ka BP) and forest expansion into abandoned land after the Spanish colonial period (Lara et al., 2012) were other major, but anthropogenic, vegetation changes which were followed by intense land clearing for lumber extraction and farming (Armesto et al., 2010). Although significant for many areas, we opted to exclude anthropogenic land use in our analysis as spatial intensities and general utilization are not well understood and are hard to qualify; furthermore, we see this as beyond the scope of this particular study.

Figure 11(a) The effect of paleoclimate input data used for LPJ-GUESS simulations on the spatial and altitudinal biome distribution and (b) the effect on simulated foliage projected cover (given as the difference between TraCE-21ka simulation and ECHAM5 simulation results, Mutz et al., 2018). For a comparison of the differences between average temperatures and precipitation between the two climate datasets see Figs. 3 and S3. The number 1 represents Pan de Azúcar, 2 represents Sta. Gracia, 3 represents La Campana, and 4 represents Nahuelbuta.

Download

5.3 Sensitivity of simulation to paleoclimate input data

As demonstrated in this study, climatic forcing of the ecosystem model is crucial for the simulated vegetation composition, biome establishment, and associated vegetation cover and surface runoff, which are key controls of catchment denudation rates (Jeffery et al., 2014). The TraCE-21ka dataset was chosen since, to our knowledge, it is the only available dataset providing continuous transient monthly data. However, the low spatial resolution (T31/$\sim \mathrm{3.75}{}^{\circ }$) is problematic in resolving regional-scale heterogeneity; therefore, we compare our LGM simulations with a regional paleoclimate model simulation (ECHAM5, resolution: T159/ $\mathrm{0.75}{}^{\circ }\times \mathrm{0.75}{}^{\circ }$, Mutz et al., 2018). Although the general latitudinal patterns of temperature deviations from PD are similar to TraCE-21ka, precipitation gradients and regional anomalies during the LGM are stronger (see Fig. 3 for LGM anomaly plot of TraCE-21ka data and Fig. S1a, b for ECHAM5). While LGM precipitation anomalies in TraCE-21ka do not exceed 650 mm, precipitation rates in the ECHAM5 data vary by well over 1500 mm for the Andean highlands from 27 to 35^∘ S and the coastal areas from 36 to 53^∘ S but are markedly lower for the lower Andes and Tierra del Fuego (Fig. S1b). This leads to a different biome distribution pattern for the LGM (Fig. 11). Due to lower temperatures, the cold deserts extend further north and a small zone of temperate rainforest exists at latitudes from 40 to 42^∘ S, whereas Magellanic woodlands are simulated with TraCE-21ka. A zone classified as Magellanic woodland stretches from 30 to 42^∘ S along the Andean range which gives way to cold deciduous forest at 28 to 32^∘ S at altitude. Temperate deciduous “Maule” forest exists north of 40^∘ S and transitions into sclerophyllous woodland extending to 30^∘ S, whereas TraCE-21ka results lead to steppe at the lowlands and mesic woodlands at higher altitudes. Higher precipitation levels north of 30^∘ S also lead to higher vegetation cover, larger areas of matorral and arid shrubland, and a reduction of desert conditions. This results in substantial differences in simulated foliar projected cover (Fig. 11b). While cover of steppe and mesic woodlands simulated with TraCE-21ka is generally 10 %–20 % higher versus the ECHAM5 simulations, the cover from 35^∘ S to the Atacama Desert is substantially lower than ECHAM5 (higher precipitation rates for this region). Cover south of 42^∘ S is substantially higher for TraCE-21ka as the colder conditions in ECHAM5 simulations generally prohibit vegetation establishment.

The results show that the choice of paleoclimate data clearly has an influence on the simulated vegetation composition, but the impact on vegetation cover depends on the type and location of change. For instance, forest-to-forest transitions due to changes in temperature can have little effect on FPC, while a different annual precipitation rate in Mediterranean to semi-arid conditions can lead to substantial changes in FPC (Sta. Gracia: +200–500 mm precipitation in ECHAM5 LGM data results in > +20 % FPC; Fig. 8a). Given the results from our study we might severely underestimate episodes of elevated erosion in future coupled model exercises – thus, thorough sensitivity studies will be required.

5.4 Impact of atmospheric CO₂ on vegetation cover

The effectiveness of photosynthesis, and the related ability of plants to build up biomass, has a large dependence on [CO₂] (Farquhar et al., 1980; Hickler et al., 2015). Increases in vegetation biomass (expressed in LAI) were observed in our simulations but coincide with simultaneous rises in temperature and [CO₂] (Fig. 9b–d). As observed in the CO₂ sensitivity runs for LGM climate conditions, higher [CO₂] concentrations lead to an expansion of vegetated areas at higher elevations and an expansion of forest biomes at the expense of herbaceous vegetation. Such changes and the magnitude of the CO₂ effect are consistent with earlier LGM simulations (Harrison and Prentice, 2003; Bragg et al., 2013). These changes suggest that [CO₂] could play crucial role in understanding landscape evolution over longer timescales – a factor that again is not considered in previous studies of landscape evolution (i.e., Istanbulluoglu and Bras, 2005; Collins et al., 2004).

The strong limiting role of [CO₂] during the LGM is generally accepted, but the extent to which substantial CO₂ fertilization effects still occur as concentrations increase and rise beyond the current stage is highly debated (Hickler et al., 2015). However, the LPJ-GUESS model with the enabled nitrogen cycle, and nitrogen limiting CO₂ effects, has been shown to generally reproduce experimental observations from “Free Air Carbon Dioxide Enrichment” (FACE) experiments (e.g., Zaehle et al., 2014; Medlyn et at al., 2015; De Kauwe et al., 2017). Enabling the nitrogen cycle in LPJ-GUESS is crucial if one aims to understand vegetation response and landscape evolution at [CO₂] levels above those of the present day, as have occurred in various episodes of Earths' history (i.e., the Miocene).

In general, we were able to show a temporal compatibility of paleovegetation state, vegetation transitions, and simulated vegetation composition using our model (Sect. 5.2). Our results also indicate that the simulated vegetation is strongly influenced by climatic model drivers (see Sect. 5.3) and [CO₂], and by feedback mechanisms occurring within the model (i.e., fire return intervals as controlled by climate and the production of fuel). Given the sensitivity of the vegetation cover to the climate forcing, careful consideration must be made regarding paleoclimate data and its characteristics and uncertainties. The TraCE-21ka data does not show high levels of variability, instead exhibiting only gradual changes over time (temperature, atmospheric CO₂ concentration); furthermore, the TraCE-21ka data does not display strong, centennial to millennial scale trends (precipitation). In contrast, available proxy data suggest stronger variability, at least locally, although this is difficult to generalize for larger areas (see Sect. 5.2). Regional paleoclimate simulations with a higher spatial resolution might be required to obtain better landscape-scale model forcing. Another approach may be to modify the climate model derived forcing based on proxy data, such as increasing inter- and intra-annual variability within reasonable ranges (e.g., Giesecke et al., 2010).