Radu Cârpa
email: radu.carpa@ens-lyon.fr

Olivier Glück
email: olivier.gluck@ens-lyon.fr

Laurent Lefevre
email: laurent.lefevre@inria.fr

Segment Routing based Traffic Engineering for Energy Efficient Backbone Networks

Keywords: energy efficiency, backbone networks, SDN, traffic engineering, SPRING, segment routing, MPLS

Energy consumption has become a limiting factor for deploying large-scale distributed infrastructures. This work 1 seeks to improve the energy efficiency of backbone networks by providing an intra-domain Software Defined Network (SDN) approach to selectively turn off a subset of links. We propose the STREETE framework (SegmenT Routing based Energy Efficient Traffic Engineering) that dynamically adapts the number of powered-on links to the traffic load. The core of the solution relies on SPRING, a novel protocol being standardized by IETF. It is also known under the name of Segment Routing. The algorithms have been implemented and evaluated using the OMNET++ simulator. Experimental results show that the consumption of 44% of links can be reduced while preserving good quality of service.

I. INTRODUCTION

Networks play an important role in today's highly connected world, and their traffic is expected to increase by a factor of three in the next five years 2 . Although they are responsible for only a fraction of the energy the IT sector consumes, networks can be made more efficient as their components often consume a constant amount of power regardless of their utilization. Backbone networks, for example, comprise devices that consume several KWh of energy even when idle.

Techniques to reduce the energy consumption of backbone networks include traffic rerouting and the use of low power consumption modes. It is possible to reroute data over alternative paths to offload a subset of links because operators over-provision their networks to handle demands during peakperiods and to offer high Quality of Service (QoS) to their clients. Evidence shows, however, that even during peak hours links are rarely used more than 50% [START_REF] Hassidim | Network utilization: The flow view[END_REF]. Fig. 1c demonstrates the utilization of links computed from traffic matrices of the Géant network where most links are less than 25% utilized even, at peak hours.

In this work, we tackle the problem of energy consumption of backbone networks by changing the status of router ports and transponders on the two extremities of a link. The status of these components is set to sleep mode whenever a link is not required to transfer data, and brought back to operational state when needed. Hereafter we term this process as switching off/on links.

We analyzed the implementation issues of an energyefficient SDN-based traffic engineering in core networks. Most works which seek to adapt the number of switched on links to the amount of traffic in the network rest at the conceptual level. They search to minimize the number of active links via solving MILP formulations or heuristics. Few went deeper into implementation details, but still miss-observed or ignored some race conditions. Our team fully implemented the proposed solution in the OMNET++ packet-based discrete event simulator. Experiments considering real network topologies (Germany50 and Géant, see Fig. 1) and real dynamic traffic matrices allowed us to quantify the trade-off between energy saving and impact of our solution on network performance.

As mean to reroute the traffic we use a promising new protocol, SPRING. This comes in contrast with other works, which use classical IP link weights changes or MPLS+RSVP-TE for this purpose. SPRING proved itself well suited for dynamic reconfiguration of the network. As a bonus, the usage of SPRING makes the implementation easier. To the best of our knowledge, SPRING has never been considered in the literature for reducing the energy consumption of networks.

The implemented solution is proposed under a form of a framework which we name STREETE (SegmenT Routing based Energy Efficient Traffic Engineering). It represents an online method to switch some links off/on dynamically according to the network load. The generecity of our framework allows gradual enhancements, like adding protection and more intelligent flow placement computations by a simple update of its composing algorithms.

The rest of this paper is structured as follows. Section II discusses related work, whereas Section III describes our assumptions and concepts such as SDN and segment routing/SPRING protocol. After that, we present our solution in Section IV and analyze simulation results in Section V. Finally, Section VI concludes the paper.

II. RELATED WORK

Gupta et al. [START_REF] Gupta | Greening of the Internet[END_REF] were among the first to consider putting interfaces into sleep mode to save energy. In backbone networks, however, links cannot be simply switched off without losing data. It has been shown that even micro-second sleeps are difficult to achieve due to the small inter-packet arrival time [START_REF] Bolla | Energy efficiency in optical networks[END_REF]. To cope with this challenge, data transmission in bursts has been considered [START_REF] Nedevschi | Reducing Network Energy Consumption via Sleeping and Rate-adaptation[END_REF]; a technique nearly impossible under modern 100Gbps speeds.

A distributed network-wide solution has been provided by Vasic et al. in order to change flow paths to reduce link utilization, enabling 21% of links of real networks to be put in sleep mode [START_REF] Vasić | Energy-aware Traffic Engineering[END_REF]. Such a solution, however, is a concept whose implementation details are not fully disclosed. Another distributed approach considers that every node locally monitors the utilization of adjacent links and decides whether to switch them off or not [START_REF] Bianzino | GRiDA: A green distributed algorithm for backbone networks[END_REF]. Machine learning mechanisms are used to avoid choosing links whose extinctions have already provoked a congestion. This solution heavily relies on the Interior Gateway Protocol Traffic Engineering extensions (IGP-TE) and the communication overhead is underestimated. After detailed analysis, we discovered that due to bad estimations of the frequency of the IGP-TE exchanges [START_REF] Salsano | Traffic Engineering with OSPF-TE and RSVP-TE: Flooding Reduction Techniques and Evaluation of Processing Cost[END_REF], the announced overhead of 0.52% is actually a 30% increase of IGP-TE flooding. This is not negligible because these floodings are costly both in number of messages and required processing power. Moreover, the authors used classical IP routing instead of MPLS, which can induce instability periods while the network is re-configuring.

Zhang et al. proposed a centralized approach based on MPLS+RSVP-TE [START_REF] Zhang | GreenTE: Power-aware traffic engineering[END_REF] that relies on estimation of traffic matrices, proven not very efficient [START_REF] Gunnar | Traffic Matrix Estimation on a Large IP Backbone: A Comparison on Real Data[END_REF]. They also opt for a hybrid approach where MPLS+RSVP-TE is used in parallel with shortest path routing. Our work uses the SPRING protocol, which yields the same advantages with less complexity.

Existing work also concentrates on searching the best links to switch off assuming full knowledge of network state. The results on energy efficiency of the proposed solutions varies drastically from one work to another. For example, at times the performance of greedy approaches like "selecting the least congested link" is claimed to be close to that of optimal algorithms in terms of energy savings (Kamola et al. [START_REF] Kamola | Shortest path green routing and the importance of traffic matrix knowledge[END_REF]). In contrast, stochastic optimization techniques, like genetic algorithms, are said to outperform greedy algorithms when considering stability of both network and energy savings (Bonetto et al. [START_REF] Bonetto | Algorithms for the multi-period power-aware logical topology design with reconfiguration costs[END_REF]). We attempt to explain this variation in our work.

[4] presents ESTOP, an algorithm which is similar to STREETE. They also test their solution on the Germany50 network. However, they do not exploit the knowledge of the network wise traffic matrices and obtain worst results in terms of energy efficiency.

III. ASSUMPTIONS AND BACKGROUND

We consider that the variation of traffic in backbone network is slow and follows a diurnal pattern due to a high level of flow aggregation [START_REF] Juva | Traffic characterization for traffic engineering purposes: analysis of Funet data[END_REF]. Table I provides figures on the power consumed by backbone-network equipments [START_REF] Van Heddeghem | Power consumption modeling in optical multilayer networks[END_REF]. Our model considers that switching a port off corresponds to powering off part of its integrated circuit and the transponder. Switching a link off means to switch off the two ports at the link's extremities. Hence, if a 100Gbps link is switched off, 2 • (135 + 150) = 570 Watts are saved.

Table I

A. SDN

In traditional networks, control and data planes are located on network nodes, and adjustments to network-wide policies often require access to, and configuration of, all devices. SDN aims to simplify management by separating the control and data planes. The control plane contains the network intelligence responsible for calculating paths to be followed by data flows and programming them into the data plane. The data plane is responsible for passing packets from an incoming port of a device to the right outgoing port. We refer to data plane devices as SDN switches.

SDN controllers exploit network-wide knowledge to execute all required computations and apply changes into the data Compared to the data-center networks, the optical backbone networks face problems (such as lightpath provisioning), for which the current SDN implementation, OpenFlow, is not ready. However, researchers succeed to bring the SDN concept into the backbone [START_REF] Siqueira | An optical SDN Controller for Transport Network virtualization and autonomic operation[END_REF]. In this paper we assume the availability of such a control plane, whose comportment we can change for our needs: for example, by modifying its source code.

We focus on a centralized SDN solution. While it is possible to chose the best links to be shut down in a distributed fashion [START_REF] Bianzino | GRiDA: A green distributed algorithm for backbone networks[END_REF], more complex and intelligent decisions can be taken if done in a centralized manner.

B. SPRING protocol

SPRING [START_REF] Filsfils | Segment Routing Architecture[END_REF], also known as Segment Routing, is an IETF draft since 2013 and has the goal to substitute MPLS+RSVP-TE for traffic engineering. It combines the power of source routing (flexible traffic engineering) with shortest path routing (less signaling/header overhead compared to MPLS/Myrinet).

The data plane of the SPRING protocol uses the same concept of label switching of MPLS and can even use existing MPLS-capable devices without modification. The control plane though, suffered a complete re-design. Called Segment IDentifiers (SID), labels are global identifiers, unlike classical MPLS where they are local. Here, we are interested in 2 types of SIDs, namely "nodal" and "adjacency" (Fig. 2a).

• A nodal SID is globally unique and identifies a node (a, b, c, ..., h). • An adjacency SID is local to a node and uniquely identifies an outgoing interface (node b has the Adjacency SIDs L1, L2, L3, L4 and L5).

Compared to MPLS, the distribution of labels is done via an extension to the IGP instead of using special protocols (LDP/RSVP-TE). After the network discovery, sending a packet to the node a through the shortest path resumes to encapsulating it into a packet with destination a. However, unlike in IP, much more flexible traffic engineering is possible:

• If node h wants to send a packet to node a while forcing it over link L1, it adds the header [b, L1] (Fig. 2b). As a result, the packet will take the shortest path to b, and b will force it through link L1 to a.

• In Fig. 2c, if h wants to send a packet to a via f , it will add the header [f, a]. As result the packet will take the shortest path to f , and afterward continue to a. SPRING is a good choice for dynamic rerouting regardless of opting for an SDN approach or not. Being a source routing protocol, it enables fast flow setup and easy reconfiguration of virtual circuits with minimum overhead. Changes must be applied only on the ingress devices and no time and signaling is lost re-configuring the midpoint devices. This is important when frequent routing policy changes are needed.

We also use SPRING protocol because MPLS+RSVP-TE faces convergence issues in case of link failures 3 . With small numbers of RSVP-TE virtual tunnels and the refresh reduction extension added to RSVP-TE, the risk has been reduced. However, in our case, due to rerouting the virtual tunnels over a subset of paths and switching other links off, we increase the load on some routers and the number of LSPs per router. For this reason, SPRING suits better our work because it is more stable with a larger number of dynamically changing LSPs.

IV. THE STREETE FRAMEWORK

A. Analysis

In this section we formalize the problem of the energy efficient traffic engineering. The goal is to avoid losing connectivity between any two nodes in the network. The term Connectivity Constraint is used to refer to satisfying the integrity of the network.

We divided our solution into the following three steps:

• Selecting links to switch off/on.

• Computing new routes to avoid/re-use these links.

• Pushing these routes to the network devices and switching off/on the links. 1) Selecting links to switch off/on: This step selects candidate links to be switched off/on. We implement the "least congested link" technique because it is the approach against which most of existing work compare their results. This step, quite straightforward, first sorts links by their instant transmission speed. The complexity of this step is

O(E • log(E))
, where E is the number of links (Edges). It then tries to switch off links that transmit less data by testing, for each link, whether shutting it down does not violate the connectivity constraint. At most it will take O(E •(E +V )) : E times breadth-first search to test the integrity of the network, where V is the number of nodes (Vertices). The algorithm that searches for good links to be switched off/on is hereafter referred to as <SelectLinksToOff/On>.

2) Computing new routes for data flows: Prior to actually switching off/on the links selected by <SelectLinksToOff/On>, we must compute new paths for every flow traversing them. It may happen that the network is not capable of routing all data flows if all the selected links are switched off. Moreover, letting every node recalculate new shortest paths throughout the network is not efficient and may produce congestion. In fact, congestion may occur even when a link is turned on, such as when the turned-on link provides a shortest path for a large number of flows.

To compute the new routes, the SDN controller exploits global knowledge about the network traffic to solve the "Multicommodity Flow Problem" which has been proven to be NPcomplete, but for which approximation techniques exist [START_REF] Madry | Faster Approximation Schemes for Fractional Multicommodity Flow Problems via Dynamic Graph Algorithms[END_REF]. As a prerequisite, the SDN controller has to know the used bandwidth between any head-end to any tail-end (the traffic matrix). The statistics must be collected on the ingress SDN switches. Fortunately, the availability of this kind of metric is already considered in the OpenFlow SDN specification. At this step, the algorithm may also compute protection paths [START_REF] Victor | Spare Capacity Allocation Using Shared Backup Path Protection for Dual Link Failures[END_REF].

In our simulation we implemented a technique that computes the all-pairs shortest path routes and tests the possibility to route all the available traffic over these paths. The complexity of the computation is O(V • (E + V • logV )) if Dijkstra from every node is computed. The complexity may decrease with the use of an optimized all-pairs shortest path algorithm.

In the case of switching on, our algorithm is let as simple as possible. It tests the utilization of the most utilized link, and if it is greater than 75%, everything is switched on the shortest path routes are recalculated. We refer to this algorithm as <ComputeNewRoutes>.

3) Rerouting and switching links off/on: Once new routes for flows are computed, the controller directly informs each node about the explicit routes to be introduced in the forwarding table and the links will be switched off or on.

The steps taken to switch off the links are described in the Algorithm 1. For example, here the links between nodes A and B, C and D, etc were selected for shutdown. The algorithm presents the steps to be executed to switch these links off.

The overhead of this solution in terms of exchanged messages is smaller than for distributed solutions. In distributed cases, the nodes need to synchronize their actions to avoid switching off a link which disconnects the network. In the centralized SDN case, if network failures do not occur, there is no risk to violate the connectivity constraint. The cases of network failures may be treated by adding protection constraints to <SelectsLinksToOff> and <ComputeNewRoutes>. We do not address this issue in this paper.

4) Main Loop: <MainLoop> demonstrates how the presented algorithms interact. This algorithm searches links that can be switched on or off. If it finds such links, it will actually execute <RerouteAndSwitchOn> or <RerouteAndSwitchOff>. The controller tells every router the new paths to be used for virtual circuits;

3 3
The controller tells nodes adjacent to links (A,B,...) to turn them off; The SDN controller sends IGP flooding : link A-B is down / * flooding is needed, because at line 2 we may not contact every node * /

6

Every router applies the changes to the forwarding table to reroutes over the new paths;

A. Environment Setup

The STREETE approach was implemented in the OM-NET++ simulator. We also developed a draft version of the SPRING protocol to dynamically reroute data in the network to avoid turned-off links. At the border of the network domain, during the classification of each packet and assignment to a virtual circuit, the ingress nodes collect the statistics about the amount of data sent to all the other nodes inside the domain. The SDN controller retrieves this information from every router and generates the network-wide traffic matrix to be used by the algorithm presented in the previous section.

We tested our solution considering two different networks, summarized in Table II. The Germany50 network is a real backbone, i.e. the traffic is the aggregation of a large number of flows and is steady over the time. Although Geant is interconnecting research networks, the aggregation degree is high enough and the traffic pattern is similar to other backbone networks. We still observe slightly higher traffic variations, but that did not have any impact on our solution.

We tested STREETE by simulating real networks with real traffic matrices from SNDLib4 , which were scaled to represent today's faster network speeds. The links in the network are switched off one by one while no link reaches the threshold utilization of 60%. Over this value, we stop switching off links. When the utilization of any link goes beyond 75% for a couple of minutes, the network is switched back on.

The propagation delay was calculated as the time needed for light to travel through the fiber. The distance between nodes is the bird-fly distance. Queuing delay is simulated by OMNET++.

B. Results

1) Energy consumption:

We simulated a day of operation for each network. In Fig. 3, green lines represent the % of links left on by our algorithm. Using the values of Table I for 100Gbps links, this translates into reductions in energy consumed during the simulated day as shown in Table III.

At the beginning of the day, when our algorithm starts to execute, both networks start to shut down links one by one every couple of minutes and converge very fast to a spanning tree. The further extinction of links is not possible without violating the connectivity constraint. In the case of the Germany 50 network, we observe a turn-on at 23:00h.

It is due to a large increase in traffic demand between two network nodes. We do not know the reason for this unusual increase in the original traffic matrices, but this demand was handled well by our algorithm which turned on sleeping links and restored the full capacity of the network. We discovered that even greedy techniques for selecting links to switch off provide good energy savings due to high over-provisioning of the analyzed networks. An a posteriori detailed analysis of the work which affirms the contrary ( [START_REF] Bonetto | Algorithms for the multi-period power-aware logical topology design with reconfiguration costs[END_REF]) on the same traffic matrices showed that they compare a greedy approach acting at the physical link granularity with complex techniques acting on the light path granularity.

2) Impact on link loads: The fact of rerouting data over paths which are not the shortest possible, induces overhead in the network load. The overhead of our solution is visible in Figures 3a and3c. The maximum overhead observed during these simulations was 18.56%. It may be noticed that even with this overhead, the networks converge to a minimum energy state.

At the beginning of the simulations and the moment when the network was turned on, in the case of Germany50 we observe that the overhead is null, as expected.

3) Impact on end-to-end delays: Figures 3b and3d show the impact of our solution on the end-to-end delay. For every source and destination node, we calculated the mean end-toend delay per intervals of one hour. The box plots represent the distribution of the means of all pairs (source,destination) during the previous hour. The blue ones describe the delays in the original network and the red ones the delay in the network with switched off links. We also represent the mean delay for the cases with switch off and without switch off as continuous lines of the same colors.

In the case of the Germany50 network, the delay increases by 37% in the worst case, but usually increases by around 20%. A much more interesting result can be seen in the Géant network. While in most cases the delay still slightly increases, it can be clearly seen that the maximum outlier suffered a delay reduction. It is due to the minimum hop count shortest path routing. The Géant network has a complex physical topology and it occurs that the minimum hop count routing is not well suited for this network: in some cases traffic between two European nodes will pass through New York. Turning off the 2 links forced the flows to pass through multiple hops, but remaining in Europe.

4) Impact on packet loss: The simulation revealed absolutely no packet loss due to our algorithms. In the case of backbone networks, the aggregated traffic flows do not change fast enough to produce congestion implied by our solution. The algorithm wakes up the links and serves the increased demand.

VI. CONCLUSION

In this work, we implemented an energy-aware traffic engineering technique for reducing energy consumption in backbone networks. Energy efficient traffic engineering was analyzed in previous work, but none addressed implementation challenges of their solutions. We showed that ignoring to test the feasibility of techniques can lead to bad estimations and unstable solutions. We proposed and implemented a working prototype in the OMNET++ simulator.

Nowadays networks are changing towards a centralized architecture and SDN is increasingly utilized in data center production networks. We believe that the SDN philosophy may be extended to the backbone networks. The simplicity of the implemented solution shows that SDN may also be a good path for reducing the energy consumption of the network devices.

Compared to previous work, in this paper we used the SPRING transport protocol to improve the stability of energy efficient traffic engineering solutions. To the best of our knowledge, this is the first work proposing the use of SPRING to improve the energy efficiency of backbone networks. The flexibility of this routing protocol is well suited to frequent route changes that happen when we switch links off and on. Moreover, this protocol is very SDN-friendly.

Via simulations, we showed that as much as 44% of links can be switched off to save energy. We also discovered by using real traffic matrices that backbone networks are very over-provisioned. Even greedy techniques can easily approach the maximum reduction in the amount of energy consumed. In fact, the bottleneck in terms of energy efficiency in energy aware traffic engineering solutions is the connectivity constraint. While this may not hold true in all networks, the analyzed ones confirmed the assumption.

In future work, we will evaluate the solution on a real NetFPGA testbed simulating a backbone network.

Figure 1 :

 1 Figure 1: Considered network topologies and Géant usage

Figure 2 :

 2 Figure 2: SPRING protocol

Data:

  List-Of-Links :[(A-B),(C-D),...] and the new routes / * given by <SelectLinksToOff> and <ComputeNewRoutes>

4 4

 4 Nodes A,B,C,... ACK to the controller ;5 5

7 8 Algorithm 1 : 7 if

 817 Nodes A,B,C,... suspend the IGP HELLO messages ; Let the links (A-B,C-D,...) go to a sleep mode / negotiate the shutdown ; 9 end Centralized SDN: <RerouteAndSwitchOff> 1 while Energy Efficient Traffic Engineering is active do 2 List-Of-Links = <SelectLinksToOn>(); 3 if <ComputeNewRoutes>(All-links-that-are-ON ∪ List-Of-Links) then 4 <RerouteAndSwitchOn>(List-Of-Links); <ComputeNewRoutes>(All-links-that-are-ON \ List-Of-Links) then 8 <RerouteAndSwitchOff>(List-Of-Links);

  our solution's overhead Initial network load % of links in Minimum Spanning Tree % of links ON (a) Impact on link loads (b) Impact on end-to-end delay our solution's overhead Initial network load % of links in Minimum Spanning Tree % of links ON (c) Impact on link loads (d) Impact on end-to-end delay

Figure 3 :

 3 Figure 3: Simulation results in Germany50 (top) and Géant (bottom) networksTable II: Evaluated network topologies.NetworkNumber of nodes Number of links Géant (Fig.1b) 22 * 36 Germany50 (Fig.1a) 50 88 *the New York node is not visible in Fig.1b, but the two links connecting it to the rest of the network are shown.

Table III :

 III Energy savings during the simulated day.

	Network	Average off time of links	Energy economy (kWh)
	Géant	14.88 out of 36	204
	Germany50	38.18 out of 88	522

The work is done as part of the CHIST-ERA STAR "SwiTching And tRansmission" European project (http://www.chistera.eu/projects/star), which aims to decrease the consumption of backbone networks by a factor of 100 until

2 http://ciscovni.com/forecast-widget/index.html

http://www.ietf.org/mail-archive/web/bmwg/current/msg01502.html

http://sndlib.zib.de/


Pavel Sakov
email: pavel.sakov@bom.gov.au

Jean-Matthieu Haussaire

Marc Bocquet

J.-M Haussaire

An iterative ensemble Kalman filter in presence of additive model error

Keywords: ensemble Kalman filter, model error, Gauss-Newton minimisation, iterative ensemble Kalman filter

The iterative ensemble Kalman filter (IEnKF) in a deterministic framework was introduced in Sakov et al. (2012) to extend the ensemble Kalman filter (EnKF) and improve its performance in mildly up to strongly nonlinear cases. However, the IEnKF assumes that the model is perfect. This assumption simplified the update of the system at a time different from the observation time, which made it natural to apply the IEnKF for smoothing. In this study, we generalise the IEnKF to the case of imperfect model with additive model error.

The new method called IEnKF-Q conducts a Gauss-Newton minimisation in ensemble space. It combines the propagated analysed ensemble anomalies from the previous cycle and model noise ensemble anomalies into a single ensemble of anomalies, and by doing so takes an algebraic form similar to that of the IEnKF. The performance of the IEnKF-Q is tested in a number of experiments with the Lorenz-96 model, which show that the method consistently outperforms both the EnKF and the IEnKF naively modified to accommodate additive model noise.

Introduction

The analysis step in the Kalman filter (KF, [START_REF] Kalman | A new approach to linear filtering and prediction problems[END_REF]) can be seen as a single iteration of the Gauss-Newton minimisation of a nonlinear cost function [START_REF] Bell | The iterated Kalman smoother as a Gauss-Newton method[END_REF]. It yields an exact solution in the linear case and works well in weakly nonlinear cases, of the system and improves the forecast (the background) and forecast covariance (the background covariance) used in calculating the analysis at the next iteration. Despite processing the same observations multiple times, the IEnKF maintains the balance between the background and observation terms in the cost function: each next iteration represents a correction to the previous one rather than a new assimilation of the same observations. It is different in this respect from the Running in Place scheme (RIP, [START_REF] Kalnay | Accelerating the spin-up of Ensemble Kalman Filtering[END_REF][START_REF] Yang | Handling nonlinearity in the Ensemble Kalman Filter: Experiments with the three-variable Lorenz model[END_REF], which adopts the latter approach. The RIP also has a stochastic implementation [START_REF] Lorentzen | An iterative ensemble Kalman filter[END_REF]. A Bayesian derivation of the IEnKF, which suggests its optimality for nonlinear chaotic models, has been given in section 2 of [START_REF] Bocquet | An iterative ensemble Kalman smoother[END_REF].

The perfect model framework makes it possible to extend the IEnKF for assimilating future observations, or smoothing. The 40 corresponding method is known as the iterative ensemble Kalman smoother (IEnKS, [START_REF] Bocquet | Localization and the iterative ensemble Kalman smoother[END_REF]Sakov 2014, 2013). The IEnKF can also be enhanced to accommodate the inflation-less EnKF (IEnKF-N, [START_REF] Bocquet | Combining inflation-free and iterative ensemble Kalman filters for strongly nonlinear systems[END_REF], and the ensemble space formulation of the IEnKF algorithm makes it possible to localise it [START_REF] Bocquet | Localization and the iterative ensemble Kalman smoother[END_REF] with the localisation method known as the local analysis [START_REF] Evensen | The Ensemble Kalman Filter: theoretical formulation and practical implementation[END_REF]). Moreover, it is possible to base the iterative EnKF on minimisation methods other than the Gauss-Newton, e.g., on the Levenberg-Marquardt method [START_REF] Bocquet | Combining inflation-free and iterative ensemble Kalman filters for strongly nonlinear systems[END_REF][START_REF] Chen | Levenberg-Marquardt forms of the iterative ensemble smoother for efficient history matching and uncertainty quantification[END_REF].

Along with the listed above single-cycle iterative schemes, there are also a variety of multi-cycle iterative EnKF methods, emerging mostly from applications with static or quasi-static model state, such as oil reservoir modelling (e.g., [START_REF] Li | Iterative ensemble Kalman filters for data assimilation[END_REF]. Such methods involve re-propagation of the system from the initial time using the last estimation of the static parameters of the model. They are less suitable for applications with chaotic (e.g., with atmospheric or oceanic) models, when the divergence at any single cycle can be typically considered as a crash of the system.

The additive model error can be straightforwardly included into Monte-Carlo, or stochastic formulations of either iterative or non-iterative EnKF schemes. While it has not been formally considered in the original EnRML [START_REF] Gu | An iterative ensemble Kalman filter for multiphase fluid flow data assimilation[END_REF], it was a part of the iterative ensemble smoother by [START_REF] Mandel | Hybrid Levenberg-Marquardt and weak-constraint ensemble Kalman smoother method[END_REF].

Despite the intensive developments of the deterministic iterative EnKF schemes, so far they have not rigorously included the model error. One reason for that is the simplicity of the asynchronous DA in the perfect model framework [START_REF] Evensen | An ensemble Kalman smoother for nonlinear dynamics[END_REF][START_REF] Hunt | Four-dimensional ensemble Kalman filtering[END_REF][START_REF] Sakov | Asynchronous data assimilation with the EnKF[END_REF]. The other reason is that the model error increases the dimension of the minimisation problem: if in the perfect model case the optimal model state at any particular time defines the whole optimal model trajectory, with a non-perfect model the global in time optimal solution represents a set of optimal model states at each DA cycle. This complicates the problem even in the simplest case of sequential DA considered in this study.

The development of an IEnKF framework with imperfect model can have a number of important theoretical and practical implications. Firstly, it can help understand limits of applicability of the perfect-model framework and limitations of empirical treatments of the model error. Further, it would be interesting to see whether/when adding empirically the model error term to the cost function can have a regularising effect similar to that of the transition from the strong constraint 4D-Var to the weakconstraint 4D-Var. Including the model error has the potential to successfully address situations when a large model error can be expected, such as of a probable algal bloom in biogeochemical models, or of a rain event in land models.

This study develops an iterative method called IEnKF-Q based 90 on the Gauss-Newton minimisation in the case of a system with additive model error. In the following, the non-Gaussianity of the data assimilation system originates from the nonlinearity of the model dynamics and the observation frequency. This data assimilation system is assumed to lie in a range from 95 a weakly nonlinear to a strongly nonlinear regime, where the EnKF might fail but where multimodality of the underlying cost function is still not prominent. Nonetheless, by construction, the IEnKF-Q could also accommodate nonlinear observation operators and to some extent, which is context-dependent, non-100 Gaussian variables due to its variational analysis. The strongly nonlinear regime where multimodality becomes prominent and where the iterative ensemble Kalman filter and smoother could fail has been discussed in the conclusions of [START_REF] Bocquet | An iterative ensemble Kalman smoother[END_REF]. We refer to [START_REF] Fillion | Quasi static ensemble variational data assimilation[END_REF] for a more complete study 105 of this strongly nonlinear regime but in a perfect model context.

The outline of the study is as follows. The IEnKF-Q method is formulated in section 2. If the observation operator is linear, an alternative formulation resulting in the decoupling into a smoothing and a filtering analysis is discussed in section 3. A 110 pseudo-code for the IEnKF-Q algorithm is presented in section 4, and its performance with the Lorenz-96 model is tested in section 5. These tests include a preliminary study of a local IEnKF-Q. The results are discussed in section 6 and summarised in section 7.

Formulation and derivation

The IEnKF-Q method is introduced from the more general context of the following global in time cost function J K :

{x i } K i=1 = arg min {xi} K i=1 J K (x 1 , . . . , x K ), (1a)
J K (x 1 , . . . , x K ) = 1 2 x 1 -x f 1 2 (P f 1 ) -1 + 1 2 K i=2 y i -H i (x i ) 2 R -1 i 120 + 1 2 K i=2 x i -M i (x i-1 ) 2 Q -1 i . ( 1b
)
Here i is the cycle number associated with In the case of linear M and H the problem (1) becomes quadratic and has recursive solutions. The last component of the solution x K is known as the filtering analysis and is given by the KF, while the whole analysis {x i } K i=1 is given by the Kalman smoother. In the nonlinear case, it is essential for the applicability of recursive, or sequential, methods based on the KF, such as the EnKF, that the nonlinearity of the system is weak.

The rationale for iterative methods such as the IEnKF is that the weak nonlinearity needs to be achieved only in the course of minimisation; then the final analysis is calculated essentially for a linear system.

Therefore, the focus of an iterative method is a single analysis cycle (i.e. K = 2) and the associated problem that arises in the course of the iterative solution of

{x 1 , x 2 } = arg min {x1,x2} J(x 1 , x 2 ), (2a)
J(x 1 , x 2 ) = 1 2 x 1 -x a 1 2 (P a 1 ) -1 + 1 2 y 2 -H(x 2 ) 2 R -1 + 1 2 x 2 -M(x 1 ) 2 Q -1 . (2b)
Here we have dropped the absolute time indices and use relative indices 1 and 2, which refer to analysis times x = E1/m, (3a)

P = AA T /(m -1), (3b)
A

≡ E -x1 T = E (I -11 T /m), ( 3c
)
165
where m is the ensemble size, and 1 is a vector with all components equal to 1.

The problem formalised by (2) can be solved by finding zero gradient of the cost function (2b):

     ∇x 1 J(x 1 , x 2 ) = 0, ∇x 2 J(x 1 , x 2 ) = 0, (4)
similarly to the approach in [START_REF] Sakov | An iterative EnKF for strongly nonlinear systems[END_REF]. However, for an ensemble-based system the derivation becomes simpler if the solution is sought in ensemble space. Let us assume

x 1 = x a 1 + A a 1 u, ( 5a
) 175 A a 1 (A a 1 ) T = P a 1 , (5b)
A a 1 1 = 0, (5c)
where A a 1 is defined as the matrix of the centred anomalies resulting from a previous analysis at t 1 , and

180 x 2 = M(x 1 ) + A q 2 v, (6a)
A q 2 (A q 2 ) T = Q, (6b)
A q 2 1 = 0.

(6c)

We seek solution in (u, v) rather than in (x 1 , x 2 ) space. Note that 185 for convenience we use a different normalisation of the ensemble anomalies in (5b) than in (3b). After substituting (5,6) into (2) the problem takes the form

{u , v } = arg min {u,v} J(u, v), (7a)
J(u, v) = 1 2 u T u + 1 2 y 2 -H(x 2 ) 2 R -1 + 1 2 v T v, ( 7b
) 190 or, concatenating u and v, w ≡    u v    , (8a)
w = arg min w J(w), (8b)
J(w) = 1 2 w T w + 1 2 y 2 -H(x 2 ) 2 R -1 , (8c)
where according to (5), ( 6) and (8a)

x 2 (w) = M(x a 1 + A a 1 w 1:m ) + A q 2 w m+1:m+mq , (9)
mq is the size of the model noise ensemble A q 2 , and wn 1:n2 denotes a subvector of w formed by elements from n 1 to n 2 .

Condition of zero gradient of the cost function (8c) yields

w -(HA) T R -1 [y 2 -H(x 2 )] = 0, (10)
where

A ≡ [MA a 1 , A q 2 ], (11)
H ≡ ∇H(x 2 ), (12)
M ≡ ∇M(x 1 ).

(13) Equation ( 10) can be solved iteratively by the Newton method:

w i+1 = w i -D i ∇J(w i ), (14)
where D i is the inverse Hessian of the cost function (8c), and hereafter index i denotes the value of the corresponding variable at the ith iteration. We ignore the second-order derivatives in calculating the Hessian, which corresponds to employing the Gauss-Newton minimisation, so that

D i ≈ I + (H i A i ) T R -1 H i A i -1 , (15)
and ( 14) becomes

w i+1 -w i = I + (H i A i ) T R -1 H i A i -1 × (H i A i ) T R -1 y 2 -H(x i
This equation, required to obtain the analysis state, is the core of the IEnKF-Q method.

The other two necessary elements of the IEnKF-Q are the computations of the smoothed and filtered ensemble anomalies

A s
1 and A a 2 . Knowledge of A s 1 is needed to reduce the ensemble spread at t 1 in accordance with the reduced uncertainty after assimilating observations at t 2 ; and A a 2 is needed to commence the next cycle.

To find the analysed ensemble anomalies at t 2 , we first define the perturbed states at t 1 and t 2 :

δx 1 = A a 1 δu, δx 2 = MA a 1 δu + A q 2 δv, (17)
in terms of the perturbed δu and δv. We note that in the linear case, D approximated by (15) represents the covariance in ensemble space:

E[w (w ) T ] = D , (18)
where E is the statistical expectation and index denotes the value of the corresponding variable after convergence, so that, using (17):

A a 2 (A a 2 ) T = E[δx 2 (δx 2 ) T ] = A E[w (w ) T ](A ) T = A D (A ) T , (19)
and

A a 2 = A (D ) 1/2 = A I + (H A ) T (R) -1 H A -1/2 , ( 20
)
where D 1/2 is the unique symmetric positive (semi-)definite square root of a positive (semi-)definite matrix D.

Similarly,

A s 1 (A s 1 ) T = E[δx 1 (δx 1 ) T ] = A a 1 E[u (u ) T ](A a 1 ) T ; (21) therefore A s 1 = A a 1 (D 1:m,1:m ) 1/2 , ( 22
)
where Dn 1:n2 ,m1:m2 denotes a submatrix of D formed by rows from n 1 to n 2 and columns from m 1 to m 2 . It can be verified using (32) that the smoothed error covariance P s 1 = A s 1 (A s 1 ) T matches the Kalman smoother solution (Rauch et al. 1965, eq. 3.31).

Equations ( 16), ( 20) and ( 22) constitute the backbone of the IEnKF-Q.

Decoupling of u and v in the case of linear observations

In the case of a linear observation operator H, it is possible to decouple the solution for u and v. This is shown below by 265 transforming ( 16) to an alternative form.

Re-writing ( 16) as

w i+1 = I + (H i A i ) T R -1 H i A i -1 × (H i A i ) T R -1 y 2 -H(x i 2 ) + H i A i w i (23)
270 and using the identity

I + B T R -1 B -1 B T R -1 = B T (BB T + R) -1 , ( 24
)
where R is positive definite, we get

w i+1 =(H i A i ) T (H i A i ) T H i A i + R -1 275 × y 2 -H(x i 2 ) + H i A i w i , (25)
or, decomposing w i and A i ,

                       u i+1 =(H i M i A a 1 ) T (H i M i A a 1 ) T H i M i A a 1 + R i u -1 × y 2 -H x i 2 + H i M i A a 1 u i + H i A q 2 v i , v i+1 =(H i A q 2 ) T (H i A q 2 ) T H i A q 2 + R i v -1 × y 2 -H x i 2 + H i M i A a 1 u i + H i A q 2 v i , (26) (27)
where

R i u ≡ H i A q 2 (H i A q 2 ) T + R, (28)
R i v ≡ H i M i A a 1 (H i M i A a 1 ) T + R.
(29) 280 c 2017 Royal Meteorological Society Prepared using qjrms4.cls

Focusing on the increments of the iterates, we equivalently obtain

                       u i+1 -u i =D i u (H i M i A a 1 ) T (R i u ) -1 × y 2 -H(x i 2 ) + H i A q 2 v i -u i , v i+1 -v i =D i v (H i A q 2 ) T (R i v ) -1 × y 2 -H(x i 2 ) + H i M i A a 1 u i -v i , (30) (31)
where

D i u ≡ I + (H i M i A a 1 ) T (R i u ) -1 H i M i A a 1 -1 , (32)
D i v ≡ I + (H i A q 2 ) T (R i v ) -1 H i A q 2 -1 . ( 33
) 285 It is straightforward to verify that D i u = D i 1:m,1:m , and D i v = D i m+1:m+mq,m+1:m+mq .
Equations ( 30) and ( 31) represent an alternative form of equation ( 16) that makes it easy to see the decoupling of u and v in the case of linear H. In this case

H i = H = Const and H(x i 2 ) = H[M(x i 1 ) + A q 2 v i ] = H • M(x i 1 ) + HA q 2 v i , ( 34
)
so that (30) becomes

u i+1 -u i =D i u (HM i A a 1 ) T (R i u ) -1 × y 2 -H • M(x a 1 + A a 1 u i ) -u i . ( 35
)
295

It follows from ( 35) that in the case of linear observations, u can be found by the IEnKF algorithm (i.e, assuming perfect model) with modified observation error (28): R i u = R + HQH T . After that, v can be found from (31):

v = (HA q 2 ) T (R v ) -1 y 2 -H(x 2 ) + HM A a 1 u , ( 36
)
300 which can further be simplified using

x 2 = M(x 1 ) + A q 2 v and u = (HM A a 1 ) T (R u ) -1 y 2 -H • M(x 1 ) (37)
obtained from ( 35), finally yielding the non-iterative estimator

v = (HA q 2 ) T (R + HQH T ) -1 y 2 -H • M(x 1 ) . ( 38
)
Computationally, the decoupling reduces the size of D i , i.e.

(m + mq) × (m + mq), to that of D i u , i.e. m × m; however, it involves the inversion of a p × p matrix R i u , where p is the number of observations, which in large-scale geophysical systems can be expected to be much larger than the ensemble sizes m and mq.

The decoupling of u and v can be analysed in terms of probability distributions. This allows one to understand it at a more fundamental level and to connect the IEnKF-Q to the particle filter with optimal proposal importance sampling [START_REF] Doucet | On sequential Monte Carlo sampling methods for Bayesian filtering[END_REF], which is an elegant particle filter solution of our original problem with applications to the data assimilation in geosciences [START_REF] Bocquet | Beyond Gaussian statistical modeling in geophysical data assimilation[END_REF][START_REF] Snyder | Performance bounds for particle filters using the optimal proposal[END_REF][START_REF] Slivinski | Exploring practical estimates of the ensemble size necessary for particle filters[END_REF].

The posterior probability density function (PDF) of the analysis

p(x 1 , x 2 |y 2 ) is related to the IEnKF-Q cost function (2b) through J(x 1 , x 2 ) = -ln p(x 1 , x 2 |y 2 ). (39)
In all generality, the posterior PDF can be decomposed into

p(x 1 , x 2 |y 2 ) = p(x 2 |x 1 , y 2 )p(x 1 |y 2 ). (40)
It turns out that when H is linear both factors of this product have an analytic expression. This simplification is leveraged over when defining a particle filter with an optimal importance sampling [START_REF] Doucet | On sequential Monte Carlo sampling methods for Bayesian filtering[END_REF]. For our problem, one can show after some elementary but tedious matrix algebra that

-ln p(x 1 |y 2 ) = 1 2 x 1 -x a 1 2 (P a 1 ) -1 + 1 2 y 2 -H • M(x 1 ) 2 (R+HQH T ) -1 + c 1 , (41)
and

-ln p(x 2 |x 1 , y 2 ) = 1 2 x 2 -M(x 1 ) -QH T (R + HQH T ) -1 × [y 2 -H • M(x 1 )] 2 Q -1 +H T R -1 H + c 2 , (42)
x 2 =M(x 1 ) + QH T R + HQH T -1 × y 2 -H • M(x 1 ) . ( 43
)
It is simple to check that this expression, albeit written in ensemble space, is consistent with (38).

This decomposition explains at a fundamental level why the computation of the MAP of the IEnKF-Q were to decouple in (35,38) when the observation operator H is linear. However, this decoupling, valid for the MAP, does not immediately convey to the computation of the posterior perturbations.

The base algorithm

In this section, we put up an IEnKF-Q algorithm based on equations ( 16), ( 20) and ( 22). We refer to it as the base algorithm, because there are many possible variations of the algorithm based on different representations of these equations, including using decoupling of u and v in the case of linear observations described in section 3. Further, we do not include localisation, which is a necessary attribute of large-scale systems. The localisation of the IEnKF and IEnKS has been explored in [START_REF] Bocquet | Localization and the iterative ensemble Kalman smoother[END_REF]. In this paper, an implementation based on the local analysis method [START_REF] Evensen | The Ensemble Kalman Filter: theoretical formulation and practical implementation[END_REF][START_REF] Sakov | Relation between two common localisation methods for the EnKF[END_REF] has been proposed and may require the use of a surrogate model, typically advection by the fluid, to propagate a dynamically covariant localisation over long data assimilation windows. Such an implementation is actually rather straightforward for the IEnKF-Q since it is already formulated in ensemble space. Even though this is not the focus of this study, we will make preliminary tests of a local variant of the IEnKF-Q at the end of section 5 and provide its algorithm in Appendix A.

While the EnKF is a derivative-less method, it is possible to vary the type of approximations of Jacobians M and H with the ensemble used in the algorithm. In various types of the EnKF, it is common to use approximations of various products of H and M using ensemble of finite spread set based on statistical estimation 380 (3b) for sample covariance:

Hx ← H(E) 1/m, (44a)
HA ← H(E)(I -11 T /m), (44b)
Mx ← M(E) 1/m, ( 44c
) MA ← M(E)(I -11 T /m), ( 44d
) 385 HMA ← H • M(E)(I -11 T /m). (44e)
However, as pointed in [START_REF] Sakov | An iterative EnKF for strongly nonlinear systems[END_REF], it is also possible to use finite difference approximations:

Hx ← H(x1 T + εA) 1/m, ( 45a
) 390 HA ← H(x1 T + εA)(I -11 T /m)/ε, (45b)
Mx ← M(x1 T + εA) 1/m, (45c)
MA ← M(x1 T + εA)(I -11 T /m)/ε, (45d)
HMA ← H • M(x1 T + εA)(I -11 T /m)/ε, (45e)
where ε 1. Using these approximations results in methods of derivative-less state-space extended Kalman filter (EKF) type.

The difference in employing approximations ( 44) and ( 45 It was found later [START_REF] Bocquet | Combining inflation-free and iterative ensemble Kalman filters for strongly nonlinear systems[END_REF]) that performance of schemes based on finite difference approximations can be improved by conducting a final propagation with a finite 410 spread ensemble. The corresponding schemes were referred to as "bundle" variants, while the schemes using finite spread approximations -as "transform" variants. The algorithm 1 is a transform variant of the IEnKF-Q method.

Algorithm 1 A "transform" variant of the IEnKF-Q. The pieces of pseudo-code highlighted in red show changes relative to the IEnKF algorithm in absence of model error. "SR(A, m)" denotes ensemble size reduction from m + mq to m.

1: function [E 2 ] = ienkf cycle(E a 1 , A q 2 , y 2 , R, M, H) 2:
x a 1 = E a 1 1/m 3:

A a 1 = (E a 1 -x a 1 1 T )/ √ m -1 4: D = I, w = 0 5: repeat 6: x 1 = x a 1 + A a 1 w 1:m 7: T = (D 1:m,1:m ) 1/2
8:

E 1 = x 1 1 T + A a 1 T √ m -1 9: E 2 = M(E 1 )
10:

HA 2 = H(E 2 )(I -11 T / m) T -1 / √ m -1 11: HA q 2 = H(E 2 11 T / m + A q 2 m q -1) ×(I -11 T / m q )/ m q -1 12: HA = [HA 2 , HA q 2 ]
13:

x 2 = E 2 1/m+A q 2 w m+1:m+mq 14: until ∆w < ε 19:

∇J = w -(HA) T R -1 [y 2 -H(x 2 )]
A 2 = E 2 (I -11 T / m) T -1 20: A = [A 2 / √ m -1, A q 2 ] D 1/2 21: A 2 = SR(A, m) √ m - 1
22: 2 , before applying the nonlinear map from the ensemble space to the observation space:

E 2 = x 2 1 T + A 2 23
w → H M(x a 1 + A a 1 w 1:m ) + A q 2 w m+1:m+mq .

(

) 46
The implementation of this joint transform is less simple than that offered by the one in line 7, which merely amounts to using the smoothing anomalies marginalised at t 1 . Another simple possibility is to choose T = [D 1/2 ] 1:m,1:m , which remains a positive definite matrix. We have checked that these three approaches yield the same quantitative results for all the experiments reported below, except for those on localisation.

However, we expect that the optimal joint transform mentioned above could make a difference in the presence of a significantly nonlinear observation operator (not tested).

Because the IEnKF-Q uses augmented ensemble anomalies A possible way of reducing the ensemble size to m is to keep the m -1 largest principal components of A a 2 and use the remaining degree of freedom to centre the reduced ensemble to zero. In practice the magnitude of ensemble members produced by this procedure can be quite non-uniform, similar to that of the SVD spectrum of the ensemble. This can have a detrimental effect on performance in a nonlinear system; therefore, one may need to apply random mean-preserving rotations to the ensemble to render the ensemble distribution more Gaussian. * The augmentation of the propagated state error anomalies and model error anomalies has also been used in the reduced rank square root filter by Verlaan and Heemink (1997, eq. ( 28)).

c 2017 Royal Meteorological Society Prepared using qjrms4.cls

This reduction, based on the SVD, actually represents a marginalisation, in a probabilistic sense, over all the remaining degrees of freedom. Assuming Gaussian statistics of the perturbations, the marginalisation can be rigorously performed this way as the excluded modes are orthogonal to the posterior ensemble subspace. In the limit of the Gaussian approximation, this guarantees that the reduction to the posterior ensemble space accounts for all information available in this ensemble space.

In the IEnKF-Q algorithm, the computational cost induced by mq is due to the cost of the observation operator to be applied to the mq additional members of the ensemble in the analysis, as seen in ( 10), ( 11) and ( 12). In contrast, the cost associated to the m members of the ensemble is due to the application of both the evolution model and observation operators, which is potentially much greater, as seen in ( 10), ( 11), ( 12) and ( 13). A large mq also potentially increases the computational cost of the nonlinear optimisation of cost function (8c), which is nonetheless expected to be often marginal compared to the computational cost of the models. For realistic applications, mq could be chosen to be reasonably small by pointing A q 2 to the most uncertain, possibly known a priori, directions of the model, such as the forcings.

These are often called stochastic perturbations of the physical tendencies, see [START_REF] Buizza | Stochastic representation of model uncertainties in the ECMWF ensemble prediction system[END_REF], section 2.5 of Wu et al.

(2008) and section 5.c of [START_REF] Houtekamer | Model error representation in an operational ensemble Kalman filter[END_REF]. Because they are randomly selected, these perturbations are actually meant to explore a number of independent model error directions greater than mq but over several cycles of the DA scheme.

Numerical tests

This section describes a number of numerical tests of the IEnKF-Q with the Lorenz-96 model [START_REF] Lorenz | Optimal sites for suplementary weather observations: simulation with a small model[END_REF] to verify its performance against the IEnKF and EnKF.

The model is based on 40 coupled ordinary differential equations in a periodic domain:

ẋi = (x i+1 -x i-2 ) x i-1 -x i + 8, i = 1, . . . , 40;
(47)

x 0 = x 40 , x -1 = x 39 , x 41 = x 1 . (48)
Following Lorenz and Emanuel (1998), this system is integrated with the fourth order Runge-Kutta scheme, using a fixed time step of δt = 0.05, which is considered to be one model step.

The model noise is added after integrating equations ( 47 averaged over 10 5 cycles after a spinup of 5000 cycles. For each run the optimal inflation is chosen out of the following set: 520 {1, 1.02, 1.05, 1.1, 1.15, 1.2, 1.25, 1.3, 1.4, 1.5, 1.75, 2, 2.5, 3, 4}. In the following experiments, we choose mq = 41 for the IEnKF-Q so that A q 2 can span the whole range of Q whatever its actual form, which ensures that (6b) is exactly satisfied. This should highlight the full potential of the IEnKF-Q.

525

As justified in section 4, random mean-preserving rotations of the ensemble anomalies are sometimes applied to the IEnKF-Q, typically in the very weak model error regime.

The performance of the IEnKF-Q is compared to that of the EnKF using the ensemble transform Kalman filter scheme (ETKF, 530 unit of time is kept constant. Even though this value of Q is two orders of magnitude smaller than R when T = 1, we found it to be realistic. Indeed, the standard deviation of the perturbation added to the truth of the synthetic experiment is 0.1 in that case, to be 615 compared to a root mean square error of about 0.2 obtained for the analysis of the EnKF in a perfect model experiment with T = 1,

the Lorenz-96 model, and the data assimilation setup subsequently described. the increase in the model error results in a more pronounced advantage of the IEnKF-Q over the other methods. The relative performance of the non-iterative schemes is better than in Figure 1 because of the increased ensemble size.

Test 2: model noise magnitude

This test investigates the relative performance of the schemes depending on the magnitude of model error both in a weakly nonlinear and significantly nonlinear case. The tests for all schemes involved are conducted with ensemble size m = 20.

Moreover, results with full-rank ensemble m = 41 are shown for IEnKF-Det and IEnKF-Q.

Figure 3 shows the performance of the schemes in a weakly nonlinear case T = 1. For small model error q 3 • 10 -3 , all schemes perform similarly; for q 3 • 10 -3 the IEnKF-Q starts to outperform other schemes; and from q 5 • 10 -2 the iterative schemes IEnKF-Rand and IEnKF-Det start to outperform their non-iterative counterparts EnKF-Rand and EnKF-Det.

Interestingly, the IEnKF-Q does not show advantage over IEnKF-Det when both use full-rank ensembles m = 41 (except, perhaps, some very marginal advantage for larger model error q 0.1). At a heuristic level, both schemes indeed explore the same model error directions. At a mathematical level, there are indications in favour of this behaviour, including the marginal advantage. Thanks to the decoupling analysis, we see that when m = 41 the smoothing analysis at t 1 of the IEnKF-Q and that of the IEnKF-Det become equivalent. However, the filtering analysis c 2017 Royal Meteorological Society Prepared using qjrms4.cls at t 2 of the IEnKF-Q, (43), is different albeit close to that of the IEnKF-Det, which is just the forecast of x 1 , i.e., the first term of (43). Concerning the update of the anomalies of the ensemble, we have seen that ( 52) is the P a 2 of the IEnKF-Det when m = 41, and is very close to the expression of P a 2 for the IEnKF-Q, (53), except maybe when Q is large.

10 -4 10 -3 10 -2 10 -1 q 0.17 In the significantly nonlinear case (T = 10, Figure 4), the noniterative schemes are no longer able to constrain the model. The empirically modified iterative schemes IEnKF-Rand and IEnKF-Det yield performance similar to the IEnKF-Q up to q 2 • 10 -3

with the ensemble size m = 20, and up to q 10 -2 with the ensemble size m = 41; however, apart from underperforming the IEnKF-Q for larger model errors, they also lose stability and are unable to complete 10 5 cycles necessary for completion of these runs. Interestingly, for very large model error q = 0.5 (Q = 5I)

the IEnKF-Q yields similar performance with the ensemble size of m = 20 and m = 41, with the RMSE ∼ 0.94 much smaller than the average magnitude of the model error ∼ 2.2. In this regime, the analysis RMSE remains smaller than that obtained with the sole observations (estimated to be 0.994 1 by A. Farchi, personal communication). This could be due to the variational analysis which spans the full state space since mq = 41, and, in this regime, little depends on the prior perturbations A a 1 .

Test 3: ensemble size

This test investigates the performance of the methods depending on the ensemble size both in a weakly nonlinear (T = 1, Figure 5) In line with the results of Test 2, we observe that the non-iterative schemes do not perform well in the significantly nonlinear case. The IEnKF-Q outperforms the other schemes when using smaller ensembles, but yields a performance similar to that of the IEnKF-Det with a full-rank (or nearly full-rank) ensemble. This is mainly due to its search for the optimal analysis state over a large subspace. Likewise, the performance of the IEnKF-Q degrades when restricting the model error directions to that of the ensemble space (yielding mq = m), and yet remains slightly better than the IEnKF-Det (not shown).

Test 4: localisation

A couple of numerical experiments are carried out to check that the IEnKF-Q can be made local. However, a detailed discussion of the results is out of scope, since our primary concern is only to confirm the feasibility of a local IEnKF-Q. To this end, we have merged the local analysis as described in [START_REF] Bocquet | Localization and the iterative ensemble Kalman smoother[END_REF] 10 -4 10 -3 10 -2 10 -1 q 0.15 0.17 Both tests show that a local IEnKF-Q is not only feasible but also yields very accurate results, with a local 10-735 member implementation outperforming a global 20-member implementation.

Discussion

The presence of model error in a DA system causes lossy transmission of information in time. Therefore, there may be potential for empirical use of the additive model error in the EnKF to improve numerics. It indeed can often be perfectly feasible to specify some sort of additive model error as a tuneable parameter of a suboptimal system, similarly to the common use of inflation. In fact, a number of studies found that using empirical additive model error in EnKF systems, alone or in combination with inflation, can yield better performance than using inflation only (e.g., [START_REF] Whitaker | Ensemble data assimilation with the NCEP global forecast system[END_REF]. Another example of employing model error in a suboptimal system is using (so far without marked success) the hybrid covariance factorised by 790 ensemble anomalies augmenting a small (rank deficient) dynamic ensemble and a large static ensemble [START_REF] Counillon | Application of a hybrid EnKF-OI to ocean forecasting[END_REF].

In this study, we have assumed that model error statistics are known. In some simple situations, these could be estimated online with techniques such as those developed by Todling 795 (2015). Nonetheless, using these empirical Bayesian estimation techniques here would have obscured the methodological introduction to the IEnKF-Q.

Summary

This study proposes a new method called IEnKF-Q that extends 800 the iterative ensemble Kalman filter (IEnKF) to the case of additive model error. The method consists of a Gauss-Newton minimisation of the nonlinear cost function conducted in ensemble space spanned by the propagated ensemble anomalies and anomalies of the model error ensemble. To simplify the 805 algebraic form of the Gauss-Newton minimisation, the IEnKF-Q concatenates the expansion coefficients u and v into a single vector w, and augments ensemble anomalies MA a 1 and A q 2 into a single ensemble A. After that, the minimisation takes the form (16) similar to that in the perfect-model case.

810

Algorithmically, the method can take many variations including "transform" and "bundle" versions, and various localisation approaches. An example algorithm suitable for low dimensional systems is presented in section 4. Using this algorithm, the method is tested in section 5 in a number of experiments with the Lorenz-815 96 model. In all experiments, the IEnKF-Q outperforms both the EnKF and IEnKF adapted for handling model error either in a "stochastic" or "deterministic" way, except in situations with full-rank ensemble and weak to moderate model error, where it performs equally with the IEnKF-Det. Surprisingly, it also 820 outperforms these methods in weakly nonlinear situations, when the solution is essentially found at the very first iteration, and iterative schemes should not have any marked advantage over noniterative schemes. This is caused by using full-rank (augmented)

forecast ensemble anomalies in the analysis, and only reducing the ensemble size back to the initial one at the very end of the cycle.

Note that in practice the cost of the IEnKF-Q in high-dimensional systems can be expected to be similar to that of the IEnKF because both methods use ensembles of the same size in propagation.

One interesting feature of the IEnKF-Q is the decoupling of iterations over u and v made possible in presence of a linear observation operator. In this case, u can be found using the (perfect-model) IEnKF with an increased observation error covariance, followed by obtaining v in a single iteration. The decoupling can be the underlying reason why in certain situations the IEnKF-Q and IEnKF-Det show equal performance.

Algorithm 2 A local analysis and transform variant of the IEnKF-Q, using a Gauss-Newton minimisation. "SR(A, m)" denotes ensemble size reduction from m + m q to m. x a 1 = E a 1 1/m 3:

A a 1 = (E a 1 -x a 1 1 T )/ 11:

[E 1 ] i = [x 1 1 T ] i + [A a 1 T i ] i √ m -1 12:
end for

13:

E 2 = M(E 1 )

14:

for i = 1, . . . , n do 15:

[x 2 ] i = [E 2 1/m] i + [A q i w i m+1:m+mq ] i 16:
end for 17:

for i = 1, . . . , n do 18:

HA i 2 = H(E 2 )(I -11 T / m) (T i ) -1 / √ m -1 19:
HA q i = H(E 2 11 T / m + A q i mq -1) ×(I -11 T / mq)/ mq -1 20:

HA i = [HA i 2 , HA q i ]
21:

∇J i = w i -(HA i ) T (R i ) -1 [y 2 -H(x 2 )]
22:

D i = [I + (HA i ) T (R i ) -1 HA i ] -1
23:

∆w i = -D i ∇J i
24:

w i := w i + ∆w i 25:

end for 26:

until n i=1 ∆w i < nε 27:

for i = 1, . . . , n do 28:

A i 2 = E 2 (I -11 T / m) (T i ) -1 29:

A i = [A i 2 / √ m -1, A q i ] (D i ) 1/2
30:

A i 2 = SR(A i , m) √ m -1

  115

  schemes.

  : end function Line 6 of the algorithm corresponds to (5a); line 7 calculates the 415 ensemble transform in (22); multiplication by √ m -1 on line 8 restores normalisation of ensemble anomalies before propagation to statistically correct magnitude; division by √ m -1 on line 10 changes it to the algebraically convenient form P = AA T used in (5). The observation ensemble anomalies of the model 420 noise ensemble H i A q 2 are calculated on line 11. This involves adding ensemble mean and re-normalisation before applying the observation operator. In the case of linear observations this line would reduce to HA q 2 = H(A q 2 ). Line 13 corresponds to (6a), and line 20 to (20). 425 The ensemble transform applied in line 7 is actually a bit restrictive, though it is sleek and convenient. Its potential suboptimality is obvious in that the transform correctly applies to the evolution model propagation of the ensemble, but not to the observation operator. In this context, faithfully enforcing the transform principle proposed in Sakov et al. (2012) would imply applying (on the right) the transform matrix T = D 1/2 to the joint anomaly matrix A a 1 , A q

( 11 )

 11 * , it increases the ensemble size from m to m + mq. Consequently, to return to the original ensemble size one needs to conduct ensemble size reduction at the end of the cycle. If the ensemble size is equal to or exceeds the dimension of the model subspace, such a reduction can be done losslessly; otherwise it is lossy. The reduction of the ensemble size is conducted on line 21 of Algorithm 1. Multiplication by √ m -1 performs renormalisation of A 2 back to the standard EnKF form (3b).

  ) for the time length of each DA cycle. An alternative would be to gradually add model noise at each model time step, which would 505 be consistent if the original model was based on a continuous stochastic differential equation. Even though less elegant, we chose the former approach because, in that case, the actual model error covariance is guaranteed to match the assumed model error covariance. However, an implication of such approach is that the 510 properties of the resulting stochastic model depend on the length of the cycle. This applies to the true model state as well as to the ensemble members of the EnKF-and IEnKF-based methods to be defined later.In the following twin experiments, each variable of the model 515 is independently observed once per DA cycle with Gaussian observation error of variance 1: R = I. The performance metric we use is the filtering analysis root mean square error (RMSE)

Figure 1 .

 1 Figure 1. Test 1: dependence of the mean analysis RMSE on the time interval between observations T in units of δt. m = 20, Q = 0.01T I.

Figure 1 Figure 2 Figure 2 .

 122 Figure 1 compares the performance of the EnKF, the IEnKF

Figure 3 .

 3 Figure 3. Test 2: dependence of the mean analysis RMSE on the magnitude of model error in a weakly nonlinear case T = 1; Q = qT I, m = 20.

Figure 4 .Figure 5 .Figure 6 .

 456 Figure 4. Test 2: dependence of the mean analysis RMSE on the magnitude of model error in a significantly nonlinear case T = 10; Q = qT I, m = 20.

Figure 7 .

 7 Figure 3, we transfer the RMSE curves of the global IEnKF-Q with ensemble size m = 20 and m = 41. In addition, we consider the local IEnKF-Q with an ensemble size of m = 10, which requires the use of localisation. Again, the localisation length is 10 grid points. The RMSEs are plotted in Figure 8. The local IEnKF-Q RMSE curve lies in between those for the global IEnKF-Q with

Figure 8 .

 8 Figure 8. Test 4: dependence of the mean analysis RMSE on the magnitude of model error in a weakly nonlinear case T = 1; Q = qT I.

  ([D i ] 1:m,1:m ) 1/2

  P f 1 -initial state error covariance, y i -observations, H i -observation operator, R i -observation error covariance, M i -model operator, and Q i -model error covariance; and the norm notation x 2

time, K -number of cycles plus one, x i -(model) state at cycle i, x i -state estimate (analysis) at cycle i, x f 1 -initial state estimate, B ≡ x T Bx is used. The cost function is assumed to be, generally, nonlinear due to nonlinear operators M and H.

  t 1 and t 2 ; x a 1 and

	model states E. These two representations are related as follows:
		P a 1 are the filtering analysis and filtering state error covariance at
		time t 1 , which have been obtained in the previous cycle; all other
		variables have direct analogues in formulation (1) of the global
		problem. The function (2b) should be seen as the state-space cost
		function associated with the analysis of the IEnKF-Q. It is the key
		to the method's derivation.
		The main difference between the KF and the EnKF is their
		representation of the state of the DA system (SDAS). In the
		KF, the SDAS is carried by the state estimate x and state error
		covariance P. In the EnKF, the SDAS is carried by an ensemble of
	c 2017 Royal Meteorological Society	Prepared using qjrms4.cls

  c 1 and c 2 are constants that neither depend on x 1 norx 2 . Note that p(x 1 |y 2 ) is non-Gaussian while p(x 2 |x 1 , y 2 ) is a Gaussian PDF thanks to the linearity of H.This decomposition enables to minimise J(x 1 , x 2 ) in two steps. First, one can minimise -ln p(x 1 |y 2 ) over x 1 yielding the maximum a posteriori (MAP) solution x 1 . This identifies with the smoothing analysis of a perfect model IEnKF but with R replaced with R + HQH T . Second, the MAP solution of the minimisation of -ln p(x 2 |x 1 , y 2 ) is directly given by

c 2017 Royal Meteorological Society Prepared using qjrms4.cls where

1 :

 1 function [E 2 ] = lienkf cycle(E a 1 , {A q i } i=1,...,n , y 2 , {R i } i=1,...,n , M, H)

	2:

c 2017 Royal Meteorological Society Prepared using qjrms4.cls

) -w i . (16) c 2017 Royal Meteorological Society Prepared using qjrms4.cls

Acknowledgements

The authors would like to thank two anonymous reviewers for their useful comments and suggestions. They are grateful to T. Janjić for her invitation to submit this work. M. Bocquet and J.-M. Haussaire acknowledge the contribution of INSU via the LEFE/MANU DAVE project. CEREA is a member of the Institut Pierre Simon Laplace (IPSL).

P. Sakov, J.-M. Haussaire and M. Bocquet

The stochastic approach is

where Ξ is an n × m matrix whose columns are independently sampled from N (0, I). With these anomalies, one has

When applying this approach 545 to the EnKF, we refer to it as EnKF-Rand. Be wary that the EnKF-Rand is not the original stochastic EnKF; its analysis step is deterministic.

The deterministic approach is to substitute the full covariance matrix Q with its projection Q onto the ensemble subspace: Q = 550 Π A QΠ A , where Π A = AA † is the projector onto the subspace generated by the anomalies (the columns of) A = MA a 1 , and

Hence, the anomalies that satisfy (50) are [START_REF] Raanes | Extending the square root method to account for additive forecast noise in ensemble methods[END_REF])

560

When applying this approach to the EnKF, we refer to it as EnKF-Det.

Those two simple ways to add model noise to the analysis of the EnKF can also be applied to the standard IEnKF. At each iteration, the IEnKF smoothing analysis at t 1 , yielding A s 1 , is followed by 565 either ( 49) or (51) with A a 1 replaced with A s 1 , which yields the IEnKF-Rand and IEnKF-Det methods, respectively.

These heuristic methods are fostered by the decoupling analysis in section 3 when H is linear. This analysis suggests that the IEnKF smoothing analysis at t 1 should actually be performed 570 with an observation error covariance matrix R + HQH T in place of R. Note, however, that we did not observe any significant differences in the performance of the IEnKF-Rand and IEnKF-Det with or without this correction.

Further, it can be shown that the heuristic IEnKF-Rand yields

whereas the rigorous IEnKF-Q yields

These posterior error covariance matrices ( 52) and ( 53

Test 1: nonlinearity

This test investigates the performance of the schemes depending on the time interval between observations, covering DA regimes from weakly nonlinear to significantly nonlinear. The ensemble size is m = 20, chosen so that it is greater than the dimension of the unstable-neutral subspace (which is here 14 and to which we where the localisation length is defined to be their c parameter.