Modeling the multiple effects of temperature and radiation on rice quality

Ongoing climate change is likely to enhance the deterioration of rice quality that has been observed in western Japan, especially in Kyushu, since the 1990s. Therefore, it is important to examine the response of rice quality to environmental variation over a wide geographical domain. To that end, the aims of this study were (i) to propose a statistical model to predict rice quality based on temperature, total radiation during the ripening period, and their multiple effects; and (ii) to evaluate the model validity and uncertainty in prediction. A Bayesian calibration was adopted to account for uncertainty in the parameter values associated with non-climatic factors. The validation results showed that the model performed well in capturing the temporal trend and interannual variation in observed rice quality in all prefectures of Kyushu. We then performed the prediction experiment for rice quality in the extremely hot summer of the year 2010, which was omitted from the model calibration data. The results showed that the predictive capability of the statistical model is somewhat dependent on the calibration data, but this dependency does not necessarily mean that useful predictions for climates not in the calibration data are impossible.


Introduction
Declines in rice quality have been observed in western Japan, especially in Kyushu, since the 1990s (Morita 2008, Okada et al 2009. Such declines are likely to lower the eating quality of rice (Terao et al 2005, Wakamatsu et al 2007 and reduce farm income and consumer utility in Japan and other countries where the demand for high-quality rice has been increasing. The major reason for the decline in rice quality is the occurrence of chalky grains, especially milky white grains (Morita 2008). Chalky grains sharply increase when the mean daily minimum temperature for the 20 days after heading (ear emergence) exceeds 22 • C (Tsukimori 2003). The underlying mechanisms for the occurrence of chalky grains in rice plants are: reduced allocatable carbohydrates in the plant associated with an increased night-time respiration rate (Vong andMurata 1977, Hirai et al 2003); reduced capacity of stems and leaves for assimilation (Kobata et al 2004, Morita et al 2005; insufficient solar radiation during the ripening period (Matsushima and Manaka 1957); and hits of typhoons during the ripening period (Wakamatsu et al 2007).
Ongoing climate change may reduce rice quality in the near future. No studies have assessed the possible impact of climate change on rice quality, although some studies have proposed process-based models to predict rice quality based on field experimental results (Nagahata et al 2006, Nakagawa et al 2008. However, these models are designed for prediction at the field scale, and a large gap exists between the spatial scale at which these models operate and the scale at which climate projections are developed. Furthermore, it is difficult to obtain detailed information on cultivars and management practices over large areas, which is essential for a process-based model to simulate the complicated biochemical processes that govern rice quality.
Another important issue for impact assessment is uncertainty about the model's applicability to accommodate unprecedented climates, because all impact models are developed and calibrated on the basis of historical data. This corresponds to the uncertainty of future impacts associated with the extrapolation of current knowledge to future unprecedented climates. Therefore, the central objective of impact assessment model validation should be evaluation of the predictive capability of impact models under unprecedented climates.
In this study, we propose a statistical model that has a medium level of complexity to predict rice quality at broad spatial scales, that is, the model is less complex than fieldscale process-based models but more complex than simple regression models. A Bayesian calibration method ) was adopted to account for the uncertainty of non-climatic factors (e.g. cultivar and management) in model parameter values. To evaluate the model's capability and applicability for impact assessment, we conducted two types of uncertainty analyses using this model: (i) the sensitivity of the modeled rice quality to temperature increases; and (ii) prediction experiments for the extremely hot summer of 2010, the data from which were not incorporated in the calibration data. The summer of 2010 was the hottest summer in Japan since 1898, and the mean temperature anomaly for June-August in that year was 1.64 • C (Japan Meteorological Agency 2011). This resulted in the lowest recorded rice quality in western Japan since collection of comparable statistics was commenced in 1999 (MAFF 2010a). All analyses were carried out for the Kyushu region in western Japan (figure 1).

Data
According to Japan's Agricultural Products Inspection Act, harvested rice grains are categorized into four grades: first grade, second grade, third grade and irregular. The major criterion for assigning rice to lower grades is the percentage of chalky grains. Chalky grains are immature and the entire endosperm has a chalky texture, whereas refined whole grains are translucent in appearance (Tashiro and Wardlaw 1991). For rice quality data, we used the proportion of first-grade rice for seven prefectures in Kyushu for the period 1979-2007 from government statistics (MAFF 2010a). The data for 2010 were obtained from a rapid assessment released by the government (MAFF 2010a). Data on heading and harvest dates were obtained from MAFF (2010b).
Daily minimum temperatures and accumulated solar radiation for the same period were obtained from a grid dataset developed at the National Institute for Agro-Environmental Science (called Mesh-AMeDAS; Seino 1993). The grid interval of the dataset is 30 × 45 in latitude and longitude (about 1 km × 1 km). Land-use data on the same grid interval were obtained from the same dataset. To combine the weather and rice quality data, daily values of the climate variables for grid cells that contained paddy fields ( 20% of a grid cell) were spatially averaged.
For model calibration, the data on typhoon track and damage in paddy rice production (MAFF 2010b) were used to exclude rice quality data from years in which severe damage occurred in the study areas during the ripening period. The data from 1991, 1993, 2004 and 2006, just four out of the 29 years , were removed from the calibration data for Fukuoka. This treatment avoided overfitting the model, considering that typhoon damage to rice quality was not considered in the model.

Rice quality model
We preliminary examined the relationship between rice quality, temperature and total radiation during the ripening period using the rice quality data from the governmental crop statistics in all prefectures of Kyushu. The relationships between total radiation for the ripening period and rice quality at three temperature levels (represented by the mean daily minimum temperature for the 20 days after heading) are shown in figure 2. The data show that as temperature increases rice Figure 2. Relationships for observed rice quality (represented by the proportion of first-grade rice) versus total radiation during the ripening period at three temperature levels (represented by the mean daily minimum temperature during the 20 days after heading, T 20 ). Curves indicate the logistic regressions fitted to the data at each temperature level annotated with their correlation coefficients (R) and statistical significance (***, P < 0.001; **, P < 0.01; and *, P < 0.05). quality tends to decline. At temperatures <21 • C, most rice quality data have a high percentage of first-grade rice across the range of radiation. At each higher-temperature level (i.e. 21-22.9 and 23 • C), the decline in rice quality at lower radiation levels becomes increasingly pronounced, showing that the sensitivity of rice quality to insufficient radiation increases as temperature increases.
We formulated the statistical relationships between rice quality, temperature, and radiation from the logistic function: where Q is the proportion of first-grade rice (%), Q max and Q min (T ) are the upper and lower limit of Q (%), respectively, T is the mean daily minimum temperature for n days after heading ( • C), S is the total radiation during the ripening period (MJ m −2 ), f a (T ) is the sensitivity coefficient of Q to S, and f b (T ) is the value of S at which rice quality becomes halfway between the upper and lower limits (i.e. (Q max + Q min (T ))/2). The variables Q min (T ), f a (T ), and f b (T ) were assumed to be linear functions of T to account for the multiple effects of temperature and radiation on rice quality: where p i (i = 1-6) are parameters. The mean daily minimum temperature for n days after heading, T , and total radiation during the ripening period, S, were represented by: and where T min i is the daily minimum temperature on the i th day after heading, n is the period after heading in which the temperature has a negative impact on rice quality (days), S i is the daily total radiation on the i th day after heading (MJ m −2 day −1 ), and m is the period from heading to maturity (days). The variable n depends upon non-climatic factors such as the ripening ability of the cultivar, fertilization and irrigation during the ripening period, water temperature, and other factors. These sources of variation were accounted for by Bayesian calibration.

Bayesian calibration
Rice quality nonlinearly responds to climate conditions during the ripening period, and a large amount of variation exists that is not explained by climatic factors (figure 2). To deal with such variation, we adopted Bayesian calibration for the estimation of the parameter values, p i (i = 1-6) and n for each of the seven prefectures of Kyushu. The general procedure for Bayesian calibration begins by quantifying the known uncertainty of a parameter value in the form of a prior distribution. Observed data corresponding to model outputs are then used to update the posterior distribution of the parameters by means of Bayes' theorem: where p(θ |D) is the posterior distribution of the parameter θ for given data D, π(D|θ) is the likelihood function, p(θ ) is the prior distribution of parameter θ , and the denominator of the right-hand side of the equation (7) is the normalizing constant. The non-informative uniform distributions were used here for the prior distributions of all parameters. The likelihood function was developed on the assumption that errors were distributed normally: where σ 2 is the variance of the error, N is the sample size, and Y andŶ are vectors of the observed and modeled rice quality, respectively. We used the Metropolis-Hastings algorithm to estimate a high-dimensional posterior distribution of parameters via a sampling procedure using the Markov chain Monte Carlo technique (MCMC; Metropolis et al 1953, Hastings 1970. We applied the Metropolis-Hastings algorithm following the procedure described in Iizumi et al (2009).

Posterior distributions of parameters
The convergence of the Markov chains to a stationary distribution was examined by checking the Gelman-Rubin statistic (Gelman and Rubin 1992) on the basis of three parallel chains and visually checking the chains. The total number of MCMC iterations was 100 000. Once the chains had reached convergence with reference to the Gelman-Rubin statistic (<1.2), the last 10 000 samples per chain (i.e. 30 000 samples in total) were used to obtain the posterior distribution (see supplementary information available online at stacks.iop.org/ ERL/6/034031/mmedia).
The posterior distributions of parameters for Fukuoka estimated from the full dataset and from two subsets of the calibration data are shown in figure 3. In particular, these subsets excluded the data from the year with the hottest (1999, +1.8σ , where σ represents the standard deviation) or coldest (1980, −2.6σ ) summers (represented by the mean daily minimum temperature for n days after heading, T ) for the 25 yr period to examine the sensitivity of posterior parameter distributions to a particular data set from years with an extremely hot or cold temperature condition. For the parameters p 3 , p 4 and n, little difference was found between the locations of the posterior distributions from the subsets and that from the full set of calibration data, indicating a comparatively low sensitivity of these parameter values to data for a particular year.
For the other parameters, the locations of the posterior distribution varied between the subsets and the full dataset, indicating a comparatively high dependency of these parameter values on the particular set of calibration data. The parameters p 1 and p 2 correspond to the slope and intercept term, respectively, to express the linear effect of temperature on the lower limit of rice quality (equation (2)). For these parameters, the posterior distributions from the subset that excluded data from years with cold summers are very close to that from the full dataset, whereas the posterior distributions from the subset that excluded the data from years with hot summers shifted remarkably away from the distribution for the full set of calibration data. This shows that the presence of the data from years with hot summers in the calibration data is essential to precisely determine the lower limit of rice quality in this model.
The parameters p 5 and p 6 are the slope and intercept terms, respectively, of f b (T ). f b (T ) denotes the temperature dependence on the threshold value of total radiation for the ripening period, S, that results in the value of rice quality halfway between the upper and lower limits. For these parameters, the posterior distributions from all data were different from those from both subsets. These differences show that the values of these parameters are sensitive to the particular set of calibration data. Both data from years with hot summers and those with cold summer are essential to precisely determine the f b (T ) because both upper and lower values of rice quality are definitely important to determine their half value. Therefore, the lack of either type of data in the calibration data could lead to bias in the parameter values of p 5 and p 6 .
Of the seven parameters, the one that can be directly compared with the results of previous studies is the length of the period after heading in which temperature conditions negatively impact rice quality, n. The posterior mean value of n was 30 days, and this is close to other reported values (Nagato and Ebata 1960, Terashima et al 2001, Kondo et al 2006.

Model validation
We validated the capability of the model to simulate observed rice quality for each prefecture from years that were not included in the calibration data. The leave-one-out crossvalidation method (Stone 1974, Geisser 1975) was used.  Specifically, we first removed the sample data from one year of the calibration data and estimated the posterior distributions. Then the model was used to simulate the removed data. We repeated these steps for all years.
A comparison of the observed and simulated rice quality in Fukuoka for years in which data were removed is shown in figure 4. Most observed rice quality was distributed within the range of the ensemble mean ± 1 standard deviation (σ ) calculated by perturbing the parameter values within the posterior distributions. Pearson's correlation coefficient between the simulated and observed rice quality data for the 29 yr period was 0.86 (P < 0.001). The corresponding rootmean-square error was 12.33%. For other prefectures, the calculated goodness-of-fit statistics were somewhat worse than those for Fukuoka, but showed good correspondence between model simulations and corresponding observations (table 1). These results indicate a high capability for the model to capture temporal trends and interannual variation in observed rice quality from climatic factors.
Relatively large discrepancies between simulated rice quality and sample data were found in some years, for example, 2005 and 2007 in Fukuoka. These discrepancies can be attributed to factors such as pests, which are not accounted for in the model. Larger than normal outbreaks of brown planthoppers occurred in 2005 and 2007 in Kyushu Figure 5. Elasticity of rice quality to temperature or radiation at various temperature levels (represented by the mean daily minimum temperature for 30 days after heading). , Watanabe et al 2007, Kajisa et al 2008. This suggests that the model will not perform well in years in which non-climatic factors (e.g. pests) are the dominant cause of rice quality decline.

Relative impacts of climatic factors on rice quality
To quantify the relative impacts of climatic factors on rice quality, we performed sensitivity analysis using artificial increases in temperature for Fukuoka. More specifically, we calculated the change in rice quality per unit change in climatic factor as follows (referred to as the elasticity of rice quality to temperature or radiation): We used the posterior mean parameter values for the calculations. A positive sign for elasticity means a positive correlation between rice quality and the climatic factor and vice versa.
In this study, we focus only on elasticity of rice quality to temperature and radiation with change in temperature level because little information on the likely effect of climate change on radiation is available. We calculated the elasticities numerically with artificial temperature data. The artificial temperature data were obtained by adding anomalies to the baseline calculated from the calibration data (=20.3 • C). The anomalies ranged from −2 to +3 • C in intervals of 0.1 • C. The total radiation during the ripening period was kept constant (=688.9 MJ m −2 ). This corresponds to the baseline radiation.
The calculated elasticity of rice quality to temperature or radiation at various temperature levels is shown in figure 5. The sign of elasticity is always negative for temperature and positive for radiation, suggesting that a reduction in rice quality is caused by temperature increase, radiation decrease, or a combination of both. This finding agrees with the results of previous studies (Matsushima andManaka 1957, Kawatsu et al 2007). Under current temperature conditions, the negative impact of temperature and positive impact of radiation have roughly the same level of elasticity (−5.45 for temperature and 4.78 for radiation), but the negative impact from temperature is slightly stronger with increasing temperature than the positive impact from increasing radiation. For current radiation and management conditions, the tipping point at which the negative impact from increasing temperature becomes 1.5 times larger than the positive impact from radiation is 21.8 • C (corresponds to a temperature increase by +1.5 • C).

Prediction of rice quality under unprecedented climates
The capability of the model to predict rice quality under unprecedented climates was evaluated, taking the year 2010 with an extremely hot summer as an example. For each of the seven prefectures of Kyushu, we estimated the posterior distributions of parameter from the calibration data without the data from 2010 and then simulated the rice quality from the weather data for that year.
The correspondence between the observed and simulated data in 2010 is good in most areas ( figure 6). Although the model overestimated rice quality in 2010 in most prefectures, a similar tendency was also observed in other years. A comparatively large discrepancy between the observed and simulated data for 2010 appeared in Saga. This is likely due to the rapid introduction of a high-temperature tolerant cultivar 'Sagabiyori' in this area since 2009. Indeed, the proportion of first-grade rice in 2010 was 14.6% for the conventional cultivar 'Hinohikari', but 79.1% for 'Sagabiyori' (MAFF 2010a). No persistence of relevant information on a rapid change in the predominant cultivar in the calibration data explains the inaccurate simulation in this area.
We further examined the sensitivity of the predictive capability of the model to the calibration data. The Saga data were omitted from this analysis because of the inaccuracy introduced by the change in cultivar. We first obtained the frequency distribution for mean daily minimum temperature for n days after heading from the calibration data and calculated the standard deviation (σ ) on the assumption that the frequency distribution could be approximated by a normal distribution. The model was then calibrated for each of four subsets of the calibration data, referred to as CTL, +1.5σ , +1.0σ , and +0.5σ . These subsets except CTL were excluded from years in which the mean daily minimum temperature for n days after heading was greater than the calibration data mean + 1.5σ , mean + 1.0σ , and mean + 0.5σ , respectively, although the CTL subset was not excluded. This meant that in the +1.5σ subset, no data from years with very hot summers were included in the calibration data. Therefore, there is no information on the effect of very high-temperature conditions on rice quality in the model calibration. As the sample size effects calibration results, for a fair comparison the sample size was set to be the same among the subsets by removing samples less than +0.5σ . Finally, we compared the simulation results from the different calibration data with the observed data.
The prediction error for 2010 from each subset of calibration data is shown in figure 7. In Kagoshima, the model accurately simulated the rice quality in 2010 even when data from years with very hot summers (>mean + 1.5σ ) were removed from the calibration data. The correspondence between observed and simulated data deteriorated if we removed data from years with hot summers (>mean + 1.0σ ) or slightly hot summers (>mean + 0.5σ ). Therefore, the predictive capability of the model is somewhat dependent on the calibration data, but this dependency does not necessarily mean that rice quality cannot be predicted for years with extremely hot summers that are not in the calibration data.

Conclusion
We propose a statistical model to predict rice quality from climatic factors at large spatial scales. The model accounts for the multiple effects of temperature and radiation during the ripening period. Bayesian calibration was adopted to account for uncertainty due to non-climatic factors in the model prediction. The model accurately reproduced the temporal trend and interannual variation in observed rice quality. However, the model was inaccurate in the occasional years in which non-climatic factors dominated the quality results.
The sensitivity analysis showed that an increase in temperature has a negative effect on rice quality, whereas an increase in radiation has a positive effect. Under present climate conditions, these two climatic factors affect rice quality Figure 7. Absolute prediction errors between the observed and simulated ensemble mean rice quality for three areas in Kyushu in 2010, a year with an extremely hot summer. For each area, the horizontal axis indicates subsets of the calibration data in which data from very hot (+1.5σ ), hot (+1.0σ ), or slightly hot (+0.5σ ) summers were removed. For reference, the result from calibration data including data in very hot summer years is also shown (CTL). Error bars indicate the range of the ensemble mean ± 1 standard deviation produced by perturbing the parameter values.
to a similar extent. However, the negative effect from temperature becomes larger compared to the positive effect from radiation as average temperature during the ripening period increases. This suggests that climate change will cause a decline in rice quality, all other things being equal.
The predictive capability of the model is somewhat dependent on the calibration data. However, the model is still reliable even when data from years with very hot summers were not included in the calibration data, indicating that at least a modest level of extrapolation for future climate is possible. Some projection results for regional climate change impacts based on climate model projection have been reported in a separate paper (Okada et al 2011).
Future research should examine the impacts of atmospheric CO 2 concentration on rice quality. Such information is currently scarce, and datasets of the spatial distribution of atmospheric CO 2 concentration do not exist in the same way as they do for temperature and radiation. Additional systematic exploration of the sensitivity of the predictive capability of the model to the calibration data would also assist in determining the applicability of the model to a wider temporal domain.