Khalid Addi

Samir Adly

Daniel Goeleven
email: daniel.goeleven@univ-reunion.fr

Michel Thera
email: michel.thera@unilim.fr

A NONSYMMETRIC LINEAR COMPLEMENTARITY PROBLEM TO SOLVE A QUASISTATIC. ROLLING FRICTION AL CON TACT PROBLEM

Keywords: Complementarity Problems, Frictional contact problems, Rolling con tact, Boundary Element Method, Mathematical programming

A simple approach and an algorithm are proposed to solve the qua sistatic rolling frictional contact problem between an elastic cylinder and a fiat rigid body. The discretization is based on the boundary element method. The unilateral frictional contact problem (nonsmooth but monotone) is formulated in a compact form as a nonsymmetric linear complementarity problem which is solved using Lemke's algorithm.

INTRODUCTION

In the quasistatic case of a rolling elastic cylinder in frictional contact with a flat rigid body, the relative slip velocity is a very important parameter of the problem� To compute it, a direct method is presented in Abascal et al. [6, 7, 8J, where two approaches are used considering the velocity as unknown variable in the LCP solving the frictional contact.

Here, we propose• another approach based in solving the static frictional contact problem and, afterwards, computing the velocity using an explicit appropriate re lation. The boundary element method (BEM) is used to discretize the mechanical problem. Then the elastostatic; frictional contact, rolling problem is formulated as a compact LCP . using the boundary element method and the rigid displacement approach.

MODELLING, DISCRETIZATION AND CONDENSATION OF THE PROBLEM.

Following the approach used in [START_REF] Addi | On solving a rolling fr ictional contact problem using BEM and Mathematical Programming[END_REF] the rolling problem is reduced to a plane strain state (see, e.g., [START_REF] Gonzalez | An algorithm to solve coupled 2D rolling contact problems[END_REF][START_REF] Gonzalez | An algorithm to solve 2D transient rolling contact problems[END_REF][START_REF] Gonzalez | Solving 2D rolling problems using the NORM-TANG iteration and mathematical programming[END_REF]) where the Coulomb's law describes the friction.

The friction is assumed to follow the dry Coulomb's law where normal and tan gential tractions on the boundaries of the contact zone are related via a coefficient of friction. Under this assumption, two points belonging to the •cylinder A and _the• rigid body B fall in three different states relative to each other (Fig.

( 1) Stick state

(1) r� < 0, r� < 0, (2)
(2) Slip state

� �n= O , r O j � n >O , r � >=O • �+-Stick/Slip�• . � X Sepa r ation Separation FIGURE 1. Potential contact zone.
r� < 0, r� < 0, r� = r�, ! r �l = �� r �l, l r ¥1 == �� r� l , On= 0, sgn ( st) = -s g n ( r #),

(3) Separation state

(3)
rJt == 0 , r� == 0, r: = 0, r!J. == 0, On > 0. Here, 1-L is the constant coefficient of friction, r�, r�, r� and r !J. are the normal and tangential tractions of cylinder A and rigid body B, respectively, at contacting points, U N and ur are the normal and the tangential cylinder surface displacements, St is the slip velocity and On is the normal separation given by: (4)

On== Ono-UN.

The initial separation, Ono, for a cylinder and a flat body in contact can be approx imated by [START_REF] Goeleven | Solvability Theory and Projection Methods for a Class of Singular Variational Inequalities. Elastostatic Unilateral Contact Applications[END_REF] where R is the cylinder radius and x is the Eulerian coordinate along the contact zone used to position each pair of points at each time relative to a rigid body position of the cylinder.

When applying an Eulerian description of particles moving through the contact area, the relative tangential slip velocity of each cylinder surface point is defined as

(6)
St == b t = d O t( X, T ) dr
where T is the time coordinate, X = x (r) is the Cartesian coordinate of each point relative to fixed axes and varying time r, and Ot is the tangential separation given by [START_REF] Gonzalez | An algorithm to solve 2D transient rolling contact problems[END_REF] . Ot == ( X A -X B ) + U� + u ¥ where xA and x8 are the cartesian coordinates of contacting points of the cylinder and the rigid respectively. u � and u¥ describe the tangential displacements of at contacting points of the cylinder and the rigid body, respectively. Substituting ( 7) into (6) leads to (8 )

_ VA V B VA A V B B A B 8 t - - + u T, x + UT,x + UT,r + U T,r
where va is the velocity of a point of the rigid body a (a == A, B) in the x-direction.

u � , x = 8;[. Under steady-state rolling conditions, the variation with time vanishes, so.

(

) S t = VA -V B + VA uA + V3 u 3 . �X �X 9
which . is usually approximated [START_REF] Kalker | Three dimenensional elastic bodies in rolling contact[END_REF] by

(10) Bt = l V I ( e � sgn( V ) (u � , x + u ¥, x ) )
where V is given by (V A + v B )/2, and e is the normalized relative rigid slip ve-vA-v B .

locity (creepage) , defined as e == I V [

. In fact, we consider St, in practice, as dimensionless variable denoted by s; [START_REF] Gonzalez | Solving 2D rolling problems using the NORM-TANG iteration and mathematical programming[END_REF]:

(11) s; == � + sgn( V )(u � ,x + u ¥, x )
When the displacement derivatives are approximated using a finite difference scheme, the tangential slip velocity for a point located at the coordinate Xi is, under steady state rolling conditions, expressed as

(12)
where hi denotes the distance between two adjacent boundary points Xi+l and Xi• To analyze coupled 2D rolling contact between two cylinders, Abascal et al. have formulated the problem considering the slip velocity as unknown variable. Indeed, in [START_REF] Gonzalez | Solving 2D rolling problems using the NORM-TANG iteration and mathematical programming[END_REF], the NORM-TANG iteration [START_REF] Panagiotopoulos | Hemivariational inequalities. Applications in mechanics and engineering[END_REF] to solve two complementarity problems. In [START_REF] Gonzalez | An algorithm to solve coupled 2D rolling contact problems[END_REF], the problem is analyzed by minimizing a function representing the equilibrium equation and the contact restrictions. In our case, the static contact-friction problem between an elastic cylinder and a flat rigid body is first solved and, subsequently, the slip velocity is computed explicitly using the above formula. The model is discretized using the boundary element method which is more suitable for this class . of problem where the nonlinearity is limited to the boundary of the body. The here adopted formulation follows the lines of [START_REF] Stavroulakis | Hard and Soft fingered robot grippers[END_REF], [START_REF] Antes | The boundary integral approach to static and dynamic contact problems[END_REF], [START_REF] Klarbring | A mathematical programming approach to three-dimensional contact problems with friction[END_REF], [12], [START_REF] Stavroulakis | Nonlinear equation approach for inequality elastostatics: a two dimensional BEM implementation[END_REF] and [START_REF] Stavroulakis | On the rigid body displacements and rotations in unilateral contact problems and applications[END_REF] where the matrix formulation of displacement boundary element in elastostatics is given by [START_REF] Brebbia | Boundary elements. An introductory course[END_REF] : [START_REF] Klarbring | A mathematical programming approach to contact problems with friction and varying contact surface[END_REF] Hu == Gt.

Here, u is the vector of the nodal boundary displacements, t is the vector of element boundary tractions, and H; G are the appropriate matrices.

The number of equations in (13) depends on the number of nodes on the dis-. cretized boundary. Let us note that the number of boundary element tractions (i.e., the dimension oft) depends on the nature of the boundary elements used. • The classical approach for the solution of the bilaterally constrained structures goes through the specification of appropriate, known boundary displacements or tractions, the rearrangement of the system (13), and finally, the formulation of a nonsymmetric system of equations, after reordering:

(14) Ay = b.
The 2n-dimensional vector y contains all the unknown boundary displacements or tractions of the problem. In the contact area, both displacements and tractions are unknown. They must be kept in the formulation and connected with the inequality and complementarity relations of the unilateral contact mechanism. After solving the arising LCP, one knows which of these . variables. vanishes. Thus, we proceed with condensation and, then, formulate the linear complementarity problem.

Let us consider Ue and te as being the boundary nodal displacements and trac tions, respectively, at the frictional unilateral contact boundary of the cylinder. After partitioning the boundary of the cylinder, the equation [START_REF] Kwack | A complementarity problem formulation for two dimensional frictional contact problems[END_REF] gives: [START_REF] Lemke | Some pivot scheme for the linear complementarity problem[END_REF] When ne is the number of the nodes at the unilateral contact boundary, then Ue and te have 2ne elements, Hee and Gee have 2ne x 2ne elements, Hff has 2(n-ne) X 2(n-ne) elements, Gje, Hje and Hef have 2(n-ne) X 2nc elements. u denotes the nodal boundary displacements of the cylinder outside the contact area. • The next step is to perform a local coordinate transformation so that normal and tangential to the unilateral boundary quantities appear in the formulation. There fore, let us consider w and r as the natural local coordinates (normal and tangential coordinates) of the displacements and. the tractions in the contact boundary, re spectively:

(16) u� = [�t], t�= [�t], w i = [��], r i = [��].
The transformation for � single unilateral boundary node i reads: [START_REF] Stavroulakis | Nonlinear equation approach for inequality elastostatics: a two dimensional BEM implementation[END_REF] with

C . i-i tUc -W, (18)
e i = [ CO� tP i sin tP i ] .

-Sill tP i COS </Ji

Since e i -1 = e[, relations (17) can be inverted:

i T i T (19) U c == c i w , tc = -C i r.
Taking into account these transformations, one has (20)

Then, the problem will be condensed, i.e., in order to get the flexibility matrix, the unknown variable u will be reduced by considering the relation (20): We follow the formulation of [START_REF] Klarbring | A mathematical programming approach to contact problems with friction and varying contact surface[END_REF] and [START_REF] Kwack | A complementarity problem formulation for two dimensional frictional contact problems[END_REF]. Let the normal forces and the friction forces be assembled in vectors rN = {rNl, ... , rNn} T and ry == {rr1, rr2, ... , rrn} T , respectively, w�ere, e.g., rri is the frictional force of the i-th contact node.

Coulomb's law of dry friction connects the tangential (frictional) forces with the normal (contact) forces by the relation

• (25) 'Yi == fl i T Nil -/rril , i == 1, . . . , n, li � 0.
Here J.J denotes the absolute value and JL is the friction coefficient. The friction mechanism is considered to work in the following way: If lrril < fliT Nil (i.e., 'Yi > 0),. the slipping value liT must be equal to zero, and if /rr i / = JL/rN i / (i.e., 'Yi � 0), then we have slipping in the opposite direction of rr i:

if "'/i > 0, then YTi = 0 (26) if 'Yi == 0,
.

[ 1 -1]

T � == 0 0 ' Finally, the slip value in (25), ( 26) is written as

(30) T \ .,>_0 YT == TA, /\.
where A is a vector of nonnegative slipping parameters. Then, 1 and .X fulfil the following complementarity condition:

(31) ')'TA== 0 FIGURE 2.
The mesh around the cylinder contour.

3.1. The Linear Complementarity Problem formulation. To formulate the linear complementarity problem the rigid body displacements approach is adopted. The slipping value ,\ and the tangential displacements ur are then related by the compatibility relation:

Tr.-\ur == dr where dr denotes the initial tangential distance.

The general decomposition scheme with slack variables can be found in more details for two-dimensional friction problems in [14}.

The structure is assumed linear elastic which, on the assumption that everything outside of the frictional contact interfaces has been condensed out, reads:

(33)
ii == Fr

•where

F == [FNN FNT] . FrN Frr ' ---_ [ rN ] r- .
rr

Here, F is the symmetric flexibility matrix where F N N is an n x n nonsingular matrix with the mechanical meaning of being the normal flexibility matrix, Frr is a 2n x 2n nonsingular matrix (the tangential flexibility) and F N T, FrN are the corresponding couple flexibility matrices. By using the previous relations, the unilateral kinema�ic conditions normal and tangential to the interface, take the form: Note also that (37) can be written as the variational inequality:

(39) z > 0 : (Mz + b) T ( v -z) == 0, V v � 0 4. NUMERICAL RESULTS
To discretize the rolling problem with BEM (Fig. 2), constant elements are used. In fact, Karami [START_REF] Karami | A boundary element method for two dimenensional contact problems[END_REF] and Stavroulakis et al. [START_REF] Stavroulakis | On the rigid body displacements and rotations in unilateral contact problems and applications[END_REF] remarked the existence of oscillations in the traction behaviour with quadratic elements. Young modulus and Poisson coefficient are taken to be 94500 and 0.1, respectively. The rigid body normal and tangential displacement used were dN = R/2 and dr == 0.1d N where R is the radius of the cylinder. The LCP (3 7 is solved using Lemke's algorithm [15) The results shown in the Figures 3,4,5, and 6 are obtained by applying a normal and a tangential rigid body displacements approach. Figure 3 shows the existence of a separation area in the potential contact zone. The stick/slip area is not centered because of the tangential effect. In the Figures 4 and5, the normal and the tangential tractions in the frictional contact case are plotted. Figure 6 shows the velocity behaviour only in the effective contact zone.

u[

  == Hj} f1 -Hj}G fee T r-Hj} H1cCT w Then, (22) [HccCT-HcJHj}HJcCT] w = [!c-HcJHj}fJ] -[ccccr-HcJH j }GJcCT] r Thus, w can be written as: Wn] == [Wno] + [Fnn Fnt] [r N] . . Wt ww Ftn Ftt rr 3. THE FRICTIONAL CONTACT FORMULATION

  then there exists a> 0 such that YTi = -aT T i• By assembling the contributions of all (n) unilateral nodes, relation (26) reads: (27) with the matrices Tr and T N (28) Tr == diag [T}, Tf, ... , Tr] , • T N == diag [T}v, TJv, . .• . , TN] . These matrices are obtained from the linearized friction law considered in 2D [2, 13].

FIGURE 3 .

 3 FIGURE 3. Distribution of stick/slip and normal separation in the potential contact zone.

FIGURE 4 .FIGURE 5 .

 45 FIGURE 4. Normal and tangential tractions in the frictional contact case (!-l == 0.1).

5 FIGURE 6 .

 56 FIGURE 6. Tangential velocity in the effective contact zone: Depen dence on the friction coefficient

Acknowledgements. The first author expresses his thanks to the Alexander von Humboldt-Stiftung, Bad Godesberg, Bonn, for the financial support and to Prof. H. Antes and Prof. G.E.• Stavroulakis for their guidance.

The research of Samir Adly and Michel Thera has been supported by the "Fan dation EADS" and the ANR project "Guidage".


Qinggang Wu
email: qgwu@nuaa.edu.cn

Minghua Hu
email: minghuahu@nuaa.edu.cn

Xiaozhen Ma
email: xzma@nuaa.edu.cn

Yanjun Wang
email: ywang@nuaa.edu.cn

Wei Cong
email: congwei@variflight.com

Daniel Delahaye
email: daniel@recherche.enac.fr

Modeling flight delay propagation in airport and airspace network

Keywords: airport and airspace network, flight delays, delay propagation

An Airport-Sector Network Delays model is developed in this paper for flight delay estimation within air transport network. This model takes both airports and airspace capacities into account by iterating among its three main components: a queuing engine, which treats each airport in the network as a queuing system and is used to compute delays at individual airport, a Link Transmission Model, which computes delays at individual sector and transmits all air delays into ground delays, and a delay propagation algorithm that updates flight itineraries and demand rates at each airport on the basis of the local delays computed by the queuing engine and flow control delays computed by the Link Transmission Model. The model has been implemented to a network consisting of the 21 busiest airports in China and 2962 links that represent to 151 enroute control sectors in mainland China, and its performance is evaluated by comparing with the actual delay data and results of Airport Network Delays model. It is found that the proposed model is well-suited for simulating delays in air transport system where either airports or airspace could be the bottleneck of the system.

I. INTRODUCTION

Minimizing flight delay is one of the major challenges in the air transportation industry today. Flight delay is usually defined as the difference between scheduled and actual arrival time. A complete flight consists of several segments: pushing back from origin airport gate, taxiing out to the runway, taking off, climbing, cruising through several sectors, landing, and finally taxiing to the arrival gate. At the gate, the aircraft waits for turn-around, after which it will continue to its next flight leg. Thus, the determination of delay is not easy since it can be caused by many factors during the whole operation such as airport capacity, severe weather conditions, sector capacity, maintenance, etc. This paper presents a new model, Airport-Sector Network Delays model (ASND), which is a combination of Airport Network Delays (AND) model and Link Transmission Model (LTM) [START_REF] Pyrgiotis | Modelling delay propagation within an airport network[END_REF][START_REF] Cao | Link Transmission Model for Air Traffic Flow Management[END_REF]. The model can approximately compute flight delays at each individual airport and sector in a network. More important, the model can help us to understand how flight delays that occur in a network will impact the performance of the whole system. Three biggest contributors to flight delays in China [START_REF] Shao | The Analysis of Flight Delay in Airport Based on Bayesian Networks[END_REF], predecessor flight delays, airport capacity constraints, and sector capacity constraints, are taken into account in this model.

Early research work focuses on flight delays in a single airport [START_REF] Hansen | Micro-level analysis of airport delay externalities using deterministic queuing models: a case study[END_REF], while it may neglect the impact of predecessor flights from upstream airports. To minimize the effect of upstream flight delays, airlines normally insert some "buffer" in their flight schedule. For instance, when an aircraft is scheduled to be on the ground for turn-around process (after it arrives at gate/apron, there will be a series of process including unloading passengers/goods, refueling, cleaning, boarding etc.) at the airport for 45 minutes, it may actually require only 35 minutes to complete all the turnaround process, thus providing 10 minutes buffer. However, these buffers are unable to absorb longer delays that typically occur on a daily basis, thus leading to the propagation of delay in the network. In order to explore how flight delays propagate within a large network system, Nikolas Pyrgiotis et al. developed an analytical queuing-based network decomposition model [START_REF] Pyrgiotis | Modelling delay propagation within an airport network[END_REF]. The AND model operates by iterating between its two main components: a queuing engine (QE), which treats each airport in the network as a queuing system and is used to compute delays at individual airport, and a delay propagation algorithm (DPA) that updates flight itineraries and demand rates at each airport on the basis of the local delays computed by the QE. AND is well-suited for simulating delays in the United States as it has been validated by operating within a large network consisting of the 34 busiest commercial airports in the continental United States and comparing the system-wide results from AND with real observations of delays [START_REF] Pyrgiotis | Modelling delay propagation within an airport network[END_REF].

However, AND model does not consider sector capacity or compute flight delays caused by airspace constraints. In the United States it is a justified assumption: it has been reported that nearly 84% of air traffic delays are generated by airports in 2001 and 95% in 2005 [START_REF] Chatterji | National Airspace System Delay Estimation Using Weather Weighted Traffic Counts[END_REF][START_REF]Congress Joint Economic Committee. Your Flight Has Been Delayed Again: Flight Delays Cost Passengers, Airlines, and the U[END_REF]. Instead, in China, airspace capacity is a main factor that causes flight delays [START_REF] Pan | Characterizing air traffic networks via large-scale aircraft tracking data: A comparison between China and the US networks[END_REF]. In general, limited available airspace capacity is responsible for up to 20% of the total delays in 2005. By 2015, this proportion has increased to nearly 30% [8]. The main reason could be that air traffic volumes continuously grow at 10% per year, while available airspace for civil aviation is limited. The imbalance between airspace capacity and traffic becomes more and more serious. Therefore, airspace capacity must be taken into account when modeling flight delays propagation in China. P.K.Menon's Eulerian model seems to be an effective tool to model realistic airspace involving multiple traffic stream [START_REF] Menon | Computer-Aided Eulerian Air Traffic Flow Modeling and Predictive Control[END_REF]. Unfortunately, although it can simulate and adjust traffic in airspace, it can hardly perform the adjustment strategies on individual flight. Taking this into account, the combination of Menon's model and AND is extremely difficult. We then turn to the Eulerian-Lagrangian large-capacity cell transmission model(CTM(L)) for air traffic flow developed by Sun et al. [START_REF] Sun | Multicommodity Eulerian-Lagrangian Large-Capacity Cell Transmission Model for En Route Traffic[END_REF], which is control volume based(Eulerian) and takes the origin-destination information of the aircraft(Lagrangian) into account. CTM(L) has three progressive levels. The primary level is link level, which can be viewed as the connection between the entry point and exit point in a specific sector. The number of cells in one link is scaled by a unit time. The secondary level is sector level model, and since there is no interconnection between different links in one sector, the link level model can be spontaneously extended to obtain the sector level model by gathering all links in the same sector. After merge/diverge nodes are added into the network, the tertiary level, multicommodity network level model is created by putting all sector level models together. Then high-altitude routes can be combined by several links from different sectors. The overall number of cells accounting for the length of a route is obtained by averaging flight time of historical trajectories. Since CTM(L) is based on the origindestination information of flights, every individual flight can be distinguished when operating in the airspace, which prompts us with the possibility of a combination of CTM(L) with AND. Since almost all links have length that are longer than 1-minute interval in CTM(L) model and at the next time instant, only a portion of the aircraft in an upstream link will fly to the downstream link, Cao et al. then use links instead of cells as the most basic element to develop a new model named Link Transmission Model(LTM) [START_REF] Cao | Link Transmission Model for Air Traffic Flow Management[END_REF]. By reducing the number of state variables to one tenth of CTM(L), LTM is about six times faster than CTM(L) while computing [START_REF] Cao | Comparison of the Performances of Two Aggregate Air Traffic Models AIAA Guidance, Navigation, and Control Conference[END_REF]. In addition, compared with CTM(L), another advantage of the LTM is that LTM is able to simulate and predict ground delays, by using a special link, defined as the first link of each path to represent the airport. But ground delays simulated by LTM are all caused by sector capacity. Although these delays occur on the ground, they are all transmitted from the air delays. LTM does not model the airport capacity and the influence among airports.

As we can see, AND and LTM seem to be complementary. AND can calculate flight delays due to airport capacity constraints and flight schedule then provide transmission coefficients as parameters into LTM, while LTM can compute delays due to sector capacity constraints between each pair of airports in AND for individual flights. Thus, it is advantageous to combine AND and LTM together. In this paper, a new model is developed taking both airports and airspace capacity into account, the Airport-Sector Network Delays model(ASND). ASND operates by iterating among its three main components: a QE, which is the same as that in AND, a LTM, which computes delays at individual sector and transmits all air delays into ground delays following the control strategy of China, and a delay propagation algorithm (DPA) that updates flight itineraries and demand rates at each airport on the basis of the local delays computed by the QE and flow control delays computed by the LTM.

The main contributions of the present paper are several. First, a unique model is constructed by combining AND and LTM. In the model, both airport capacity and sector capacity are taken into account. Second, A network consisting of the 21 busiest airports in China and 2962 links that refer to 151 continental en route control sectors is constructed. ASND is successfully validated in that network by comparing with the results of AND and the empirical flight delays.

The rest of paper is organized as follows. Section Ⅱ introduces ASND model and main algorithms when running the model. Section Ⅲ describes the simulation and validation of ASND for 21 busiest airports and related sectors. The comparisons of ASND and AND for computing flight delays in two typical days are carried out. Section Ⅳ provides a summary of conclusions and future work with ASND model.

II. MODEL DESCRIPTION

ASND is generally used to simulate one-day operation of flights, beginning at a time when there are little air traffic activities within the whole network (for example, 6 a.m. in China). An entire day is subdivided into 𝑚 sub-periods, 𝑇 1 , 𝑇 2 , ⋯ , 𝑇 𝑚 of equal length. Here, we use 𝑚 = 96, with each sub-period is 15 min. The five main steps of ASND are described as follows.

Step 1. At first, QE is run for all airports respectively in all sub-periods Each airport is modeled as a queuing system and its runway system is modeled as a single server that serves both arrivals and departures. At each airport, demand rates can be obtained from flights itineraries and service rates can be obtained from the capacity of the runway system. Let 𝑝 𝑎,𝑗 (𝑡) denote the state probability that there are 𝑗 aircraft in airport 𝑎 at time 𝑡 .The expected waiting time in queue, 𝑊 𝑎,𝑞 (𝑡), is determined as follows:

𝑊 𝑎,𝑞 (𝑡) ≈ 𝐿 𝑎,𝑞 (𝑡) 𝜇 𝑎 (𝑡) = ∑ (𝑗-1)𝑝 𝑎,𝑗 (𝑡) 𝑁 𝑗=1 𝜇 𝑎 (𝑡)



where 𝐿 𝑎,𝑞 (𝑡) is the expected number of aircraft in airport 𝑎 at time 𝑡, 𝑁 is a sufficiently integer that denotes the capacity of the queue, and 𝜇 𝑎 (𝑡) is the service rate in airport 𝑎 at time 𝑡 . The detail information for computing 𝑝 𝑎,𝑗 (𝑡) can be found in [START_REF] Pyrgiotis | A stochastic and dynamic model of delay propagation within an airport network for policy analysis[END_REF].

Step 2. For flight 𝑓, we define the following variables. Step 3. Update demand rates of each airport by using the adjust departure and arrival times that are calculated above. Then, QE is performed to recalculate the expected waiting time after the time 𝑡 𝑑 , which is the earliest time that demand rate changes. The adjust departure and arrival time after 𝑡 𝑑 should also be update by using (3) and (4) since the expected waiting time may have changed. This step will repeat until no demand rates change during the whole day. So far we have run Pyrgiotis' AND completely. The flight delays obtained can be regarded as under un-capacitated sector capacity.

Step 4. Let 𝑘 be the index of path, and 𝑠 be the index of sector. We define 𝑙 to be the index of link, which can be understood as the connection between the entry point and exit point inside a sector [START_REF] Sun | Multicommodity Eulerian-Lagrangian Large-Capacity Cell Transmission Model for En Route Traffic[END_REF]. The length of a link is scaled by expected travel time of a flight through it. Then we define 𝑇 𝑙 𝑘 to denote the length of link 𝑙 on path 𝑘. In particular, 𝑇 0 𝑘 denotes the length between the flight's takeoff and the entrance time of first en route control sector on path 𝑘. Since most of the links are longer than 1-minute interval and at the next minute, only a proportion of the aircraft in an upstream link will fly to the downstream link, ℎ 𝑙 𝑘 (𝑡) is used to indicate the number of aircraft that transit from link 𝑙 to link 𝑙 + 1 on path 𝑘 at time 𝑡, which can be determined by adjusted flights itineraries that are calculated above [START_REF] Cao | Link Transmission Model for Air Traffic Flow Management[END_REF]. The number of aircraft 𝑥 𝑙 𝑘 (𝑡) in link 𝑙 on path 𝑘 at time 𝑡 can be determined as follows (no flight in the links at 𝑡 = 0):

𝑥 𝑙 𝑘 (𝑡 + 1) = { 𝑥 𝑙 𝑘 (𝑡) -ℎ 𝑙 𝑘 (𝑡 -∑ 𝑇 𝑖 𝑘 𝑙 𝑖=0 ) + ℎ 𝑖 𝑘 (𝑡 -∑ 𝑇 𝑖 𝑘 𝑙-1 𝑖=0 ) , 𝑡 ≥ ∑ 𝑇 𝑖 𝑘 𝑙 𝑖=0 𝑥 𝑙 𝑘 (𝑡) + ℎ 𝑖 𝑘 (𝑡 -∑ 𝑇 𝑖 𝑘 𝑙-1 𝑖=0 ), ∑ 𝑇 𝑖 𝑘 𝑙-1 𝑖=0 ≤ 𝑡 < ∑ 𝑇 𝑖 𝑘 𝑙 𝑖=0 0, 0 ≤ 𝑡 < ∑ 𝑇 𝑖 𝑘 𝑙-1 𝑖=0   𝑥 0 𝑘 (𝑡 + 1) = { 𝑥 0 𝑘 (𝑡) + ℎ 0 𝑘 (𝑡) -ℎ 0 𝑘 (𝑡 -𝑇 0 𝑘 ), 𝑡 ≥ 𝑇 0 𝑘 𝑥 0 𝑘 (𝑡) + ℎ 0 𝑘 (𝑡), 0 ≤ 𝑡 < 𝑇 0 𝑘 
Then, the number of aircraft in sector 𝑠 can be counted as:

𝐸 𝑠 (𝑡) = ∑ 𝑥 𝑙 𝑘 (𝑡) (𝑙,𝑘)∈𝑠



Let 𝐶 𝑠 (𝑡) denote the capacity of sector 𝑠 at time 𝑡. The number of aircraft in the sector should not exceed the sector capacity, this is:

𝐸 𝑠 (𝑡) ≤ 𝐶 𝑠 (𝑡) (7)
 By using LTM, the earliest time 𝑡 𝑟 can be determined when sector overload occurs, then all overflow aircraft should be adjusted to their previous sectors, as well as the previous links. In order to reduce the influence on the demand rates, the flight with later adjusted departure time will be pushed back preferentially. When there is no overload sector, we then calculate air delays of flights which have been pushed back as:

𝑊 𝑓,𝑠 = 𝑡 𝑟 -(𝐴𝐷(𝑓) + 𝑊 𝑜(𝑓),𝑞 (𝐴𝐷(𝑓)) -∑ 𝑇 𝑖 𝑘 𝑙 𝑖=0
)  where 𝑙 is the link that 𝑓 has been adjusted into. Then all air delays should be transmitted into ground delays by performing further adjustment on the adjusted departure time of those flights by: 𝐴𝐷 ′ (𝑓) = 𝐴𝐷(𝑓) + 𝑊 𝑓,𝑠 

Step 5. Update demand rates of each airport by using the adjust departure and arrival time that are changed due to sector capacity constraints. Then from the period 𝑇 𝑟 , run AND again. 𝑇 𝑟 is the period that contains the earliest adjust departure time of flight being pushed back, for it is the earliest period that some other flights' adjusted departure or arrival times in it may also change. After that, we can run LTM again to check the requirement of sector capacity. This loop "LTM-AND-LTM-AND-LTM" continues until no overload sector can be found at any time.

Input aircraft itineraries for a single day demand, capacity profile at each airport and sector, and length of each link.

Run QE for every airport: Calculates the expected delay on landing and takeoff for every T.

Start at T0 (first period of day). To satisfy all airports' and sectors' constraints, the calculations of "LTM-AND-LTM-AND-LTM" may run many times. Some tricks can be used to improve its computing efficiency. For example, from the second time we run LTM, its start time can be set on the earliest adjusted time of flights which have been readjusted, for flights before that time are not changed and the airspace condition will stay the same. When we use DPA to adjust flights after running LTM, we can only run DPA on flights within two hours from the start time, since it is unlikely that there will be no overload sectors during the next two hours. These tricks will greatly reduce running time of ASND.

Run

III. SIMULATION AND RESULTS

Detailed, auto recorded flights data were obtained from Operations Management Center of Civil Aviation Administration of China (CAAC). The database details the flights departure and landing in every Chinese mainland airport, providing a comprehensive picture of air transport in China. Each flight record reports the flight number, execute date, schedule/actual departure (arrival) time, and the unique aircraft register number (tail number). Such data allows us to easily reconstruct the path of each aircraft flying in the network. Additional information about airport runway information, airspace structure, and aircraft performance are also provided, to determine airport capacity, flying time, and turnaround time.

The ASND model then is implemented with Python 3.5. The investigated air transport network consists of the 21 busiest airports in China and 2962 links that refer to 151 mainland en route control sectors (89.35% of total). In addition, a 22nd "virtual airport" acts as an un-capacitated source to generate and "absorb" flights. All flights between the 21 busiest airports and some other airports can be taken into account by having the relevant aircraft fly to/from the virtual airport, which include most of daily flights.

The following five sets of inputs are critical to the ASND model.

 Aircraft itineraries, including aircraft tail number, the aircraft type, the airline, the origin and destination airports and the scheduled departure and scheduled arrival times;

 Expected service rate at each airport, which indicates the expected number of arrivals and departures that can be served per sub-period at each airport;

 Minimum turnaround time, which is the shortest necessary time between a flight and its immediate predecessor flight (unload, clean, refuel, load, etc.);

 Link length, which is the expected travel time between the entry point and exit point inside a sector, can be computed from historical traffic data of all flights in the network (one year in the present work from 1 January to 31 December in 2016);

 Sector capacity, which is the maximum number of aircraft that can exit in the sector at the same time, is determined by experienced controllers, referring to several indicators such as maximum capacity of sector in history, weather condition and so on.

To illustrate and validate the ASND, the results from the model of a day should be compared with the actual operation data. In particular, since sector capacity is the bottleneck of civil aviation in China, we would like to select actual operation data of two contrasting days from 2016 for our demonstration. On one of the two days, the volume of flights is more than average throughout the year, while on the other day, it is less than the average. Ideally, the service rates and sector capacities depend largely on some factors, such as weather, and special use of airspace etc. The two days we choose are both have no severe weather or large-scale military activities, so that they can be considered to approximately have the same condition. Therefore, the input of expected service rate at each airport and sector capacities are the same when using ASND to model for the two days, so that it will make it easier for us to analyze our model. We also compare the results obtained from ASND with AND, to demonstrate the improvement our model when investigating flight delays in China. Finally, as ASND is a macroscopic model of airport delays, rather than to predict individual flight delays, we then validate our model by comparing airport hourly delay obtained from our model and real data. A delayed flight is defined by CAAC as one whose actual arrival time is later than its schedule arrival time by more than 15 minutes [START_REF]Statistics methods of civil aviation flight normality[END_REF].

The flight data tested here include all the flights within China mainland from 0600 CST on 11 July 2016 to 0559 CST on 12 July 2016, containing 10765 flights. Figure 2 shows flight delay rates at the main airports estimated by AND and ASND, and the actual delay rates (the ratio of delay flights to the total arrival flights) at these airports reported by CAAC. We can see that ASND displays a better result than AND, for flight delay rates are closer to the actual delay rates at most airports, especially the hub ones such as Beijing Capital International Airport (ZBAA), Guangzhou Baiyun International Airport (ZGGG) and Shanghai Pudong International Airport (ZSPD). These results may reveal the actual air transportation in China, that nearly 30 percent of delays are caused by sector limitation. Figure 3 shows the hourly delay distributions from 0800 to 2300 on 11 July computed by the two models and the actual hourly delay distribution at ZSPD. Compared with AND, ASND follows a pattern more similar to the observed data.

The results of Wuhan Tianhe International Airport (ZHHH) is another evidence for the validity of ASND. Figure 4 shows the hourly delay distributions from 0800 to 2300 on 11 July computed by the two models and the actual hourly delay distributions at ZHHH. Compared with AND, in some hour periods the delay aircraft estimated by ASND become fewer. Although we just added sector capacity constraints into the network, it can be found that not all airports' delays have become more serious and delays at some airports like ZHHH have even eased and become closer to the real data. It demonstrates that ASND does not simply increase the delays at all airports, but more appropriately simulates the actual operation. The second validation test is based on the flights operated between 0600 CST on 1 January 2016 and 0559 CST on 2 January 2016, containing 9138 flights (less than the daily average). Because of fewer flights on that day, sector capacity constraints should have less influence on flight delays compared with 11 July. Figure 5 shows the flight delay rates at airports estimated by AND and ASND, and the actual delay rates at these airports reported by CAAC. It can be seen that the difference between results of the two models are not obvious and they are both consistent with the actual operation. In fact, the amount of delayed flights obtained by ASND is 2022, 320 more than the result computed by AND, while on 11 July, the increase number is 1572. Figure 6 shows the hourly delay distributions from 0800 to 2300 on 1 January computed by the two models and the actual hourly delay distributions at ZBAA. We can see that expect for the first two hours, ASND and AND obtain similar results, which reveals that sector compacity constraints have a limited impact on delays, even at the busiest airport in China. These results above conform to our expectations. Performance of the proposed model is evaluated in this paper by comparing with the actual delay data and results of AND. It is found that ASND is well-suited for simulating delays in air transport system where either airports or airspace could be the bottleneck of system, such as the system in China. It can be used to explore critical sectors or airports within air transport network and simulate the scenario when some sectors' or airports' capacities decrease suddenly. Since some potential delays can be identified by ASND in air transport network, it is able to support the pre-warning system for serious flight delays and strategies concerning demand management.

Ongoing work involves further validation of ASND. More major airports and sectors can be added into the model to make it more complete. In addition, uncertain capacities can be considered, which means the sector capacity will not be fixed during a piece of time, but randomly fluctuation within a certain range, which will approach to the actual operation better.

  𝑜(𝑓): origin airport of 𝑓 𝑑(𝑓): destination airport of 𝑓 𝑆𝐷(𝑓): scheduled departure time of 𝑓 𝐴𝐷(𝑓): adjusted departure time of 𝑓 𝑆𝐴(𝑓): scheduled arrival time of 𝑓 𝐴𝐴(𝑓): adjusted arrival time of 𝑓 F𝑇 𝑜(𝑓),𝑑(𝑓) : flying time between 𝑜(𝑓) and 𝑑(𝑓) of 𝑓 𝑓′: the predecessor flight to 𝑓, if there is one 𝑡𝑢𝑟𝑛(𝑓 ′ , 𝑓): turnaround time, which is the time between the scheduled arrival time of 𝑓′ and the scheduled departure time of 𝑓 𝑚𝑖𝑛𝑡𝑢𝑟𝑛(𝑓 ′ , 𝑓): a "minimum turnaround time" which is required to complete all the turn-around processes Then the "buffer" time associated with 𝑓 defined as 𝑏𝑢𝑓𝑓𝑒𝑟(𝑓), can be computed as follow: 𝑏𝑢𝑓𝑓𝑒𝑟(𝑓) = 𝑡𝑢𝑟𝑛(𝑓 ′ , 𝑓) -𝑚𝑖𝑛𝑡𝑢𝑟𝑛(𝑓 ′ , 𝑓)  We initially set 𝐴𝐷(𝑓) = 𝑆𝐷(𝑓) and 𝐴𝐴(𝑓) = 𝑆𝐴(𝑓) for all 𝑓. Then if there is a predecessor flight associated with 𝑓, DPA is performed to calculate the adjusted departure time of 𝑓 as follow: 𝐴𝐷(𝑓) = max [𝑆𝐷(𝑓), 𝑆𝐷(𝑓) + 𝐴𝐴(𝑓 ′ ) + 𝑊 𝑑(𝑓 ′ ),𝑞 (𝐴𝐴(𝑓′)  -𝑆𝐴(𝑓 ′ ) -buffer(𝑓)]  And the adjusted arrival time can be calculated for all flights as follows: 𝐴𝐴(𝑓) = max [𝑆𝐴(𝑓), 𝐴𝐷(𝑓) + 𝑊 𝑜(𝑓),𝑞 (𝐴𝐷(𝑓))  +𝐹𝑇 𝑜(𝑓),𝑑(𝑓) ] 

Fig. 1 .

 1 Fig. 1. The schematic of ASND model

Fig. 2 .Fig. 3 .Fig. 4 .

 234 Fig. 2. Delay rates at major airports on 11 July

Fig. 5 .Fig. 6 .

 56 Fig. 5. Delay rates at major airports on 1 January

ACKNOWLEDGEMENT

This research was supported by the National Natural Science Foundation of China (Grant Nos. 61773203, U1833126), and the State Key Laboratory of Air Traffic Management System and Technology (Grant No. SKLATM201707).


Magdalena Oryaëlle Chevrel

Jérémie Labroquère

Andrew Harris

Scott Rowland

PyFLOWGO: an open-source platform for simulation of channelized lava thermo-rheological properties

Keywords: FLOWGO, Python, lava flow, heat budget, rheology

published or not. The documents may come

Introduction

The thermo-rheological properties of lava flowing in a channel depend on the evolution of the moving volume, where viscosity and yield strength are increasing due to cooling and crystallization (e.g. Lipman and Banks 1987, Crisp et al.

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

1994, [START_REF] Cashman | Cooling and crystallization of lava in open channels, and the transition of pahoehoe lava to aa[END_REF]. Harris and Rowland (2001) produced a 1-D model called FLOWGO in which velocity of a lava control volume flowing down a channel is computed via the Jeffreys (1925) equation as modified for a Bingham fluid by Moore (1987). In this approach velocity depends on the lava rheological properties computed according to the cooling and crystallization path of the control volume as estimated via a heat balance box model (Fig. 1). FLOWGO is thus a framework within which thermo-rheological models can be integrated to test fits between output parameters and natural data. By selecting appropriate models to place within this framework Harris and Rowland (2001) succeeded in simulating the down flow heat budget, cooling, crystallinity, viscosity, yield strength, velocity, channel width and maximum length of several lava flows including those of Mauna Loa 1984, Pu'u ' Ō 'ō 1997and Etna 1998. Rowland et al. (2004) later adapted the models contained within FLOWGO to run in a Martian environment allowing cooling-limited, channelized lava flows on Mars to be simulated and their emplacement properties to be inferred [START_REF] Ramsey | What can thermal infrared remote sensing of terrestrial volcanoes tell us about processes past and present on Mars?[END_REF].

During the 14 years since the inception of this approach, the basic physical principles on which FLOWGO is based have not changed. But recently, Harris and Rowland (2015) and Harris et al. (2015) incorporated an alternative model to compute the melt phase viscosity that is based on lava composition, rather than on a given assumed viscosity as originally proposed. They also introduced a three phase rheological model to estimate the effect of crystals and bubbles on viscosity.

To correctly simulate the evolution of thermo-rheological parameters down flow using FLOWGO the user thus is allowed a degree of flexibility so as to best-fit the natural cases, while changing thermo-rheological models and variables within plausible limits (e.g. Harris et al. 2007). Originally, Harris and Rowland (2001)

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

wrote FLOWGO in the programming language IDL (Interactive Data Language) but due to license price and other computing issues this code was set aside and an Excel version was written (officially published in Harris et al. 2015). This was freely shared when needed by other scientists. Although Excel is a convenient tool and is easily and widely used by geologists, it has limited applications, a poor flexibility for model evolution, and when many equations and input parameters are stacked in sequence, it becomes too easy to key in a hidden (or very-hard to find) error. Besides, it cannot be easily incorporated into other software. Lava modeling capabilities and computer processing power has improved over the past decade, and FLOWGO remains often cited, being recognized as the only thermorheological-based model. Some authors have therefore used it to reproduce natural flow evolution of past (e.g. [START_REF] Riker | The length of channelised lava flows: insight from the 1859 eruption of Mauna Loa Volcano, Hawaii[END_REF]Wantim et al. 2013) or ongoing eruptions (Harris et al. 2011;Wright et al. 2008), as well as producing hazards maps (Rowland et al. 2005) and applying FLOWGO as a reference to compare results of other models against (e.g. [START_REF] Cordonnier | Benchmarking lava-flow models[END_REF] or as input to develop new probabilistic models (e.g. QLAVAH, Mossoux et al. 2016).

The focus of the present work is thus to provide FLOWGO in a modern and flexible language. We chose Python because it provides useful libraries, is opensource, and its object-oriented approach allows for great flexibility. Python also has been widely adopted in scientific computing during the recent years and has been described as "the next wave in Earth Sciences Computing because it simply enables users to do more and better science" [START_REF] Lin | Why python is the next wave in earth sciences computing[END_REF]. Furthermore, Python can be run on any operating system which guarantees portability. Here, we describe

Model architecture in Python

PyFLOWGO is developed in Python v3 which is an object-oriented programming language. The code has been designed to allow the user to switch between any existing models and add new models as they become available, without modifying the architecture of the code. The software acts as a framework that provides interfaces to implement multiple models, and calls them in the correct sequence to build the lava flow differential equations and solve them using a numerical approach. The interfaces basically define the methods necessary for the solver to work and can be implemented with specific models depending on the desired simulation. The top level of the architecture is the integrator which solves the differential equations depending on heat fluxes and on input physical characteristics of the lava (described by the material lava class), terrain conditions and a given crystallization rate model (Fig. 2). The integrator solves the differential equations and updates accordingly the current lava state (temperature, crystallization, position, etc.) which is then used for the next integration step. This process is iterated until termination conditions are reached. The material lava class is composed of multiple models such as the melt viscosity model, the relative viscosity model, the yield

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT
strength model and the vesicle fraction model. Each model is defined by the same interface that governs inputs and outputs delivered to and from the model. As an example, to compute the melt viscosity, all the models available to the user share a common interface called base melt model viscosity (Fig. 3). This interface makes sure that the model receives the state of the lava in order to deliver the viscosity value in the expected unit, that is Pa s. In the same way, all heat fluxes that compose the differential equation share the same interface called base flux (Fig. 4). In this case, the interface provides a unique method to compute and return the flux in W/m based on the state and channel dimensions as input parameters. With this architecture, new physical models or fluxes can then easily be added by implementing the given interface it depends from. Communication is carried out only between the interfaces, and models can be switched from one to another with no modification of the code structure, thus avoiding implementation errors and allowing a great flexibility.

Modeling

Differential equation for heat budget and crystallization down flow

PyFLOWGO is built around the main differential equation established in FLOWGO which is based on the heat budget for a control volume of lava within a channel (e.g. [START_REF] Danes | Dynamics of lava flows[END_REF][START_REF] A N U S C R I P T A C C E P T E D Accepted Manuscript | Dynamics of lava flow: Thickness growth characteristics of steady two-dimensional flow[END_REF][START_REF] Crisp | Influence of crystallization and entrainment of cooler material on the emplacement of basaltic 'a'a lava flows[END_REF]Keszthelyi 1995b; [START_REF] Keszthelyi | Some physical requirements for the emplacement of long basaltic lava flows[END_REF]

∆H = Q rad + Q conv + Q rain + Q cond -Q cryst -Q visc (1)
where Q rad , Q conv and Q rain represent the heat loss from the surface due to radiation, forced atmospheric convection and rain vaporisation; Q cond is the heat loss by conduction through the flow base and levées; and Q cryst and Q visc are the heat gains due to crystallization and viscous dissipation. Note that heat loss due to entrainment of cold material from the crust into the hotter flow interior could also be involved [START_REF] Crisp | Influence of crystallization and entrainment of cooler material on the emplacement of basaltic 'a'a lava flows[END_REF]. This case is not treated here but all details may be found in Harris and Rowland (2001). Following Keszthelyi (1995b), heat gain from crystallization, Q cryst can be written as:

Q cryst = ∂ T ∂ x ρ bulk L cryst E r ∂ φ ∂ T cool (2)
where ∂ T /∂ x is the cooling per unit length (in K/m), ρ bulk (kg/m 3 ) is lava bulk density, L cryst (J/kg) is latent heat of crystallization, E r (m 3 /s) is effusion rate, and

∂ φ /∂ T cool (K -1 )
∂ T ∂ x = -Q rad -Q conv -Q cond -Q rain + Q visc E r ρ bulk L cryst ∂ φ /∂ T cool (3) M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT
where the fluxes are expressed in W/m and detailed in the supplementary material (Appendix A) together with the thermal conditions of the lava (crust temperature and coverage). Each heat flux (Q i ) can be written as a function of the state (Y ) and the position (x), such that: Q i (Y, x). This allows rewriting Eq.3 in the following form:

∂ T ∂ x = F Y, x, ∂ φ ∂ T cool with: F Y, x, ∂ φ ∂ T cool = 1 E r ρ bulk L cryst 1 ∂ φ /∂ T cool ∑ f luxes Q i (Y, x) (4)
The increase of down flow crystal fraction (∂ φ /∂ x) is then computed via:

∂ φ ∂ x = ∂ φ ∂ T cool ∂ T ∂ x (5)
where the crystallization rate per degree of cooling, ∂ φ /∂ T cool , can be calculated from one of the models presented in the supplementary material (Appendix A).

Conservation of volume

Based on mass conservation, if the effusion rate and channel depth are kept fixed down flow, then width can be computed at each down flow step from:

w = E r d V mean (6)
with w and d being the channel width and depth (in m), and V mean being the flow velocity (in m/s) which in turn is controlled by the underlying slope and the rheological properties of the lava. The effusion rate is computed at the vent from initial

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT
flow geometry (width x depth), rheology and velocity.

Flow velocity

The velocity model originally used by FLOWGO is derived from the Jeffreys (1925) equation that was later adapted for a Bingham rheology by Moore (1987).

Although this equation is based on the mean of the velocity gradient inside the channel, here it is used as to estimate a single value to best characterize the velocity of the control volume at each step. Note that this is not the rate of advance of the flow front and is expressed as:

V mean = ρ bulk g d 2 sin θ n η bulk 1 - 3 2 τ 0 τ b + 1 2 τ 0 τ b 3 (7)
where n is the channel shape factor, θ is the underlying slope in radians, g (m/s 2 ) is acceleration due to gravity, ρ bulk (kg/m 3 ) is the lava bulk density, η bulk (Pa•s) is the bulk viscosity of the lava mixture (considering the melt phase and the effect of fraction equations, i.e., Eq. 4 and Eq. 5, respectively. Theses two equations are discretized in the space domain using a linearization at the position x i . By denoting

δ 2 1 Λ = Λ 2 -Λ 1
as the variation in variable Λ from state 1 to state 2, and

Λ i = Λ(x i )
as the variable value at position x i , we obtain the following first order integration (neglecting higher order components):

δ i+1 i T δ i+1 i x ≈ F Y i , x i , ∂ φ ∂ T cool i (8)
which leads to the common Euler scheme for temperature:

T i+1 = T i + (δ i+1 i x)F Y i , x i , ∂ φ ∂ T cool i (9) with δ i+1 i
x being the step distance between i and i + 1. Then the second equation, for the crystal fraction, is solved using the same approach:

φ i+1 = φ i + (δ i+1 i x) ∂ φ ∂ T cool i F Y i , x i , ∂ φ ∂ T cool i (10)
The values of the state (Y i+1 ) now allow computation of the rheology, and hence the velocity of the control volume, at position x i+1 , which is then used to estimate channel width by considering a constant effusion rate E r and a constant depth in the conservation of volume equation (Eq. 6).

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

The thermal conditions and the crystal fraction are thus integrated by propagating the initial state conditions at x 0 down-channel, and by updating velocity and control volume thermo-rheological properties at each step. In practice the step size is fixed such that ∀i, δ i+1 i x = ∆x. Determination of this step is model and case dependent, and must be chosen to be small enough to remove any numerical error and provide enough accuracy in the variables of interest. A convergence study must therefore be performed for every new studied case.

Data interpolation

At multiple places in PyFLOWGO, interpolation is needed when data are provided as discrete values. In particular, the line-of-steepest-descent down which the control volume is moved is extracted from a Digital Elevation Model (DEM) where data are usually provided for every 1, 10, 50 or even every 100 meters (depending on the spatial resolution of the DEM). The same holds for the MELTS data used to estimate crystallization as a function of temperature (see crystallization rate model in supplementary material). In this case, data are usually given with a temperature step of 0.1, 1 or 10 °C. To be independent of data discretization, and thus able to read any type of data, the data has to be interpolated to be set at a common step value. In PyFLOWGO, we use a linear interpolation to reconstruct any missing data during the integration of the differential equations.

Description of the software package

PyFLOWGO is packaged as a module (or library) to be imported and reused in any Python program. The models and simulation configurations are chosen using a single configuration file (json format) that contains all the necessary numerical and As described above (section 2), any new model can be added by implementing the interface it depends from.

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT model parameters (
The line-of-steepest-descent down which the control volume is moved has to be previously extracted from a DEM and input as a x, y text file where x is distance down flow and y is the slope (in degrees) at that point. All computed variables (as a function of distance) can be stored for every step in an output file (CSV format).

Some tools are provided to plot graphs such as crystallization rate, velocity, evolution of the bulk viscosity (interstitial melt + effect of particles), crystal fraction, yield strength, channel width as a function of the distance. Ground-truth data can also be added and plotted with the model results. 2015), respectively for the three test cases (see table A.3 in appendix).

Model Verification

Convergence analysis

Numerical convergence is a necessary step to verify of the numerical solution.

Indeed, the integration scheme propagates a numerical error which can then exceed the model error. The numerical error of an Euler integration scheme is controlled by the integration step size. Convergence analysis consists of reducing the step size and then tracking the value returned for the variable of interest. When the variable seems to stabilize within an acceptable range, then the step size is sufficiently small.

A convergence analysis is performed here for each case by reducing the iteration step size from 100 m down to 1 m. Figure 5 shows the maximum distance attained by the control volume (i.e. the point where mean velocity equals 0 m/s) and the lava core temperature near the end of the flow with respect to the step size. For the three cases the convergence for distance is reached at step size less than 10 m with an error of less than 50 m. For ML84 and MU74, the temperature convergence is reached at step sizes of less than 10 m, with an error within ∼ 2-3

°C, but for PdF2010, temperature convergence seems to be reached only at smaller step sizes (Fig. 5). For these three test cases, we conclude that a 10 m step is a maximum step size that must be used to run PyFLOWGO in order to guarantee small numerical errors on simulated lava flow properties and dimensions. Initial the previous section a convergence study must be performed for every new studied case. Here one can see that using an appropriate step size is necessary to avoid large errors (in this case the distance reached differs by 1 km on a total distance of 6.5 km).

Results and validation against FLOWGO

Conclusion

This paper describes PyFLOWGO, a software written in Python to run FLOWGO, a thermo-rheological framework for lava flowing in a channel as originally presented by Harris and Rowland (2001). PyFLOWGO is constructed in a similar 2001). The lava viscosity and yield strength are estimated within the control volume according to the lava state within the box (including thermal state: heat budget, temperature of core, base, surface, crust; and physical state: crustal coverage, crystallinity, vesicularity) in order to compute the velocity and corresponding channel width (for a fixed effusion rate) used for the next step. This model assumes a "cooling limited" lava flow behavior: the lava stops flowing because it has cooled to such an extent that its rheological behavior impede motion. Q rad , Q conv and Q rain are heat losses into the atmosphere due to radiation, forced convection due to heating of the air above the lava surface and effect of rain, respectively. Q cond is the heat loss by conduction into the cooler base and levees. Q cryst and Q visc are the heat gain due to crystallization and viscous dissipation, respectively. Dimensions d and w are the channel depth and width, and θ is the underlying ground slope.

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

Einstein-Roscoe model

This model calculates the effect of crystals on viscosity according to the Einstein-Roscoe relationship, as first introduced by Shaw (1965) and as used in the original FLOWGO version: .19) where here R = 1.51 for spherical solid particles, as suggested by [START_REF] Pinkerton | Methods of determining the rheological properties of magmas at sub-liquidus temperatures[END_REF], this equation is therefore only applicable for spherical particles and for a volume fraction maximum of 0.66 (i.e., 1/R).

η r = (1 -Rφ ) -2.5 (A

Krieger-Dougherty model

This model calculates the effect of crystals on viscosity according to the Krieger-Dougherty relationship [START_REF] Krieger | Rheology of monodispersed latices[END_REF][START_REF] Krieger | A mechanism for non-Newtonian flow in suspensions of rigid spheres[END_REF], Pabst 2004):

η r = (1 -φ /φ m ) -bφ m (A.20)
where b is the Einstein coefficient (also termed intrinsic viscosity) and φ m is the crystal maximum packing, both being fitting parameters that depend on particle shape. In theory, for spherical particles (aspect ratio of 1) this relationship reduces to Eq. A.19. For elongated particles of aspect ratio of approx. 9, Mueller et al. to Maron and Pierce (1956):

η r = (1 -φ /φ m ) -2 (A.21)
where φ m is a fitting parameter that depends on particle shape (e.g. Mueller et al. for spherical particles and φ m = 0.339 for elongated particles with an aspect ratio of approx. 9.

Costa model

The Costa model allows the effect of crystal fraction in an intermediate range of crystallinity (30 to 80 vol % crystals) to be calculated by taking into account applied deformation (strain rate) following [START_REF] Costa | A model for the rheology of particle-bearing suspensions and partially molten rocks[END_REF]:

η r = 1 + φ φ * δ (1 -F) bφ *
in which:

F = (1 -ξ )erf √ π 2(1 -ξ ) φ φ * 1 + φ φ * γ (A.22)
here, φ * is the critical solid fraction that is present at the onset of the exponential increase in η r with φ ; γ is the slope of the relation between η r and φ as the crystal fraction approaches φ * , and δ is the slope of the relations for values of φ greater than φ * . ξ , γ and δ are all empirical parameters that depend on particle shape and applied strain rate. Two default models are offered in PyFLOWGO: costa1 that

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT
is only applicable for spherical particles (aspect ratio of 1) and costa2 that is for elongated particles (aspect ratio of 9). Both models can be used for strain rate set either at 1s -1 or 10 -4 s -1 and use the values given in [START_REF] Cimarelli | Rheology of magmas with bimodal crystal size and shape distributions: Insights from analog experiments[END_REF]. For example, at 1s -1 , costa1 produces: φ * = 0.67, ξ = 0.01, γ =1.6 and δ =11.4, while costa2 produces: φ * = 0.28, ξ =0.001, γ =8.55 and δ =4.45. More examples can be found in [START_REF] Costa | A model for the rheology of particle-bearing suspensions and partially molten rocks[END_REF][START_REF] Cimarelli | Rheology of magmas with bimodal crystal size and shape distributions: Insights from analog experiments[END_REF] and [START_REF] Chevrel | Lava flow rheology: A comparison of morphological and petrological methods[END_REF].

Phan-Thien and Pham model

PyFLOWGO offers one model that allows the treatment as a three-phase mixture comprising a suspension of rigid spherical particles (φ ) and bubbles (φ b ) following [START_REF] Phan-Thien | Differential multiphase models for polydispersed suspensions and particulate solids[END_REF]. This model is applicable only for φ +φ b < 1.

One of the three following cases can be applied:

Case ptp1, crystals are smaller than bubbles:

η r = 1 - φ 1 -φ b -5/2 (1 -φ b ) -1 (A.23)
Case ptp2, crystals and bubbles are the same size:

η r = (1 -φ -φ b ) 5φ -2φ b 2φ -φ b (A.24)
Case ptp3, crystals are larger than bubbles: strain rate, only for "costa1" and "costa2" relative viscosity models s -1 0.0001 or 1 *used in case the limiting conditions (v mean = 0 or φ = φ max or T core = T solid ) are not reached.

η r = 1 - φ b 1 -φ -1 (1 -φ ) -5/2 (A.
** C H = (U /U) 2 where U is the fraction of wind speed according to [START_REF] Keszthelyi | The initial cooling of pahoehoe flow lobes[END_REF]. Q cond "yes" / "no" convection Q conv "yes" / "no" rain Q rain "yes" / "no" viscous heating Q visc "yes" / "no" models crystallization rate model ∂ φ /∂ T "basic" / "bimodal" / "bimodal f temp" / "melts" melt viscosity model η melt "basic" / "dragoni" / "shaw" / "vft" relative viscosity model η r "er" / "mp" / "kd" / "costa1" / "costa2" / "ptp1" / "ptp2" / "ptp3" yield strength model τ 0 "basic" / "dragoni" crust temperature model T crust "basic" / "hon" / "bimodal" effective cover crust model

  the architecture of our new open-source code, named PyFLOWGO, explaining the various models (heat flux, rheology, crystallization rate, crust temperature, crust M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT cover fraction) that can currently be chosen to set-up a lava flow simulation. So that the model can be trusted as an operational tool with known uncertainty we tested the output against previous iterations of the model. As validation, we followed three cases for which appropriate natural data are available and have been previously tested: Mauna Loa 1984 (Hawaii), Mauna Ulu 1974 (Kilauea, Hawaii) and Piton de la Fournaise 2010 (La Réunion), and results are discussed in terms of model convergence and error.

  . The thermal budget (∆H) represents the balance of the heat fluxes flowing in (gain) and out (loss) of a box model as illustrated in Figure 1. The change in heat content of a unit length per unit time (∆H in J s -1 m -1

  particles, see supplementary for details), and τ 0 (Pa) and τ b (Pa) are respectively the lava yield strength and the basal shear stress. Channel shape factor can be obtained via 3(1 + d/w) 2 (Wilson and Parfitt 1993) and reduces to 3 when the channel is much wider than deeper. Solution of the velocity model now requires definition of ρ bulk , η bulk , τ 0 and τ b where η bulk , τ 0 are functions of ∂ φ /∂ T cool and ∂ T /∂ x and hence dependent on the thermal box model. Details about the models to calculate these variables are given in Appendix A. PyFLOWGO uses the FLOWGO solution (Harris and Rowland 2001; Rowland et al. 2005; Harris and Rowland 2015) to solve the lava flow temperature and crystal

  Three lava channels, Mauna Loa 1984 (ML84), Mauna Ulu 1974 (MU74) and Piton de la Fournaise 2010 (PdF2010) have been chosen as test cases from previously published FLOWGO papers. These flows were chosen because their input parameters are well constrained and because they cover a large range of lava flow characteristics: ML84 being a long channel of more than 25 km, MU74 being a mature channelized flow of about 8 km in length and PdF2010 being slower and cooler than the Hawaiian cases and of only 1-2 km in length. For these cases, using the same input parameters as provided inHarris and Rowland (2015),Robert et al. (2014) and Harris et al. (

  convergence tests by Harris and Rowland (unpublished) and by us here indicate 10 m as being an optimum step size in terms of errors and run time. As mention in

Figure 6

 6 Figure6plots lava core temperature, bulk viscosity and mean velocity computed with PyFLOWGO against the results from FLOWGO Excel spreadsheets prepared for the same test cases byHarris and Rowland (2015),Robert et al. (2014) and Harris et al. (2015), using identical input parameters (see tableA.3 in appendix). For all the cases, PyFLOWGO reproduces the FLOWGO results perfectly. Note that the oscillations in mean velocity for ML84 and Pdf2010 come from the small spatial resolution of the line of steepest descent. To obtain less noisy results, one could filter the DEM data. For Mauna Ulu 1974, we used the same slope path as the one provided byRobert et al. (2014), specifically one value every 200 m. The PyFLOWGO run at this 200 m step size thus reproduces very well the results ofRobert et al. (2014), but a comparison with PyFLOWGO run at a step size of 10 m shows that convergence was not actually reached using the 200 m step size. Here one can see that using an appropriate step size is necessary

  to allow estimation of all parameters involved in the thermorheological evolution of a control lava volume flowing down a channel. We present here the architecture of the code, as well as the discretized formulation of the channelized lava flow problem and the various models that can be selected according to the study case. This new code is written with the object-oriented programming language Python v3, and offers more flexibility while reducing the risk of making error when changing models in comparison to the previous FLOWGO version which was written in Excel. The user can run PyFLOWGO using already implemented models, or extend the code with new models by simply implementing the base classes. Communication through interfaces allow tests of different models on the same study case, without modifying the code architecture. This software, used as a model testing platform, also allows the user to easily and quickly set up new complex cases of lava flow simulation to test.PyFLOWGO has been successfully validated against FLOWGO via three test cases (Mauna Loa 1984, Mauna Ulu 1974 and Piton de la Fournaise 2010). For each test case, a convergence study has been performed, which is an essential practice that must be conducted for each new study. PyFLOWGO has also been unit tested and packaged in a Python library form to allow ease of installation. This software is open-source, thus available at any location and institution. .,Mouginis-Mark, P., 1994. Temperature of an active lava channel from spectral measurments, kilauea volcano, hawaii. Bulletin of Volcanology 56, 297-301. Fulcher, G., 1925. Analysis of recent measurements of the viscosity of glasses. J. Am. Ceram. Soc. 8, 339-355. Ghiorso, M.S., Sack, O., 1995. Chemical mass transfer in magmatic processes IV.A revised and internally consistent thermodynamic model for the interpolation and extrapolation of liquid-solid equilibria in magmatic systems at elevated temperatures and pressures. Contrib. Mineral. Petrol. 119, 197-212. Giordano, D., Russell, J.K., Dingwell, D., 2008. Viscosity of magmatic liquids: A model. Earth and Planetary Science Letters 271, 123-134. Greeley, R., Iversen, J.D., 1987. Measurements of wind friction speeds over lava surfaces and assessment of sediment transport. Geophysical Research Letters 14, 925-928. Harris, A., Dehn, J., Patrick, M., Calvari, S., Ripepe, M., Lodato, L., 2005. Lava effusion rates from hand-held thermal infrared imagery: an example from the june 2003 effusive activity at stromboli. Bulletin of Volcanology 68, 107-117. Harris, A., Favalli, M., Mazzarini, F., Pareschi, M., 2007. Best-fit results from application of a thermo-rheological model for channelized lava flow to high spatial resolution morphological data. Geophysical research letters 34, L01301. Harris, A., Favalli, M., Wright, R., Garbeil, H., 2011. Hazard assessment at Mount hybrid lava flow inundation model and satellite-based land classification. Natural Hazards 58, 1001-1027. Harris, A., Flynn, L., Keszthelyi, L., Mouginis-Mark, P., Rowland, S., Resing, J., 1998. Calculation of lava effusion rates from landsat tm data. Bulletin of Volcanology 60, 52-71. Harris, A.J.L., Rhéty, M., Gurioli, L., Villeneuve, N., Paris, R., 2015. Simulating the thermorheological evolution of channel-contained lava: FLOWGO and its implementation in EXCEL. Detecting, Modelling and Responding to Effusive Eruptions. Eds, Harris AJL, De Groeve T, Garel F and Carn SA, Geological Society, London, Special Publications 426. Harris, A.J.L., Rowland, S.K., 2001. FLOWGO: a kinematic thermo-rheological model for lava flowing in a channel. Bulletin of Volcanology 63, 20-44. Harris, A.J.L., Rowland, S.K., 2015. FLOWGO 2012: An updated framework for thermorheological simulations of channel-contained lava. Hawaiian Volcanoes: From Source to Surface, Geophysical Monograph 208, Eds, Carey R, Cayol V, Poland M, and Weis D, American Geophysical Union .Hon, K., Kauahikaua, J., Denlinger, R., Mackay, K., 1994. Emplacement and inflation of pahoehoe sheet flows: Observations and measurements of active lava flows on Kilauea Volcano, Hawaii. Geol. Soc. Am.Bull. 106, 351-370. Hulme, G., 1974. The interpretation of lava flow morphology. Geophys. J.

Figure 1 :

 1 Figure 1: Schematic view of the thermo-rheological model FLOWGO illustrating the heat box model of the control volume of lava advancing through a channel (modified from Harris and Rowland,2001). The lava viscosity and yield strength are estimated within the control volume according to the lava state within the box (including thermal state: heat budget, temperature of core, base, surface, crust; and physical state: crustal coverage, crystallinity, vesicularity) in order to compute the velocity and corresponding channel width (for a fixed effusion rate) used for the next step. This model assumes a "cooling limited" lava flow behavior: the lava stops flowing because it has cooled to such an extent that its rheological behavior impede motion. Q rad , Q conv and Q rain are heat losses into the atmosphere due to radiation, forced convection due to heating of the air above the lava surface and effect of rain, respectively. Q cond is the heat loss by conduction into the cooler base and levees. Q cryst and Q visc are the heat gain due to crystallization and viscous dissipation, respectively. Dimensions d and w are the channel depth and width, and θ is the underlying ground slope.

Figure 2 :Figure 3 :

 23 Figure 2: PyFLOWGO UML class diagram -top level. The interfaces (labeled "I" and with the prefix "base") provide parameters to main classes (labeled "C") that enable the Integrator to update the flow State at discrete positions along a slope.

Figure 4 :Figure 6 :

 46 Figure 4: Example of PyFLOWGO UML class diagram describing the interface, base flux, for the various heat fluxes (see Appendix A for details about the fluxes). The user is free to chose which fluxes to consider (see Table A.2 in appendix for the available fluxes at this date) or implement a new flux.

( 2010 )

 2010 give b = 6.07 and φ m = 0.343. See for more examples Mueller et al. 2010,[START_REF] Cimarelli | Rheology of magmas with bimodal crystal size and shape distributions: Insights from analog experiments[END_REF] Mader et al. 2013. This model calculates the effect of crystal cargo on relative viscosity according

  2010 and Mader et al. 2013). For example, Mueller et al. (2010) use φ m = 0.633

  25) Appendix A.5. Vesicle fraction models, φ b Two vesicle state models are available. The first is a simple model whereby the vesicle fraction is held constant down flow and is equal to the initial (at vent) value.

  Table A.1 in appendix). The user can therefore select which heat fluxes to consider and the associated models for lava thermal condition down flow including effective crust cover fraction, crust and uncrusted surface temperature (see TableA.2 in appendix). The package also include models for crystallization per degree of cooling, density, melt viscosity, relative viscosity, vesicle fraction and for yield strength and basal shear stress (see Appendix A for details).

  Table A.3: Input parameters of the json files used for the test cases

	"lava name"	ML84 HR2015 cold MU74 Robertetal2014 PdF2010 srtm
	"slope file"	ML84-slope file.txt	DEM maunaulu74.txt	DEM pdf2010 srtm.txt
	"step size"	10.0	200.0	10.0
	"models"
	"crystallization rate model"	basic	bimodal	basic
	"melt viscosity model"	shaw	vft	vft
	"relative viscosity model"	er	ptp2	er
	"yield strength model"	dragoni	basic	basic
	"crust temperature model"	hon	constant	constant
	"effective cover crust model"	basic	basic	basic
	"vesicle fraction model"	constant	constant	constant
	"heat budget models"
	"radiation"	yes	yes	yes
	"conduction"	yes	yes	yes
	"convection"	yes	yes	yes
	"rain"	no	no	no
	"viscous heating"	no	no	no
	"terrain conditions"
	"width"	21	2.475	4.5
	"depth"	3	2.475	1.4
	"gravity"	9.8	9.8	9.8
	"max channel length"	50000	10000	4000
	"eruption condition"
	"eruption temperature"	1410.15	1438.15	1387.15
	"lava state"
	"critical distance"	0	4190	0
	"position"	0	0	0
	"time"	1	0	1
	"crystal fraction"	0.25	0.1097	0.104
	"density dre"	2724	2900	2970
	"vesicle fraction"	0.140	0.4866	0.64
	"liquidus temperature"	1473.15	1473.15	0
	"radiation parameters"
	"stefan-boltzmann sigma"	5.67e-8	5.67e-8	5.67e-8
	"emissivity epsilon"	0.95	0.98	0.95
	"conduction parameters"
	"basal temperature"	773.15	1273.15	773.15
	"core base distance"	19	19	19
	"convection parameters"
	"wind speed"	5.12259623	5.12259623	5
	"ch air"	0.0036	0.0036	0.0036
	"air temperature"	283.15	293.15	293.15
	"air density"	0.4412	0.4412	0.4412
	"air specific heat capacity"	1099	1099	1099
	"thermal parameters"
	"buffer"	140	0	140
	"crust cover fraction"	0.9023	0.5	1
	"alpha"	-0.1601	0	-0.0076
	"crust temperature"	698.15	1273.15	773.15
	"melt viscosity parameters"
	"shaw slope"	2.36	0	0
	"a vft"	0	-4.7	-4.52
	"b vft"	0	5429.7	5558
	"c vft"	0	595.5	582.9
	"crystals parameters"
	"crystals grown during cooling" 0.45	0.89	0.89
	"solid temperature"	1243.15	1268.15	1237.15
	"latent heat of crystallization"	350000	350000	350000
	"crystallization rate 1"	-	0.003	-
	"crystallization rate 2"	-	0.025	-

Acknowledgements

The authors greatly acknowledge Simone Tarquini and an anonymous reviewer for revisions that improve the clarity of the manuscript. M.O.C acknowledges the Auvergne fellowship program for full support. This is ClerVolc publication number 270.

ACCEPTED MANUSCRIPT

Appendix A. Models

A C C E P T E D ACCEPTED MANUSCRIPT

Following Harris and Rowland (2001), PyFLOWGO also allows a model that considers the eruption temperature (T erupt ) instead of T 0 in Eq.A.27. In the PyFLOWGO framework these are the dragoni and basic model, respectively.

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

Textural and rheological evolution of basalt flowing down a lava channel. Bulletin of Volcanology 76, 824. Rowland, S.K., Garbeil, H., Harris, A., 2005. Lengths and hazards from channelfed lava flows on Mauna Loa, Hawai-i, determined from thermal and downslope modeling with flowgo. Bulletin of Volcanology 67, 634-647. Rowland, S.K., Harris, A., Garbeil, H., 2004. Effects of martian conditions on numerically modeled, cooling-limited, channelized lava flows. Journal of Geophysical Research 109, E100101.

Ryerson, F.J., Weed, H.C., Piwinskii, A.J., 1988. Rheology of subliquidus magmas: I Picritic compositions. J. Geophys. Res. 93, 3421-3436. Shaw, H., 1965. Comments on viscosity, crystal settling, and convection in granitic magmas. Am. J. Sci. 263, 120-152. Shaw, H.R., 1972. Viscosities of magmatic silicate liquids: An empirical method of prediction. Am. J. Sci. 272, 870-893. Tammann, G., Hesse, W., 1926. Die Abhängigkeit der Viskosität von der Temperatur bei unterkühlten Flüssigkeiten. Z. Anorg. Allg. Chem. 156. Vogel, D., 1921. Temperaturabhängigkeitsgesetz der Viskosität von Flüssigkeiten. Phys. Zeit 22, 645-646. Wantim, M., Kervyn, M., Ernst, G., del Marmol, M., Suh, C., Jacobs, P., 2013. Numerical experiments on the dynamics of channelised lava flows at mount The basic model proposed by Harris and Rowland (2001) takes into account the amount of crystallization during flow (φ grown ). This is defined as the amount of crystallization occurring between the eruption temperature (T erupt ) and the temperature at which the lava cannot flow anymore (T solid ):

Bimodal model as function of distance

This bimodal model was proposed by Robert et al. (2014) and Harris and Rowland (2015) and allows the crystallization rate to be changed after a given distance

where x critic , and the constants C 1 and C 2 are of the user's choice. This bimodal model allows the crystallization rate to be changed after a given temperature (T critic ) is reached.

where T critic , and the constants C 1 and C 2 are of the user's choice.

Appendix A.1.4. MELTS model

The MELTS model allows the crystallization rate per degree of cooling to be set from a MELTS-based look-up table as suggested by Harris and Rowland (2001), Harris and Rowland (2015) and [START_REF] Riker | The length of channelised lava flows: insight from the 1859 eruption of Mauna Loa Volcano, Hawaii[END_REF]. The look-up table is a file containing the amount of crystals (fraction) as a function of temperature (in °C) that must be previously built using the MELTS software of Ghiorso and Sack (1995). A linear interpolation of these data is computed by PyFLOWGO and gives a function (φ interp ) that represents the fraction of crystals grown as a function of temperature.

The fraction of crystals grown per degree of cooling is then computed using the finite differences via the interpolated function:

with ∆T being the temperature step that is chosen to be small enough (for example 10 -6 ). Note that this model considers that crystallization happens under equilibrium conditions, which is probably not the case during lava emplacement (e.g., [START_REF] Chevrel | Lava flow rheology: A comparison of morphological and petrological methods[END_REF][START_REF] Kolzenburg | In situ thermal characterization of cooling/crystallizing lavas during rheology measurements and implications for lava flow emplacement[END_REF].

Heat loss due to radiation from the lava surface to the atmosphere is expressed as:

where σ (W/m 2 K 4 ) is the Stefan -Boltzmann constant, ε is emissivity, w is the channel width and T e f f (K) is the effective surface temperature, which is calculated using a two-component model for the lava surface [START_REF] Pieri | Eruption rate area, and length relationships for some hawaiian lava flows[END_REF][START_REF] Crisp | A model for lava flows with two thermal components[END_REF][START_REF] Pieri | Thermal radiance observations of an active lava flow during the june 1984 eruption of mount etna[END_REF]:

where T atmo is the temperature of the surrounding atmosphere, f crust is the fraction of crusted lava, T crust is the cool crust temperature, 1f crust represents the fraction of exposed uncrusted hot lava and T hot is the hot component temperature. The different models used to calculate f crust , T crust and T hot are described in sections Appendix A.7, Appendix A.8 and Appendix A.9, respectively.

Appendix A.2.2. Forced convection heat flux

Heat loss due to forced atmospheric convection from the lava surface is calculated via (e.g. [START_REF] Keszthelyi | Observations of the effect of wind on the cooling of active lava flows[END_REF]:

where h conv is the convective heat transfer (in W/m 2 K) and T conv (K) the characteristic surface temperature. The convective heat transfer depends on atmospheric conditions and can be defined as:

where U is wind speed (m/s), C H the wind friction factor as defined by Greeley and Iversen (1987), ρ atmo (kg/m 3 ) is atmospheric density and Cp atmo the heat capacity of the air (J/kg K) in contact with the lava surface. The characteristic surface temperature is calculated via:

Heat flux due to rain

The heat flux due to vaporization of rainwater falling onto the lava surface is expressed by:

where ∂ R/∂t (m/s) is the rainfall rate and ρ H 2 O (kg/m 3 ) and L H 2 O (J/kg) are, respectively, the density and latent heat of vaporisation of water.

Appendix A.2.4. Conductive heat flux

The heat flux through the base and the levées of the flow occur via conduction and is expressed as (after Keszthelyi 1995a):

where κ lava is the thermal conductivity of the lava (in W/mK), T core the lava core temperature (in K), T base (K) the temperature at the base of the basal layer and h base (m) the thickness of the basal layer that is defined between the underlaying surface and the thermal boundary when T core ) is reached. It is usually calculated via:

where H b is the proportion occupied by the basal layer in respect to the entire flow thickness (d, in m).

Appendix A.2.5. Viscous heating

Viscous heating in the lava channel is expressed here, for a channel that is wider than it is deep (w > d) following [START_REF] Costa | Viscous heating in fluids with temperaturedependent viscosity: implications for magma flows[END_REF]:

where η bulk (Pa•s) is the bulk viscosity of the molten lava as calculated in section Appendix A.4 and V mean is the mean velocity of the lava as calculated in section 3.3.

Appendix A.3. Density model

PyFLOWGO provides one model to calculate the bulk density:

where ρ DRE is the density of the dense rock equivalent and φ b is the volume fraction of bubbles in the lava obtained via the methods described in section Appendix A.5.

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

Appendix A.4. Viscosity models

Magma and lava are complex systems composed of a polydispersed particle mixture of crystals and bubbles of various shapes and sizes in a liquid phase (the silicate melt). The viscosity of this mixture may be defined as:

where the viscosity of the interstitial melt, η melt (Pa s) is Newtonian and depends on temperature and composition, and the relative viscosity, η r (dimensionless) is obtained by the ratio η bulk /η melt , and depends on the volumetric abundance and aspect ratio of the particles (bubbles and crystals) in the mixture as well as on the strain rate of the flow. PyFLOWGO offers the possibility of calculating the bulk viscosity of the lava using one of four melt viscosity models which can be combined with one of five relative viscosity models.

Appendix A.4.1. Melt viscosity models, η melt

Dragoni and basic model

The Dragoni model calculates the viscosity of the melt at the lava temperature (T core ) using the relation proposed by [START_REF] Dragoni | A dynamical model of lava flows cooling by radiation[END_REF]:

where η 0 (Pa s) is the viscosity of the lava at the liquidus temperature T 0 .

The basic model, as proposed in the original FLOWGO version, is adapted from [START_REF] Dragoni | A dynamical model of lava flows cooling by radiation[END_REF] where instead of liquidus viscosity and temperature, it is the eruption viscosity (η erupt ) and temperature (T erupt ) that are used in Eq.A.16.

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

Shaw model

This model calculates the melt viscosity according to the Arrhenian relationship proposed by Shaw (1972) and reformulated here as:

where s is the characteristic slope of the η melt versus T core relationship that needs to be computed from the melt chemical composition using Shaw (1972).

VFT model

This model is based on the Vogel-Fulcher-Tammann equation (Vogel 1921, Fulcher 1925, Tammann and Hesse 1926) and takes into account the non-Arrhenian behavior of the melt viscosity [START_REF] Dingwell | Volcanic dilemma flow or blow[END_REF] and allows η melt to be calculated via: .18) where A (Pa s), B (J/mol) and C (K) are fitting parameters that depend on chemical composition. These fitting parameters need to have been previously determined either from viscosity measurements at high and low temperature or from the melt chemical composition using for example the model proposed by Giordano et al.

(2008).

Appendix A.4.2. Relative viscosity models

The first four relative viscosity models given here take into account the affect of crystals whereas the fifth model considers those of both crystals and bubbles. More complex formulations may take into account bimodal particle size distribution and

The second is a model as proposed by Harris and Rowland (2015), which allows the vesicle fraction to be changed after a given distance, x critic , and is intended to take into effect down flow degassing. In the PyFLOWGO framework these are the constant and bimodal model, respectively. In the bimodal model:

where x critic , φ b1 and φ b2 are the proximal and distal vesicularities, respectively and can be set using down flow assessments of lava density (e.g. Robert et al. 2014).

Appendix A.6. Yield strength and shear stress model

Velocity depends also on the yield strength of lava and on the basal shear stress (Eq. 7). PyFLOWGO provides one basal shear stress (τ b ) model (Hulme 1974):

where τ b is in Pa.

Lava yield strength can instead be calculated as a function of temperature and crystallinity following [START_REF] Dragoni | A dynamical model of lava flows cooling by radiation[END_REF], and [START_REF] Pinkerton | Methods of determining the rheological properties of magmas at sub-liquidus temperatures[END_REF] as proposed in the original version of FLOWGO. PyFLOWGO uses this approach of allowing τ 0 to be calculated as function of lava temperature (T core ) using the liquidus temperature (T 0 ) and the lava crystal content (φ ) as proposed by Ryerson et al. (1988) in:

. Effective crust cover fraction model

The upper surface of the lava is partially covered by a cooler crust. The fraction of this crusted lava is termed as the effective crust cover fraction, f crust , and varies between zero (crust free, no insulation: rare in nature) and one (complete crust coverage, well-insulated; but not equivalent to a lava tube). Effective crust cover fraction directly affects the effective surface temperature (Eq. A.6) and the characteristic surface temperature (Eq. A.9) which, in turn, influences the heat fluxes due to radiation and forced convection. PyFLOWGO offers two models to calculate

The basic model, as proposed in the original version of FLOWGO, allows f crust to vary down flow as function of velocity:

where f init is the initial (at vent) crust fraction and α is a coefficient that varies crust cover as function of V mean : crust cover increases as flow velocity decreases.

Based on examination of aerial photographs of active channels flowing at known velocities, Harris and Rowland (2001) derived f init of 0.9 and α of -0.16, for poorly insulated flow, and f init of 1.0 and α of -0.00756 for more heavily crusted flow.

Alternatively, f crust can be held constant down flow, equals f init at all down flow location when α = 0.

The second model, named bimodal in PyFLOWGO, allows the dependence of

effective crust cover fraction with velocity to be changed after a given distance,

x critic , as proposed by Harris and Rowland (2015):

where α 1 and α 2 are the crust cover growth coefficients for proximal and distal In the PyFLOWGO framework, this is the honbimodal model.

Appendix A.9. Uncrusted surface temperature model

The temperature of the uncrusted lava surface, T hot , will be lower than of the flow core (e.g. [START_REF] Calvari | The 1991-1993 etna eruption: chornology and lava flow-field evolution[END_REF], Flynn and Mouginis-Mark 1994, Harris et al. 1998). PyFLOWGO accounts for this difference via:

T hot = T corebuffer (A.34)

where buffer is the temperature difference between the maximum surface temperature and the core temperature, and is set by the user. Based on field measurements using thermocouples and radiometers, Harris and Rowland (2001) give buffer = 140. The buffer value may also be lower; [START_REF] Bailey | The changing morphology of an open lava channel on mt. etna[END_REF]