A land surface model combined with a crop growth model for paddy rice ( MATCRO-Rice Ver . 1 ) – Part II : Model validation

A land surface model combined with a crop growth model for paddy rice (MATCRO-Rice Ver. 1) – Part II: Model validation Yuji Masutomi1, Keisuke Ono2, Takahiro Takimoto3, Masayoshi Mano4, Atsushi Maruyama2, and Akira Miyata2 1College of Agriculture, Ibaraki University, 3-21-1, Chuo, Ami, Inashiki, Ibaraki 300-0393, Japan 2Intstitute for Agro-Environmental Sciences, NARO, 3-1-3, Kannondai, Tsukuba, Ibaraki 305-8604, Japan 3Institute for Global Change Adaptation Science, Ibaraki University, 3-21-1, Chuo, Ami, Inashiki, Ibaraki 300-0393, Japan 4Graduate School of Horticulture, Chiba University, 648 Matsudo, Matsudo-shi, Chiba 271-8510, Japan Correspondence to: Yuji Masutomi (yuji.masutomi@gmail.com)


Introduction
It has been recognized that crop growth and management in agricultural land are important factors that affect climate at various spatial and temporal scales via exchange of heat, water, and gases (Tsvetsinskaya et al., 2001;Bondeau et al., 2007;Osborne et al., 2009;Levis et al., 2012).Betts (2005) pointed out that integration of crop growth models (CGMs) into climate models is needed for accurate climate simulations by climate models.To consider the influence of agricultural land on climate in climate simulations, several land surface models (LSMs) or dynamic vegetation models (DVMs) incorporated with a CGM have been developed (Tsvetsinskaya et al., 2001;Kucharik, 2003;Gervois et al., 2004;Bondeau et al., 2007;Osborne et al., 2007;Lokupitiya et al., 2009;Maruyama and Kuwagata, 2010;Levis et al., 2012;Osborne et al., 2015).Masutomi et al. (2016) have developed a new LSM-CGM combined model, called MATCRO-Rice, by incorporating a CGM into an LSM, MATSIRO (Takata et al., 2003).The most important feature of the model is that it can consistently simulate latent heat flux (LHF), sensible heat flux (SHF), net carbon uptake by crop, and crop yield in paddy rice fields by exchanging variables between an LSM and a CGM.The consistency among model outputs enable us to apply the model to a wide range of integrated issues.For example, the model can investigate the interaction between climate and paddy rice fields, consistently considering impacts of climate on rice productivity and impacts of paddy rice fields on climate.Osborne et al. (2009) showed that this interaction can affect variability in climate and crop production.Therefore, -degree latitude of the simulation site W glu,0 0.5 kg ha −1 dry weight of glucose reserve at emergence W lef,0 1.0 kg ha −1 dry weight of leaf at emergence W rot,0 1.0 kg ha −1 dry weight of root at emergence W stm,0 1.0 kg ha −1 dry weight of stem at emergence z a 3 m reference height at which wind speed is observed z max 4 m depth of soil layer z t 0.05 m depth of top soil layer z b 2 m depth from the soil surface to the upper bound of the bottommost layer of soil δt 1800 s time resolution the understanding of the interaction is important for securing food security.However, little is known about the interaction.MATCRO-Rice can be a useful tool for studying the interaction between climate and paddy rice fields.
The objective of the present paper is to present the results of the comprehensive validation of MATCRO-Rice and to show the effects of modifications from the original LSM, MATSIRO.Before presenting the results of the validation and the effects of modification, we first show the numerical method (Sect.2) and the method and results of parameterization for model parameters (Sect. 3).The results of model validation and the effects of modifications are shown in Sects.4 and 5, respectively, followed by concluding remarks in Sect.6.

Numerical setting and method
All simulation setting parameters are shown in Table 1.We set the time resolution of the simulation to half an hour, i.e. δ t = 1800 s.For time discretization, the forward difference method was used.
To simulate soil water and heat transfer (Sect.3.5 in Masutomi et al., 2016), we spatially discretized soil into five layers with thickness of 0.05, 0.2, 0.75, 1.0, and 2.0 m, resulting in z max = 4.0 m, z t = 0.05 m, and z b = 2.0 m.To simulate soil water content for each soil layer (w s ), we replaced the gradient of water flux by net water fluxes between layers.In the calculation for water fluxes between layers, we used the hydraulic conductivity that is smaller among soil layers and the difference in water potentials between soil layers.After the calculation for soil water content for each layer, water content beyond saturation was taken out to base flow.
To simulate soil temperature for each soil layer, we solved the system of equations for soil layers by using the Gauss-Jordan method.In the calculation of soil temperatures, we replaced the gradient of heat flux by net heat fluxes between layers.In the calculation of heat fluxes between layers, we used thermal conductivity averaged between soil layers and soil temperatures for each layer.
The downhill simplex method (Nelder and Mead, 1965) was used to simulate temperatures of the canopy and surface (T c and T g ; Sect.3.1 in Masutomi et al., 2016), bulk transfer coefficients (C E g , C E c , C H g , C H c , C M , and C Mg ; Sect.3.3 in Masutomi et al., 2016), and variables related to carbon assimilation (A n,x , c i,x , and g st,x ; Sect.4.1 in Masutomi et al., 2016).
We set z a = 3 m.CO 2 concentration (C a,ppm ) and the depth of surface water (d w ) were set at 390 ppm and 0.025 m, respectively.The initial dry weight of each organ was set at 1 kg ha −1 for leaf (W lef,0 ), stem (W stm,0 ), and root (W rot,0 ), and at 0.5 kg ha −1 for glucose reserve in the leaf (W glu,0 ).D oy,Ie , D oy,Is , D oy,sw , and L t depend on the simulations.Values for these parameters are shown in the sections of each simulation.

Parameterization
Table 2 shows model parameters parameterized in the present paper.All parameters are parameterized using observations, the literature, and assumptions.The method of the parameterization is explained in this section.

Parameterization site and observation data
Table 3 shows the observational data used for parameterization.The data were observed from 2003 to 2006 at a site which is located in Tsukuba, Japan (lat: 36 • 03 14.3 N; long: 140 • 01 36.9E), at 13 m above sea level.The climatic zone of the site is temperate, with the mean annual air tempera-Geosci.Model Dev., 9, [4155][4156][4157][4158][4159][4160][4161][4162][4163][4164][4165][4166][4167]2016 www.geosci-model-dev.net/9/4155/2016/Biomass for each organ (W lef , W pnc , W rot , and W stm ) and leaf area index (L) were measured nearly every 2 weeks.At each measuring time, 10 stands were sampled from the fields.Yield (Y ld ) and phenological dates including trans-planting (D oy,tr ), heading (D oy,hd ), and harvest (D oy,hv ) were observed every year.The values of observed yield are the husked rice yield with 15 % water content.The rice grains for measuring yield were sampled from the whole fields of the observational site.The crop height (h gt ) was measured on average every 5 days.

Phenology
Phenological parameters that represent development stages (D vs,e , D vs,h , G ds,m , D vs,tr , and D vs,te ) were parameterized.First, we calculated D vs values at heading and G ds,m values from 2003 to 2006 using the phenological model given by Masutomi et al. (2016).The mean values were set to D vs,h and G ds,m , resulting in D vs,h = 0.616 and G ds,m = 167759940 K s. Figure 1 compares the heading and harvest dates between observations and simulations from 2003 and those from 2006.The simulated heading and harvest dates were in good agreement with the observations.The average errors were 2.25 and 4.5 days for heading and harvest, respectively.
D vs,e , D vs,tr , and D vs,te were determined so that the duration from sowing to emergence, transplanting, and the end of transplanting shock was 5, 20, and 25 days, respectively.Thus, D vs,e = 0.012, D vs,tr = 0.06, and D vs,te = 0.08.

Partitioning
MATCRO partitions carbohydrates in leaves, in the form of glucose, into each organ, according to MACROS (Penning de Vries et al., 1989).Parameters related to glucose partitioning (D vs,rot1 , D vs,rot2 , D vs,lef1 , D vs,lef2 , D vs,pnc1 , D vs,pnc2 , P rot , and P lef ) were parameterized as follows: (i) we calculated the ratio of glucose partitioned to each organ (leaf, stem, root, panicle) during the growing period using the observed biomass for each organ; (ii) we conducted the curve fitting of the calculated ratios in (i). Figure 2 shows the calculated ratios of glucose partitioned to each organ and the fitting curves for the ratios.
To determine the ratio of dead leaves at harvest (r d1,lef ), we first calculated the observational ratios of dead leaves during growing periods by dividing the decrease in leaf biomass between observational dates by the duration among the observational dates.Then, by graphically fitting a curve to the calculated ratios of dead leaves, we determined r d1,lef .Figure 3 shows the calculated ratios of dead leaves and the fitted curve.
The fraction of glucose allocated to starch reserve (f stc ) is determined as follows: (i) we first calculated the ratios of stem biomass at harvest to maximum stem biomass for each year from 2003 to 2006 (Bouman et al., 2001); (ii) then, a 4-year average was calculated for f stc .

LAI, crop height, and specific leaf weight
To obtain the parameters for the relationship between LAI and crop height (h aa , h ab , h ba , and h bb ), we conducted linear regressions of the data before and after heading using observations for LAI and crop height from 2003 to 2006.Thus, h aa = 0.439, h ab = 0.675, h ba = 0.366, and h bb = 0.318.Figure 4 compares the LAI-height relation between observations and simulations.
To obtain parameters for specific leaf weight (k S lw , S lw,mn , and S lw,mx ), we plotted observations for specific leaf weights during growing periods from 2003 to 2006 and conducted the curve fitting of the plotted data.Thus, k S lw = 3.5, S lw,mn = 350, and S lw,mx = 600.Figure 5 shows the specific leaf weights and the fitted curve.

Crop yield
To determine the ratio of crop yield to dry weight of panicle at harvest (k yld ), we calculated the dry weight of panicle at harvest, because the weight was not observed.By assuming linear increase of dry weight from the last date in which dry weight of the panicle was measured, we calculated the dry weight of the panicle at harvest from 2003 to 2006.The median value among the ratios of observed yields to the calculated dry weight of panicle produced k yld .

RuBisCO-limited photosynthesis rate
Parameters related to RuBisCO-limited photosynthesis rate (V max (0), s 1 , and s 2 ) were parameterized using the values obtained from the literature.In this parameterization, we adjusted the parameters so that the RuBisCO-limited photosynthesis rate (ω c,x ) simulated by MATCRO agrees with the observational value reported by Borjigidai et al. (2006).In the simulations, CO 2 concentration in the leaf was fixed to c i,x = 30 Pa. Figure 6, showing the comparison of RuBisCOlimited photosynthesis rates among MATCRO, those reported by Borjigidai et al. (2006), and MATSIRO (Takata et al., 2003), on which MATCRO is based, indicates that there is a good agreement in the photosynthesis rate between the simulations of MATCRO and the observational value in Borjigidai et al. (2006); the simulations for the photosynthesis rate of MATCRO were significantly improved compared to those of MATSIRO.

Validation
We conducted two types of validation.The first validation was conducted at the parameterization site as explained in Sect.3.1.In the validation, the simulated LHF, SHF, carbon uptake by crop, and crop yields were compared with the observations from 2003 to 2006.The second validation was conducted for a territory across Japan.The simulated crop yields for Japan were compared with the national statistics from 1991 to 2010.

Input and validation data
Table 4 shows the observational data used for the validation.Information on the instruments used for the observations is available from the AsiaFlux website (AsiaFlux, 2016).The height at which the wind speed was measured was different each year.Assuming the logarithmic vertical profile of the wind, we transformed the observed wind speed to that at 3 m       7).These RMSE values are comparable to those reported in earlier studies (Kimura and Kondo, 1998;Maruyama and Kuwagata, 2010).One of the major reasons for the errors of LHF and SHF between simulations and observations is thought to be a problem in flux observations.Aubinet et al. (1999) reported that the energy balance in observations is not closed.In contrast, the energy balance simulated by MATCRO is completely   closed.Therefore, the energy imbalance in flux observations can cause errors between simulations and observations.El Maayar et al. (2008) suggested to test the degree of energy imbalance in observations before comparing the observations with simulations.This degree is generally evaluated by where H , λE, R n , and G are the daily averages for SHF, LHF, net radiation, and heat flux into ground, respectively; d indicates a day; and N is the number of days.The observation values for R n and G in this equation are expected to be sufficiently accurate (Twine et al., 2000;Wilson et al., 2002).The values of I m in the observations from 2003 to 2006 were 0.79, 0.77, 0.78, and 0.74, with the average of 0.78.In other words, these results imply that the total flux of observed LHF and SHF can be smaller than a true value.The ratios of the total flux of observed LHF and SHF to those of simulated LHF and SHF from 2003 to 2006 were 0.84, 0.79, 0.80, and 0.83, with the average of 0.82.This suggests that the errors of LHF and SHF between observations and simulations can be largely attributed to the energy imbalance in observations.

Comparison of net carbon uptake by crop
In this section, we tested the accuracy of MATCRO for simulating net carbon uptake by crops during growing periods by comparing the changes in total biomass between simulations and observations.Figure 11    Hence, we conclude that the model has high accuracy for simulating net carbon uptake by crops during growing periods.

Method, input, and validation data
Rice yields for all prefectures in Japan were simulated from 1991 to 2010, and then the average national rice yield was calculated for each year from the prefectural yields weighted by the prefectural planting areas.The simulated rice yields for Japan were compared with the national crop statistics (MAFF, 1991(MAFF, -2010)).
The Global Meteorological Forcing Dataset (GMFD) for land surface modelling (Sheffield et al., 2006) was used for meteorological input data.Because the spatial resolution of the GMFD is 1 • , the average meteorological values for each prefecture were calculated from the gridded values weighted by the prefectural planting areas.The time resolution of the GMFD is 3 h.We used the same values every 3 h for halfhourly simulations.The comparison with the agrometeorological dataset in Japan (Ohno et al., 2012) revealed that shortwave radiation of the GMFD has a large error; we corrected the error using the agrometeorological dataset, so that the daily values for shortwave radiation of the GMFD could agree with those of the agrometeorological dataset.In addition, we changed the ratio of photosynthetic active radiation to shortwave radiation from 0.5 (default) to 0.423, according to the value observed at the parameterization site (Sect.3.1).
The sowing and harvesting dates (D oy,sw and D oy,hv ) for each prefecture were obtained from the national crop statistics (MAFF, 1991(MAFF, -2010)).In this validation, simulated harvesting dates were not used because the results in the previous section (Sect.4.1.4)showed that simulated harvesting dates may cause errors in yield simulations.
For simplicity, land surface was assumed to be annually flooded and irrigated.Hence, we set D oy,Is = 1 and D oy,Ie = 365.L t for each prefecture was set to the latitude of the centre of each prefecture.The other simulation settings and parameters for crops and soil were set using the values shown in Tables 1, 2, and 5.

Comparison of yield
Figure 13 shows the comparison of the observed and simulated yields from 1991 to 2010.The average error between simulations and observations was 16.7 %.Although the simulated yields tended to overestimate the observations, the correlation between simulations and observations was significant at 0.663 (P < 0.01).Hence, we conclude that MATCRO can reproduce substantial parts of weatherinduced variability in yields.One of the reasons for the simulation errors is thought to be the method of parameterization.For the simulations over Japan, we used parameters parameterized for a variety Koshihikari only at one site in Japan (Sect.3), although there is a large diversity in agricultural management and techniques, and rice varieties planted throughout Japan.Therefore, parameterization on a large scale is necessary for better large-scale simulations.

Effects of modifications
There are two major modifications of MATCRO from the original LSM (MATSIRO).The first one is the dynamic calculation of LAI, crop height, and root length.The other is the consideration of flooded surface and irrigation.We quantify the effects of the two major modifications on the simulation www.geosci-model-dev.net/9/4155/2016/Geosci.Model Dev., 9, 4155-4167, 2016 of LHF and SHF.Both simulations are conducted at the parameterization site from 2003 to 2006 (Sect.3.1).

Effect of dynamic calculation of LAI
The original LSM (MATSIRO) uses the monthly constant LAI, which is given in gridded form.The default gridded LAI data of MATSIRO were obtained from the Global Soil Wetness Project 2 (Dirmeyer et al., 2006).Figure 14 20.78, 19.28, and 18.72 W m −2 , respectively, with the 4year average of 19.54 W m −2 (Table 6).The RMSEs of daily SHF from 2003 to 2006 were 14.69, 20.30, 16.93, and20.68W m −2 , respectively, with the 4-year average of 18.15 W m −2 (Table 7).These errors are comparable to those of MATCRO.Hence, we conclude that the effects of the dynamic calculation of LAI, crop height, and root length on LHF and SHF are small.

Effect of flooded and irrigated surface
We simulated LHF and SHF from 2003 to 2006 without flooded surface and irrigation.The simulations are called MATCRO-RF, which denotes simulations under rain-fed conditions.Figures 15 and 16 show the comparison of daily LHF and LHF between observations and simulations obtained by MATCRO and MATCRO-RF.The LHF and SHF simulated by MATCRO-RF were not in agreement with the observations.The RMSEs of daily LHF simulated by MATCRO-RF from 2003 to 2006 were 16.63, 36.90, 29.32, and24.93W m −2 , respectively, with the 4-year average of 26.95 W m −2 (Table 6).The RMSEs of daily SHF from 2003 to 2006 were 16.34, 42.02, 34.16, and31.56W m −2 , respectively, with the 4-year average of 31.02W m −2 (Table 7).These errors in MATCRO-RF are considerably larger than those simulated by MATCRO.Figures 15 and 16 show that MATCRO-RF tends to underestimate daily LHF and to overestimate daily SHF.The underestimation of daily LHF can be attributed to lower evaporation from the soil due to the absence of irrigation and flooding.The overestimation of daily SHF can be caused by high surface temperature due to lower evaporation from the soil.Hence, we conclude that flooded surface and irrigation have large effects on simulations of LHF and SHF.

Concluding remarks
In this paper, we presented the results of the validation of MATCRO-Rice and the effects of the modification of the original LSM (MATSIRO), and the numeric and parameterization methods.First, the comparison of the LHF and SHF between simulations and observations at the paramerization site confirmed that the model can reproduce the observed LHF and SHF data well.The accuracy of the simulations for LHF and SHF was comparable to those obtained in earlier studies.Second, we showed that the simulated growth of the total biomass was in good agreement with the observations at the parameterization site.This indicates that the model can simulate the net carbon uptake by crops during a growing period at paddy rice fields.Last, we demonstrated that the model has high ability to simulate crop yield by comparing the simulated and observed yields at the parameterization site and over Japan.
The validation results suggest that MATCRO-Rice has high ability to accurately and consistently simulate LHF, SHF, net carbon uptake by crop, and yield.There have been many models that simulate some of the four variables with high accuracy, but a few models can accurately and consistently simulate all four of them.This point is the most important feature of MATCRO-Rice.The model can be applied to a wide range of issues, including climate change impact (e.g.Masutomi et al., 2009), and it will facilitate the scientific research especially on the climate-crop interactions (Osborne et al., 2009).We validated LHF, SHF, and carbon flux simulated by this model with observations from only one site.The model should be further validated at multiple sites in order to enforce the reliability and applicability of the model.However, since there are a few flux sites on agriculture land worldwide, it will be necessary to increase their number on agricultural land to promote climate-crop modelling studies.
We assessed the effects of the dynamic simulation of LAI, crop height, and root length on LHF and SHF and the effect of flooded surface and irrigation on LHF and SHF.The results show that the effects of the dynamic simulation on LHF and SHF are small, whereas the flooded surface and irrigation have large effects on LHF and SHF.These results suggest that climate-crop modelling should incorporate flooded surface and irrigation.

Code and data availability
The source code of MATCRO will be distributed upon request to the corresponding author (Yuji Masutomi: yuji.masutomi@gmail.com).The website for MATCRO-Rice will be developed in the near future.

Figure 1 .
Figure 1.Comparison of heading and harvest dates.SIM: simulations; OBS: observations; DOY indicates the number of days from 1 January.

Figure 2 .
Figure 2. Partitioning ratio of glucose.Red lines are fitted.D vs indicates the development stage.

Figure 3 .
Figure 3. Ratio of dead leaves.A red line is fitted.D vs indicates the development stage.

Figure 4 .
Figure 4. Relationship between the leaf area index (LAI) and crop height (black curve: h gt = h aa L h ab ; red curve: h gt = h ba L h bb ).

Figure 5 .
Figure 5. Relationship between specific leaf weight and development stage (D vs ).A red curve denotes the fitted curve used in the model.

Figure 7 .
Figure 7.Comparison of daily LHF between simulations and observations.DOY indicates the number of days from 1 January.

Figure 8 .
Figure 8.Comparison of half-hourly LHF between simulations and observations.
compares the growth of the total biomass between simulations and observations from 2003 to 2006.As indicated by the figure, the simulated total biomass was in good agreement with the observations.

Figure 9 .
Figure 9.Comparison of SHF between simulations and observations.DOY indicates the number of days from 1 January.

Figure 10 .
Figure 10.Comparison of half-hourly SHF between simulations and observations.

Figure 11 .
Figure 11.Comparison of total biomass between simulations and observations during growing periods from 2003 to 2006.Circles indicate mean values of observations and the ranges indicate standard deviation of observations.Red lines denote simulations.DOY indicates the number of days from 1 January.

Figure 12 .
Figure 12.Comparison of yields between simulations and observations.

Figure 13 .
Figure 13.Comparison of yields over Japan between simulations and observations.

Figure 14 .
Figure 14.Comparison of LAI between observations, simulations by MATCRO, and the default value of MATSIRO.

Figure 15 .
Figure 15.Comparison of daily LHF between observations and simulations by MATCRO and MATCRO-RF.DOY indicates the number of days from 1 January.

Figure 16 .
Figure 16.Comparison of SHF between observations and simulations by MATCRO and MATCRO-RF.DOY indicates the number of days from 1 January.

Table 3 .
Observational data used for parameterization.
• C and precipitation of 1200 mm.The soil type is clay loam.The variety planted at the site is Koshihikari, which is the most planted variety in Japan.

Table 4 .
Observational data used for validation at the parameterization site.
shows the comparison of LAI between observations, simulations by MATCRO, and the default values of MATSIRO.We can see that MATCRO adequately reproduces seasonal changes in LAI, although the default LAI are not in agreement with the observed LAI.MATSIRO also uses constant crop height and root length, which are vegetation-specific parameters.The default values of crop height and root length for crops are 1 m.Using the default data for LAI and the default values for crop height and root length, we simulated LHF and SHF from 2003 to 2006.In the simulations, we also used the original equation of MATSIRO for calculating the maximum canopy water (w cap ).The RMSEs of daily LHF from 2003 to 2006 were 19.4,