Modeling Potential Distribution and Carbon Dynamics of Natural Terrestrial Ecosystems: A Case Study of Turkey.

We derived a simple model that relates the classification of biogeoclimatezones, (co)existence and fractional coverage of plant functional types (PFTs), and patternsof ecosystem carbon (C) stocks to long-term average values of biogeoclimatic indices in atime- and space-varying fashion from climate-vegetation equilibrium models. ProposedDynamic Ecosystem Classification and Productivity (DECP) model is based on the spatialinterpolation of annual biogeoclimatic variables through multiple linear regression (MLR)models and inverse distance weighting (IDW) and was applied to the entire Turkey of780,595 km² on a 500 m x 500 m grid resolution. Estimated total net primary production(TNPP) values of mutually exclusive PFTs ranged from 108 26 to 891 207 Tg C yr-1under the optimal conditions and from 16 7 to 58 23 Tg C yr-1 under the growth-limiting conditions for all the natural ecosystems in Turkey. Total NPP values ofcoexisting PFTs ranged from 178 36 to 1231 253 Tg C yr-1 under the optimalconditions and from 23 8 to 92 31 Tg C yr-1 under the growth-limiting conditions. Thenational steady state soil organic carbon (SOC) storage in the surface one meter of soil wasestimated to range from 7.5 1.8 to 36.7 7.8 Pg C yr-1 under the optimal conditions andfrom 1.3 0.7 to 5.8 2.6 Pg C yr-1 under the limiting conditions, with the national range of 1.3 to 36.7 Pg C elucidating 0.1% and 2.8% of the global SOC value (1272.4 Pg C), respectively. Our comparisons with literature compilations indicate that estimated patterns of biogeoclimate zones, PFTs, TNPP and SOC storage by the DECP model agree reasonably well with measurements from field and remotely sensed data.


Introduction
Understanding biogeoclimatic controls and its spatio-temporal variability is essential to the quantification of the dynamics of biological productivity under a changing environment at the local, regional and global scales [1][2][3][4][5]. Our current understanding of the seasonal and geographical distribution of the global terrestrial net primary productivity (NPP) estimated at 56.4 Pg carbon (C) yr -1 (on average, 426 g C m -2 yr -1 ) [6] and 59 Pg C yr -1 is based on the extrapolation of local and regional studies to the global scale (1 Pg = 10 15 g) [7,8]. Dynamic classification of regional plant functional types (PFTs) in response to changes in forcing biogeoclimate variables such as elevation, geographical position, moisture index, biotemperature, and growing season precipitation is needed for a better estimation of the global NPP, sustainable management of natural resources, and modelling of biogeochemical cycles [9][10][11]55]. Changes in the predominant PFTs are primarily determined using the analysis of time series datasets derived from one or the combination of the following sources: (1) atlases [12], (2) remote sensing [13], and (3) biogeoclimate relationships [14][15][16]. The use of atlases and satellite images serves to describe the actual distribution of natural and cultivated vegetation patterns, while (geo)statistical and process-based models based on biogeoclimatic controls serve not only to describe but also to predict the potential distribution of natural PFTs and changes in patterns of plant and soil C storage in a changing global climate.
Flux rates, storage sizes, and residence times of C cycle are intimately coupled with the distribution patterns of biogeoclimate zones. A various number of spatial interpolation techniques, multiple regression models, biogeoclimate indices, and land-cover classification systems have been implemented to quantify dynamics of terrestrial biogeochemical metabolism including NPP and net ecosystem productivity (NEP) [14,15,[17][18][19]. Due to the impracticality of continuous field NPP measurements of all ecosystem types at the regional and global scales by harvest or 14 C-based methods, various algorithms have been devised to spatially interpolate biogeoclimatic datasets for each pixel of a gridded digital elevation model (DEM) [20] and to use them as inputs into processed-based ecosystem models [21][22][23][24][25][26][27].
In this study, we present a simple algorithm of Dynamic Ecosystem Classification and Productivity (called DECP hereafter) to quantify the dynamics of potential natural PFTs, NPP, and soil organic carbon (SOC) as a function of biogeoclimatic determinants that reflect the geo-referenced long-term mean bioclimate and apply it to the entire Turkey of 780,595 km 2 .

Description of Study Region
Turkey (36-42°N and 26-45°E) is located where Asia, Europe, and the Middle East meet, with an average altitude of 1250 m. The temperature reaches 45 °C in July in the southeastern region, and falls to -30 °C in February in the eastern regions, with a mean annual temperature of around 13 °C. Annual precipitation ranges from 258 mm in the central and southeastern regions to 2220 mm in the northeastern Black Sea coasts, with a mean annual precipitation of around 634 mm. Annual evapotranspiration varies from 624 mm in the eastern region to 2400 mm in the southeastern region, with a mean annual evapotranspiration of 1280 mm according to the long-term mean climate data between 1968 and 2004 [32].

Derivation of Bioclimatic Indices
Monthly climate data in Turkey were obtained from 269 meteorological stations for the period of 1968 to 2004 and included solar radiation (SR, MJ m -2 ), mean, minimum and maximum air temperature (T mean , T min and T max , o C), cloudiness (CLD, %), potential evapotranspiration (PET, mm), precipitation (PPT, mm), and soil temperature (0 to 10 cm in depth) (ST 10 , o C) [32]. The three Holdridge life zone (HLZ) classification indices of biotemperature (BT), potential evapotranspiration ratio (PER), and potential evapotranspiration (PET HLZ ) [15], and the following three bioclimatic indices used by Box [18] were derived from the climate data for each data point as follows: where mean annual biotemperature (BT, o C) was calculated by substituting monthly mean temperatures (MMT i ) both above 30 °C and below 0 °C with 0 °C [15].
where MMT coldest is mean monthly temperature of the three coldest months (MMT min , o C).
where GSP is growing season precipitation calculated as the mean monthly precipitation (mm) of the three warmest months (MMP warmesti) .

MI = HLZ
PPT PET (4) where annual moisture index (MI) refers to the ratio of the mean annual precipitation (PPT, mm yr -1 ) to mean annual potential evapotranspiration (PET HLZ , mm yr -1 ).

Mapping Biogeoclimate Zones and Potential Natural Plant Functional Types
A national map of potential natural vegetation was derived from (biogeo)climatic variables and indices, assuming a natural vegetation distribution is in equilibrium with the mean long-term climate, without human interference and cultivated vegetation types. Natural vegetation is thus separated into three broad PFTs: trees, shrubs, and grass. Three life-forms (evergreen needleleaf, EN; deciduous broadleaf, DB; and evergreen broadleaf, EB) were distinguished for tree and shrub covers in terms of leaf phenology (evergreen vs. deciduous) and leaf shape (needleleaf vs. broadleaf), and one grass lifeform (C3) in terms of photosynthetic pathway since the life-forms of deciduous needleleaf, and C4 grass species are not found in the prevailing natural vegetation of Turkey.
The fractional cover of trees (f T ) was estimated from the annual MI and assumed to linearly decrease from 100% at MI = 1.0 to 0% at MI = 0.6 as follows [18]: for MI < 0.6 (MI-0.6)/0.4 for 0.6 MI 1.0 Treeline at altitudes beyond which tree cover fraction becomes zero was assumed to occur at a growing degree-day sum (GDD 5 ) ≤ 350 according to Woodward [33]. GDD 5 was calculated thus: where GDD 5 refers to the annual sum of monthly mean temperatures above the base temperature (T b = 5 o C) multiplied by the number of days in the month i (m i ).
According to the 17-class land cover classification system of the International Geosphere-Biosphere Programme (IGBP) and the IGBP DISCover definition of forest (greater than 60% canopy cover), four tree cover classes of <l0%, 10-30%, 30-60% and 60-100% were assumed to translate to barren or sparsely vegetated, grassland (steppe), woodland/shrubland (W/S) and forest cover types, respectively [34]. Unlike the mutually exclusive four cover categories derived from the IGBP system, fractional shrub and grass covers compatible with the fractional forest cover were also mapped directly based on the assumption that the 500 m x 500 m pixels do not all consist of a single PFT. The cover fractions of shrubs and C3 grass were estimated based on MLR models developed by Paruelo and Lauenroth [35] as follows: f S = 1.7105 + 1.5451PPT DJF -0.2918 ln PPT (n = 70; R 2 adj. = 0.62; P < 0.001) (7) f C3 = 1.1905 -0.02909MAT + 0.1781 ln PPT DJF -0.2383BIOME (n = 69; R 2 adj. = 0.37; P < 0.001) (8) where f S and f C3 refer to the cover fractions of shrubs and C3 grass, respectively, and were linearly normalized based on 1-f T . PPT DJF is the ratio of PPT during the three winter months (December to February) to annual PPT. The value of BIOME is an indicator (dummy) variable representing grassland (1) vs. shrubland (2) and 1 at f S ≤ 0.2 and 2 at f S > 0.2. Assignment of f T to one of PFTs (evergreen needleleaf, EN; deciduous broadleaf, DB; evergreen broadleaf, EB) was based on MMT coldest and GSP according to Box [18]. The mixture (100 to 0%) of the PFTs (f PFTs ) was linearly interpolated according to Box [18] as follows: The (bio)climatic variables of BT, MMT coldest , and GSP for the classification of biogeoclimate zones; MAT, PPT, PET HLZ , PPT DJF , GDD 5 , and soil temperature at the depth of 0 to 10 cm during the three warmest months (June to August) of growing season (ST JJA10 ) for the cover fraction of PFTs; and MMT coldest for the relative dominance of the life-form composition were mapped using geostatistical interpolation based on inverse distance weighted (IDW) method and multiple linear regression (MLR) models in ArcGIS 9.1 [36]. The spatial interpolation of the long-term mean annual bioclimatic variables was carried out for each (about 500 m x 500 m) of 3,182,206 grid cells (ca. 780,595 km 2 ). IDW was used to create accurate surface maps of PPT and GSP, while the remaining bioclimatic variables were mapped based on their best MLR models using digital layers for each explanatory variable of DEM (m), latitude (LAT, m), longitude (LON, m), distance to sea (DtS, m), and aspect (Asp, compass degree) ( Figure 1). Best MLR models were selected based on lowest Cp and highest R 2 adj. , and significant P values (<0.05) for all the explanatory variables. Digital elevation model was constructed from a 1:250,000 scale topographic map of Turkey generated by the Turkish General Command of Mapping, projected to the Universal Transverse Mercator (UTM) coordinates of the World Geographic System (WGS 1984) and re-sampled to a grid size of 500 m. Digital maps of distance to sea (m) and aspect ( o ) were derived from DEM using ArcGIS 9.1. The classification and distribution patterns of major biogeoclimate zones and PFTs were revealed overlaying the maps of the biogeoclimate variables, and detailed classification algorithm is presented in Figure 1.

Quantifying Total Net Primary Productivity and Soil Organic Carbon Density
The annual quantification of TNPP was based on the approach of light use efficiency (LUE) by Monteith [37,38], and Kumar and Monteith [39] as follows: where TNPP is total net primary productivity (g C m -2 yr -1 ). LUE T (g DM MJ -1 ) is light use efficiency of intercepted photosynthetically active radiation (PAR) into total dry matter (DM) of aboveground and belowground biomass. The value of 0.45 is a conversion coefficient for C content per unit DM biomass [40]. Biome-specific minimum, mean and maximum LUE T values were obtained from Ruimy et al. [40] ( Table 1). The fraction of incident PAR in incoming solar radiation was assumed to be 0.48 according to McCree [41]. The amount of intercepted PAR (IPAR) was assumed to be a function of fractional canopy covers of PFTs (fPFTs). The length of growing season period (GSL) (in days) was estimated as the number of days between the first and last months when the monthly minimum temperature (T min ) < 0 o C. LCF refers to growth-limiting climate factors and is a reduction factor consisting of cloudiness (CLD) and humidity (MI) indices combined and normalized to a scale of 0 to 1 (LCF is 0 for completely growth-limiting conditions and 1 for optimal conditions). Aboveground NPP (NPP A ) and belowground NPP (NPP B ) values were estimated based on the ratios of NPP B to NPP A compiled by the literature review of Ruimy et al. [40] (Table 1). Net ecosystem productivity (NEP) can be quantified as follows: where NEP represents the net ecosystem exchange of CO 2 between the atmosphere and an ecosystem after accounting for the rate of C release from the soil to the atmosphere by respiration of microorganisms and roots (R h ). NEP reveals whether a terrestrial ecosystem is acting as a sink (net ecosystem sequestration) or a source (net ecosystem emission) for atmospheric CO 2 .
Based on the empirical Miami model modified by Friedlingstein et al. [42], originally developed by Lieth [43], TNPP values were also estimated as a function of PPT as follows: in nontropical regions (12) where TNPP (PPT) represents dependence of annual net primary production (kg C m -2 yr -1 ) on precipitation (PPT) (mm yr -1 ). Dai and Fung [44] reported that the modified Miami model estimated TNPP ranges of major terrestrial ecosystems closer to observations than the original Miami model. The dynamics of SOC pools can be expressed in a single compartment by the following simple form of the first-order ordinary differential equation [45]: where SOC is the quantity of soil organic carbon aggregated as one compartment, A the total mean annual litter (organic carbon) input to the soil, and k the annual decomposition rate of SOC. At steady state, annual NPP should equal annual R h which can be expressed as follows: Based on natural vegetation and undisturbed sites of the global observed SOC density dataset (100 cm in depth) [46], Yang et al. [47] reported a simple model for steady state k values as follows: where MAT is mean annual temperature, k w and k d decomposition rates for wet and dry soil conditions, respectively. PER provides an aridity index that represents the interactive impact of BT and PPT on decomposition rates [15,48,54]. Soil organic carbon density (kg C m -3 ) at steady state was directly derived from Eqn (14), based on the DECP model, the modified Miami model, and steady state k values. As with the classification of biogeoclimate zones and PFTs, all the variables and parameters were mapped to a regular 500 m x 500 m pixel for the quantification of TNPP and SOC, based on their best MLRs and IDW interpolations, using ArcGIS 9.1.

Biogeoclimate Zones
Classification of biogeoclimate zones was based on the overlay of annual surface maps created by geographical position (latitude, longitude, elevation, and distance to sea)-sensitive MLR models of MMT coldest and BT, and by the IDW interpolation of GSP. Best MLR models elucidated 89.8% of variation in BT as a function of LAT, LON, and DEM (R 2 adj. = 89.8%; P < 0.001; n = 265) and 91.2% of variation in MMT coldest as a function of LON, DEM, and DtS (R 2 adj. = 91.2%; P < 0.001; n = 265) ( Figure 1). The power that determines how influential neighboring points are to the point whose value is being interpolated (the higher the power, the lesser the influence from distant points) was used as the default optimized power value of 1.6274 in the IDW interpolation of GSP (n = 269). Cross-validation between predicted and observed values of GSP as a consequence of the IDW interpolation resulted in R 2 of 78.9%, mean prediction error of -0.0721 mm, and root mean square prediction error of 9.58 mm. The classification algorithm resulted in a total of 14 biogeoclimate zones (P < 0.001) ( Figure 2). The resulting biogeoclimate zones were distinguished due to what appeared to be completely distinct or disjunct bioclimatic and geographical ranges. About 82% and 18% of Turkey appeared to comprise dry and moist ecosystems, respectively. (Sub)nival, alpine, boreal, cool temperate, warm temperate and Mediterranean ecosystems covered 0.19%, 0.64%, 5.54%, 51.72%, 33.89, and 8.03% of Turkey, respectively ( Table 2).

Potential Natural Land Cover and Plant Functional Types
The existence and fraction of tree cover drive potential distribution of natural land cover types and were derived from the MLR-based surface maps of GDD 5 , and MI (the ratio of PPT to PET HLZ ) (  According to the IGBP land-cover classification system, the rule-based generation of potential natural land cover map showed that forest (fT ≥ 60%), W/S (30% ≤ fT < 60%), grassland (steppe) (10% ≤ fT < 30%) and barren or sparsely vegetated (fT < 10%) ecosystems covered 371,097 km 2 (47.5%), 150,608 km 2 (19.3%), 96,533 km 2 (12.4%), and 162,357 km 2 (20.8%) of Turkey, respectively (Table 3). Once the digital map layers of the biogeoclimate zones and the IGBP land-cover classification were overlaid, biome-specific land cover types were determined ( Figure 4). Existence, fractional cover, and spatial extent of PFT mixtures (EN versus DB forests and W/S) were determined using the gradient maps of MMT coldest and GSP ( Figure 5) (Table 3). According to the relative dominance of EN versus DB trees, 63%, 26%, and 11% of the forests appeared to be mixed (EN or DB < 70%), EN and DB forests, respectively. Pure EN forests (EN or DB ≥ 70%) range from boreal EN forests (e.g., Abies nordmanniana, and Picea orientalis) to Mediterranean EN forests (e.g., Pinus brutia, and Pinus pinea) along the MMT coldest and GSP gradients across Turkey.
The MLR models by Paruelo and Lauenroth [35] of fS and fC3 as a function of PPT and MAT that were normalized to fT were used to estimate the coexistence and fractional cover of trees, shrubs, and C3 grass in a single pixel based on the IGBP land-cover classification (Figure 1). MLR-based spatial interpolations of PPT DJF and MAT had R 2 adj. values of 73.1% as a function of LAT, LON, DEM, and DtS and 90.1% as a function of LAT, LON, DEM, DtS, and aspect, respectively (P < 0.001; n = 265). Areas with fS of <20%, 20% to 50%, 50% to 80%, and >80% comprise 302,331 km 2 (39%) of the eastern and northern Anatolia and southwestern Mediterranean regions; 120,030 km 2 (15%) of the northwestern and central Anatolia regions; and 154,259 km 2 (20%) and 203,975 km 2 (26%) of the central, southeastern and western Anatolia and southeastern Mediterranean regions, respectively ( Figure 6). The distribution map of C3 grass in Turkey revealed that lowland steppes are concentrated in the central and southeastern Anatolia regions, and the Iğdır plain in the eastern border of Turkey, while highland steppes mostly take place in the northeastern and eastern Anatolia regions (Figure 7).

Total Net Primary Production under Growth-Limiting and Optimal Conditions
For each land cover type of the IGBP classification system, minimum, mean and maximum TNPP values of both mutually exclusive and inclusive PFTs under the growth-limiting and optimal environmental conditions were predicted for each pixel of 500 m x 500 m resolution based on the DECP model ( Table 4) The DECP model revealed that TNPP values of mutually exclusive PFTs ranged from 108 + 26 to 891 + 207 Tg C yr -1 under the optimal conditions and from 16 + 7 to 58 + 23 Tg C yr -1 under the growth-limiting conditions for all the natural ecosystems in Turkey (Table 4). Total NPP values of coexisting PFTs ranged from 178 + 36 to 1231 + 253 Tg C yr -1 under the optimal conditions and from 23 + 8 to 92 + 31 Tg C yr -1 under the growth-limiting conditions. Forest TNPP varied between 147 + 41 g C m -2 yr -1 in boreal EN forest and 3431 + 709 g C m -2 yr -1 in warm temperate DB forest under the optimum conditions and between 26 + 11 g C m -2 yr -1 in warm temperate DB forest and 240 + 62 g C m -2 yr -1 in Mediterranean EN forest under the growth-limiting conditions according to the mutually exclusive IGBP classification (Table 4). Minimum and maximum TNPP values of W/S ecosystems were estimated at 87 + 19 g C m -2 yr -1 and 1689 + 371 g C m -2 yr -1 for warm temperate DB W/S under the optimal conditions and at 9 + 3 g C m -2 yr -1 in warm temperate DB W/S and 116 + 18 g C m -2 yr -1 in Mediterranean EN W/S under the growth-limiting conditions, respectively. Steppe ecosystems in Turkey were estimated to have TNPP values ranging from 57 + 17 g C m -2 yr -1 (5 + 2 g C m -2 yr -1 under the limiting conditions) in the cool temperate zone to 1023 + 319 g C m -2 yr -1 (51 + 18 g C m -2 yr -1 under the limiting conditions) in the Mediterranean zone under the optimal conditions (Table 4).
Our DECP model did not estimate any productive natural area for the mutually exclusive IGBP land cover of barren or sparsely vegetated land. However, the coexistence of PFT mixtures revealed TNPP values by the shrub and grass life forms of 39 + 5 Tg C yr -1 to 211 + 26 Tg C yr -1 under the optimal conditions and of 3 + 1 Tg C yr -1 to 17 + 4 Tg C yr -1 under the limiting conditions for the barren or sparsely vegetated ecosystems. Shrub TNPP values under the coexistence of PFTs ranged from 26 + 8 g C m -2 yr -1 in the (sub)nival zone to 1475 + 202 g C m -2 yr -1 in the warm temperate steppe under the optimal conditions and from 8 + 2 g C m -2 yr -1 in the warm temperate forest to 153 + 29 g C m -2 yr -1 in the warm temperate steppe under the limiting conditions. The algorithm of the DECP model for the coexisting PFTs led to a 2.5-to 4.6-fold increase in the estimation of steppe TNPP values under the (non)-limiting conditions, relative to the algorithm for the mutually exclusive PFTs.   TNPP: total net primary productivity, B: boreal, CT: cool temperate, WT: warm temperate, M: Mediterranean, EN: evergreen needleleaf, DB: deciduous broadleaf, W/S: woodland/shrubland, BSVL: barren or sparsely vegetated land. Rows designated with and without the term "IGBP" indicate values estimated according to mutually exclusive and inclusive algorithms using the IGBP land-cover classification and MLR models by Paruelo and Lauenroth [35], respectively. Values given as mean + standard deviation may not total due to rounding.
The DECP model estimated the areal extent of productive forest ecosystems as 349,084 km 2 (44.7% of the entire country). Total spatial extent of productive steppe ecosystems was estimated as about 12% of the total land area based on the IGBP classification of mutually exclusive and inclusive PFTs with the fractional cover of 10% to 30% and as 94% of the total land area based on the MLR model [35] with the fractional cover less than 10%. Productive W/S ecosystems occupied 18% to 19% of the total land area with the fractional cover rule (30% to 60%) of the IGBP classification and 60% the total land area with the fractional cover range of 0% to 100% (Table 4).
Based on the IDW interpolation of PPT, the modified Miami model [42] was used to estimate mean TNPP values of biogeoclimatic land cover types under the optimal conditions (Table 5). Mean TNPP values estimated by the modified Miami model across Turkey ranged from 427 + 100 g C m -2 yr -1 in boreal forest to 653 + 49 g C m -2 yr -1 in Mediterranean forest; from 359 + 30 g C m -2 yr -1 in cool temperate W/S to 564 + 30 g C m -2 yr -1 in Mediterranean W/S; and from 335 + 20 g C m -2 yr -1 in cool temperate steppe to 499 + 28 g C m -2 yr -1 in Mediterranean steppe (Table 5). Sub(nival) vegetation type was estimated to have the second highest TNPP of 584 + 200 g C m -2 yr -1 after that of Mediterranean forest.

Soil Organic Carbon Density
Steady state SOC density for a depth of 1 m based on the potential distribution patterns of TNPP and biogeoclimate zones, and the gradients of BT, MAT, and PPT was lowest in warm temperate forest (3.5 + 0.7 kg C m -3 yr -1 ) and W/S (1.9 + 0.4 kg C m -3 yr -1 ), and cool temperate steppe (1.3 + 0.4 kg C m -3 yr -1 ) and highest in warm temperate forest (31 + 7 kg C m -3 yr -1 ) and W/S (17 + 4 kg C m -3 yr -1 ), and Mediterranean steppe (11 + 3 kg C m -3 yr -1 ) under the optimal conditions. Under the limiting conditions, SOC storage ranges from 0.5 + 0.2 kg C m -3 yr -1 in warm temperate forest to 4.7 + 0.9 kg C m -3 yr -1 in Mediterranean forest; from 0.2 + 0.07 kg C m -3 yr -1 in warm temperate W/S to 2.2 + 0.6 kg C m -3 yr -1 in Mediterranean W/S; and from 0.1 + 0.04 kg C m -3 yr -1 in cool temperate steppe to 1.2 + 0.4 kg C m -3 yr -1 in Mediterranean steppe (Table 6). Values given as mean + standard deviation may not total due to rounding.
Estimates of the modified Miami model for steady state amount of SOC storage range from 9.5 + 3.8, 8.1 + 0.6 and 7.9 + 0.4 kg C m -3 yr -1 in cool temperate forest, W/S and steppe to 13.1 + 1, 12.8 + 0.6 and 11.9 + 0.6 kg C m -3 yr -1 in Mediterranean forest, W/S and steppe, respectively ( Table 5). The national steady state SOC pool in the surface one meter of soil was estimated to range from 7.5 + 1.8 to 36.7 + 7.8 Pg C yr -1 under the optimal conditions and from 1.3 + 0.7 to 5.8 + 2.6 Pg C yr -1 under the limiting conditions, based on the DECP model. The modified Miami model estimated, on average, potential SOC storage of 40.9 + 14 Pg C yr -1 including the amount of SOC stored in barren or sparsely vegetated land cover.

Discussion
The proposed DECP algorithm links patterns of land cover and potential natural vegetation dynamics to changes in NPP and steady state SOC density on the national scale. Biogeoclimatic classification of the DECP model appears to represent distribution, spatial extent, and mosaics of potential natural PFTs of Turkey in harmony with evident differences in long-term mean climate data rendered sensitive to the geographic gradients of elevation, latitude, longitude, aspect, and distance to sea through MLR-based spatial interpolations. Classification accuracy of the DECP model is a function of the classification categories as well as quality of the source data used to derive biogeoclimate indices and MLR models. Spatial extent of potential natural forest ecosystems of 349,084 km 2 estimated by the DECP model is about 64% greater than the actual forest area of 212,000 km 2 reported for Turkey [53]. Values given as mean + standard deviation may not total due to rounding.
Literature compilations by Ruimy et al. [40] of biome-specific TNPP values indicated ranges of 100 to 2000 g C m -2 yr -1 for Mediterranean forest, 100 to 900 g C m -2 yr -1 for temperate forest, 100 to 800 g C m -2 yr -1 for boreal forest, 20 to 1500 g C m -2 yr -1 for W/S of all latitudes, 50 to 2000 g C m -2 yr -1 for temperate grassland, and 20 to 300 g C m -2 yr -1 for tundra grassland. The range of TNPP estimates, particularly under the optimal conditions, by the DECP model and the modified Miami model agrees well with the minimum and maximum values of TNPP estimated in the literature compilations of TNPP by Ruimy et al. [40]. Our comparisons for different land cover types show that the DECP model estimates of TNPP fall within the ranges of the MODIS-derived computations of NPP by Running et al. [49] and existing field data by Kucharik et al. [50] and Zheng et al. [11]. On the national scale, minimum TNPP (16 Tg C yr -1 ) of mutually exclusive PFTs under the limiting conditions and maximum TNPP (1231 Tg C yr -1 ) of mutually inclusive PFTs under the optimal conditions according to the DECP model, and mean TNPP (1.5 Pg C yr -1 ) of mutually exclusive PFTs under the optimal conditions according to the modified Miami model elucidate 0.02%, 2.0%, and 2.5% of the global mean TNPP value of 59 Pg C yr -1 .
The DECP model predicted SOC storage patterns across the biogeoclimate zones and PFTs of Turkey similar to those revealed by global SOC dataset derived from potential natural vegetationsupporting soil profiles and classified according to HLZ [48]. Global SOC data by Post et al. [48] showed 9.2 + 4.5 kg C m -3 for Mediterranean moist forest, 11.5 + 13.9 kg C m -3 for Mediterranean dry forest, 5.4 + 2.2 kg C m -3 for Mediterranean W/S, 9.3 + 7.4 kg C m -3 for warm temperate moist forest, 8.3 + 4.6 kg C m -3 for warm temperate dry forest, 7.6 + 6.8 kg C m -3 for warm temperate steppe, 12 + 8.2 kg C m -3 for cool temperate moist forest, 13.3 + 9.5 kg C m -3 for cool temperate steppe, and 15.5 + 30 kg C m -3 for boreal moist forest. The estimates of SOC values under the optimal conditions by the DECP and modified Miami models are well within one standard deviation of the global SOC dataset by Post et al. [48]. Our national range of estimates for SOC storage of 1.3 Pg C to 36.7 Pg C by the DECP model and 40.9 Pg C by the modified Miami model accounts for 0.1%, 2.8%, and 3.2% of the global SOC value (1272.4 Pg C) according to Post et al. [48], respectively.
According to DECP, mean SOC carbon density appeared to decrease for forest and W/S ecosystems and increase for steppe ecosystems in transition from cool temperate to Mediterranean zones. This pattern for forest and W/S reveals that SOC storage decreases in response to an increase in BT and PER (the ratio of PET HLZ to PPT), and a decrease in elevation and GSP. When PER (1.1 + 0.35) approached unity in the cool temperate zone, SOC pool under the cool temperate forest peaked at 12.4 + 3.3 kg C m -3 which shows a close agreement with mean SOC pool value of 10 kg C m -3 at PER of 1 derived from the global soil dataset [48]. The reverse pattern for steppe ecosystems appears to be due to the increase in TNPP values estimated in transition from cool temperate to Mediterranean zones. On the other hand, the modified Miami model estimated increases in TNPP and SOC for forest, W/S and steppe ecosystems in transition from cool temperate to Mediterranean zones and failed to capture detailed patterns among biogeoclimate variables, SOC storage and TNPP revealed by the DECP algorithm.
Values of vegetation-specific parameters used in the DECP model approximately correspond to those of the biogeochemical models TEM [23], CASA [51] and SILVAN [52]. As with all regression models, MLR models should be used cautiously to estimate response variables at large spatial scales beyond the region from which they are derived. The proposed DECP model does not account for certain areas of a grid cell that may be unsuitable for establishment of natural vegetation (including rivers, lakes, roads, and rock and stony outcrops), and for changes caused by human disturbances (deforestation, urban sprawl, and land-use changes). The model must be further validated by other means available such as satellite images before its use in predicting spatio-temporal responses of biogeoclimate zones, PFTs, TNPP, and SOC to the projected scenarios of global climate change [56]. According to available field and remotely sensed data, the DECP model based on the 37-year mean climatic data was able to generate realistic distributions of potential natural land cover and estimates of TNPP and steady state SOC in Turkey under the limiting and non-limiting conditions in mutually exclusive and inclusive ways of PFTs. The difference between the steady state estimates of net C stocks for potential natural ecosystems under the non-limiting and limiting conditions may reveal the magnitude of historic national C loss as well as the national potential for C sequestration if 50% recovery of C loss is assumed over the next 50-100 years as a reasonable upper limit in constant climate.
This research provides the first extensive quantification of potential conditions for NPP and SOC in Turkey. Coupling this information with current land-use/land-cover and related NPP would lead to important implications about C loss due to human-induced disturbances across Turkey since most of the C loss is most likely to result from alteration of land-use/land-cover, and land management practices. Our model potentially estimated ca. 45% of the land area of Turkey to be productive forest ecosystems, while only ca. 25% of the country is currently forested, with more than half being unproductive. This reveals the significant potential for Turkey's forests to sequester C as well as the importance of sustainable land management practices to achieve C-sequestration potential given rapid rates of population growth and conversions of forests and grasslands into urban-industrial and cropland areas.