René Aïd

Luciano Campi
email: campi@ceremade.dauphine.fr

Adrien Nguyen

Nizar Touzi

A structural risk-neutral model of electricity prices *

Keywords: energy markets, electricity prices, fuel prices, risk-neutral probability, no-arbitrage pricing, forward contract. JEL Classification: D41, G13. AMS Classification (2000): 91B24, 91B26

The objective of this paper is to present a model for electricity spot prices and the corresponding forward contracts, which relies on the underlying market of fuels, thus avoiding the electricity nonstorability restriction. The structural aspect of our model comes from the fact that the electricity spot prices depend on the dynamics of the electricity demand at the maturity T , and on the random available capacity of each production means. Our model explains, in a stylized fact, how the prices of different fuels together with the demand combine to produce electricity prices. This modeling methodology allows one to transfer to electricity prices the risk-neutral probabilities of the market of fuels and under the hypothesis of independence between demand and outages on one hand, and prices of fuels on the other hand, it provides a regression-type relation between electricity forward prices and forward prices of fuels. Moreover, the model produces, by nature, the well-known peaks observed on electricity market data. In our model, spikes occur when the producer has to switch from one technology to the lowest cost available one. Numerical tests performed on a very crude approximation of the French electricity market using only two fuels (gas and oil) provide an illustration of the potential interest of this model.

Introduction

In security markets, the following relationship between spot and forward prices of a given security holds:

F (t, T ) = S t e r(T -t) , t ≤ T.
As usual, T is the maturity of the forward contract, S t is the spot price at t and r is the interest rate which is assumed constant for simplicity. We also assume no dividends. The no-arbitrage arguments usually used to prove such an equality lie heavily upon the fact that securities are storable at zero cost. For storable commodities (oil, soybeans, silver...), the former relation has been extended by including storage costs and and an unobservable variable called convenience yield (see Schwartz [START_REF] Schwartz | The stochastic behavior of commodity prices: Implications for valuation and hedging[END_REF], [START_REF] Routledge | Equilibrium forward curves for commodities[END_REF], and Geman [START_REF] Geman | Commodities and commodity derivatives: Modelling and pricing for agriculturals, metals and energy[END_REF], sec. 3.7). But, when one considers electricity markets (see Burger et al. [START_REF] Burger | Managing energy risk: An integrated view on power and other energy markets[END_REF] or Geman and Roncoroni [START_REF] Geman | Understanding the fine structure of electricity prices[END_REF] for an exhaustive description), such a property does not hold anymore: Once purchased, the electricity has to be consumed, so that the above relation does not make sense. This fact is very well documented in electricity market literature (see, e.g., Clewlow & Strickland [12]) but has not prevented the development of many electricity spot price models following the Black & Scholes framework [START_REF] Benth | A note on arbitrage-free pricing of forward contracts in energy markets[END_REF][START_REF] Benth | Extracting and applying smooth forward curves from average-based commodity contracts with seasonal variation[END_REF][START_REF] Benth | Stochastic modeling of electricity and related markets[END_REF][START_REF] Benth | Stochastic modeling of financial electricity contracts[END_REF][START_REF] Burger | A spot market model for pricing derivatives in electricity markets[END_REF][START_REF] Cartea | Pricing in electricity markets: a mean reverting jump diffusion model with seasonality[END_REF][START_REF] Eydeland | Energy and power risk management: New developments in modeling, pricing and hedging[END_REF] (see Benth [START_REF] Benth | Stochastic modeling of electricity and related markets[END_REF] for a survey of the literature).

Nevertheless, the non-storability of electricity is not enough to claim that no relation holds between spot and forward prices and that no arbitrage constraint affects the term structure of electricity prices, except the constraints coming from overlapping forward contracts. Indeed, one could argue that even if electricity can not be stored, the fuels that are used to produce electricity can. To see that this observation leads to constraints on the term structure of electricity prices, let us consider a fictitious economy in which power is produced by a single technology -coal thermal units with the same degree of efficiency -and that the electricity spot market is competitive. Then, the electricity price should satisfy the following relation :

F e (t, T ) = q c F c (t, T ), t ≤ T,
where the subscript e stands for electricity, c stands for coal, and q c denotes the heat rate. If there is t < T such that F e (t, T ) > q c F c (t, T ), then one can at time t:

• Sell a forward on electricity at F e (t, T ) and buy q c coal forward at F c (t, T ) and, at time T :

• Sell q c coal at S c (T ), buy electricity at S e (T ) = q c S c (T ).

One can check that this strategy is indeed an arbitrage. Moreover, the opposite relation can be obtained in a similar way. Here, in this fictitious economy, the important feature is not that electricity can be produced by coal, but that the relation between spot prices of coal and electricity is known. Furthermore, it extends directly to forward prices.

In real economies, similar no-arbitrage relations between electricity and fuel prices can not be identified so easily. The reason for this is that electricity can be produced out of many technologies with many different efficiency levels: Coal plants more or less ancient, fuel plants, nuclear plants, hydro, solar and windfarms, and so on. Generally, the electricity spot price is considered to be the day-ahead hourly market price. At that time horizon, any producer will perform an ordering of its production means on the basis of their production costs. This process refers to a unit commitment problem and one can find a huge literature on this optimization problem in power systems literature (see, e.g., Batut and Renaud [START_REF] Batut | Daily generation scheduling with transmission constraints: A new class of algorithms[END_REF] and Dentcheva et al. [START_REF] Dentcheva | Solving the unit commitment problem in power generation with primal and dual methods[END_REF]). Depending on the market prices of fuels and on the state of the power system (demand, outages, inflows, wind and so forth), this ordering may vary through time. Hence, when the forward contract is being signed, the ordering at the contract maturity is not known.

The objective of this paper is to build a model for electricity spot prices and the corresponding forward contracts, which relies on the underlying markets of fuels, thus avoiding the non-storability restriction. The structural aspect of our model comes from the fact that the electricity spot prices depend on the dynamics of the electricity demand at the maturity T , and on the random available capacity of each production means. Our model allows one to explain, in a stylized fact, how the prices of different fuels together with the demand combine to produce electricity prices. This modeling methodology allows us to transfer to electricity prices the risk-neutral probabilities of the market of fuels, under a certain independence hypothesis (see Assumption 2.2). Moreover, the model produces, by nature, the well-known peaks observed on electricity market data. In our model, spikes occur when the producer has to switch from one technology to the lowest cost available one. Moreover, the dynamics of the demand explain this switching process. Then, one easily understands that the spikes result from a high level of the demand which forces the producer to use a more expensive technology.

Our model is close to Barlow's [START_REF] Barlow | A diffusion model for electricity prices[END_REF], since the electricity spot price is defined as an equilibrium between demand and production. But, in our model, the stack curve is described by the different available capacities and not a single parametrized curve. Moreover, this model shares some ideas with Fleten and Lemming forward curve reconstruction method [START_REF] Fleten | Constructing forward price curves in electricity markets[END_REF]. But, whereas the authors methodology relies on an external structural model provided by the SINTEF, our methodology does not require such inputs.

The article is structured in the following way: Section 2 is devoted to the description of our model; Section 3 describes the relation between future prices of electricity and fuels; Section 4 presents the model in the case of only two fuels; Section 5 presents numerical results showing the potential of the model on the two technologies case of the preceeding section; finally, Section 6 provides some future research perspectives.

The Model

Let (Ω, F, P) be a probability space sufficiently rich to support all the processes we will introduce throughout this paper. Let (W 0 , W ) be an (n + 1)-dimensional standard Wiener process with W = (W 1 , . . . , W n ), n ≥ 1. In the sequel, we will distinguish between the filtration F 0 = (F 0 t ) generated by W 0 and the filtration

F W = (F W t ) generated by the n-dimensional Wiener process W = (W 1 , . . . , W n ).
Commodities market. We consider a market where agents can trade n ≥ 1 commodities and purchase electricity. We consider only commodities that can be used to produce electricity. With a slight abuse of language, we will always identify in this paper any given production technology with the corresponding commodity (also called fuel) used. For i = 1, . . . n, S i t denotes the price of the quantity of commodity i necessary to produce 1 KWh of electricity and is assumed to follow the following SDE:

dS i t = S i t   µ i t dt + n j=1 σ ij t dW j t   , t ≥ 0, (2.1)
where µ i and σ ij are F W -adapted processes suitably integrable (see Assumption 2.1). We also assume that the market contains a riskless asset with price process

S 0 t = e R t 0 rudu , t ≥ 0,
where the instantaneous interest rate (r t ) t≥0 is an F W -adapted non-negative process such that t 0 r u du is finite a.s. for every t ≥ 0. As a consequence, (r t ) is independent of the Brownian motion W 0 . We will frequently use the notation Xt := X t /S 0 t for any process (X t ). We make the following standard assumption (see, e.g. Karatzas [START_REF] Karatzas | Lectures on the mathematics of finance[END_REF], Section 5.6).

Assumption 2.1 The volatility matrix σ t = (σ ij t ) 1≤i,j≤n is invertible and both matrices σ and σ -1 are bounded uniformly on [0, T * ] × Ω. Finally, let θ denote the market price of risk, i.e.

θ t := σ -1 t [µ t -r t 1 n ], t ≥ 0,
where 1 n is the n-dimensional vector with all unit entries. We assume that such a process θ satisfies the so-called Novikov condition

E exp 1 2 T * 0 ||θ t || 2 dt < ∞ a.s.
Remark 1 Imposing the Novikov condition on the commodities market price of risk ensures that the minimal martingale measure we will use for pricing in Section 3 is well defined. The reader is referred to Section 5.6 in Karatzas's book [START_REF] Karatzas | Lectures on the mathematics of finance[END_REF].

Market demand for electricity. We model the electricity market demand by a real-valued continuous process D = (D t ) t≥0 adapted to the filtration F 0 = (F 0 t ) generated by the Brownian motion W 0 . Observe that, under our assumptions, the processes S i (i = 0, . . . , n) are independent under P of the demand process D. To be more precise, the process D models the whole electricity demand of a given geographical area (e.g. U.K., Switzerland, Italy and so on). In this respect, it must be strictly positive. Nevertheless, in Section 5, where the empirical analysis is performed, it is more convenient to use a residual demand to reduce the number of possible technologies. A residual demand is the whole demand less the production of some generation assets (e.g. nuclear power, run of the river hydrolic plants, wind farms). It is clear that the residual demand can be negative.

Electricity spot prices. We denote by P t the electricity spot price at time t. At any time t, the electricity producer can choose among the n commodities which is the most convenient to produce electricity at that particular moment and the electricity spot price will be proportional to the spot price of the chosen commodity. We recall that the proportionality factor is already included in the definition of each S i so that, if at time t the producer chooses commodity i then P t = S i t , 1 ≤ i ≤ n.

How does the electricity producer choose the most convenient commodity to use? For each i = 1, . . . , n, we denote ∆ i t > 0 the given capacity of the i-th technology for electricity production at time t. (∆ i t ) is a stochastic process defined on (Ω, F, P) and assumed independent of (W 0 , W ). We denote F ∆ = (F ∆ t ) its filtration. Moreover, we assume that each ∆ i t takes values in [m i , M i ] where 0 ≤ m i < M i are the minimal and the maximal capacity of i-th technology, both values being known to the producer. In reality, the producer has to fill capacity constraints, so he faces demand variability, security conditions and failures risk. Thus, if one wants to consider capacity management and partial technology failures in the model, the production capacity has to be modelled as a stochastic process.

For every given (t, ω), the producer performs an ordering of the commodities from the cheapest to the most expensive. The ordered prices of commodities are denoted by

S (1) t (ω) ≤ • • • ≤ S (n) t (ω).
This order induces a permutation over the index set {1, . . . , n} denoted by π t = {π t (1), . . . , π t (n)}. Notice that π t defined an F W -adapted stochastic process, and we follow the usual probabilistic notation omitting its dependence on ω.

Given a commodities order π t at time t, we set

I πt k (t) := k-1 i=1 ∆ πt(i) t , k i=1 ∆ πt(i) t , 1 ≤ k ≤ n,
with the convention

0 i=1 ≡ 0.
For the sake of simplicity, we will assume from now on that the electricity market is competitive and we will not take into account the short term constraints on generation assets as well as start-up costs. Hence, the electricity spot price equals the cost of the last production unit used in the stack curve (marginal unit). Thus, if the market demand at time t for electricity D t belongs to the interval I πt k (t), the last unit of electricity is produced by means of technology π t (k), when available. Otherwise, it is produced with the next one with respect to the time-t order π t . This translates into the following formula:

P t = n i=1 S (i) t 1 {Dt∈I π t i (t)} , t ≥ 0. (2.2)
Let T * > 0 be a given finite horizon, in the sequel we will work on the finite time interval [0, T * ]. Typically, all maturities and delivery dates of forward contracts we will consider in the sequel, will always belong to [0,

T * ]. Assumption 2.2 Let F t = F 0 t ∨ F W t ∨ F ∆ t , t ∈ [0, T * ],
be the market filtration. There exists an equivalent probability measure Q ∼ P defined on F T * , such that the discounted prices of commodities S = ( S1 , . . . , Sn ) (i.e. without electricity!) are local Q-martingales with respect to (F t ). This hypothesis is equivalent to assuming absence of arbitrage in the market of fuels (see [START_REF] Delbaen | A general version of the fundamental theorem of asset pricing[END_REF]). Notice that we are not making this assumption on the electricity market, as announced in the introduction. Thanks to relation (2.2), any electricity derivative can be viewed as a basket option on fuels. Hence, Assumption 2.2 allows us to properly apply the usual risk neutral machinery to price electricity derivatives.

The market of commodities and electricity is clearly incomplete, due to the presence of additional unhedgeable randomness source W 0 driving electricity demand D. Thus, in order to price derivatives on electricity we have to choose an equivalent martingale measure among infinitely many to use as a pricing measure. One possible choice is the following: Let Q = Q min denote the minimal martingale measure introduced by Föllmer and Schweizer [START_REF] Föllmer | Hedging of contingent claims under incomplete information[END_REF], i.e.

dQ dP = exp - T * 0 θ u • dW u - 1 2 T * 0 ||θ u || 2 du (2.3)
where we recall that θ t = σ -1 t (µ t -r t 1 n ) is the market price of risk for the commodities market (S 1 , . . . , S n ). In the previous formula as well as in the sequel of this paper x • y denotes the scalar product between two vectors x, y.

Notice that, due to Assumption 2.1, such a measure is well defined, i.e. (2.3) defines a probability measure on F T * , which is equivalent to the objective measure P.

Remark 2 Furthermore, it can be easily checked that under Q the laws of processes W 0 and ∆ i (1 ≤ i ≤ n) are the same as under the objective probability P and the independence between the filtrations F 0 , F ∆ and F W is preserved under Q. This property will be very useful in the proof of Proposition 3.1.

Under such a probability Q the prices of commodites S i , 1 ≤ i ≤ n, satisfy the SDEs

dS i t = S i t   r t dt + d j=1 σ i,j t dW Q,j t   , S i 0 > 0,
whose solutions are given by

S i t = S i 0 exp t 0 r u - 1 2 ||σ i u || 2 du + t 0 σ i u • dW Q u , t ≥ 0,
where

W Q = (W Q,1 , . . . , W Q,d
) is an n-dimensional Brownian motion under Q, and σ i = (σ i,1 , . . . , σ i,n ).

The measure Q will be used as pricing measure in the rest of the paper. We recall that in the literature, such a measure Q is related to locally risk minimization procedure, in the sense that, given a contingent claim H with some maturity T > 0, E Q [ H] is the minimum price allowing an agent to approximately (and locally in L 2 ) hedge the claim (see Schweizer's survey [START_REF] Schweizer | A guided tour through quadratic hedging approaches, Option pricing, interest rates and risk management[END_REF] for further details).

Remark 3 Notice that including storage costs c i and convenience yields δ i changes only the drift coefficients in commodities dynamics from r t to r t + c i -δ i .

Electricity forward prices

We now consider a so-called forward contract on electricity with maturity T 1 > 0 and delivery period [T 1 , T 2 ] for T 1 < T 2 ≤ T * , i.e. a contract defined by the payoff

(T 2 -T 1 ) -1 T2 T1 P T dT (3.4)
at the maturity T 1 , whose time-t price F t (T 1 , T 2 ) is to be paid at T 1 .

The following observation is crucial: According to formula 2.2, the payoff (3.4) can be expressed in terms of prices of fuels, so that in our model the forward contract on electricity can be viewed as a forward contract on fuels and since the classical no-arbitrage theory makes sense on the market of fuels, it can also be used to price electricity derivatives such as (3.4). In other terms, our production-based structural model relating electricity and fuels allows us to transfer the whole no-arbitrage classical approach from fuels to electricity market, so overcoming the non-storability issue.

By Assumption 2.2 and classical result on forward pricing (see [START_REF] Björk | Arbitrage theory in continuous time[END_REF] Chapter 26), it immediately follows that:

F t (T 1 , T 2 ) = 1 T 2 -T 1 T2 T1 E Q t e - R T t rudu P T E Q t e - R T t rudu dT, (3.5)
E Q t denoting the conditional Q-expectation given market's filtration F t , for t ≥ 0. Let T ∈ [T 1 , T 2 ].
It is convenient for the next calculations to introduce the forward measure Q T defined by the density

dQ T dQ := e - R T t rudu B t (T ) on F W T , where B t (T ) := E Q t e - R T t rudu
is the time-t price of a zero-coupon bond with maturity T . Then:

F t (T 1 , T 2 ) = 1 T 2 -T 1 T2 T1 E Q T [P T |F t ] dT (3.6) = n i=1 1 T 2 -T 1 T2 T1 E Q T S (i) T 1 {D T ∈I π T i (T )} |F t dT. (3.7)
We denote by Π n the set of all permutations over the index set {1, . . . , n}. Let π ∈ Π n be a given non-random permutation. Under the assumption S i t ∈ L 1 (Q t ) for any t ≥ 0 and 1 ≤ i ≤ n, we can define the following changes of probability on F W T :

dQ i T dQ T = S i T E Q T [S i T ] , 1 ≤ i ≤ n, T ≤ T * . Proposition 3.1 If our model assumptions hold and if S i T ∈ L 1 (Q T ) for all T ∈ [T 1 , T 2 ] and 1 ≤ i ≤ n, we have F t (T 1 , T 2 ) = 1 T 2 -T 1 n i=1 π∈Πn T2 T1 F π(i) t (T )Q π(i) T [π T = π|F W t ]Q T [D T ∈ I π i (T )|F 0,∆ t ]dT, (3.8)
for t ∈ [0, T 1 ], where F i t (T ) denotes the price at time t of forward contract on the i-th commodity with maturity T and F 0,∆ t is the natural filtration generated by both W 0 and ∆.

Proof. Notice first that

F t (T 1 , T 2 ) = 1 T 2 -T 1 T2 T1 F t (T )dT, where F t (T ) = E Q T [P T |F t ]
can be interpreted as the t-price of a forward contract with maturity T and instantaneous delivery at maturity. By the definition of electricity forward price F t (T ), we have

F t (T ) = n i=1 E Q T S (i) T 1 {D T ∈I π T i (T )} |F t = n i=1 π∈Πn E Q T S π(i) T 1 {D T ∈I π i (T )} 1 {π T =π} |F t .
If we use the mutual (conditional) independence between W , W 0 and ∆ as in Remark 2, we get

F t (T ) = n i=1 π∈Πn E Q T S π(i) T 1 {π T =π} |F W t Q T [D T ∈ I π i (T )|F 0,∆ t ].
Using the change of probability dQ

π(i)
T /dQ T yields

E Q T S π(i) T 1 {π T =π} |F W t = F π(i) t (T )Q π(i) T [π T = π|F W t ],
so giving, after integrating between T 1 and T 2 and dividing by T 2 -T 1 , the announced formula.

The main formula (3.8) provides a formal expression to the current intuition of electricity market players that the forward prices are expected to be equal to a weighted average of forward prices of fuels. Such weights are determined by the crossing of the expected demand with the expected stack curve of the technologies. We will see in Section 5 that this model is able to explain the spikes of electricity. Nonetheless, we can already observe that formula (3.8) reproduces the stylized fact that the paths of electricity forward prices are much smoother than those of spot prices. This is due to the averaging effect of the conditional expectation on the indicator functions appearing in formula (2.2), even in the degenerate case when the delivery period reduces to a singleton.

In the next section, we will perform some explicit computations of the conditional probabilities involved in the previous formula for electricity forward prices, under more specific assumptions on the dynamics of prices and demand.

A model with two technologies and constant coefficients

In order to push further the explicit calculations, we assume now that the volatilities of fuels are constant, i.e. σ i,j t = σ i,j for some constant numbers σ i,j > 0, 1 ≤ i, j ≤ n, and that the interest rate is constant r t = r > 0. Under the latter simplification, the forward-neutral measures Q T all coincide with the minimal martingale measure Q = Q min . Similar closed-form expressions can be obtained by assuming a Gaussian Heath-Jarrow-Morton model for the yield curve.

Let us assume from now on that only two technologies are available, i.e. n = 2.

Dynamics of capacity processes ∆ i . In order to get explicit formulae for forward prices we have to specify the dynamics of each capacity process ∆ i . We assume that the probability space (Ω, F, P) supports four (independent) standard Poisson processes N 1,u t , N

d∆ i t = (m i -M i )1 (∆ i t =Mi) dN i,d t + (M i -m i )1 (∆ i t =mi) dN i,u t , ∆ i 0 = M i (4.9)
Remark 4 Basically we are assuming that each capacity i can take only two values M i > m i and it switches from m i to M i (resp. from M i to m i ) when the process N i,u (resp. N i,d ) jumps. Each capacity evolves independently of each other. At t = 0 both technologies have maximal capacity M i . The fact that the intensities of upside and downside jumps of ∆ i are not necessarily equal introduces a skewness in the probability of being at capacity M i or m i .

Let T be any time in the delivery period [T 1 , T 2 ]. First observe that, since ∆ is independent of W 0 and its law is invariant under the probability change from

P to Q = Q T as in Remark 2, we have Q T [∆ π(1) T = x 1 |F 0,∆ t ] = P[∆ π(1) T = x 1 |∆ t ] as well as Q T [∆ π(1) T = x 1 , ∆ π(2) T = x 2 |F 0,∆ t ] = P[∆ π(1) T = x 1 , ∆ π(2) T = x 2 |∆ t ] for x 1 ∈ {m 1 , M 1 } and x 2 ∈ {m 2 , M 2 }.
As a consequence of the previous assumption on the dynamics of capacities ∆ i , the conditional probabilities

Q T [D T ∈ I π k (T )|F 0,∆ t
] appearing in the main formula (3.8) can be decomposed as follows

Q T [D T ∈ I π 1 (T )|F 0,∆ t ] = Q T D T ≤ ∆ π(1) T |F 0,∆ t = P[∆ π(1) T = m 1 |F ∆ t ]Q T D T ≤ m 1 |F 0 t +P[∆ π(1) T = M 1 |F ∆ t ]Q T D T ≤ M 1 |F 0 t A similar decomposition for Q T [D T ∈ I π 2 (T )|F 0,∆ t
] holds too. It is clear now that the building blocks appearing in such formulae are the probabilities

P[∆ k T = x|∆ k t ] and Q T D T ≤ y|F 0 t . It remains to compute P[∆ k T = x|F ∆ t ] for k = 1, 2 and x = M k , m k .
As an example, we will compute

P[∆ k T = m k |∆ 0 = M k ].
For the sake of simplicity, we will drop for a while the index k from the notation, that is we will write ∆ T for ∆ k T , M for M k , and so on. Let τ d be the last jump time of the process N d t before T , i.e. τ d = sup{t ∈ [0, T ] : ∆N d t = 1} with the convention that sup ∅ = 0. Notice that on the event {τ d > 0} we have

{∆ T = m} = {N u τ d = N u T }.
On the other hand, on the set {τ d = 0} the process ∆ has no jump downwards over the time interval [0, T ], so that P(∆ T = m, τ d = 0|∆ 0 = M ) = 0. Using the independence between N d and N u and the stationarity of N u , one has

P[∆ T = m|∆ 0 = M ] = E[P(N u τ d = N u T |τ d )1 τ d >0 ] = E[P(N u T -τ d = 0|T -τ d )1 T -τ d <T ] = E[e -λ u (T -τ d ) 1 T -τ d <T ].
By the time-reversal property of the standard Poisson process1 , the random variable T -τ d has the same law as T d 1 ∧ T , where T d 1 is the first jump time of (N d t ) t≥0 . We recall that T 1 has exponential law with parameter λ d . Thus we have

P[∆ T = m|∆ 0 = M ] = E[e -λ u (T d 1 ∧T ) 1 T d 1 <T ] = E[e -λ u T d 1 1 T d 1 <T ] = λ d λ d + λ u (1 -e -(λ d +λ u )T )
The general result follows by stationarity :

P[∆ k T = m k |∆ k t = M k ] = λ d k λ d k + λ u k (1 -e -(λ d k +λ u k )(T -t) ), k = 1, 2. (4.10)
Using the same arguments, one can obtain similar expressions for the remaining probabilities

P[∆ k T = x|F ∆ t ] for k = 1, 2 and x = M k , m k .
Dynamics of the electricity demand D. We also assume that the residual demand is defined by the a mean-reverting Ornstein-Uhlenbeck process. It is well-known that this process has a positive probability to be negative. Nonetheless, in the empirical study, it will applied to a residual demand, which can be negative (see Section 2).

dD t = a(b(t) -D t )dt + δdW 0 t , D 0 > 0, (4.11)
for given strictly positive constants a and δ, and a long-term mean b(t) which can vary with time, to incorporate annual seasonal effects as in [START_REF] Barlow | A diffusion model for electricity prices[END_REF] :

b(t) = b 0 + b 1 cos(2πt -b 2 ) - 2π a sin(2πt -b 2 ) ,
[D T ≤ x 1 |F 0 t ] and Q[x 1 < D T ≤ x 1 + x 2 |F 0 t ],
for any 0 ≤ t ≤ T and x 1 , x 2 ∈ R, given by

Q[D T ≤ x 1 |F 0 t ] = Φ   x 1 -b(T ) -(D t -b(t))e -a(T -t) δ 1 2a 1 -e -2a(T -t)   (4.12) Q[x 1 < D T ≤ x 1 + x 2 |F 0 t ] = Φ   (x 1 + x 2 ) -b(T ) -(D t -b(t))e -a(T -t) δ 1 2a 1 -e -2a(T -t)   (4.13) -Φ   x 1 -b(T ) -(D t -b(t))e -a(T -t) δ 1 2a 1 -e -2a(T -t)   ,
where Φ denotes the cumulative distribution function of an N (0, 1) random variable.

Let T ∈ [T 1 , T 2 ]
. The next step consists in computing the law of the couple (S 1 T , S 2 T ) under each probability Q

π(i) T
for any permutation π ∈ Π 2 and any i = 1, 2, in order to get an explicit expression for the conditional probability

Q T [π T = π|F W t ] = Q[π T = π|F W t ]
appearing in formula (3.8). It can be easily done in this setting by using multidimensional Girsanov's theorem (see, e.g., Karatzas and Shreve's book [START_REF] Karatzas | Brownian motion and stochastic calculus[END_REF], Theorem 5.1 in Chapter 3). Indeed, if we denote σ i the 2-dimensional vector (σ i,1 , σ i,2 ) and we set

Z i t := dQ i T dQ | F W t ,
we get that

Z i t = exp σ i • W Q t - 1 2 ||σ i || 2 t , t ∈ [0, T ].
A simple application of Girsanov's theorem provides the following Q i T -dynamics of each price process S j for j = 1, 2 :

S j t = S j 0 exp r - 1 2 ||σ j || 2 + σ j • σ i t + σ j • W t , t ∈ [0, T ],
where

W = ( W 1 , W 2 ) is a 2-dimensional Brownian motion under Q i T .
The following result follows from direct calculation. Proposition 4.1 Let T 2 > T 1 > 0. Under our model assumptions, the price at time t of an electricity forward contract with maturity T 1 and delivery period [T 1 , T 2 ], denoted by F t (T 1 , T 2 ), is given by the following formula:

F t (T 1 , T 2 ) = π∈Π2 1 T 2 -T 1 T2 T1 (A 1 (t, T ) + A 2 (t, T ))dT, ( 4

.14)

where

A 1 (t, T ) := {x1=m π(1) ,M π(1) } F π(1) t (T )Q π(1) T [π T = π|F W t ]P[∆ π(1) T = x 1 |∆ t ]Q[D T ≤ x 1 |F 0 t ] A 2 (t, T ) := {x 1 =m π(1) ,M π(1) ; x2=m π(2) ,M π(2) } F π(2) t (T )Q π(2) T [π T = π|F W t ]P[∆ π(1) T = x 1 , ∆ π(2) T = x 2 |∆ t ] ×Q[x 1 < D T ≤ x 1 + x 2 |F 0 t ]
where, for any π ∈ Π 2 and i = 1, 2, the conditional probabilities

Q[D T ≤ x 1 |F 0 t ] and Q[x 1 < D T ≤ x 1 +x 2 |F 0 t ]
are given by (4.12) and (4.13), and

Q π(i) T [π T = π|F W t ] = 1 -Φ(m(t)/γ(t)),
where m(t) and γ(t) are defined as follows:

m(t) = ln S π(1) t S π(2) t - 1 2 ||σ π(1) -σ π(2) || 2 -(σ π(1) -σ π(2) ) • σ π(i) (T -t) γ 2 (t) = ||σ π(1) -σ π(2) || 2 (T -t).
Proof. It suffices to combine the different formulae obtained in this section and observe that for any π ∈ Π 2 and i = 1, 2 we have

Q π(i) T [π T = π|F 0 t ] = Q π(i) T [S π(1) T ≤ S π(2) T |F W T ] = Q π(i) T [X ≤ 0|F W t ]
where X := ln(S π(1) T

/S π(2) T

). Under Q

π(i) T , X = ln S π(1) t S π(2) t + 2 j=1 (σ π(1),j -σ π(2),j )( W j T -W j t ) - 2 j=1 1 2 ((σ π(1),j ) 2 -(σ π(2),j ) 2 ) -(σ π(1),j -σ π(2),j )σ π(i),j (T -t).
Thus, conditioned to F W t , the random variable X is normal with mean m(t) and variance γ2 (t), where

m(t) = ln S π(1) t S π(2) t - 2 j=1 1 2 ((σ π(1),j ) 2 -(σ π(2),j ) 2 ) -(σ π(1),j -σ π(2),j )σ π(i),j (T -t)
and

γ 2 (t) = 2 j=1 (σ π(1),j -σ π(2),j ) 2 (T -t).
Notice that only the mean m(t) depends on π(i). Finally, we have

Q π(i) T [π T = π|F W t ] = Q π(i) T [X ≤ 0|F W t ] = Q π(i) T [(X -m(t))/γ(t) ≤ -m(t)/γ(t)|F W t ] = Φ(-m(t)/γ(t)) = 1 -Φ(m(t)/γ(t)),
where Φ is the c.d.f. of a standard gaussian random variable. The proof is complete.

Numerical results

To provide a coherent and tractable framework for numerical examples, we follow the two fuels model of the previous section and we push further the simplification.

Data choice. We test the model on the French deregulated power market. The data cover the period going from January, 1 st , 2007 to December, 31 st , 2008. For the demand process (D t ), we used the data provided by the French TSO, RTE, on its web site. 2 The hourly demand can be retrieved. The two technologies we have chosen are natural gas plants and fuel combustion turbines. They are known to frequently determine the spot price during peaking hours, since they are the most expensive ones. Moreover, a decomposition of the production is provided by RTE for each type of generation asset (nuclear, hydrolic plants, coal and gas, fuels, peak). Hence, it allowed us to deduce the residual demand addressed to gas and fuels technologies by substracting the nuclear and hydrolic production from the demand. Since these two technologies are setting the price during peaking hours, we focus our analysis on one particular hour of the day. We have chosen the 12 th hour, which is usualy the first peaking hour of the day (the next one being 19 th hour). The electricity spot and future prices are provided by Powernext. The CO 2 prices are provided by PointCarbon data. For fuel and gas prices, we used Platt's data. Gas prices are quoted in GBP and fuel prices en USD. We used the daily exchange rate to convert GBP into EUR.

Reconstruction of S 1

t and S 2 t . In our model, we need to rebuild the spot prices of the two technologies S 1 t and S 2 t . To tackle with the problem of aggregating the numerous gas and fuel power plants into only two technologies, we used the information provided by the French Ministry of Industry on electricity production costs 3 . It gives an average heat rate for each techology. We use also an average emission rate for CO 2 emissions of each technology. Furthermore, for fuel power plants production costs, one need to take into account the transportation cost from ARA zone to the location of the plants. We use an average fixed cost. Thus, we obtain the following expressions for the prices of the two technologies :

S 1 t = 101.08 • S g t + 0.49 • S co2 t S 2 t = 0.38 • S f t + 0.88 • S co2 t + 13.44
where S g , S f and S co2 denote respectively gas (A C/therm), fuel and carbon emission prices (A C/ton).

Remark 5 One can observe that the ordering between the two technologies never changes on historical data. Fuel combustion turbines are known to be more expensive than gas plants. If the prices of technologies follow the dynamics given by (2.1), the probability to have different orders π(t) ∈ Π can be positive. Nevertheless, for a reasonable choice of parameters, this probability can be made sufficiently small. Hence, we make the rough approximation that P(S 1 t < S 2 t ) = 1 for all t.

Estimation of electricity demand. The demand process given by expression (4.11) is estimated via the Maximum Likelihood Principle. Let's remind that the demand process is given by :

D t = b(t) + X t = b 0 + b 1 cos(2πt -b 2 ) + X t
where X t is an Ornstein-Uhlenbeck process with a known Likelihood expression (see [START_REF] Ait-Sahalia | Maximum likelihood estimation of discretely sampled diffusions: a closed-form approximation approach[END_REF], Section 5). For a discrete sample (D t1 , . . . , D tn ) observed at fixed times with a constant time step (t i -t i-1 ) = ∆t, i = 1 . . . n, an expression of the Likelihood is

L(b 0 , b 1 , b 2 , a, δ, D t1 , . . . , D tn ) = 1 ( √ 2πv) n exp - 1 2v n-1 i=1 (D ti+1 -b(t i+1 )) -e a∆t (D ti -b(t i )) 2 ,
where v = δ 2 e 2a∆t -1 2a and b(t) is the same as above. We maximize numerically this expression to obtain an estimation for the set of parameters. We then test the hypothesis that each parameter is null and finally obtain the set given in Table 1. The parameter b2 is not significantly different from 0 with threshold 99 %, thus it is taken to be zero.

Estimation of capacity process. For two technologies, the implementation of formula (2.2) is very simple. We define the following variables: where here D t is the sum of residual demands for the two technologies. The electricity spot price is defined by the following rule: If R 2 is positive, then we take P = S 2 , and if it is zero, P = S 1 . However, in our electricity spot market model, applying this rule to estimate the capacity processes ∆ 1 and ∆ 2 would lead to claiming that only the second technology (the most expensive one) is being used. Hence, to take into account all the complexity of the short-term bidding process involving production constraints (start-up cost, ramp constraints, minimal runtime...), we introduce a threshold ∆1 such that the price is given by the second technology althought

R 1 = min(D + t , ∆ 1 t ), R 2 = min((D t -∆ 1 t ) + , ∆ 2
R 1 = ∆1 < ∆ 1 .
Noting that the inequality on R 1 is equivalent to R 2 > (∆ 1 -∆1 ), the threshold ∆1 is obtained by solving the following program: min

(∆ 1 -∆1 ) n i=1 R P ti -S 1 ti 1 {R 2 t i ≤(∆ 1 -∆1 )} -S 2 ti 1 {R 2 t i >(∆ 1 -∆1 )} .
The function R is a risk criterion: we tested two cases, the L 1 and the L 2 norms. The absolute error (L 1 ) showed a global minimum and the quadratic error (L 2 ) showed a local minimum on a reasonable interval (very high price peaks disturb the convergence). Thus, we use the L 1 criterion to determine that the intermediate parameter ∆ 1 -∆1 equals 610 MW. Eventually, we have new values for (D t -∆ 1 t )1 {Dt>∆ 1 t } and since we know exactly when P t = S i t , for i = 1, 2, the estimation of the model on historical data is straightforward (see Figure 2). Finally, we can estimate parameters for the capacity process ∆ 1 t as D t = R 1 t + R 2 t is available. Theoretically, capacity thresholds m i and M i are structural and are known to producers. But, since they vary over time due to maintenance scheduling and weather conditions, we estimate their constant counterparts. Moreover, we had to deal with the fact that in our model ∆ 1 does take two values. Thus, we proceed in two steps. First, we filter the data to define a ∆ 1 t taking only two values. Second, we estimate the free parameters λ u 1 and λ d 1 using that filtered time series. The capacity process ∆ 1 is partially hidden, since it is observed only if D t > ∆ 1

t . Thus, we suppose that we observe data at discrete times t i , and we calibrate the capacity levels by minimizing the quadratic error between the series (∆ 1 ti 1 {Dt i >∆ 1 t i } ) i=1...n and two constant values, under the following structural constraints:

M 1 ≥ sup t∈[0,T ],Dt≤∆ 1 t D t ; m 1 ≥ inf t∈[0,T ],Dt>∆ 1 t D t .
Solving this calibration problem, we deduce the transformed serie ∆1 which takes two values:

∆ti = m 1 1 |∆t i -m1|<|∆t i -M1| + M 1 1 |∆t i -m1|≥|∆t i -M1| , i = 1 . . . n.
On that series, we estimate λ u 1 and λ d 1 by observing the series ( ∆1 ti 1 {Dt i > ∆1 t i } ) i=1...n . We denote (t k(i) ) i=1...n the subgrid of the discrete times where t k(i) is the last time before t i when we observe (∆ 1 ti ) i=1...n . Then, by the Bayes rule and the independence between D t and ∆1

t , the probability Q ∆1 ti = x|D ti > ∆1 ti , ∆1 t k(i) , for i = 1 . . . n, is given by:

Q i [x] := Q ∆1 ti = x|D ti > ∆ti , ∆1 t k(i) = P ∆1 ti = x| ∆1 t k(i) Q [D ti > x] Q D ti > ∆1 ti | ∆1 t k(i)
.

If follows that:

Q i [x] ≡ P ∆1 ti = x| ∆1 t k(i) Q [D ti > x] P ∆1 ti = M 1 | ∆1 t k(i) Q [D ti > M 1 ] + P ∆1 ti = m 1 | ∆1 t k(i) Q [D ti > m 1 ]
.

An expression of the Likelihood for the given sample is:

L(λ u 1 , λ d 1 , ∆t1 , . . . , ∆tn , D t1 , . . . , D tn ) = n i=1 Q i [x] 1 { ∆1 t i =x} (1 -Q i [x]) (1-1 { ∆1 t i =x} ) 1 {D t i > ∆1 t i } .
We maximize this expression to obtain the intensities. The values of the parameters of the capacity process are summarized in Table 2. We notice that λ u

1 > λ d 1 means that P[ ∆1 T = M 1 ] > P[ ∆1 T = m 1 ]
for a sufficiently long maturity T . A comparison with a naive econometric model. To evaluate the benefit of adding the demand and production capacity to the market model, we make a comparison between a simple econometric approach and ours. We consider the alternative linear model:

M 1 (MW) m 1 (MW) λ u 1 (y -1 ) λ d 1 (y -
P t = α 0 + α 1 S 1 t + α 2 S 2 t + t , (5.15)
where t is a Gaussian white noise. We compare the linear model (5.15) to ours where we add some free linear parameters and a Gaussian noise to facilitate the comparison:

P t = β 0 + i=1,2 β i S i t 1 {Dt∈I π t i (t)} + t .
In both cases, we estimate the parameters using a quadratic loss minimization. The Table 3 as well as Figure 3 shows that including explicitly demand and production capacity in the model, produces a better fit. Forward prices computation. Following the approximation given in Remark 5, in the case of two fuels, the expression (3.8) becomes

Price

F t (T 1 , T 2 ) = 1 T 2 -T 1 T2 T1 x1=m1,M1 P[∆ 1 T = x 1 |∆ t ] F 2 t (T ) + (F 1 t (T ) -F 2 t (T ))(Q[D T ≤ x 1 |F 0 t ]) dT. (5.16)
We do not have forward prices F i t (T ) at our disposal but only swap prices, i.e., values of 1 T2-T1 T2 T1 F i t (T )dT for delivery periods [T 1 , T 2 ]. Nevertheless, we make the approximation that:

F i t (T ) ≈ 1 T 2 -T 1 T2 T1 F i t (T )dT, T ∈ [T 1 , T 2 ] .
One may think that this approximation is quite rough for forward gas prices, since the spot market has daily granularity; but, for the prices of fuels, it is quite reasonable since spot prices take only one value per month.

We calibrate the spot price model on the former period, till June 2008, and then backtest it on future prices from July 2008 to February 2009. On that sufficiently wide interval, we focus on two assets: The two quarters ahead and three quarters ahead futures, covering Spring 2009 (April, May, June) and Summer 2009 (July, August and September). The results are illustrated on Figure 4 and Figure 5. We see that, as expected, the predicted price overestimates the real price. Indeed, we estimated the model on high peak hours of each day, which is over the mean price most of the time. However we observe strong correlation between predicted and historical prices as shown in Table 4.

Asset

Corr Calibration on forward prices. The model gives two relations between the price of power and the prices of commodities. As we estimated the parameters on spot prices, we will now do the same on forward prices. Using formula (5.16), and under the previous assumptions on prices F i t (T ), i = 1, 2, the model can be calibrated directly on forward prices. However, given the great number of parameters, we must fix some of them in order to solve the identification problem: The capacity levels M 1 and m 1 , and the parameters of the demand D t are now fixed. Thus, the probability P ∆ 1 T = x|∆ t for x = m 1 , M 1 , which is integrated on the period [T 1 , T 2 ], is the only free variable. The goal is to calibrate numerically this variable on the following expression :

E [∆F t (T 1 , T 2 )] V [∆F t (T 1 , T 2 )]
t (T 1 , T 2 ) = f 1 (λ, ∆ t , D t )F 1 t (T 1 , T 2 ) + (1 -f 1 (λ, ∆ t , D t ))F 2 t (T 1 , T 2 ) where f 1 (λ, ∆ t , D t ) = x=m1,M1 1 T 2 -T 1 T2 T1 P ∆ 1 T = x|∆ 1 t Q [D T = x|D t ] dT.
These expressions depend on ∆ t and D t via the formulae (4.12) and (4.10). Thus, f 1 (λ, ∆ t , D t ) actually depends on t in an explicit manner. For the sake of simplicity, we make a few more approximations. Indeed, calibration can be difficult because of the fact that e -(λ d 1 +λ u 1 )(T -t) is very small when T t. Hence, if T t or the parameters λ (relation (4.10)) and a (relation (4.12)) are large enough, we can make the following approximations:

P [∆ T = x|∆ t ] ∼ = lim T ↑∞ P [∆ T = x] and Q [D T > x|D t ] ∼ = lim T ↑∞ Q [D T > x].
Thus, the calibration is equivalent to a linear model estimation under constraints, whose coefficients are f 1 (λ) and 1 -f 1 (λ).

Under that approximation, we obtain P [∆ T = M 1 ] and P [∆ T = m], giving the expected failure probabilities for the cheapest technology on the delivery period [T 1 , T 2 ]. The computation gives a sound result for calibration on Summer 2009 Future price (P [∆ T = M 1 ] = 0.865), but not for Spring 2009 Future, which is clearly overestimated. We explain this drawback by the fact that we used the two most expensive technologies to price electricity.

Spot price simulations. Our model can be easily improved to obtain trajectories with high spikes. If the residual demand D t is negative, it corresponds to the case when nuclear power is the marginal unit of the system. Its cost is well-known to be constant over time ( ∼ = 15A C/MWh). On the other hand, if the residual demand D t exceeds the total capacity ∆ 1 t + ∆ 2 t of our two technologies, it corresponds to situations when electricity has to be imported. In the French market, which is a structural exporter, it corresponds to tension on the system and electricity is bought at high cost. This high cost is arbitrarily fixed to a constant value (500A C/MWh). In order to simulate the prices of commodities, we quickly estimate on our first sample of data (January 2007 to December 2008) the multivariate diffusion process given by the relation (2.1). The Figure 6 shows that this simple device allows us to get visible spikes.

Conclusion and perspectives

By building a market model for electricity and fuels, we provided a possible answer to the issue of pricing electricity-based derivatives using a risk-neutral approach, as in security markets. This model should be considered more as a methodology than as a definitive model for electricity spot and forward prices. Indeed, we think it may offer many perspectives for further developments. We see three different areas to explore. First, we assumed competitive equilibrium on the spot market; this assumption could be changed to take into account possible strategic bidding, so quantifying the possible deviation of forward electricity prices from their equilibrium due to frictions on the spot. Second, the spot market could be extended to a multizonal framework to take into account the fact that electricity is exchanged between different countries and that a spot price is formed in each country. Finally, the relation linking forward electricity prices to forward fuels prices could be extended to a wider class of contingent claims. We hope to develop these points in future papers.

where b 0

 0 , b 1 and b 2 are (positive) constants. Then we set b(t) = b 0 + b 1 cos(2πt -b 2 ). In this case, there are explicit formulae for Q

Figure 1 :

 1 Figure 1: Midday daily demand (day-ahead peakload demand from 01/01/2007 to 31/12/2008, RTE) and simulation with fitted parameters. In black line, we showed the long trend b(t).

Figure 2 :

 2 Figure 2: Midday daily prices and model fitted on historical data (POWERNEXT R day-ahead peakload prices from 01/01/2007 to 31/12/2008).

Figure 3 :

 3 Figure 3: Prices and econometric estimation of our model and a linear model (POWERNEXT R day-ahead prices from 01/01/2007 to 31/12/2008).

Figure 4 :

 4 Figure 4: Forward prices : model anticipations and market data (POWERNEXT R Future prices on peak load from 01/07/2008 to 27/02/2009, 169 obs.). Left = Spring 2009; right = Summer 2009.

F

Figure 5 :

 5 Figure 5: Forward yields : model anticipations and market data (POWERNEXT R Future yields on peak load from 01/07/2008 to 27/02/2009, 169 obs.). Up = Spring 2009; down = Summer 2009.

Figure 6 :

 6 Figure 6: Spot price simulation. Parameters calibrated on the period 01/2007 -12/2008. We use two thresholds for very high price peaks (when Dt > 8500MWh, the price is fixed to 500A C) and low demand prices (when Dt < 0MWh, the price is fixed to 15A C). The process is simulated on 780 points (3 years).

Table 1 :

 1 Parameters estimation for the demand process.

	t ),

Table 2 :

 2 Parameters for the capacity process.

	1 )

Table 3 :

 3 Model comparison. Corr := correlation with historical price; MaxE := maximum error; MAE := mean absolute error; MSE := mean square error; MPE=Mean percentage error. Errors are calculated w.r.t. historical data (POWERNEXT R day-ahead prices from 01/01/2007 to 31/12/2008).

		Corr MaxE MAE	MSE	MPE
	Linear model	0.756 406.96 18.35 919.53 23.734%
	Structural Model 0.702 385.23 17.54 786.20 23.956%

Table 4 :

 4 Model anticipations results. Corr = correlation with historical price; E = yield mean (in parenthesis the real asset value); V = yield variance; ME = maximum price error; MAE = mean absolute error; MSE = mean squared error; MPE = mean percentage error. Errors are calculated w.r.t. historical data.)

	MaxE MAE	MSE	MPE

The process (N d T -N d (T -t)-) t≥0 as the same law as (N d t ) t≥0 .

RTE: www.rte-france.fr

Ministère de l'Industrie et des Finances, www.energie.minefi.gouv.fr/energie/electric/f1e elec.htm, see "Les coûts de référence de la production électrique"

The authors thank the "Chair Finance and Sustainable Development" sponsored by EDF and Calyon for their support.


Guillaume Charrier
email: guillaume.charrier@inra.fr

Isabelle Chuine

Marc Bonhomme

Thierry Ameglio

Charrier Guillaume

Bonhomme Marc

Améglio Thierry

Assessing frost damages using dynamic models in walnut

Keywords: Frost hardiness, abiotic stress, temperature, photoperiod, dormancy, phenology

Introduction

Frost damages are of major importance for perennial crops. Across weather hazards impacting agriculture, frost induces the most important economic losses [START_REF] Snyder | Frost protection: fundamentals, practice and economics. Volume 2. Environment and Natural Resources Series -Food and Agriculture Organization of the United Nations[END_REF] e.g. damages can cost up to several hundreds of millions dollars in fruit and tree production [START_REF] Attaway | A history of Florida citrus freezes[END_REF]. Frost damages develop when the exposure (i.e. freezing event) overtakes the vulnerability (i.e. frost hardiness) of the plant [START_REF] Williams | Towards an integrated framework for assessing the vulnerability of species to climate change[END_REF]Charrier et al., 2015a).

Freezing events mainly occur when the heat energy balance is negative (i.e. radiative frost; [START_REF] Lindkvist | A frost assessment method for mountainous areas[END_REF]. Exposure to radiative frosts, which is observed under clear sky and low wind speed, exhibits temporal pattern in frequency and intensity. Although substantial freezing events could occur at different periods along the year, most studies focused on spring period when vulnerability is high (close to leaf unfolding and blooming dates (cf. [START_REF] Rigby | Spring frost risk in a changing climate[END_REF][START_REF] Kubin | Predicting spatial and temporal patterns of bud-burst and spring frost risk in north-west Europe: the implications of local adaptation to climate[END_REF][START_REF] Augspurger | Reconstructing patterns of temperature, phenology, and frost damage over 124 years: Spring damage risk is increasing[END_REF]. To assess risks of frost damages, three different approaches are commonly used:

(i) Probability of freezing event at a given date in the recorded temperature time series from meteorological stations [START_REF] Cooter | Recent frost date trends in the northeastern United States[END_REF][START_REF] Easterling | Recent changes in frost days and the frost-free season in the United States[END_REF][START_REF] Robeson | Increasing growing-season length in Illinois during the 20th century[END_REF].

(ii) Comparison between budburst dates, either observed or simulated by a phenological model and last freezing dates [START_REF] Cannell | Climatic warming, spring budburst, and frost damage on trees[END_REF][START_REF] Scheifinger | Trends of spring time frost events and phenological dates in Central Europe[END_REF][START_REF] Schwartz | Onset of spring starting earlier across the Northern Hemisphere[END_REF][START_REF] Kubin | Predicting spatial and temporal patterns of bud-burst and spring frost risk in north-west Europe: the implications of local adaptation to climate[END_REF][START_REF] Ferguson | Dynamic thermal time model of cold hardiness for dormant grapevine buds[END_REF][START_REF] Marino | Reconstruction of false spring occurrences over the southeastern United States, 1901-2007: an increasing risk of spring freeze damage?[END_REF] (iii) Estimation of frost damage using phenological models dynamically coupled to frost hardiness models [START_REF] Leinonen | A simulation model for the annual frost hardiness and freezing damage of Scots pine[END_REF][START_REF] Morin | Will tree species experience increased frost damage due to climate change because of changes in leaf phenology?[END_REF].

In the context of global change (IPCC, 2014), these different approaches led to contradictory results, independently of the approach used. Some studies forecasted an increasing vulnerability [START_REF] Hanninen | Does climatic warming increase the risk of frost damage in northern trees?[END_REF][START_REF] Inouye | Effects of climate change on phenology, frost damage, and floral abundance of montane wildflowers[END_REF]Auspurger, 2013), whereas it was supposed to decrease in other ones [START_REF] Murray | Date of budburst of fifteen species in Britain following climate warming[END_REF][START_REF] Scheifinger | Trends of spring time frost events and phenological dates in Central Europe[END_REF]Eccel et al., 5 2009;[START_REF] Baraer | Evaluation of winter-Freeze damage risk to apple trees in global warming projections[END_REF][START_REF] Kubin | Predicting spatial and temporal patterns of bud-burst and spring frost risk in north-west Europe: the implications of local adaptation to climate[END_REF][START_REF] Dai | The decreasing spring frost risks during the flowering period for woody plants in temperate area of eastern China over past 50 years[END_REF][START_REF] Morin | Will tree species experience increased frost damage due to climate change because of changes in leaf phenology?[END_REF].

Finally, some authors highlighted the fact that the pattern is complex [START_REF] Linkosalo | Predicting spring phenology and frost damage risk of Betula spp. under climatic warming: a comparison of two models[END_REF]2006;[START_REF] Gu | The 2007 eastern US spring freezes: Increased cold damage in a warming world?[END_REF] due to antagonistic effects of increasing temperature in mean and variance [START_REF] Rigby | Spring frost risk in a changing climate[END_REF].

One important missing point in most of these studies is that frost hardiness is dynamic exhibiting successive stages of hardening and dehardening (Charrier et al., 2013a), including re-hardening [START_REF] Kalberer | Deacclimation and reacclimation of cold-hardy plants: Current understanding and emerging concepts[END_REF]. Frost hardiness varies along the year and within the plant (Charrier et al., 2013b), which therefore drives the plant's frost vulnerability.

At the end of the growing season, when photoperiod, and temperature decrease, bud meristems get into endodormancy [START_REF] Welling | Molecular control of cold acclimation in trees[END_REF][START_REF] Druart | Environmental and hormonal regulation of the activity-dormancy cycle in the cambial meristem involves stage-specific modulation of transcriptional and metabolic networks[END_REF], and frost hardening is induced (Charrier et al., 2011). Endodormancy, which is a temporary suspension of meristem activity controlled by intrinsic factors, is released by chilling exposure [START_REF] Lang | Endo-, para-, and ecodormancy: physiological terminology and classification for dormancy research[END_REF]. Finally, when endodomancy is released, bud goes through the ecodormancy phase during which growth is temperature-limited. Increase of temperature in late winter-early spring promotes both bud growth and dehardening (Greer & Warrington, 1984).

Photoperiodic control of the ecodormancy stage has also been evidenced in evergreen conifers, and some deciduous species (e.g. Fagus sylvatica; Basler & Körner, 2012).

Environmental conditions control both exposure and vulnerability to freezing event, which can therefore explain the complexity in frost risk assessment. Exposure and vulnerability diverging across organs (Charrier et al., 2013b;2015), bud and shoot frost hardiness indicate both short and long term risks, respectively. Furthermore, frost hardiness being dynamic along the frost exposure period in relation with annual phenological cycle, delayed frost hardening in autumn and premature dehardening in spring can subsequently damage trees [START_REF] Bannister | Frost resistance and the distribution of conifers[END_REF]Aitken & Hannertz, 2001). An efficient process-based model for 6 frost damage should thus integrate processes, revealing indirect and direct effects of environmental factors (Hänninen & Kramer, 2007;[START_REF] Hänninen | Boreal and temperate trees in a changing climate: modelling the ecophysiology of seasonality[END_REF]:

-Indirect effects mediated by the annual phenological cycle i.e. time course of endo-and ecodormancy, bud growth and budset, -Direct effects i.e. modulation of frost hardiness depending on temperature, and eventually photoperiod.

Among existing dynamic models [START_REF] Fuchigami | A degree growth stage (degree GS) model and cold acclimation in temperate woody plants[END_REF][START_REF] Leinonen | A simulation model for the annual frost hardiness and freezing damage of Scots pine[END_REF][START_REF] Jonsson | Climate change and the effect of temperature backlashes causing frost damage in Picea abies[END_REF][START_REF] Ebel | A temperature index model to estimate long-term freeze-risk of Satsuma mandarins grown on the northern coast of the Gulf of Mexico[END_REF][START_REF] Ferguson | Dynamic thermal time model of cold hardiness for dormant grapevine buds[END_REF], we selected two which were compatible with these pre-requisites: the one from [START_REF] Leinonen | A simulation model for the annual frost hardiness and freezing damage of Scots pine[END_REF], derived from Kellomäki et al., 1992;1995, thereafter called Photothermal model), including thermal and photoperiodic control of frost hardiness, has already been used to assess frost vulnerability in future climatic conditions of several North American deciduous species [START_REF] Morin | Will tree species experience increased frost damage due to climate change because of changes in leaf phenology?[END_REF]). An alternative model from Ferguson et al. (2011, thereafter called Thermal model) predicts frost hardiness under a strict thermal control.

We hypothesized that critical periods facing frost damages could be (i) in autumn, when trees have not acclimated yet, (ii) in the core of the winter, when temperature are minimum, or (iii) in spring when trees have deacclimated. Finally, to identify such critical periods, we calibrated Photothermal and Thermal models to test the relevancy of including photoperiod to model frost vulnerability in a deciduous species. Calibration was performed using a large amount of ecophysiological data (e.g. frost hardiness measured over five years in two locations and three genotypes, chilling and forcing requirements for bud dormancy release measured over approx. twenty years), on the walnut tree, J. regia, which is the second perennial crop in France [START_REF] Rodier | Exploitations arboricoles et vergers: Les noyers sont désormais le deuxième verger français[END_REF], with an extensive leafless period. We selected three walnut genotypes that differed either in phenology or in maximal frost hardiness to perform pair-wise comparison and compare situations with similar phenological behavior, but contrasted maximal frost hardiness (Juglans regia L. Cv Franquette and Juglans regia x nigra cv 'Early'), and contrasted phenological behavior, but similar maximal frost hardiness (Juglans regia x nigra cv 'Early', and 'Late'). Finally, we compared the predicted damages at different periods in buds and branches with the temperature of the freezing event, the actual frost hardiness, the ratio of phenological stage completion, the date and also the annual yield of fruits the subsequent growing season.

Material and Methods

Models description

The computation of frost damages is based on three sub-models: (i) annual phenological cycle, (ii) frost hardiness and (iii) frost damage. For the reasons detailed above, frost hardiness was modeled according to [START_REF] Leinonen | A simulation model for the annual frost hardiness and freezing damage of Scots pine[END_REF] or [START_REF] Ferguson | Dynamic thermal time model of cold hardiness for dormant grapevine buds[END_REF] which shares similar philosophy i.e. frost hardening competence of a particular organ (H C ) depends on direct effect (i.e. temperature and, eventually, photoperiod) and indirect effect of environmental conditions (i.e. phenological stage). For the sake of clarity, details of frost hardiness models are provided in supplementary material, although they did not differ from the previously published versions.

Annual cycle

The annual cycle is divided into four different phases: endodormancy, ecodormancy, growth and lignification, which occur sequentially [START_REF] Hanninen | A framework for modelling the annual cycle of trees in boreal and temperate regions[END_REF].

Endodormancy

Endodormancy is induced by decreasing photoperiod and temperature, and released by chilling temperatures [START_REF] Lang | Endo-, para-, and ecodormancy: physiological terminology and classification for dormancy research[END_REF]. Chilling accumulation was modeled according to the inverse of the Richardson function (originally used to predict ecodormancy release; [START_REF] Richardson | A model for estimating the completion of rest for Redhaven and Elberta peach trees[END_REF] which best described endodormancy release dates in walnut trees and other fruit species [START_REF] Chuine | Can phenological models predict tree phenology accurately in the future? The unrevealed hurdle of endodormancy break[END_REF]. When CU(t) reaches parameter CU crit (arbitrary chilling units, CU), endodormancy is released and ecodormancy stage initiates.

ݐܷሺܥ + 1ሻ = ܷܥሺݐሻ + Max൫Min൫T ୦୧୦ -ߠሺݐሻ; T ୦୧୦ -T ୪୭୵ ൯; 0൯ (1)
with CU(t), the chilling unit at day t, T high , both the temperature above which CU equals 0 and the amount of CU when temperature equals T low or lower; CU being linear between T low and T high .

Ecodormancy

Ecodormancy stage was modeled according to a sigmoid function, which is one of the most relevant for this particular stage (Caffarra et al., 2011). When FU(t) reaches parameter FU crit (arbitrary forcing units, FU), budburst occurs and apparent growth initiates.

ݐܷሺܨ + 1ሻ = ܷܨሺݐሻ + ଵ ଵା ష౩ౢ౦ሺ ഇሺሻష ఱబ ሻ (2)
with FU(t), the forcing unit at day t, slp, the slope of the function at the temperature inducing half of the maximal apparent growth rate T 50 .

Growth

Growth stage was modeled with a temperature sum in Degree Days (de Réaumur, 1735).

When GDD(t) reaches GDD crit (1200 DD), growth stage stops and lignification stage initiates [START_REF] Dreyer | Comportement d'une plante pérenne soumise à des contraintes hydriques: réponses physiologiques de jeunes Noyers à des périodes de sécheresse[END_REF].

ݐሺܦܦܩ + 1ሻ = ሻݐሺܦܦܩ + Maxሺ0; ߠሺݐሻ -T ୲୦ୋ ሻ (3)
with GDD t , the growing degree day at day t, and T thG the threshold temperature for growth.

Lignification 9

Lignification stage was considered under photoperiodic control. This stage constitutes a transition to endodormancy stage, which was arbitrary set to 1 st of September [START_REF] Chuine | Can phenological models predict tree phenology accurately in the future? The unrevealed hurdle of endodormancy break[END_REF].

Frost hardiness

Frost [START_REF] Leinonen | A simulation model for the annual frost hardiness and freezing damage of Scots pine[END_REF][START_REF] Ferguson | Dynamic thermal time model of cold hardiness for dormant grapevine buds[END_REF].

Frost damage

Based on the relation between frost hardiness (FH, temperature inducing 50% relative electrolyte leakage REL, see below) and frost sensitivity (FS, slope at this temperature; Fig. 1B), potential frost damages were calculated on a daily basis. Potential frost damages (FD, in % of electrolyte leakage) induced by a freezing event (θ min < -4°C) were therefore calculated as:

ܦܨ = ଵ ଵା ಷೄሺಷಹషഇ ሻ (4)
Quantification of damages could therefore be calculated using minimal temperature on a daily basis. For each year, the maximal predicted damages over a specific period were also calculated (from September to November, December to February and March to May, respectively called autumn, winter and spring period in the subsequent parts) as well as the damages induced by the first, maximal (i.e. absolute lowest temperature), and last freezing events. The predicted damages were related to temperature, frost hardiness, date, year or ratio of phenological stage completion (i.e.

ሺ௧ሻ େ ౙ౨౪

during endo-and

ிሺ௧ሻ ౙ౨౪
eco-dormancy, respectively).

Data

Annual cycle

Parameters used for endodormancy release stage were taken from [START_REF] Chuine | Can phenological models predict tree phenology accurately in the future? The unrevealed hurdle of endodormancy break[END_REF], which was fitted on endodormancy release dates from Juglans regia cv Franquette in orchards from France from 1975 until 2013 (18 dates). Parameters driving endodormancy release was fixed similar across genotypes since chilling requirements are similar in these genotypes (Charrier et al., 2011). From this latter study, sigmoid functions describing the response to temperature during the ecodormancy phase was fit for each genotype ('Franquette', 'Early', and 'Late', called Fra, E-Hyb, and L-Hyb, respectively). Finally, GDD crit was taken from [START_REF] Dreyer | Comportement d'une plante pérenne soumise à des contraintes hydriques: réponses physiologiques de jeunes Noyers à des périodes de sécheresse[END_REF].

Frost hardiness

Frost hardiness was measured from September until budburst on one-year-old branches along the period 2007-2012 (58 to 61 dates with n=5 replicates per date; supplementary materials) in two different orchards (low : 45° 46' N 03°08' E, 340m a.s.l. and high elevation: 45°43'N 03° 01' E 880m a.s.l.). We used the electrolyte leakage method that provides accurate estimation of frost hardiness along winter, and across organs (Zhang and Willison 1987;[START_REF] Sutinen | Seasonal differences in freezing stress resistance of needles of Pinus nigra and Pinus resinosa: evaluation of the electrolyte leakage method[END_REF]Charrier & Améglio, 2011). Frost hardiness was concomitantly measured on buds and bearing branch during winter 2008-2009 (10 sampling dates) in low elevation orchard.

Samples were cut into six 5-cm long segments without buds. Buds were kept intact attached to a small piece of branch. Samples were exposed to four different freezing temperatures among this set of temperatures: -5, -10, -15, -20, -30°C, -40°C. Depending on the season either the highest or the lowest temperatures were not used. Two supplementary subsamples were exposed to control temperature (+5°C) and maximal freezing temperature (- 75°C). Frozen and thawed rate was set to 5K.h -1 . Details are provided in Charrier and Améglio (2011).

Relative electrolytic leakage (REL) was calculated as (C 1 /C 2 )*100 as described in Zhang & Willison (1987). We assumed the following relationship between REL and temperature (θ) for each sample:

( ) ( ) d + e + a = REL θ c b - 1 (5)
where parameters a and d define asymptotes of the function, and b is the slope at the inflection point c.

Frost hardiness level was estimated as the temperature of the inflection point (c) of the adjusted logistic sigmoid function Eq. ( 1) [START_REF] Repo | Estimation of standard error of impedance-estimated frost resistance[END_REF], whereas frost sensitivity is considered to be estimated by the parameter b in % Damage.°C -1 . Parameter estimation was performed by nonlinear regression using ExcelStat ver. 7.5.2.

Calibration and validation

Parameters driving the annual phenological cycle were taken from literature data and therefore not fitted. In the frost hardiness and frost damage sub-models, some parameters were fixed according to observed values (FH MinBch , FH MinBud , FH MaxBch , α, β, and γ; Table I).

However

Statistical analysis

Statistical analyses were performed using R software (R Development Core Team, 2015).

The normality of the potential damage distribution was tested using Shapiro-Wilk test. Most of the distributions of predicted damages being not normally distributed, the comparison 13 between seasons, genotypes, location, organs were performed using the Kruskall-Wallis non parametric test and using Spearman's rank correlation coefficient ρ.

Results

Parameters value estimated from data

Frost hardiness and frost damage sub-models included parameters that have been estimated from the dataset (i.e. FH MinBud , FH MinBch , FH MaxBch, α, β, and γ; Tab. 1). Frost hardiness measured for buds (FH Bud ) was strongly and significantly correlated to frost hardiness measured in the bearing branch (FH Bch ) through a logarithmic function (R² = 0.904; P<0.001; Fig. 1A). Frost hardiness data were mostly measured on branches (56 different dates). We therefore calibrated the frost hardiness model on branch data and extrapolated to buds (10 common dates) such as:

ܪܨ ௨ௗ = ܽ × ܪܨ-݊ሺܮ ሻ -ܾ (6)
with a = -14.4, and b =-28.5°C.

Frost sensitivity (FS in % Damage.°C -1 ) was related to frost hardiness in all genotypes and organs (R² = 0.517; P<0.0001; Fig. 1B) such as:

FH e γ β α ⋅ + = FS (7)
with α= 0.106, β=-0.135, and γ=-13.57 (Tab. 1).

Frost sensitivity was thus high for low frost hardiness, and almost constant for the maximum frost hardiness. Minimum frost hardiness in branch (FH MinBch = -7.7 ± 0.2°C, mean ± SE) were not significantly different similar across genotypes (P= 0.104) or locations (P=0.766). For buds, frost hardiness was only measured in Franquette and was minimum at budburst date FH MinBch = -4.4 ± 0.2°C, mean ± SE). These parameters were therefore fixed to -7.7°C and -4.4°C for branches and buds, respectively (Tab. 1). However, maximal resistances (FH MaxBch ) significantly differed between Franquette (Fra, FH MaxBch = -32.2 ± 1.1°C, mean ± SE) and hybrids (Early: FH MaxBch = -36.6 ± 0.7°C, mean ± SE and Late:

FH MaxBch = -37.6 ± 1.0°C, mean ± SE; P<0.01), but not across hybrids (Early vs Late; P=0.725). Different values were therefore fixed: -32.2, and -37.1°C, for Fra and hybrids, respectively (Tab. 1).

Model validation

The calibrated Photothermal model predicted accurate values of frost hardiness in the three genotypes on the calibration dataset (Efficiency = 0.951, 0.909, 0.948, and RMSE = 1.45, 2.64, and 2.12°C in Fra, E-Hyb, and L-Hyb, respectively; Fig. 2A; Fig. S1). On the validation dataset, accuracy was also high (Efficiency = 0.855, 0.879, 0.904, and RMSEP = 2.28, 2.91, and 2.47°C in Fra, E-Hyb and L-Hyb, respectively; Fig. 2B).

The Photothermal model split maximal frost hardiness (FH Max ) into two components: one driven by photoperiod (FH Pmax ) and another one driven by temperature (FH Tmax ). Whereas Fra and E-Hyb predominantly responded to temperature (δ = 0.517 and 0.578, respectively), in L-Hyb, frost hardiness seemed to be under higher photoperiodic control than other genotypes (δ = 0.293). Differences across genotypes were also observed in the extent of effective nyctiperiod, i.e. NL 2 -NL 1 : from 2h in E-Hyb to 4h in L-Hyb. Extent of effective temperature, i.e. T 2 -T 1, also exhibited genotypic differences: from 18.9°C (-6.9 to 12.0°C) in in Fra and E-Hyb, respectively). However, in L-Hyb, this model more accurately predicted frost hardiness in the calibration dataset (Efficiency = 0.980 and RMSE = 1.30°C), but not in the validation dataset (Efficiency = 0.891 and RMSEP = 2.63°C).

Predicting frost damages

The frost damages were, finally, predicted using the Photothermal model because, in the Thermal model, the decrease in number of parameters did not compensate the loss in accuracy in predicting frost hardiness (BIC = 383.9 and 418.3, in Photothermal and Thermal version, respectively). However, dynamics were qualitatively similar in both versions of the model (Fig. S1). Maximal potential frost damages per season predicted for the three genotypes (Fra, E-Hyb, and L-Hyb) were calculated for buds and branches, at low and high elevation in autumn, winter, and spring (Fig. 3). Approximately two times less damages were predicted in branches compared to buds: maximum potential damages were relatively moderate in branches (median ca. 10-20% and lower than 35%, even in extreme cases), and higher in buds (median ca. 10-30%, reaching levels higher than 50% during coldest winters). Significantly higher damages were predicted at high than in low elevation (P<0.001; e.g. 8.9 and 4.1% of days inducing more than 10% damage in buds from Fra, in high and low elevation orchard, respectively), and in winter than in spring and autumn, for all genotypes, organs and locations (P<0.001; e.g. 2.2, 14.6 and 2.8% of days inducing more than 10% damage in buds from Fra, in autumn, winter, and spring, respectively). In spring compared to autumn, significantly higher damages were only predicted in branches from Fra at low elevation.

Frost risk is dynamic along the year and differs across orchards. Along the year, the safety margin (i.e. difference between actual frost hardiness and minimal temperature) calculated on a daily basis was relatively wide for all genotypes and locations i.e. between 10 and 20°C along the whole year for buds and between 15 and 30°C for branches (Fig. 4). However, risky periods for buds (gray areas in Fig. 4) are expected when frost exposure (frost probability > 0) could meet frost vulnerability at least once per decade (SM 10 > -10°C). Frost risk was almost negligible in low elevation orchard except in spring for early genotypes (Fra and E-Hyb), and autumn and winter in Fra only. In high elevation orchard, frost risk was higher in October (only in Fra), December and January (Fran and E-Hyb) and in February for all genotypes.

Higher frost risk was thus predicted in high elevation for all genotypes because the maximal level of bud frost hardiness achieved during winter could be overtaken by minimal temperature, even though maximum frost hardiness observed during leafless period was tightly correlated to the absolute minimal temperature observed along winter (Fig. 5A-C).

Although branches showed much lower frost hardiness than the minimal temperature they were exposed to, buds can suffer more than 50% cell lyses (4.3% of years in Fra and 2.2% in Hybrids during the considered period). Finally, early genotypes (Fra and E-Hyb) may suffer frost damage during spring in low elevation orchard, whereas Fra may be affected by autumnal freezing events, especially at high elevation.

Among the numerous freezing events recorded along the 1970-2015 period (n= 7630), predicted frost damages were primarily and negatively correlated to freezing temperature (ρ = -0.983; P < 0.0001) through a polynomial function (R² = 0.979; Fig. 6A). Relatively weaker correlations, although significant, were observed with actual frost hardiness (ρ = -0.515; P < 0.0001; Fig. 6B) and ratio of phenological stage completion (ρ = 0.336; P < 0.0001; Fig. 6C).

The predicted damages induced by the first and last freezing events with potential injuring temperature (i.e. θ min < -4°C occurred) exhibited different correlations. The correlation with temperature of the freezing event was very similar and also highly significant (ρ = -0.939; P < 0.0001; Fig. 6D), whereas, reversely oriented correlations were observed with actual frost hardiness (ρ = 0.243; P = 0.0009; Fig. 6E) or ratio of phenological stage completion (ρ = -0.268; P = 0.0002; Fig. 6F). Also, first frosts occurred on average, ten days earlier at high (Nov. 13 th ) than in low elevation (Nov. 23 rd ; P= 0.012), and more than three weeks later for last freezing event e.g. Mar. 13 th and Apr. 5 th in low and high elevation respectively (P< 0.001). The potential frost damages generated by the first and the last freezing events were similar in other genotypes, although exhibiting lower absolute values. Finally, no significant trend were observed with years (ρ = -0.037; ρ = 0.032, for the complete set or only first and last freezing events, respectively).

Yield in nuts production, obtained in an independent orchard, exhibited a wide range over a ten years period (2005-2014) from 1 to 25 kg per tree (Fig. 7). Damages predicted during the previous autumn were lower (9.1%) than during the previous winter, or spring periods (16.2, and 9.6%, respectively). However, these autumn damages were more strongly correlated to fruit production the subsequent year (R² = 0.381), although only significant at the 10% level (P = 0.057). On the contrary, no significant correlation was observed with yield in winter (R² = 0.005; P = 0.841), nor spring (R² = 0.089; P = 0.402).

Discussion

Frost damage, under field conditions, develops when frost exposure overtakes frost hardiness. Although frost exposure is relatively easy to measure using, as a proxy, air temperature from weather station, the access to accurate frost hardiness data is rather limited.

Here, a large dataset of frost hardiness seasonality over five years in two locations was used to calibrate integrated models predicting frost damages for three cultivars of walnut trees.

Three sub-models interact, namely the annual phenological cycle, the frost hardiness and the frost damage sub-models. The frost damage sub-model was based on empirical observations of the sigmoidal function between frost hardiness and frost sensitivity, over a wide range (from -39.2 to -3.9°C). The annual phenological cycle sub-model was mostly based on experimental observations of temperature response during endodormancy, ecodormancy (Charrier et al., 2011) and growth stages [START_REF] Dreyer | Comportement d'une plante pérenne soumise à des contraintes hydriques: réponses physiologiques de jeunes Noyers à des périodes de sécheresse[END_REF]. The parameter estimates of the response function to temperature during endodormancy were statistically inferred using a unique dataset of 20 years of endodormancy release date. The onset of endodormancy was arbitrarily set to the beginning of September. Using a dynamic date instead of a constant date e.g. calculated from leaf fall model (Delpierre et al., 2007) might have delayed frost hardening and, consequently, over-estimated potential damages. However, the annual cycle sub-model has proven robust in previous dormancy and budburst simulations [START_REF] Chuine | Can phenological models predict tree phenology accurately in the future? The unrevealed hurdle of endodormancy break[END_REF].

The frost hardiness sub-model rested upon parameter optimization and comparison across two alternative models: Photothermal, using photoperiod and temperature as input variables [START_REF] Leinonen | A simulation model for the annual frost hardiness and freezing damage of Scots pine[END_REF], and Thermal, using only temperature [START_REF] Ferguson | Dynamic thermal time model of cold hardiness for dormant grapevine buds[END_REF]. The photothermal model used here was based on many experimental results on Pinus sylvestris [START_REF] Valkonen | Frost hardiness of Scots pine seedlings during dormancy[END_REF][START_REF] Repo | Modelling frost resistance of trees. Modelling to understand forest functions[END_REF][START_REF] Leinonen | A simulation model for the annual frost hardiness and freezing damage of Scots pine[END_REF][START_REF] Repo | Rehardening potential of Scots pine seedlings during dehardening[END_REF][START_REF] Kellomäki | A simulation model for the succession of the boreal forest ecosystem[END_REF]1995;[START_REF] Leinonen | A 2 nd -order dynamic model for the forst hardiness of trees[END_REF][START_REF] Leinonen | A simulation model for the annual frost hardiness and freezing damage of Scots pine[END_REF]. The Thermal model, developed on Vitis sp, shares similar philosophy i.e. environmental conditions influence frost hardiness through direct and indirect (i.e. via annual cycle) effects. These models thus seem relevant for other perennial species, such as deciduous angiosperms, because of the similarity in controlling factors and processes involved: frost hardiness and annual cycle are related (Charrier et al., 2011), hardening and dehardening are under temperature control (Charrier & Améglio, 2011) and reversible (e.g. dehardening after warm spells, rehardening after cold period: [START_REF] Kalberer | Deacclimation and reacclimation of cold-hardy plants: Current understanding and emerging concepts[END_REF][START_REF] Augspurger | Reconstructing patterns of temperature, phenology, and frost damage over 124 years: Spring damage risk is increasing[END_REF][START_REF] Arora | Increased risk of freeze damage in woody perennials visà-vis climate change: importance of deacclimation and dormancy response[END_REF], and the variability in frost hardiness is observed at the intra-individual scale (e.g. buds vs branches; Charrier et al., 2013b).

Even though Thermal model more accurately predicted frost hardiness in the calibration dataset (RMSE = 1.79, and 2.07°C, in the Thermal and Photothermal model, respectively), its predictive value was lower (RMSEP = 3.24, and 2.55°C, in the Thermal and Photothermal model, respectively). Furthermore, it does not properly simulate frost dehardening and rehardening (Fig. S1; [START_REF] Kalberer | Deacclimation and reacclimation of cold-hardy plants: Current understanding and emerging concepts[END_REF][START_REF] Arora | Increased risk of freeze damage in woody perennials visà-vis climate change: importance of deacclimation and dormancy response[END_REF]. Due to its higher number of parameters, Photothermal model more accurately predicts changes in frost hardiness. Photoperiodic control has been presented as an adaptative strategy to prevent frost damages and budburst to occur during atypically warm periods (Basler & Körner, 2012).

While the influence of photoperiod has been found relevant to predict phenological cycle in some species, including Juglans regia [START_REF] Laube | Chilling outweighs photoperiod in preventing precocious spring development[END_REF], the additive effect of photoperiod and temperature on frost hardening has not been observed during controlled experiments [START_REF] Zhang | Quantification of additive response and stationarity of frost hardiness by photoperiod and temperature in Scots pine[END_REF]. However, frost hardening has been observed under short day lengths in evergreen species (Schwartz, 1970) and, even in Juglans regia (Charrier & Améglio, 2011).

Although recent studies tend to include photoperiodic control in phenological stages (Caffarra & Donelly, 2011;Caffarra et al., 2011), the exact mechanism involved for frost hardening is unclear, even though it might involve photoreceptors in buds [START_REF] Zohner | Perception of photoperiod in individual buds of mature trees regulates leaf-out[END_REF] or cambium in species performing cortical photosynthesis such as beech trees [START_REF] Berveiller | Interspecific variability of stem photosynthesis among tree species[END_REF]). If photosystems or phytochroms are probably involved in this response, the exact signaling pathway still has to be clarified.

Parameter estimates of the Photothermal model suggest that both temperature and photoperiod do explain the variability in frost hardiness (Tab. 1). The L-Hyb genotype exhibits lower δ parameter (i.e. FH Max mainly under photoperiod control) and, consequently, a shallower slope than other genotypes in the relation between FH Max and absolute minimal temperature of the winter (Fig. 5). Furthermore, branches also exhibit steeper slopes than buds. This would imply that branches would support even lower temperature respectively), whereas buds already reached their theoretical limit respectively). High temperature dependence (δ > 0.5) could be critical under sudden temperature changes when trees are not yet hardened (autumn), or already dehardened (spring). This is especially true for Fra which presents a lower FH Max and τ parameters (i.e. more reactive to environmental changes). The slope of the relation between minimum temperature and FH Max thus highlights the plasticity of frost hardiness under changing environment and would imply a trade-off between higher safety at the margin of frost-exposed period (L-Hyb) and higher plasticity to respond to rapid temperature changes (Fra and E-Hyb). Early hybrid (E-Hyb) indeed mainly differed from late hybrids (L-Hyb) during dehardening period in relation with their faster rate of ontogenic development and therefore exhibit higher spring risks (Charrier et al., 2011).

However, higher FU critR (i.e. ability to reacclimate for a longer period after budburst), as observed in Fra, may compensate too high plasticity and prevent this genotype to endure too much late frost damages, especially under warmer springs.

According to the model, maximum frost hardiness was considered to be achieved every winter [START_REF] Larcher | Das kälteresistenzverhalten von Quercus pubescens, Ostrya carpinifolia und Fraxinus ornus auf drei thermisch unterschiedlichen standorten[END_REF]Charrier et al., 2011). This was probably the case during the considered period, but may not always happen with warmer winter temperatures, as observed, for instance, during Spanish warm autumns (Guardia et al., 2013). However, this parameter may be shunt by the photoperiodic control of frost hardening (Schwartz et al., 1970;Charrier & Améglio, 2011). This parameter should be useful to determine potential altitudinal limit (Charrier et al., 2013). In conclusion, parameter estimates for the three genotypes differed from each other (Tab. 1) and differed from that found for Pinus sylvestris. However, this highlights the necessity to use specific parameter sets for each particular species and/or genotypes (as discussed in [START_REF] Morin | Will tree species experience increased frost damage due to climate change because of changes in leaf phenology?[END_REF], especially given the potential plasticity affecting thermal prerequisites involved in frost hardiness and phenological cycles [START_REF] Kramer | Chilling and forcing requirements for foliage bud burst of European beech (Fagus sylvatica L.) differ between provenances and are phenotypically plastic[END_REF].

The Photothermal and Thermal models both accurately predict frost hardiness and sensitivity of the calibration dataset. Potential frost damages were consequently assessed using the Photothermal model, because of its better accuracy (crucial for a predictive purpose), its ability to simulate dehardening and rehardening and the sensitivity of J. regia to photoperiod for annual cycle [START_REF] Laube | Chilling outweighs photoperiod in preventing precocious spring development[END_REF], even though any causal involvement of photoperiod for frost hardiness has been proved yet.

The factors inducing frost damages are manifold, i.e. related to environmental conditions (frost exposure), to genotypic differentiation (maximal frost hardiness, sensitivity to photoperiod), and their interaction (acclimation, chilling and forcing accumulation). However, our results highlighted the major role of frost exposure (using minimum temperature as a proxy) on frost damage throughout the whole winter leafless period (ρ = -0.983; Fig. 6A), compared to frost vulnerability, using actual frost hardiness (ρ = -0.515; Fig. 6B) or phenological stage (ρ = 0.336; Fig. 6C), as proxies. These latter correlations indeed indicate that more damage are predicted when the plant is fully hardened, at the end of endodormancy, or beginning of ecodormancy. According to the hypothesis (ii), maximal damages are predicted during winter months, for colder location (Fig. 3), when minimal temperature could overtake bud frost hardiness (Fig. 45). Potential damages generated by the first and last freezing events are much lower and also depend on frost exposure (ρ = -0.939; Fig. 6D)

whereas not to vulnerability (frost hardiness: ρ = -0.243, phenological stage completion ρ = -0.268; Fig. 6E-F). Finally, frost damages are more likely to happen during autumn and winter season in colder orchards, whereas in spring in warmer orchards, especially in early genotypes. Among the parameters of this model, the difference in ∂ and τ revealed the plasticity of frost vulnerability, which could be critical after warm spells. However, FU critR may mitigate frost risks in early genotypes.

Surprisingly, even though predicted frost damages generated by first freezing events were lower than in other periods, it was significantly correlated to next growing season's fruit production (Fig. 7). Whereas, many studies focused either on winter or spring frost hardiness's (e.g. [START_REF]Spring frosts in deciduous fruit trees-morphological damage and flower hardiness[END_REF][START_REF] Kenney | Change and evolution in the plant hardiness zones of Canada[END_REF], only a few of them characterized autumnal damages which can be of similar intensities and frequencies (Charrier & Améglio, 2011;Guardia et al., 2013). Under extended warm conditions during autumn, frost damages may be critical (Guardia et al., 2013;[START_REF] Hänninen | Boreal and temperate trees in a changing climate: modelling the ecophysiology of seasonality[END_REF]. Frost hardening, which is tightly connected to chilling exposure, may be subsequently delayed until the moment of first freezing event [START_REF] Kalberer | Deacclimation and reacclimation of cold-hardy plants: Current understanding and emerging concepts[END_REF][START_REF] Arora | Increased risk of freeze damage in woody perennials visà-vis climate change: importance of deacclimation and dormancy response[END_REF]. However, the correlation between frost damage and fruit yield were relatively weak because of the numerous abiotic and biotic factors influencing fruit production, especially during the growing period (Fig. 7). In mountain plants, reproductive buds are more vulnerable than vegetative buds [START_REF] Ladinig | How endangered is sexual reproduction of high-mountain plants by summer frosts? Frost resistance, frequency of frost events and risk assessment[END_REF], and leave primordia than xylem parenchyma or cambium [START_REF] Andergassen | Frost hardiness in the juvenile and adult life phase of ivy (Hedera helix L.)[END_REF]. The exact temporal pattern of frost hardiness in flower buds would deserve further studies during autumn period using methods that could quantify frost damages in different cellular types (e.g. Infra-Red differential analysis; [START_REF] Hacker | Ice propagation in plants visualized at the tissue level by infrared differential thermal analysis (IDTA)[END_REF]Charrier et al., 2015b;in press).

At the organ scale, meristem temperature may also diverge from air temperature of several degrees [START_REF] Savvides | Meristem temperature substantially deviates from air temperature even in moderate environments: is the magnitude of this deviation species-specific?[END_REF]. Although damages on vegetative bud may be compensated by the development of vegetative latent buds, flower buds are generated the previous growing season. Resiliency of the vegetative system would therefore insure survival of the tree, if these events do not occur every year. Frost damages on buds still accounts for major losses for crop production, more than any other abiotic stress [START_REF] Rieger | Freeze protection for horticultural crops[END_REF][START_REF] Snyder | Frost protection: fundamentals, practice and economics. Volume 2. Environment and Natural Resources Series -Food and Agriculture Organization of the United Nations[END_REF], also during autumn period (Fig. 7).

Under future climate, many disturbances are expected to affect the biology of trees during winter, even though such a trend has not been observed in the conditions of the study. Global warming may displace the risk at a given location [START_REF] Kenney | Change and evolution in the plant hardiness zones of Canada[END_REF], through both exposure and vulnerability modulations. In relation with minimal temperature, frost damages should decrease if the number and intensity of freezing events decrease. Winter damages might remain similar, whereas autumn or spring period could become frost-free for a longer extent. However, many uncertainties remain about the frequency and intensity of future freezing events when microclimatic variability is not taken into account into the projections i.e. actual organ temperature (Charrier et al., 2015a). Autumn period with warmer conditions therefore deserve a particular attention in upcoming studies (Charrier & Améglio, 2011;Guardia et al., 2013;[START_REF] Hänninen | Boreal and temperate trees in a changing climate: modelling the ecophysiology of seasonality[END_REF].

Frost damage occurs at the crossing of high frost exposure and vulnerability. The modeling approach used here highlighted the major role of frost exposure rather than vulnerability on frost damages along the dormant periods for walnut trees. However, flower buds, which could endure >50% cellular damages during the harshest winters, may be more impacted by the first freezing events regarding fruit production the subsequent growing season. Charrier et al. (2011) and [START_REF] Dreyer | Comportement d'une plante pérenne soumise à des contraintes hydriques: réponses physiologiques de jeunes Noyers à des périodes de sécheresse[END_REF], for endodormancy, ecodormancy and growth stages respectively. The parameters of the frost damage sub-model were fixed according to our own observations, as well as FH MinBch , FH MinBud , and FH MaxBch (frost hardiness sub-model). Other parameters driving the frost hardiness sub-models were fitted.

Sub-model

  , in the Photothermal model, FH Max has two components i.e. FH PMax and FH TMax corresponding to photoperiodic and thermal control, respectively. The parameter δ was introduced to split FH Max into its thermal (δ) and photoperiodic (1-δ) components, for each genotype. Finally, 9 and 6 parameters were fitted in the Photothermal and the Thermal model, respectively. Optimization of these parameters was assessed by minimization of the mean squared error between simulated values and observed data from the calibration dataset. We used Generalized Reduced Gradient algorithm with different initial values of parameters to avoid local minima problems. The complete dataset was split into two datasets of similar size, one for calibration (low elevation from autumn 2007 to spring 2009 and high elevation from autumn 2008 to spring 2009: 34, 32 and 33 sampling dates for Fra, E-Hyb and L-Hyb, respectively), and the other one for validation (low elevation from autumn 2009 to spring 2012: 26, 26 and 28 sampling dates for Fra, E-Hyb and L-Hyb, respectively; Tab. S1).

Fra

  to 35.2°C (-9.3 to 25.9°C) in L-Hyb. Another difference lied in the time constant (τ= 5 days in Fra; 16 and 15 days in E-Hyb and L-Hyb, respectively).The Thermal model (under strict temperature control), more accurately predicted frost hardiness than the Photothermal model on calibration dataset (Efficiency = 0.928, and 0.929 and RMSE = 1.75 and 2.33°C in Fra and E-Hyb, respectively; Fig.S1), but accuracy was lower on the validation dataset (Efficiency = 0.760 and 0.751 and RMSEP = 2.94 and 4.16°C

Figure 1 .

 1 Figure caption:

Figure 2 .

 2 Figure 2. A) Simulated versus observed frost hardiness (°C) in branches from different walnut genotypes: J. regia cv Franquette (Fra), early (E-Hyb) and late hybrids J. regia x nigra (L-Hyb) in the calibration (A) and validation (B) datasets. Symbols and bars represent mean and standard errors from 5 replicates, respectively. Root Mean Standard Errors are indicated for each genotype and dataset.

Figure 3 .

 3 Figure 3. Distribution of maximal frost damages per season (A: autumn, W: winter and S: spring) in buds (upper panel) or branches (lower panel); in low (left side) or high elevation orchards (right side) for three genotypes of walnut trees (J. regia cv Franquette (Fra), early (E-Hyb) and late hybrids J. regia x nigra (L-Hyb)) during the 1970-2015 period. Thick bar represents the median, extremities of the box 1 st and third quartile, extremities of the dotted lines 1 st and 9 th decile, and open circles outliers. For each organ, same letters indicate nonsignificantly different distribution according to Kruskall-Wallis test (P>0.05).

Figure 4 .

 4 Figure 4. Safety margin (i.e. difference between bud frost hardiness and minimum temperature) for buds: median (green line), 10 th (red line) and 90 th centile (blue line) and branches: median only (brown line) and probability of freezing event (black line), along the year in buds in J. regia cv. Franquette (A-B), early (C-D) and late hybrids (E-F), from low (left) and high elevation orchards (right). Gray areas indicate when safety margin for buds falls below 10°C with a probability higher than 10% (SM 10 ) when probability of freezing exposure is not null (i.e. frost risk).

Figure 5 .

 5 Figure 5. Maximum frost hardiness predicted per year depending on the absolute minimum temperature observed during the same period for buds and branches in low and high elevation orchard for different walnut genotypes: J. regia cv Franquette (A), J. regia x nigra Early (B)and Late (C). Slopes did not differ across location (P = 0.139), but across organs and genotypes (P < 0.001): Branches: 0.287, 0.252, 0.108°C.°C -1 in Fra, E-Hyb and L-Hyb,

Table 1

 1

presents the values of the different parameters. Frost damage predicted at different period were finally compared to the yield measured in the Walnut research facility in Creysse (44° 52' N, 0°30' W, 45 m a.s.l.).

Climate input data

For frost risk simulation, meteorological data (daily mean and minimal temperatures, photoperiod) from the closest weather station i.e. 2 km away from the orchards were used. Data were spread from September 1970 until August 2015, from September 1970 until August 2015, and from September 1990 until August 2015 for low elevation, high elevation and research facility orchard, respectively.

Table 1 .

 1 Parameters of the different sub-models of the study. Parameters driving the annual cycle sub-model were taken from Chuine et al. (2016),

Plant, Cell & Environment

Acknowledgements:

The authors thank Prof. Christine Foyer, Prof. Keith Mott and two anonymous reviewers for their helpful comments and critical readings on the initial version of the manuscript. We would like to acknowledge Anne-Emilie Bouchardeau, Christophe Serre and Christian Bodet for technical help. Plant, Cell & Environment Predicted damage (%) Yield (kg dry nuts.tree -1 )

respectively; Buds: 0.142, 0.113, 0.045°C.°C -1 in Fra, E-Hyb and L-Hyb, respectively. R² of the linear regressions are indicated.