C N Abadie

B W Byrne

G T Houlsby

MODELLING OF MONOPILE RESPONSE TO CYCLIC LATERAL LOADING IN SAND

HAL is

Introduction

Offshore wind turbines are large structures, for which the foundation is critical and costly. The most common foundation is the monopile, which is a single stiff tubular pile of up to 8m diameter and perhaps 45m long. During their lifetime, monopiles are exposed to a range of loads of very many cycles. Because of an inadequate understanding of the effects of these loading conditions (e.g. DNV, 2014), currently installed monopiles may be over-designed [START_REF] Kallehave | Modification of the API formulation of initial stiffness of sand[END_REF], causing excessive manufacturing, transportation and installation costs. While significant progress has been achieved on the design of monopiles for ultimate capacity [START_REF] Byrne | New design methods for large diameter piles under lateral loading for offshore wind applications[END_REF], understanding and modelling the effects of cyclic loading is still developing. It is nonetheless recognized that, as the monopiles are made shorter and with smaller diameters, following optimization, cyclic loading could cause progressive tilting of the foundation along with changes to stiffness, damping and capacity. This paper outlines experimental and theoretical research, aimed at providing further insight, and a new modelling technique (HARM -"Hyperplastic Accelerated Ratcheting Model") for piles in sand subjected to lateral cyclic loading. The experimental work highlights the important effects that need to be accounted for in detailed modelling. The theoretical development then elaborates on these key findings and leads to a constitutive framework that describes the soil-pile interaction response in a rigorous and concise manner. The final part of this paper demonstrates how the experimental results are used as source data for calibration of the model and demonstrates the suitability of the methodology for complex load scenarios with large cycle number.

Model test experiments 2.1 Background and motivation

Pile response to lateral loading is commonly assessed using the so-called p-y method (API, 2010, DNV, 2014). The soil-pile interaction is represented by a series of independent springs down the pile:

  y y x E Ap y K Ap p py u u            , tanh 0 , 0
Where p is the local soil reaction at a specific depth x and y is the corresponding pile deflection. A is an empirical parameter, pu is the ultimate lateral pile resistance and K0 is the initial soil reaction.

When considering the case of cyclic loading, it is recommended (API, 2010, DNV, 2014) to reduce the parameter A to a constant value (A=0.9 for sands). This simple method was originally derived from investigation on long slender piles used in the oil and gas industry [START_REF] O'neill | An evaluation of p-y relationships in cohesionless soils[END_REF] for jacket structures, for which down-rating the ultimate lateral capacity, following cyclic loading, is appropriate. It is questionable whether this approach, which provides the same predicted response irrespective of the number of cycles applied, is relevant to offshore wind monopiles. Accordingly, alternative approaches have been proposed recently in the literature, often making use of small-scale experimental tests. The pile deflection at peak load is measured and a direct fitting of the results provides an empirical relationship between pile deflection and cycle number (e.g. LeBlanc et al., 2010b[START_REF] Peralta | Investigations on the Behavior of Large Diameter Piles under Long-Term Lateral Cyclic Loading in Cohesionless Soil[END_REF][START_REF] Klinkvort | Centrifuge modelling of drained lateral pile-soil response[END_REF].

These laws can be applied to the monotonic response (e.g. LeBlanc et al., 2010b), or can be integrated within Equation 1to "degrade" the monotonic p-y modulus with cycle number:

   N E E py N py 0 , ,
 is an empirical parameter directly derived from the fitted empirical law [START_REF] Abadie | Cyclic Lateral Loading of Monopile Foundations in Cohesionless Soils[END_REF].

Despite providing a fast methodology for estimating the cyclic pile deflection, these semi-empirical models lack rigorous description of the fundamental behaviour of the soil-pile response. They do not account for change in capacity, damping or stiffness, and indeed could be misleading for the latter, as they predict a degraded soil modulus, while increased foundation stiffness has been observed (e.g. LeBlanc et al., 2010b[START_REF] Klinkvort | Centrifuge modelling of drained lateral pile-soil response[END_REF][START_REF] Abadie | Cyclic loading response of monopile foundations in cohesionless soils[END_REF]. Finally, they are difficult to adapt for load packets at different amplitudes and, to do so, requires strong assumptions such as linear cumulative damage [START_REF] Miner | Cumulative damage in fatigue[END_REF]) (e.g. LeBlanc et al., 2010a). To address this issue, an experimental investigation was performed, analysing qualitatively the key mechanisms driving the global pile response.

Experimental techniques

The experimental tests were performed using a rig developed by [START_REF] Rovere | Cyclic loading test machine for caisson suction foundations[END_REF], of which the set-up is fully described in [START_REF] Abadie | Cyclic Lateral Loading of Monopile Foundations in Cohesionless Soils[END_REF]. It involved a small rigid pile, of 77mm diameter and 360mm length. A thorough consideration of scaling, applicable to this set of tests, including pile dimensions and stress level, is described in LeBlanc et al. (2010b), [START_REF] Abadie | Cyclic Lateral Loading of Monopile Foundations in Cohesionless Soils[END_REF] and [START_REF] Abadie | Cyclic plastic response of stiff piles in sand. Submitted to Géotechnique. API[END_REF]. The load was applied at an appropriate eccentricity above ground level using a suspended mass system and a driving motor of frequency 0.106Hz (Figure 1). The pile was driven into loose dry Yellow Leighton Buzzard sand samples.

A load cell and a pair of LVDTs enabled accurate recording of the macro-response of the pile. In the following, the test results are shown as lateral load (H) vs. displacement at the point of load application (yt), as defined in Figure 1. For comparison with larger scale piles, both parameters are normalised using the corresponding values at ultimate capacity, designated by the subscript R and established in [START_REF] Abadie | Cyclic Lateral Loading of Monopile Foundations in Cohesionless Soils[END_REF]. For consistency with Section 3, the load ratio is called  and the displacement ratio : The tests are described in Table 1, with the load scenarios specified in terms of number of cycles N, maximum load magnitude b and amplitude c:

max max     R b H H max min max min      H H c
Even though c is not zero, it is sufficiently low that the cyclic tests are considered as 1-way loading.

Monotonic response

The initial monotonic response, called the backbone curve, is observed through a simple load-reversal test (H0). The test result (Figure 2(a)) proves a hysteresis loop that conforms closely to the Masing rules [START_REF] Masing | Eiganspannungen und Verfestigung beim Messing[END_REF]: the reversed loading (unload [START_REF] Pyke | Nonlinear soil models for irregular cyclic loadings[END_REF]: (i) when exceeding the maximum load, the unloading or reloading curve follows that of the backbone curve and (ii) each time an unloading or reloading path intersects a curve from a previous cycle, the stress-stain response follows that curve.

Continuous cyclic loading

This is extended to long-term cyclic loading at nonzero mean load, for which a typical response is displayed in Figure 3 (test CMLT3). The graph shows accumulation of permanent deformation with cycle number that, despite reducing, does not decay to zero. This phenomenon is called ratcheting.

The graph also shows a tightening of the hysteresis loop during cycling, with an increase in stiffness and decrease in loop area. This effect was quantified by [START_REF] Abadie | Cyclic Lateral Loading of Monopile Foundations in Cohesionless Soils[END_REF] and [START_REF] Abadie | Cyclic plastic response of stiff piles in sand. Submitted to Géotechnique. API[END_REF], and is of second order compared to ratcheting. For instance, the secant stiffness typically increases by only 15% between 100 and 10,000 cycles.

Although the tests were performed in loose dry sand, for which larger ratcheting deformation might be expected, the underlying response of Figure 3 is in good agreement with other test results from the literature, performed at higher relative densities and

Multi-amplitude cyclic loading

The effect of cyclic loading on capacity is explored through a series of cyclic loading tests of the same load magnitude but different cycle number, immediately followed by a monotonic reloading (CMC0-5, CMLT2, Figure 4). For clarity, the past cyclic load history is not plotted on the graph but is indicated by the dashed line at maximum peak load. Each reloading curve is plotted from the final cyclic deformation to put each into context. The graph shows that, despite the effects of ratcheting, on exceeding the maximum cyclic load, the curve tends to rejoin the initial backbone curve. This conforms to the second observation made by [START_REF] Pyke | Nonlinear soil models for irregular cyclic loadings[END_REF].

Finally, a series of cyclic load packets at different magnitudes was applied, with a typical loaddeflection response displayed in Figure 5. The standard backbone curve is also plotted. The results conform to the above, showing that the key features characterising the pile response are: (i) For each load packet, on initial reloading exceeding the maximum past cyclic load, the response tends towards re-joining the backbone curve. This means that the extended Masing rules still apply after cyclic loading. (ii) However, this mechanism is competing with ratcheting, which moves the global response away from the backbone curve.

Constitutive modelling 3.1 Model description

The model described here is derived through the hyperplasticity framework, a novel approach for modelling plasticity behaviour that is more concise and rigorous than the standard methods widely used to date. The mathematical background of the approach is thoroughly described in [START_REF] Houlsby | Principles of Hyperplasticity: an Approach to Plasticity Theory Based on Thermodynamic Principles[END_REF] and full details of the model presented below are provided in [START_REF] Houlsby | A model for nonlinear hysteretic and ratcheting behaviour[END_REF]. For consistency, the notation adopted here are the same as in [START_REF] Houlsby | A model for nonlinear hysteretic and ratcheting behaviour[END_REF]. The model is expressed in terms of stress ( and strain ( for application to the experimental results through Equations 1 and 2.

The model provides a 0-D representation of the soilpile macro-response outlined experimentally, which could be used in a structural analysis of a wind turbine. This justifies the consideration of H and y at the location of load application as this bypasses the need for a two-variable model to capture both rotation and displacement at ground level. [START_REF] Masing | Eiganspannungen und Verfestigung beim Messing[END_REF] and [START_REF] Pyke | Nonlinear soil models for irregular cyclic loadings[END_REF] showed that models employing pure kinematic hardening give rise to the extended Masing rules. Hence, the starting point of the model developed below is kinematic hardening.

A typical multi-surface kinematic hardening model is described using the schematic representation of Figure 6(a) [START_REF] Houlsby | Principles of Hyperplasticity: an Approach to Plasticity Theory Based on Thermodynamic Principles[END_REF]. On loading the system, the springs (Hn) of each individual unit contribute to the global response when the corresponding slider reaches the threshold value (kn). The typical response is shown in Figure 6(b). A choice of a large number of spring-slider units, or yield surfaces NS, enables computation of a smooth continuous curve (Figure 6(c)). For a load-controlled case, the incremental strain response would be defined according to: where the incremental plastic strains are defined according to the family of yield surfaces:

    N n n d E d d 1   
  n n n S k H N n       , 1
As illustrated in Figure 6(c), this model provides a stabilised loop when predicting cyclic loading. This is the case of most standard models and further development is required to implement the second key feature of the monopile response: ratcheting.

Using this base model, a new ratcheting element is added (Figure 7(a)). This creates an additional increment of strain dr for each increment of plastic strain, in the direction of loading:

   N n n n r d R S d 1 ) (   
The amount of ratcheting strain generated by each hardening surface is characterised by the ratcheting rate Rn within that surface and now:

r N n n d d E d d         1
with the increment of plastic strain still defined by Equation 6. This modified model, called HARM, provides the incremental response of Figure 7(b), featuring both key mechanisms of Section 2.5.

Computation of large number of cycles

To model large cycle numbers (typically for monopiles, above 10 8 cycles), incremental calculations such as the above are computationally prohibitive. However, because of the definition of the ratcheting strain, it is possible to accelerate the effects of ratcheting by multiplying the rate Rn by a factor Rfac corresponding to the number of cycles that are accelerated through. This is illustrated in Figure 7(b), where the programme of Table 2 was used. This is a powerful feature of the model, which is demonstrated on small numbers of cycles for clarity, but can be extended to a larger number of cycles (see later Table 3).

Application to Test Results

Kinematic hardening variables: backbone curve

The values for the stiffness Hn and strength kn are calibrated in order to fit the kinematic hardening model to the desired shape of the backbone curve. For a range of model tests, a close fit is achieved with [START_REF] Abadie | Cyclic Lateral Loading of Monopile Foundations in Cohesionless Soils[END_REF]: Where pU, mh, E and kU are found by direct fitting of the backbone curve. For test CMC0, this gives:

h m U pU k E              
pU=1, mh=3, E=84 and kU= 1. Based on Equation 11, it can be shown that [START_REF] Houlsby | A model for nonlinear hysteretic and ratcheting behaviour[END_REF]:

S S U n N n N n k k  1 , 1 1                               0 2 1 2 1 1 2 1 1 S h h h S h N m m m N n m S pU Uo H n n n H H N k H 
This provides the kinematic hardening values.

Ratcheting behaviour

The ratcheting behaviour is specified by the rate of ratcheting, Rn. First, the amount of ratcheting deformation decreases with cycle number, which means that the ratcheting rate (Rn) decreases with cyclic load history. The implications of the above are twofold: (a) an additional variable that records the cyclic history is needed and (b) Rn will be a decreasing function of this new variable. The accumulated ratcheting strain  is hence defined:

    increments r r d d d    
The choice of function for Rn is based on the experimental results of LeBlanc et al. (2010b) and [START_REF] Abadie | Cyclic Lateral Loading of Monopile Foundations in Cohesionless Soils[END_REF]. They showed that, for the pile geometry and soil conditions explored, and in the case of 1-way cyclic loading, the accumulated deformation can be described by:

          m m m m b p pN p pN N T N T max 0 0 0 0     
with subscript p for the peak cyclic load. Because there is no displacement created by kinematic hardening (Masing rule), the accumulated deformation (pN-p0) is entirely due to ratcheting (pN-p0). m and m are empirical exponents found to be equal to m=4 and m=0.31 and T0 is a dimensionless parameter of 0.5 [START_REF] Abadie | Cyclic Lateral Loading of Monopile Foundations in Cohesionless Soils[END_REF].

According to the above equation, it is expected that the rate of ratcheting will also depend on the stress level. For both cases, a power law function is chosen so that Rn is expressed as:

s r m m n n k R R                           0 0 0 0     
R0 is the initial rate of ratcheting and mr and ms are empirical exponents. 0 is the ultimate stress level, which by definition of , is equal to 1. 0 is the initial value of hardening strain, which is an arbitrarily small value compared with p0, introduced for normalisation purposes. The term

  0  n k
enables to adapt the ratcheting rate to each hardening surface.

Finally, the equation of any unload-reload loop provides an analytical expression for the ratcheting strain at peak load pN:

 

1 max 0 1 0 1 1                h s r r m m U m m pU m m pN k N R           1 1 1 0      h s r h m m m m m      1 1 , 1 ) 1 ( 2 1      

Comparison with experimental results

The calibrated model is used to simulate the tests described in Table 1. Figure 8 shows the predicted response of test HIM, demonstrating that the model captures the Masing behaviour and that ratcheting deformation is negligible during a perfectly symmetric load reversal.

Figure 9(a) displays the stress-strain curve obtained for test CMLT3, computed using the accelerated programme of Table 3. The results exhibit the key features outlined in Section 2 and the graph compares favourably with that of Figure 3. The model is applied to the multi-amplitude cyclic load cases of Table 1. First, tests CMC0-5 and CMLT2 were simulated and the re-load curves are shown in Figure 10. Comparison with Figure 4 shows good agreement with the experimental trends. Likewise, the predicted response of test MALL3 displayed in Figure 11 shows that the resulting stress-strain curve represents the experimental response fairly well (Figure 5). In both cases, the responses feature a competition between (i) Masing behaviour and (ii) ratcheting; this is the key modelling objective of HARM.

Conclusions

This paper outlines a new modelling approach for rigid monopiles subjected to lateral cyclic loading. The methodology builds on experimental test results, identifying that the key mechanisms driving the pile response are a competition between (i) Masing behaviour and (ii) accumulated ratcheting deformation with cycle number.

These results underpin the development of a constitutive model, HARM. A methodology for calibration has also been highlighted and the predicted responses demonstrate close match with experimental results, with an accurate prediction of the key mechanisms that were initially identified. Since the complete load-deflection response of any unload-reload loop is captured, the change in stiffness and damping can also be predicted. This point is not developed further here but is explored in more detail through [START_REF] Abadie | Cyclic Lateral Loading of Monopile Foundations in Cohesionless Soils[END_REF].

This paper provides the basis for a more extensive study, aimed at providing rigorous theoretical models, properly calibrated against experimental data, to be used with confidence by industry for design of the next generation of offshore foundations. On-going and future research explores extension of this modelling concept to Winkler-type approaches and 3D Finite Element Analyses, as well as more complex soil and loading conditions (e.g. multi-directional, partially-drained to undrained sands and clays, layered soils).

Figure 1 :

 1 Figure 1: Lateral loading test equipment

  Figure 3: Test CMLT3: Load-deflection curve

Figure 4 :

 4 Figure 4: Tests CMC0-5, CMLT2: Reload response following cyclic loading at different cycle number

  Typical kinematic hardening and ratcheting model: (a) schematic representation, (b) typical incremental and accelerated response Table 2: Accelerated programme of Figure 7(b)

  to the experimental data CMLT1-3 is obtained with mr=3, ms=12 and R=0.8.

Figure 8 :

 8 Figure 8: Predicted response for test HIM

Finally,

  Figure 11(b) shows the evolution of predicted and measured displacement at peak load, for tests MALL1-3. The results demonstrate that the trends are properly captured, while the quantitative values over-predict the final deformation.

Figure 10 :

 10 Figure 10: Predicted monotonic re-loading response, tests CMC0-5+CMLT2

Table 1 :

 1 Test programme

	Monotonic
	H0	Symmetric reversed loading, 1.5 cycle
	HIM	(b=0.70) Symmetric reversed loading at increasing mag-nitudes, 4.5 cycles (b= .310.440.560.69)
	Continuous cyclic
		b	c	N
	CMLT1	0.31	0.24	100, 000
	CMLT2 * 0.42	0.18	100, 000
	CMLT3	0.47	0.13	100, 000
	Continuous cyclic + Monotonic reload ( * Includes CMLT2)
	CMC0	-	-	0
	CMC1	0.42	0.18	1
	CMC2	0.42	0.18	10
	CMC3	0.42	0.18	100
	CMC4	0.42	0.18	1000
	CMC5	0.42	0.18	10000
	Load packets of various magnitudes
		Load History		Load Cases
	MASL1	1000×A -100×C

-1×E A: (b=0.30; c=0.11)

Acknowledgements

We are grateful to EDF Energies Nouvelles and EDF R&D as well as more recently DONG Energy for their funding of this work, including the financial support of the first Author.


Manuel Lafond

Anne Bergeron

Krister M Swenson

Two Strikes Against the Phage Recombination Problem

The recombination problem is inspired by genome rearrangement events that occur in bacteriophage populations. Its goal is to explain how to transform a bacteriophage population into another using the minimum number of recombinations. Here we show that the general combinatorial problem is NP-Complete, both when the target population contains only one genome of unbounded length, and when the size of the genomes is bounded by a constant. In the first case, the existence of a minimum solution is shown to be equivalent to a 3D-matching problem, and in the second case, to a satisfiability problem. These results imply that the comparison of bacteriophage populations using recombinations will have to rely on heuristics that exploit biological constraints.

Introduction

Genetic recombinations or, more generally, the exchange of DNA material between organisms, have been a source of computational problems since the 1865 report of Gregor Mendel on plant hybridization [START_REF] Abbott | Experiments on Plant Hybrids by Gregor Mendel[END_REF]. Recombinations occur in the reproduction of all living organisms, including asexual reproduction, and are fundamental producers of diversity. In this paper, we study the computational complexity of problems related to modular recombination, which is a form of exchange pervasive in viruses that infect bacteria, called phages.

The biological theory of modular recombination was proposed a few decades ago by Botstein [START_REF] Botstein | A theory of modular evolution for bacteriophages[END_REF], who envisioned " ... viruses as belonging to large interbreeding families, members of which share only a common genome organization consisting of interchangeable genetic elements each of which carries out an essential biological function." The common genome organization that Botstein refers to is the preservation of the order of biological functions, called modules, along the virus genome, although the actual sequences that carry the function may diverge substantially.

The computational models were slower to emerge, since genomic data about "large interbreeding families" were not commonplace until a few years ago. In 2010 a study of a few dozen sequenced strains of Staphyloccoccus aureus was conducted [START_REF] Kahánková | Multilocus pcr typing strategy for differentiation of Staphylococcus aureus siphoviruses reflecting their modular genome structure[END_REF], and a scenario of interbreeding was inferred on the population [START_REF] Swenson | Reconstructing the modular recombination history of staphylococcus aureus phages[END_REF]. The recent availability of other datasets monitoring phage populations evolving through time [START_REF] Kupczok | Rates of Mutation and Recombination in Siphoviridae Phage Genome Evolution over Three Decades[END_REF][START_REF] Lavelle | A decade of streptococcus thermophilus phage evolution in an irish dairy plant[END_REF] or space [START_REF] Chmielewska-Jeznach | Molecular, physiological and phylogenetic traits of lactococcus 936-type phages from distinct dairy environments[END_REF][START_REF] Murphy | Comparative genomics and functional analysis of the 936 group of lactococcal siphoviridae phages[END_REF] suggested the problem of computing the minimum number of recombination events that transforms one population of phages into another. In a previous paper [START_REF] Bergeron | On the comparison of bacteriophage populations[END_REF] we developed a heuristic with approximation bounds based on certain properties of the input and found that, on phages infecting bacteria responsible for cheese fermentation, our heuristic performed well. The question remained, however, as to the computational complexity of the optimization problem.

We answer that question in this article, showing that two basic problems related to the comparison of phage populations are computationally difficult. The first one reduces the problem of finding a perfect 3D-matching to reconstructing a single phage, from a population of phages that represents the triples of the 3D matching instance, with a minimum number of recombinations. The second one reduces a variant of a classic satisfiability problem to the reconstruction of a population of phages, with only 4 modules that represent variables and clauses, with a minimum number of recombinations.

Basic definitions and properties

Phage genomes can adopt either a circular or linear or shape during their life cycle. Genomic data found in databases are linearized by choosing, as a starting point, one module shared by all members of a family, yielding the following representation of phages.

Given an alphabet A, a phage p with n modules can be represented by p = p[0..n -1] where p[a] ∈ A. The recombination operation at positions a and b between two phages p and q :

p = p[0..a -1]|p[a..b -1]|p[b..n -1] q = q[0..a -1]|q[a..b -1]|q[b..n -1]
yields new phages c and d:

c = p[0..a -1]|q[a..b -1]|p[b..n -1] d = q[0..a -1]|p[a..b -1]|q[b..n -1].
Positions a and b are called the breakpoints of the recombination. The recombining phages are called parents, and the newly constructed phages, their children. This relation allows us, when several recombinations are considered, to refer to descendants and ancestors, of both phages and positions; each recombination creates two descendants to the two parents, while the each character in each of the children has exactly one ancestral character from the parents. Note that, naturally, ancestor and descendant relationships are transitive through the generations.

A recombination scenario S from P to Q is a sequence of recombinations that constructs all phages of Q using phages of P and their descendants. Note that no phage is discarded in the process, in the sense that P grows until it is a superset of Q.

The problem that we address in this article is the following: Input: Populations P and Q of equal-length phages, and an integer r. Question: Does there exist a recombination scenario S from P to Q of length at most r?

Minimum Phage Population Reconstruction (MinPPR)

A break in a phage q with respect to the set of phages P is a position b such that for all parents p in P,

p[b -1..b] ̸ = q[b -1..b]. A recombination heals a break b of a phage q if it creates a child c such that c[b -1..b] = q[b -1..b].
In order to be healed, a break b must be one of the breakpoints of the recombination.

Since a recombination can heal at most two breaks in a single phage, if a phage q has n breaks with respect to the set of phages P, then the minimum number of recombinations to construct q is n+1 2 . A crucial remark is that, even if all the breaks are healed, the reconstruction of a phage q with n = 2r breaks with respect to a set of parents might require more than r recombinations. This is the case, for example, if two parents p 1 = 10111 and p 2 = 11101 are used to reconstruct q = 11111: phage q has no break with respect to the set {p 1 , p 2 }, but one recombination is necessary to reconstruct q. This recombination must cut an already healed break in p 1 or p 2 , and we say that the break is reused. Definition 1. In a recombination scenario, a break is said to be reused if it is a breakpoint of more than one recombination in the scenario.

Finally, there is an easy upper bound for the number of recombination necessary to reconstruct a phage: Proposition 1. If there exists a scenario that reconstructs a phage q with n modules from a population P, then there exists one of length at most n -1.

Proof. A scenario exists if, for each position b, there exists a phage p b ∈ P such that p b [b] = q[b], otherwise no recombination can produce the value q[b] at position b. We first recombine p 0 and p 1 using breakpoints 1 and 2, to produce a child that equals q on its first 2 positions, and proceed in a similar way up to position n -1.

⊓ ⊔

Here we study the decision problem where one asks if Q can be generated from P using at most r recombinations, for some given r. Let us first argue that the problem is in NP. A given scenario of r recombinations can be verified in time proportional to r, |P|, and |Q|, but this is not polynomial if r is not polynomial in |P| and |Q| (e.g. if r is exponential). However, Proposition 1 gives an upper bound on the number of required recombinations based on the number of modules. Hence, we may assume that r is bounded by a polynomial in |P| and |Q| and a scenario can be verified in polynomial time, and thus the problem is in NP.

Reconstructing one target genome

We first consider the case in which the population Q consists of a single phage of unbounded length. We reduce the 3D-Perfect-Matching problem to it, where we receive a set of triples T = {(i 1 , j 1 , k 1 ), . . . , (i n , j n , k n )} ⊆ [1..m] 3 , where m ≥ 2 is an integer [START_REF] Karp | Reducibility among combinatorial problems[END_REF]. The goal is to find a subset T ′ ⊆ T of size m such that for any two distinct (i, j, k), (i ′ , j ′ , k ′ ) ∈ T ′ , we have i ̸ = i ′ , j ̸ = j ′ , and

k ̸ = k ′ . Such a set T ′ is called a perfect 3D-matching.
Since there is a single phage in Q, the alphabet is the set {0, 1}, and Q = 11111 . . . 1111 will be the only element of the target population. We consider the following phages, each of length 15m + 2, that form the input population P. See example in Figure 1.

1. For each element (i, j, k) ∈ T , we construct a phage P ijk that has three 1's in positions 5i, 5j + 5m and 5k + 10m, and 0's elsewhere. 2. For each element (i, j, k) ∈ T , we associate three phages, P ij-, P -jk , and P i-k

with two 1's respectively in positions 5i + 1 and 5j -1 + 5m, 5j + 1 + 5m and 5k -1 + 10m, 5i -1 and 5k + 1 + 10m, and 0's elsewhere. 3. P has 0's in every position in which one of the above phages has a 1.

i = 1 i = 2 j = 1 j = 2 k = 1 k = 2 P122 . . . .
Q 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 Fig. 1: Example input with T = {(1, 2, 2), (2, 1, 2), (2, 1, 1), (2, 2, 2)} and m = 2.
The 1's related to phage P 222 are in red. Dots are used to represent the value 0, in order to better highlight the relative positions of the 1's.

We show that T has a 3D-matching if and only if P can generate Q with at most 6m recombinations, implying: Theorem 1. The Minimum Phage Population Reconstruction problem is NP-complete, even when population Q has a single phage.

We have already established that the problem is in NP. For NP-hardness, we show that there exists a 3D-matching T ′ ⊆ T if and only if it is possible to reconstruct Q from P using at most 6m recombinations.

The (⇒) direction

Suppose that there is a 3D-matching T ′ ⊆ T . For each (i, j, k) ∈ T ′ , it is possible to apply three recombinations to P ijk , P ij-, P -jk and P i-k to obtain 00 . . . 01110 . . . 01110 . . . 01110, where the first 111 is centered at column i, the second 111 is centered at column j, the third 111 is centered at column k. The genome P has 000 in these three triples of positions, and so with three more recombinations we can bring in these 111 into P . Use Figure 2 as an illustration. This costs 6 events. Since T ′ is a perfect 3D-matching, we can repeat this m times to fill in all the remaining 000's in P , hence achieving cost 6m.

The (⇐) direction

We next show that if Q can be reconstructed from P using at most 6m recombinations, then T admits a 3D-matching. We first establish several properties that hold in general for scenarios that transform P into Q, before proving the main result of the section.

Given a scenario S that reconstructs phage Q, we identify the following subsets of parents:

1. S ijk contains phages of the form P ijk that belong to the scenario. 2. S xy contains phages of the form P ij-, P -jk or P i-k that belong to the scenario.

Let P = S ijk ∪S xy ∪{ P } be the set of parents that initially belong to scenario S. We prove that scenario S reconstructs phage Q in 6m recombinations only if the set S ijk corresponds to a perfect matching.

By construction, phage Q has 12m breaks with respect to the set of phages P. Since a recombination can heal at most two breaks of a single phage, we need at least 6m recombinations to reconstruct Q. In a scenario of length 6m, no break can be reused.

We distinguish two types of breaks: red breaks connect a phage in S xy to a phage in S ijk , and green breaks connect a phage in S xy to phage P , (see Figure 2). We say that a recombination is red when its two breaks are red, and green if they are green. There is an equal number of red and green breaks in Q, thus, in a scenario of length 6m, the number of red recombinations is equal to the number of green recombinations, and is at most 3m, allowing for eventual red-green recombinations. The following easy result links properties of a recombination scenario to the existence of a perfect matching:

•1• •1• •1• •1• •1• •1• •1• •1• •1• 1• 0 0 0 •1 1• 0 0 0 •1 1• 0 0 0 •1 i j k ... ...
Lemma 1. In any scenario that reconstructs Q, |S ijk | ≥ m. If |S ijk | = m, then the set {(i, j, k)|P ijk ∈ S ijk } is a perfect matching.
Proof. In order to reconstruct Q, all values of i, j and k ∈ [1..m] must appear at least once in the indices of elements

P ijk ∈ S ijk , thus |S ijk | ≥ m. If |S ijk | = m,
then all values of i, j and k ∈ [1..m] appear exactly once implying that the set {(i, j, k)|P ijk ∈ S ijk } corresponds to a perfect matching. See Figure 3 for an example.

⊓ ⊔

In order to show that |S ijk | = m, we first introduce three lemmas that constrain the order of recombinations contained in a scenario of length 6m. The first one concerns the red interval of a phage in S xy , which contains the 0's adjacent to its red breaks, along with the (circularly) intervening columns that are all 0's. See Figure 4 for an example of a red interval.

Lemma 2 (red interval). All descendants of a phage p ∈ S xy must heal red breaks shared with p before acquiring a 1 in p's red interval.

Proof. Consider a descendant of p where one of its red breaks b is not yet healed, along with the first recombination producing a child c that contains b and a 1 in p's red interval. If this recombination does not heal b, then it must first heal a break between b and the 1 in p's red interval, thereby creating a phage c containing both the 1, and the break b. This implies that the second break β of this recombination must be in p's red interval, which is a contradiction since the second break cannot be healed in an interval with all 0's. See Figure 4 for an illustration.

⊓ ⊔

The second lemma establishes the property that all four breaks spanning two consecutive groups of 1's in P will be healed in the same ancestral lineage of Q. ...0 0 0 Lemma 3 (sticky breaks). Consider a break that is healed when phage p is produced by a recombination in a scenario of length 6m. Any adjacent break must be healed in a descendant of p.

i = 1 i = 2 j = 1 j = 2 k = 1 k = 2 P122
Q 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
•1• •1•0 0•0 0 0... b 1• ...•1...
Proof. Consider the 1 in column x that is adjacent to a single healed break in phage p, so that p[x -1..x + 1] = 110, or p[x -1..x + 1] = 011. Consider the 1 that is put in column x + 1, or x -1, while healing the remaining break with column x, producing some child c. Since any break is healed exactly once, the 1's to either side of a healed break must be ancestors of 1's in Q. Therefore, the 1 in column x of both phages p and c must be the ancestor of a 1 in Q, which is only possible if c is a descendant of p. ⊓ ⊔ Lemma 4 (independent S xy ). Consider phages p 1 and p 2 in S xy . There cannot exist a descendant of both phages with an unhealed red break from both p 1 and p 2 .

Proof. Define the green interval of a phage p ∈ S xy to be the interval containing the 0's next to the green breaks in p, along with the (circularly) intervening columns that are all 0's. If the green intervals of p 1 and p 2 do not intersect, the red interval lemma gives the result, since there can be no 1 in either red interval before both of the red breaks in a phage are healed. Suppose their green intervals intersect and, without loss of generality, that the green interval for p 1 starts to the left of the green interval of p 2 . Say that there is a descendant containing the left 1 of p 1 and the left 1 of p 2 . The red interval lemma implies that a recombination happened directly to the left of the 1 in p 2 , which healed that red break. Say that there is a descendant with the left 1 of p 1 and the right 1 of p 2 . Due to the red interval lemma, this implies a recombination happened directly to the right of the 1 in p 2 , which healed that red break. By symmetry, the other cases are covered by those already listed. ⊓ ⊔

The next lemma states a desirable property of recombination scenarios of length 6m, saying that if a phage is in S xy , then its two -unique -siblings are also in S xy . We say that phages p and p ′ eventually recombine in a scenario S if there exists a recombination in S between p, or one of its descendants, and p ′ , or one of its descendants.

Lemma 5. In a scenario S of length 6m, if a phage P ijk ∈ S ijk eventually recombines using the red breaks of a phage in S xy , then all three phages P ij-, P -jk and P i-k eventually recombine with P ijk using both of their red breaks. Proof. Consider the first time that a phage P ijk appears in scenario S, and suppose that this recombination involves a phage in S xy , or its descendant. Without loss of generality we may assume this phage is P i-k , due to the circularity of the genomes and symmetry of our construction. We will show that both P ij- and P -jk must eventually recombine with P ijk in the scenario, using both of their red breaks.

By the lemma statement, both red breaks of P i-k were healed in this recombination producing some child c, having only 0's between columns i and k, with the exception of column j. See Figure 5 for an illustration. Now, the recombination healing the remaining red break adjacent to column i must be in a descendant of c, due to the sticky breaks lemma. Aside from the descendant of c, the other parent p in this recombination must be a descendant of some phage P ij ′ -∈ S xy , for some j ′ .

In the following, we show that the only possible companion breakpoint for this recombination occurs when P ij ′ -= P ij-. Since every recombination must heal two breaks, there is a companion breakpoint between columns i and j, between columns j and k, or (circularly) between columns k and i.

If the companion breakpoint lies (circularly) between k and i, this implies that p has a 1 in the red interval of P ij ′ -, which is impossible by Lemma 2.

Say the companion breakpoint lies between columns i and j. In order to heal two breakpoints, there must be a break b between columns i and j in a descendant p ′ of c. If b is adjacent to j, then it can only be healed using a 1 descending from a phage in S xy . If this phage is not P ij-, then the independent S xy lemma prohibits a common descendant between this phage and P ij ′ -, a contradiction. Say b is not adjacent to j, but rather in the zone of all 0's in c. Then the existence of b implies that there has been a recombination at the breakpoint adjacent to column j in an ancestor of p ′ . This leads to the same contradiction as in the previous case.

The same argument applies to a breakpoint occurring between j and k. Now consider the symmetric case, where a recombination heals the remaining break adjacent to column k in phage c. The same reasoning shows that both red breaks of P -jk are used to recombine with a descendant of P ijk .

⊓ ⊔ The previous lemma depends on the assumption that the first recombination with a phage in S ijk heals the two red breaks of a phage in S xy . The following lemma show that this must be the case. Lemma 6. In a scenario S of length 6m, the first recombination using P ijk ∈ S ijk must heal both red breaks of a phage in S xy , be it P ij-, P -jk , or P i-k .

1•1•0 0 0 aefewf 0 0 0•1•0 0 0 0 0 0•1•1 ... •1 ... range for j indices •1 1• •1 1• 1• j k i ... ... ... c p ... ...
Proof. Since a phage P ijk ∈ S ijk has only red breaks, it must eventually recombine using exactly two red breaks, b 1 and b 2 . Suppose that P ijk does not eventually recombine with a single phage of S xy , then b 1 and b 2 are breaks on different phages p 1 and p 2 of S xy . This implies that p 1 and p 2 eventually recombine to produce a child containing both b 1 and b 2 . Such a recombination is impossible, due to Lemma 2.

⊓ ⊔

Thus we have the result:

Proposition 2.
A scenario S reconstructs Q in 6m recombinations only if the set S ijk is a perfect matching.

Proof. Lemma 5 shows that for each P ijk in S ijk , the three corresponding P ij-, P -jk and P i-k must belong to the scenario. Since there are 3m pairs of red breaks, the maximum number of elements of S ijk is m. Lemma 1 gives the result. ⊓ ⊔

NP-hardness for genomes of length 4

In the preceding section, we showed that the Minimum Phage Population Reconstruction was hard when the length of the genomes was unbounded. Is it still the case for genome of bounded length? The answer is yes, and we dedicate the remainder of the section to the proof of the following statement:

Theorem 2. The Minimum Phage Population Reconstruction problem is NP-complete, even when the genomes of P and Q have length 4.

We reduce from the Balanced-4Occ-SAT problem, where we are given a boolean formula ϕ in conjunctive normal form, such that each variable has exactly two positive occurrences in the clauses of ϕ, and exactly two negative occurrences [START_REF] Berman | Approximation hardness of short symmetric instances of max-3sat[END_REF].

Consider an instance ϕ of Balanced-4Occ-SAT with variables x 1 , . . . , x n and clauses C 1 , . . . , C m . We construct a corresponding instance (P, Q, r) of the phage problem. See Figure 6 for an example with 3 variables and 4 clauses, and Figure 7 for a more abstract view.

The alphabet for the phages of P and Q has, for each variable x i , a corresponding symbol i ∈ [1..n], and for each clause C j , a corresponding symbol c j . We also add two unique symbols '-' and '•' to the alphabet.

Consider a variable x i , where i ∈ [1.

.n]. Let C g , C h be the clauses in which x i occurs positively, and C r , C s those in which x i occurs negatively. Add the following phages to P:

X i = i • i i X + i = c g i c h - X - i = c r i c s -

The (⇒) direction

Suppose that ϕ is satisfied by an assignment α : {x 1 , . . . , x n } → {true, f alse} of the variables. Let us produce Q from P. For each i ∈ [1.

.n], if α(x i ) = true, then recombine X i with X + i by exchanging 2nd positions:

Xi = i • i i X + i = cg i c h - X * i = i i i i p = cg • c h - This produces children X * i and p = c g •c h -,
where C g and C h are the clauses that are satisfied by setting x i to true. At this point, if D * g is not already in Q, then recombine D g,1 = -c g c g with p = c g •c hby exchanging positions 1 and 2, thereby producing D * g :

Dg,1 = --cg cg p = cg • c h - D * g = cg • cg cg --c h -
For an illustration of the previous two recombinations, see the black edges in Figure 7. In the same way, if

D * h is not already in Q, recombine D h,2 = c h --c h with c j •c h -by
exchanging positions 2 and 3, which produces D * h . If instead α(x i ) = f alse, recombine X i with X - i by exchanging 2nd positions. This creates X * i and c r • c s -, where C r and C s are satisfied by setting x i to f alse. Produce D * r and D * s , if not already there, as in the previous case. Since every clause C j is satisfied by α, for each D * j , there will be some X i in the above procedure that produces D * j . Moreover, there are exactly n + m recombinations: one to produce each X * i , and one to produce each D * j .

The (⇐) direction

Suppose that there exists a sequence S = (R 1 , . . . , R r ) of at most r ≤ n + m recombinations that reconstructs Q. Let X = {X 1 , . . . , X n } and D = {D 1 , . . . , D m } where D j is either D j,1 or D j,2 , whichever contains the character c j that is the ancestor of the c j in the 4th position of D * j . We define a function f : X ∪ D → S in the following way: set f (X i ) to the first recombination that creates X * i , or an ancestor of X * i , with [1..n] in the 2nd position and i in the 4th. Note that the parent used by f (X i ), having i in the 4th position, must also have {-, •} in the 2nd position, as otherwise there would be a previous recombination with the required properties for being f (X i ).

Set f (D j ) to the first recombination that creates D * j , or an ancestor of D * j , with [1..n] ∪ {•} in the 2nd position and c j in the 4th. As previously, note that one of the parents used by f (D j ) contains '-' in the 2nd position and c j in the 4th position. × symbol. In black, X i and X + i first recombine to produce q = c j • c h -, which then recombines with D j,1 to produce D * j (and some other unused phage). Notice that q can also be used to produce D * h in a similar manner. In dark gray, an alternate way to produce X * i and D * j is depicted. Here, X + i is first recombined with D j,1 to produce p = c j ic j c j . Phage p is not in Q, but can be recombined with X i to produce both X * i and D * j in one operation. Proposition 3. The function f is a bijection.

Proof. We prove that f is injective, and since there can be at most n + m recombinations in S, we conclude that f is a bijection.

If i ̸ = k, then f (X i ) ̸ = f (X k ) since equality implies a recombination where at least one parent has an element of {-, •} in 2nd position, and both children have an element of [1..n].

If g ̸ = h, then f (D g ) ̸ = f (D h ) since equality implies a recombination where at least one parent has a '-' in 2nd position, and both children do not.

Finally, f (X i ) ̸ = f (D j ) since equality implies a recombination where one parent has a '-' in 2nd position, and both children do not.

⊓ ⊔

The crucial consequences of Proposition 3 are the three following results:

Proposition 4. There is exactly one element of X ∪ D in each recombination of S.

Proof. Since all phages of X ∪ D are necessary to produce Q, then each one must be used by at least one recombination. We show that the image of f contains no recombination between two of these phages.

Suppose X i and X k recombine. Then this recombination is not f (X i ) or f (X k ) since both parents have a '•' in 2nd position, implying that none of the children have an element of [1..n].

Suppose D g and D h recombine. Then this recombination is not f (D g ) or f (D h ) because both parents have a '-' in 2nd position, implying that both children have a '-' in 2nd position.

Suppose X i and D j are used in a recombination R. Then R is not f (X i ) because one parent has a '•' in 2nd position, and the other has a '-', implying that none of the children have an element of [1..n]. So R must be f (D j ), implying that one of the children p is an ancestor of D * j , having • in 2nd position and c j in 4th position. In the remaining paragraph we show that there is no recombination from S that makes the c j in 4th position of p the ancestor of the c j in 4th position of D * j , contradicting the existence of R.

Consider any subsequent recombination R ′ that uses p as a parent. Each R ′ cannot be a f (X i ) since p cannot contribute an element of [1..n] in 2nd position, so it must be f (D h ), for D h ∈ D and h ̸ = j. Note that the other parent q, used in R ′ , must have '-' in 2nd position and c h in 4th position, so the children are then q ′ with • and c h in 2nd and 4th positions, and p ′ with '-' and c j in 2nd and 4th positions. But, while p ′ is the child that contains the c j that could be ancestral to the 4th position of D * j , it cannot be used as a parent in a recombination of f (X ∪ D), since it has '-' in 2nd position, and f (D j ) has already been applied. Therefore, p ′ is not an ancestor of D * j . This is true for any such child p ′ produced by a subsequent recombination, which contradicts that p is an ancestor of D * j . This contradict our initial supposition that R is f (D j ).

⊓ ⊔ Corollary 1. The recombination that uses X i produces X * i . Proof. Note that X i is the only phage in X ∪ D that shares any character with X * i . Since both parents of any recombination creating X * i must share at least one character with X * i , the recombination in S that uses X i is the only one that can create X * i . ⊓ ⊔ Corollary 2. The recombination that uses D j produces either D * j or c j ic j c j . Proof. By definition D j is an ancestor of D * j , and by Corollary 4 we know that any descendant of D j must recombine with an element of X ∪ D. Since no other member of X ∪ D has c j in positions 1, 3 or 4, a child of D j must be an ancestor of D * j and have that character c j in those positions. This child must be either D * j , p = c j ic j c j , or q = c jc j c j . A subsequent recombination using child q = c jc j c j does not exist since it would either recombine with an X k , but not produce X * k in contradiction with Corollary 1, or with a D h , whose 2nd position is also '-'. Therefore, q has no descendant, contradicting that it is an ancestor of D * j .

⊓ ⊔

We now establish that a recombination scenario S of length n + m implies a valid, and satisfiable truth assignment for ϕ.

For an X i ∈ X , we say that X i chose X + i if the only recombination that uses X i is with an X + i or its descendant, and we say that it chose X - i if X i recombined with X - i or its descendant. If X i chose X + i , we set x i = true, and if X i chose X - i we set x i = f alse. Let us call the resulting assignment α, which we claim is satisfying. We first argue that each x i is assigned only one value, and then show that each clause is satisfied. Proposition 5. A recombination scenario S of length n + m implies a valid, and satisfiable truth assignment for ϕ.

Proof. We first show that the assignment α is well-defined, i.e. each variable x i is assigned either true or f alse, but not both. Since each X i chooses at least one of X + i or X - i , we know that x i is assigned true or f alse. Assume now that x i is assigned both. Then X i chose both X + i and X - i , meaning that X i recombines with a phage that descends from both X + i and X - i . But the existence of this phage requires a recombination between descendants of both X + i and X - i , which contradicts Corollary 4.

We now show that each clause is satisfied by α. By Corollary 2, we know that D j,1 and D j,2 do not recombine, and only one of the two, that we called D j , appears in S. Since D j contributes at most two c j characters to D * j , the third c j can only be a descendant of an X + i or X - i , since they are the only other phages in P that may contain c j . By construction, this means that the clause C j is satisfied by the variable x i being set to the truth assignment implied by the corresponding X + i or X - i . ⊓ ⊔

Conclusion

The notion of recombination used in this article is the same as the two-point crossover function [START_REF] Gitchoff | Recombination induced hypergraphs: a new approach to mutation-recombination isomorphism[END_REF] used for characterizing fitness landscapes for the exploration of genotypes. This two-point crossover has since been studied in a general form as a k-point crossover [START_REF] Changat | Transit sets of-point crossover operators[END_REF]. While, to the best of our knowledge, there is no work directly linking the area of fitness landscape exploration to the minimization problem discussed in this article, we hope that these related areas can be fused in the future.

In this paper, we showed that the Minimum Phage Population Reconstruction is NP-Complete in two extreme cases: bounded length, and a single target phage. Although negative, such results may come as a relief, since we can turn our focus to algorithms that, by accounting for biological constraints, could provide drastically reduced search spaces for parsimonious solutions. For example, the use of other measures of evolution, or information about community structure [START_REF] Chevallereau | Interactions between bacterial and phage communities in natural environments[END_REF] might play a significant role in reducing the complexity of the problem: after all, phages recombine all the time, and thrive doing so.

  1 . . . . . . . . . . . . . . 1 . . . . . . . . . 1 . . P212 . . . . . . . . . 1 . . . . 1 . . . . . . . . . . . . . . 1 . . P211 . . . . . . . . . 1 . . . . 1 . . . . . . . . . 1 . . . . . . . P222 . . . . . . . . . 1 . . . . . . . . . 1 . . . . . . . . . 1 . . P22-. . . . . . . . . . 1 . . . . . . . 1 . . . . . . . . . . . . . P12-. . . . . 1 . . . . . . . . . . . . 1 . . . . . . . . . . . . . P21-. . . . . . . . . . 1 . . 1 . . . . . . . . . . . . . . . . . . P-22 . . . . . . . . . . . . . . . . . . . . 1 . . . . . . . 1 . . . P-12 . . . . . . . . . . . . . . . 1 . . . . . . . . . . . . 1 . . . P-11 . . . . . . . . . . . . . . . 1 . . . . . . . 1 . . . . . . . . P2-2 . . . . . . . . 1 . . . . . . . . . . . . . . . . . . . . . 1 . P1-2 . . . 1 . . . . . . . . . . . . . . . . . . . . . . . . . . 1 . P2-1 . . . . . . . . 1 . . . . . . . . . . . . . . . . 1 . . . . . . P 1 1 1 . . . 1 1 . . . 1 1 . . . 1 1 . . . 1 1 . . . 1 1 . . . 1

Fig. 2 :

 2 Fig.2: Red and green breaks. Red breaks connect a phage in S xy to a phage in S ijk , and green breaks connect a phage in S xy to phage P . Phage Q has 6m red breaks and 6m green breaks.

  . . . . 1 . . . . . . . . . . . . . . 1 . . . . . . . . . 1 . . P211 . . . . . . . . . 1 . . . . 1 . . . . . . . . . 1 . . . . . . . P12-. . . . . 1 . . . . . . . . . . . . 1 . . . . . . . . . . . . . P21-. . . . . . . . . . 1 . . 1 . . . . . . . . . . . . . . . . . . P-11 . . . . . . . . . . . . . . . 1 . . . . . . . 1 . . . . . . . . P-22 . . . . . . . . . . . . . . . . . . . . 1 . . . . . . . 1 . . . P1-2 . . . 1 . . . . . . . . . . . . . . . . . . . . . . . . . . 1 . P2-1 . . . . . . . . 1 . . . . . . . . . . . . . . . . 1 . . . . . P 1 1 1 . . . 1 1 . . . 1 1 . . . 1 1 . . . 1 1 . . . 1 1 . . . 1

Fig. 3 :

 3 Fig. 3: A possible output for the example of Figure 1, with the sets S ijk and S xy used by a recombination scenario of length 6m = 12. Here |S ijk | = m, and {(1, 2, 2), (2, 1, 1)} is a perfect matching.

Fig. 4 :

 4 Fig.4: A recombination between a descendant d of p ∈ S xy and a phage p ′ . The recombination creates phage c, in blue, which contains both break b and a 1 from a phage p ′ , where the 1 is located in a column from p's red interval. This recombination must have one breakpoint between b and the other red break of p, and one breakpoint β after the location of the 1 in phage p ′ . Breakpoint β cannot heal a break because it is between two 0's.

Fig. 5 :

 5 Fig.5: At the creation of child c, it has only 0's between columns i and k, with the exception of column j. The gray 1's in p are possible breaks that can be healed.

Fig. 7 :

 7 Fig. 7: An illustration of how phages X i , X + i , D j,1 , and D h,2 can recombine to produce X * i , D * j , and D * h . Recombinations are represented by pairs of arrows that meet in the middle at a box with a ■

Acknowledgments

ML was partially supported by Canada NSERC Grant number 2019-05817. AB was partially supported by Canada NSERC Grant number 05729-2014. KMS was partially supported by the grant ANR-20-CE48-0001 from the French National Research Agency (ANR).

and add the following to Q:

Now for each clause C j , j ∈ [1.

.m], add the following to P:

and add the following to Q:

We show that ϕ is satisfiable if and only if Q can be reconstructed from P with at most r = n + m recombinations.

Coding the clauses:

Other phages in P:

Fig. 6: In this example there are n = 3 variables and m = 4 clauses, thus n+m = 7 phages in Q. One possible recombination scenario of length 7 is depicted. The three recombinations in Group 1 first construct * 1 using X + that asserts that clauses 1 and 3 are satisfied when variable x 1 is true. The resulting phage q is then used to generate both D * 1 and D * 3 . In 2, X * 2 and D * 2 are constructed using X - 2 that asserts that clause 2 is satisfied when when variable x 2 is false. Group 3 shows an alternative strategy that first constructs phage p = c 4 3 c 4 c 4 , and uses it to simultaneously construct X * 3 and D * 4 .


Quentin Febvre
email: quentin.febvre@imt-atlantique.fr

Julien Le Sommer

Clément Ubelmann

Ronan Fablet

J Le Sommer

Training neural mapping schemes for satellite altimetry with simulation data

la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Plain Language Summary

For an artificial intelligence (AI) to learn, one need to describe a task using data and an evaluation procedure. Here we aim at constructing images related to the ocean surface currents. The satellite data we use provide images of the ocean surface with a lot of missing data (around 95% of missing pixels for a given day), and we aim at finding the values of the missing pixels. Because we don't know the full image, it is challenging to train an AI on this task using only the satellite data. However, today's physical knowledge makes it possible to numerically simulate oceans on big computers. For these simulated oceans, we have access to the gap-free image, so we can train AI models by first hiding some pixels and checking if the model fill the gaps with the correct values. Here, we explore under which conditions AIs trained on simulated oceans are useful for the real ocean. We show that today's simulated oceans work well for training an AI on this task and that training on more realistic simulated oceans improve the performance of the AI!

Introduction

Satellite altimeters have brought a great leap forward in the observation of sea surface height on a global scale since the 80's.

Altimetry data have greatly contributed to the monitoring and understanding of key processes such as the sea-level rise and the role of mesoscale dynamics. The scarce and irregular sampling of the measurements presents a challenge for training deep neural networks. The retrieval of mesoscale-to-submesoscale sea surface dynamics for horizontal scales smaller than 150 km however remains a challenge for operational systems based on optimal interpolation [START_REF] Taburet | DUACS DT2018: 25 years of reprocessed sea level altimetry products[END_REF] and data assimilation [START_REF] Lellouche | The Copernicus Global 1/12°Oceanic and Sea Ice GLORYS12 Reanalysis[END_REF] schemes. This has motivated a wealth of research to develop novel mapping schemes [START_REF] Ballarotta | Dynamic Mapping of Along-Track Ocean Altimetry: Performance from Real Observations[END_REF][START_REF] Ubelmann | Reconstructing Ocean Surface Current Combining Altimetry and Future Spaceborne Doppler Data[END_REF][START_REF] Guillou | Mapping Altimetry in the Forthcoming SWOT Era by Back-and-Forth Nudging a One-Layer Quasigeostrophic Model[END_REF].

In this context, data-driven and learning-based approaches (Alvera [START_REF] Azcarate | Reconstruction of incomplete oceanographic data sets using empirical orthogonal functions: Application to the Adriatic Sea surface temperature[END_REF][START_REF] Barth | DINCAE 2.0: Multivariate convolutional neural network with error estimates to reconstruct sea surface temperature satellite and altimetry observations[END_REF][START_REF] Lguensat | The Analog Data Assimilation[END_REF][START_REF] Fablet | END-TO-END PHYSICS-INFORMED REPRESENTATION LEARNING FOR SATELLITE OCEAN REMOTE SENSING DATA: APPLICATIONS TO SATELLITE ALTIMETRY AND SEA SURFACE CURRENTS[END_REF][START_REF] Martin | Synthesizing Sea Surface Temperature and Satellite Altimetry Observations Using Deep Learning Improves the Accuracy and Resolution of Gridded Sea Surface Height Anomalies[END_REF] appear as appealing alternatives to make the most of the available observation and simulation datasets. Especially, Observing System Simulation Experiments (OSSE) have stressed the potential of neural schemes trained through supervised learning for the mapping of satellite-derived altimetry data [START_REF] Fablet | END-TO-END PHYSICS-INFORMED REPRESENTATION LEARNING FOR SATELLITE OCEAN REMOTE SENSING DATA: APPLICATIONS TO SATELLITE ALTIMETRY AND SEA SURFACE CURRENTS[END_REF][START_REF] Beauchamp | [END_REF]. Their applicability to real datasets has yet to be assessed and recent studies have rather explored learning strategies from real gappy multi-year altimetry datasets [START_REF] Martin | Synthesizing Sea Surface Temperature and Satellite Altimetry Observations Using Deep Learning Improves the Accuracy and Resolution of Gridded Sea Surface Height Anomalies[END_REF]. Despite promising results, schemes trained with unsupervised strategies do not reach the relative improvement of the operational processing suggested by OSSEbased studies.

Here, we go beyond using OSSEs as benchmarking-only testbeds. We explore their use for the training of neural mapping schemes and address the space-time interpolation of real satellite altimetry observations. Through numerical experiments on a Gulf Stream case-study with a 5-nadir altimeter constellation, our main contributions are three-fold. We demonstrate the relevance of the simulation-based learning of neural mapping schemes and their generalization performance for real nadir altimetry data. We benchmark the proposed approach with state-of-the-art operational products as well as neural schemes trained from real altimetry datasets. We also assess how the characteristics of the training datasets, especially in terms of resolved ocean processes, drives the mapping performance. To ensure the reproducibility of our results, our code is made available through an open source license along with the considered datasets and the trained models [START_REF] Febvre | Code and Data release Simulation-based 4dVarNets for satellite altimetry[END_REF].

The content of this paper is organized as follows. Section 2 offers background information on related work, Section 3 presents our method, Section 4 reports our numerical experiments, and Section 5 elaborates on our main contributions.

Background

Gridded satellite altimetry products

The ability to produce gridded maps from scattered along-track nadir altimeter measurements of sea surface height is key to the exploitation of altimeter data in operational services and science studies [START_REF] Abdalla | Altimetry for the future: Building on 25 years of progress[END_REF]. As detailed below, we can distinguish three categories of approaches to produce such maps: reanalysis products (Lellouche et al., 2021) using data assimilation schemes, observation-based products [START_REF] Taburet | DUACS DT2018: 25 years of reprocessed sea level altimetry products[END_REF] and learning-based approaches [START_REF] Fablet | END-TO-END PHYSICS-INFORMED REPRESENTATION LEARNING FOR SATELLITE OCEAN REMOTE SENSING DATA: APPLICATIONS TO SATELLITE ALTIMETRY AND SEA SURFACE CURRENTS[END_REF].

Reanalysis products such as the GLORYS12 reanalysis [START_REF] Lellouche | The Copernicus Global 1/12°Oceanic and Sea Ice GLORYS12 Reanalysis[END_REF] leverage the full expressiveness of state-of-the-art ocean models. They aim at retrieving ocean state trajectories close to observed quantities through data assimilation methods including among others Kalman filters and variational schemes [START_REF] Carrassi | Data assimilation in the geosciences: An overview of methods, issues, and perspectives[END_REF]. Such reanalyses usually exploit satellite-derived and in situ data sources. For instance, GLO-RYS12 reanalysis assimilates satellite altimetry data, but also satellite-derived observations of the sea surface temperature, sea-ice concentration as well as in situ ARGO data [START_REF] Wong | Argo Data 1999-2019: Two Million Temperature-Salinity Profiles and Subsurface Velocity Observations From a Global Array of Profiling Floats[END_REF].

The second category involves observation-based products. In contrast to reanalyses, they only rely on altimetry data and address a space-time interpolation problem. They usually rely on simplifying assumptions on sea surface dynamics. In this category, optimal-interpolation-based product DUACS (Data Unification and Altimeter Combination System) [START_REF] Taburet | DUACS DT2018: 25 years of reprocessed sea level altimetry products[END_REF]) exploits a covariance-based prior, while recent studies involve quasi-geostrophic dynamics to guide the interpolation scheme [START_REF] Guillou | Mapping Altimetry in the Forthcoming SWOT Era by Back-and-Forth Nudging a One-Layer Quasigeostrophic Model[END_REF][START_REF] Ballarotta | Dynamic Mapping of Along-Track Ocean Altimetry: Performance from Real Observations[END_REF].

Data-driven and learning-based approaches form a third category of SSH mapping schemes. Similarly to observation-based methods, they are framed as interpolation schemes.

Especially deep learning schemes have gained some attention. Recent studies have explored different neural architectures both for real and OSSE altimetry datasets [START_REF] Archambault | Multimodal Unsupervised Spatio-Temporal Interpolation of satellite ocean altimetry maps[END_REF][START_REF] Beauchamp | Data-driven and learning-based interpolations of along-track Nadir and wide-swath SWOT altimetry observations[END_REF][START_REF] Martin | Synthesizing Sea Surface Temperature and Satellite Altimetry Observations Using Deep Learning Improves the Accuracy and Resolution of Gridded Sea Surface Height Anomalies[END_REF]. These studies investigate both different training strategies as well as different neural architectures from off-the-shelf computer vision ones such as convolutional LSTMs and UNets [START_REF] Ronneberger | U-Net: Convolutional Networks for Biomedical Image Segmentation[END_REF] to data-assimilation-inspired ones [START_REF] Beauchamp | Data-driven and learning-based interpolations of along-track Nadir and wide-swath SWOT altimetry observations[END_REF]Fablet, Chapron, et al., 2021).

Ocean Modeling and OSSE

Advances in modeling and simulating ocean physics have largely contributed to a better understanding of the processes involved in the earth system and to the development of operational oceanography [START_REF] Barnier | Impact of partial steps and momentum advection schemes in a global ocean circulation model at eddy-permitting resolution[END_REF][START_REF] Ajayi | Spatial and Temporal Variability of the North Atlantic Eddy Field From Two Kilometric-Resolution Ocean Models[END_REF]. High-resolution simulations used in Observing System Simulation Experiments (OSSE) also provide a great test-bed for the design and evaluation of new of ocean observation systems [START_REF] Benkiran | Assessing the Impact of the Assimilation of SWOT Observations in a Global High-Resolution Analysis and Forecasting System Part 1: Methods[END_REF]. The availability of numerical model outputs enables the computation of interpretable metrics directly on the quantities of interest. This avoids challenges met when working solely with observation data that may be incomplete, noisy or indirectly related to the desired quantity. For example, in the case of the recently launched SWOT mission, OSSEs combined ocean and instrument simulations to address calibration issues and interpolation performance for SWOT altimetry data [START_REF] Dibarboure | Data-Driven Calibration Algorithm and Pre-Launch Performance Simulations for the SWOT Mission[END_REF]. Such OSSEs have also promoted novel developments for the interpolation of satellite altimetry such as the BFN-QG and 4DVarNet schemes [START_REF] Guillou | Mapping Altimetry in the Forthcoming SWOT Era by Back-and-Forth Nudging a One-Layer Quasigeostrophic Model[END_REF][START_REF] Beauchamp | [END_REF].

In OSSE settings, we can train learning-based mapping schemes in a supervised manner using model outputs as the "ground truth" during the training phase. Nonetheless, these training methods cannot be straightforwardly applied to Observing System Experiments (OSEs) due to a lack of comprehensive groundtruthed observation datasets. Applied machine learning practitioners often grapple with insufficient amount of labelled data during the training of supervised learning schemes, as the collection of large annotated datasets for a specific task can be costly or unattainable. Proposed solutions includes the exploitation of large existing datasets (such as ImageNet [START_REF] Deng | Ima-geNet: A large-scale hierarchical image database[END_REF]) to train general purpose models (like [START_REF] He | Deep Residual Learning for Image Recognition[END_REF]). Another approach involves the generation of synthetic datasets to facilitate the creation of groundtruthed samples (Gomez [START_REF] Gomez Gonzalez | VIP: Vortex Image Processing Package for High-contrast Direct Imaging[END_REF][START_REF] Dosovitskiy | FlowNet: Learning Optical Flow with Convolutional Networks[END_REF]. OSSEs, which combine ocean model outputs and observing system simulators [START_REF] Boukabara | Community Global Observing System Simulation Experiment (OSSE) Package (CGOP): Assessment and Validation of the OSSE System Using an OSSE-OSE Intercomparison of Summary Assessment Metrics[END_REF], can deliver such large synthetic groundtruthed datasets. We propose to investigate how OSSE-based training strategies apply to the analysis of real satellite altimetry datasets. Recent results of SSH superresolution model trained on simulation datasets and evaluated on real ones [START_REF] Buongiorno Nardelli | Super-Resolving Ocean Dynamics from Space with Computer Vision Algorithms[END_REF] support the relevance of such strategies.

Physics-aware deep-learning

In the last decades, DL advances combined with the rise in computational resources and amount of data have shown the power of extracting knowledge from data in domains ranging from computer vision to language processing [START_REF] Lecun | Deep learning[END_REF]. Yet, despite to the universality of DL architectures [START_REF] Hornik | Multilayer feedforward networks are universal approximators[END_REF], a central challenge persists in learning from data: the generalization performance beyond the distribution of the training data. To tackle this problem, the literature includes a variety of strategies such as data augmentation [START_REF] Shorten | A survey on Image Data Augmentation for Deep Learning[END_REF] and regularization techniques, including dropout layers [START_REF] Srivastava | Dropout: A Simple Way to Prevent Neural Networks from Overfit-ting[END_REF] and weight decay schemes [START_REF] Krizhevsky | ImageNet Classification with Deep Convolutional Neural Networks[END_REF]. This is of critical importance for physical systems, where models trained on past data will be challenged when the system evolves and reaches dynamics absent from the training data. We can see evidence of this shortcoming in the instability challenges faced by neural closures for climate models [START_REF] Brenowitz | Interpreting and Stabilizing Machine-Learning Parametrizations of Convection[END_REF].

There have been a variety of approaches to harness physical priors within learning schemes to address this issue. Some injects trainable components in classical integration schemes of physical models such as [START_REF] Yin | Augmenting physical models with deep networks for complex dynamics forecasting*[END_REF]. Others leverage physical priors within their learning setups which can been used in the training objective [START_REF] Raissi | Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations[END_REF][START_REF] Greydanus | Hamiltonian Neural Networks[END_REF], as well as in the architecture [START_REF] Li | Fourier neural operator for parametric partial differential equations[END_REF][START_REF] Wang | Towards physics-informed deep learning for turbulent flow prediction[END_REF]. However most of these works have focused on relatively simple physical models and it remains challenging to combine current state-of-the-art ocean models with such methods. Obstacles include the complexity and cost of running the physical models, the differences in programming tools and the computing infrastructures used in each domain, as well as the availability of automatic differentiation tools for state-of-the-art ocean models.

The proposed simulation-based training strategy offers another way to benefit from the advances in high-resolution ocean modeling in the design of deep neural models for ocean reanalysis problems.

Method

Overview

We designate our approach as "simulation-based", it consists in leveraging ocean models and simulations of observing systems to design supervised training environments. In this section, we describe the proposed method for assessing the potential of simulationbased neural schemes for the mapping real altimetry tracks. We describe the architecture considered in our study, as well as the different datasets used for training purposes. We also detail our simulation-based training setup and the proposed evaluation framework on real altimetry.

Figure 1: Overview of the experimental setup. On the left side we display the simulation-based training strategy based on an ocean simulation which will be used for 1) generating synthetic observation and 2) computing the training objective of the neural mapping scheme. On the right side we show the evaluation principle of splitting the available satellite observations to evaluate the method on data that were not used for the inference.

Neural mapping scheme

The neural mapping scheme considered for this study is the 4DVarNet framework [START_REF] Fablet | END-TO-END PHYSICS-INFORMED REPRESENTATION LEARNING FOR SATELLITE OCEAN REMOTE SENSING DATA: APPLICATIONS TO SATELLITE ALTIMETRY AND SEA SURFACE CURRENTS[END_REF]. We choose this scheme due to the performance shown in the OSSE setup. As reported in [START_REF] Beauchamp | [END_REF], it significantly outperforms the DUACS product [START_REF] Taburet | DUACS DT2018: 25 years of reprocessed sea level altimetry products[END_REF] in the targeted Gulf stream region. 4DVarNet relies on a variational data assimilation formulation. The reconstruction results from the minimization of a variational cost. This cost encapsulates a data fidelity term and a regularization term. It exploits a prior on the space-time dynamics through a convolutional neural network inspired from [START_REF] Fablet | Bilinear Residual Neural Network for the Identification and Forecasting of Geophysical Dynamics[END_REF], and an iterative gradient-based minimization based on a recurrent neural network as introduced by [START_REF] Andrychowicz | Learning to learn by gradient descent by gradient descent[END_REF]. The overall architecture and components are similar to those presented in [START_REF] Beauchamp | [END_REF]. We adapt some implementation details based on cross-validation experiments to improve the performance and reduce the training time. We refer the reader to the code for more details [START_REF] Febvre | Code and Data release Simulation-based 4dVarNets for satellite altimetry[END_REF].

SSH Data

Resolution Reanalysis Tide DAC NATL60 [START_REF] Ajayi | Spatial and Temporal Variability of the North Atlantic Eddy Field From Two Kilometric-Resolution Ocean Models[END_REF] 1/60 The last column indicate whether the Dynamic Atmospheric Correction was applied on the synthetic SSH. It justify the presence of both eNATL60-0 and NATL60 to isolate the impacts of resolution and tide.

We use numerical simulations of ocean general circulation models (OGCM) to build our reference SSH datasets. Such simulations involve a multitude of decisions that affect the resulting simulated SSH. Here we consider NEMO (Nucleus for European Modelling of the Ocean) (Gurvan et al., 2022) which is among the state-of-the art OGCM in operational oceanography [START_REF] Ajayi | Spatial and Temporal Variability of the North Atlantic Eddy Field From Two Kilometric-Resolution Ocean Models[END_REF] as well as in climate studies [START_REF] Voldoire | The CNRM-CM5.1 global climate model: Description and basic evaluation[END_REF]. The selected SSH datasets reported in Table 1 focus on three main aspects: the added-value of high-resolution eddy-rich simulations, the impact of reanalysis datasets and the relevance of tide-resolving simulations.

In order to evaluate the impact of eddy-rich simulations, we consider NATL60, GLORYS12f and ORCA025 free runs, respectively with a horizontal grid resolution of 1/60 • , 1/12 • , and 1/4 • . Finer grids allow for more processes to be simulated. We therefore expect higherresolution simulations to exhibit structures closer to the real ocean and the associated trained deep learning model to perform better. Regarding the impact of reanalysis data, we compare numerical experiments with the GLORYS12-r reanalysis and the associated free run GLORYS12-f. This reanalysis dataset relies on the assimilation of temperature, sea level and sea ice concentration observations. Besides, the recent eNATL60 twin simulations eNATL60-t and eNATL60-0 allow us to evaluate the impact of tide-resolving simulations. We summarize in Table 1 the characteristics of the different datasets.

OSSE-based training setup

We sketch the proposed OSSE-based training setup on the left side of the Figure 1. In order to fairly evaluate the datasets' quality as a training resource, we standardize the training procedure. We regrid all simulations to the same resolution (1/20°) and we use daily-averaged SSH fields as training targets. We generate noise-free pseudo-observations by sampling values of the daily-averaged fields corresponding to realistic orbits of a 5 altimeterconstellation. We train all models from a one-year dataset in a Gulfstream domain from (66°W, 32°N) to (54°W, 44°N) in which we keep the same two months for validation. The hyper-parameters of the model and training procedure such as the number of epoch, learning rate scheduler are the same for all the experiments. The detailed configuration can be found by the reader in the available implementation. As training objective, we combine the mean square errors for the SSH fields and the amplitude of the gradients as well as a regularization loss for the prior model.

OSE-based evaluation setup

As sketched on the right side of the Figure 1, the evaluation setup relies on real altimetry data from the constellation of 6 satellites from 2017 (SARAL/Altika, Jason 2, Jason 3, Sentinel 3A, Haiyang-2A and Cryosat-2 ). We apply the standardized setup presented in a data-challenge https://github.com/ocean-data-challenges/2021a SSH mapping OSE. We use the data from the first five satellites as inputs for the mapping and the last one (Cryosat-2) for computing the performance metrics. We compute these metrics in the along-track geometry. The evaluation domain spans from (65°W, 33°N) to (55°W, 43°N) and the evaluation period from January 1 st to December 31 st 2017. Given η c2 and η the measured SSH and the reconstructed SSH respectively, we compute the following two metrics:

• µ ssh is a score based on the normalized root mean squared (nRMSE) error computed as 1 -RM S(η -η c2 ) RM S(η c2 )

• λ x is the wavelength at which the power spectrum density (PSD) score 1-P SD(η -η c2 ) P SD(η c2 ) crosses the 0.5 threshold, which characterize the scales resolved by the reconstruction (the error below that wavelength makes up for more than half of the total signal)

In Table 3, we also consider the root mean square error (RMSE) as well as the nRMSE score of the sea level anomaly µ sla obtained by subtracting the mean dynamic topography to the SSH. Lastly, we assess the performance degradation resulting from the transition from simulated to real data by quantifying the improvement relative to DUACS in the resolved scale λ x on our OSE setup as well as on the OSSE benchmarking setup proposed in [START_REF] Guillou | Mapping Altimetry in the Forthcoming SWOT Era by Back-and-Forth Nudging a One-Layer Quasigeostrophic Model[END_REF]. This benchmarking setup relies on NATL60-CJM165 OSSE dataset. We refer the reader to https://github.com/ocean-data-challenges/ 2020a SSH mapping\ NATL60 for a detailed description of this experimental setup.

Results

This section details our numerical experiments for the considered real altimetry casestudy for a Gulf Stream region as described in Section 3.5. We first report the benchmarking experiments to assess the performance of the proposed learning-based strategy with respect to (w.r.t.) state-of-the-art mapping schemes. We then analyse how the characteristics of the training datasets drive the mapping performance.

Benchmarking against the state of the art

We report in Table 2 the performance metrics of state-of-the-art approaches including both operational observation products [START_REF] Taburet | DUACS DT2018: 25 years of reprocessed sea level altimetry products[END_REF][START_REF] Ubelmann | Reconstructing Ocean Surface Current Combining Altimetry and Future Spaceborne Doppler Data[END_REF], deep-learning-based schemes trained on observation data [START_REF] Archambault | Multimodal Unsupervised Spatio-Temporal Interpolation of satellite ocean altimetry maps[END_REF][START_REF] Martin | Synthesizing Sea Surface Temperature and Satellite Altimetry Observations Using Deep Learning Improves the Accuracy and Resolution of Gridded Sea Surface Height Anomalies[END_REF] as well as methods using explicitly a model-based prior on sea surface dynamics [START_REF] Guillou | Mapping Altimetry in the Forthcoming SWOT Era by Back-and-Forth Nudging a One-Layer Quasigeostrophic Model[END_REF][START_REF] Ballarotta | Dynamic Mapping of Along-Track Ocean Altimetry: Performance from Real Observations[END_REF][START_REF] Lellouche | The Copernicus Global 1/12°Oceanic and Sea Ice GLORYS12 Reanalysis[END_REF]. We compare those methods with a 4DVarNet trained on eNATL60-0 OSSE dataset. The latter outperforms all other methods on the two metrics considered (22% improvement in RMSE w.r.t. the DUACS product and 33% improvement in the resolved scale). We report a significantly worse performance for GLORYS12 reanalysis. This illustrates the challenge of combining large ocean general circulation models and observation data for the mapping of the SSH.

The last column indicates that the 4DVarNet scheme leads to the best mapping scores for both the OSE and OSSE setups. For the latter, the reported improvement of 47% is twice greater than the second best at 22%. The performance of the 4DVarNet drops by 11% when considering the former. By contrast, other methods do not show such differences between the OSE and OSSE case-studies. This suggests that the finer-scale structures that are well reconstructed in the OSSE setup are not as beneficial in the OSE setup. While one could question the representativeness of the OSSE datasets for the finescale patterns in the true ocean, real nadir altimetry data may also involve multiple processes which could impede the reconstruction and evaluation of horizontal scales below 100km. 1; the subsequent columns indicate the reconstruction metrics described in Section 3.5. Note that the NATL60 could not be evaluated on the OSSE setup since the evaluation data were used for validation during the training stage. We analyse here in more detail the impact of the spatial resolution of the training dataset onto the reconstruction performance. In Table 3, as expected, the higher resolution grid in the ocean run simulation leads to a better mapping with a 22% improvement in λ x and a 17% improvement in the RMSE score between the experiments with the coarsest (ORCA025) and finest (NATL60) resolutions. We also observe qualitative differences in the relative vorticity fields in Figure 2. Residual artifacts due to the altimetry tracks appear (60°W, 39°N) for the two lower-resolution training datasets. They are greatly diminished when considering the NATL60 dataset. Despite these differences, the reconstructed vorticity and kinetic energy fields in Figure 2 look very similar for the different 4DVarNet schemes, whatever the training datasets. By contrast, the vorticity and kinetic energy fields in the training datasets clearly depict fewer fine-scale structures and weaker gradients for the lower-resolution simulation datasets, namely ORCA025 and GLORYS12-f. These results support the generalization skills of 4DVarNet schemes to map real altimetry tracks despite being trained on SSH sensibly different from the reconstruction.

Eddy-present datasets versus eddy-rich ones

We draw similar conclusions from the analysis of the spectral densities shown in Figure 4. The differences in the energy distribution of the training data significantly reduce in the reconstructions. 4DVarNet schemes trained from higher-resolution datasets however result in more faithful reconstruction at all scales. The patterns observed for the temporal PSD are slightly different in Figure 3. We do not observe the same homogenization as for the spatial PSD. Lower-resolution training datasets involve a significant drop of an order of magnitude for periods greater than 10 days and wavelength greater than 200km.

Forced simulation datasets versus reanalysis ones

Looking in more specifically at the effect of ocean reanalysis between the two experiments GLORYS12-f and GLORYS12-r. We can first note the impact of observation data assimilation in Figure 3 where we see how the power spectrum of the reanalysis is significantly raised compared to the free run. The spectrum is closer to ones of the higher resolution simulations. Visually we also clearly see stronger gradients in the kinetic energy in Figure 2. Quantitatively in Table 1 we see an improvement of 11% in both the RMSE and the scale resolved, besides training on a reanalysis increase the relative gain w.r.t. DU-ACS significantly more on real data (+9%) than on simulated data (+5%) as we can see in the right most column. This suggests that the reanalysis dataset conveys information on real world observations which improves the generalization performance.

Tide-free datasets versus tide-resolving ones

We assess here the impact of tide-resolving simulation used as training data. We use the twin eNATL60 runs eNATL60-t and eNATL60-0. Contrary to other runs, those simulations contain barometric and wind forcing, we therefore remove the Dynamic Atmospheric Correction [START_REF] Carrere | Major improvement of altimetry sea level estimations using pressure-derived corrections based on ERA-Interim atmospheric reanalysis[END_REF] from the SSH fields. Additionally since the barotropic tide signals are removed from real altimetry tracks prior to interpolation, we also remove the signal from the training data by subtracting the spatial mean over the training domain for each hourly snapshot before calculating the daily averages.

Given those processing steps, the two training datasets exhibit very similar wavenumber spectra as shown in Figures 3. We also find that training on those two datasets produce little differences in the reconstructions both quantitatively (see Table 3) and qualitatively (Fig. 2). The resulting performance is comparable to that of the NATL60 experiment.

We identify two hypotheses for explaining why tide-resolving simulation do not lead to better mapping schemes:

• The preprocessing applied on the training field remove the main tide signals. We therefore effectively measure the impact of tide modeling on other ocean processes that may be less significant; • The evaluation procedure applied on altimetry tracks on which the barotropic tide has been filtered may not be interpretable enough to measure the reconstruction of residual tide signals. New instruments like the KaRIN deployed in the SWOT mission may provide new ways to better quantify those effects.

These findings provide motivation for carefully considering the purpose of the learningbased model when making decisions about the training data. In our case, explicitly modeling tide processes that are removed from the observations in the evaluation setup added overheads in the computational cost of running the simulation as well as in the preprocessing of the training data. Additionally given the considered evaluation data and metrics, we were not able to quantify any significant differences between the two trained mapping schemes.

Discussion

This study has been greatly facilitated by the standardized tasks and evaluation setups proposed in data-challenges https://ocean-data-challenges.github.io/. Datachallenges are used to specify a targeted problem of interest to domain experts through datasets and relevant evaluation metrics. This preliminary work have been instrumental in constituting the comprehensive benchmark and combining methods from different teams and institution around the world. Additionally, it also constitutes a strong basis for a trans-disciplinary collaboration between the ocean and machine learning research communities.

Moreover, the results presented in this study introduce a use of ocean simulations for developing altimetry products. This opens new ways for ocean physicist, modelers and operational oceanographers to collaborate. In order to assess the range of these new synergies, it would be interesting to explore if the approach proposed here of training neural schemes using simulation data would generalize to other tasks such as forecast or sensor calibration and to other quantities like surface temperature, currents, salinity or biochemical tracers.

If the simulation-based training approach introduced here is successfully extended to other ocean problems, one could envision training large foundation deep learning models (Brown et al., n.d.) capturing the inner structure of high resolution ocean simulations which could then be used in many downstream applications. This could be the way to capitalize on all the advancement in ocean modeling without having to run OGCM numerical simulation for each downstream products. Furthermore, we would like to highlight the cost consideration when running numerical simulation intended for training learning based schemes. Indeed given that the eNATL60 run took 2700x CPU hours and 350x memory compared to the ORCA025 run for a smaller domain, a trade-off arises between generating multiple "cheap" trajectories versus generating a single more realistic trajectory.

To conclude, we have shown in this study that training machine learning models on simulations datasets leads good performance on real altimetry data mapping and outperforms current state of the art approaches. The model trained on NATL60 reduces the RMSE by 18% compared neural schemes trained on observation data and improves the scales resolved by 33% compared to the DUACS operational product. Even the coarsest simulation considered ORCA025 provides competitive results with current operational methods. We have shown that using a more realistic SSH fields using reanalysis or higher resolution simulations increases the performances of the trained model. This is an exciting result that shows the potential for training operational products from ocean simulations and how advances in ocean modeling in operational oceanography can be beneficial. The results shown here are limited to the interpolation problem on a regional domain but the robustness of the performance shown are encouraging for further developing these results using a larger domain.

Figure 2 :

 2 Figure 2: Samples Kinetic energy and relative vorticity of the training and reconstruction data of January 6 t h . The reconstructed year is 2017 while the training year vary depending on the simulation. The first two columns (a) and (b) show the training data while columns (c) and (d) show the associated 4DVarNet reconstruction. The kinetic energy is displayed in columns ((a) and (c)) and the relative vorticity normalized by the local Coriolis parameter in columns ((b) and (d)). Each row shows the experiment using respectively: ORCA025 (I), GLORYS12-f (II), GLORYS12-r (III), NATL60 (IV), eNATL60-t (V) and eNATL60-0 (VI)

Figure 3 :

 3 Figure 3: Space-time spectral densities of the training datasets (first row) and of their associated reconstruction (second row). Darker blue in the lower left corner indicates higher energy at larger wavelength and periods. The different SSH fields exhibit different energy cascades when moving to finer temporal (upward) or spatial (rightward) scales.

Figure 4 :

 4 Figure 4: Spectral analysis of the training and reconstructed SSH datasets. We display the PSD of the training dataset (left plot), reconstructed SSH field (center plot) as well as the associated PSD score (right plot)

Figure 5 :

 5 Figure 5: Spectral impact of model reanalysis. We display the PSD of the training dataset (left plot), reconstructed SSH field (center plot) as well as the associated PSD score (right plot)

Table 1 :

 1 Summary table of the different synthetic SSH fields used for training.

	•	No	No	No

Table 2 :

 2 SSH reconstruction performance of the benchmarked methods (a) 4DVarNet from this study trained on eNATL60-0 (b) Archambault et al. (2023), (c and d) ConvLstm-SST and ConvLstm from Martin et al. (2023), (e) DYMOST from Ballarotta et al. (2020), (f) MIOST from Ubelmann et al. (2021), (g) BFN-QG from Guillou et al. (2021), (h) DUACS from Taburet et al. (2019), (i) GLORYS12 fromLellouche et al. (2021. The columns indicate from left to right: whether athe mapping schemes rely only on SSH data or also exploit additional data such as gap free SST products; if the method uses deep learning architectures; the data used to calibrate (or train) the mapping scheme; the numerical model of the ocean used for the mapping if any (QG stands for quasi-geostrophic); µ and λ x are the metrics as described in Section 3.5

		SSH Deep	Calibrated on Physical rmse µ ssh	λ x	1 -λx λ ref
		Only Learning data from	Model	(cm)	() (km) (% ose, osse)
	(a) 4DVarNet	Yes	Yes	Simulation	-	5.9 0.91	100	33, 47
	(b) MUSTI	No	Yes	Satellite	-	6.3 0.90	112	26, 22
	(c) ConvLstm-SST No	Yes	Satellite	-	6.7 0.90	108	28, -
	(d) ConvLstm	Yes	Yes	Satellite	-	7.2 0.89	113	25, -
	(e) DYMOST	Yes	No	Satellite	QG	6.7 0.90	131	13, 11
	(f) MIOST	Yes	No	Satellite	-	6.8 0.90	135	11, 10
	(g) BFN-QG	Yes	No	Satellite	QG	7.6 0.89	122	19, 21
	(h) DUACS	Yes	No	Satellite	-	7.7 0.88	151	0, 0
	(i) GLORYS12	No	No	Satellite	NEMO	15.1 0.77	241	-60, -

Table 3 :

 3 Performance of 4DVarNet mapping schemes trained on different simulated datasets. The first column shows the source of the training dataset as described in Table

Acknowledgments

This work was supported by ANR Projects Melody and OceaniX and CNES. It benefited from HPC and GPU resources from GENCI-IDRIS (Grant 2020-101030) and Ifremer.

Open Research Section

The authors provide the training data, source code, reconstructed maps and trained model for each experiments of the manuscript at https://doi.org/10.5281/zenodo.8064114.