Xuhui Wang
email: xuhui.wang@pku.edu.cn

Philippe Ciais

Laurent Li

Francois Ruget

Nicolas Vuichard

Nicolas Viovy

Feng Zhou

Jinfeng Chang

Xiuchen Wu

Hongfang Zhao

Shilong Piao
email: slpiao@pku.edu.cn

Management outweighs climate change on affecting length of rice growing period for early rice and single rice in China during 1991-2012

Keywords: length of growing period, rice, climate change, crop management, ORCHIDEE-crop, China

Whether crop phenology changes are caused by change in managements or by climate change belongs to the category of problems known as detection-attribution. Three type of rice (early, late and single rice) in China show an average increase in Length of Growing Period (LGP) during 1991-2012: 1.0±4.8 day/decade (±standard deviation across sites) for early rice, 0.2±4.5 day/decade for late rice and 2.0±6.0 day/decade for single rice, based on observations from 141 long-term monitoring stations. Positive LGP trends are widespread, but only significant (P<0.05) at 25% of early rice, 22% of late rice and 38% of single rice sites. We developed a Bayes-based optimization algorithm, and optimized five parameters controlling phenological development in a process-based crop model (ORCHIDEE-crop) for discriminating effects of managements from those of climate change on rice LGP. The results from the optimized ORCHIDEE-crop model suggest that climate change has an effect on LGP trends dependent on rice types. Climate trends have shortened LGP of early rice (-2.0±5.0 day/decade), lengthened LGP of late rice (1.1±5.4 day/decade) and have little impacts on LGP of single rice (-0.4±5.4 day/decade). ORCHIDEE-crop simulations further show that change in transplanting date caused widespread LGP change only for early rice sites, offsetting 65% of climate change induced LGP shortening. The primary drivers of LGP change are thus different among the three types of rice. Management are predominant driver of LGP change

for early and single rice. This study shows that complex regional variations of LGP can be reproduced with an optimized crop model. We further suggest that better documenting observational error and management practices can help reduce large uncertainties existed in attribution of LGP change, and future rice crop modeling in global/regional scales should consider different types of rice and variable transplanting dates in order to better account impacts of management and climate change.

Introduction

The Length of the Growing Period (LGP), defined as the interval in days from the day of planting/transplanting to the day of maturity, is an integrated indicator of crop development that has been related to production [START_REF] Bassu | How do various maize crop models vary in their responses to climate change factors?[END_REF], Zhang & Tao, 2013). Shortening LGP caused by warmer climate is recognized as a key emerging response through which climate change may impact agricultural production [START_REF] Bassu | How do various maize crop models vary in their responses to climate change factors?[END_REF][START_REF] Estrella | Trends and temperature response in the phenology of crops in Germany[END_REF][START_REF] Lin | Climate change impacts on crop yield and quality with CO2 fertilization in China[END_REF], Porter et al., 2014). However, historical change in LGP has been reported diversely across different crops and regions. Some studies found shortening LGP over the past decades [START_REF] Chmielewski | Climate changes and trends in phenology of fruit trees and field crops in Germany, 1961-2000[END_REF][START_REF] He | Impacts of recent climate warming, cultivar changes, and crop management on winter wheat phenology across the Loess Plateau of China[END_REF][START_REF] Siebert | Spatio-temporal patterns of phenological development in Germany in relation to temperature and day length[END_REF], Tao et al., 2014b[START_REF] Xiao | Observed changes in winter wheat phenology in the North China Plain for 1981-2009[END_REF]. For example, oat in Germany was found to have shorter LGP over the past five decade with rates of change ranging from -0.1 to -0.4 day/decade [START_REF] Siebert | Spatio-temporal patterns of phenological development in Germany in relation to temperature and day length[END_REF]. On the other hand, there are also studies finding little change or even a lengthening in LGP [START_REF] Liu | Contrasting effects of warming and autonomous breeding on single-rice productivity in China. Agriculture, Ecosystems &amp[END_REF][START_REF] Liu | Contributions of climatic and crop varietal changes to crop production in the North China Plain, since 1980s[END_REF][START_REF] Sacks | Crop management and phenology trends in the U.S. Corn Belt: Impacts on yields, evapotranspiration and energy balance[END_REF], Tao et al., 2013, Zhang et al., 2013).

For example, maize in the US Corn Belt shows lengthening LGP during 1981-2005 with an average positive trend of 5 day/decade [START_REF] Sacks | Crop management and phenology trends in the U.S. Corn Belt: Impacts on yields, evapotranspiration and energy balance[END_REF].

The LGP change of China's rice (Oryza sativa), which is the staple food resource for more than half of Chinese population and the crop with the largest growing area in the country, has attracted research interest. Observed trends of rice LGP across different stations vary largely from -2 day/decade to more than 7 day/decade over the past 2-3 decades, the majority of the field-scale observations showing either non-significant change or a lengthening of LGP [START_REF] Liu | Contributions of climatic and crop varietal changes to crop production in the North China Plain, since 1980s[END_REF][START_REF] Tao | Climate changes and trends in phenology and yields of field crops in China, 1981-2000[END_REF], Tao et al., 2013). One hypothesis explaining the lack of evidence for shortening trend of rice LGP was that management practices has counterbalanced the effects of climate change (e.g. [START_REF] Liu | Contrasting effects of warming and autonomous breeding on single-rice productivity in China. Agriculture, Ecosystems &amp[END_REF], Tao et al., 2013, Zhang et al., 2013).

However, large uncertainties remain on the relative contributions of climate change, shifts in transplanting date and other management practices (e.g. use of longer-duration cultivar), which limits our ability to understand the past trends and project the near term evolution of

LGP and its possible consequences for future crop production. Attribution of the observed trend of LGP from past observations remains challenging because both changes in climate and in management practices have taken place simultaneously. Recent studies used statistical models to characterize the interannual sensitivity of rice LGP to temperature and to separate the contribution of the temperature trend to LGP trend for rice and maize crops over the period 1981-2009 (Tao et al., 2014a, Tao et al., 2013, Zhang et al., 2013). This approach has some limitations: first, statistical models built from interannual LGP variations cannot isolate the impact of changing planting dates from the effects of climate change; second, statistical analyses usually assume linear and constant response to climatic variations (Zhang et al., 2013), but several studies showed that the response is neither linear [START_REF] Lobell | The critical role of extreme heat for maize production in the United States[END_REF] nor constant with time [START_REF] Lobell | Greater Sensitivity to Drought Accompanies Maize Yield Increase in the U[END_REF][START_REF] Burke | Adaptation to climate change: evidence from US agriculture Stanford University[END_REF]. On the other hand, crop models can provide an alternative mean to further understand mechanisms and quantify the attributions of different drivers (e.g. [START_REF] Lobell | Extreme heat effects on wheat senescence in India[END_REF]. Therefore, a question to ask in complement of the statistical models is whether crop models can be used as an independent method to separate climate change impacts from management. Using crop models factorial simulations where each driver is varied at a time, or combined, instead of statistical models based on historical data can overcome the limitations by having mechanistic representation of climate change impacts [START_REF] Gregory | Attribution of climate change: a methodology to estimate the potential contribution to increases in potato yield in Scotland since 1960[END_REF], but earlier application of crop models for the attribution of rice LGP trends were criticized for lack of validation for the study region [START_REF] Tao | Single rice growth period was prolonged by cultivars shifts, but yield was damaged by climate change during 1981-2009 in China, and late rice was just opposite[END_REF].

Methods

Rice phenology observations from Agrometeorological stations

Transplanting and maturity date of rice in China during 1991-2012 were recorded over 287 agro-meteorological field stations by the Chinese Meteorological Administration, covering the entire rice growing area, from the northeast to the southwest and Hainan Island (Fig. 1). The length of These observations were made following a standardized protocol across sites (CMA, 1993). The dataset includes single rice (177 stations), early rice (110 stations) and late rice (110 stations). Early rice and late rice have the same number of stations because they are two consecutive crops on the same site comprising the double rice cropping system (i.e. rotation between early rice and late rice [START_REF] Tao | Single rice growth period was prolonged by cultivars shifts, but yield was damaged by climate change during 1981-2009 in China, and late rice was just opposite[END_REF]). 80% of the 287 stations are used to optimize ORCHIDEE-crop model parameters. Time coverage of the stations ranges from few years to 21 years (Fig. 1) with 141 stations (88 for single rice and 53 for early/late rice) having records longer than 15 years, which are the long-term stations used for the detection and attribution of LGP trends (Figure S3). (Frolking et al., 2002) that darker pixel has larger area of rice croplands. [START_REF] Krinner | A dynamic global vegetation model for studies of the coupled atmosphere-biosphere system[END_REF], simulating carbon, water and energy fluxes (e.g. photosynthesis, respiration and evapotranspiration) and includes an agronomical module simulating crop phenology, leaf area dynamics, growth of reproductive organs, carbon allocations and management impacts [START_REF] Wu | ORCHIDEE-CROP (v0), a new process based Agro-Land Surface Model: model description and evaluation over Europe[END_REF]. The formula for crop phenology, leaf area dynamics, growth of reproductive organs were originated from a generic crop model STICS [START_REF] Brisson | An overview of the crop model stics[END_REF]. Compared with ORCHIDEE-STICS [START_REF] Gervois | Including Croplands in a Global Biosphere Model: Methodology and Evaluation at Specific Sites[END_REF], an earlier version of the crop model, which chained the ORCHIDEE model with STICS only through leaf area dynamics, ORCHIDEE-crop has a complete coupling between crop growth and physiology of carbon and water exchanges in soil-vegetation-atmosphere continuum. ORCHIDEE-crop calculates thermal unit accumulation, photosynthesis and energy exchange on a half-hourly time step, while leaf area dynamics, carbon allocation and biomass and soil organic carbon change are simulated on a daily time step.

Like most crop models, the crop growth cycle in ORCHIDEE-crop is divided into several stages with the developments driven by accumulated thermal unit. Since simulation of rice growth starts from transplanting (LEV), the growth cycle is divided into only three phases, which are divided by the onset of grain filling (DRP) and the physiological maturity (MAT).

The thermal unit (gdd) needed to grow from transplanting to maturity are prescribed parameters (GDD LEVDRP and GDD DRPMAT ). Accumulation of thermal unit (gdd) is calculated at each half-hourly time step following Eq. 1:

= × × × × min , + 1 - . 1
Where f(t) is a tri-linear function of temperature (T) following Eq. 2, δ p (δ v , δ n , δ w ) are crop-specific scalars for photo-period (vernalization, nitrogen, water) regulations respectively.

ε is a scalar parameter describing the sensitivity of the crop to nitrogen and water stress.

= 0, < ! "# > %& -! , ! < < ' ( ' ( -! ' ( -%& × -%& , ' ( < < %& . 2
As described above, temperature change has a first-order control over gdd (Fig. S1).

Therefore, the most important parameters for accumulations of gdd are GDD LEVDRP , GDD DRPMAT , T min , T opt and T max (Table 1), which are to be optimized in the parameter optimization. Details of the regulation scalars can be found in Brisson et al. (2008). In our study, δ v =1 because transplanted rice require no vernalization to develop; we assumed that δ n =1 and δ w = 1 because 93% of rice cropland in China is irrigated (http://www.knowledgebank.irri.org/country-specific/asia/rice-knowledge-for-china/2013-06-03-07-15-17, [START_REF] Salmon | Global rain-fed, irrigated, and paddy croplands: A new high resolution map derived from remote sensing, crop inventories and climate data[END_REF], and the nitrogen fertilizer application rate is higher than 100 kgN ha -1 [START_REF] Zhou | A New High-Resolution N2O Emission Inventory for China in 2008[END_REF]. In this study, we also assumed δ p =1, which indicates that photoperiodic constraint on the phenology is minimal for rice. This is probably true for early and single rice, because varieties insensitive to day-length change are commonly used (Cao et al., 2011). There are, however, cases for late rice, where day-length sensitive varieties are used (Cao et al., 2011), but we cannot account it due to lake of information on the extent for application of day-length sensitive varieties. Further details on ORCHIDEE-crop structure and parameters can be found in [START_REF] Wu | ORCHIDEE-CROP (v0), a new process based Agro-Land Surface Model: model description and evaluation over Europe[END_REF]. It should be noted that rice phenology development is modelled mostly by temperature driven processes in almost all rice models [START_REF] Li | Uncertainties in predicting rice yield by current crop models under a wide range of climatic conditions[END_REF], so the parameter we chose here represent the main processes driving the phenology development. [START_REF] Liu | Contrasting effects of warming and autonomous breeding on single-rice productivity in China. Agriculture, Ecosystems &amp[END_REF], Tao et al., 2013, Zhang et al., 2013), but it could cover other changes in agronomic practice, such as fertilization change.

Parameter optimization with particle filter

We used a particle filter method with sequential importance resampling (PFSIR) to optimize the ORCHIDEE-crop parameters for early, late and single rice phenology respectively over China. Particle filter is a Monte-Carlo implementation of recursive Bayesian theorem to estimate the posterior probability density of a state-space (here is the parameter set of the model) [START_REF] Van Leeuwen | Particle Filtering in Geophysical Systems[END_REF]. A set of ensemble members of the parameter set called "particles" hereafter, are used as a discrete approximation of the multi-dimensional probability density function (PDF) of the model parameters. The PDF evolves by propagating all particles forward in space or time through the ORCHIDEE-crop model. Each step when observations become available, each particle is assigned a weight (or importance) according to the model-observation differences. A new set of particles is generated after each iteration by resampling the weighted particles (sequential importance resampling). The optimized parameter sets for early rice, late rice and single rice are obtained through applying PFSIR to ORCHIDEE-crop model respectively. Particle filters has been found to have broader suitability than traditional variational methods [START_REF] Chorin | Conditions for successful data assimilation[END_REF], in particular for non-linear cases. Thus, variant forms of particle filter have become growingly popular when applying in earth system models (e.g. [START_REF] Bilionis | Crop physiology calibration in CLM[END_REF], Yu et al., 2014). Further details of PFSIR used in this study can be found in the Appendix.

Advantages of using the PFSIR method are multiple: First, unlike error minimization methods or manual adjustments previously adopted (e.g. [START_REF] Gregory | Attribution of climate change: a methodology to estimate the potential contribution to increases in potato yield in Scotland since 1960[END_REF], Zhang et al., 2014a), PFSIR not only provides a best (maximum likelihood) estimate, given an observation probability, according to the Bayes theorem, but also the uncertainties of the optimized parameters; Second, unlike variational methods (e.g. 4D-Var) assuming Gaussian distributions of the parameters, no assumptions are necessary for the posterior parameter distribution of the parameters in the particle filter, which makes it suitable for a model like ORCHIDEE-crop that uses some non-Gaussian and threshold-like parameters; Third, particle filter does not assume linearity of the state-space, which overcomes some of the limitations of methods based upon linearization of the state-space such as ensemble Kalman filter (van Leeuwen, 2010); Fourth, when being fed with large dataset, the Bayes-based particle filter is less sensitive to data outliers than error minimization methods (e.g. [START_REF] Kersebaum | Analysis and classification of data sets for calibration and validation of agro-ecosystem models[END_REF], which also make it suitable for application in crop models and over regional scale; Fifth, the particle filter does not require the effort of constructing the tangent linear model of the original model for calculating sensitivities of the model output to its parameters, and tends to have higher efficiency than other Monte-Carlo methods [START_REF] Gaucherel | Parameterization of a process-based tree-growth model: Comparison of optimization, MCMC and Particle Filtering algorithms[END_REF]. The particle filter is thus recommended for non-linear data assimilation, though has limitations for high-dimensional system [START_REF] Van Leeuwen | Particle Filtering in Geophysical Systems[END_REF]. With growing number of parameters (dimension of the parameter space), the filter may become less efficient and required a huge number of computing resources in order to obtain satisfactory estimates. Some improvements to the particle filter would be needed in such high-dimensional cases (e.g. [START_REF] Van Leeuwen | Nonlinear data assimilation in geosciences: an extremely efficient particle filter[END_REF]. Given the relatively small dimension of the parameter set (Table 1), this poses little threats to our study.

To evaluate the robustness of the optimized model, we randomly selected 20% of the sites It should be noted that we infer only one set of optimized parameter for each rice type over China, which is consistent with our intention to use a generic model across large regions, but could be a limitation in cases when local varieties within the same rice type have very different parameters. Separating the rice growing area into finer zones and producing multiple parameter sets for each rice type (Zhang et al., 2014a) may yield smaller errors due to increased degree of freedom and a potentially better calibration reflecting the diversity of local varieties. But doing this would also increase the risk of over-fitting and would require a detailed zoning map of rice crop varieties instead of zoning map of climate. In addition, there are growing requests for assessing climate change impacts over regional and global scales [START_REF] Rosenzweig | Assessing agricultural risks of climate change in the 21st century in a global gridded crop model intercomparison[END_REF] asking for robust parameter sets representing a broad scale of the growing area.

Trend analyses

We calculated the trend of rice LGP from the observations, the simulations S0 and S1, and for the residual ∆ by regressing time series of LGP at each station against year using least square regression. The trend estimates were compared with non-parametric test (Sen's slope) (Fig. S2). The similar estimates between least square regression slope and Sen's slope indicate robustness of the trend estimates to potential outliers. Statistical significance was reported based on two-tailed t-test. Only stations with more than 15 years of observations during 1991-2012 are used in the trend analyses (Fig. S3).

Results

Simulated LGP with prior and posterior parameters

Fig. 2 shows the histogram of the simulated bias of LGP (difference between observed

LGP and simulated LGP) for simulation S0 before and after optimization, and for the three rice types. Over site-years used in optimization, the posterior model largely reduces the root mean square error (RMSE) from 32.7 days (prior) to 14.8 days for early rice (optimized) (Fig. 2a), from 108.9 days to 12.4 days for late rice (Fig. 2b), and from 73.7 days to 24.4 days for single rice (Fig. 2c). When we only look at spatial variations across sites (Fig. S4), we found that the posterior model not only reduces the absolute errors (indicated by the vicinity to 1:1 line), but also better reproduces the spatial LGP gradient among the sites used for optimization. The R 2 for the spatial gradient improves from 0.41 (P<0.01) to 0.55 (P<0.01) for early rice (Fig. S4a), from 0.00 (P=0.91) to 0.33 (P<0.01) for late rice (Fig. S4b), and from 0.21 (P<0.01) to 0.47 (P<0.01) for single rice (Fig. A2c). Interannual variations of LGP at the long-term sites used for optimization also show significant improvement for all three rice types (P<0.05) (Fig. S5). These show that given the structure of the ORCHIDEE-crop model, with the PFSIR optimization method, it is possible to find a set of parameters for each of the three rice types, which provides an improved fit to the LGP observations across sites and years.

To test whether the improvements due to optimization is limited to the sites chosen for optimization, we also use the prior and posterior model parameters in ORCHIDEE-crop runs at the cross-validation sites. The RMSE of LGP for the simulation S0 with prior parameters are 33.9 day for early rice, 113.0 day for late rice and 74.5 day for single rice, respectively.

The RMSE of LGP with posterior parameters at the cross-validation sites are 16.3 day for early rice, 10.2 for late rice and 19.2 for single rice, which are close to that over the optimization sites (Fig. 2d-f). Therefore, the RMSE reduction over the validation sites is similar to that over the optimization sites (Fig. 2d-f). The improved spatial gradients (Fig. S4d-f) and interannual correlation between observed and modeled LGP (Fig. S5d-f) also hold for the validation sites. Indeed, when we re-selected the sites used for optimization and running the particle filter once again for testing, we obtain a similar set of parameter set than the one presented in The optimization of ORCHIDEE-crop parameters not only significantly reduced the misfit with site observations but also significantly changed the simulated trend in LGP (Fig S4). For early and single rice, the trend in optimized LGP (-0.7±2.7 day/decade (mean ± standard deviation across sites) for early rice and -0.5±5.2 day/decade for single rice) differs by more than 60% (P<0.01) from the prior modeled LGP trend (-1.7±4.8 day/decade for early rice and -1.5±18.4 day/decade for single rice)(Fig. S6a andc). For late rice, the optimization even changes the sign of the simulated LGP trend and largely reduced the spatial variations of the trend (Fig. S6b). The average LGP trend for late rice is changed from -7.5±16.7 day/decade to 1.0±3.0 day/decade (Fig. S6b). The optimized model thus produces lengthening instead of shortening LGP for late rice. The LGP trend simulated by the optimized model is further analyzed in the section "attribution of LGP trends to climate change, transplanting date change and other management factors".

Optimized parameter sets

Fig. 3 shows the probability distribution of the five optimized parameters (see Methods section for descriptions of the parameters) of the ORCHIDEE-crop simulation for LGP before (prior) and after (posterior) optimization for early rice, late rice and single rice, respectively.

Optimized (posterior) parameters of thermal unit requirements (GDD LEVDRP and GDD DRPMAT )

show largest uncertainty reduction (UR) with a 90% error reduction in the standard deviation after optimization (Fig. 3a andb, Table 1) ±0.5 o C for single rice, Fig. 3c) with UR ~85%. The higher optimal T opt and T min values found for early rice, compared to single and late rice suggest that early rice must be more acclimated to the high temperature in spring and summer over southern China, which matches its geographical distributions (Fig. 1) and was not accounted in the prior values of these parameters. For all the three rice types, the posterior probability distribution of T max shows a large range (Fig. 3e) indicating that this temperature threshold that corresponds to the stop of phenology development is less well constrained from the LGP observations, likely because T max is a high-end threshold, which is not frequently reached in the historical period 1991-2012 (4 site-days for early rice, no site-day for late rice and 7 site-days for single rice).

Fig. 3. Histogram of the prior and posterior parameter distribution for early rice, late rice and single rice. The optimized parameters include (a) GDD LEVDRP , (b) GDD DRPMAT , (c) T opt , (d)

T min , and (e) T max (see Methods section for definitions and descriptions of the parameters). At country scale, observations show an average lengthening of LGP for all three types of rice (Fig. 4). Single rice sites show the largest lengthening rate of 2.0±6.0 day/decade (mean ± standard deviation in spatial variations), followed by early rice (1.0±4.8 day/decade) and late rice (0.2±4.5 day/decade). But there are large site-to-site variations in the observed LGP trend (Fig. S7). For single rice, 61% of the sites show a trend towards longer LGP, 50% of which are statistically significant (Fig. s7c). For early and late rice, the percentage of sites showing longer LGP is similar (58% and 55% for early and late rice respectively), but the percentage of significant positive trends was smaller than that for single rice (27% and 19% for early and late rice respectively). There is a large proportion of sites showing no significant change of

Attribution of LGP trends to climate change, transplanting date change and other

LGP (more than 50% for all three types of rice), indicating that LGP change is either weakly sensitive to climate change or compensated by effects of change in climate and managements.

To further understand the drivers of the LGP trends, we estimated the contribution of climate change alone from simulation S1, the contribution of transplanting date from the difference between simulation S0 and S1, and interpreted the contribution of all other management (OM)

as being caused by a non-modeled residual term ∆, as explained in the Method section. As Fig. 4 and Fig. 5 shows, the impacts of climate change on LGP change differs between the three rice types. For early rice sites using the simulation S1 with the optimized model, we infer an average shortening of LGP induced by climate change alone of -2.0±5.0 day/decade (Fig. 4). Except for sites in Hainan and Guangxi, the shortening of LGP in simulation S1 is widespread (71%) and significant at 41% of the early rice sites (Fig. S7j). However, for late rice, climate change alone leads to an average lengthening of the LGP of 1.1±5.4 day/decade, with 16% of the sites having a significant lengthening mostly in Hunan, Jiangxi and Fujian provinces (Fig. S7k). This positive LGP trend for late rice in response to climate change occurs in ORCHIDEE-crop because temperature during the growing season is reaching the optimum temperature of phenology development for late rice in southern China (Table 1). For single rice, the contribution of climate change to LGP trends shows regional differences.

Climate change is modeled to shorten LGP over northeastern China and high-altitude Yungui plateau over southwestern China, but to lengthen LGP in the middle and lower reach of Yangtze River basin (Fig. S7l). These regional contrasts for single rice LGP trends leads to a near neutral average impact of climate change on LGP trend across China (-0.4±5.4 day/decade, Fig. 4). Among all the sites, climate change is the dominant factor contributing to the observed LGP trend for 26% of early rice sites, 28% of late rice sites and 19% of single rice sites (Fig. 5).

We found that 66% of the early rice sites experienced earlier transplanting date during 1991-2012 (Fig. S8). From the difference between modeled LGP in simulation S0 and S1, we infer that the earlier shift of the transplanting date (-2.0±4.8 day/decade) alone, has lengthened the LGP of early rice by 1.3±5.5 day/decade (Fig. 4). But earlier transplanting practice have not been adopted widely for late rice and single rice sites, and the observation sites showing positive and negative trends in transplanting dates are of similar proportion for late rice and single rice (Fig. S8b andc). The magnitude of the average trend in transplanting date is also small for these two types of rice (-0.3±3.4 day/decade for late rice and 0.1±4.1 day/decade for single rice), which has minor impacts on the simulated LGP change in the S0-S1 difference (-0.1±5.0 day/decade for late rice and -0.1±1.7 day/decade for single rice, Fig. 4). Therefore, the earlier shift of transplanting date is the dominant factor contributing to the trend of LGP at 17% of early rice sites (Fig. 5a), and a minor driver of LGP trends for other rice types, being dominant at only 7% of the late rice sites (Fig. 5b) and 2% of the single rice sites (Fig. 5c). On average across sites, the role of other management practices (OM), inferred from the residual trend not explained by transplanting date and climate change, is found to be the predominant factor for LGP change for early and single rice. OM are identified to be responsible for a lengthening of LGP by 2.1±3.9 day/decade for early rice and 2.8±7.6 day/decade for single rice (Fig. 4). A great majority of the early rice sites (71%) and single rice sites (64%) show positive contributions of OM trends. Over 20% of early rice sites sand 27% of single rice sites, the OM induced LGP trend is statistically significant (P<0.05, Fig. S7d-f). On the contrary, OM contributes to a shortening of LGP for late rice by -0.8±5.8 day/decade (Fig. 4), with a significant LGP shortening in Hunan, Jiangxi, Guangdong and Fujian provinces (Fig. S7e). The dominant role of OM is prevalent in southern China provinces, such as Guangdong, Guangxi and Yunnan for both early rice and late rice (Fig. 5a-b). For single rice sites, OM is the predominant driver of the LGP trend from the northeast to the southwest at 78% of the sites (Fig. 5c).

Discussion

Our analyses of a large network of rice phenological observations with more than 100 long-term stations across rice growing area in China indicate that the LGP of single rice has become longer over the past two decades, which is consistent with a recent study focused on Northeast China and Yangtze River basin during 1980-2009 [START_REF] Tao | Single rice growth period was prolonged by cultivars shifts, but yield was damaged by climate change during 1981-2009 in China, and late rice was just opposite[END_REF]. Although site-to-site variations are large, all three rice types exhibit an average trend towards longer

LGP. The ORCHIDEE-crop model optimized upon observed LGP was run using factorial simulations, with either climatological (fixed) or observed transplanting dates, and variable climate. The results suggest that the primary factors driving the LGP trends are not the same among the three rice types.

We found that recent climate change considered as a single driver in the model, shortened the LGP of early rice (Fig. 4 & Fig. S7j), which is consistent with previous statistical modelling (Zhang et al., 2013) and process modeling based on four sites [START_REF] Liu | Contrasting effects of warming and autonomous breeding on single-rice productivity in China. Agriculture, Ecosystems &amp[END_REF].

For late rice, climate change appears to have induced little change or a lengthening of LGP, which is different from early rice [START_REF] Liu | Contrasting effects of warming and autonomous breeding on single-rice productivity in China. Agriculture, Ecosystems &amp[END_REF], Tao et al., 2013) and from some other temperate crops [START_REF] Lobell | Extreme heat effects on wheat senescence in India[END_REF]. This is because the optimized parameter values indicate a lower optimum temperature (23.4 ± 0.6 o C) for phenology development of late rice than for early rice. Late rice sites are mainly located in southern China where temperature after transplanting (around July and August) is higher than this optimal temperature of phenology development of late rice [START_REF] Li | Impacts of climatic condition and its change on the temperature suitability and potential productivity of rice in the Yangtze River delta of China[END_REF]. Thus, further warming beyond the temperature optimum will not accelerate the phenology development and cause a lengthening of LGP (Fig. S1, [START_REF] Yin | A Nonlinear Model to Quantify Temperature Effect on Rice Phenology and It's Application[END_REF]. It should be noted that the optimum temperature that we determined from PFSIR is consistent with statistical analyses of rice phenology observations in southern China [START_REF] Xie | Index Identification of Suitable Temperature at the Booting Stage and Accumulated Temperature over 10℃ during the Whole Growth Period in Rice in South China[END_REF] and with the incubation study [START_REF] Summerfield | Photothermal Responses of Flowering in Rice (Oryza sativa)[END_REF], but lower than that used in previous models [START_REF] Liu | Contrasting effects of warming and autonomous breeding on single-rice productivity in China. Agriculture, Ecosystems &amp[END_REF], Zhang et al., 2014b), parameters of which may have originally derived from earlier studies based on assumptions or rice varieties in Southeast Asia (e.g. [START_REF] Kropff | ORYZA1, a basic model for irrigated lowland rice production[END_REF]. Our capability to further assess this parameter is rather limited since field trials determining the optimum temperature of phenology development are rarely available, requiring more data and future studies to refine this key parameter in order to more accurately project climate change impacts on LGP change. It should also be noted that, although high temperature stress did not necessarily shorten LGP, it may still adversely affect rice yields as it stresses photosynthesis [START_REF] Yin | C3 and C4 photosynthesis models: An overview from the perspective of crop modelling[END_REF], and thus reduce biomass accumulation for the harvest.

By comparing the simulations driven by fixed transplanting dates (S1) and by variable transplanting dates (S0), we can separate the contribution of transplanting date trends on LGP trends. Although an earlier transplanting date is a pragmatic autonomous adaptation through which farmers adapt to climate change [START_REF] Olesen | Impacts and adaptation of European crop production systems to climate change[END_REF], its effect on the regional trends of LGP was not separated by previous statistical models [START_REF] Tao | Single rice growth period was prolonged by cultivars shifts, but yield was damaged by climate change during 1981-2009 in China, and late rice was just opposite[END_REF], Zhang et al., 2013), probably due to its co-variations with climate [START_REF] Tao | Climate changes and trends in phenology and yields of field crops in China, 1981-2000[END_REF]. It may also be related with the linear assumption of previous statistical analyses (e.g. [START_REF] Tao | Single rice growth period was prolonged by cultivars shifts, but yield was damaged by climate change during 1981-2009 in China, and late rice was just opposite[END_REF]Zhang et al., 2013), which can be improved using recent progresses in statistical analyses including non-linear or threshold like equation (e.g. [START_REF] Burke | Adaptation to climate change: evidence from US agriculture Stanford University[END_REF][START_REF] Solomon | Climate econometrics[END_REF]. We found that changes in transplanting date were widespread over the last 20 years for early rice sites in southern China, and that they contributed to lengthen LGP, whereas climate change has the opposing effect of shortening LGP. This suggests that the adoption of earlier transplanting date has partly mitigated climate change impacts on early rice growth over the past two decades.

However, the same adaptation strategy is probably not possible for late rice because earlier transplanting and lengthening of LGP nearly compensate for each other for early rice, leaving no more time during the season available for earlier transplanting of late rice (MOA, 2014). In addition, advancing transplanting dates for late rice to mitigate climate change will have limitation due to frost risk and photo-period constraints in the autumn. The same reason may also explain why single rice sites show large site-to-site variations on the sign of change in transplanting date (Fig. S8).

Other management practices were found to be the dominant driver of LGP trends for early rice and single rice across the country (Fig. 5), which is about one magnitude larger than the contribution of transplanting date and climate trends for early rice and single rice, though with large site-to-site variations (Fig. 4). Previous studies usually interpreted this residual contribution not explained by climate change as the contribution of cultivar change, in particular the adoption of long-duration cultivars [START_REF] Liu | Contrasting effects of warming and autonomous breeding on single-rice productivity in China. Agriculture, Ecosystems &amp[END_REF], Tao et al., 2013, Zhang et al., 2013), which was supported by the empirical assessment of change in thermal requirements (Zhang et al., 2014b). This hypothesis is reasonable, since use of longer-duration cultivars is one of the most commonly used practices to achieve higher yields and mitigate the impacts of climate change [START_REF] Aggarwal | Climate change and rice yields in diverse agro-environments of India. II. Effect of uncertainties in scenarios and crop models on impact assessment[END_REF], Porter et al., 2014).

However, there are other management practices that could also impact LGP trends. For example, foliage nitrogen fertilizer spraying on leaf in the late growing season, can also lead to increase of leaf longevity and the growing season [START_REF] Russell | Plant Canopies: Their Growth, Form and Function[END_REF]. Future studies should account for these effects with spatially and temporally explicit datasets in order to more accurately attribute and project LGP change. In addition, OM trends may not necessarily induce longer LGP. Local agronomists in China have been studying and recommending the combination of rice varieties with shorter-duration and longer-duration cultivars in order to improve yield and to minimize risk of exposure to climate extremes (e.g. [START_REF] Ai | Ecological suitability of different combinations of hybrid rice varieties for double rice cropping in Hunan province[END_REF][START_REF] Mao | Effects of different ecological conditions on yield of different high quality rice varieties[END_REF][START_REF] Li | Adaptability test report of different double-cropping late Japonica varieties and their sowing dates in Miluo[END_REF] Shorter-LGP induced by OM seems to be widespread for late rice in southern China. These efforts were taken likely because shorter LGP for late rice can have the advantage to avoid the damage induced by cold weather events during anthesis and grain filling, known as the "cold dew wind" [START_REF] Huo | Agriculture and Meteorological Disasters[END_REF][START_REF] Wu | Geographical distribution of risk factors on cold dew wind of late rice in southern China[END_REF]. The risk of late rice exposure to cold damage can be more than 30% for some regions in southern China according to [START_REF] Wu | Geographical distribution of risk factors on cold dew wind of late rice in southern China[END_REF], and warming over past decades does not alleviate the risk of the weather events and reduce late rice production when it occurs [START_REF] Huo | Agriculture and Meteorological Disasters[END_REF], Ministry Of Agriculture, 2014).

Unlike previous studies identifying climate change as the dominant driver of rice phenology change, using field trials [START_REF] Vries | Adaptability of irrigated rice to temperature change in sahelian environments[END_REF], statistical models (Zhang et al., 2013) or crop model simulation [START_REF] Yao | Assessing the impacts of climate change on rice yields in the main rice areas of China[END_REF] not only can reproduce the phenology of the sites used for optimization, but also remains robust when applied to validation sites (Fig. 3). Therefore, the optimized model provides some confidence in the attribution, compared to models not optimized for rice croplands in China (e.g. [START_REF] Liu | Contrasting effects of warming and autonomous breeding on single-rice productivity in China. Agriculture, Ecosystems &amp[END_REF]. Indeed, the posterior model largely differs from the prior model in the estimated climate change impacts on LGP change (Fig. S6), further highlighting the necessity of optimizing crop models for regional studies. Admittedly, the optimized model simulations still cannot perfectly reproduce spatiotemporal variations in LGP, which may introduce uncertainties in the attribution of LGP trends to climate trends, but this should not largely impact our conclusions because we found no significant correlation between trend in the residual LGP (difference between observations and simulation S0) and the trend in growing season temperature (Fig. S9). This indicates that the trend attributed to OM is probably not biased by climate trend unexplained by ORCHIDEE-crop. On the other hand, in addition to optimizing the parameters of a single model against observations to reduce parameter uncertainties, recent studies indicate that multiple models can perform better than one model [START_REF] Li | Uncertainties in predicting rice yield by current crop models under a wide range of climatic conditions[END_REF][START_REF] Martre | Multimodel ensembles of wheat growth: many models are better than one[END_REF], due to the consideration of structural uncertainties. Although there are many difficulties in coordinating multiple models, promising future studies using model ensembles with the same protocol can improve our understanding regarding the structural uncertainties (e.g. [START_REF] Elliott | The Global Gridded Crop Model Intercomparison: data and modeling protocols for Phase 1 (v1.0)[END_REF]. It should also be noted that almost all current rice models, including ORCHIDEE-crop predict phenology development based on variations in temperature. The physiological impacts of non-temperature drivers should be further explored in future studies. Finally, observational error may also play an important role in the attribution to OM, which have largely been neglected both in our modelling study and previous statistical attribution (e.g. Zhang et al., 2013). Since the observation over all the stations followed the same protocol (CMA, 1993), it is often assumed that the observational error is uniform across sites and years. Thus, it would not impact the trend estimates and therefore attribution of the LGP trends. Although the assumption is reasonable, the reliability of this assumption remains uncertain. For better data-model integration, we recommend future data collection efforts to further report the error term together with the observations, which will help minimize potential biases in model parameterization and attribution efforts.

Conclusions

In

Appendix: Particle filter with sequential importance resampling

The basic idea of the particle filter is to represent the probability distribution function (PDF) of the parameters through an ensemble of parameters, each set of which is called a particle. At each step of the particle filter, the relative importance of the particle, or weight (w) is given by Eq. A1:

w ! = + ,|. ! ∑ + ,|. 0 1 023
. 41

where N is the number of particles, y is the observation and p(y|x i ) is probability density of the observations given the simulation with the particle x i (M(x i )) following Eq. A2: where x 1:N is the particle after N iterations. This Markov process makes the entire calculation iterative. When there is no observation in site i, the Markov process can still evolve by adding a random term to the particle in site i-1, but what we aim is to obtain final posterior PDF of the parameters given the observations over N sites (y 1:N ):

+(. 3:1 |, 3:1 = +(, 3:1 |. 3:1 +(. 3:1 +(, 3:1 ( . 44

Using Eq. A3 to further break down Eq. A4, we obtain Eq. A5:

+(. Applying Eq. A2 to Eq. A5, we obtained the numerical solution for all terms from 1 to N. For each step i, importance resampling is taking place to randomly redraw a new ensemble of particles from the weighted old ensemble to represent p(x i ), which leads to disregard particles that have very small weights and thus refine the ensemble. Sometimes the importance resampling may disregard some locally low probably particles but having global significance.

Therefore, we usually perform twice of the entire PFSIR process with different re-order observations to test its convergence in order to minimize the potential error due to this. More details and illustration of the particle filter can be found in van Leeuwen (2010). To adapt ORCHIDEE-crop model to different cropping systems, single rice and double rice (early rice and late rice) in China, we adapted particle filter with sequential importance resampling (van Leeuwen, 2009) separately for the three rice types (Table 1).

Prior parameters of the ORCHIDEE-crop was obtained from [START_REF] Irfan | Adaptation of the generic crop model STICS for rice (Oryza sativa L.) using farm data in Camargue[END_REF]. The range of prior parameters were obtained from Sanchez et al. ( 2014), which synthesized experiment knowledge on the range of basal, optimal and maximum temperature thresholds of rice development, and [START_REF] Valade | Modeling sugarcane yield with a process-based model from site to continental scale: uncertainties arising from model structure and parameter values[END_REF], which contains modeller's prior knowledge for the range of the parameters. Since we knew little about the prior probability distribution of the parameters, we assumed the prior parameter equally distributed within its range in order to guarantee a well sampling of the entire parameter space. Another issue may limit the capability of PFSIR is the error in the observation data. Unfortunately, accuracy description of the phenology observations are not available except that observations were made following the same standard protocol. However, the dataset is being treated as reliable data source without the need to do further filtering (e.g. [START_REF] Tao | Single rice growth period was prolonged by cultivars shifts, but yield was damaged by climate change during 1981-2009 in China, and late rice was just opposite[END_REF]Zhang et al., 2013).

  The first objective of this study is to optimize a process-based crop model to represent rice phenology in China. The second objective is to run the optimized model for attributing LGP change to climate change and change in various management practices during the last two decades. To achieve these goals, we first collected and harmonized observations of the rice LGP during 1991-2012 from an extensive station network in China (287 sites). Then, a random set of 80% of the sites is used to optimize the process-based crop model (ORCHIDEE-crop) under a Bayesian framework, by calibration of the parameters controlling rice phenology. The optimized model results are then evaluated against the remaining 20% of the site observations. Finally, contributions to LGP trends from climate change, transplanting date change and other management practices (including cultivar change) are separated by comparing the LGP observations and simulations of the optimized model driven by observed and fixed transplanting date.

Fig. 1 .

 1 Fig. 1. Spatial distribution of agrometeorological stations in China for (a) early rice, (b) late

  rice phenology with ORCHIDEE-crop model ORCHIDEE-crop model (svn version no. 2409) is a process-based crop model, which is based on the generic vegetation model ORCHIDEE

(

  22 sites of early rice, 21 sites of late rice and 35 sites of single rice, see Fig.1for its spatial distribution) as validation sites. The validation sites are not used into the PFSIR, providing independent evaluation measurements of the performance for the optimized model. It should be noted that the probability of posterior parameter distribution usually reflects the strength of constraint from the observation data, thus the spread of posterior probability distribution is also a metric to evaluate the performance of the particle filter. Larger spread of posterior probability distribution would indicate loose constraint from the observations.

Fig. 2 .

 2 Fig. 2. Histogram of the differences between observed length of rice growing period (LGP)

  , indicating strong observational constraints on these parameter values. Early, late and single rice have their posterior thermal unit requirements (GDD LEVDRP and GDD DRPMAT ) concentrated in a narrow range of values, which are significantly different from each other (P<0.05). On the other hand, the temperature threshold parameters regulating phenological development (T min , T opt , and T max in Eq. 2) show different values after optimization among the three rice types. For early rice, T min for phenology development is well constrained with a UR of 78% (9.9±0.5 o C, Fig.3d), while T opt has a large posterior range between 29 o C and 35 o C (32.3±1.9 o C, Fig.3c) and a UR of 55%. For late and single rice, optimized T min are slightly lower than early rice (9.2±1.1 o C for late rice and 9.4±0.5 o C for single rice, Fig.3d) and UR of 52% and 78%. On the contrary, optimized T opt for late and single rice are much lower than early rice (23.4±0.6 o C for late rice and 22.8

Fig. 4 .

 4 Fig. 4. Box plot of change in the length of rice growing period length (LGP) over the past two

Fig. 5 .

 5 Fig. 5. Spatial distribution of the controlling factors on change in the length of growing

  , our analyses combining phenology observations and optimized crop model simulations indicate that management practices (including both change in transplanting date and changes of OM) probably outweigh the impact of climate change on LGP change for early rice and single rice in China during the past two decades. However, we are only able to separate the effects on LGP trends of trends transplanting date from other management practices, such as cultivar change, due to limited data on spatio-temporal variations of other management practices. On the other hand, attribution of LGP trends to OM has the largest uncertainty in this analysis since the role of OM is inferred from the misfit of model runs driven by climate change and observed transplanting date and the observations. Errors in the attribution of LGP trends to climate or transplant date trends, depends largely on the crop model used, a structural bias in this model, and non-unified observational error across sites and years will translate into an erroneous attribution to OM. Through the Bayesian optimization framework (particle filter with sequential importance resampling), we optimized the ORCHIDEE-crop model to fit the spatio-temporal variations of LGP for the three rice types across China. The optimized model

  where δ is the observation error. In this study, we assume observational error is uniform across sites and years, since the observations over the network were made by trained staff following the same protocol(CMA, 1993), which are designed to unify and minimize the observational error across the network. Theoretically, it is possible to analytically have the PDF of the particles by putting all observations into the equation in one time. However, in practice, over a large number of sites/time steps, it requires a large number of particles to well sample the entire parameter space and computationally inefficient by wasting time in barely possible particles. Therefore, the Markov process (filter) to realize recursive Bayesian theorem is applied here (Eq. A3): +(. 3:1 = +(. 1 |. 163 +(. 163 |. 16: … +(. : |.3 (Eq. A3

Table 1 ,

 1 further indicating the robustness of the optimized models in reproducing the spatiotemporal variations of rice LGP in China during 1990-2012, for the three rice types.

Table 1 .

 1 Prior and posterior parameters for early rice, late rice and single rice.

	Prior		Posterior
	Generic rice	Early rice	Late rice	Single rice

  this study, we calibrated ORCHIDEE-crop model to represent spatio-temporal variations of rice LGP for three different types of rice in China, and applied this model forced by historical change in climate and transplanting date to attribute the trend in rice LGP. On one hand, we showed that, technically, 1) using Bayes-based particle filter, a generic process-based crop model can be objectively parameterized to represent spatio-temporal variations in rice LGP over China and 2) attribution of LGP trend based on calibrated model provides an alternative to statistical attribution previously used. On the other hand, through factorial simulations, we found that LGP change for different rice types show contrasting dominant drivers. Managements outweighs climate change in affecting LGP of early and single rice, but not for late rice. This suggests that future modelling efforts at global and regional scale should consider various types of rice grown and time-varying management

practices. Since large uncertainties still remain in understanding change in LGP, improving documentation of management practices in addition to transplanting date, better description of observational error and ensemble crop modelling can further reduce uncertainties in attributing climate change impacts in future studies.

Acknowledgements

We thank the editor and two anonymous reviewers for comments and help in improving the manuscript. This work was jointly supported by National Natural Science Foundation of China (NSFC) and Agence Nationale de la Recherche (ANR) (41561134016). The work of X. W. was supported by AXA foundation. The work of P.C. and X.W. was also supported by Imbalance-P (ERC Synergy Grant 610028). The work of F. Z. was supported by National Natural Science Foundation of China (41201077).

Supplementary Information Captions

Supplementary information can be found, in the online version of this article.