Yasmine Briefs
email: ybriefs@mpi-inf.mpg.de

Hendrik Leidinger
email: hleiding@mpi-inf.mpg.de

Christoph Weidenbach
email: weidenbach@mpi-inf.mpg.de

KBO Constraint Solving Revisited

Keywords: KBO Constraint Solving, NP-complete problem, Weight Ordering Constraint Solving

KBO constraint solving is very well-known to be an NPcomplete problem. Motivated by the needs of the family of SCL calculi, we consider the particular case where all terms occurring in a constraint are bound by a (single) ground term. We show that this problem and variants of this problem remain NP-complete even if the form of atoms in the constraint is further restricted. In addition, for a non-strict, partial term ordering solely based on symbol counting constraint solving remains NP-complete. Nevertheless, we provide a new simple algorithm testing KBO constraint solvability that performs well on benchmark examples.

Introduction

The family of SCL calculi (Clause Learning from Simple Models) [START_REF] Bromberger | Deciding the Bernays-Schoenfinkel fragment over bounded difference constraints by simple clause learning over theories[END_REF][START_REF] Bromberger | SCL(FOL) revisited[END_REF]13] perform reasoning on a set of first-order clauses. They develop a trail of ground literals with respect to a ground term (atom) bound β and an ordering ≺. All ground literals on the trail are ≺ (or ) smaller than the ground term (atom) β and ≺ should in particular have the property that for any term t there are only finitely many literals s such that s ≺ t. In case SCL does not detect a conflict with respect to a finite, exhaustive trail of ground literals, they constitute a model candidate for the clause set [START_REF] Bromberger | Exploring partial models with SCL[END_REF]. If SCL detects a conflict it learns a new first-order non-ground clause. It is derived by resolution and factoring with guidance from the trail. A natural choice for the ordering ≺ is the Knuth-Bendix (KBO) ordering [START_REF] Knuth | Simple word problems in universal algebras[END_REF]. For the ground case, a KBO relation can be efficiently computed [14]. All SCL calculi propagate literals from clauses with respect to the trail. For example, given a trail [P (a)] and a clause ¬P (x) ∨ R(x, y) the literal R(a, y) could be propagated. The SCL theory only enables ground literals on the trail, however, in practice it is not affordable to put all groundings of R(a, y) on the trail that are ≺ smaller than β. Therefore, we already considered trail literals with variables when we developed a two-watched literal scheme for SCL [START_REF] Bromberger | A two-watched literal scheme for first-order logic[END_REF]. Recall that this propagation situation is not exceptional as typically not all literals in a clause carry all occurring variables. The consequence of this extension is that for SCL we now need to decide solvability of conjunctions of inequations t i ≺ β where the t i may contain (shared) variables, i.e., we have to decide solvability of a particular form of KBO constraints if ≺ is the KBO.

For the SCL(EQ) calculus [13] the requirements on constraint solving get more sophisticated. Now the trail is a sequence of unit (in)equalities and propagation and conflicting clauses are decided with respect to the resulting congruence. For an extended congruence closure algorithm [START_REF] Downey | Variations on the common subexpression problem[END_REF][START_REF] Hurd | Congruence classes with logic variables[END_REF]15,16] we need now in addition to inequations t i ≺ β to consider inequalities t i = s i in order to separate congruence classes. In its simplest form, constraints consist of inequations t i ≺ β and inequalities t i = s i where β and the s i are ground, so called simple rightground constraints, Definition 4. In a more general setting, the s i carry variables and then a quantifier alternation on variables occurring in s i but not in t i needs to be considered. Such constraints are called alternating, Definition 27.

In this paper we investigate the complexity of all these variants with respect to a KBO <, Definitions 3, 4, 25, 27, but also a weaker non-strict ordering based on pure symbol counting, Definition 22. Except for constraints bound by a single ground term, Proposition 26, all problems are NP-hard, Propositions 5, 21, 24, 28.

Korovin and Voronkov developed a decision procedure [START_REF] Korovin | A decision procedure for the existential theory of term algebras with the Knuth-Bendix ordering[END_REF] for KBO constraints consisting of inequations s j < t j only and refined it to an NP algorithm [START_REF] Korovin | Knuth-Bendix constraint solving is NP-complete[END_REF]. According to Löchner [14], these results are "of more theoretical interest" because they are "too involved to be implemented with reasonable effort". In fact, to the best of our knowledge we present the first implemented algorithm for KBO constraint solving in this paper. Later, Korovin and Voronkov [START_REF] Korovin | Orienting rewrite rules with the Knuth-Bendix order[END_REF] showed that checking satisfiability of a KBO constraint consisting of a single inequation s < t can be done in polynomial time. For the special case of a right-ground constraint consisting of a single inequation s < t, what their algorithm essentially does is assigning the minimal constant to every variable.

To the best of our knowledge the problem of simple right-ground KBO constraints has never been studied before. We are also not aware of any implementation of a KBO constraint solving algorithm. The paper is now organized as follows: In Sect. 3 we prove the NP-completeness of this problem and present an algorithm to solve it. In Sect. 4 we study the complexity of variants of this problem including alternating constraints. We also consider a non-strict, partial ordering based on symbol counting and weaker than a KBO. The algorithm for right-ground constraints is extended to alternating constraints. In Sect. 5 we put the algorithm developed in Sects. 3 and 4 to practice and end the paper with a discussion of the obtained results, Sect. 6.

Preliminaries

In the following let Σ be a signature, i.e., a finite set of function symbols. Every function symbol f has an associated arity which we denote by arity(f ). Function symbols c with arity(c) = 0 are called constants. We denote the set of all terms by T (Σ, X ) where X is an infinite set of variables. Vars(t) denotes the set of variables occurring in the term t. A term t is called ground if it contains no variables, i.e., Vars(t) = ∅. The set of all ground terms is denoted by T (Σ). We assume that Σ contains at least one non-constant function and at least one constant, i.e., that T (Σ) is infinite. For otherwise, constraint solving becomes trivial. A substitution is a mapping σ : X → T (Σ, X ) such that σ(x) = x for only finitely many x ∈ X. The application tσ of a substitution σ to a term t ∈ T (Σ, X ) is defined in the usual way. We call a substitution grounding for some term t ∈ T (Σ, X ) if tσ is ground. A substitution σ is a matcher from s to t if sσ = t. We consider the following version of the Knuth-Bendix ordering (KBO) on ground terms: Definition 1 (KBO on Ground Terms [START_REF] Knuth | Simple word problems in universal algebras[END_REF]). Let be a strict total ordering (a precedence) on Σ, and w : Σ → N + a weight function. w is extended to terms recursively by w(f In particular, the precedence is strict and total, no unary function f with w(f ) = 0 is allowed and all weights are natural numbers. It can be shown that > KBO is a strict, total and well-founded ordering on ground terms. In the following, we simply write > for > KBO . Definition 2. A KBO constraint C is a finite set of atoms t#s where t, s ∈ T (Σ, X ) and # ∈ {<, >, =, ≤, ≥, =}. We say that C = {t 1 # 1 s 1 , . . . , t n # n s n } is satisfiable if there exists a substitution σ that is grounding for all t j , s j such that n j=1 t j σ # j s j σ. Such a grounding substitution σ is called a solution. Definition 3. A right-ground KBO constraint C is a KBO constraint where s 1 , . . . , s n ∈ T (Σ), i.e., only the t j may contain variables.

(t 1 , . . . , t n )) = w(f )+ n i=1 w(t i ).
Definition 4. A simple right-ground KBO constraint C is a right-ground KBO constraint where # ∈ {<, =}.
For simple right-ground KBO constraints, we prefer more explicit notation: We now assume t 1 , . . . , t n , l 1 , . . . , l m ∈ T (Σ, X ), s 1 , . . . , s n , r 1 , . . . , r m ∈ T (Σ) and call C satisfiable if there exists a substitution σ that is grounding for all t j , l j such that ⎛ ⎝ n j=1

t j σ < s j ⎞ ⎠ ∧ ⎛ ⎝ m j=1 l j σ = r j ⎞ ⎠ .
3 Simple, Right-Ground KBO Constraints

We start by investigating the complexity of simple, right-ground KBO constraint solving.

Proposition 5. Checking satisfiability for simple right-ground KBO constraints is NP-hard.

Proof. We reduce from MONOTONE 3SAT which is NP-complete by [START_REF] Garey | Computers and Intractability: A Guide to the Theory of NP-Completeness[END_REF]. Let N M be a set of clauses where N consists of the clauses with only positive literals and M consists of the clauses with only negative literals. We consider a signature with a constant a, a ternary function f and a unary function g. We use a KBO instance where all weights are 1 and f g a. For every propositional variable P occurring in N M , we introduce a variable x P . Then the equation x P = a stands for P is true and x P = a stands for P is false. Now every positive clause

(P ∨ Q ∨ R) ∈ N is encoded as an inequation f (x P , x Q , x R ) < f(g(a)
, g(a), g(a)). Obviously, this inequation can only be satisfied by a grounding that maps at least one of these variables to a, i.e., that sets at least one of P, Q, R to true.

Every negative clause (¬P ∨ ¬Q ∨ ¬R) ∈ M is encoded as an inequality f (x P , x Q , x R ) = f (a, a, a). Obviously, this can only be satisfied if not all of these variables are mapped to a, i.e., if at least one of P, Q, R is false. Now the clause set has a solution iff there is a solution to the constructed simple right-ground KBO constraint. Assume N M is satisfiable by a valuation β. Then for every propositional variable P map x P to a if β(P ) = 1 and to g(a) otherwise. As explained above, this grounding will satisfy the constraint. Now let σ be a solution to the constraint. Then the valuation β where β(P ) = 1 if σ(x P ) = a and β(P ) = 0 otherwise satisfies N M .

We have added |M | inequalities and |N | inequations which can be constructed in polynomial time, so the reduction works in polynomial time. Proposition 6. Checking satisfiability for simple right-ground KBO constraints is in NP.

Proof. Let C = {t 1 < s 1 , . . . , t n < s n , l 1 = r 1 , . . . , l m = r m } be a constraint. If for some inequality l j = r j , there is no matcher from l j to r j , we can ignore this inequality since it is true for every grounding. If for some inequality l j = r j , it actually holds that l j = r j , then this inequality is impossible to satisfy, so we are done. After sorting out these two cases, as r j is ground, every inequality l j = r j has a unique matcher τ j which has linear size with respect to r j . In the following, we say that the term τ j (x) is restricted by the inequality l j = r j . The inequality l j = r j then signifies

x∈Vars(lj ) σ(x) = τ j (x).
For the inequations t j < s j , it is obviously optimal to assign the smallest possible term to every variable. Larger terms only have to be considered due to the inequalities l j = r j . If there is a grounding σ that satisfies t j < s j , then any grounding σ with σ (x) ≤ σ(x) for all variables x satisfies t j < s j . Hence, if there exists a solution, then there also exists a solution that only uses the m + 1 smallest terms for every variable. This is because every inequality l j = r j only restricts at most one term for every variable, so for every variable the m + 1 smallest terms contain the smallest term that is not restricted for that variable.

As we only have to consider the m + 1 smallest terms for every variable, the size of the groundings we have to consider is polynomially bounded by the input size. Let f be the function with the maximal arity and let p = arity(f ). Let a be the smallest constant. We claim that every of the m + 1 smallest terms has at most mp + 1 symbols. Proof by contradiction: Assume t 0 is one of the m + 1 smallest terms with #t 0 > mp + 1. Perform the following m times: Obtain t i+1 by replacing any subterm g(s 1 , . . . , s n ), where the s i are constants, by a. The number of symbols decreases by at most p, so #t i > (mi)p + 1. As none of the t i is a constant, such a subterm always exists. After m steps, we obtain terms

t 0 > t 1 > • • • > t m with #t m > (m -m)p + 1 = 1, i.e., t m is not a constant, so t m > a.
This contradicts the fact that t 0 was one of the m + 1 smallest terms since at least m + 1 terms are smaller than t 0 . Thus, we can guess a grounding and check in polynomial time whether it is a solution.

Next we propose an algorithm for testing satisfiability of simple right-ground KBO constraints. Of course, by Proposition 6, there already exists an algorithm, but we expect that the following algorithm performs better in practice. Let C be a simple right-ground KBO constraint with n inequations t j < s j and m inequalities l j = r j .

Assume that Vars ({t

j | 1 ≤ j ≤ n} ∪ {l j | 1 ≤ j ≤ m}) = {x 1 , . . . , x k }.
As explained in the proof of Proposition 6, we only have to consider the m + 1 smallest terms for the grounding, so to begin, we generate an ordered list S of the m + 1 smallest terms. This way, a grounding substitution σ corresponds to a vector v ∈ N k where v

i < m + 1 is the index of the term σ(x i ) in S, i.e., S[v i ] = σ(x i ). Let σ( v) with σ( v)(x i ) := S[v i ]
denote the grounding corresponding to the vector v. Later on, we give a dynamic programming algorithm to compute the k smallest terms for some number k. Actually, we do not directly generate the m + 1 smallest terms, but start with a constant number of terms and generate more terms as needed.

The algorithm is given by three inference rules that are represented by an abstract rewrite system. They operate on a state which is either ⊥ or a four-tuple (T ; v; F ; C) where T is a sequence of variables, the trace; v ∈ N k is a grounding substitution in vector notation, the current grounding; F is a set of forbidden groundings; and C is a simple right-ground KBO constraint. The initial state for a constraint C is (ε; (0, . . . , 0); ∅; C), i.e., the trace is empty, every variable is mapped to the smallest constant and there are no forbidden groundings.

We use the following partial ordering ≤ F on groundings: v ≤ F u iff for all i ∈ {1, . . . , k} we have v i ≤ u i . By inc( v, i) we denote the grounding v with v i = v i + 1 and v l = v l for all l ∈ {1, . . . , k} with l = i, i.e., the grounding where we increase the term for the variable x i by one. Analogously, we define dec( v, i), where we instead decrease the term for the variable x i by one, i.e., v i = v i -1. The two operations inc and dec are only used when they are well-defined, i.e., they yield a grounding v ∈ N k where v i < m + 1. The operation inc is only used when an inequality l j = r j is not satisfied, and this can happen at most m times without intermediate Backtrack steps. The operation dec( v, i) is only used for Backtrack, and by Lemma 15, in this case v i > 0.

The role of F is that we want to keep the algorithm from considering wrong groundings again. For all u ∈ F , we do not visit states with grounding v if v ≥ F u. When we Backtrack, we insert the current grounding into F . The trace T records the last updated variables so Backtrack is able to undo the last Increase operation. As will be proven in Theorem 18, the algorithm terminates in ⊥ iff there exists no solution, and if there exists a solution, then it terminates in a state where the current grounding v is a solution.

Increase (T ; v; F ; C) ⇒ KCS (T x i ; v ; F ; C) provided v = inc( v, i), l j σ( v) = r j for some l j = r j ∈ C, l j σ( v ) = r j and there is no u ∈ F with v ≥ F u Backtrack (T x i ; v; F ; C) ⇒ KCS (T ; v ; F ∪ { v}; C) provided v = dec( v, i
) and either 1. l j σ( v) = r j for some l j = r j ∈ C, but for all l ∈ {1, . . . , k}, we have that

l j σ(inc( v, l)) = r j implies that there is a u ∈ F with inc( v, l) ≥ F u, or 2. t j σ( v) ≥ s j for some t j < s j ∈ C Fail (ε; v; F ; C) ⇒ KCS ⊥
provided either 1. l j σ( v) = r j for some l j = r j ∈ C, but for all l ∈ {1, . . . , k}, we have that

l j σ(inc( v, l)) = r j implies that there is a u ∈ F with inc( v, l) ≥ F u, or 2. t j σ( v) ≥ s j for some t j < s j ∈ C
Informally, Increase is applicable if some inequality l j = r j is not fulfilled and we can fix this with the new grounding inc( v, i) which is not forbidden by F . Backtrack undoes an operation and is applicable if either some inequality l j = r j is not fulfilled, but Increase is not applicable, or if some inequation t j < s j is not fulfilled. Fail is applicable if Backtrack would be applicable on an empty trace, i.e., there is no operation to undo.

Obviously, there is no state on which we can apply both Backtrack and Fail.

Definition 7. A reasonable strategy is a strategy that prefers Backtrack and Fail over Increase.

Example 8. Consider a signature with constants a, b, c and a binary function f . We set w(a) = 1;

w(b) = w(c) = 2; w(f ) = 3 and a ≺ b ≺ c ≺ f . We consider the constraint C = {x 1 = a, f (x 1 , x 2 ) < f(a, c)}.
The m +1 smallest terms, where m = 1, are a, b. This is the unique execution of the algorithm. In order to increase readability, for v, we write the terms instead of the indices.

(ε; (a, a); ∅; C) ⇒ Increase KCS (x 1 ; (b, a); ∅; C) ⇒ Backtrack KCS (ε; (a, a); {(b, a)}; C) ⇒ Fail KCS ⊥
The algorithm terminates in ⊥, so there is no solution.

Example 9. Consider a signature with constants a, b, a binary function g and a ternary function f . Let

w(a) = 1, w(b) = w(f ) = w(g) = 2 and a ≺ b ≺ g ≺ f . The constraint is C = {x 1 < b, g(x 2 , a) < g(b, b), f(x 1 , x 2 , x 3 ) = f (a, a, a), g(x 1 , x 2 ) = g(a, b)}.
The m + 1 smallest terms, where m = 2, are a, b, g(a, a). The algorithm has found a solution, so no rule is applicable and it terminates. Note that after the third and fifth operation, we cannot increase x 1 because (b, b, a) ≥ F (b, a, a) ∈ F .

Next we prove the correctness of the algorithm.

Lemma 10. If (ε; (0, . . . , 0); ∅; C) ⇒ l KCS (T ; v ; F ; C), then there is no u ∈ F with v ≥ F u.
Proof. We prove this by induction on l. For l = 0, this holds since F = ∅. For l > 0, the last applied rule must have been either Increase or Backtrack. If the last applied rule was Increase, then there cannot be such a u because Increase does not modify F and because this is part of the condition of the Increase rule. Now assume the last applied rule was Backtrack, so the previous state was (T x i ; v; F ; C) with v = dec( v, i) and F = F ∪ { v}. If there was some u in F such that v ≥ F u, then, since v < F v, we have v > F u. Hence, by the induction hypothesis, u / ∈ F , so as

F = F ∪ { v}, it must hold that u = v, contradiction to v > F u. Lemma 11. If (ε; (0, . . . , 0); ∅; C) ⇒ i KCS (T ; u; F ; C) ⇒ l KCS (T ; u ; F ; C) for l > 0, then u = u or F = F .
Proof. If all l rule applications are applications of the Increase rule, then clearly u < F u , so in particular, u = u . There is no rule that removes elements from F , so F ⊆ F . If there is at least one application of the Backtrack rule among the l rule applications, the current assignment v is added to F , and by Lemma 10, v / ∈ F , so F is modified and F = F .

Proposition 12. ⇒ KCS is well-founded, i.e., the algorithm always terminates.

Proof. By Lemma 11, we can reach every combination of v and F at most once.

For v, there are (m + 1) k possibilities. We only add occurring groundings to F , so the number of possibilities for F is upper bounded by the number of subsets of all possible groundings which is 2 (m+1) k . Thus, the number of reached states is finite (it is at most (m + 1) k 2 (m+1) k ), so the algorithm terminates.

Of course, the upper bounds in the proof of Proposition 12 are far too high and the algorithm will run much faster in practice.

Lemma 13. If (ε; (0, . . . , 0); ∅; C) ⇒ l KCS (T ; v; F ; C) and u ∈ F , then for all u ≥ F u it holds that u cannot be a solution.

Proof. The proof is by induction on l. For l = 0, we have F = ∅, so this holds. For l > 0, if the last applied rule was Increase, the statement follows by the induction hypothesis since F is not modified. Now assume that the last applied rule was Backtrack. Let (T ; v ; F ; C) be the previous state. We only have to show that all u ≥ F v cannot be solutions, for all other elements of F = F ∪ { v }, this follows by the induction hypothesis. First assume that Backtrack is applicable because of condition [START_REF] Barrett | The Satisfiability Modulo Theories Library (SMT-LIB)[END_REF]. Then v cannot be a solution since l j σ( v ) = r j . For u > F v , if l j σ( u ) = r j , then u clearly cannot be a solution. Otherwise, there is a variable x i such that u ≥ F inc( v , i) and l j σ(inc( v , i)) = r j . However, it is part of condition (1) that then, there is an element u ∈ F with u ≤ F inc( v , i) ≤ F u , so by the induction hypothesis, u cannot be a solution. If Backtrack is applicable because of condition (2), then t i σ( v ) ≥ s i for some i ∈ {1, . . . , n}. Clearly, if u ≥ F v , then also t i σ( u ) ≥ s i , so u cannot be a solution.

Corollary 14. If (ε; (0, . . . , 0); ∅; C) ⇒ l KCS (T ; v; F ; C) and condition [START_REF] Barrett | The Satisfiability Modulo Theories Library (SMT-LIB)[END_REF] or condition [START_REF] Bromberger | Deciding the Bernays-Schoenfinkel fragment over bounded difference constraints by simple clause learning over theories[END_REF] of Fail is fulfilled for (T ; v; F ), then for all u ≥ F v, u cannot be a solution.

Proof. The conditions for Fail are the same as the conditions for Backtrack, so this follows by the proof of Lemma 13.

Lemma 15. If (ε; (0, . . . , 0); ∅; C) ⇒ l KCS (T ; v; F ; C), then for all i ∈ {1, . . . , k} the number of occurrences of x i on the trace T equals v i .

Proof. In the following, we denote the number of occurrences of x i in T by C(T, x i ). The proof is by induction on l. If l = 0, the statement trivially holds. For l > 0 let (T ; v ; F ; C) be the previous state. If the last applied rule was Increase, then

T = T x i and v = inc( v , i), so C(T, x i ) = C(T , x i ) + 1 IH = v i + 1 = v i .
For j = i, C(T, x j ) = C(T , x j ) and v j = v j , so the statement follows by the induction hypothesis. If the last applied rule was Backtrack, then

T = T x i and v = dec( v , i), so C(T, x i ) = C(T , x i ) -1 IH = v i -1 = v i .
Again, for j = i, C(T, x j ) = C(T , x j ) and v j = v j , so the statement follows by the induction hypothesis. Proof. Since Fail is applicable on (T ; v; F ; C), T = ε, so by Lemma 15, v = (0, . . . , 0). Hence, by Corollary 14, for all u ≥ F (0, . . . , 0), u cannot be a solution, so there exists no solution.

Lemma 17. If (ε; (0, . . . , 0); ∅; C) ⇒ l KCS (T ; v; F ; C) and no rule is applicable on (T ; v; F ; C) then v is a solution.

Proof. Assume that for some j ∈ {1, . . . , n}, we had t j σ( v) ≥ s j . Then, either Backtrack or Fail would be applicable. Now assume that for some j ∈ {1, . . . , m}, we had l j σ( v) = r j . Then, either Increase or Backtrack would be applicable.

Theorem 18. The algorithm is correct: If there exists a solution, then starting from (ε; (0, . . . , 0); ∅; C), the algorithm terminates in a state (T ; v; F ; C) where v is a solution. If there is no solution, the algorithm terminates in ⊥.

Proof. Follows by Proposition 12, Lemma 16 and Lemma 17.

We have implemented the above algorithm in the context of the SPASS reasoning workbench. The efficiency of the algorithm depends on the respective variables we choose for Increase. If there exists a solution, then there exists an execution using only the rule Increase. The following criteria might be useful to select the best variable for Increase:

-We prefer variables that do not occur in "critical" inequations, or in a minimal number of inequations. A "critical" inequation is one where the weight difference is 0 or close to 0. -We prefer variables x i for which the next term is not restricted by any inequality l j = r j . -We prefer variables x i for which the next term does not have a larger weight, or for which the increase in weight is minimal. -We prefer variables that fix multiple inequalities l j = r j instead of just one.

It is possible to calculate and maintain some score for every variable here and decide based on this score. The exact selection criteria still need to be further explored.

A remaining problem from the presentation of the algorithm is how to compute the k smallest terms. If the occurring weights are rather small, the following dynamic programming algorithm might be useful in practice. The idea is to compute all terms of a specific weight for increasing weights until we generated at least k terms. Unfortunately, there may be exponentially many terms of a specific weight where the exponent is the maximal arity of a function and the base is the number of terms of smaller weights. However, k is bounded above by the number of inequalities m, the number of terms with smaller weights is bounded above by k and the maximal arity is probably small, so it is to be expected that this is not a big problem.

As it is probably hard to find the next possible weight, we simply always increase the weight by 1 starting by the weight of the smallest constant. Our DP array is two-dimensional, one dimension having the weight and the other dimension having the size of the tuple from 1 to max_arity. Actually, it is fourdimensional since every entry is a list of tuples of terms and every tuple is a list of its entries. A tuple of size 1 is just a term of the specific weight. The tuples of larger size are needed for the DP transitions where they serve as argument tuples for the functions. We maintain an array smallest_terms that will in the end contain the at least k smallest terms.

We iterate over the weights starting at the weight of the minimal constant. Let curweight denote the current weight. The idea is to compute all terms of weight curweight, sort them, add them to smallest_terms, and proceed with weight curweight + 1 if |smallest_terms| is still smaller than k. To do so, if curweight is not the smallest weight, we first compute the tuples of size 2 to max_arity for the previous weight. This is done via DP: For tuple size i we iterate over the terms s ∈ smallest_terms. Then we iterate over the tuples t of size i -1 and weight curweight -1w(s) using the DP array and add (s, t) to the current DP entry. Afterwards, we calculate all terms of weight curweight by iterating over all symbols f and all tuples t of size arity(f ) and weight curweightw(f ) using the DP array. Then, the term f (t) has weight curweight.

We finish this section by a discussion of potential heuristics, sufficient conditions for a simple right-ground KBO constraint to have a solution. As explained before, every inequality l j = r j rules out any assignment that satisfies τ j , the matcher from l j to r j . Now assume we have m inequalities and know that there are more than m solutions for the inequation t < s, then one might think that there is a grounding that solves all inequalities l j = r j and the inequation t < s. However, this is not true.

C = {x = a, f (x, y) < f(a, c)}.
The inequation has two solutions, namely {x → a, y → a} and {x → a, y → b}. However, it has no solution where x is not mapped to a, so for the overall problem, there is no solution.

So the above sufficient condition needs to refined in order to be correct. However, calculating the number of solutions is again NP-hard.

Proposition 20. Calculating the number of solutions σ for some right-ground inequation t < s is NP-hard.

Proof. We reduce from the Unbounded Subset Sum Problem (USSP) which is NP-complete by [START_REF] Garey | Computers and Intractability: A Guide to the Theory of NP-Completeness[END_REF]. Let s 1 , . . . , s n , T ∈ N + . We have to find out whether there are x 1 , . . . , x n ∈ N such that n i=1 x i s i = T , i.e., whether there is a multiset of values from {s 1 , . . . , s n } that sums up to T . Assume we had an oracle that could compute the number of solutions for any inequation l < r where r is ground. We will use this oracle twice.

For both uses, we use a signature with constants c and d and unary functions f 1 , . . . , f n . We have

w(c) = 1, w(f i ) = s i for i ∈ {1, . . . , n} and d ≺ c ≺ f 1 ≺ • • • ≺ f n .
For the first case, set w(d) = T +2. Using the oracle with the inequation x < d, we get the number of terms smaller than d. Since d is the smallest term of weight T + 2, this is exactly the number of terms with weight ≤ T + 1. For the second case, set w(d) = T + 1. Again, using the oracle with the inequation x < d, we get the number of terms smaller than d. This time, this is the number of terms with weight ≤ T . If we now subtract those values, we get the number of terms with weight exactly T + 1.

Now the USSP has a solution iff the number of terms with weight exactly T + 1 is not 0. Every term t of weight T + 1 must have the constant c as subterm since the weight of d is too large. The rest of t must consist of the unary functions. Hence, the weights of the unary functions used sum up to T + 1 -1 = T . Since the weights of the unary functions correspond to the numbers from the USSP, this yields a solution for the USSP. Conversely, given a solution to the USSP, we can construct a term of weight T + 1 analogously.

The problem with the aforementioned insufficient condition is that an inequality l j = r j does not necessarily rule out only one grounding, but possibly infinitely many groundings. This happens if there are variables that are not restricted by the matcher τ j of l j and r j . However, the criterion can be refined to a correct sufficient condition. If we restrict ourselves to the m + 1 smallest terms again, we would again at least have a finite number of groundings that l j = r j rules out. If we now sum up these numbers over all inequalities, we have an upper bound on the total number of ruled out groundings. For the inequation t < s, the same problem with variables that do not occur arises (there may be infinitely many solutions), so here, we restrict ourselves to the m + 1 smallest terms again. If now, the number of solutions for t < s is larger than the upper bound on the total number of ruled out groundings, we can actually be sure that there is a solution. However, this correct sufficient condition is hard to compute and therefore seems to be not very useful in practice.

Further Constraint Variants and Ordering Relaxation

In this section we study further variants of constraint problems and eventually extend the algorithm of Sect. 3 to alternating KBO constraints.

Proposition 21. Checking satisfiability for right-ground KBO constraints restricted to strict inequations is NP-hard.

Proof. The proof strategy is the same as the one used in the proof of Proposition 5. The encoding for positive clauses stays the same as < is still allowed. For negative clauses ¬P ∨ ¬Q ∨ ¬R we encode them as f (x P , x Q , x R ) > f(a, a, a). This inequation can only be satisfied by a grounding that does not map all of these variables to a, and is trivially satisfied by any such grounding.

In particular, we have seen that only having constraints of the form t i < s i and t i > s i suffices to make the problem NP-hard. Next we turn to a weaker term ordering ≤ sym solely based on symbol counting. Even for this ordering constraint solving remains NP-hard.

Definition 22. For ground terms t, s ∈ T (Σ), we define t ≤ sym s : ⇐⇒ |sym(t)| ≤ |sym(s)|, i.e., t does not contain more symbols than s. Definition 23. A right-ground symbol constraint C is a finite set of atoms t#s with t ∈ T (Σ, X ), s ∈ T (Σ) and # ∈ {≤ sym , =}. Satisfiability is defined analogously to the satisfiability of KBO constraints.

Proposition 24. Checking satisfiability for right-ground symbol constraints is NP-hard.

Proof. The proof strategy is the same as the one used in the proof of Proposition 5. We encode positive clauses

P ∨ Q ∨ R as f (x P , x Q , x R ) ≤ f (g(a)
, g(a), a). The only way to satisfy this inequation is to map at least one of these variables to a.

Negative clauses ¬P ∨ ¬Q ∨ ¬R are encoded as f (x P , x Q , x R ) = f (a, a, a).
In particular, the NP-hardness of these problems is not caused by the complicated structure of the KBO since the problem is already NP-hard for a comparison as simple as counting the number of symbols.

Our next variants are motivated by the definition of congruence classes with respect to terms with variables. For the first variant, all instances of the defining term t have to be smaller than a single ground term β and different from ground terms s 1 , . . . , s n . Definition 25. A simple, single ground KBO constraint C consists of terms t ∈ T (Σ, X ) and s 1 , . . . , s n , β ∈ T (Σ). We say that C is satisfiable if there exists a substitution σ that is grounding for t such that

⎛ ⎝ n j=1 tσ = s j ⎞ ⎠ ∧ tσ < β.
Proposition 26. Assuming that we are given the n+1 smallest terms, checking satisfiability of simple, single right-ground KBO constraints is in P.

Proof. Actually, for this problem, if a reasonable strategy (Definition 7) is used, the algorithm from Sect. 3 runs in polynomial time. The key difference to the other problems is that here, every variable occurs in every inequality, so every inequality rules out at most one grounding. We first show that we can only reach polynomially many states. First, consider states (T ; v; F ; C) where tσ( v) < β. If v does not violate any inequality t = s i , then the algorithm terminates, so there is at most one such state. For every inequality t = s i , there is at most one grounding v that violates it. We claim that we reach at most k + 1 states with current grounding v where k is the number of variables. v is reached at most once using Increase because otherwise, there must be an intermediate Backtrack, so v would be inserted into F . If we reach v using Backtrack for some variable x i , then inc( v, i) was inserted into F , so for every variable x i , we can reach v using Backtrack at most once. Hence, v is reached at most k + 1 times. Now, consider states (T ; v; F ; C) where tσ( v) ≥ β. Since a reasonable strategy is used, we must have reached this state from a state that does not violate t < β, so by the argumentation before, at most k such states can be reached for every inequality, so at most n • k in total where n is the number of inequalities.

Hence, in total, there are at most (k + 1) • n + n • k + 1 states. The state transitions can be done in polynomial time because we only need to iterate over all inequalities and inequations and over all entries in F . Since there are only polynomially many states and every rule application inserts at most one element into F , F has polynomially many entries.

Definition 27 (Alternating KBO Constraint). An alternating KBO constraint C consists of terms t, s 1 , . . . , s n ∈ T (Σ, X ) and β ∈ T (Σ). We say that C is satisfiable if there exists a substitution σ that is grounding for t such that for all substitutions τ that are grounding for all s j we have

⎛ ⎝ n j=1 tσ = s j τ ⎞ ⎠ ∧ tσ < β.
Proposition 28. Checking satisfiability for alternating KBO constraints is NPhard.

Proof. We reduce from SAT. Let N be a set of clauses and X 1 , . . . , X k be the variables occurring in N . We use a signature with a k-ary function f , two constants a and b, variables x 1 , . . . , x k and y 1 , . . . , y k . Set t = f (x 1 , . . . , x k ). Now for every clause C j ∈ N we introduce an inequality f (x 1 , . . . , x k ) = f (s j,1 , . . . , s j,k ) where we set s j,i = b if X i occurs positively in C j , s j,i = a if X i occurs negatively in C j and s j,i = y i if X i does not occur in C j . The idea is that x i = a stands for X i is set to and x i = b stands for X i is set to ⊥. ∀τ.x i σ = y i τ is obviously impossible to satisfy, so the inequality must be made true by setting some positive variable to a or some negative variable to b.

To ensure that the x i are only mapped to a or b, we do the following: We first introduce a new constant c and set β = c. Then we set w(f ) = w(a) = w(b) = 1, w(c) = k + 2 and c ≺ a ≺ b ≺ f . If all variables x i are mapped to a or b, we have w(tσ) = k + 1, i.e., tσ < β. Any grounding where some x i is not mapped to a or b results in tσ ≥ β. Now there is a solution σ iff there is a satisfying valuation for N .

If a reasonable strategy is used, satisfiability of alternating KBO constraints can be checked using the algorithm from Sect. 3. Any solution σ must be such that tσ < β, so we only have to consider instances of the s j τ with s j τ < β. What we can now do is to calculate for all s j all groundings τ with s j τ < β and add the inequality t = s j τ to the constraint. There are only finitely many such groundings because we did not allow unary functions f with w(f ) = 0. This way, we obtain a simple right-ground KBO constraint, so we can apply the algorithm. A more efficient possibility to do this is to add the groundings of the s j implicitly, i.e., to change the condition of Increase (and the first case of Backtrack and Fail) to whether there exists a matcher τ such that l j σ( v) = r j τ . Also, the condition for the next grounding for Increase changes: It is not that we fix the inequality anymore, but that we change a variable that occurs on the left side of the inequality.

Example 29. Consider the signature Σ = {f (2) , g (1) , a (0) }, where the superscript numbers denote the function arities, together with the following alternating KBO constraint C:

t = f (x 1 , x 2 ) s 1 = f (g(y 1 ), y 2 ) β = f (f (a, a), a) s 2 = f (a, a)
We set w(a) = w(g) = w(f ) = 1 and a ≺ g ≺ f . The few smallest terms are a, g(a), g(g(a)), f(a, a).

Note that for alternating KBO constraints, it does not suffice anymore to consider the n + 1 smallest terms only since an inequality may rule out more than one term for a variable. However, as mentioned in Sect. 3, we calculate the smallest terms as needed, so this is not a problem. For shorter notation, for F , we omit groundings u if there is a grounding v ∈ F with v < F u. A possible run of the algorithm looks as follows: We implemented the algorithm of Sect. 3 and its extension to constraints with right hand side variables, Definition 27, and tested it in the context of an extended congruence closure (CC) algorithm with variables [START_REF] Downey | Variations on the common subexpression problem[END_REF][START_REF] Hurd | Congruence classes with logic variables[END_REF]15,16]. We implemented a rather naive variant of [START_REF] Hurd | Congruence classes with logic variables[END_REF] with the only goal to generate KBO constraints in order to test our new algorithm on KBO constraints. In contrast to [START_REF] Hurd | Congruence classes with logic variables[END_REF] our algorithm considers a finite signature, as usual for first-order logic problems. All experiments were carried out on a Debian Linux server equipped with AMD EPYC 7702 64-Core CPUs running at 3.35GHz and an overall memory of 2TB. The result of all runs as well as all input files and binaries can be found at https://nextcloud.mpi-klsb.mpg.de/index.php/s/BAwd99cxFpSJmSp. As a first test case we considered all eligible UEQ problems from CASC-J11 [17]. We consider equations and all inequalities for the congruence closure algorithm. The equations generate the congruence and for the inequalities we compute the congruence classes for the respective right and left side term of the inequality. For each example, the KBO function weight was always set to one and the precedence is generated with respect to the occurrence of symbols in the input file in ascending order. For β we chose a fixed nesting depth of 4 and build for each input file a nested term of exactly this depth using function symbols in the order of occurrence in the input, starting with a non-constant function symbol. Out of all eligible problems our CC algorithm terminated on 186 problems within a time limit of 30 min. Please note that although our CC implementation is rather naive, in contrast to the classical ground CC algorithm it does not need a complete grounding; for the examples where our naive algorithm runs out of time a complete grounding is not affordable. The below table shows some typical runs on the UEQ domain. All timings are presented in hundredths of a second and if they take less than one hundredth of a second we write zero. The below table shows the problem name, the number of ground terms smaller than β indicating the solution space for the constraint, the summed up time of all calls to the KBO constraint solver during the CC run, the number of calls to the KBO constraint solver, and the results of these calls. The three selected examples are typical: most of the problems are satisfiable and the constraint solving algorithm needs almost no time. Note that for the first example all 8014 calls to the constraint solver needed in sum 3 hundreds of a second. The LAT143-1 is the example showing the worst constraint solving performance, i.e., still less than a hundreth of a second per call. For the SMT-LIB examples of the UF domain [START_REF] Barrett | The Satisfiability Modulo Theories Library (SMT-LIB)[END_REF], we expanded let operators, removed the typing, coded predicates as equations, did a CNF transformation and then took the first literal of each clause as input for the CC algorithm. Nesting depth was set to 2, the rest done as for the UEQ examples. Removing types means that the number of smaller terms increases, i.e., the problems get potentially more difficult for the constraint solver, in particular for unsatisfiable constraints. The below Here uf.555113 is the worst example on constraint solving time with 1.34 s for 5120 calls. Although alternating KBO constraint solving is NP-hard, in practice there are typically only a few inequalities meaning that out of the overall number of terms smaller β, only a few need to be considered.

Discussion

We have studied a number of specific KBO constraint solving problems motivated by the SCL calculus and established their complexity. Except for simple, single right-ground KBO constraints all studied problems are proven NP-hard. We propose an algorithm that eventually runs for alternating KBO constraints which include a quantifier alternation. The algorithm shows nice performance on benchmark problems. Our next step is to turn our naive CC implementation with variables into a robust algorithm. The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

  (ε; (a, a, a); ∅; C) ⇒ Increase KCS (x 1 ; (b, a, a); ∅; C) ⇒ Backtrack KCS (ε; (a, a, a); {(b, a, a)}; C) ⇒ Increase KCS (x 2 ; (a, b, a); {(b, a, a)}; C) ⇒ Increase KCS (x 2 x 2 ; (a, g(a, a), a); {(b, a, a)}; C) ⇒ Backtrack KCS (x 2 ; (a, b, a); {(b, a, a), (a, g(a, a), a)}; C) ⇒ Backtrack KCS (ε; (a, a, a); {(b, a, a), (a, g(a, a), a), (a, b, a)}; C) ⇒ Increase KCS (x 3 ; (a, a, b); {(b, a, a), (a, g(a, a), a), (a, b, a)}; C)

Lemma 16 .

 16 If (ε; (0, . . . , 0); ∅; C) ⇒ l KCS (T ; v; F ; C) ⇒ Fail KCS ⊥, then there exists no solution.

Example 19 .

 19 Consider a signature with constants a, b and c and a binary function f . The weights are w(a) = 1; w(b) = w(c) = 2; w(f ) = 3 and we use a ≺ b ≺ c ≺ f as a precedence. Now consider the constraint

1 ⇒ 1 ⇒

 11 (ε; (a, a); ∅; C) ⇒ Increase KCS (x 1 ; (g(a), a); ∅; C) s 2 , τ = {} ⇒ Increase KCS (x 1 x 1 ; (g(g(a)), a); ∅; C) s 1 , τ = {y 1 → a, y 2 → a} ⇒ Increase KCS (x 1 x 1 x 1 ; (f (a, a), a); ∅; C) s 1 , τ = {y 1 → g(a), y 2 → a} ⇒ Backtrack KCS (x 1 x 1 ; (g(g(a)), a); {(f (a, a), a)}; C) β ⇒ Backtrack KCS (x 1 ; (g(a), a); {(g(g(a)), a)}; C) s Backtrack KCS (ε; (a, a); {(g(a), a)}; C) s Increase KCS (x 2 ; (a, g(a)); {(g(a), a)}; C) s 2 , τ = {} 5 Experiments

  The Knuth-Bendix ordering > KBO induced by and w is defined by s >

KBO t iff 1. w(s) > w(t), or 2. w(s) = w(t), and (a) s = f (s 1 , . . . , s m ), t = g(t 1 , . . . , t n ) and f g, or (b) s = f (s 1 , . . . , s m ), t = f (t 1 , . . . , t m ) and (s 1 , . . . , s m ) > lex KBO (t 1 , . . . , t m ).

  table again shows some typical results. 1112 examples could be performed by the CC algorithm inside 30 min. The UF domain contains larger examples compared to the UEQ domain, but the characteristics remain. Constraint solving itself takes almost no time. Again all timings are presented in hundredths of a second.

	Problem < β	Time KBO Constraint Solver #calls # true # false
	00336	2120806	0	131	131	0
	uf.926761 138397692 0	3882 1988	1894
	uf.555113 254939	134	5120 3306	1814

  13. Leidinger, H., Weidenbach, C.: SCL(EQ): SCL for first-order logic with equality. In: Blanchette, J., Kovács, L., Pattinson, D. (eds.) IJCAR 2022. LNCS, vol. 13385, pp. 228-247. Springer, Cham (2022). https://doi.org/10.1007/978-3-031-10769-6_14 14. Löchner, B.: Advances in Equational Theorem Proving-Architecture, Algorithms, and Redundancy Avoidance. Dissertation, Fachbereich Informatik, TU Kaiserslautern (2005) 15. Nelson, G., Oppen, D.C.: Fast decision procedures based on congruence closure. J. ACM 27(2), 356-364 (1980) 16. Shostak, R.E.: Deciding combinations of theories. J. ACM 31(1), 1-12 (1984) 17. Sutcliffe, G.: The CADE ATP system competition -CASC. AI Mag. 37(2), 99-101 (2016) Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

Acknowledgments. We thank our reviewers for their constructive comments that helped us improve the paper.


Constantin Ardilouze

Aaron A Boone

Impact of initializing the soil with a thermally and hydrologically balanced state on subseasonal predictability

Keywords: Soil enthalpy, Tibetan plateau, Subseasonal predictability, initialization, land-atmosphere coupling, S2S, LS4P

Accurate soil moisture initial conditions in dynamical subseasonal forecast systems are known to improve the temperature forecast skill regionally, through more realistic water and energy fluxes at the land-atmosphere interface. Recently, results from a multi-model coordinated experiment have provided evidence of the primal contribution of the initial surface and subsurface soil temperature over the Tibetan Plateau for capturing a hemispheric scale atmopsheric teleconnection leading to improved subseasonal forecasts. Yet, both the soil temperature and water content are key components of the soil enthalpy and we hypothesize that properly initializing one of them without modifying the other in a consistent manner can alter the soil thermal equilibrium, thereby potentially reducing the benefit of land initial conditions on subsequent atmospheric forecasts. This study builds on the protocol of the above-mentioned multi-model experiment, by testing three different land initialization strategies in an Earth system model. Results of this pilot study suggest that a better mass and energy balance in land initial conditions of the Tibetan Plateau triggers a wave train which propagates through the northern hemisphere mid-latitudes, resulting in an improved large scale circulation and temperature anomalies over multiple regions of the globe. While this study is based on a single case, it strongly advocates for enhanced attention towards preserving the soil energy equilibrium at initialization to make the most of land as a driver of atmospheric extended-range predictability.

Introduction

At the frontier between meteorological and climatic time scales, subseasonal forecasting, often referred to as S2S (Subseasonal to seasonal), remains a challenge for the scientific community. While ocean interannual variability is the dominant source of CNRM, Université de Toulouse, Météo-France, CNRS, Toulouse, France seasonal predictability at longer lead times, a substantial number of studies have emphasized the key role of continental surfaces in predictability at time scales between two weeks and two months (e.g. [START_REF] Dirmeyer | Land-atmosphere interactions and the water cycle. Seamless prediction of the Earth system: from minutes to months[END_REF]. Soil moisture, in particular, has received special attention since the seminal work of [START_REF] Koster | Regions of strong coupling between soil moisture and precipitation[END_REF] which identified certain regions of the globe as hotspots of soil moisture-atmosphere coupling. Predictability experiments using numerical models that followed this work such as GLACE2 multi-model experiment [START_REF] Koster | Contribution of land surface initialization to subseasonal forecast skill: First results from a multi-model experiment[END_REF] and the pilot study from [START_REF] Douville | Relative contribution of soil moisture and snow mass to seasonal climate predictability: a pilot study[END_REF] confirmed that a realistic initialization of soil moisture in dynamic forecasting systems could improve 2-meter temperature forecasts, particularly in socalled "water-limited" semi-arid regions. The snow-atmosphere coupling, which plays a role in the surface radiation and water balance [START_REF] Xu | Snow-atmosphere coupling strength. Part II: Albedo effect versus hydrological effect[END_REF], explains the importance of snow cover initialization for subseasonal forecasting [START_REF] Orsolini | Impact of snow initialization on sub-seasonal forecasts[END_REF]. More recently, a teleconnection has been identified between the surface and subsurface ground temperature of high altitude mountains and precipitation downstream of the flow [START_REF] Xue | Spring land temperature anomalies in northwestern us and the summer drought over southern plains and adjacent areas[END_REF][START_REF] Xue | Spring land surface and subsurface temperature anomalies and subsequent downstream late spring-summer droughts/floods in north america and east asia[END_REF]. These findings set the ground for the international initiative 'Impact of Initialized Land Surface Temperature and Snowpack on Subseasonal to Seasonal Prediction Project' (LS4P) [START_REF] Xue | Impact of initialized land surface temperature and snowpack on subseasonal to seasonal prediction project, phase I (LS4P-I): organization and experimental design[END_REF], gathering over 40 partners and more than 20 earth system models around a common experimental forecasting protocol. The phase 1 of LS4P focused on a case study, namely the impact of land surface temperature over the Tibetan Plateau (TP) in May 2003 on the subsequent drought (flood) occurrence south (north) of the Yangtze river basin in June 2003. The experimental setup consisted in adjusting the TP land temperature at the model initialization step of the subseasonal forecast systems, so as to mitigate the May 2003 surface temperature bias analyzed in a previous baseline control forecast. In this protocol, as well as in other multi-model forecast experiments implying changes in initial soil prognostic variables such as soil moisture [START_REF] Koster | The second phase of the global land-atmosphere coupling experiment: soil moisture contributions to subseasonal forecast skill[END_REF][START_REF] Ardilouze | Multimodel assessment of the impact of soil moisture initialization on mid-latitude summer predictability[END_REF], the physical consistency between soil prognostic variables at initialization has been overlooked. Inconsistencies may trigger dramatic model-dependent adjustements and partly obscure the outcome of the numerical experiments. Moreover, it prevents disentangling the respective contributions of soil temperature, soil moisture and soil enthalpy (i.e. latent energy) on land-atmosphere coupling and subsequent atmospheric response (Zhao et al, 2018). To tackle this issue, we reproduce the protocol of the LS4P prediction case study with three different initialization strategies detailed in section 2, including their impact on soil variables (subsection 2.4). Results on atmospheric forecasts are commented in section 3, before a discussion section and concluding remarks.

Data and methods

Subseasonal forecast setup

The four numerical experiments of this study each consist in a 30-member ensemble re-forecast, initialized on May 1st, 2003, and run over a 61-day forecast period, hence until June 30th. They are carried out with the CNRM-CM6-1 earth system model at the standard horizontal resolution which is approximately 150 km [START_REF] Voldoire | Evaluation of CMIP6 DECK experiments with CNRM-CM6-1[END_REF]. The land component, Interactions between the Soil Biosphere Atmosphere-CNRM Total Runoff Integrating Pathways (ISBA-CTRIP: [START_REF] Decharme | Recent changes in the ISBA-CTRIP land surface system for use in the CNRM-CM6 climate model and in global off-line hydrological applications[END_REF], is included in the Externalized Surface (SURFEX: [START_REF] Masson | The SURFEXv7.2 land and ocean surface platform for coupled or offline simulation of earth surface variables and fluxes[END_REF] model and the OASIS-MCT (Ocean Atmosphere Sea Ice Soil -Model Coupling Toolkit) coupler designed to process heat, water and momentum exchanges at the surface [START_REF] Voldoire | SURFEX v8. 0 interface with OASIS3-MCT to couple atmosphere with hydrology, ocean, waves and sea-ice models, from coastal to global scales[END_REF]. Atmospheric initial conditions are directly derived from the ERA-Interim re-analysis (Dee et al, 2011) while land initial conditions result from an offline 1993-2016 SURFEX simulation forced by 3-hourly ERA-Interim atmospheric fields, as described in [START_REF] Boisserie | Land surface initialization strategy for a global reforecast dataset[END_REF]. Similarly, the 1 • ocean and sea-ice initial conditions are derived from an ERA-Interim forced ocean run constrained towards the Mercator Ocean International Glorys 2V4 reanalysis [START_REF] Ferry | Mercator global Eddy permitting ocean reanalysis GLORYS1V1: Description and results[END_REF]. The ensemble members are generated through the addition of perturbations to the atmospheric initial conditions.These perturbations are designed as corrections of model drift errors, following the method described in [START_REF] Batté | Randomly correcting model errors in the ARPEGE-Climate v6. 1 component of CNRM-CM: applications for seasonal forecasts[END_REF].

Soil initial state adjustments

In the soil moisture perturbation experiments referred to in the introduction, neither water mass nor enthalpy of the soil are conserved during the initialization process. In this pilot study, an effort is made to use a very simple adjustment methodology to conserve these quantities as best as possible, or to at least have a more physically consistent initial state.

The enthalpy for a soil layer (J m -2 ) can be defined following (Zhao et al, 2018) as:

ℎ = [𝜌 𝑠 𝑐 𝑠 (1 -𝑤 𝑠𝑎𝑡 ) + 𝜌 𝑙 𝑐 𝑙 𝑤 𝑙 + 𝜌 𝑖 𝑐 𝑖 𝑤 𝑖 ] (𝑇 -𝑇 𝑓 ) -𝜌 𝑖 𝐿 𝑓 𝑤 𝑖 (1)
where 𝑇 represents the temperature of the soil layer (K), 𝐿 𝑓 is the specific latent heat of fusion (J kg -1 ), and 𝑤 𝑙 , 𝑤 𝑖 and 𝑤 𝑠𝑎𝑡 represent the liquid volumetric water content of the soil, the frozen water and the soil porosity, respectively (m 3 m -3 ). The densities for liquid water, ice and soil solids are represented by 𝜌 𝑙 , 𝜌 𝑖 and 𝜌 𝑠 , respectively (kg m -3 ), while 𝑐 𝑠 , 𝑐 𝑤 and 𝑐 𝑖 represent the specific heat capacity of the soil constituents, liquid water, and ice, respectively (J kg -1 K -1 ). The reference temperature corresponds to the freezing point of liquid water, 𝑇 𝑓 (273.15 K). In Eq. 1, the contribution of air in the soil pores to ℎ has been neglected.

To compute soil water phase changes, ISBA uses the relation between the soil water potential and temperature for sub-freezing conditions from [START_REF] Fuchs | An Analysis of Sensible and Latent Heat Flow in a Partially Frozen Unsaturated Soil[END_REF]:

𝜓 * = 𝜓 𝑠𝑎𝑡 + 𝐿 𝑓 𝑇 -𝑇 𝑓 𝑔 𝑇 (2)
where 𝑔 is the gravitational constant (m s -2 ), and 𝜓 𝑠𝑎𝑡 represents the soil matric potential at saturation (m). The potential 𝜓 * can be substituted in the expression for the soil matric potential in order to obtain the maximum unfrozen (liquid) water content at a given soil temperature. ISBA uses the pedotransfer functions by [START_REF] Clapp | Empirical equations for some soil hydraulic properties[END_REF] so that

𝑤 𝑙,max = 𝑤 𝑠𝑎𝑡 𝜓 * 𝜓 𝑠𝑎𝑡 -1/𝑏 (3)
where 𝑏 is related to the slope of the water retention curve. This method is described in more detail in [START_REF] Decharme | Impacts of snow and organic soils parameterization on northern Eurasian soil temperature profiles simulated by the ISBA land surface model[END_REF]. Finally, a minimum liquid water content, 𝑤 𝑙,min must be defined. Here, it is simply computed as a vegetation-cover weighted average between wilting point (corresponding to a matric potential of -300 m) and a value about half of wilting point. The latter value is an approximation of the minimum value in ISBA generally found in simulations owing to baresoil evaporation.

For LS4P, we first apply a temperature perturbation, Δ𝑇, to soil layers in a masked region. Then the soil water and ice contents are adjusted to either conserve enthalpy which corresponds to the newly imposed temperature or the soil moisture components are adjusted in order to result in a state which is closer to a thermal equilibrium (or at least what is deemed to be more physically realistic in terms of wetting-cooling or drying-warming). The four possible responses of ℎ to the application of a Δ𝑇 using the methodology herein can be summarized as follows:

-1 A potentially frozen soil thaws owing to the application of a sufficiently large Δ𝑇. Any soil ice is simply removed thus a warming soil is becoming potentially more dry. Mass is not conserved if the original state had soil ice. Note that using the definition of enthalpy in Eq. 1, ℎ conservation is not possible in this pathway as it changes sign. -2 A soil with a temperature initially above 𝑇 𝑓 cools to below freezing. As in case 1, ℎ conservation is not possible, but mass is potentially conserved: the partitioning between frozen and unfrozen soil water is made by imposing 𝑤 𝑙 ≤ 𝑤 𝑙,max where the maximum liquid water content for 𝑇 > 𝑇 𝑓 is defined using Eq. 3. -3 Both the initial and final soil temperatures are above freezing. There is no soil ice initially so that warming/cooling correspond to drying/wetting by simple inversion of Eq. 1 for 𝑤 𝑙 using the new 𝑇. Mass is not conserved, but enthalpy is unless the final adjusted 𝑤 𝑙 value either exceeds saturation, 𝑤 𝑠𝑎𝑡 , or falls below some minimum value, 𝑤 𝑙,min (described in more detail below). -4 Both the initial and final soil temperatures are below freezing. In a first step, the total mass is conserved but the partitioning between ice and liquid soil water is modified to conserve enthalpy. The ice content is first diagnosed under this assumption by solving Eq. 1 for 𝑤 𝑖 . The soil ice is then checked to see if it exceeds physical limits (totally frozen soil or negative soil ice). In this case, the ice is limited and liquid water is adjusted to conserve enthalpy (and mass). As in case 3, if after adjustment the liquid water exceeds physical limits, it is limited thus both mass and enthalpy conservation imbalances can occur.

These pathways are represented by arrows in the left-hand panel of the schematic shown in Fig. 1. Indeed, when the imposed temperature change Δ𝑇 leads to the crossing of the freezing temperature threshold 𝑇 𝑓 (meaning going from sub-freezing temperatures to freezing or vice versa), the soil enthalpy can no longer be conserved using the definition of enthalpy in Eq. 1 (i.e. using a reference temperature of 𝑇 𝑓 ) together with the fact that we are attempting to preserve the temperature anomalies. Note that in reality, latent heat changes (going between frozen and unfrozen soil states) would balance input heating or cooling. The Table in Fig. 1 summarizes the four possible pathways for the LS4P phy experiment and which quantities are conserved (or not).

Soil temperature Soil liquid water Soil ice

Ctrl ✗ ✗ ✗ LS4P t ✓ ✗ ✗ LS4P sm ✗ ✓ ✓ LS4P phy ✓ ✓ ✓
Table 1: List of experiments and prognostic variables which are potentially adjusted at initialization over the TP.

Sensitivity experiments and reference data

The baseline (Ctrl) experiment corresponds to the setup described in Sec. 2.1, and the 3 additional experiments only differ from the Ctrl case through changes in land initial conditions over the TP. In LS4P t , a temperature correction mask (Fig. 2c) is applied to the 14 soil layers of ISBA-CTRIP, corresponding to the default 10-meter soil depth.

The underlying idea of this correction is to reduce the average cold bias of the model affecting the TP. This positive temperature mask was generated following the method described in the section 3.2.3 of [START_REF] Xue | Impact of initialized land surface temperature and snowpack on subseasonal to seasonal prediction project, phase I (LS4P-I): organization and experimental design[END_REF]. It was calculated by taking into account both the average cold bias over the entire TP (Fig. 2b) and the cold observed anomaly (Fig. 2a). The LS4P phy experiment consists in adjusting soil liquid water and ice content while applying the prescribed temperature correction, in order to impose physical consistency and/or enthalpy conservation, to the maximum possible extent.

Finally, the last sensitivity experiment, LS4P sm , consists in adjusting the soil moisture over the TP as in LS4P phy , but without applying the soil temperature correction. For example, soil is dried (wetted) if a positive (negative) temperature change was prescribed by LS4P. Thus, this experiment is somewhat akin to a GLACE type experiment in that only the soil water content is modified, but in a manner consistent with the temperature corrections (thus differing from GLACE in this respect). The aim of LS4P sm is to provide insight on the respective contributions of the initial soil temperature and soil moisture over the TP in subseasonal forecasts. Table 1 recapitulates, for each experiment, the soil prognostic variables that are changed at initialization. In this study, three reference datasets are used to evaluate our simulations. The 2meter temperature (t2m) is assessed against version 4 of the monthly gridded CRU-TS dataset [START_REF] Harris | Version 4 of the CRU TS monthly high-resolution gridded multivariate climate dataset[END_REF] and precipitation is evaluated using the monthly multisource MSWEP V2 [START_REF] Beck | MSWEP V2 global 3-hourly 0.1 precipitation: methodology and quantitative assessment[END_REF] dataset. The ERA5 reanalysis [START_REF] Hersbach | The ERA5 global reanalysis[END_REF] is used as the reference for atmospheric fields on standard pressure levels. Simulation output and reference data have been re-gridded onto a common 1x1 • regular rectangular grid.

The comparison of 2003 subseasonal re-forecasts with observational counterparts is based on anomalies relative to the 1993-2016 period rather than the 2003 full fields, which is a standard bias correction strategy. In that respect, the model climatology is computed from an extension of the Ctrl re-forecast to all the other years of the 1993-2016 period, with the same initialization day (May 1st), and following the setup described in Sec. 2.1.

Mass conservation

Enthalpy conservation If adjustments are indeed required (for example, to avoid super saturation or overdrying resulting from the first pass at enthalpy conservation), enthalpy might not be conserved.

1 ✓ ✗ 2 ✓ ✗ 3 ✗ ✓ 4 ✓ ✓

Impact on land prognostic variables

To verify the validity of our setup, we first analyze the vertical profile of land prognostic variables after one day of model integration. In Fig. 3, only the profile for the top first meter of soil is displayed since it has an impact on land-atmosphere interactions over the S2S timescale. Although initial soil temperature is unchanged in LS4P sm (green line) with respect to the Ctrl (black), the temperature profile appears slightly shifted towards higher values (Fig. 3c) due to the latent heat release owing to soil water freezing occurring during the first steps of the model integration. This is supported by the slightly higher ice content and lower liquid water content profiles (Fig.s 3b and 3a, respectively) in LS4P sm , as compared to LS4P phy (orange line), despite their identical change in initial soil water content. The difference in the temperature profile ). Thus, there is a large physical adjustment owing to an initial state which is not in thermal equilibrium. Note that the four temperature profiles are relatively close to each other near the surface, because of the rapid adjustment to the overlying atmospheric conditions after the first day of model integration. Not surprisingly, the evolution of the land prognostic state variables throughout the 2-month forecast time (Fig. 4) shows a progressive decrease of the soil ice content in the uppermost m of the soil for the Ctrl and LS4P t , in agreement with the seasonal atmospheric warming occurring from late spring to early summer (Fig. 4b). In response to ice thawing, soil liquid water in LS4P t and Ctrl progressively increases with time (Fig. 4a). While the soil moisture content remains nearly identical for LS4P sm and LS4P phy , the latter is substantially warmer at all times, despite a decreasing trend during the first 15 days (Fig. 4c). LS4P sm starts as cool as in Ctrl, by construction, but ends up slightly warmer than LS4P t . This relates to the abundant soil ice in the However, although the soil becomes progressively warmer than LS4P t , the seasonal decline of the snow cover does not catch up, because the amount of accumulated snow at the beginning of the melting season (∼ May 20th) is higher than in the other sensitivity experiments. Note that in the current study, the soil states are modified while holding the snow properties (mass, snow pack enthalpy, liquid content...) constant between the original and perturbed setups. But owing to the feedbacks between the soil and snow pack, we plan to investigate perturbing the snow pack properties along with the soil in the next phase of LS4P.

Results

In this section, we analyze the atmospheric response to the land initialization strategy over the TP. For each simulation, the ensemble mean is considered. The statistical sig-nificance of the difference between two simulations is assessed by means of a two-sided Student's t-test for paired samples. Here, the samples correspond to the 30 ensemble members of the simulations. The null hypothesis is rejected for p-values exceeding 0.05. The comparison of the simulations with observational references is based on pattern anomaly correlation calculated over the entire northern hemisphere and 3 sub-regions: North America (140

• W-60 • W,25 • N-60 • N), Europe (15 • W-45 • E,30 • N- 70 • N) and East Asia (80 • E-140 • E,20 • N-55 • N).
The boundaries of these regions are depicted in Fig. 5(a-1). The observed anomalies are relative to the observed climatology over the 1993-2016 period, while the simulated anomalies are relative to a model climatology, derived from a full hindcast initialized on May 1st for each year of the 1993-2016 period, following the setup described in Sec. 2.1. This manner of comparing anomalies rather than full fields is a standard bias-correction procedure in the evaluation of subseasonal forecasts (e.g. [START_REF] Dirmeyer | Land-atmosphere interactions and the water cycle. Seamless prediction of the Earth system: from minutes to months[END_REF]. Note that for assessing the statistical significance of pattern correlations, we take into account the spatial autocorrelation between adjoining grid points [START_REF] Wilks | Statistical methods in the atmospheric sciences[END_REF]. As a consequence, the effective sample size of grid points retained to compute the the t-statistic is higher for spatially heterogeneous fields (precipitation), than for smoother fields (temperature or geopotential height). Thus, for similar correlation values, the precipitation pattern correlation is more likely to be significant than the geopotential height at 500 hPa (z500) or temperature counterparts (see Table 2).

Evaluation of the Ctrl forecast for May and June 2003

The Ctrl re-forecast simulation is compared to observations for t2m, precipitation and mid-troposphere circulation patterns (Fig. 5). As expected, the magnitude of the June 2003 Ctrl anomaly patterns (sub-figures a-4, b-4 and c-4) are weaker than the May counterpart (a-3, b-3 and c-3) : indeed, due to the chaotic nature of the atmosphere, the forecast members tend to diverge with time and therefore the ensemble spread is larger in June than in May. The warm anomaly observed in May and June over Europe corresponds to the early stage of the extremely hot summer of 2003 [START_REF] Black | Factors contributing to the summer 2003 european heatwave[END_REF]. The anomaly pattern correlations reported in Table 2 confirm that Ctrl compares remarkably well with observations in May, while to a lesser extent in June, where a spurious cool anomaly is forecast over the Iberian peninsula. The wet-North dry-South pattern, as well as the southern Europe anticyclonic blocking, are also well captured by the model in May, resulting in very high correlations for precipitation and z500. In June, a weaker, but significant signal is still found over central Europe for precipitation. More generally, the correlations drop substantially in June for most of the considered variables and regions, which is not surprising, given that the May correlations are augmented by the higher predictability of the first 10 days after initialization belonging to the meteorological predictability horizon. when Ctrl succeeds in simulating the positive anomaly near the West coast of North America (Fig. 5c-4).

Impact of the TP land initialization strategy

In May 2003 (Fig. 6a), the three perturbed simulations show a significant impact on t2m over North America. While they all tend to amplify the Ctrl warm bias from the Hudson Bay to New England, only LS4P sm and LS4P phy mitigate the cold bias over the western US. Additionally, LS4P phy simulates a cooler pattern over much of the south-east US affected by a pronounced warm bias in Ctrl. This is confirmed by the increased pattern correlation in Figure 7a. Elsewhere, none of the simulations particularly outperforms the other ones except for LS4P t which improves temperature patterns over East Europe and West Russia. At the northern hemisphere scale, all the simulations are fairly even in terms of statistically-based performance. In June 2003 (Fig. 6b), the signal is weaker over America, and the simulations do not reduce the strong warm bias over Central US, a classic summer pattern in many climate models (e.g. Zhang et al, 2018). The bias is even accentuated in LS4P t and LS4P phy to a lesser extent. Over Eurasia, the three simulations show warmer t2m north east of the Caspian Sea but only LS4P phy partly compensates the cool bias over the Iberian peninsula and Central Siberia. This simulation also performs better over the TP, North China and Lapland. The t2m pattern correlation (Fig. 7b) further underscores the good performance of LS4P phy over these regions in June.

A similar analysis for precipitation also indicates a strong impact of our experimental setup in May 2003 over North America (Fig. 8a), with LS4P sm and LS4P phy showing enhanced mitigation of the wet bias over the Rockies and Great Plains. In June 2003, the signal is weak overall (Fig. 8b). However, LS4P phy indicates several local improvements. More precisely, this simulation captures a dryer pattern over central Siberia, and a wetter pattern over northern India, western Sahel and north of the Yangtze, keeping in mind that the latter regions are affected by a strong dry bias in Ctrl. The markedly wet bias over Myanmar is mitigated in both LSP t and LS4P phy . For all of the large-scale focus regions, the pattern correlations do not reveal any substantial difference between simulations (not shown). However, the aforementioned regional improvements in LS4P phy suggest a rather beneficial impact of this land initialization strategy also for the precipitation forecast. The local improvements captured nearby the Yangtze river are further commented at the end of Sec. 3.3. The t2m and precipitation patterns described above are partly driven by changes in the atmospheric circulation The maps for z500 (Fig. 9) reveal differences with a larger amplitude in May, especially over the western north hemisphere. Differences are weaker for LS4P t , though. The LS4P phy map suggests a wave-like pattern. This feature will be commented in the following section. For this simulation, the regional correlations indicate substantially improved circulation anomalies across North America (Fig. 10a), but no particular impact elsewhere.

In June (Fig. 9b), the sensitivity experiments exhibit overall less widespread differences with respect to Ctrl. However, here again, LS4P phy distinguishes itself from the other simulations owing to a relatively marked geopotential anomaly tripole composed by a central positive anomaly stretching across southern Siberia and two negative anomalies over Northern Europe and, more weakly, from northern China to north-west Pacific. This feature is also visible in LS4P sm and LS4P t , but with barely, if at all, any significant patterns. The pattern correlations reveal that LS4P phy considerably improves the circulation anomaly over both Europe, thanks to lower geopotential height over Scandinavia, and East Asia. This result is consistent with the t2m analysis.

The relatively higher performance of the Ctrl simulation over North America in June is largely due to the successful simulation of the z500 positive anomaly located near the west coast of America. A possible reason for the poorer results of the sensitivity experiments in this respect is discussed in Sec. 4.

Description of the mechanisms involved

In the previous section, we identified a number of common features between the sensitivity experiments: the adjustement of land initial conditions over the TP affects mainly the western half of the northern hemisphere during the first month of the As suggested earlier, a striking feature of the change in z500 circulation over North America in May 2003 for LS4P phy is a wave-like pattern. Additionally, it is remarkably out-of-phase from the Ctrl bias pattern (the pattern correlation value between the z500 bias and the LS4P phy z500 anomaly over the North America domain is -0.74). This is also true for LS4P sm but to a lesser extent, especially for south-East US (spatial correlation value : -0.48). The anomaly pattern found in LS4P phy matches very well with the Tibetan Plateau-Rocky Mountains Circumglobal (TRC) wave train identified by Xue et al (2022) (Their Fig. 4B) and named as such in Xue et al (2023, this issue). In June (Fig. 9b), the relative agreement between sensitivity simulations on a geopotential anomaly tripole described earlier and stretching from Northern Europe to Northern China suggests that the TRC, identified as a May teleconnection pattern, can propagate in June across Eurasia. In the upper troposphere, the geopotential height anomaly is affected consistently by the mid-troposphere (Fig. 11), and consequently, the jet stream is shifted northwards in May over North America and in June over Asia, with a stronger amplitude in LS4P phy .

In order to illustrate the dynamic propagation, we show the time evolution of 250 hPa meridional wind anomalies from Ctrl, averaged out over the 30 • N-60 • N meridional cross-section (Fig. 12). The three perturbed simulations clearly show a wave train emerging from the TP longitudes and propagating eastwards during the first 30 days. This wave train is very weak in LS4P t and rapidly fades out while on the other hand, it is particularly strong in LS4P phy . The attenuated wave train re-emerges past day 40 in longitudes stretching eastwards from Central Asia to the North Pacific. Again, it is the strongest and best established in LSP phy . The choice of a cross section bounded between 30 • N and 60 • N was a trade-off to include most of the TP and mid-latitudes, without averaging out too many opposite anomalies, that could mask the signal. The drawback is that the wave propagation at higher latitude across northern Scandinavia and Lapland, if any, is not detectable. This could explain the discontinuity found in wave propagation illustrated in the Hovmöller diagram between day 25 and day 40 after initialization. A similar Hovmöller analysis has thus been achieved for a cross section between 55 • N and 75 • N (Fig. 13). It reveals an eastward wave propagation between day 25 and day 45 after initialization from the Greenwich meridian to 80 • E. This wave propagation is more pronounced in LS4P phy than in other perturbed experiments.

The way land surface is initialized over the TP in early May is therefore susceptible to trigger an extratropical circumglobal atmospheric wave train. This wave train considerably modifies the atmospheric flow in our model for several weeks.

Finally, even if the impact of our experiments on continental south East Asia precipitation remains inconclusive, the potential interaction with the onset and evolution of the East-Asian monsoon cannot be discarded, and in particular one of its main components known as the Meiyu-Baiu front. This quasi-stationary area of intense rainfall forms in spring at the confluence of the summer monsoon moist air south of China and the warm dryer air to the North (e.g. Zhou et al, 2004). The location and evolution of the Meiyu-Baiu front is governed by a wide range of factors covering multiple scales, including the thermal forcing of the TP (e.g. [START_REF] Liu | Impact of Tibetan orography and heating on the summer flow over Asia[END_REF], as well as the mid-latitude westerly jet [START_REF] Sampe | Large-scale dynamics of the meiyu-baiu rainband: Environmental forcing by the westerly jet[END_REF]. In this study, both factors are impacted by the TP land initialization, and could thus change the location and intensity of this front. In Fig. 14, we show the anomaly of vertically integrated moisture divergence in ERA5 and Ctrl, in order to localize the Meiyu-Baiu front. In May, Ctrl captures relatively well the divergence (convergence) anomalies in red (blue) shades, corresponding to dry (wet) spells, especially over the China sea. The strip of moisture convergence simulated south of Japan is mislocated and too zonally tilted with respect to observations. The location is partly improved in LS4P phy and LS4P sm with a south-eastward shift of this convergence front (Fig. 15a). This shift results from the change in large-scale circulation, and in particular the meridional displacement of the subtropical westerly jet [START_REF] Sampe | Large-scale dynamics of the meiyu-baiu rainband: Environmental forcing by the westerly jet[END_REF]. More precisely, in these two simulations, the wave propagation triggers the development of a barotropic low pressure anomaly between China and Southern Japan at 500 hPa (fig. 9a), and to a lesser extent at 250 hPa, in the lee of the TP (two bottom-left maps of fig. 11). Accordingly, the subtropical westerly jet is accelerated south of this low anomaly, which favors the southwards advection of cool dry air from the mid-latitudes. Eventually, the convergence with northward monsoonal flow is shifted further south than in the Ctrl experiment. This result is consistent with the study from [START_REF] Volonté | Characterising the interaction of tropical and extratropical air masses controlling East Asian summer monsoon progression using a novel frontal detection approach[END_REF] highlighting the sensitivity of the Meiyu-Baiu front location to the transient mid-latitude weather disturbances steered by the westerly jet.

In June (Fig. 15b), the observed moisture divergence anomaly is rather patchy, suggesting a little organised Meiyu-Baiu front. Over China, the dry (wet) anomalies located south (north) of the Yangtze river materialize the exceptional drought (precipitation) that motivated the choice of the year 2003 as the initial case study of LS4P [START_REF] Xue | Impact of initialized land surface temperature and snowpack on subseasonal to seasonal prediction project, phase I (LS4P-I): organization and experimental design[END_REF]. Ctrl misses these anomalies and simulates a convergence anomaly north-east of the Philippines, where observations indicate a anomalously dry pattern (red shades). This is partly improved in LS4P t and LS4P phy where dryer conditions are simulated. It is noteworthy that LS4P phy simulates a strip of enhanced moisture convergence north of the Yangtze and conversely further south. This could be related to the improved z500 circulation pointed out in Sec. 3.2, but considering the small scale of these features, verification of any relationship would require further analysis beyond the scope of the current study.

Discussion and conclusion

We have revisited the LS4P exercise [START_REF] Xue | Impact of initialized land surface temperature and snowpack on subseasonal to seasonal prediction project, phase I (LS4P-I): organization and experimental design[END_REF] which consists in adjusting the land initialization over the Tibetan Plateau (TP) in ESM-based late spring subseasonal forecasts. The rationale of LS4P is to mitigate the initial subsurface temperature bias affecting this high altitude region in order to make the most out of this source of extended range atmospheric predictability. In this pilot study focusing on a 2-month forecast initialized on May 1st 2003, we compared three initial adjustment strategies over the TP to assess the impact of the physical equilibrium of land initial conditions.

Adjusting the land initialization over the TP with more thermodynamicallybalanced state variables leads to marked changes in the subsequent atmospheric response. In particular, the strategy that adjusts the soil liquid water and ice content while applying the prescribed temperature correction triggers the most effective and long-lasting atmospheric response at the hemispheric scale among the three setups. This setup also better preserves the soil physical consistency and/or enthalpy at initialization. The atmospheric response translates into a an extratropical tropospheric wave train emerging from the TP propagating eastwards towards north America during the first month, and more weakly to Europe and Asia during the second month. In this study, our balanced initialization strategy improves the atmospheric circulation pattern over North America in May and Europe and East-Asia in June, with respect to an observational reference. It also interacts with the East-Asian monsoon flow by changing the location and intensity of the Meiyu-Baiu front but without any noticeable improvement.

Conversely, we find that neither the correction of initial land temperature alone nor that of initial soil moisture alone has a major impact on the atmospheric response within our model. Only the joint correction of soil moisture and temperature with physical coherence succeeds in improving the forecast in this case study. A coordinated multi-model exercise such as GLACE-2 or LS4P, but focused on this issue of land physical consistency at initialization, would be required to verify if our findings can be generalized to other models. Similarly, this paper only reports on a single case study, and thus a way forward could be to extend this exercise over multiple consecutive years, with an ad hoc temperature correction mask for the TP each year. Such a multi-year reforecast is needed to evaluate the screen-level temperature and circulation prediction skill in different regions of the Northern hemisphere. Yet, our work clearly suggests the importance of preserving the land physical equilibrium in coupled land-atmosphere numerical experiments. This question is especially acute in subseasonal and seasonal dynamical forecast systems where the land component is initialized in a fashion to approach observed states. This implies a direct or indirect data assimilation system fed by superficial observations. Therefore, a joint assimilation of both observed soil temperature and moisture while imposing a thermally-consistent state upon initialization of the land surface could lead to improved coupled subseasonal to seasonal forecasts. Finally, our numerical simulations have been carried out with an interactive ocean model in order to mimick the setup of the majority of S2S forecast systems [START_REF] Vitart | The subseasonal to seasonal (S2S) prediction project database[END_REF]. However, it is an exception among other contributing models in the LS4P exercise. A similar set of experiments with prescribed sea-surface temperature would be worth implementing to verify the robustness of our conclusions and to better compare our results with those of other modeling systems. The limited success of the June 2003 forecast over North America in the sensitivity experiments (Sec. 3.2) could originate from a poor simulation of the seesaw atmospheric linkage between the Tibetan plateau and the Rocky Mountains, as identified by [START_REF] Xue | Spring land surface and subsurface temperature anomalies and subsequent downstream late spring-summer droughts/floods in north america and east asia[END_REF]. In the second phase of the LS4P project, we plan to apply our new initialization method over both mountain ranges simultaneously, to try and improve the model atmospheric response.

Considering the increasing refinement of the physical processes simulated in land surface models, our study makes the case for particular care to be taken in the initialization of the land component in weather-to-seasonal forecast systems, preserving the land physical consistency to the best possible extent. Additionally, this works sheds light on the soil enthalpy which may well be a serious candidate to help reconcile the tug of war between supporters of land subsurface temperature versus root zone soil moisture as the dominant driver stemming from land in subseasonal forecasts.

Fig. 1 :

 1 Fig.1: Schematic of the four possible scenarios for soil enthalpy, ℎ, change when initial soil temperature 𝑇 is modified through the application of a temperature change in LS4P phy . 𝑇 𝑓 represents the freezing point temperature. The dashed quadrants are deemed to be physically impossible. The physically-based limits prevent enthalpy conservation in cases 1 and 2: liquid water mass is conserved for case 1 (ice is removed if initially present), and total mass is conserved for Case 2 unless physical limits of the final adjusted moisture state are exceeded. Cases 3 and 4 allow enthalpy conservation with adjustments (to soil water content) if physical limits are exceeded. If adjustments are indeed required (for example, to avoid super saturation or overdrying resulting from the first pass at enthalpy conservation), enthalpy might not be conserved.

Fig. 2 :Fig. 3 :

 23 Fig. 2: (a) May 2003 observed monthly mean 2-m temperature anomaly, (b) difference between the Ctrl ensemble mean and the observation, (c) and soil temperature correction mask applied at initialisation. Unit: K

Fig. 4 :

 4 Fig. 4: Ensemble mean evolution of top 1 meter liquid water (a), ice (b), temperature (c) and snow water equivalent departure from the initial Ctrl value (d) over the TP

Fig. 5 :Fig. 8 :

 58 Fig. 5: Comparison of May and June 2003 observed and Ctrl anomalies for (a) 2-meter temperature in K, (b) precipitation in mm/day and (c) 500 hPa geopotential height (m)

Fig. 9 :

 9 Fig. 9: As figure 6 for geopotential height at 500 hPa (unit: m).

Fig. 11 :

 11 Fig. 11: May (left column) and June (right column) 2003 250 hPa geopotential height (shade) and zonal wind (contour) anomaly

Fig. 12 :

 12 Fig. 12: Hovmöller diagrams of 250 hPa meridional wind anomaly averaged between 30 • N and 60 • N. Values not significant at the 95% have been set to 0.

Fig. 13 :

 13 Fig. 13: As figure 12 between 55 • N and 75 • N.

Fig. 14 :

 14 Fig. 14: Comparison of May and June 2003 observed and Ctrl anomalies for vertically integrated moisture divergence in kg m -2 day -1

  Fig. 15: (a) May and (b) June 2003 vertically integrated moisture divergence (shade) in kg.m -2 .day -1 and 850 hPa wind anomalies (vectors).

  Over the rest of the northern hemisphere, the temperature patterns are less well captured by Ctrl with a few local exceptions like over north-west America, Asian steppes and western Africa in May. The shift from dry May to wet June anomalies north of the Gulf of Guinea is very well simulated by Ctrl. Finally, the z500 patterns are very well captured over Europe and East Asia in May but less so over North America, and conversely in June,

		2-m temperature	precipitation	z500
	Northern hemisphere	M:0.44 J:0.31	M:0.23 J:0.17	M:0.49 J:0.23
	North America	M:0.42 J:-0.04	M:0.28 J:0.14	M:0.29 J:0.46
	Europe	M:0.59 J:0.28	M:0.62 J:0.23	M:0.89 J:-0.1
	East Asia	M:0.50 J:0.08	M:0.14 J:0.16	M:0.81 J:0.03
	Table 2: Spatial anomaly correlation of Ctrl with Obs in May (M) and June (J) 2003.
	The bold font highlights positive Values which are significant at the 95% confidence
	level.

sea surface temperature effect -The GEWEX/LS4P Phase I experiment. Climate Dynamics Zhang C, Xie S, Klein SA, Ma Hy, Tang S, Van Both authors contributed to the study conception and design. Material preparation, data collection and analysis were performed by Constantin Ardilouze. Both authors wrote the first draft of the manuscript and commented on previous versions of the manuscript. Both authors read and approved the final manuscript.

ERA5 data (Hersbach, H. et al., 2018) was downloaded from the Copernicus Climate Change Service (C3S) (2023). The datasets generated during the current study are not publicly available but are available from the corresponding author on reasonable request