Yohei Hosoe
email: hosoe@kuee.kyoto-u.ac.jp

Dimitri Peaucelle
email: peaucelle@laas.fr

Static Output Feedback Stabilization of Discrete-Time Linear Systems with Stochastic Dynamics Determined by an i.i.d. Process

This paper discusses an extension of our earlier results about state feedback stabilization toward static output feedback stabilization for discrete-time linear systems with stochastic dynamics determined by an independent and identically distributed (i.i.d.) process. A completely-linearized necessary and sufficient inequality condition for static output feedback synthesis is known to be difficult to obtain even in the case with deterministic systems. Hence, this paper derives with the so-called S-variable approach a bilinear matrix inequality condition that can be effectively solved through exploiting the result of state feedback synthesis. An illustrative numerical example is also given for such synthesis.

I. INTRODUCTION

Randomness is regarded as one of the important concepts in modern science and technology, with which various kinds of phenomenon have been interpreted and described. Control technology is also no exception, and the demand for approaches that can take into account randomness behind systems is increasing more and more. To make a contribution in this direction of studies, this paper discusses static output feedback stabilization for discrete-time linear systems with dynamics determined by a stochastic process. Such systems are called discrete-time linear random dynamical systems in the field of analytical dynamics [START_REF] Arnold | Random Dynamical Systems[END_REF]. In this paper, we particularly focus on the case that the process determining the system dynamics is independent and identically distributed (i.i.d.) with respect to the discrete time, as in the earlier studies [START_REF] Koning | Compensatability and optimal compensation of systems with white parameters[END_REF], [START_REF] Hosoe | Robust stability analysis and state feedback synthesis for discrete-time systems characterized by random polytopes[END_REF], [START_REF] Hosoe | Equivalent stability notions, Lyapunov inequality, and its application in discrete-time linear systems with stochastic dynamics determined by an i.i.d. process[END_REF].

Markov jump systems [START_REF] Costa | Discrete-Time Markov Jump Linear Systems[END_REF] are one of the most famous subclasses of random dynamical systems. Although the dynamics of a standard Markov jump system is described with a finite-mode Markov chain, the earlier study [START_REF] Costa | Stochastic stability of jump discrete-time linear systems with Markov chain in a general Borel space[END_REF] succeeded in alleviating this part of assumption so that a more general stationary Markov process can be dealt with in stability analysis. Since i.i.d. processes are a special case of stationary Markov processes, one might consider our arguments would be covered by the results for Markov jump systems. However, the above results on Markov jump systems were actually obtained under an additional assumption that the supports of random coefficient matrices (depending on a Markov process) are bounded, which is needless to use in the i.i.d. process case; instead, a milder assumption is used just for defining second-moment stability (see Assumption 2 introduced later). Hence, the arguments for the above two subclasses do not cover each other, in general. Regarding the relationship of the projection-based control approach [START_REF] Koning | Compensatability and optimal compensation of systems with white parameters[END_REF] with our study for essentially the same systems, the latter can be seen as a redevelopment of the theoretical basis first developed in the former with the recent results on linear matrix inequality (LMI)-based control. A difference can be found, e.g., in the way of dealing with stability; stability is characterized by the spectral radius of some mapping in [START_REF] Koning | Compensatability and optimal compensation of systems with white parameters[END_REF] while it is by a Lyapunov inequality in our study that can be solved as a standard LMI through some pre-calculation of expectation (details will be clearer later).

In [START_REF] Hosoe | Robust stability analysis and state feedback synthesis for discrete-time systems characterized by random polytopes[END_REF], [START_REF] Hosoe | Equivalent stability notions, Lyapunov inequality, and its application in discrete-time linear systems with stochastic dynamics determined by an i.i.d. process[END_REF], we discussed state feedback stabilization as a part of main results. The purpose of this paper is to extend those results toward static output feedback stabilization (under the assumption that the output equation of a system is described with a constant C matrix). It is considered to be unusual in practice that all the system states can be directly measured, and one may have to design an output feedback controller or to use a state observer together with a state feedback controller. This paper is concerned with the former theme (i.e., output feedback synthesis). In the case of deterministic systems, complete linearization of the inequality condition for static output feedback stabilization controller synthesis is known to be difficult [START_REF] Arzelier | Robust static output feedback design with deterministic and probabilistic certificates[END_REF]. Hence, a similar difficulty finally remains also in the present arguments for stochastic systems. However, at least the same level of results as the case of deterministic systems is obtained even in such a situation through using the S-variable approach [START_REF] Ebihara | S-Variable Approach to LMI-Based Robust Control[END_REF], which are expected to be sufficiently practical as is the case with deterministic systems.

The results in this paper could be used, e.g., for taking into account the effect of randomness of packet interarrival times in control for networked systems; an associated numerical example about random aperiodic sampling (i.e., a random version of aperiodic sampling [START_REF] Hetel | Recent developments on the stability of systems with aperiodic sampling: An overview[END_REF]) was dealt with in [START_REF] Hosoe | Equivalent stability notions, Lyapunov inequality, and its application in discrete-time linear systems with stochastic dynamics determined by an i.i.d. process[END_REF]. Another prospective application would be the combination of our theory with statistical approaches for parameter estimation such as the ensemble Kalman filter [START_REF] Evensen | The ensemble Kalman filter: Theoretical formulation and practical implementation[END_REF], [START_REF] Evensen | Data Assimilation: The Ensemble Kalman Filter[END_REF]. As is well known, the estimated parameters in modeling more or less fluctuate in practice. Our framework has a potential of dealing with such fluctuations more effectively in control than deterministic approaches.

This paper is organized as follows. Section II describes the system to be dealt with in this paper, and introduces some assumptions and the definition of exponential secondmoment stability for the system. Section III briefly reviews a Lyapunov inequality characterizing the stability as well as an inequality condition for synthesis of stabilizing state feedback. Then, Section IV discusses main results about static output feedback synthesis. Since the technique of change of variables used in the earlier study on state feedback synthesis is difficult to employ in the case of output feedback directly, we consider introducing an auxiliary variable called an Svariable for deriving a new inequality condition. Although we will not obtain a completely-linearized inequality condition even with the S-variable, the new condition can be effectively solved by exploiting the result of state feedback synthesis appropriately. We also describe the details of this, and propose a two-step synthesis method using our inequality conditions. An illustrative numerical example of our synthesis is given in Section V.

We use the following notation in this paper. The set of real numbers, that of positive real numbers and that of nonnegative integers are denoted by R, R + and N 0 , respectively. The set of n-dimensional real column vectors and that of m × n real matrices are denoted by R n and R m×n , respectively. The set of n × n positive definite matrices is denoted by S n×n + . The identity matrix of size n is denoted by I n . The Euclidean norm of the vector (•) is denoted by ||(•)||. The vectorization of the matrix (•) in the row direction is denoted by row(•), i.e., row(

•) := [row 1 (•), . . . , row m (•)],
where m is the number of rows of the matrix and row i (•) denotes the ith row. The Kronecker product is denoted by ⊗. The (block) diagonal matrix is denoted by diag(•). For the real square matrix M , He(M ) := M + M T , where M T denotes the transpose of M . The expectation (i.e., the expected value) of the random variable (•) is denoted by E[(•)]; this notation is also used for the expectation of a random matrix. If s is a random variable obeying the distribution D, then we represent it as s ∼ D.

II. DISCRETE-TIME LINEAR STOCHASTIC SYSTEMS AND SECOND-MOMENT STABILITY

A. Discrete-Time Linear Systems with Stochastic Dynamics Determined by an i.i.d. Process

Let us consider the Z-dimensional discrete-time stochastic process ξ = (ξ k ) k∈N0 satisfying the following assumption.

Assumption 1: ξ k is independent and identically distributed (i.i.d.) with respect to the discrete time k ∈ N 0 .

This assumption naturally makes ξ stationary and ergodic [START_REF] Klenke | Probability Theory: A Comprehensive Course[END_REF]. The support of ξ k is denoted by Ξ. By definition, Ξ ⊂ R Z , and Ξ corresponds to the set of values that ξ k can take (at each k).

With such a process ξ, consider the discrete-time linear random dynamical system

x k+1 = A(ξ k )x k + B(ξ k )u k , (1)
y k = Cx k , (2)
where

x k ∈ R n , u k ∈ R mu , y k ∈ R my , A : Ξ → R n×n , B : Ξ → R n×mu , C ∈ R my×n
, and the initial state x 0 is assumed to be deterministic. The dynamics of this system is determined by the random coefficient matrix A(ξ k ), which is i.i.d. with respect to k because of Assumption 1. Hence, we call the system having such dynamics a discrete-time linear system with stochastic dynamics determined by an i.i.d. process. The representation of this A matrix is indeed general in the sense that any random matrix (denoted by A k ) can be represented in that form with appropriate ξ; we can always take A(•) and ξ such that A(ξ k ) = A k (under the i.i.d. assumption on A k and ξ k ). A similar comment also applies to B(ξ k ). Although C in ( 2) is confined to be a deterministic constant matrix, this is natural, e.g., when a deterministic filter is connected to the output of the system.

B. Closed-Loop System and Second-Moment Stability

As a controller for the system described in the preceding subsection, let us consider the output feedback

u k = Ky k (3)
with static time-invariant gain K ∈ R mu×my . Then, the closed-loop system becomes

x k+1 = A of (ξ k )x k , A of (ξ k ) = A(ξ k ) + B(ξ k )KC. ( 4
)
To define second-moment stability for this system, we introduce the following assumption (which is actually a minimal requirement for the stability definition, although the details are omitted). Assumption 2: The squares of elements of A(ξ k ) and B(ξ k ) are all Lebesgue integrable, i.e.,

E[A ij (ξ k ) 2 ] < ∞, E[B ij (ξ k ) 2 ] < ∞, (5)
where A ij (ξ k ) and B ij (ξ k ) represent the (i, j)-entries of A(ξ k ) and B(ξ k ), respectively. For each fixed K, this assumption ensures the squares of elements of A of (ξ k ) are also Lebesgue integrable. Under Assumptions 1 and 2, we define exponential second-moment stability (i.e., exponential mean square stability) [START_REF] Kozin | A survey of stability of stochastic systems[END_REF] as follows.

Definition 1: The system (4) with a fixed K satisfying Assumptions 1 and 2 is said to be exponentially stable in the second moment if there exist a ∈ R + and λ ∈ (0, 1) such that

E[||x k || 2 ] ≤ a||x 0 ||λ k (∀k ∈ N 0 , ∀x 0 ∈ R n ). (6)
In this paper, we discuss a method of designing K that stabilizes the closed-loop system (4) in the sense of Definition 1 under Assumptions 1 and 2.

III. EARLIER RESULTS ON STABILIZING STATE

FEEDBACK SYNTHESIS If all the system states can be measured at each time step (i.e., m y = n and C = I n in (2)), the present output feedback stabilization problem reduces to that of designing a stabilizing state feedback gain. Such a special case was already studied in [START_REF] Hosoe | Equivalent stability notions, Lyapunov inequality, and its application in discrete-time linear systems with stochastic dynamics determined by an i.i.d. process[END_REF], and this section reviews the earlier results before proceeding to the main discussions on static output feedback stabilization; those results themselves will also play an important role in the output feedback stabilization discussed later.

Let us consider the state feedback

u k = F x k (7)
with static time-invariant gain F ∈ R mu×n , and the associated closed-loop system

x k+1 = A sf (ξ k )x k , A sf (ξ k ) = A(ξ k ) + B(ξ k )F. (8)
We see that the output feedback stated in the preceding section reduces to the above state feedback under K = F if m y = n and C = I n . This implies that Assumption 2 also ensures the squares of elements of A sf (ξ k ) are Lebesgue integrable for each F , and hence, the second-moment stability can be defined for system [START_REF] Ebihara | S-Variable Approach to LMI-Based Robust Control[END_REF] in a fashion similar to Definition 1. Then, the following theorem is known to hold (see Theorems 2 and 3 in [START_REF] Hosoe | Equivalent stability notions, Lyapunov inequality, and its application in discrete-time linear systems with stochastic dynamics determined by an i.i.d. process[END_REF]), which gives a Lyapunov inequality condition for second-moment stability.

Theorem 1: Suppose the open-loop system (1) satisfies Assumptions 1 and 2. For given F ∈ R mu×n , the following two conditions are equivalent.

1) The closed-loop system ( 8) is exponentially stable in the second moment. 2) There exist P ∈ S n×n + and λ ∈ (0, 1) such that

E[λ 2 P -A sf (ξ 0 ) T P A sf (ξ 0 )] > 0.
(9) With the Lyapunov inequality (9), we can further obtain the following theorem [START_REF] Hosoe | Equivalent stability notions, Lyapunov inequality, and its application in discrete-time linear systems with stochastic dynamics determined by an i.i.d. process[END_REF], which gives an inequality condition for synthesis of stabilizing state feedback.

Theorem 2: Suppose the open-loop system (1) satisfies Assumptions 1 and 2. There exists a gain F such that the closed-loop system (8) is exponentially stable in the second moment if and only if there exist X ∈ S n×n + , Y ∈ R mu×n and λ ∈ (0, 1) satisfying

λ 2 X * Ḡ A X + Ḡ B Y X ⊗ I n > 0 (10)
for Ḡ A and Ḡ B given by

Ḡ A := [ ḠT A1 , . . . , ḠT An ] T ∈ R nn×n , (11)
Ḡ B := [ ḠT B1 , . . . , ḠT Bn ] T ∈ R nn×mu , (12)
Ḡ =: ḠA1 , . . . , ḠAn , ḠB1 , . . . , ḠBn ( ḠAi ∈ R n×n , ḠBi ∈ R n×mu (i = 1, . . . , n)) (13)
with a matrix

Ḡ ∈ R n×(n+mu)n (n ≤ (n+m u )n) satisfying ḠT Ḡ = E [row(A(ξ 0 )), row(B(ξ 0 ))] T • [row(A(ξ 0 )), row(B(ξ 0 ))] . (14)
In particular, F = Y X -1 is one such stabilizing gain. While the Lyapunov inequality (9) involved decision variables contained in the expectation operation, the decision variables in [START_REF] Evensen | The ensemble Kalman filter: Theoretical formulation and practical implementation[END_REF] are all uncontained in the expectation operation. Hence, once Ḡ A and Ḡ B are calculated, (10) can be solved as a standard LMI for each fixed λ. By minimizing λ through a bisection method with respect to λ 2 under such an LMI, we can obtain a stabilizing gain that is optimal in the sense of exponential second-moment stability; the minimal λ corresponds to that in [START_REF] Costa | Stochastic stability of jump discrete-time linear systems with Markov chain in a general Borel space[END_REF], which further corresponds to the convergence rate of E[||x k || 2 ] with respect to k.

Remark 1: In [START_REF] Hosoe | Equivalent stability notions, Lyapunov inequality, and its application in discrete-time linear systems with stochastic dynamics determined by an i.i.d. process[END_REF], the Lyapunov inequality was shown with the non-strict inequality sign, unlike [START_REF] Hetel | Recent developments on the stability of systems with aperiodic sampling: An overview[END_REF]. However, this difference does not make any problem since for fixed P , there exists λ ∈ (0, 1) satisfying ( 9) if and only if there exists λ ∈ (0, 1) satisfying its non-strict version; this can be confirmed with P > 0 and the gap between λ and 1. Since we need to deal with strict inequalities in the output feedback synthesis later, we introduced the strict version [START_REF] Hetel | Recent developments on the stability of systems with aperiodic sampling: An overview[END_REF] of Lyapunov inequality in this section.

IV. STATIC OUTPUT FEEDBACK STABILIZATION

As reviewed in the preceding section, we can easily design a stabilizing state feedback gain F by using Theorem 2. One of the implicitly used tricks in deriving Theorem 2 was the change of variables Y = F X for X = P -1 in [START_REF] Evensen | The ensemble Kalman filter: Theoretical formulation and practical implementation[END_REF], which has been used also in the case of deterministic systems. Such conventional techniques and some other key ideas associated with ( 11)-( 14) led us to Theorem 2 in the earlier study. However, that change of variables is difficult to use in the case of the present static output feedback stabilization. This section first briefly describes this issue, and then shows an inequality condition for static output feedback stabilization that is derived with a different approach.

A. Difficulty in Change of Variables

Since the squares of elements of A of (ξ k ) are Lebesgue integrable under Assumptions 1 and 2, the stability of system (4) can be characterized by the following Lyapunov inequality as in Theorem 1.

E[λ 2 P -A of (ξ 0 ) T P A of (ξ 0 )] > 0. (15)
Then, it can be shown from the arguments in [START_REF] Hosoe | Equivalent stability notions, Lyapunov inequality, and its application in discrete-time linear systems with stochastic dynamics determined by an i.i.d. process[END_REF] that

E[A(ξ 0 ) T P A(ξ 0 )] = Ḡ T A (P ⊗ I n) Ḡ A , (16)
E[A(ξ 0 ) T P B(ξ 0 )] = Ḡ T A (P ⊗ I n) Ḡ B , (17)
E[B(ξ 0 ) T P B(ξ 0 )] = Ḡ T B (P ⊗ I n) Ḡ B (18)
(with which Theorem 2 was derived); although n was fixed at (n+m u )n in [START_REF] Hosoe | Equivalent stability notions, Lyapunov inequality, and its application in discrete-time linear systems with stochastic dynamics determined by an i.i.d. process[END_REF], it may actually be arbitrary as in Theorem 2 of the present paper as long as ( 14) is satisfied. By using ( 16)-( 18), the Lyapunov inequality ( 15) can be equivalently rewritten as

λ 2 P -( Ḡ A + Ḡ B KC) T (P ⊗ I n)( Ḡ A + Ḡ B KC) > 0. ( 19
)
This leads us to the following lemma. Lemma 1: Suppose the open-loop system (1) and (2) satisfies Assumptions 1 and 2. For given K ∈ R mu×my , the closed-loop system (4) is exponentially stable in the second moment if and only if there exist P ∈ S n×n + and λ ∈ (0, 1)

satisfying (19) for Ḡ

A and Ḡ B given by ( 11)-( 14). Hence, if (19) can be linearized as in Theorem 2, we can easily design a static output feedback. However, the Schur complement technique [START_REF] Boyd | Linear Matrix Inequalities in System and Control Theory[END_REF] and the congruence transformation using diag(X, X ⊗ I n) for X = P -1 only leads us to the following inequality from (19).

λ 2 X * Ḡ A X + Ḡ B KCX X ⊗ I n > 0 (20)
Since C is generally not a nonsingular square matrix in the present output feedback control problem, we cannot apply the change of variables to the product KCX in the above inequality. This is the issue in the output feedback synthesis, which does not occur in the state feedback synthesis (see, e.g., [START_REF] Arzelier | Robust static output feedback design with deterministic and probabilistic certificates[END_REF] for the associated discussions on deterministic systems). Even if we try to solve the bilinear matrix inequality (BMI) (20) iteratively by fixing a part of decision variables (i.e., K or X) at each iteration, the way of selecting the initial value for the variables becomes a problem; we have no reasonable way of the selection for K and X in (20). Since synthesis results obtained with BMIs do generally depend on initial values, the above issue is significant. To alleviate this issue, the next subsection shows an inequality condition with what we call an S-variable [START_REF] Ebihara | S-Variable Approach to LMI-Based Robust Control[END_REF].

B. Inequality Condition Derived with S-variable Approach and Two-Step Synthesis

We first show the following lemma, which gives another inequality condition for stability of closed-loop system (4).

Lemma 2: Suppose the open-loop system (1) and ( 2) satisfies Assumptions 1 and 2. For given K ∈ R mu×my , P ∈ S n×n + and λ ∈ (0, 1), the following two conditions are equivalent.

1) The inequality (19) is satisfied.

2) There exists S ∈ R (n+mu)×mu satisfying

Ḡ A Ḡ B I n 0 T P ⊗ I n 0 0 -λ 2 P Ḡ A Ḡ B I n 0 <He S KC -I mu .
(21) Proof: 2⇒1: Let us denote the left-hand side of (21) by Q. Then, by pre-and post-multiplying [I n , C T K T ] and its transpose on (21), we have

I n KC T Q I n KC < 0, (22)
which is equivalent to (19) (regardless of S). 1⇒2: Since (19) can be rewritten as ( 22), it follows from Finsler's lemma [START_REF] Ebihara | S-Variable Approach to LMI-Based Robust Control[END_REF] that (19) holds if and only if there exists τ ∈ R + such that

Q < τ C T K T -I mu KC -I mu . ( 23
)
This is nothing but (21) with

S = τ 2 C T K T -I mu . ( 24
)
This completes the proof. The auxiliary variable S in ( 21) is called an S-variable [START_REF] Ebihara | S-Variable Approach to LMI-Based Robust Control[END_REF]. By using S-variables, various kinds of inequality condition have been derived for control of deterministic systems; roughly speaking, the S-variable approach is an approach to introducing auxiliary variables in matrix inequalities based on Finsler's lemma (or the elimination lemma). Lemma 2 implies that this approach is useful even for the present stochastic systems control (in which the expectation operation has to be dealt with appropriately). For F ∈ R mu×n and S 2 ∈ R mu×mu , let us represent the S-variable in ( 21) by

S = F T S 2 -S 2 (25)
without loss of generality; S 2 is nonsingular when ( 21) is satisfied, which can be confirmed from the lower-right block and P > 0 (to ensure this nonsingularity, the inequality sign in the Lyapunov inequality needed to be strict). Then, (21) can be rewritten as

Ḡ A Ḡ B I n 0 T P ⊗ I n 0 0 -λ 2 P Ḡ A Ḡ B I n 0 <He F T -I mu S 2 KC -S 2 . ( 26
)
Hence, this, together with Lemmas 1 and 2 and the change of variables S 1 = S 2 K, leads us to the following theorem. Theorem 3: Suppose the open-loop system (1) and ( 2) satisfies Assumptions 1 and 2. There exists a gain K such that the closed-loop system (4) is exponentially stable in the second moment if and only if there exist

F ∈ R mu×n , S 1 ∈ R mu×my , S 2 ∈ R mu×mu , P ∈ S n×n + and λ ∈ (0, 1) satisfying Ḡ A Ḡ B I n 0 T P ⊗ I n 0 0 -λ 2 P Ḡ A Ḡ B I n 0 <He F T -I mu S 1 C -S 2 (27)
for Ḡ A and Ḡ B given by ( 11)- [START_REF] Boyd | Linear Matrix Inequalities in System and Control Theory[END_REF]. In particular, S 2 is nonsingular, and K = S -1 2 S 1 is one such stabilizing gain. As we can see, the inequality condition (27) is still bilinear. However, if we fix F , the inequality becomes an LMI. The important point here is that the variable F is not a Lyapunov matrix but an auxiliary variable, and hence, the influence of fixing it is relatively limited. In addition, a reasonable way of selecting F actually exists, which is related to the state feedback synthesis as indicated by the used symbol for the variable. To see this, we show the following lemma.

Lemma 3: Suppose the open-loop system (1) and (2) satisfies Assumptions 1 and 2. For given F ∈ R mu×n , P ∈ S n×n + and λ ∈ (0, 1), there exist S 1 ∈ R mu×my and S 2 ∈ R mu×mu satisfying (27) only if (9) holds.

Proof: By pre-and post-multiplying [I n , F T ] and its transpose on (27), we have ( 9) under ( 16)-( 18). This implies that (27) can have a solution only if (9) holds.

From this lemma, it is meaningless in (27) to consider F that does not satisfy [START_REF] Hetel | Recent developments on the stability of systems with aperiodic sampling: An overview[END_REF], which in turn implies that selecting a stabilizing state feedback gain as the variable F is reasonable. Although such selection of F does not necessarily lead to a globally optimal (or near) solution, a similar idea is empirically known to be effective in the case of deterministic systems [START_REF] Arzelier | Robust static output feedback design with deterministic and probabilistic certificates[END_REF]. Hence, this paper also exploits this idea and proposes the following two-step synthesis method for static output feedback stabilization controller synthesis.

Step 1: Minimize λ such that there exist X ∈ S n×n + and Y ∈ R mu×n satisfying (10) through a bisection method with respect to λ 2 . If λ < 1, take F = Y X -1 , and proceed to the next step; if not, the problem is unsolvable by Lemma 3.

Step 2: For fixed F , minimize λ such that there exist S 1 ∈ R mu×my , S 2 ∈ R mu×mu and P ∈ S n×n + satisfying (27) through a bisection method with respect to λ 2 . If λ < 1, then K = S -1 2 S 1 is a stabilizing output feedback gain. Although we did not refer to the iterations about the products of F , S 1 and S 2 in (27) in the above method, we can try them if the minimal λ does not become sufficiently small with the fixed F .

C. About Assumption on Constant C Matrix

Before closing this section, we would like to make remarks on the inequality condition (27), which facilitates understanding of the difference between the present study on stochastic systems and those for deterministic systems. With ( 16)-( 18), we can rewrite (27) as

E A(ξ 0 ) B(ξ 0 ) I n 0 T P 0 0 -λ 2 P A(ξ 0 ) B(ξ 0 ) I n 0 <E He F T -I mu S 1 C -S 2 . ( 28
)
This is actually a straightforward representation of the inequality condition that can be first anticipated from the form of the expectation-based Lyapunov inequality [START_REF] Löfberg | YALMIP: A toolbox for modeling and optimization in MATLAB[END_REF] with knowledge of LMIs for deterministic systems. By pre-and post-multiplying [I n , C T K T ] and its transpose on (28) (as in the proof of Lemma 2), we obtain

E A of (ξ 0 ) I n T P 0 0 -λ 2 P A of (ξ 0 ) I n <E He (F T -C T K T )(S 1 C -S 2 KC) (29)
with [START_REF] Hosoe | Equivalent stability notions, Lyapunov inequality, and its application in discrete-time linear systems with stochastic dynamics determined by an i.i.d. process[END_REF]. Since the right-hand side of this inequality becomes 0 for K = S -1 2 S 1 , the inequality reduces to [START_REF] Löfberg | YALMIP: A toolbox for modeling and optimization in MATLAB[END_REF]. Although this conclusion itself was already led in Theorem 3, the above arguments provide us with more insight into the influence of the presence of ξ in our system. To see this, let us temporarily consider the case where the C matrix in ( 2) is given by C(ξ k ) (i.e., a random matrix depending on ξ k ). Then, one might consider (28) with C replaced by C(ξ 0 ) could be the corresponding inequality condition for output feedback synthesis. Unfortunately, however, the arguments for stochastic systems are not so simple, and we soon notice that the dependence of A(ξ 0 ), B(ξ 0 ) and C(ξ 0 ) becomes a problem in obtaining (29) (i.e., [START_REF] Löfberg | YALMIP: A toolbox for modeling and optimization in MATLAB[END_REF]) with C replaced by C(ξ 0 ); for random variables s 1 and s

2 , E[s 1 ]E[s 2 ] is not equal to E[s 1 s 2 ],
in general. This is actually the reason why we had to confine ourselves to the case where the C matrix is given by a constant matrix in this paper.

In the case of deterministic systems, it is known that stability of the system x k+1 = A d x k with the A matrix A d is equivalent to that of the system η k+1 = A T d η k (i.e., the transposed system). Hence, a similar idea could be used also in the case of our stochastic systems; in that case, we might be able to deal with the random C matrix C(ξ 0 ) through alternatively confining the B matrix to a constant matrix. However, the arguments for this direction is also not simple, and the expectation operation in inequality conditions causes an issue when we transform them with conventional LMI techniques; this issue is related to the non-square Ḡ

A and the extended Lyapunov matrix P ⊗ I n in (27). Resolving this issue is left for future work.

V. NUMERICAL EXAMPLE

This section illustrates our synthesis approach through a numerical example. Let us consider the two-dimensional random process ξ k = [ξ 1k , ξ 2k ] T satisfying Assumptions 1 and 2 whose distribution is given by

ξ 1k ∼ N (0, 0.2 2 ), ξ 2k ∼ U(-0.5, 0.5), (30)
where N (0, 0.2 2 ) and U (-0.5, 0.5) denote the normal distribution with mean 0 and standard deviation 0.2 and the continuous uniform distribution with minimum -0.5 and maximum 0.5, respectively. We assume that ξ 1k and ξ 2k are also independent of each other at each k. Then, let us further consider the system (1) and ( 2) with coefficients

A(ξ k ) =   1.3 + ξ 2k 0.8 + ξ 1k -0.5 0.5 0.3 + ξ 1k ξ 2k -1.2 + ξ 2 1k -0.2 0 .8 0 .6   , B(ξ k ) = 0 ξ 1k 1 T , C = 1 0 0 0 0 1 . (31)
This system can be confirmed to be unstable by Theorem 1 with F = 0. For the above system, we first design a state feedback gain based on Theorem 2, which corresponds to Step 1 of our synthesis method. The expectation in ( 14) can be calculated with MATLAB and Symbolic Math Toolbox. Since n is arbitrary, taking it as the rank of the expectation (i.e., constant matrix) is reasonable. The corresponding Ḡ (and hence, Ḡ

A and Ḡ B ) can be constructed with the singular value decomposition. In addition, LMIs can be solved with MATLAB, YALMIP [START_REF] Löfberg | YALMIP: A toolbox for modeling and optimization in MATLAB[END_REF] and SDPT3 [START_REF] Tütüncü | Solving semidefinitequadratic-linear programs using SDPT3[END_REF]. In our calculation, n in ( 11) and ( 12) was 5. Through minimizing λ with respect to [START_REF] Evensen | The ensemble Kalman filter: Theoretical formulation and practical implementation[END_REF], we obtained λ = 0.8174, and X and Y leading to

F = Y X -1 = 2.3120 0.4735 -2.2271 . ( 32
)
Since the above λ is less than 1, the closed-loop system [START_REF] Ebihara | S-Variable Approach to LMI-Based Robust Control[END_REF] with this gain is stable by Theorem 2. Since the condition in the theorem is necessary and sufficient, this value of λ is optimal in the sense of state feedback stabilization using a static gain. This further ensures that the above value is a lower bound of λ that can be achieved by static output feedback (recall Lemma 3), which we next design. The matrices Ḡ A and Ḡ B were already calculated in the above state feedback synthesis, and we can use them also for Overlays of initial value response of output feedback system generated with 100 sample paths of ξ and x 0 = [1, 0, 0] T . the output feedback synthesis based on Theorem 3, which corresponds to Step 2 of our synthesis method. We fixed F obtained above, and minimized λ with respect to (27). Then, we obtained λ = 0.8614, and S 1 and S 2 leading to K = S -1 2 S 1 = 1.8093 -1.8487 .

Since the above λ is less than 1, the closed-loop system (4) with this gain is stable by Theorem 3. Hence, our two-step synthesis succeeded in stabilization for the present example.

The initial value response of state of the closed-loop system with x 0 = [1, 0, 0] T became as in Fig. 1; this figure shows the overlays of the response obtained with 100 sample paths of ξ.

Since the result λ = 0.8614 of output feedback synthesis is relatively close to the result λ = 0.8174 of state feedback synthesis, the room for further improvement about λ (i.e., the convergence rate of E[||x k || 2 ]) in the output feedback synthesis is considered to be not large in the present example. For reference, with an additional iteration for solving (27), we obtained λ = 0.8587 with a slightly different gain.

VI. CONCLUSIONS

In this paper, we derived an inequality condition for static output feedback stabilization of discrete-time linear systems with stochastic dynamics determined by an i.i.d. process. Since the complete linearization of the inequality condition was difficult as is the case with deterministic systems, we propose a two-step synthesis method in which a state feedback gain is used as the fixed parameter in the inequality condition. A lemma was also provided for showing reasonability of such synthesis.

As a future issue, the case of the random C matrix is remaining. In addition, the extension of the results in this paper toward robust output stabilization might be possible through exploiting the idea of random polytopes [START_REF] Hug | Random polytopes[END_REF] as in the state feedback synthesis in [START_REF] Hosoe | Robust stability analysis and state feedback synthesis for discrete-time systems characterized by random polytopes[END_REF]; random polytopes can be used for representing uncertainties, e.g., in mean and variance of random matrices.

Fig. 1 .

 1 Fig. 1.

This work was partially supported by JSPS KAKENHI Grant Number 17K14700.


Benoit Hingray
email: benoit.hingray@ujf-grenoble.fr

Juliette Blanchet

Guillaume Evin

Jean-Philippe Vidal

J.-P Vidal

Uncertainty component estimates in transient climate projections Precision of estimators in a single time or time series approach

Keywords: Uncertainty Sources, Climate Projections, ANOVA, Internal Variability, Model Uncertainty, Scenario Uncertainty, Precision of Estimates

Quantifying model uncertainty and internal variability components in climate projections has been paid a great attention in the recent years. For multiple synthetic ensembles of climate projections, we compare the precision of uncertainty component estimates obtained respectively with the two Analysis of Variance (ANOVA) approaches mostly used in recent works: the popular Single Time approach (STANOVA), based on the data available for the considered projection lead time and a time series based approach (QEANOVA), which assumes quasi-ergodicity of climate outputs over the available simulation period.

We show that the precision of all uncertainty estimates is higher when more members are used, when internal variability is smaller and/or the response-touncertainty ratio is higher. QEANOVA estimates are much more precise than STANOVA ones: QEANOVA simulated confidence intervals are roughly 3 to 5 times smaller than STANOVA ones. Except for STANOVA when less than 3 members is available, the precision is rather high for total uncertainty and moderate for internal variability estimates. For model uncertainty or responseto-uncertainty ratio estimates, the precision is low for QEANOVA to very low for STANOVA. In the most unfavorable configurations (small number of members, large internal variability), large over-or underestimation of uncertainty components is thus very likely. In a number of cases, the uncertainty analysis should thus be preferentially carried out with a time series approach or with a

Introduction

A critical issue in climate change studies is the estimation of uncertainties in projections along with the contribution of the different uncertainty sources, including scenario uncertainty, the different components of model uncertainty, and internal variability (e.g. [START_REF] Hawkins | The potential to narrow uncertainty in regional climate predictions[END_REF].

Scenario uncertainty is related to the poorly known future of greenhouse gas emissions. Model uncertainty corresponds to the dispersion between the different climate responses obtained with different models for the same forcing configuration. Model uncertainty concerns Global Climate Models (GCMs) but also all subsequent models of the climate change impact modeling chain, such as Regional Downscaling Models (regional climate models and/or statistical downscaling methods) and impact models like Hydrological Models (HMs).

Internal variability first originates from the chaotic variability of the climate at the global scale (e.g. [START_REF] Räisänen | CO2-induced climate change in CMIP2 experiments: Quantification of agreement and role of internal variability[END_REF][START_REF] Deser | Uncertainty in climate change projections: the role of internal variability[END_REF]. At the local scale, it also results from the fact that very similar large scale atmospheric circulation configurations can lead to very different meteorological observations [START_REF] Braun | Internal variability of the Canadian RCMs hydrological variables at the basin scale in Quebec and Labrador[END_REF][START_REF] Lafaysse | Internal variability and model uncertainty components in a multireplicate multimodel ensemble of hydrometeorological projections : the Alpine Durance basin[END_REF].

Estimating and partitioning uncertainty in future climate projections is first intended to help evaluating the significance of estimated changes for adaptation purposes. Besides, this is intended to highlight the most important uncertainty sources. This thus allows estimating the fraction of total uncertainty that could be narrowed via scenario refinement and model improvement, and its irreducible fraction pertaining to natural variability (e.g. [START_REF] Hawkins | The potential to narrow uncertainty in projections of regional precipitation change[END_REF][START_REF] Lafaysse | Internal variability and model uncertainty components in a multireplicate multimodel ensemble of hydrometeorological projections : the Alpine Durance basin[END_REF].

Over the recent years, uncertainty has been mostly explored and partitioned based on Multiscenarios Multimodel Multimember Ensembles (MMEs) of transient climate projections. Various methods have been proposed for this, most of them based on an Analysis of Variance (ANOVA) of projections available for the specific projection lead time considered [START_REF] Hingray | Development of probability distributions for regional climate change from uncertain global mean warming and an uncertain scaling relationship[END_REF][START_REF] Yip | A simple, coherent framework for partitioning uncertainty in climate predictions[END_REF][START_REF] Bosshard | Quantifying uncertainty sources in an ensemble of hydrological climate-impact projections[END_REF][START_REF] Giuntoli | Future hydrological extremes: the uncertainty from multiple global climate and global hydrological models[END_REF][START_REF] Van Pelt | Uncertainty in the future change of extreme precipitation over the Rhine basin: the role of internal climate variability[END_REF][START_REF] Paeth | Quantifying the evidence of climate change in the light of uncertainty exemplified by the Mediterranean hot spot region[END_REF]. In this single time approach -called STANOVA in the following -, and provided multiple members are available for each modeling chain, the model uncertainty components are estimated from the dispersion between the climate responses of the different modeling chains, obtained for each chain from the multimember mean of the projections. Similarly, the internal variability component for each modeling chain is estimated from the inter-member variance of the projections.

The different uncertainty components can also be estimated from a time series approach, benefiting from the long time series now available for most climate experiments (e.g. [START_REF] Johns | Climate change under aggressive mitigation: the ensembles multi-model experiment[END_REF][START_REF] Jacob | EURO-CORDEX: new high-resolution climate change projections for European impact research[END_REF]. A quasi-ergodic assumption for climate simulations in transient climate is often used in this context [START_REF] Hingray | Partitioning internal variability and model uncertainty components in a multimodel multimember ensemble of climate projections[END_REF]. It considers that the temporal variations of the climate response of a particular simulation chain are necessarily gradual and smooth, the higher frequency variations of the time series being due to internal variability alone. It considers next that the internal variability can be assumed to remain roughly constant over the considered period or to vary as a gradual and smooth function of the climate response of the chain. These assumptions were used for instance by [START_REF] Räisänen | CO2-induced climate change in CMIP2 experiments: Quantification of agreement and role of internal variability[END_REF]; Hawkins andSutton (2009, 2011); [START_REF] Hawkins | The potential to narrow uncertainty in projections of stratospheric ozone over the 21st century[END_REF]; [START_REF] Lafaysse | Internal variability and model uncertainty components in a multireplicate multimodel ensemble of hydrometeorological projections : the Alpine Durance basin[END_REF]; [START_REF] Bracegirdle | Sources of uncertainty in projections of twenty-first century westerly wind changes over the Amundsen Sea, West Antarctica, in CMIP5 climate models[END_REF]; [START_REF] Reintges | Uncertainty in twentyfirst century projections of the Atlantic Meridional Overturning Circulation in CMIP3 and CMIP5 models[END_REF]. In this time series approach called QEANOVA in the following, the climate change response of each simulation chain and its possible evolution with time is obtained from the long term trend estimated from the time series of the chain. The variance over time of the deviation from the climate response defines the internal variability of each chain.

The climate responses of the different chains can be used to estimate the components of model uncertainty for any projection lead time, thanks again to a usual ANOVA.

In both approaches, a potentially critical step is the possibility to identify a precise estimate of the climate response of each modeling chain. A number of previous works advocated the use of multiple runs to improve this identification (e.g. [START_REF] Kendon | Robustness of future changes in local precipitation extremes[END_REF][START_REF] Deser | Uncertainty in climate change projections: the role of internal variability[END_REF][START_REF] Kew | Robust assessment of future changes in extreme precipitation over the Rhine basin using a GCM[END_REF]. Unfortunately a limited number of members have until now been usually available for most modeling chains, logically deriving from the small number of GCM runs (e.g. [START_REF] Johns | Climate change under aggressive mitigation: the ensembles multi-model experiment[END_REF][START_REF] Jacob | EURO-CORDEX: new high-resolution climate change projections for European impact research[END_REF][START_REF] Paeth | Quantifying the evidence of climate change in the light of uncertainty exemplified by the Mediterranean hot spot region[END_REF][START_REF] Reintges | Uncertainty in twentyfirst century projections of the Atlantic Meridional Overturning Circulation in CMIP3 and CMIP5 models[END_REF]. In a number of analyses, the discrimination between internal climate variability and the climate response of each modeling chain is thus expected to be quite inaccurate owing to the limited ensemble size especially when the internal variability is non-negligible compared to the chain's climate responses (e.g. [START_REF] Deser | Uncertainty in climate change projections: the role of internal variability[END_REF][START_REF] Van Pelt | Uncertainty in the future change of extreme precipitation over the Rhine basin: the role of internal climate variability[END_REF].

These inaccuracies translate to inaccuracies in estimates of uncertainty components and of derived characteristics (e.g. response to uncertainty ratio, significance of projected changes, time of emergence of climate change, etc.).

In specific analysis configurations, the precision of uncertainty component estimates is potentially low to very low, but to our knowledge this issue has never been considered.

In the following, we characterize the precision of uncertainty components estimates when obtained either with the single time or the time series approaches introduced previously. The precision is obviously expected to depend on the multimodel multimember ensemble of projections (MME) available for the analysis and on the climate variable considered for the analysis. In the following, both ANOVA approaches are applied to multiple synthetic MMEs obtained via Monte Carlo simulations. We consider a simplified configuration where MME account for two uncertainty sources: model uncertainty and internal variability. The analysis framework we follow here is rather crude but it allows comparing both ANOVA approaches for a large number of MME config-urations with different characteristics. We especially discuss how the precision depends on the number of members used in the analysis, and we characterize the gain in precision obtained with the time series based approach, owing to the larger size of the data set accounted for. We also discuss the results for different levels of internal climate variability contribution to total uncertainty variance. We first focus on the precision of estimates obtained for model uncertainty and internal variability. We next consider the precision of estimates obtained on the one hand for the mean climate response from the model ensemble and on the other hand for the response to uncertainty ratio frequently used to estimate the significance of estimated changes.

All recent analyses of uncertainty sources in climate projections make use of the moment estimators of variance components. The present study also makes use of such estimators in order to put the precision of uncertainty estimates derived from these analyses into perspective.

The statistical bases for the estimation of uncertainty components in a given MME are summarized in section 2 for both the single lead time ANOVA and the quasi-ergodic ANOVA approaches. The Monte Carlo experiment for the simulation of multiple synthetic MMEs is presented in section 3. The precision of uncertainty estimates is presented in section 4, and discussed in section 5. Section 6 concludes.

ANOVA methods

Transient ensembles of climate experiments

Let us consider a hypothetical ensemble of climate experiments where M different members are available for each of G different climate modeling chains.

For the sake of simplicity, we assume that the same number of members is available for each chain. A chain refers for instance to a given GCM and the members to the different runs available for each GCM. A chain could also refer to a given GCM/RDM combination or to a given GCM/RDM/HM model combination. Members would respectively refer to the potentially multiple realizations obtained with the different runs for each GCM/RDM combination or for each GCM/RDM/HM combination (e.g. [START_REF] Lafaysse | Internal variability and model uncertainty components in a multireplicate multimodel ensemble of hydrometeorological projections : the Alpine Durance basin[END_REF][START_REF] Vidal | Hierarchy of climate and hydrological uncertainties in transient low-flow projections[END_REF].

Note Y (g, m, t) the climate projection obtained for the m th member of chain g for time t. Y (g, m, t) can be expressed as:

Y (g, m, t) = λ(g, t) + ν(g, m, t), (1)
where λ(g, t) is the climate response of model g for year t and ν(g, m, t) is the deviation from the climate response obtained with the member m for this year as a result of internal variability (IV).

As in most climate impact studies, the uncertainty analysis can be carried out for the change variable X. In the following, we consider absolute changes, but a similar analysis could be obtained for relative changes. Absolute changes are expressed as:

X(g, m, t) = Y (g, m, t) -Y (g, m, t C ), (2)
where X(g, m, t) are the experiment outputs of the change variable for a given member m of each chain g, each year t and where t C refers to a given reference year (or n-yr reference period centered on year t C ). These change outputs can be written as:

X(g, m, t) = ϕ(g, t) + η(g, m, t), (3)
where ϕ(g, t) is the climate change response of chain g for t and where η(g, m, t)

are the residuals obtained for the m th member of chain g as a result of IV.

In the following, we assume that the climate change response can be further decomposed as:

ϕ(g, t) = µ(t) + α(g, t), (4)
where µ(t) is the mean climate change response from all available chains at projection lead time t, α(g, t) is the deviation from µ(t) for the climate change response of chain g.

By assumption, no correlation is expected between the change responses and the residuals η s. If the residuals are additionally independent and identically distributed (i.i.d.), the total uncertainty of the change variable at t is

σ 2 X (t) = s 2 α (t) + σ 2 η (t)
where s 2 α (t) is the sample variance of the α's in equation (4) and σ 2 η (t) the variance of the η s in equation (3). They respectively correspond to the model uncertainty and IV components of total uncertainty variance for X. Note that, for any given chain g, no correlation is expected between the values of the residuals ν(g, m, t) obtained from the raw variable Y in the reference and in a future climate. The IV of the change variable X is thus the sum of the IV estimates for Y obtained for the reference and the future periods respectively.

In the following, we detail the two ANOVA approaches used to estimate (1) the mean climate change response µ(t), (2) the deviations α(g, t) for the different models, (3) the total uncertainty variance σ 2 X (t), (4) its model uncertainty s 2 α (t), and (5) its IV component σ 2 η (t). We will also consider the Response-to-Uncertainty ratio R2U(t) defined as R2U(t) = µ(t)/σ X (t), the fraction F α (t) of total variance explained by model uncertainty and the fraction F η (t) of total variance explained by internal variability. Moment estimators of variance components are used here for both ANOVA approaches.

Single Time ANOVA method

The Single Time ANOVA method (STANOVA in the following) is by far the most frequently used approach for the last years. It consists in estimating an ANOVA model on the projected changes X(g, m, t) for each time step t in turn. For any given t, climate projections data available for the time steps that precede and follow t are disregarded.

The STANOVA method uses an ANOVA model with fixed effects of the form of equations ( 3) and (4). The residuals are assumed to be i.i.d. In this configuration, an unbiased estimator of the residual's variance σ 2 η at time t is the multimodel mean of the inter-member variance at t, expressed as (Montgomery, 2012):

σ 2 η,ST (t) = 1 G(M -1) G g=1 M m=1 {X(g, m, t) -ϕ ST (g, t)} 2 (5)
where ϕ ST (g, t) is an estimate of the mean climate change response for model

g at t estimated over all M members as ϕ ST (g, t) = 1 M M m=1 X(g, m, t).
The model uncertainty variance is the inter-model sample variance of the deviations α(g, t) for all g, namely s 2 [START_REF] Montgomery | Design and Analysis of Experiments[END_REF]:

α (t) = 1 G-1 G g=1 {α(g, t)} 2 . An unbiased estimator of s 2 α (t) under constraint G g=1 α ST (g, t) = 0 is
s 2 α,ST (t) = 1 G -1 G g=1 { α ST (g, t)} 2 - 1 M σ 2 η,ST (
(t) = 1 G G g=1 ϕ ST (g, t).
An unbiased estimator of the total variance for X(g, m, t) finally reads:

σ 2 X,ST (t) = s 2 α,ST (t) + σ 2 η,ST (t)
.

Quasi-Ergodic ANOVA method

A time series ANOVA approach consists first in estimating for each chain g the climate response from the raw projections Y (g, m, t). This is typically achieved by fitting a trend model to one of the members available for each chain (e.g. [START_REF] Hawkins | The potential to narrow uncertainty in regional climate predictions[END_REF][START_REF] Reintges | Uncertainty in twentyfirst century projections of the Atlantic Meridional Overturning Circulation in CMIP3 and CMIP5 models[END_REF]. In the Quasi-Ergodic ANOVA method (QEANOVA in the following) considered here, the trend model is not obtained from one given member but from all the M members Y (g, m, t)

available for the chain [START_REF] Hingray | Partitioning internal variability and model uncertainty components in a multimodel multimember ensemble of climate projections[END_REF]. The trend model is assumed to give an unbiased estimator of the climate response function λ(g, t) of the chain (the limitations linked to the choice of the trend model are discussed in section 5.2). In this case, the trend model for Y gives an unbiased estimator of the climate change response function of the chain ϕ(g, t), as:

ϕ QE (g, t) = λ QE (g, t) -λ QE (g, t C ) (7)
where λ QE (g, t) and λ QE (g, t C ) are the trend estimates of the raw projections Y for the future and reference periods respectively.

Under the quasi-ergodic assumption, the IV of chain g for the raw data Y is estimated as the variance over time of the residuals from the climate response estimated for the chain, namely {Y (g, m, t)-λQE (g, t)}. Note that all residuals from all members are used here. If the experiment period covers T time steps, the IV of chain g is thus estimated from a sample of size T × M . Assuming no correlation between residuals obtained for the reference and future periods and assuming that the IV variance for Y is constant over the entire experiment period, the IV of chain g for the change variable X is twice that of the raw variable Y [START_REF] Hingray | Partitioning internal variability and model uncertainty components in a multimodel multimember ensemble of climate projections[END_REF]).

An unbiased estimator of the IV component for X is derived from the multimodel mean of the IV estimates obtained for the G different chains respectively [START_REF] Hingray | Uncertainty components estimates in transient climate projections. bias of moment-based estimators in the single time and time series approaches[END_REF], see also Appendix A, Eq.A.3 and Eq.A.6). When for each chain g, the trend model for Y can be expressed as a linear combination of L functions of time, it reads :

σ 2 η,QE = 2 G(T M -L) G g=1 M m=1 T k=1 Y (g, m, t k ) -λQE (g, t k ) 2 . ( 8
)
When the trend model is a simple linear function of time, L = 2.

An unbiased estimator of the sample model uncertainty variance [START_REF] Hingray | Uncertainty components estimates in transient climate projections. bias of moment-based estimators in the single time and time series approaches[END_REF], see also Appendix A, Eq.A.5):

s 2 α (t) under constraint G g=1 α QE (g, t) = 0 is
s 2 α,QE (t) = 1 G -1 G g=1 { α QE (g, t)} 2 - A(t, C) M σ 2 η,QE , (9)
where α QE (g, t) is as previously an estimate of the deviation α(g, t) for model g

expressed as α QE (g, t) = ϕ QE (g, t) -µ QE (t) with µ QE (t) = 1 G G g=1 ϕ QE (g, t)
the overall mean climate change response at t. A(t, C) depends on the projection lead time and on the different functions of time used for the trend model of Y [START_REF] Hingray | Uncertainty components estimates in transient climate projections. bias of moment-based estimators in the single time and time series approaches[END_REF]. It is positive and smaller than 1. The expression of A(t, C) is given in Appendix A when the trend model is a simple linear function of time (Eq. A.7).

An unbiased estimator of total variance for X(g, m, t) finally also reads:

σ 2 X,QE (t) = s 2 α,QE (t) + σ 2 η,QE .
3 Simulations

In the following, we mainly focus on the precision of uncertainty estimators obtained with STANOVA and QEANOVA methods for the change variable X(g, m, t e ) obtained between some future prediction lead time t e and the ref-

erence period t C .
We construct via Monte Carlo simulations multiple synthetic MMEs where

each MME is composed of M × G times series of raw projections Y (g, m, t),
with M the number of members for each of the G chains (see Appendix C for details). For reasons of simplification, each MME is simulated based on the assumptions that the climate response function of each modeling chain g is a linear function of time and that the random deviations ν(g, m, t) due to internal variability come from independent and normally distributed random variables: ν(g, m, t) ∼ N (0, σ 2 ν ).

Each MME simulated for Y (g, m, t) allows deriving the corresponding time series for X(g, m, t). In the present analysis, for graphical simplification purposes, each MME for Y is constructed so that it leads, for X at time t = t e , to a prescribed value of the Response-to-Uncertainty ratio [R2U(t e )] and to a prescribed value of the fractional variance F η (t e ) due to internal variability (Appendix C).

In the following, simulations were produced over the 1960-2100 period for a 10-yr average climate variable (see illustrative MMEs in figure 1). Decadal time series were next aggregated to 20-yr time series for all analyses presented in section 4, in order to match with the usually preferred 20-yr temporal resolution in climate analyses. Y (g, m, t) time series are thus constituted from 7 values each (T = 140/20). The reference period is t C = 1990.

For the sake of simplicity, a unique number of chains is considered in the following section (G = 5). The influence of the number of chains on the results is discussed in section 5.1. Most results are finally presented for prediction lead time t e = 2050. t e will be thus omitted in most following notations.

[

FIGURE 1 HERE ]

For each MME, both STANOVA and QEANOVA models are applied, allowing estimating the ratio R(θ) = θ/θ th between the estimated value θ and the theoretical (prescribed) value θ th of different parameters θ: the grand ensemble mean response µ(t E ), the uncertainty variances

s 2 α (t E ), σ 2 η (t E ), σ 2 X (t E ), the
Response-to-Uncertainty ratio [R2U(t E )] and the contribution to total variance of internal variability and model uncertainty, F η (t E ) and F α (t E ).

We run the simulations for different combinations of prescribed R2U and F η values and for different numbers of members M . For each combination, we simulate 1,000 MMEs. We subsequently estimate, for each parameter θ, the mean E[R(θ)] and the standard deviation SD[R(θ)] of the corresponding 1,000 ratios R(θ).

Application of a QEANOVA method requires fitting a trend model to the raw climate projections Y . A parametric model is usually used for convenience.

For simplicity, we here also assume that the trend is a linear function of time for the whole transient period [t 1 , t T ] (the influence of the trend model on the results is discussed in section 5.2). The estimation is made using ordinary least squares.

Equations ( 5) and ( 8) give unbiased estimators of σ 2 η . Respectively, equations ( 6) and ( 9) give unbiased estimators of s 2 α . E[R(θ)] obtained for σ 2 η and s 2 α with both STANOVA and QEANOVA approaches are therefore expected to equal one whatever the (M, F η , R2U) configuration. This actually also applies to all other STANOVA and QEANOVA estimators (not shown).

We therefore focus on the precision of estimators from STANOVA and QEANOVA approaches, i.e. on results obtained for SD[R(θ)], denoted SD below for shortness. The smaller the standard deviation, the better the precision of the estimator. We will additionally look at the ratio SD ST /SD QE , denoting the gain in precision obtained with the QEANOVA approach with respect to the STANOVA approach. The larger this ratio, the larger the gain.

The uncertainty analysis with the STANOVA method requires at least two members for each chain. This is not the case for the QEANOVA method which can be carried out with no more than one single member for each chain.

In the following, we also present, for information, results obtained with the QEANOVA method in case of single member MMEs.

Results

In our case and whatever the (M, F η , R2U) configuration, the distribution of the R(θ) ratios obtained respectively from the 1000 synthetic data sets was always found to be roughly normal (not shown). In a first approximation, the SD value of R(θ) thus determines the probability that the R(θ) ratio is outside of any [1 ± e] interval, that is the probability that the estimated value of the considered parameter is at least 100e% smaller or 100e% greater than the theoretical value. For instance, the probability that the ratio R(θ) is outside the [0.5,1.5] interval is 1%, 32% or 62% when SD is 0.2, 0.5 or 1 respectively. When the SD value is greater than 0.5, the probability for a large underestimation (let say R(θ) < 0.5) or for a large overestimation (let say R(θ) > 1.5) is thus non negligible (i.e. p > 0.32). When SD is greater than 1, a large underestimation or overestimation of the considered parameter is very likely (i.e. p > 0.62).

4.1 Variance components, total and fractional variance SD results obtained for variance parameters (σ 2 η , s 2 α , σ 2 X , F η , F α ) are roughly independent on the R2U value prescribed for the simulation of the synthetic MMEs (not shown). They conversely depend on the value prescribed for F η and on the number of members M . We here therefore only present results via the SD response functions in the (M, F η ) domain.

For the different variance parameters and both ANOVA methods, SD values are logically higher for small M numbers. Except for internal variability, they are also higher for large F η values, highlighting the lower precision of the estimators in these configurations (figure 2). For all parameters, STANOVA SD values are also logically higher than QEANOVA SD values.

[ FIGURE 2 HERE ]

For internal variability estimates, σ 2 η , SD values do only depend on M .

For configurations with M ≤ 3, STANOVA SD values are greater than 0.4, highlighting the possibility for significant errors in σ 2 η estimates. Whatever the (M ,F η ) configuration, QEANOVA SD values are at least 3 times smaller than STANOVA SD values. QEANOVA estimates are thus expected to be much closer to the theoretical value.

For model uncertainty estimates s 2 α , the precision gain obtained with QEANOVA is at least 2.5. It is up to 4 for high F η values (e.g. > 90%) or even more when M is smaller than 5. Whatever the (M ,F η ) configuration, SD values are much larger than those obtained for internal variability. In the most unfavorable configurations (very low M and very high F η values), STANOVA SD values are noticeably higher than 6! s 2 α estimates are very likely to be largely over or underestimated. In those cases, a misestimation can also be obtained with QEANOVA but it would be somehow less critical.

Compared to s 2 α and σ 2 η , the potential error for total variance estimates, σ 2 X , is likely to be small for both methods. SD values are globally smaller than 0.2 for QEANOVA. For STANOVA, they are smaller than 0.3 except in the low M and high F η configurations. This suggests that for any given dataset, errors obtained for s 2 α and σ 2 η partly compensate. The gain in precision obtained with QEANOVA is around 2.5, whatever the (M ,F η ) configuration.

SD response functions obtained for fractional variances estimates, F η and F α , are very similar to those obtained for σ 2 η and s 2 α estimates respectively (not shown). Whatever the method, errors in F η estimates are thus not expected to be very large. For F α , large errors are conversely expected with STANOVA in the most unfavorable (M ,F η ) configurations. The gain in precision obtained with QEANOVA is again at least 2.5.

Mean response and Response-to-Uncertainty ratio

We focus here on the estimates obtained for the mean change response µ and for the Response-to-Uncertainty ratio R2U. For both ANOVA methods and both parameters, SD values are again logically higher for smaller M and/or higher F η values. SD values additionally depend on the theoretical value prescribed for the Response-to-Uncertainty ratio R2U. SD values decrease when the theoretical value for R2U increases.

Results are presented in Figure 3 for µ estimates and three theoretical R2U values (0.2, 1, and 5). They illustrate typical configurations encountered in climate analyses. A R2U value as high as 5 can be for instance typically obtained in the 2050s for large scale temperature variables (e.g. global mean temperature as in [START_REF] Hawkins | The potential to narrow uncertainty in projections of regional precipitation change[END_REF]). A R2U value as low as 0.2 can be conversely obtained for noisy systems, e.g. regional climate variables with a high natural variability (e.g. 20-yr average regional precipitation as in [START_REF] Hingray | Partitioning internal variability and model uncertainty components in a multimodel multimember ensemble of climate projections[END_REF]).

For µ estimates, STANOVA SD values obtained when R2U = 0.2 are often greater than 0.5 and even exceed 2.5 for the most unfavorable (M , F η ) configurations. When R2U = 1, SD values are much smaller. They are lower than 0.15 except for the low M and high F η configurations. For larger prescribed R2U values, SD values tend to be negligible: they are always smaller than 0.03 when R2U = 5. QEANOVA SD values present a similar dependence on M , F η and R2U. They are however again much smaller than STANOVA SD values (roughly 2.5 times smaller whatever the M , F η and R2U configuration).

[

FIGURE 3 HERE ]

For R2U estimates, the structure of SD response functions are very similar to those obtained for µ estimates (not shown). SD values are however systematically larger (from +20% for R2U=0.1 to more than +100% for R2U=5). For noisy systems and the most unfavorable (M ,F η ) configurations, large errors are therefore expected with STANOVA. Errors are conversely expected to be low to moderate with QEANOVA estimates, which are at least 2.5 times more precise.

Discussion

Dependence to the setup and size of the projection ensemble

For the sake of simplicity, we considered that MMEs are composed of 5 modelling chains. This is much less than the number of chains currently available in recent CMIP MMEs. We also considered that the same number of members is available for all simulations chains. This is obviously never the case in real MMEs which are mostly unbalanced. In real MMEs, a large number of model combinations is also typically missing when the projections are produced with chains composed of different models (e.g. not all GCM/RCM combinations are available when projections are produced with different RCMs from different GCM outputs).

Whatever the ANOVA method, estimates of the different uncertainty components are expected to be more precise when more members are available, when the number of chains is larger and/or when the number of missing model combinations is smaller. In most of the cases however, a time series approach is still expected to produce more precise estimates of uncertainty components than a single time approach. This obviously results from the larger sample size used for the estimation which allows for a more precise estimate of the climate response and of the internal variability of each simulation chain.

As an illustration, the previous analysis was repeated with MMEs composed of 10 simulation chains instead of 5 (see Supplementary Material, Fig-

ure A1

). As expected, the precision of all uncertainty components estimates is better with 10 chains (for both approaches) and QEANOVA estimates are more precise than STANOVA ones. In the present case, the gains in precision are actually roughly the same than those presented previously with 5 chains.

The trend estimation issue in the QEANOVA approach

In the present case and for convenience, time series of the synthetic MMEs are simulated based on the assumption that the climate response of each chain is a linear function of time. The QEANOVA trend function used to estimate the climate response of each chain is also assumed to be a linear function of time.

In real datasets of projections, the climate response likely evolves in a more complex way, perhaps differently in different chains. In practice, the precise form of the climate response is not known, and for small ensembles and noisy systems it may be quite difficult to guess and/or fit. If the chosen trend model were incorrect for some chain, the response function estimate could be locally biased. Incorrect trend models could then lead to significant errors in the estimation of the local effects of the chains, and additional uncertainty would arise, particularly if internal variability is small. Because of this trend estimation issue, the results of a time series approach can thus be sub-optimal.

The linear assumptions retained for practical reasons in Section 3, make the QEANOVA results likely too optimistic since the assumed trend model is the true one. Thus the precision gain with QEANOVA in Section 4 is likely to overestimate the gain to be obtained with this method for real MMEs.

For real MMEs, non-linear models can be considered for the estimation of the trend (e.g. 4rd order polynomial used in [START_REF] Hawkins | The potential to narrow uncertainty in regional climate predictions[END_REF] and in Hingray and Said (2014) for temperature changes) and the risk of choosing non relevant trend models can be reduced with dedicated statistical tests and a visual check of the fit (if any expert knowledge can be used for instance to assess the climatic relevance of the estimated trend).

The trend estimation issue in a time series ANOVA approach, which may be critical for some MME configurations, requires thus specific care. Notwithstanding, a time series approach is likely to present a significant advantage in many cases, even when the trend is not linear, especially when internal variability is large. For illustration, we considered the same set of MMEs as previously but, instead of fitting a linear trend for estimating the climate response of each simulation chain, we fitted a 3rd order polynomial function.

Results are presented in Figure A2 of Supplementary Material. As expected, whatever the configuration in the (M, F η ) domain, the precision of QEANOVA estimates is lower with the polynomial trend model, especially for noisy MMEs and for small ensembles. Indeed, with a polynomial trend, the degrees of freedom are larger and in case of noisy systems, the local trend is often dominated by internal variability which logically leads to a lower precision of the estimated climate response. Nevertheless, QEANOVA is still much more precise than STANOVA. Note that a more relevant experiment would have been to simulate synthetic MMEs with non linear trends for the true trends but such a simulation is unfortunately not convenient and would have required to make different arbitrary and non representative choices (e.g. shape and amplitude of non linearity to be retained for each chain). The experiment discussed here thus likely again underplays the trend estimation issue in real MMEs where the true trends are not linear. Fitting higher order polynomials when the true trend is linear is indeed expected to introduce little error especially when the sample is large (the fit essentially yields the linear trend), but, as mentioned previously, the opposite does not hold: fitting a linear trend when the true trend is non linear (e.g. quadratic or cubic) will always yield local biases which may be large.

Finally, another possible limitation of a time series ANOVA approach is when internal variability is assumed to be constant over time. This assumption may be not really valid as suggested by the recent work of [START_REF] Olonscheck | Consistently Estimating Internal Climate Variability from Climate Model Simulations[END_REF]. As presented by [START_REF] Hingray | Partitioning internal variability and model uncertainty components in a multimodel multimember ensemble of climate projections[END_REF] or [START_REF] Bracegirdle | Sources of uncertainty in projections of twenty-first century westerly wind changes over the Amundsen Sea, West Antarctica, in CMIP5 climate models[END_REF], this assumption is not necessarily required and can also be relaxed.

The internal variability can for instance be assumed to vary gradually from one lead time to the other [START_REF] Hingray | Partitioning internal variability and model uncertainty components in a multimodel multimember ensemble of climate projections[END_REF]. It can also be estimated from the variance of the residuals within a moving temporal window centered on the targeted prediction lead time [START_REF] Bracegirdle | Sources of uncertainty in projections of twenty-first century westerly wind changes over the Amundsen Sea, West Antarctica, in CMIP5 climate models[END_REF].

Toward a local-QEANOVA approach

The potential limitations discussed above could make the apparent higher attractiveness of a time series approach somehow reduced for some MMEs configurations (e.g. MMEs with complex trend shapes and low noise).

However, if difficult to exploit over the whole available simulation period, the QEANOVA principles can be at least used locally, in the temporal neighborhood Ω(t) of each time point t. In such a local QEANOVA approach, the climate response of any given chain g can be obtained for instance from local linear regression, in which the trend is assumed to be smooth so that the climate response λ(g, t) of the chain can be approximated by a straight line in the neighborhood of t [START_REF] Cleveland | Robust locally weighted regression and smoothing scatterplots[END_REF]. The residuals estimated from all data in Ω(t) can next be used to assess the local internal variability of the projections. When the uncertainty analysis is produced for the climate change variable, the local linear regression and the estimation of the local internal variability is done for both the reference time t c and the future lead time t e .

In this configuration, note that the internal variability no longer has to be assumed constant between both periods.

The theoretical developments of this local QEANOVA approach, given in [START_REF] Hingray | Uncertainty components estimates in transient climate projections. bias of moment-based estimators in the single time and time series approaches[END_REF] are summarized in Appendix B. For illustration, the local QEANOVA approach was also applied to the different simulation ensembles considered in the previous sections. The neighborhood Ω(t) was here chosen to be composed from the 3 consecutive time steps centered around t. In this context, the risk that the linear trend is dominated by natural variability is likely important for noisy systems and small ensembles. This configuration would have led to a risk of underestimating the natural variability but this is not the case here thanks to the unbiased estimator derived for this variable (see Appendix B).

Results of the local QEANOVA are presented for internal variability, model uncertainty and for the Response-to-Uncertainty ratio in figure 4. Results are not a surprise. For all uncertainty parameters, the precision obtained with the local QEANOVA is significantly higher than that obtained with STANOVA (confidence intervals are 2 to 3 times smaller for the local QEANOVA). The precision of the local QEANOVA is conversely smaller than that obtained with QEANOVA. These results again logically derive from the different sizes of the samples used by the three methods.

A local QEANOVA is thus a promising alternative to STANOVA. The risk for a large misestimation of uncertainty components is indeed much lower.

Despite its smaller precision, a local QEANOVA can be also a good alternative to a time series approach. When the time series of the simulation chain presents a low to moderate noise but has a complex and nonlinear climate response for instance, a local QEANOVA actually reduces the risk of having a biased estimate of the climate response, which may occur in a time series approach because of the choice of an inappropriate trend function (see section 5.2).

[ FIGURE 4 (white lines), obtained from Monte Carlo simulations similar to those carried out in section 4, highlight the much lower precision of STANOVA, even when the fractional variance due to internal variability is low. This is also clear from the temporal variations of uncertainty estimates. When only 3 members are available, the variations of STANOVA estimates are erratic and non-negligible whatever the value of F η . They are much smoother with a local QEANOVA and they completely vanish with QEANOVA. As expected also, the model uncertainty component estimated with QEANOVA starts from zero at the beginning of the century and then gradually increases. This is thus fully coherent with the expected behavior of this uncertainty component (limit to zero for the control period). This is not the case for the two other methods, especially for STANOVA in the case of small ensembles and noisy systems (e.g. significant non zero values of model uncertainty can be achieved at time t = 2000 in Figure 5, STANOVA, 2nd line).

[ FIGURE 5 HERE ]

[ FIGURE 6 HERE ] 6 Conclusion

Precision of different methods for partitioning uncertainty

Two ANOVA approaches are currently used for partitioning model uncertainty and internal variability components in ensembles of climate projections: the popular Single Time approach (STANOVA), based on the data available for the considered projection lead time and time series approaches such as the QEANOVA approach considered here. For multiple synthetic ensembles, we compared the precision of uncertainty component estimates obtained with both approaches. The main findings are as follows:

• In both cases, the precision is logically lower for noisy MMEs and/or small ensembles of members. The precision of the mean climate response and of the Response-to-Uncertainty ratio estimates are additionally lower when the theoretical (prescribed) Response-to-Uncertainty ratio is lower.

• In many configurations, QEANOVA estimates are expected to be more precise than STANOVA ones. This does not depend on the setup and size of the MME considered (unbalanced MMEs, MMEs with missing chains, MMEs with more GCMs,...). The larger sample size used for the estimation is indeed expected to produce a better estimate of the climate response and of the internal variability of each simulation chain.

• Whatever the estimation method, the temporal variations of the uncertainty estimates should be considered as a relevant indicator of the precision of estimates. Large temporal fluctuations should call for caution and should warn for the likely low relevance of the uncertainty analysis.

Under the specific setup and assumptions of our analysis framework (e.g. length of time series, number of simulations chains, linearity of trend, constancy of internal variability over time...), we find that:

• Simulated confidence intervals obtained with QEANOVA are at least 2.5 times smaller than those obtained with STANOVA.

• The precision of estimates is rather high for total uncertainty variance and for the multichain mean response (except for STANOVA for noisy MMEs and/or for MMEs with small ensembles). The two main features of any given ensemble of projections are thus likely to be rather correctly estimated in most cases.

• This is obviously not the case for the individual uncertainty components.

For model uncertainty variance for instance, the precision of estimates is low to very low for small ensembles and/or noisy MMEs. In the latter configurations, a large over-or underestimation of uncertainty components is thus very likely preventing any precise partition of uncertainty sources.

This problem is particularly critical for STANOVA.

• STANOVA actually appears to be inappropriate in a large number of configurations as a result of the limited number of members, e.g. GCM runs, usually available for most modeling chains up to now. In the case of small ensembles and noisy MMEs, improving the relevance of uncertainty estimates requires using more data. A time series analysis approach as the QEANOVA method offers such an opportunity.

For a time series approach, the trend estimation can be an issue. For MMEs where internal variability is small for instance, there is some risk to misrepresent the forced trend. In such a case, a local QEANOVA approach offers an interesting alternative as it is not subject to the same limitation as a full time series approach. The precision of a local QEANOVA is additionally always better than that of STANOVA. Note also that, similarly to a time series approach, a local QEANOVA can be also applied when only a single member is available for the chains. It can therefore make use of all climate experiments produced in the recent years. This is another advantage over STANOVA which requires multiple members of each chain.

Precision of uncertainty estimates for real ensembles of projections

Most ensembles of projections have potentially more uncertainty sources than those considered in the present work. They may actually include scenario uncertainty, different components of model uncertainty (e.g. associated respectively to GCMs, Hydrological Models), of internal variability (with potentially its large or small scale components) and, eventually, uncertainty associated to interactions between models and/or scenarios (e.g. [START_REF] Hawkins | The potential to narrow uncertainty in projections of regional precipitation change[END_REF][START_REF] Vidal | Hierarchy of climate and hydrological uncertainties in transient low-flow projections[END_REF][START_REF] Paeth | Quantifying the evidence of climate change in the light of uncertainty exemplified by the Mediterranean hot spot region[END_REF]. In all cases, and whatever the ANOVA method used for the analysis, an idea of the precision of uncertainty components should be provided along with estimates of uncertainty components.

Results presented in the present work give already some benchmark for simple ensembles of projections. They were obtained with moment estimators of variance components. Likelihood-based estimation methods such as maximum likelihood, restricted maximum likelihood or Bayesian methods could (should) be used instead. They indeed allow for the estimation of confidence or credibility intervals associated to each uncertainty component. They also present next a number of interesting properties (e.g. [START_REF] Hartley | Maximum-likelihood estimation for the mixed analysis of variance model[END_REF][START_REF] Gelman | Analysis of variancewhy it is more important than ever[END_REF]. Conversely to method-of-moment based ANOVA approaches, they can especially not produce negative estimates of variance components.

For the two ANOVA methods considered here, negative estimates are indeed obtained for model uncertainty variance when the second terms in the right hand sides of equations ( 6) or ( 9) are greater than the first terms. The probability of having a negative estimate for s 2 α is actually directly related to the precision that can be achieved with the method for this component. For our work, this probability is roughly 2.5%, 16% or 31% when the SD value obtained for s 2 α is 0.5, 1 or 2 respectively. It can be up to 50% in the most critical (M, F η , R2U) configurations.

Likelihood-based estimation methods are also expected to use the data more efficiently than moment-based estimation methods. They are for instance easily applicable for unbalanced MMEs and/or MMEs with missing chains.

The work of [START_REF] Northrop | Quantifying sources of uncertainty in projections of future climate[END_REF] and of [START_REF] Evin | Traiter l'incertitude des projections climatiques: l'exemple des conditions meteorologiques et d'enneigement dans les Alpes francaises[END_REF] give some nice illustrations of their uses in a Bayesian context. The use of likelihoodbased estimation methods is thus to be preferred to the use of moment-based estimation methods. Even with such methods however, our results strongly suggest that a full-or partial time series approach should be used instead of a standard single time approach.

Author contributions. BH designed the analysis, developed the local-QEANOVA and the synthetic simulations. JB derived the theoretical expressions for unbiased estimators of uncertainty components and wrote the appendixes. All authors contributed to write the manuscript and discuss results.

In particular, an unbiased estimator of Λ 2 g2 is

Λ 2 g2 = Λ2 g2 -σ 2 νg M -1 g V 22 , (A.4)
where V 22 is the element (2, 2) of V.

(t k ) = 1 G-1 G g=1 {α(g, t k )} 2 , is s 2 α (t k ) = s 2 α(t k ) - 12 T T -1 T + 1 t k -t K t T -t 1 2 1 G G g=1 σ 2 νg M g , where
s 2 α(t k ) = (t k -t 1 ) 2 G -1 G g=1 Λ2 g2 .
When all GCMs have the same number of runs (M ), this expression reduces to:

s 2 α (t k ) = s 2 α(t k ) - A(t k , C) M σ 2 η , (A.5)
where σ 2 η is an unbiased estimator of internal variability variance for X

σ 2 η = 2 G G g=1 σ 2 νg (A.6)
and where

A(t k , C) = 6(T -1) T (T + 1) t k -t K t T -t 1 2 (A.7)
Appendix B Estimates with a local QEANOVA approach

We here summarize the expressions of the different uncertainty estimators obtained with a local-QEANOVA approach. The full developments, similar to those presented in appendix A for the QEANOVA approach, are detailed in [START_REF] Hingray | Uncertainty components estimates in transient climate projections. bias of moment-based estimators in the single time and time series approaches[END_REF].

Model Unbiased estimators of model uncertainty and internal variability variance.

Following [START_REF] Hingray | Uncertainty components estimates in transient climate projections. bias of moment-based estimators in the single time and time series approaches[END_REF], an unbiased estimator of model uncertainty variance, i.e. of the sample variance of α(g, t e ), is

s 2 α (t e ) = s 2 α(t e ) - 4 T 1 G G g=1 σ 2 νg M g (B.2)
where T = 2T * is the total number of time steps considered in the analysis and where σ 2 νg is an unbiased estimator of σ 2 νg given by

σ 2 νg = 1 T M g -4 Mg m=1   T * j=1 Y (g, m, t j ) -λc (g, t j ) 2 (B.3) + 2T * j=T * +1 Y (g, m, t j ) -λe (g, t j ) 2 (B.4)
with λc (g, t j ) = Λg1,c + Λg2,c (t j -t c ) and λe (g, t j ) = Λg1,e + Λg2,e (t j -t e ) where Λg1,c , Λg2,c , Λg1,e and Λg2,e are the regression coefficients of the two linear models in Eq.B.1. When all GCMs have the same number of runs (M ), the expression of model uncertainty in Eq.B.2 reduces to

s 2 α (t e ) = s 2 α (t e ) - 1 M T * σ2 η (B.5)
where σ 2 η is an unbiased estimator of internal variability variance for X For each MME simulation, the deviations of the different chains, D(g), g = 1, . . . , G, are obtained from a sample of G realizations in a normal distribution.

σ 2 η = 2 G G g=1 σ 2 νg . (B.6) Appendix C Simulation of
These realizations are scaled so that their mean is zero and their variance corresponds to the prescribed value s 2 α (t e ).

  variations of uncertainty component estimates Whatever the estimation method, the uncertainty estimates are typically derived for the different prediction lead times of a given climate simulation period. The temporal behavior of these estimates is worth considering as it may be a relevant indicator of their precision. In most cases, the climate response of each simulation chain is indeed expected to evolve gradually with time. In most cases, the different uncertainty components of a given MME are thus also expected to evolve gradually with time. Similarly, uncertainty components for model uncertainty and scenario uncertainty should start from zero when the climate variable considered for the analysis is a climate change variable. The climate change response of any simulation chain is indeed expected to start from zero when the prediction lead time corresponds to the control period. Large temporal fluctuations of uncertainty estimates and/or large values of model uncertainty for the control period should warn for the possible lack of robustness of the analysis. The possibility to consider the temporal variations of uncertainty estimates to assess their precision is illustrated for three synthetic data sets in Figure 5 for s 2 α and σ 2 η and in Figure 6 for F η . The three MMEs have been simulated following the simulation process described in section 3. Confidence intervals (colored bounds) of uncertainty estimates

FFig. 1 Fig. 2 Fig. 3 Fig. 4 Fig. 5

 12345 Fig. 1 Three synthetic multimodel multimember ensembles of time series projections (MMEs) obtained via simulation with different Fη ratios and M numbers (left: Fη = 0.1, M = 3; middle Fη = 0.5, M = 3; right: Fη = 0.5, M = 9). Time series over the 1960-2100 period for the raw projections Y (g, m, t) of a 10-yr average climate variable.Cases on the left and in the middle are illustrative of typical configurations encountered in climate analyses: "small" number of members available for each GCM (M =3) with relatively low (resp. high) contribution of IV to total variance (case 1: Fη = 0.1; case 2: Fη = 0.5). The case on the right is similar as the middle case but in a configuration with a (unrealistic) "large" number of members (M = 9).

  Considering t 1 , . . . , t T regularly spaced on [t s ; t f ], we have V 22 = 12(T -1)/{T (T + 1)(t T -t 1 ) 2 }.

	Unbiased estimation of the sample variance of the α's in the change vari-
	able. Given (A.2) and (A.4), an unbiased estimator of the sample variance of
	α(g, t k ), i.e. of s 2 α

  . A regression model is still considered to estimate the response function λ(g, t) for Y in Eq.A.1 but λ(g, t) is assumed to be only locally linear in time, in the neighborhoods of t c and t e respectively, i.e. on [t c -ω, t c + ω] and [t e -ω, t e + ω], where t e ∈ [t c , t f ] is the future prediction lead time under consideration. λ(g, t) can thus be expressed asλ(g, t) = λ c (g, t) = Λ g1,c + (tt c )Λ g2,c for t c -ω ≤ t ≤ t c + ω, λ e (g, t) = Λ g1,e+ (tt e )Λ g2,e for t e -ω ≤ t ≤ t e + ω. (B.1) The change variable X(g, m, t), for t e -ω ≤ t ≤ t e + ω, in Eq.A.2 is such that α(g, t) = (Λ g1,e -Λ g1,c ) + (tt e )Λ g2,e . Each interval [t c -ω, t c +ω] and [t e -ω, t e +ω] is discretized into T regular periods of length dt = 2ω/(T * -1), with T odd, giving respectively the sequences t 1 , . . . , t T * and t T * +1 , . . . , t 2T * . The values of Y for these different times are further considered to estimate the regression coefficients of the linear trend models in Eq.B.1. For the illustration given in section 5.3, dt = ω = 20yrs and T * = 3.

  MMEsEach MME is simulated for Y (g, m, t) assuming that, for all t s ≤ t ≤ t f : where the climate response function λ(g, t) is a linear function of time and where the ν(g, m, t) are independent and homoscedastic random variables (with variance σ 2 νg ). For convenience we further assume that λ(g, t) can be decomposed as λ(g, t) = w(t) + d(g, t), for t = t 1 , . . . , T where the mean climate response w(t) of the G chains and the deviations d(g, t) of chain g are linear functions of time, expressed as: w(t) = B + P.(tt 1 )/(t e -t C ) and = 1, . . . , G) lead, for the change variable X at time t = t e , to a prescribed value of the Response-to-Uncertainty ratio [R2U(t e )] and to a prescribed value of the fractional variance F η (t e ) due to internal variability. For X, we have by definition ϕ(g, t) = λ(g, t) -λ(g, t C ) and η(g, m, t) = ν(g, m, t) -ν(g, m, t C ). We have thus for the change variable µ(t) = w(t)w(t C ), α(g, t) = d(g, t) -d(g, t C ) and in turn µ(t) = P.(tt C )/(t e -t C ) and α(g, t) = D(g).(tt C )/(t e -t C ). The theoretical values for µ(t e ) and s 2 α (t e ) are thus as follows : µ(t e ) = P and s 2 α (t e ) = Var(D(g)). Fixing P to 1, we thus simply require in turn

	d(g, t) = D(g).(t -t 1 )/(t e -t C ) with the constraint	G g=1 D(g) = 0.
	For graphical simplification purposes, each MME for Y is constructed so
	that the parameters (σ 2 ν , P, D(g), g σ 2 X (t e ) =	1 [R2U (t e )] 2 ;	(C.2)
	σ 2 ν (t e ) =	1 2	σ 2 η (t e ) =	1 2	F η (t e ).σ 2 X (t e );	(C.3)
	Y (g, m, t) = λ(g, t) + ν(g, m, t),	(C.1)

Var(D(g)) = s 2 α (t e ) = (1 -F η (t e )).σ 2 X (t e ). (

C

.4)

Acknowledgements We thank the three anonymous reviewers for their constructive suggestions which helped to significantly improve the content of our manuscript.

Appendix A QEANOVA Estimates of uncertainty components

We summarize the theoretical developments of the QEANOVA approach to achieve unbiased estimators of uncertainty components for a simplified configuration where the trend model is a simple linear function of time. The full developments are given in [START_REF] Hingray | Uncertainty components estimates in transient climate projections. bias of moment-based estimators in the single time and time series approaches[END_REF] for the general configuration where the trend model is a linear combination of L functions of time.

Here, we first consider the case where the number of members differs from one chain to the other; the simplified equations obtained when all chains have the same number are then given. For the sake of conciseness, we omit the subscript "QE" related to the QEANOVA approach.

Model. We first consider the raw projections Y (g, m, t) with M g members for each of the G chains, assuming that, for all t s ≤ t ≤ t f :

where λ(g, t) is the trend model expressed as a linear function of time:

where t c is the reference period. We have thus α(g, t) = Λ g2 (tt c ) and

Unbiased estimation of the model parameters for the raw variable Y . We discretize [t s , t f ] into T time steps (from t s = t 1 to t f = t T ) and write t c as the Kth time step (i.e. t c = t K ). We are interested in the future prediction lead time t k ∈ [t K , t T ]. Let consider the regression model (A.1) for a particular g. Unbiased estimators of the regression parameters (Λ g1 , Λ g2 ) are given by the least square estimates

where "'" denotes the transpose , R is the T × 2 matrix of covariates whose kth row is (1, t k -t 1 ), for 1 ≤ k ≤ T , and V = (R R) -1 . Covariance matrix of the estimators ( Λg1 , Λg2 ) is given by σ 2 νg M -1 g V where σ 2 νg is an unbiased estimator of σ 2 νg given by

where L = 2 and λQE (g, t k ) = Λg1 -Λg2 (t k -t 1 ).