Estimating spatial variation in the effects of climate change on the net primary production of Japanese cedar plantations based on modeled carbon dynamics

Abstract

[...]the PTDB enables model users to select parameters in a physiologically supported manner. Japan has the unique characteristic that a single tree species is planted across a wide climatic range, from subtropical to cool temperate [18]. [...]the parameterization and validation of process-based models for planted forests in this region remain a challenge, and clarification of the application range of such models is needed. Forest covers two-thirds of the land area at present, and a stand age of approximately 50–60 years is common for planted forests [20]; those stands are now mature and can be harvested. [...]forest resources in Japan are currently rich. [...]the effects of ongoing climate change on mature tree plantations remain unclear, and previous research has suggested the potential risk of declining stand productivity [5, 21, 22]. [...]demand for spatiotemporal prediction of the productivity of planted forests in Japan is increasing to support the implementation of adaptation strategies for climate change.

About the Authors:

Jumpei Toriyama

Roles Conceptualization, Formal analysis, Methodology, Software, Writing – original draft, Writing – review & editing

* E-mail: [email protected]

Affiliation: Kyushu Research Center, Forestry and Forest Products Research Institute (FFPRI), Kumamoto-city, Kumamoto, Japan

ORCID logo https://orcid.org/0000-0001-8061-6398

Shoji Hashimoto

Roles Conceptualization, Methodology, Project administration, Supervision, Writing – review & editing

Affiliations Department of Forest Soils, FFPRI, Tsukuba, Ibaraki, Japan, Graduate School of Agricultural and Life Sciences, The University of Tokyo, Bunkyo-ku, Tokyo, Japan

ORCID logo https://orcid.org/0000-0003-3022-7495

Yoko Osone

Roles Data curation, Formal analysis, Resources, Validation, Writing – review & editing

Affiliation: Department of Forest Soils, FFPRI, Tsukuba, Ibaraki, Japan

Naoyuki Yamashita

Roles Data curation, Resources, Validation, Writing – review & editing

Affiliation: Department of Forest Soils, FFPRI, Tsukuba, Ibaraki, Japan

ORCID logo https://orcid.org/0000-0003-3398-6825

Tatsuya Tsurita

Roles Data curation, Resources, Writing – review & editing

Affiliation: Department of Forest Soils, FFPRI, Tsukuba, Ibaraki, Japan

Takanori Shimizu

Roles Data curation, Resources, Writing – review & editing

Affiliation: Department of Disaster Prevention, Meteorology and Hydrology, FFPRI, Tsukuba, Ibaraki, Japan

Taku M. Saitoh

Roles Data curation, Methodology, Resources, Writing – review & editing

Affiliation: River Basin Research Center, Gifu University, Gifu-city, Gifu, Japan

Shinji Sawano

Roles Resources, Writing – review & editing

Affiliation: Hokkaido Research Center, FFPRI, Sapporo, Hokkaido, Japan

Aleksi Lehtonen

Roles Supervision, Writing – review & editing

Affiliation: Natural Resources Institute Finland, Helsinki, Finland

ORCID logo https://orcid.org/0000-0003-1388-0388

Shigehiro Ishizuka

Roles Resources, Writing – review & editing

Affiliation: Department of Forest Soils, FFPRI, Tsukuba, Ibaraki, Japan

Introduction

At present, spatiotemporal prediction of the response of planted forests to a changing climate is increasing in importance. The planted forest area worldwide reached 293.9 million ha, which is 7% of total forest area, by 2020 [1]. This increase has highlighted the importance of predicting the spatial responses of planted forests to changing climatic conditions in terms of material and energy flows in terrestrial ecosystems. In particular, predicting the future risk of decreasing productivity and carbon (C) sequestration in planted forests and modeling of C dynamics [2] are expected play important roles.

Prediction of a planted forest’s response to climate change through modeling of C dynamics depends strongly on the climate scenarios employed. To fully use the information in climate scenarios and quantify changes in the stock and flow of C, process-based modeling, which quantifies processes such as photosynthesis and respiration and represents them as numerical expressions, is an effective approach. Several process-based models, including ORCHIDEE [3], LPJ-GUESS [4] and 3-PG [5], have been applied to planted forests, and Biome-BGC (BBGC) has been used to model various tree species in Europe [6, 7] and East Asia [8, 9].

Appropriate parameterization is an essential step in the use of process-based models. Process-based models have been parameterized for various plant functional types [10]. In recent years, however, a powerful dataset of plant traits has been developed for each tree species [11, 12], enabling species-specific or plant trait-based parameterization of process-based models [13] with a focus on the relationship between the diversity of plant traits and the ecosystem response to climate change [14]. A plant trait database (PTDB) for Japanese cedar (Cryptomeria japonica D. Don), which is described later, and Japanese cypress (Chamaecyparis obtusa) was published recently [15]. We can refer to this PTDB for representative values of eco-physiological parameters to apply in process-based models, such as leaf nitrogen (N) content and specific leaf area (SLA). Thus, the PTDB enables model users to select parameters in a physiologically supported manner. In addition, Bayesian calibration [16, 17] or optimization (BO) is an effective approach for certain parameters with little or no data in the PTDB. These Bayesian approaches support the selection of parameters through repeated model runs, known as Markov chain Monte Carlo (MCMC) sampling, in which data are effectively used to estimate optimal parameter values and their uncertainties [16, 17].

The modeling research described above has also focused on the coastal region of East Asia, which is the region with the highest proportion of planted forest [1] and is characterized by a relatively warm and humid climate. Japan has the unique characteristic that a single tree species is planted across a wide climatic range, from subtropical to cool temperate [18]. Thus, the parameterization and validation of process-based models for planted forests in this region remain a challenge, and clarification of the application range of such models is needed. On the other hand, Japan increased its stock of planted trees at an extremely high rate after World War II [19]. Forest covers two-thirds of the land area at present, and a stand age of approximately 50–60 years is common for planted forests [20]; those stands are now mature and can be harvested. Thus, forest resources in Japan are currently rich. However, the effects of ongoing climate change on mature tree plantations remain unclear, and previous research has suggested the potential risk of declining stand productivity [5, 21, 22]. Therefore, demand for spatiotemporal prediction of the productivity of planted forests in Japan is increasing to support the implementation of adaptation strategies for climate change. Research into the productivity of Japanese cedar plantations, which cover 18% of the total forested area and 44% of the planted forest area in Japan [20], must be undertaken with the highest priority due to its importance to the timber industry.

Previous studies have shown that Japanese cedar may be vulnerable to climate change, especially drier climatic conditions. Matsumoto et al. [22] estimated the vulnerability of Japanese cedar to a warming climate at the national scale based on two indices—namely, the ratio of annual transpiration to precipitation and soil water-holding capacity. They found that the area of high vulnerability is distributed in the western part of Japan and that Japanese cedar growth may decrease extensively from western to eastern Japan under future climate scenarios with increases in mean annual temperature of 2.2–3.2 degrees C by 2081–2100 [22]. Matsumoto et al. [22] calculated transpiration [23] based on stomatal conductance [24–26] and vapor pressure deficit (VPD) using second mesh order (approximately 10-km resolution) climate scenarios at a monthly timescale. At present, third mesh order (approximately 1-km resolution) climate scenarios at a daily timescale are available [27], and prediction of the stomatal response of Japanese cedar to climate scenarios with high spatiotemporal resolution can be undertaken. In addition, quantitative prediction of the changes in stand productivity is essential for the development of forestry policy.

Based on this background, we focus on estimating changes in stand productivity, evaluated as net primary production (NPP), rather than on mortality under severe climatic conditions. In extremely severe climate scenarios, increased mortality may be caused by hydraulic failure and reduced production of non-structural C [28]; process-based models covering these processes are being developed rapidly and will become popular in the near future. In this study, our general objective is to contribute to the development of adaptation strategies for climate change for the forestry sector in East Asia. As a specific objective, we parameterize a process-based model for planted forests in Japan at the national scale to quantify the climatic limitations on stand productivity and detect regions with possible declines in NPP, with a focus on Japanese cedar as a representative tree species that grows in a warm and humid region of East Asia.

Materials and methods

Biome-BGC (BBGC)

We used the process-based ecosystem model BBGC ver.4.2 [29] as the basis for modeling C dynamics. BBGC is a two-leaf model that simulates the stocks and fluxes of C, N, water and energy. This model requires input data on daily maximum, minimum, and mean air temperatures; precipitation; daytime VPD and solar radiation; and day length. BBGC also requires eco-physiological data (see section Parameterization) and site data (see section Site data) to run its simulations. The NPP is calculated by gross primary production (GPP) minus autotrophic respiration in daily routine. GPP is calculated using Farquhar’s model [24] connected with the stomatal conductance model [30]. In this routine, the function of the maximum rate of carboxylation (Vcmax) is associated with the parameterization conducted in this study, as follows:(1)where FLNR is the fraction of leaf N in Rubisco, ACT is the activity of Rubisco scaled based on temperature and the concentrations of oxygen (O2) and carbon dioxide (CO2), and CN_leaves is the C:N ratio in leaves. The photosynthesis rate is reduced with four scaling factors—namely, air temperature, VPD (S1 Fig in S1 File), soil water potential (S1 Fig in S1 File), and solar radiation. Autotrophic respiration comprises maintenance respiration, which is a function of temperature and the N pool of the living tree, and growth respiration, which is a fixed proportion of GPP. The allocation of C to tree organs, such as leaves, stems, coarse roots, and fine roots, is determined as a fixed ratio based on eco-physiological parameters. Other processes such as the decomposition of woody debris and soil organic matter have been described in detail previously [29]. In this study, we use BBGC as applied in many previous studies of planted forests to estimate changes in NPP based on a simple physiological response without acclimation to climatic factors, i.e., increased temperature and elevated CO2 concentration.

First, we made minor changes to the BBGC code. The turnover rate of fine roots was separated from that of leaves, as in previous studies [6, 7, 31]. In addition, the function driving the soil water retention curve, which represents the relationship between soil water content and soil matric potential, was replaced; the original pedo-transfer function based on soil texture (i.e., sand, silt, and clay content) was switched to the function of van Genuchten’s model [32]. This change affected the estimation of GPP, in combination with site data (see section Soil water retention curve), due to the reduction function for stomatal conductance associated with soil water potential (S1 Fig in S1 File).

Plant trait database (PTDB)

We used a PTDB for Japanese cedar [15] for the parameterization of BBGC. The PTDB contains 16,410 data entries for 177 types of plant traits in Japanese cedar, which were collected through review of previous studies published since 1950. The PTDB was also used for the validation of modeled NPP through extraction of observed NPP values from the database (see section Statistics and validation).

Flux data for Japanese cedar plantations

In addition to the PTDB, monitored and calculated flux values for Japanese cedar at two eddy covariance monitoring sites, the Takayama coniferous forest (TKC) and Kahoku experimental watershed (KHW), were used for the parameterization of BBGC.

Takayama coniferous forest (TKC).

The TKC is located in a cool temperate region (36.140 deg N, 137.371 deg E; Fig 1A). The tree plantation was established in 1961 or later. The mean annual temperature is 10.7 deg C (from 1967 to 2006) and the mean annual precipitation is 1722 mm (from 1967 to 2006). The plant area index ranges from 4.9 to 5.2 m2 m-2. Further site information is available from the Asiaflux website (http://asiaflux.net/index.php?page_id=111). The monitoring data used in this study include daily datasets of climate, GPP, ecosystem respiration (RE), and net ecosystem change (NEE) estimated using the eddy covariance method for 2006 (stand age of 46 years). For more detailed information about the flux calculation and gap-filling processes, see Saitoh et al. [33]. Climatic input data for BBGC over 45 years (1961–2005) were obtained through correction of data from the nearest meteorological station.

[Figure omitted. See PDF.]

Fig 1. Maps of regional blocks, flux-monitoring sites, and climate in Japan.

(A) Locations of the three blocks and two flux-monitoring sites, and (B) mean annual temperature (MAT) and (C) mean annual precipitation (MAP) in the 1996–2000 period. Red, yellow, and green in map (A) indicate the southwestern (SW), central (CT), and northeastern (NE) blocks, respectively. Maps (B) and (C) were constructed using historical data from MIROC5 (Model for Interdisciplinary Research On Climate Version Five).

https://doi.org/10.1371/journal.pone.0247165.g001

Kahoku experimental watershed (KHW).

The KHW is located in a warm temperate region (33.137 deg N, 130.710 deg E; Fig 1A). The tree plantation was established in 1956. The mean annual temperature is 15.1 deg C (from 2000 to 2003) and the mean annual precipitation is 2106 mm (from 2001 to 2003). The leaf area index as measured with a LAI-2000 plant canopy analyzer (Li-Cor, Inc., Lincoln, NE) ranges from 3.6 to 5.2 m2 m-2. Further information is available from the FFPRI FluxNet Database website [34]. The monitoring data used in this study were daily datasets of climate (2000–2003) and NEE (2001–2003, stand age of 46–48 years). Because Japanese cypress grows in part of the watershed, we extracted flux mainly from the Japanese cedar stand based on wind direction (from the right side of the site), as described by Shimizu et al. [35]. Then, we conducted gap-filling through the mutual imputation method (NORM software) using in-situ measurements of meteorological factors (solar radiation, net radiation, air temperature, VPD, and wind speed; [35, 36]) and generated daily NEE for the Japanese cedar stand. As in the TKC, NEE in the KHW was separated into GPP and RE. We averaged daily GPP to generate monthly data and used these data for parameterization, as the calculated response of Japanese cedar to climate in the KHW was considered less reliable at the daily timescale according to the gap-filling procedure described above. In addition to monthly GPP, RE, and NEE, monthly soil respiration measurements were up-scaled using soil temperature and the site-specific regression model [37]. Climatic input data for BBGC over 21 years (1978–1999) were obtained and corrected using data from the nearest meteorological station.

Parameterization

We varied 12 of 44 eco-physiological parameters in BBGC in two steps, while considering the balance between simplicity and effectiveness of parameterization. STEP-0 uses the parameter set prepared for evergreen needle-leaf forest in BBGC ver.4.2 (Table 1). In STEP-1 and -2, we mainly updated the eco-physiological parameters of leaves, which are the basis of the two-leaf model BBGC and control photosynthesis at the daily timescale. On the other hand, the eco-physiological parameters of turnover, except for fine roots and the chemical composition of tree organs, which control processes at the monthly or annual timescales, remained unchanged from STEP-0 in this study, mainly due to sparse data in the PTDB. In addition, parameters related to C allocation, the mechanism of which remains unclear for Japanese cedar, were also unchanged.

[Figure omitted. See PDF.]

Table 1. Eco-physiological parameters used in BBGC.

https://doi.org/10.1371/journal.pone.0247165.t001

STEP-1.

Nine eco-physiological parameters were varied in STEP-1. The annual fire mortality was changed from 0.005 to 0. We assumed that the risk of forest fire in Japan is quite low relative to other countries based on the long-term trend of forest fires [38]. In particular, no fire history has been recorded at the TKC or KHW. Then, the turnover of fine roots was changed from 0.25 to 0.69 year-1 (value of Konôpka cited in Noguchi et al. [39]), considering its effect on biomass accumulation [31].

We also changed five parameters related to the C:N ratio, SLA, and Max_sc (Table 1). These parameters are sensitive to NPP [40] and have relatively rich data in the PTDB. We used median values in the PTDB, as some parameters were ratios or exhibited log-normal distributions in the PTDB (S2 Fig in S1 File). The C:N ratios of leaf litter after retranslocation and of live wood were scaled to those of leaves and dead wood.

STEP-2.

Three additional eco-physiological parameters were optimized in STEP-2. We applied BO using the algorithms employed for Bayesian calibration [16] with slight modification. In BO, the likelihood function p (E) was defined for each site, the TKC and KHW, as follows:(2)where E is the difference between the observed data, D, and output of BBGC, M, based on the vector of three eco-physiological parameters, θ. The term i is an index of the observed data and the output of BBGC, n is the number of observational data, φ denotes a Gaussian probability density function with a given mean and variance, and Dave is the average D for the TKC or KHW. Then, we combined the p(E) of the TKC and KHW and estimated a vector for the parameters that provide maximum likelihood. As described below, we constrained the model using GPP at the TKC and KHW simultaneously and weighted it using two likelihood functions to avoid being affected more strongly by the TKC, for which we have more data than for the KHW. The Bayesian calibration procedure was described in detail by Van Oijen et al. [16].

GPP in the TKC (daily data in 2006, 365 count) and KHW (monthly data in 2001–2003, 36 count) were used as constraints on BO. In this study, we focused on the improvement of modeled GPP as a C input to the ecosystem, and only validated the modeled RE (also, soil respiration in the KHW) and NEE values through comparisons with calculated RE and observed NEE, respectively. As RE is related to the decomposition process of organic matter over several decades and is affected by many eco-physiological parameters for which little information is available, we found that optimization using RE data complicated the procedure of parameterization, even when the effective approach of a hierarchal parameterization scheme was used [41].

We optimized the fraction of leaf N in Rubisco (FLNR) and VPD values for starting and ending the reduction of stomatal conductance (VPD_rsc_start and VPD_rsc_end, respectively). Although the PTDB contains little or no data for these parameters, we consider these parameters essential for the improvement of modeled GPP. The prior distribution in BO was uniform, and the range of FLNR was 25–200% of 0.04 (0.01–0.08), the value in STEP-0; those of VPD_rsc_start and VPD_rsc_end were 10–100% of 930 Pa (93–930 Pa) and 25–100% of 4100 Pa (1025–4100 Pa), respectively. We set these ranges of VPD_rsc_start and VPD_rsc_end for BO based on a previous report stating that canopy conductance of Japanese cedar can decrease when VPD is much lower than 1.0 kPa [42]. In MCMC sampling for BO, three chains of 100,000 iterations were calculated, starting from the 25th, 50th, and 75th percentiles of the prior distribution.

Site data

We prepared site data as an ini file in BBGC. The mesh for calculation contains forest areas larger than 0.4 km2 (approximately 40% of the total area in the mesh) based on information in the National Forest Resources Database ([43], compiled at the third mesh order). The mesh for Hokkaido, the northernmost prefecture of Japan, was excluded in this study because plantations of Japanese cedar are very scarce in that region, except for a small proportion of the southern area, and because few observational data of NPP and biomass are available for comparison with model outputs. We set shortwave albedo to a uniform value of 0.15 in the ini file.

Soil thickness.

We prepared a map of the thicknesses of the A and B soil horizons and used their sum as soil thickness in the ini file in BBGC (S3 Fig in S1 File). The data source is soil profiles at 2056 points in Japan, obtained from National Forest Soil Carbon Inventory Project [44] and a legacy dataset compiled by Morisada et al. [45]. We first extracted the maximum soil thickness within the area represented as a coarse mesh of 15-km resolution with regional-scale values of soil thickness (N = 671). The map of soil thickness at the third mesh order with 1-km resolution was estimated through ordinary kriging from regional-scale data using the Geostatistical Analyst Extension of ArcGIS Pro 10.5 software (ESRI Inc., USA). Thus, the interpolated map shows the maximum potential soil thickness at the regional scale (ca. 15 km) in Japanese forests.

Soil water retention curve.

We prepared the parameters related to the soil water retention curve of van Genuchten ([32]; called vG parameters in this study) for each soil type based on the current forest soil classification system (FSD system) for Japan [46]. The equation of van Genuchten is as follows:(3)where θ(ψ) is the volumetric soil water content (m3 m-3) as a function of ψ, the soil matric potential (MPa). θs and θr are the saturated and residual soil water contents in volume, respectively, and α, and n, and m (= 1 − 1/n) are curve-fitting parameters. We first converted the vector soil map at a scale of 1:200,000 [47] into a raster soil map of the third mesh order (approximately 1-km resolution) based on the soil type with maximum coverage. The soil type in the map of the National Land Agency, which differs from that of the FSD system, was converted to align with the FSD system using the criteria proposed by Morisada et al. [45] to allocate the representative vG parameters of the A and B horizons to the map. Representative vG parameters of the A and B horizons of eight soil types were obtained through fitting of the vG function to soil water retention curve datasets from previous studies, which were allocated to 13 soil types (S4 Table in S1 File). Then, the vG parameters of the A and B horizons were weight-averaged based on soil thickness. To express the soil water content of the entire soil profile at a certain value of soil matric potential, arithmetic and geometric averaging [48] were used for θs and θr, and for α and n, respectively. The map of effective soil water capacity in the soil profile, calculated based on soil thickness, θs and θr, is shown in S3 Fig in S1 File.

Deposition and biological fixation of nitrogen.

In this study, differences in the input of N among grids were expressed based on differences in the deposition of N (S3 Fig in S1 File). We obtained the sum of dry and wet deposition of N in a 15-km grid in 2012 from data simulated using an atmospheric chemical transport model [49], which was used in Nishina et al. [50]. Deposition data at the third mesh order with 1-km resolution were interpolated through ordinary kriging.

On the other hand, no spatial information on the biotic fixation of N was available, although it occurs in the litter layer of Japanese cedar plantations [51, 52]. Therefore, we assigned a relatively low constant rate of biotic N fixation across all regions (0.40 gN m-2 year-1) considering the ranges associated with lichen (0–0.20 gN m-2 year-1), bryophytes (0.07–1.00 gN m-2 year-1), and litter (0.10–0.20 gN m-2 year-1) in temperate forests [53].

Spin-up and normal run procedures

In spin-up runs for the TKC and KHW sites, we used the eco-physiological parameters of STEP-0, set the atmospheric CO2 concentration to the value for the first year of climate data, and fixed the N deposition rate to 0.1 gN m-2 year-1. To determine the initial conditions of the normal runs for the TKC and KHW, we cut 100% of the biomass C in existing forest generated during the spin-up run, removed above-ground biomass C, and increased leaf C by 0.1 kgC m-2. In the normal runs for the TKC and KHW sites, we set the N deposition rate to 1.1 gN m-2 year-1. For the TKC, climate data for 39 years before 2000 (1961–1999) were used in the spin-up run and those for 46 years (1961–2006) were used in the normal run. For the KHW, climate data for 21 years before 2000 (1978–1999) were used in the spin-up run, those for 26 years (1978–2003) were used in the normal run, and climate data for the 1978–1999 period were applied to the young stand (1956–1977).

For the initial conditions of the normal run at the national scale, we set the cutting (and removal) ratio of biomass C in the existing forest to 90%. The cutting ratio has been controlled in previous studies to explore the effect of forest disturbance on simulated flux [9, 29]. This ratio was determined practically in this study based on preliminary analysis, and a cutting ratio lower than 100% improved the stability of the calculation through increased leaf biomass at the initial growth stage, especially in severely cold regions, although the percentage of such regions in the total mesh area is quite small. We confirmed that changing the cutting ratio of biomass C in existing forest scarcely affected the modeled NPP and increased rate of vegetation biomass, although the baseline level of vegetation biomass increased (see section of Result, Recorded and modeled NPP and biomass C relative to stand age). Other forest management options, such as thinning, were not considered in the normal run because they were too spatially diverse to capture with 1-km resolution and were outside of the objectives of this study.

Climate scenarios

We used climate scenarios that were statistically downscaled, i.e., Inverse Distance Weighted method, by NARO [27] (Nishimori et al., 2019; NAROv2.7r). The original dataset is the product of phase 5 of the Coupled Model Intercomparison Project (CMIP5). These scenarios were generated using five global climate models (GCMs; Table 2) under two representative concentration pathways, RCP2.6 and RCP8.5, following the conventions of the Intergovernmental Panel on Climate Change [54]. RCP2.6 is a representative scenario based on low emissions of greenhouse gasses (GHGs) that aims to keep global warming below 2 degrees C above pre-industrial temperatures [54]. Meanwhile, RCP8.5 is based on very high GHG emissions, higher than those in scenarios without additional efforts to constrain GHG emissions (baseline scenarios). The downscaled dataset contains climate scenarios based on historical trends (1981–2005) and future predictions (2006–2100) [27]. The atmospheric CO2 concentration for each RCP is shared among GCMs and was obtained from the RCP Database [55, 56]. The atmospheric CO2 concentration in 2000 was 368.9 ppm and those in 2100 under RCP2.6 and RCP8.5 are 420.9 and 935.9 ppm, respectively.

[Figure omitted. See PDF.]

Table 2. The five GCMs used for climate prediction in this study.

https://doi.org/10.1371/journal.pone.0247165.t002

In BBGC, some climate input data are based on daytime-average values rather than daily-average values. Therefore, we used daily maximum and minimum air temperatures (Tmax and Tmin, deg C) and precipitation (mm day-1) in the climate scenarios in this study, and employed MT-CLIM module ver.4.3 [57, 58] to generate the remaining input data for BBGC, i.e., daily mean air temperature (deg C), daytime VPD (Pa), solar radiation (W m-2), and day length (sec). In our exploratory analysis, the spatiotemporal trend of daytime VPD from the MT-CLIM module was similar to that of daily VPD estimated from daily relative humidity in the climate scenarios of MIROC5 (Model for Interdisciplinary Research On Climate Version Five), with both VPDs increasing in summer, especially in the western part of Japan.

Here, we mainly summarize the characteristics of the temperature increases according to five GCMs for the target area (see section Target stand age and area for prediction). In this study, we considered the outputs for 1996–2000 as the historical trend (H_2000) and those for 2096–2100 as future predictions (F_2100). The mean annual temperature (MAT) for forested areas from the five GCMs for H_2000 was 11.3 deg C, with high values in the southwestern region (Table 3, Fig 1B) and little variation among the GCMs. The increases in MAT (ΔMAT) from H_2000 to F_2100 were 1.8 and 5.3 deg C for RCP2.6 and RCP8.5, respectively, averaged over the five GCMs (Table 3). The ΔMAT from MIROC5, which is well-documented in this study, was the second highest for RCP2.6 and fourth highest for RCP8.5 (Fig 2). The mean annual precipitation (MAP) from the five GCMs for H_2000 was 2127 mm year-1 (Table 3, Fig 1C) and increased along with MAT in future scenarios. The mean ΔMAP in MIROC5 was relatively low for RCP2.6 and high for RCP8.5 relative to those of the other GCMs (Table 3). The increasing trends for daytime VPD from the five GCMs, calculated using MT-CLIM, were similar to those for MAT, and are shown in S5 Table in S1 File along with the trends for solar radiation.

[Figure omitted. See PDF.]

Fig 2. Trends in ΔMAT according to five GCMs.

Comparison between H_2000 and F_2100.

https://doi.org/10.1371/journal.pone.0247165.g002

[Figure omitted. See PDF.]

Table 3. MAT and MAP and their changes in future scenarios.

https://doi.org/10.1371/journal.pone.0247165.t003

Target stand age and area for prediction

We estimated the average NPP for 36–40-year-old stands to reduce yearly fluctuations caused by climate data, as the age of 40 years is commonly used as a site index in Japan. In the calculation for H_2000 (1996–2000 for 36–40-year-old stands), climate data for 1981–2000 were used twice to simulate a 40-year-old stand, whereas in the calculation for F_2100 (2096–2100 for 36–40-year-old stands), climate data for 2061–2100 were used once, except in HadGEM2-ES, which contains climate data for 2060–2099. The ΔNPP of each point of mesh was calculated as follows:(4)where NPPF_2100 and NPPH_2000 are mean NPP for F_2100 and H_2000, respectively. In total, 196,928 points of third mesh order were used for prediction (BBGC_full), covering all of Japan except Hokkaido. We divided this area into three blocks (Fig 1A)—namely, the southwestern (SW; 33% of total forest area, covering the Kyushu, Okinawa, Chugoku, and Shikoku Japanese administrative districts), central (CT; 35%, Kinki, Tokai, and Hokuriku districts), and northeastern (NE, 32%, Kanto and Tohoku districts) blocks. In addition, we used a dataset of 1967 points (BBGC_sample) that were systematically sampled from BBGC_full to analyze the differences between relatively warm (May to October) and cool (November to April) seasons, those among the five GCMs, and the effects of climate factors on GPP and NPP. We confirmed that the probabilistic distribution of NPP in BBGC_sample was approximately the same as that in BBGC_full (S6 Fig in S1 File), indicating that the use of BBGC_sample had little impact on the results and conclusions of this study.

Statistics and validation

The normalized root mean square error (NRMSE) for the model output based on eco-physiological parameters in STEP-0, -1, and -2 was calculated as follows:(5)where, as in Eq (2), D and M are the observed value and modeled output from BBGC, respectively. The term i is an index of the observed data and the output of BBGC, and n is the number of observations. Dmax and Dmin are the maximum and minimum observed values, respectively, across the TKC or KHW datasets.

The NPP of planted forests, not only for Japanese cedar but in general, can be affected by growth stage [59]. Therefore, we validated the stand age used for estimating NPP, i.e., 36–40 years, by extending the stand age from 40 to 60 years in H_2000 using climate data from MIROC5 for 1981–2000 three times. We compared the trend in mean NPP with age between the model output and observations. We extracted observed NPP values from the PTDB [15] by filtering for a stand age of 31–50 years. We also analyzed the increases in biomass C with age and extracted wood volume data from the National Forest Inventory (NFI) for 2009–2013 by filtering for a stand age of 30–60 years without thinning history. Then, we converted the volume to biomass C using a biomass expansion factor (BEF; [60]). In addition, as one step in the validation of the estimated range of NPP, we compared the NPP from H_2000 in STEP-2 between cedar-dominated and non-cedar-dominated mesh points based on data in the National Forest Resources Database [43] (see S7 Fig in S1 File), although the target of this study is the total forested area in Japan. This comparison was undertaken based on the idea that cedar-dominated areas are currently under environmental conditions that favor the growth of Japanese cedar plantations relative to those of non-cedar-dominated areas.

In this study, we did not conduct uncertainty analysis, which would require the probabilistic density of all parameters but instead analyzed the sensitivity of modeled NPP to selected eco-physiological and site parameters. The result of sensitivity analysis is described in supporting information files (S8 Text in S1 File).

Results

Simulated flux before and after parameterization

The modeled fluxes, including GPP and others, improved after two steps of parameterization. The parameters before and after optimization are summarized in Table 1. Three parameters in BO—namely, FLNR, VPD_rsc_start, and VPD_rsc_end—converged (S9 Fig in S1 File). The NRMSE values of modeled GPP were 0.25 and 0.53 in STEP-0 for the TKC and KHW, respectively, which increased slightly in STEP-1 relative to STEP-0 and then decreased in STEP-2 to lower than half the STEP-0 value (Table 4, Figs 3A and 4A). The other model outputs, namely RE, soil respiration, and NEE, which is the only flux with observational data, also improved in STEP-2 (Table 4, Fig 3B and 3C, and Fig 4B–4D).

[Figure omitted. See PDF.]

Fig 3. Comparison of simulated flux at the TKC before and after parameterization.

(A) Gross primary production (GPP), (B) ecosystem respiration (RE), and (C) net ecosystem exchange (NEE). Gray solid, gray dashed, and red lines are modeled values based on eco-physiological parameters from STEP-0 (parameter set for evergreen needle-leaf forest), and those modified in STEP-1 and STEP-2, respectively.

https://doi.org/10.1371/journal.pone.0247165.g003

[Figure omitted. See PDF.]

Fig 4. Comparison of simulated flux at the KHW before and after parameterization.

(A) Gross primary production (GPP), (B) ecosystem respiration (RE), (C) net ecosystem exchange (NEE), and (D) soil respiration. Gray solid, gray dashed, and red lines are modeled values based on eco-physiological parameters from STEP-0 (parameter set for evergreen needle-leaf forest) and those modified in STEP-1 and STEP-2, respectively.

https://doi.org/10.1371/journal.pone.0247165.g004

[Figure omitted. See PDF.]

Table 4. Normalized RMSE for modeled flux for the TKC and KHW sites.

https://doi.org/10.1371/journal.pone.0247165.t004

The modeled GPP in the TKC exhibited a sharp decrease followed by an increase in the daily time step during summer in STEP-0, and these variations were removed in STEP-2 (Fig 3A). The modeled RE in the TKC was underestimated in summer, even in STEP-2 (Fig 3B). As for the TKC, the modeled GPP and RE for the KHW exhibited sharp changes with variations in monthly climatic conditions in STEP-0 (Fig 4A and 4B). The effects of parameterization in STEP-1 on GPP and RE in the KHW were small, whereas modeled soil respiration improved. The modeled GPP and RE values became closer to estimates in the KHW in STEP-2, although they were slightly overestimated in summer.

NPP and ΔNPP

Estimated historical trends in NPP.

The annual NPP from the five GCMs for H_2000 was 0.77 ± 0.10 kgC m-2 year-1 (mean ± standard deviation [s.d.]; BBGC_sample) with the highest values in the SW block (Table 5). The annual NPP from MIROC5 for H_2000 was 0.76 ± 0.10 kgC m-2 year-1, and the trends in regional blocks were quite similar to those generated in the five GCMs (Table 5, Fig 5A). The range of annual NPP from MIROC5 at STEP-2 was significantly higher in currently cedar-dominated mesh points (0.80 ± 0.06 kgC m-2 year-1, 26% of total area) than in non-cedar-dominated (0.75 ± 0.11 kgC m-2 year-1, 74% of total area) forests (p <0.001; Fig 6).

[Figure omitted. See PDF.]

Fig 5. Estimated NPP and ΔNPP of Japanese cedar.

(A) NPP for H_2000 from MIROC5, (B) ΔNPP from H_2000 to F_2100 from MIROC5 under RCP2.6, (C) ΔNPP from H_2000 to F_2100 from GFDL-CM3 under RCP2.6, and (D) the NPP-decline area extracted from maps (B) and (C). The blue area of map (D) is based on ΔNPP < −0.02 kgC m-2 year-1 as indicated in map (B), the yellow area is based on ΔNPP < −0.075 kgC m-2 year-1 as indicated in map (C), and the red area indicates the overlap area of these two conditions.

https://doi.org/10.1371/journal.pone.0247165.g005

[Figure omitted. See PDF.]

Fig 6. Distribution of modeled NPP in cedar-dominated and non-cedar-dominated areas.

https://doi.org/10.1371/journal.pone.0247165.g006

[Figure omitted. See PDF.]

Table 5. NPP, ΔNPP, and proportion of NPP-decline area.

https://doi.org/10.1371/journal.pone.0247165.t005

The seasonal NPP estimates from the five GCMs for H_2000, in daily terms, were 3.21± 0.47 and 0.98 ± 0.50 gC m-2 day-1 during the warm and cool seasons, respectively, and the mean values in three regional blocks, SW, CT, and NE, were 3.05, 3.24, and 3.35 gC m-2 day-1 (mean from five GCMs) for the warm season and 1.40, 0.93, and 0.59 gC m-2 day-1 for the cool season, respectively. These seasonal trends in NPP from MIROC5 were similar to those obtained from the five GCMs.

Estimated ΔNPP under the RCP2.6 and RCP8.5 scenarios.

The changes in NPP (ΔNPP) according to the five GCMs were 0.04 ± 0.07 and 0.14 ± 0.11 kgC m-2 year-1 (mean ± s.d.; BBGC_sample) under the RCP2.6 and RCP8.5 scenarios, respectively, and the ΔNPPs were lower in the SW region than in the CT and NE regions (Table 5). Decreases in NPP (negative ΔNPP) in RCP2.6 were predicted more frequently with MIROC5 (27% of total area) and GFDL-CM3 (67%, Fig 7), and the NPP-decline area was distributed primarily in the western part of Japan (Fig 5B–5D). Negative ΔNPP under RCP8.5 was predicted solely with GFDL-CM3 (Fig 7).

[Figure omitted. See PDF.]

Fig 7. Trends in ΔNPP according to five GCMs.

Comparison between H_2000 and F_2100. The outliers of ΔNPP are not shown in the figure, as the resulting wide ranges obscure the differences among GCMs.

https://doi.org/10.1371/journal.pone.0247165.g007

The mean ΔNPP from the five GCMs under RCP2.6 was negative and positive in the warm and cool seasons, respectively, and was especially low in the warm season in the SW block (Table 6). The mean ΔNPP from the five GCMs under RCP8.5 was much higher than that under RCP2.6 during the cool season.

[Figure omitted. See PDF.]

Table 6. Seasonality of ΔNPP in F_2100.

https://doi.org/10.1371/journal.pone.0247165.t006

Recorded and modeled NPP and biomass C relative to stand age

A total of 48 records of NPP for 31–50 year-old stands, aged 40.3 years on average, were extracted from among 116 records for 5–70 year-old stands in the PTDB. The NPP for 31–50-year-old stands was 0.78 ± 0.24 kgC m-2 year-1 and exhibited a decreasing trend with stand age (p <0.05; Fig 8A). Meanwhile, the mean value of modeled NPP in STEP-2, based on the historical trend estimated in MIROC5, ranged from 0.75 to 0.78 kgC m-2 year-1 for every 5-year bin from 31 to 50 years of age, and its change with stand age was unclear (Fig 8B).

[Figure omitted. See PDF.]

Fig 8. Changes in observed and modeled NPP with stand age.

(A) Observed NPP. Open circles indicate records in the PTDB, and the solid line represents the regression with stand age. (B) Modeled NPP. The solid lines show the mean and standard deviation in the total area for every 5-year bin from 31 to 50 years. SW, CT, and NE indicate mean values of the southwestern, central, and northeastern blocks, respectively.

https://doi.org/10.1371/journal.pone.0247165.g008

In total, 1094 records of wood volume for 30–60-year-old stands were extracted from the NFI. The temporal trends in stand biomass C, which was converted from wood volume based on a BEF, were determined for every 10-year period as the coefficient of a single regression with stand age; these values were 0.279, 0.303, and −0.01 kgC m-2 year-1 for stand ages of 30–40, 40–50, and 50–60 years, respectively (Fig 9A). The differences in stand biomass C among regional blocks were unclear in this dataset, except that higher biomass was observed in the SW region relative to the CT and NE regions at age 30–40 years (p <0.001). Meanwhile, the mean increments of modeled stand biomass C were 0.388, 0.355, and 0.350 kgC m-2 year-1 for stand ages of 30–40, 40–50, and 50–60 years, respectively, indicating a relatively monotonic increase in stand biomass (Fig 9B). The modeled stand biomass C was generally higher in the SW block compared to the CT and NE blocks.

[Figure omitted. See PDF.]

Fig 9. Changes in estimated and modeled biomass C with stand age.

(A) Estimated biomass C based on wood volume in the NFI. The solid line represents the regression of biomass C with age for 10-year bins. (B) Modeled biomass C. The black solid line indicates the mean and standard deviation of biomass C in the total area at 10-year intervals. The gray dotted line indicates mean biomass C based on a cutting ratio of 99% in existing forests. SW, CT, and NE indicate mean values of the southwestern, central, and northeastern blocks, respectively.

https://doi.org/10.1371/journal.pone.0247165.g009

Discussion

Validity and limitations of estimating NPP through two-step parameterization

The parameterization conducted in this study improved the estimation of the physiological responses of Japanese cedar plantations growing under humid climatic conditions in East Asia. We confirmed that changes in eco-physiological parameters in STEP-1, which were mainly based on PTDB data, including the C:N ratio of leaves, SLA, and turnover of fine roots, rarely contributed to improving the simulated seasonality of GPP in daily or monthly outputs for the TKC and KHW. Meanwhile, changes in FLNR and the VPD-stomatal relationship due to BO in STEP-2 sharply improved simulated GPP at those sites. As a result, our estimation of annual NPP in the TKC in 2006 approached 0.79 kgC m-2 year-1, which is the biometric estimate over 4 years (2005–2009) as reported by Yashiro et al. [61], rising from 0.65 kgC m-2 year-1 in STEP-0 to 0.73 kgC m-2 year-1 in STEP-2. The VPD-stomatal relationship may be the key factor for simulating seasonality rather than FLNR, which instead altered the overall range of NPP [40]. The range of the VPD-stomatal relationship in STEP-2 was 135–1824 Pa, which included the entire ranges of mean daytime-VPD in the warm seasons of 302–1257 Pa for H_2000 and 428–1497 Pa for F_2100 from MIROC5 under RCP8.5. This range is quite low compared to that presented in BBGC for Japanese cypress (610–4100 Pa; [9]), which exhibits more water-saving characteristics than Japanese cedar when the soil water supply is limited [62], making the modeled stomatal behavior of Japanese cedar more sensitive to increases in VPD within the climatic range of Japan. The observed range of the VPD-stomatal relationship at the canopy level ([42, 63–65] supports adjustment to a substantially lower range of VPD compared to that in STEP-0 (930–4100 Pa). We optimized these parameters by simultaneously improving modeled GPP at two sites under different climatic conditions and timescales, the TKC and KHW, and showed that the weighting of the likelihood function in BO was effective.

The modeled NPP results for H_2000 were considered valid from the perspective of their geographical distribution. The mean NPP for H_2000 for stands aged 36–40 years was similar to that for stands aged 31–50 years in the PTDB (Fig 8A and 8B). The NPP and stand biomass C for 36–40-year-old stands for H_2000 were higher in the SW and lower in the NE region (Figs 5A, 8B and 9B). This result was consistent with statistical trends showing that the growth rate in terms of height of Japanese cedar is high in the southern and Pacific Ocean coastal regions of Japan, especially during the initial stage of growth [18]. Ito [10] estimated NPP in various biomes in East Asia using a C cycle model and obtained high values for evergreen needle-leaf forests in the southern and coastal regions of Japan. Hashimoto et al. [66] inversely estimated NPP in forested areas of Japan, not only for Japanese cedar forests, using a soil C stock dataset and the RothC model, and obtained higher values in southwestern regions compared to northeastern regions. Mitsuda [5] estimated NPP in Japanese cedar plantations on Kyushu Island (covering approximately 40% of the SW block) using the 3-PG model and obtained high values in the southern area, in accordance with this study. Regardless of the spatial validity of estimated NPP, it is important to note that we might not have sufficiently simulated the age dependency of NPP. Although the observed NPP values in the PTDB are not successional data for a single stand, they exhibited a decreasing trend with stand age from 31 to 50 years (Fig 8A), whereas the average NPP values for H_2000 and F_2100 at 5-year intervals exhibited a relatively monotonous trend over that period (Fig 8B). The estimated biomass C from NFI exhibited an increasing trend with stand age for at least up to 50 years, which probably includes the negative effect of harvesting at higher stand ages (Fig 9A). Meanwhile, the modeled biomass C in this study increased monotonically from 30 to 60 years of stand age (Fig 9B). In summary, we could not validate the modeled high growth rate of old Japanese cedar in stands >50 years of age, and it would be inappropriate to apply the model parameterized in this study to stands much older or younger than 40 years. The age dependency of NPP is a sensitive factor that complicates the modeling and assessment of climate impacts on planted forests. In this study, the weak dependency of NPP on stand age enabled comparison of NPP at a fixed age, namely 36–40 years, between current and future climatic conditions.

Difference in estimates of ΔNPP between the RCP2.6 and RCP8.5 scenarios

The mean ΔNPP for F_2100 was generally lower under the RCP2.6 scenario than under RCP8.5. This difference was likely due to mitigation of the decrease in modeled GPP during warm seasons related to the high CO2 concentration in RCP8.5. The modeled relationship between GPP and Tday (mean air temperature in daytime) during the cool season and that between NPP and Tday exhibited similar trends of temperature dependency for both H_2000 and F_2100 (based on MIROC5; Fig 10C and 10D). On the other hand, the GPP-Tday and NPP-Tday relationships during the warm season shifted upward from the historical climate to RCP2.6 and further to RCP8.5 (Fig 10A and 10B). In an additional test of the RCP8.5 climate scenario with the CO2 concentration reduced to the level of RCP2.6, the modeled GPP and NPP for F_2100 were sharply reduced, and the percentage of the total area with negative ΔNPP increased from 0% under the original RCP8.5 conditions to 45%. Similar results have been reported in previous studies of coniferous forests in East Asia, mainly in China. When the effect of elevated CO2 was considered, increased NPP was predicted more frequently under RCP8.5 than under RCP2.6 for tree species such as Pinus tabulaeformis [67] and Abies fabri [68], although if that effect was neglected, the risk of decreased productivity became high under RCP8.5, as noted for Larix olgensis [69] and A. fabri [68]. A positive effect of elevated CO2 on GPP could be detected at the global scale, even under the current climatic conditions [70]. Down-regulation of photosynthesis under elevated CO2 conditions, leading to an observed rate lower than that estimated by Farquhar’s model, has not been confirmed for mature stands of Japanese cedar. However, Hiraoka et al. [71] reported that the photosynthetic rate of 1-year-old cuttings of Japanese cedar increased under elevated CO2 conditions, whereas the maximum carboxylation rate decreased. Thus, the increase in NPP under highly elevated CO2 conditions may be lower than that predicted in this study. To account for this down-regulation in future modeling of photosynthesis, reduction of Max_sc through scaling of CO2 concentration ([72, 73]) may be an effective strategy, although the sensitivity analysis in this study indicated that reducing Max_sc did not necessarily decrease modeled NPP under the RCP8.5 scenarios (S8 Text in S1 File).

[Figure omitted. See PDF.]

Fig 10. Relationships of GPP and NPP with daytime-mean air temperature.

Tday, daytime-mean air temperature. (A) GPP in the warm season, (B) NPP in the warm season, (C) GPP in the cool season, and (D) NPP in the cool season.

https://doi.org/10.1371/journal.pone.0247165.g010

Geographic characteristics of the modeled decrease in NPP in Japanese cedar plantations

Here, we focus on the area of decreased NPP (NPP-decline area) under the RCP2.6 scenario, considering the potential risk of decreased productivity of Japanese cedar plantations. The NPP-decline area in this study was concentrated in the SW region (Fig 5D), in accordance with previous studies reporting high vulnerability of Japanese cedar in western Japan to a slowly warming climate [5, 22]. The geographic distribution of negative annual ΔNPP was derived from the decrease of NPP in the warm season, rather than its increase in the cool season (Table 6). The area that is sensitive to increasing temperature, evaluated based on the ratio of ΔNPP to ΔMAT, was characterized mainly by MAT in the historical dataset. Negative values of ΔNPP/ΔMAT were predicted more frequently by two GCMs, MIROC5 and GFDL-CM3, in this study, and were dominant when MAT for H_2000 exceeded approximately 13 deg C (Fig 11), whereas their relationship with MAP was unclear. Risk assessment of the NPP-decline area at the continental scale showed that the risk tended to be high in the southern area, where the level of NPP in the current climate is high [74–76].

[Figure omitted. See PDF.]

Fig 11. Distribution of the ΔNPP-to-ΔMAT ratio in relation to MAT and MAP under the current climate.

Results from (A) MIROC5 and (B) GFDL-CM3, based on the RCP2.6 scenario. R represents the ΔNPP-to-ΔMAT ratio. Red (R < −0.015) and yellow (−0.015 ≤ R <0) indicate decreases in NPP with increased temperature, whereas green (0 ≤ R < 0.015), light blue (0.015 ≤ R < 0.03), and blue (R ≥ 0.03) indicate increases in NPP.

https://doi.org/10.1371/journal.pone.0247165.g011

In summary, we first evaluated the vulnerability of Japanese cedar plantations to a changing climate by modeling stand productivity at the national scale in Japan. Despite geographic variation in ΔNPP, the mean values of ΔNPP under RCP2.6 and RCP8.5 were positive for all but one GCM, and the overall vulnerability to decreasing NPP was predicted to be relatively low. These results may help with recommending the best management options for planted forests in other regions with warm and humid climates. However, the consequences of decreasing NPP in the warm season were not fully examined in this study, such as the possible linkage between reduced non-structural C and mortality [28]. Monitoring of planted forests should be enhanced in southwestern Japan. The part of this region with high MAT under the current climate may require the development of new varieties of Japanese cedar or a shift toward other tree species such as Japanese cypress to properly address the vulnerability of planted forests to warming summers in the future.

Supporting information

S1 File.

https://doi.org/10.1371/journal.pone.0247165.s001

(PDF)

S2 File.

https://doi.org/10.1371/journal.pone.0247165.s002

(ZIP)

Acknowledgments

Biome-BGC v4.2 and MT-CLIM v4.3 were developed by the Numerical Terradynamic Simulation Group, University of Montana. The maps were created using the Generic Mapping Tools version 5 (http://gmt.soest.hawaii.edu/). The flux data at KHW were provided by the FFPRI FluxNet Database. The climate scenario at 1-km resolution was created and provided by the Institute for Agro-Environmental Sciences, National Agriculture and Food Research Organization. Use of the National Forest Inventory (the third term, ver. 1.0) and National Forest Resources Database (compiled in the third mesh order) was permitted by the Japanese Forestry Agency. Soil thickness data were extracted from the National Forest Soil Carbon Inventory, managed by the Japanese Forestry Agency, and legacy data accumulated by soil survey staff of the local government in the 1950s–1980s. We are grateful to Drs. Yu Morino and Kazuya Nishina of the National Institute for Environmental Studies for the use of 15-km grid nitrogen deposition data, Drs. Naoki Kabeya and Akira Shimizu of the Kyushu Research Center, FFPRI, for providing site information on the KHW, Dr. Kyotaro Noguchi of the Tohoku Research Center, FFPRI, for information on fine roots of Japanese cedar, Dr. Kazuhito Morisada of FFPRI for the use of the soil profile database, and all those who helped to develop the database of soil water retention curves. We also thank the following individuals for providing fruitful discussion of the project: Drs. Satoshi Saito, Masatake G Araki, Yuta Inoue, Tanaka Kenzo, Hiroyuki Tobita, Tadashi Sakata, Koji Shichi, Shin-ichi Iida, and Yasumasa Hirata of FFPRI, Prof. Koh Yasue of Shinshu University, Prof. Satoshi Ito of the University of Miyazaki, Prof. Yasushi Morikawa of Waseda University, Prof. Yutaka Maruyama of Nihon University, and Prof. Raisa Mäkipää and colleagues of the Natural Resources Institute, Finland.

Citation: Toriyama J, Hashimoto S, Osone Y, Yamashita N, Tsurita T, Shimizu T, et al. (2021) Estimating spatial variation in the effects of climate change on the net primary production of Japanese cedar plantations based on modeled carbon dynamics. PLoS ONE 16(2): e0247165. https://doi.org/10.1371/journal.pone.0247165

References

1. FAO. Global Forest Resources Assessment 2020: Main report. Rome. [accessed on 14th Oct., 2020]. Available from: https://doi.org/10.4060/ca9825en.

2. Van Oijen M, Balkovic J, Beer C, Cameron D, Ciais P, Cramer W, et al. Impact of droughts on the carbon cycle in European vegetation: a probabilistic risk analysis using six vegetation models. Biogeosciences. 2014;11(22): 6357–6375.

3. Bellassen V, Le Maire G, Dhote JF, Ciais P, Viovy N. Modelling forest management within a global vegetation model—Part 1: Model structure and general behaviour. Ecol. Mod. 2010;221(20): 2458–2474.

4. Jönsson AM, Lagergren F, Smith B. Forest management facing climate change-an ecosystem model analysis of adaptation strategies. Mitig. Adapt. Strat. Glob. Chan. 2015;20(2): 201–220.

5. Mitsuda Y. Evaluating the effects of climate change on the potential site productivity of Sugi (Cryptomeria japonica) planted forests in Kyushu Island, Japan. J. For. Plan. 2018;22(2): 47–53.

6. Pietsch SA, Hasenauer H, Thornton PE. BGC-model parameters for tree species growing in central European forests. For. Ecol. Manag. 2005;211: 264–295.

7. Cienciala E, Tatarinov FA. Application of BIOME-BGC model to managed forests. 2. Comparison with long-term observations of stand production for major tree species. For. Ecol. Manag. 2006;237: 252–266.

8. Ueyama M, Ichii K, Hirata R, Takagi K, Asanuma J, Machimura T, et al. Simulating carbon and water cycles of larch forests in East Asia by the BIOME-BGC model with AsiaFlux data. Biogeosciences. 2010;7(3): 959–977.

9. Ueyama M, Kai A, Ichii K, Hamotani K, Kosugi Y, Monji N. The sensitivity of carbon sequestration to harvesting and climate conditions in a temperate cypress forest: observations and modeling. Ecol. Mod. 2011;222(17): 3216–3225.

10. Ito A. The regional carbon budget of East Asia simulated with a terrestrial ecosystem model and validated using AsiaFlux data. Agric. For. Met. 2008;148(5): 738–747.

11. Kattge J, Diaz S, Lavorel S, Prentice IC, Leadley P, Bönisch G, et al. TRY–a global database of plant traits. Glob. Chan. Biol. 2011;17(9): 2905–2935.

12. Kattge J, Bönisch G, Díaz S, Lavorel S, Prentice IC, Leadley P, et al. TRY plant trait database–enhanced coverage and open access. Glob. Chan. Biol. 2020;26(1): 119–188. pmid:31891233

13. Hickler T, Vohland K, Feehan J, Miller P A, Smith B, Costa L, et al. Projecting the future distribution of European potential natural vegetation zones with a generalized, tree species-based dynamic vegetation model. Glob. Ecol. Biogeo. 2012;21(1): 50–63.

14. Anderegg WR, Konings AG, Trugman AT, Yu K, Bowling DR, Gabbitas R, et al. Hydraulic diversity of forests regulates ecosystem resilience during drought. Nature. 2018;561(7724): 538–541. pmid:30232452

15. Osone Y, Hashimoto S, Kenzo T, Araki MG, Inoue Y, Shichi K, et al. Plant trait database for Cryptomeria japonica and Chamaecyparis obtusa (SugiHinoki DB): Their physiology, morphology, anatomy and biochemistry. Eco. Res. 2020;35(1): 274–275.

16. Van Oijen M, Rougier J, Smith R. Bayesian calibration of process-based forest models: bridging the gap between models and data. Tree Physiol. 2005;25(7): 915–927. pmid:15870058

17. Minunno F, Peltoniemi M, Härkönen S, Kalliokoski T, Makinen H, Mäkelä A. Bayesian calibration of a carbon balance model PREBAS using data from permanent growth experiments and national forest inventory. For. Ecol. Manag. 2019;440: 208–257.

18. Nishizono T, Kitahara F, Iehara T, Mitsuda Y. Geographical variation in age–height relationships for dominant trees in Japanese cedar (Cryptomeria japonica D. Don) forests in Japan. J. For. Res. 2014;19(3): 305–316.

19. Fang J, Oikawa T, Kato T, Mo W, Wang Z. Biomass carbon accumulation by Japan’s forests from 1947 to 1995. Glob. Bio. Cycles. 2005;19(2): GB2004.

20. Forestry Agency. Statistical Handbook of Forest and Forestry (in Japanese). Ministry of Agriculture, Forestry and Fisheries. 2019. Available from: https://www.rinya.maff.go.jp/j/kikaku/toukei/youran_mokuzi2019.html

21. Matsumoto Y, Maruyama Y, Morikawa Y. Some aspects of water relations on large Cryptmeria japonica D. Don trees and climatic changes on the Kanto plains Japan in relation to forest decline. Jpn. J. For. Env. 1992;34: 2–13. (in Japanese with English summary)

22. Matsumoto Y, Shigenaga H, Miura S, Nagakura J, Taoda H. Mapping of Japanese cedar (Cryptomeria japonica) forests vulnerable to global warming in Japan. Glob. Env. Res. 2006;10(2): 181–188.

23. Shigenaga H, Matsumoto Y, Taoda H, Takahashi M. The potential effect of climate change on the transpiration of Sugi (Cryptomeria japonica D. Don) plantations in Japan. J. Agric. Met. 2005;60(5): 451–456.

24. Farquhar GD, von Caemmerer SV, Berry JA. A biochemical model of photosynthetic CO2 assimilation in leaves of C3 species. Planta. 1980;149(1): 78–90. pmid:24306196

25. Collatz GJ, Ball JT, Grivet C, Berry JA. Physiological and environmental regulation of stomatal conductance, photosynthesis and transpiration: a model that includes a laminar boundary layer. Agric. For. Met. 1991;54(2–4): 107–136.

26. Harley PC, Thomas RB, Reynolds JF, Strain BR. Modelling photosynthesis of cotton grown in elevated CO2. Plant Cell Env. 1992;15(3): 271–282.

27. Nishimori M, Ishigooka Y, Kuwagata T, Takimoto T, Endo N. SI-CAT 1km-grid square regional climate projection scenario dataset for agricultural use (NARO2017). J. Japan Soc. Sim. Tech. 2019;38: 150–154. (in Japanese)

28. Anderegg WR, Hicke JA, Fisher RA, Allen CD, Aukema J, Bentz B. et al. Tree mortality from drought, insects, and their interactions in a changing climate. New Phytol. 2015;208(3): 674–683. pmid:26058406

29. Thornton PE, Law BE, Gholz HL, Clark KL, Falge E, Ellsworth DS. et al. Modeling and measuring the effects of disturbance history and climate on carbon and water budgets in evergreen needle forests. Agric. For. Met. 2002;113: 185–222.

30. Jarvis PG. The interpretation of the variations in leaf water potential and stomatal conductance found in canopies in the field. Phill. Trans. R. Soc. Lond. B. 1976;273: 593–610.

31. Tatarinov FA, Cienciala E. Application of BIOME-BGC model to managed forests: 1. Sensitivity analysis. For. Ecol. Manag. 2006;237(1–3): 267–279.

32. van Genuchten MT. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 1980;44(5): 892–898.

33. Saitoh TM, Tamagawa I, Muraoka H, Lee NY, Yashiro Y, Koizumi H. Carbon dioxide exchange in a cool-temperate evergreen coniferous forest over complex topography in Japan during two years with contrasting climates. J. Plant Res. 2010;123: 473–483. pmid:20157756

34. Forestry and Forest Products Research Institute. FFPRI FluxNet Database. 2014. [data published on 31st Jul., 2014]. Available from: https://www2.ffpri.affrc.go.jp/labs/flux/index.html

35. Shimizu T, Kumagai T, Kobayashi M, Tamai K, Iida S, Kabeya N, et al. Estimation of annual forest evapotranspiration from a coniferous plantation watershed in Japan (2): comparison of eddy covariance, water budget and sap-flow plus interception loss. J. Hydrol. 2015;522: 250–264.

36. Iida S, Ohta T, Matsumoto K, Nakai T, Kuwada T, Kononov AV, et al. Evapotranspiration from understory vegetation in an eastern Siberian boreal larch forest. Agric. For. Met. 2009;149: 1129–1139.

37. Ishizuka S, Sakata T, Sawata S, Ikeda S, Takenaka C, Tamai N, et al. High potential for increase in CO2 flux from forest soil surface due to global warming in cooler areas of Japan. Ann. For. Sci. 2006;63(5): 537–546.

38. Agency Forestry. Trend in the number of forest fire in Japan since 1947. (in Japanese) [accessed on 2nd Oct., 2020]. Available from: https://www.rinya.maff.go.jp/j/hogo/yamakaji/con_1.htm

39. Noguchi K, Konôpka B, Satomura T, Kaneko S, Takahashi M. Biomass and production of fine roots in Japanese forests. J. For. Res. 2007;12(2): 83–95.

40. White MA, Thornton PE, Running SW, Nemani RR. Parameterization and sensitivity analysis of the BIOME-BGC terrestrial ecosystem model: net primary production controls. Earth Interact. 2000;4: 1–85.

41. Wang W, Ichii K, Hashimoto H, Michaelis AR, Thornton PE, Law BE, et al. A hierarchical analysis of terrestrial ecosystem model Biome-BGC: Equilibrium analysis and model calibration. Ecol. Mod. 2009;220: 2009–2023.

42. Komatsu H, Shinohara Y, Kumagai T, Kume T, Tsuruta K, Xiang Y. et al. A model relating transpiration for Japanese cedar and cypress plantations with stand structure. For. Ecol. Manag. 2014;334: 301–312.

43. Matsumoto M, Awaya Y, Iehara T, Takahashi M, Fujiwara T, Hosoda K. et al. Development of National Forest Resources Database for accounting and reporting under the Kyoto Protocol. FORMATH. 2006;6: 141–163. (in Japanese)

44. Ugawa S, Takahashi M, Morisada K, Takeuchi M, Matsuura Y, Yoshinaga S. et al. Carbon stocks of dead wood, litter, and soil in the forest sector of Japan: general description of the National Forest Soil Carbon Inventory. Bull. FFPRI. 2012;425: 207–221.

45. Morisada K, Ono K, Kanomata H. Organic carbon stock in forest soils in Japan. Geoderma. 2004;119(1–2): 21–32.

46. Forest Soil Division. Classification of forest soils in Japan (1975). Bull. Gov. For. Exp. Stat. 1976;280: 1–28. (in Japanese with English summary)

47. National Land Agency. Digital National Land Information: land classification survey. (in Japanese) 1967–1978. Available from: https://nlftp.mlit.go.jp/kokjo/inspect/landclassification/land/l_national_map_20-1.html.

48. Zhu J, Mohanty BP. Spatial averaging of van Genuchten hydraulic parameters for steady-state flow in heterogeneous soils: A numerical study. Vadose Zone J. 2002;1(2): 261–272.

49. Morino Y, Ohara T, Kurokawa J, Kuribayashi M, Uno I, Hara H. Temporal variations of nitrogen wet deposition across Japan from 1989 to 2008. J. Geophys. Res. 2011;116: D06307.

50. Nishina K, Watanabe M, Koshikawa MK, Takamatsu T, Morino Y, Nagashima T, et al. Varying sensitivity of mountainous streamwater base-flow NO3− concentrations to N deposition in the northern suburbs of Tokyo. Scientific Reports. 2017;7(1): 1–9. pmid:28127051

51. Nioh I. Nitrogen fixation associated with the leaf litter of Japanese cedar (Cryptomeria japonica) of various decomposition stages. Soil Sci. Plant Nut. 1980;26(1): 117–126.

52. Yamanaka T, Hirai K, Aizawa S, Yoshinaga S, Takahashi M. Nitrogen-fixing activity in decomposing litter of three tree species at a watershed in eastern Japan. J. For. Res. 2011;16(1): 1–7.

53. Menge DN, Hedin LO. Nitrogen fixation in different biogeochemical niches along a 120 000-year chronosequence in New Zealand. Ecology. 2009;90(8): 2190–2201. pmid:19739381

54. IPCC. Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, Stocker TF, et al. (Eds.). Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA, 1535 pp, 2013. https://doi.org/10.1017/CBO9781107415324

55. Riahi K, Gruebler A, Nakicenovic N. Scenarios of long-term socio-economic and environmental development under climate stabilization. Tech. Forecas. Soc. Chan. 2007;74(7): 887–935.

56. van Vuuren D, den Elzen M, Lucas P, Eickhout B, Strengers B, van Ruijven B. et al. Stabilizing greenhouse gas concentrations at low levels: an assessment of reduction strategies and costs. Clim. Chan. 2007;81(2): 119–159.

57. Kimball JS, Running SW, Nemani R. An improved method for estimating surface humidity from daily minimum temperature. Agric. For. Met. 1997;85: 87–98.

58. Thornton PE, Running SW. An improved algorithm for estimating incident daily solar radiation from measurements of temperature, humidity, and precipitation. Agric. For. Met. 1999;93: 211–228.

59. Ando T, Hatiya K, Doi K, Kataoka H, Kato Y, Sakaguchi K. Studies on the system of density control of Sugi (Cryptomeria japonica) stand. Bull. Gov. For. Exp. Stat. 1968; 209: 1–76. (in Japanese with English summary)

60. Greenhouse Gas Inventory Office of Japan. National greenhouse gas inventory report of Japan. 2008. [accessed on 14th Oct., 2020] Available from: https://www.env.go.jp/earth/ondanka/ghg-mrv/unfccc/material/NIR-JPN-2008_E.pdf

61. Yashiro Y, Lee NYM, Ohtsuka T, Shizu Y, Saitoh TM, Koizumi H. Biometric-based estimation of net ecosystem production in a mature Japanese cedar (Cryptomeria japonica) plantation beneath a flux tower. J. Plant Res. 2010;123(4): 463–472. pmid:20339895

62. Nagakura J, Shigenaga H, Akama A, Takahashi M. Growth and transpiration of Japanese cedar (Cryptomeria japonica) and Hinoki cypress (Chamaecyparis obtusa) seedlings in response to soil water content. Tree Physiol. 2004;24(11): 1203–1208. pmid:15339729

63. Kumagai T, Tateishi M, Shimizu T, Otsuki K. Transpiration and canopy conductance at two slope positions in a Japanese cedar forest watershed. Agric. For. Met. 2008;148(10): 1444–1455.

64. Saitoh TM, Tamagawa I, Muraoka H, Kondo H. An analysis of summer evapotranspiration based on multi-year observations including extreme climatic conditions over a cool-temperate evergreen coniferous forest, Takayama, Japan. Hydro. Process. 2013;27(23): 3341–3349.

65. Saito T, Kumagai T, Tateishi M, Kobayashi N, Otsuki K, Giambelluca TW. Differences in seasonality and temperature dependency of stand transpiration and canopy conductance between Japanese cypress (Hinoki) and Japanese cedar (Sugi) in a plantation. Hydro. Process. 2017;31(10): 1952–1965.

66. Hashimoto S, Wattenbach M, Smith P. Litter carbon inputs to the mineral soil of Japanese Brown forest soils: comparing estimates from the RothC model with estimates from MODIS. J. For. Res. 2011;16(1): 16–25.

67. Ding Y, Liang S, Peng S. Climate change affects forest productivity in a typical climate transition region of China. Sustainability. 2019;11(10): 2856.

68. Sun XY, Wang GX, Huang M, Hu ZY, Song CL. Effect of climate change on seasonal water use efficiency in subalpine Abies fabri. J. Mount. Sci. 2017;14(1): 142–157.

69. Lei X, Yu L, Hong L. Climate-sensitive integrated stand growth model (CS-ISGM) of Changbai larch (Larix olgensis) plantations. For. Ecol. Manag. 2016;376: 265–275.

70. Ueyama M, Ichii K, Kobayashi H, Kumagai T, Beringer J, Merbold L, et al. Inferring CO2 fertilization effect based on global monitoring land-atmosphere exchange with a theoretical model. Env. Res. Let. 2020;15(8): 084009.

71. Hiraoka Y, Iki T, Nose M, Tobita H, Yazaki K, Watanabe A, et al. Species characteristics and intraspecific variation in growth and photosynthesis of Cryptomeria japonica under elevated O3 and CO2. Tree Physiol. 2017;37(6): 733–743. pmid:28369644

72. Franks PJ, Adams MA, Amthor JS, Barbour MM, Berry JA, Ellsworth DS, et al. Sensitivity of plants to changing atmospheric CO2 concentration: from the geological past to the next century. New Phytol. 2013;197: 1077–1094. pmid:23346950

73. Hidy D, Barcza Z, Dobor L, Gelybó G, Fodor N, Pintér K. et al. Terrestrial ecosystem process model Biome-BGCMuSo v4. 0: summary of improvements and new modeling possibilities. Geosci. Mod. Dev. 2016;9: 4405–4437.

74. Van Oijen M, Beer C, Cramer W, Rammig A, Reichstein M, Rolinski S, et al. A novel probabilistic risk analysis to determine the vulnerability of ecosystems to extreme climatic events. Env. Res. Let. 2013;8(1): 015032.

75. Reeves MC, Moreno AL, Bagne KE, Running SW. Estimating climate change effects on net primary production of rangelands in the United States. Clim. Chan. 2014;126(3–4): 429–442.

76. Yin Y, Ma D, Wu S. Climate change risk to forests in China associated with warming. Scientific Reports. 2018;8(1): 1–13. pmid:29311619