Thery Raphaele

Llovell Felix

Meyer Xuan

Joulia Xavier

MODELLING OF A DYNAMIC MULTIPHASE FLASH: THE POSITIVE FLASH. APPLICATION TO THE CALCULATION OF TERNARY DIAGRAMS

Keywords: Process System Engineering Distillation, Dynamic simulation, Phase Equilibrium

A general and polyvalent model for the dynamic simulation of a vapor, liquid, liquid-liquid, vapor-liquid or vapor-liquid-liquid stage is proposed. This model is based on the τ -method introduced as a minimization problem by Han & Rangaiah (1998) for steady-state simulation. They suggested modifying the mole fraction summation such that the same set of governing equations becomes valid for all phase regions. Thanks to judicious additional switch equations, the τ -formulation is extended to dynamic simulation and the minimization problem is transformed into a set of differential algebraic equations (DAE). Validation of the model consists in testing its capacity to overcome phase number changes and to be able to solve several problems with the same set of equations: calculation of heterogeneous residue curves, azeotropic points and distillation boundaries in ternary diagrams.

INTRODUCTION

In recent years, dynamic simulation has become a valuable and efficient tool for plant design, analysis and operation. Besides, many chemical processes like batch distillation columns or evaporators or continuous distillation column involve phase equilibrium that may lead to the formation of a new phase or to the disappearance of an existing one because of composition, pressure or temperature changes.

Such issues are particularly significant for the study of distillation column start up and shut down and of heterogeneous azeotropic distillation column where "… the interface between trays having one and two liquid phase(s) moves in response to disturbances [START_REF] Widagdo | Azeotropic Distillation[END_REF])". To detect, handle and manage phase pattern changes in phase equilibrium related dynamic processes, reliable dynamic models are required. As summarized in a review on azeotropic distillation by [START_REF] Widagdo | Azeotropic Distillation[END_REF], "it is important to accurately compute the phase distributions when performing dynamic simulation for heterogeneous azeotropic distillations.". Today, two kinds of approaches are usually proposed, namely the equation solving approach and the global stability analysis [START_REF] Müller | Dynamic multiple-phase flash simulation: global stability analysis versus quick phase determination[END_REF] which combined search of the number of coexisting phases and solving of the relevant set of equations. But the aim of the paper is to present an alternative approach, furthermore a dynamic one, namely a formulation of a flash able to handle any phase configuration: one (vapor or liquid), two (liquid -liquid or vapor -liquid) or three (vapor -liquid -liquid) with the same single set of equations and without a systematic phase pattern investigation.

First, we review briefly the two existing approaches and also former concepts that led us to the presented formulation. Then, we introduce the "positive flash" modelling which is an adaptation of an equation solving approach to dynamic simulation. Third, we validate the model by evaluating its capacity to handle phase number changes (scenario 1) and by addressing practical problems linked to the assessment of distillation feasibility: the computation of residue curves maps for homogeneous and heterogeneous systems (scenario 2). Further implementation in a heterogeneous dynamic simulator is in progress but out of scope of this contribution.

CALCULATION APPROACHES FOR VAPOR -LIQUID -LIQUID

MIXTURE EQUILIBRIUM CALCULATIONS.

Equation solving approaches

The equation solving approach consists in performing a stability test to know the number of phases and in solving a set of non-linear algebraic equations deduced from mass balances and equilibrium relations. This can be done either simultaneously using two imbricated loops or sequentially. As the number of phases changes, so does the set of equations used and such a switching poses serious problems for the integration methods used because the number of equations and of unknown variables change and discontinuities arise [START_REF] Widagdo | Azeotropic Distillation[END_REF]. More specifically, equation solving approaches uses physical insights of the multiphase equilibrium to speed up their calculation.

They rely on necessary equilibrium conditions for the different phase combinations. These necessary conditions can be variously formulated as fugacity-equality conditions (Withson and [START_REF] Whitson | The negative flash[END_REF], isoactivity conditions (Cairns and Furzer, 1990b) or as equilibrium conditions in term of K-values [START_REF] Eckert | Modelling of dynamics for multiple liquid-vapour equilibrium stage[END_REF]. Equation solving methods usually require reasonable computational demand makes them suitable for implementation in a dynamic simulation model of processes based on phase equilibrium (Muller and Marquardt, 1997). Yet, as the number of phase is not known in advance, it is tricky to select the phase number relevant set of equations that is to be solved and numerical problem difficulties may occur due to non-smooth problems or pseudo solutions (Buzzi Ferraris andMorbidelli, 1981, 1982 ;[START_REF] Bullard | Iterated linear programming strategies for nonsmooth simulation : a penalty based method for vapor-liquid equiibrium application[END_REF]. To overcome this difficulty, equationsolving method usually involve two sequential or imbricated loops : first the determination of the number and type of coexisting phases relying on a stability test like the tangent plane criterion [START_REF] Baker | Gibbs energy analysis of phase equilibrium[END_REF], Michelsen, 1982a[START_REF] Michelsen | The isothermal flash problem. Part II: Phase split calculations[END_REF], Cairn and Fürzer, 1990a, 1990b[START_REF] Gupta | A method for the simultaneous phase equilibria and stability calculations for multiphase reacting and non reacting mixtures[END_REF][START_REF] Widagdo | Dynamic Analysis of Heterogeneous Azeotropic Distillation[END_REF][START_REF] Sun Amy | Homotopy-continuation method for stability analysis in the global minimization of the Gibbs free energy[END_REF][START_REF] Jalali-Faharani | Use of homotopy-continuation method in stability analysis of multiphase, reacting systems[END_REF] and second, the calculation of phases composition and ratio solving the correct set of governing equations. For more than two coexisting phases like liquid-liquid-vapor equilibria, these two loops may be repeated several times until the correct phase configuration is found. [START_REF] Bullard | Iterated linear programming strategies for nonsmooth simulation : a penalty based method for vapor-liquid equiibrium application[END_REF] considered at that time that the back and forth stability calculation for multiphase equilibria was not very suitable within equation-oriented simulation programs.

Global stability analysis

Approaches based upon a global stability analysis determine the equilibrium solution corresponding to the global minimum of the Gibbs free energy: a non linear optimization method is implemented so as to minimize the Gibbs free energy of an unknown number of phases [START_REF] Castillo | Computation of chemical and phase equilibrium[END_REF][START_REF] Soares | Three phase flash calculations using free energy minimisation[END_REF], Mac Donald and Floudas, 1995, 1997[START_REF] Jalali-Faharani | Use of homotopy-continuation method in stability analysis of multiphase, reacting systems[END_REF]. This approach is known to present the most accurate results but requires a large computational demand.

To perform dynamic simulation, both equation solving and global stability analysis approaches require two calculation loops; phase number calculation one and equation solving one; to be associated to a dynamic set of equations adapted to the appropriate number of phases. Our approach does not necessitate these two loops: the positive flash model is formulated below so as to handle any number of phases during dynamic simulations. With the limitation, and its solution pointed below, of the existence of trivial solution far away from the boundary between one or two liquid phases, at no point during the simulation is the number of coexisting phases calculated independently. Appearance or disappearance of a phase occurs as a discrete event during the solving of the set of equations.

Furthermore, as a single set of equations is kept and the same number of unknown variables is used whatever the number of phase (one; Vapor or Liquid, two; VL or LL or three; VLL), the model can be easily integrated as an existing DAE solver is used without any particular modification or tuning. In the literature, papers about models for the dynamic simulation of vapor -liquid -liquid trays have expressed the need to devise an efficient solving strategy and algorithms have been shown [START_REF] Rovaglio | Dynamics of Heterogeneous Azeotropic Distillation Columns[END_REF]: a sparse Gear method, [START_REF] Wong | Simulation of Dynamics and Phase Pattern Changes for an Azeotropic Distillation Column[END_REF]: a Runge Kutta method, Widagdo et al., 1992, a Gear method and Perregaard et al., 1992: several methods including Backward Differentiation Formula method for stiff problems). Finally initialization and reinitialization of the model; another significant problem of dynamic simulation of heteroazeotropic distillation is discussed in detail in the paper.

CONCEPTS LEADING TO THE "POSITIVE FLASH" MODELLING

The "negative flash" An approach called the "negative flash" in which the vapor -liquid phase split ratio is allowed to take a value outside the physical domain range was introduced [START_REF] Neoschil | Convergence pressure concept. A key for high pressure equilibria[END_REF]Withson and Michelsen, 1989). Later, this formulation was extended to liquid -liquid dynamic equilibrium [START_REF] Heinichen | Dynamic negative liquid-liquid flash[END_REF]. The resolution of this model leads to a successful solution whatever the phase region (L or LL) as the set of equation is kept unchanged: If the phase split α takes a value between 0 and 1, two liquid phases exist. Otherwise, if α leaves the physical range, the set of equations remains the same but one of the liquid phases is non existent physically. Although this approach leads to interesting results, problems may occur due to the absence of solutions or the existence of multiple ones [START_REF] Müller | Dynamic multiple-phase flash simulation: global stability analysis versus quick phase determination[END_REF].

The Ratchford-Rice equation for multiphase systems

To obtain a model consistent whatever the phase configuration, [START_REF] Michelsen | Calculation of multiphase equilibrium[END_REF] formulated a general constrained optimization problem. The objective function is deduced from material balances, equilibrium and summation equations and includes variables defining the mole fraction and the molar amount of each phase. The latter are limited to the physical range. Later, Neoschil andLeibovici (1992, 1995) addressed the same concept as an unconstrained problem in which phase amounts may lie outside the physical domain. With this approach, the phase configuration may easily and efficiently be deduced from the solution of a set of non-linear equations solved using a Newton Raphson algorithm. Both approaches are very useful for performing oil or gas reservoir simulations with equation of state thermodynamics where millions of isothermal flash calculations are needed. For such an application, the fugacities and activity coefficients are assumed to be composition independent and the equilibrium constants values are known. Both approaches are non-dynamic formulations that enable the disappearance of any number of phases during the process of finding the solution. So far, these approaches are not extended to dynamic simulation of phase appearance or disappearance and to problems where activity coefficients are used instead of equation of state thermodynamics.

The τ-method

The τ-method proposes an alternative approach. It is a steady state model for the computation of phase equilibrium derived by [START_REF] Han | A method for multiphase equilibrium calculations[END_REF] that exploits the fact that the mole fraction summation of a non-existent phase is less than unity. Consequently, to deal with all possible regions in multiphase systems, a new variable denoted τ is defined for each phase and is used to modify the mole fraction summation. If the phase exists, the pseudo phase variable τ is equal to 0; on the contrary if the phase is non-existent, τ takes a value between 0 and 1. Consider the n th adiabatic three-phase stage with nc components as illustrated on figure 1. The relevant equations for modeling the VLL equilibrium in steady state are [START_REF] Han | A method for multiphase equilibrium calculations[END_REF]:

Mass Balances (i=1,…,nc): 0 L L V L L V F II n I n n I 1 n II 1 n 1 n n = - - - + + + + + - (1) 0 x L x L y V x L x L y V z F II i , n II n i , n I n i , n n I i , 1 n I 1 n i , 1 n II 1 n i , 1 n 1 n i , n n = - - - + + + + + + + - - (2) Energy Balance: 0 h L h L H V h L h L H V h F II n II n I n I n n n II 1 n II 1 n I 1 n I 1 n 1 n 1 n F n n = - - - + + + + + + + - - (3)
Equilibrium relations (i=1,…,nc):

I i , n I i , n i , n x K y = (4) II i , n II i , n i , n x K y =
(5)

Mole fraction summations:

∑ = τ + nc i V n i , n 0 y (6) ∑ = τ + nc i I n I i , n 0 x (7) ∑ = τ + nc i II n II i , n 0 x (8)
Because the resulting model has more unknown variables than equations, it is solved as a minimization problem:

( )

II n I n V n min τ + τ + τ (9)
subject to the former constraint and to the inequalities:

1 , , , y , x , x 0 II n I n V n n , i II n , i I n , i ≤ τ τ τ ≤ (10) 0 T , V , L , L n n II n I n ≥ (11)
The set of equations and the number of unknown variables are kept constant whatever the number of coexisting phases. The efficiency of the τ-method has been illustrated for VL, LL and VLL equilibrium calculations via several examples [START_REF] Han | A method for multiphase equilibrium calculations[END_REF][START_REF] Shyamsundar | A method for simulation and optimization of multiphase distillation[END_REF]. Yet, the optimization formulation restricts the model application to steady-state simulations. Nevertheless, a dynamic variant of the τ-method involving a pseudo vapor component has been developed for vaporliquid dynamic flash calculations [START_REF] Marci | Simulation of a VLLE stage using the positive/negative flash formulation[END_REF]. It contains an equal number of variables and equations that can be readily solved using differential algebraic equations (DAE) solvers. In this paper, we extend this dynamic VLE formulation to VLLE systems.

THE POSITIVE FLASH MODELLING

This contribution extends the τ-method to dynamic simulation by adding equations to obtain an equal number of unknown variables and equations. Furthermore, these new equations constrain the pseudo phase variables τ. The resulting model is presented below. It is made of two parts, one that describes the continuous dynamic equations to follow the dynamic evolution of the continuous variables and one that describes the discrete modeling of the phase number switches to cope with changes in the number of coexisting phases.

As suggested in the literature [START_REF] Rovaglio | Dynamics of Heterogeneous Azeotropic Distillation Columns[END_REF][START_REF] Wong | Simulation of Dynamics and Phase Pattern Changes for an Azeotropic Distillation Column[END_REF][START_REF] Widagdo | Dynamic Analysis of Heterogeneous Azeotropic Distillation[END_REF], the usual assumptions for dynamic three phases equilibrium modeling have been made:

1. negligible vapor holdup in the mass and energy balance, 2. perfect mixing in the liquid phases, 3. perfect split of the liquid phases, 4. uniform temperature on each tray, 5. negligible heat of mixing, 6. tray outlet vapor and liquid streams are in equilibrium.

The Positive Flash continuous modeling

Considering an n th multiphase tray, equations are:

Mass balance equations (i=1, nc):

0 F V L L V L L dt dU dt dU n n II n I n 1 n II 1 n I 1 n II n I n = + - - - + + + + - + + (12) 0 z F y V x L x L y V x L x L dt x dU dt x dU i , n n n n II i , n II n I i , n I n i , 1 n 1 n II 1 n II 1 n I 1 n I 1 n II i , n II n I i , n I n = + - - - + + + + - - + + + + (13) Energy balance equation: 0 ) y , P , T ( H V ) x , P , T ( h L ) x , P , T ( h L Q H F ) y , P , T ( H V ) x , P , T ( h L ) x , P , T ( h L dt ) x , P , T ( h dU dt ) x , P , T ( h dU n n n n II n II n II n I n n I n I n n Fn n 1 n n 1 n 1 n II 1 n n II 1 n II 1 n I 1 n n I 1 n I 1 n II n II n II n I n I n I n = - - - + + + + + + - - - + + + + + + (14)
Equilibrium relations (i=1, nc):

( ) I i , n i , n I i , n n I i , n i , n x y , x , P , T K y ⋅ = (15) ( ) II i , n i , n II i , n n II i , n i , n x y , x , P , T K y ⋅ = (16)
Mole fraction summations:

0 y x V n nc i i , n I n nc i I i , n = τ - - τ + ∑ ∑ (17) 0 y x V n nc i i , n II n nc i II i , n = τ - - τ + ∑ ∑ (18)
Model for the liquid Holdups:

0 ) U mod( U U n II n I n = - + (19)
Relation between the liquid flow rates and the holdups:

0 ) L L ( U ) U U ( L II n I n I n II n I n I n = + ⋅ - + ⋅ (20) Additional equations: Switch equations 0 V V n n = τ ⋅ (21) ( ) ( ) 0 dt U d dt d U I n I n I n I n = ⋅ τ + τ ⋅ (22) ( ) ( ) 0 dt U d dt d U II n II n II n II n = ⋅ τ + τ ⋅ ( 23
)
Mod U n means that a model for the total liquid holdup is used (e.g. for tray hydraulic models: Francis weir formula or Mackowiack model [START_REF] Mackowiack | Pressure drop in irrigated packed columns[END_REF]). The last three switch equations permit to obtain a system with equal numbers of equations and unknown variables. Furthermore, equation ( 21)

guarantee that either the pseudo-variable τ in the case of an existing phase, or the quantity of the considered phase in the case of a non-existent phase, is equal to 0. Equation ( 22) and ( 23) guarantee that either τ n and its derivative (resp. U n and its derivative) are equal to zero when the phase exists (resp. does not exists). Boiling (condensation) phenomena that cause (dis)appearance of the vapor phase is physically well understood: V n goes rapidly to a large positive value (resp. to zero) when this occurs. But the appearance or disappearance of a liquid phase is less easy to describe physically. In particular, the liquid holdup are likely to display less brutal variations of their magnitude than the vapor does. To us this complies with stricter equations ( 22) and ( 23) where the holdup derivative goes smootly from (or to) zero when a liquid phase (dis)appear.

Both liquid phases are treated symmetrically in the equilibrium and switch equations.

) K mod( K I i , n I i , n = - (24) 0 ) K mod( K II i , n II i , n = - (25) 0 ) h mod( h I n I n = - (26) 0 ) h mod( h II n II n = - (27) 0 ) H mod( H n n = - (28)
The first modeling requires implementing the thermodynamic functions into the mathematical model subroutine. The second one is more versatile as it enables us to use a thermodynamic interface to a thermodynamic property server. We have used it with the BibPhy thermodynamic library (Prosim, 2000).

The Positive Flash Discrete modeling

The DAE solver DISCo that is used to solve the differential algebraic set of equation is based on a Gear algorithm and is capable of handling index 1 or 2 and sparse sets of equations. It also handles the automatic detection of time or state events and bears an automatic procedure to calculate initial conditions after discontinuities [START_REF] Sargousse | Un Nouvel Environnement de Simulation Orienté Objet[END_REF]. Indeed, although the model is able to describe all phase configurations, the transition between two configurations leads to discontinuities that constitute one of the key difficulties of dynamic simulation of heterogeneous azeotropic distillation [START_REF] Widagdo | Azeotropic Distillation[END_REF]. Consequently, it is necessary to detect the phase transition events and to reinitialize the discontinuous variables accordingly. Figure 2 illustrates through a Finite State

Machine, the phase configuration transitions. Each phase configuration of the system (one liquid, one vapor, liquid -vapor, liquid -liquid -vapor) is represented by a state (a circle) and the transitions between these states are triggered by an event. Each event is described by an equation that will take the value zero when the event occurs:

( )V vap _ Flag 1 vap _ Flag G V n 1 ⋅ - + τ ⋅ = (29) ( ) I n I n 2 U 1 liq _ Flag 1 1 liq _ Flag G ⋅ - + τ ⋅ = (30) ( ) II n II n 3 U 2 liq _ Flag 1 2 liq _ Flag G ⋅ - + τ ⋅ = (31) II n I n 4 U U G + = ( 32
)
Four kinds of events are related respectively to the four equations above:

-Appearance (Flag_vap=0) or disappearance (Flag_vap=1) of the vapor phase -Appearance (Flag_liq1=0) or disappearance (Flag_liq1=1) of the liquid I phase -Appearance (Flag_liq2=0) or disappearance (Flag_liq2=1) of the liquid II phase -Tray empty or full Let's take for example, the case of the vapor phase mixture. If we are initially below the mixture boiling temperature, no vapor exists (V = 0), Flag_vap = 1 and τ v = 0. When the liquid boils, vapor rises which tricks Flag_vap = 1, τ v = 0 and V ≠ 0.

The reinitialization procedure is adapted to the number of phases existing after the phase transition detected. In particular, derivatives of the molar fractions are reinitialized to non zero values, preferably to unity. Liquid holdup derivatives are set to zero to comply with switch equations ( 22) and ( 23). If the enthalpies are considered as unknowns as well, their derivatives are also reinitialized. The integration step is also reduced at each reinitialization.

ILLUSTRATIVE EXAMPLES

The use of a dynamic multiphase flash formulation is particularly interesting in the study of distillation processes. During conception and design steps, residue curve maps are thoroughly analyzed to assess distillation feasibility and distillation column sequences; at the simulation level, the dynamic behavior of a distillation column offers insights to its practical operation, startup and shutdown and may exhibit changes in the number of phases in equilibrium on the column trays. In this paper, we validate the model by evaluating its capacity to handle phase pattern changes and by computing residue curves to draw maps for homogeneous and heterogeneous ternary systems as well as to predict azeotropes.

Scenario 1. Detection of phase configuration transitions

In this example, the positive flash model is applied to the dynamic simulation of a vapor -liquidliquid equilibrium Rayleigh distillation of the water -ethanol -cyclohexane mixture (figure 3). The NRTL thermodynamic model is used with vapor -liquid -liquid binary parameters taken from the DECHEMA (see Table 1) [START_REF] Gmehling | Vapor -Liquid Equilibrium Data Collection[END_REF]. For this system, the model is similar to the former one without equations ( 19) and ( 20) related to the holdups. A scenario is described on figure 4.

Two key parameters of the process, the heat duty Q and the pure ethanol feed flow rate F ethanol , are sequentially modified to trigger off the apparition or the disappearance of one of the phases:

-First, the mixture is heated at a constant rate from its initial liquid -liquid state until it boils and enters the vapor -liquid -liquid region.

-Further heating removes a vapor rich in ethanol that moves the mixture composition across the vapor -liquid -liquid boiling envelope and settles it in the vapor -liquid region.

-Then, the heat duty is reduced until vapor condenses once the mixture temperature reaches the boiling temperature. Further cooling causes the apparition of a second liquid phase.

-Finally, pure ethanol is added to the mixture to make disappear the aqueous liquid phase.

As illustrated in figure 5, this scenario permits to test the LL/VLL, VLL/VL, VL/L, L/LL and LL/L transitions. At the transitions, the unknown variables which present discontinuities are reinitialized (e.g. the vapor flow rate V for the LL / VLL transition). As can be seen on figure 5, whenever a phase disappears, its corresponding τ value becomes strictly positive: At the beginning, the vapor phase does not exist then its corresponding τ V value is positive. τ V becomes zero when the two liquid phases mixture boils up while V becomes positive. This event is triggered by equation ( 29) meanwhile the switch equation ( 21) is always valid. The vapor flow rate V and its derivatives show discontinuities and are reinitialized accordingly. While heating is continued, the mixture composition moves away and exits the VLL region as the phase II disappear: U II becomes zero and τ II becomes positive. Then the mixture is cooled down until the vapor is totally condensed to leave a single liquid phase I: V becomes zero and τ V becomes positive. Cooling proceeds, forcing the composition to enter again the two liquid phase region: U II becomes positive and τ II equals zero. Finally, ethanol is added to exit the LLE region and end with a single liquid mixture. Liquid phase numbering is not known in advance. It is a result of the calculation.

Two comments are worthy: first, any phase pattern change is readily assessed without any need to implement a phase stability test; second, the number of unknown variables and equations of the models remains the same whatever the phase number configuration.

Scenario 2. Solving the initialization problem.

In scenario 1, the initial mixture composition is in the liquid -liquid equilibrium region. The positive flash model readily finds correctly that two phases coexists and that the pseudo phase variable τ I and τ II equals zero. However, further tests showed that in most of the vapor -liquid region, the most attractive solution is not a solution where one of the liquid phases has disappeared and the corresponding pseudo phase variable τ has some positive value. It is rather the so-called 'trivial solution': both liquid phases have equal compositions and both τ I and τ II equals zero This readily makes the two equilibrium (equations ( 15) and ( 16)) and the two liquid switch (equations ( 22) and ( 23)) equations identical two by two and the set of equations looks like a usual VLE dynamic model.. This is annoying because once it is found by the DAE solver, the trivial solution is always found even if the composition may move into the vapor -liquid -liquid region. In that case, the mixture behaves in the whole composition simplex as a one vapor -one liquid system: residue curves, azeotropic points and distillation boundaries are not computed accurately which may lead to an erroneous assessment of the distillation process conception and design.

Such a problem was already been pointed out by [START_REF] Han | A method for multiphase equilibrium calculations[END_REF]. It has also been known to occur in equilibrium calculation problems where phases are treated symmetrically, for example using identical equation of states for both phases (Ferraris andMorbidelli, 1981, 1982). The general formulation of the positive flash model falls into this category: the model is perfectly symmetrical concerning the two liquid phases and the trivial solution seems to be the most attractive one. Another possible explanation for this problem could also reside in the absence of physical meaning we attribute to the liquid τ values. Indeed, a positive τ V value can always be assimilated to an incondensable (e.g. air) composition that resides above any liquid or liquid -liquid state mixture in consideration. On the other hand, a positive τ I or τ II value could be assimilated to a non volatile component (e.g. coke) unable to be involved in a liquid -liquid equilibrium. Yet, although the presence of a pseudo-vapor phase above a liquid phase is conceivable, the continuous existence of a second pseudo liquid phase remains difficult to imagine if coke is not present initially.

To overcome the problem a method that is able to provide the system with initial values of the liquid τ variables when one comes close to a VLLE region is needed. An adaptation of the τ method is the ideal initialization method candidate as it does consider liquid τ variables. As we do not know beforehand whether there is a VLLE region and where it is eventually located, we have to solve with each step of integration a modified τ method to make sure that one is not near or in a VLLE region until a positive pseudo phase τ variable value is found. As in [START_REF] Han | A method for multiphase equilibrium calculations[END_REF]. an Iterative Linear Programming method to solve the τ method. The ILP method can handle inequality constraints and uses Newton Raphson method concepts as it is quadratically convergent towards the solution (Llovell 2002).

The modified τ method set of equations that we solve to check any occurrence of a non trivial solution bears some modifications compared to the τ method equations written in the first section: for the purpose of residue curve calculations, the vapor phase is supposed to always exists. The initial mixture temperature is then set at its boiling point. The equations become:

( )

II n I n min τ + τ (33) II n I n Total n U U U + = (34) nc ..., , 1 i x U x U x U II i , n I n I i , n I n Total i Total n = ⋅ + ⋅ = ⋅ (35) ( ) nc ,..., 1 i x y , x , T , P K y I i , n i , n I i , n n I i , n i , n = ⋅ = (36) ( ) nc ,..., 1 i x y , x , T , P K y II i , n i , n II i , n n II i , n i , n = ⋅ = (37) ∑ = = τ + nc 1 i I n I i , n 1 x (38) ∑ = = τ + nc 1 i II n I i , n 1 x (39) ∑ = = nc 1 i i , n 1 y (40)
with the constraints:

II I n n U , U 0 ≤ (41) 1 , , y , x , x 0 II n I n i , n II i , n I i , n ≤ τ τ ≤ ( 42
)
The combination of the modified τ method VLL flash calculation and of the dynamic positive flash formulation tempers our initial ambition of developing a single model suitable for any phase regions.

But it performs quite well and consists of an efficient alternative to phase stability VLL equilibrium models usually implemented in distillation column dynamic simulators: the positive flash model handles any phase pattern change as long as it is correctly initialized by the modified τ method.

To illustrate these concepts, scenario 2 is devised for a cyclohexane -water -ethanol mixture. Figure 6 depicts its composition movements on the ternary diagram and Figure 7 shows the evolution of key variables during the scenario. At first, the initial ternary mixture is located in the liquid -vapor region of the cyclohexane -water -ethanol diagram. Then the mixture is cooled and its composition follows a residue curve in the opposite direction towards an unstable node. At some arbitrary point in the vapor -liquid -liquid region, cooling is stopped and water is added to shift the composition. Then the mixture is heated and its composition follows a residue curve ending at the pure water stable node.

At each integration step of the positive flash model until a VLLE region is encountered, the modified τ method is solved to check whether any pseudo phase τ variable has a positive value. From the initial mixture composition of scenario 2, both modified τ method and positive flash converge at the trivial solution until the limit of trivial solution attractiveness region depicted in figure 6 is reached. At this point, the modified τ method finds one positive τ variable (figure 7). This value is then used to initialize the positive flash model solving and convergence is obtained at a non trivial solution. Once a positive τ value is found, the modified τ method is no longer solved and integration of the positive flash model proceeds normally: at the vapor -liquid -liquid boiling envelope, a VL -VLL transition is triggered like in scenario 1. Integration proceeds into the VLLE region until it is arbitrarily stopped to add water. The resulting mixture is then heated and integration resumes normally: the composition follows a residue curve towards and across the VLLE envelope boundary and the VLL -VL transition is triggered at the VLLE envelope. Again, at some point located on the approximate dotted line (the limit of the trivial solution attractiveness region) on figure 6, the number of iterations required to solve the positive flash model increases exponentially (up to 1000 iterations) until it converges at the trivial solution thanks to the internal reinitialization procedure of the DAE solver DISCo (figure 7). Finally, integration proceeds correctly along the residue curve in the VL region until the mixture composition reaches the water vertex.

Scenario 2 has illustrated the initialization procedure when residue curves are integrated in the reverse way (cooling of the Rayleigh kettle). In the direct integration procedure, the same initialization procedure is implemented as it may happen that a residue curve starting in a VLE region enters a VLLE region for ternary diagrams where a binary heteroazeotrope is a saddle point (e.g. the wateracetonitrile -acrylonitrile diagram shown in figure 8).

Residue curve maps

A residue curve is a trajectory of the liquid composition versus time left in a heated kettle. This is equivalent to a simple distillation process (Rayleigh distillation -figure 3). Residue curves maps are very useful for the feasibility analysis of distillation processes: they enable in particular to display azeotropes and distillation boundaries that make some separations impossible and to hint at column liquid composition profiles in continuous packed distillation towers (Widagdo andSeider, 1996, Kiva et al., 2003). Homogeneous residue curve maps are easy to compute, but residue curve maps for heterogeneous systems that exhibit a vapor -liquid -liquid region are more complicated because the dynamic set of equation must handle phase number changes along the calculation of some residue curves and distillation boundaries are not trivial to compute. A usual solution is to solve a VL flash and run a liquid -liquid stability check at every integration step. When a liquid -liquid split is detected, one shifts to a VLL set of equations. Instead, the "positive flash model" presented here enables to keep a single set of equations along the calculation of a residue curve. For both homogeneous and heterogeneous residue curve maps, starting from an initial mixture composition at boiling temperature, the dynamic set of equations of the positive flash model is integrated in two opposite directions to calculate a whole residue curve. Direct (resp. inverse) integration corresponds to heating (resp.

cooling) the Rayleigh kettle and removes the light (resp. heavy) component and the trajectory moves towards the stable (resp. unstable) node of the basic distillation region. Unlike the heating of the Rayleigh kettle, which corresponds to evaporation, cooling has no physical meaning as it condensates some vapor overhead.

The "positive flash model" set of equations is kept all along the calculation. As stated previously, unless a VLLE region has been encountered, the modified τ method is also solved to detect any non trivial solution. As this happens, the liquid τ values then found are used to reinitialize the positive flash model. Then the τ method is no longer used and the integration resumes. For the purpose of residue curve calculation, the vapor phase existence is assumed beforehand. Then τ V = 0 and the first call to the modified τ method finds a temperature value equal to the boiling temperature.

For practical purposes, the residue curve computation is stopped when the equality of vapor and liquid compositions in equilibrium is obtained with a given tolerance. In a composition diagram, this condition is satisfied for all singular points (vertexes and azeotropes), and in particular for the unstable and stable nodes that are respectively the starting point and the ending point of residue curves. But it also applies to saddle points to which residue curve may come very close. In order to avoid an unexpected stop at these points, the convergence tolerance is lowered whenever a unary or a binary saddle point is neared during integration.

In addition to residue curves, residue curve maps may display other key features that are the knowledged of all singular points (pure components and azeotropes), their stability within a distillation region (stable -highest boiling point of the region ; unstable -lowest boiling point; saddleintermediate boiling point) and distillation boundaries existence and curvature [START_REF] Kiva | Azeotropic Phase Equilibrium Diagram: A Survey[END_REF].

Including all singular points, the topology equation [START_REF] Zharov | Physicochemical fundamentals of distillations and rectifications[END_REF][START_REF] Matsuyama | Topological and Thermodynamic Classification of Ternary Vapor -Liquid Equilibrium[END_REF][START_REF] Doherty | On the dynamics of distillation processes : III. Topological Structure of Ternary Residue Curve Maps[END_REF][START_REF] Widagdo | Azeotropic Distillation[END_REF][START_REF] Kiva | Azeotropic Phase Equilibrium Diagram: A Survey[END_REF]) must be verified. The topology equation for ternary mixtures is:

2•N 3 -2•S 3 + N 2 -S 2 + N 1 =2 (38)
Where N 3 is the number of ternary node, S 3 the number of ternary saddle points, N 2 the number of binary nodes, S 2 the number of binary saddle points and N 1 the number of pure components nodes.

Methods for the computation of all azeotropes are out of scope of this paper [START_REF] Widagdo | Azeotropic Distillation[END_REF][START_REF] Kiva | Azeotropic Phase Equilibrium Diagram: A Survey[END_REF]. In our cases, existing binary azeotropes are first computed, from residue curve calculation on binary mixtures. Their existence is evaluated from the simple following algorithm:

For each binary mixture A-B do { Evaluate quickly possible azeotropy { -Component A being the most volatile, compute

A A A y x K = at 0 x A ≈ and at 1 x A ≈ -if 1 K 0 x A A > ≈ and 1 K 1 x A A < ≈
then there exists a minimum boiling temperature azeotrope

-if 1 K 0 x A A < ≈ and 1 K 1 x A A > ≈
then there exists a maximum boiling temperature azeotrope -else no binary azeotrope exists } Calculate accurately azeotrope composition, type and stability { -compute a residue curve (direct integration for maximum boiling temperature azeotrope and inverse integration for minimum boiling temperature azeotrope). The final point is the binary azeotrope.

} }

Then unary and binary stability is obtained from an eigenvalue analysis related to the jacobian associated to the dynamic model (Doherty andPerkins, 1978, 1979). Knowing them, the topology equation tells us whether there exists any ternary saddle azeotrope or any ternary unstable azeotrope.

If so, the saddle composition and temperature is found from a procedure devised by Doherty andPerkins (1978, 1979) any consisting in successive maximization and minimization along the boiling temperature surface to follow crests and valley towards the saddle points. To find an eventual unstable azeotrope, a residue curve is computed using an inverse integration scheme from an initial composition involving all three components. Finally, as distillation boundaries have the property of being residue curves, they are computed using information provided by the ternary diagram classification of serafimov [START_REF] Kiva | Azeotropic Phase Equilibrium Diagram: A Survey[END_REF] as this classification readily indicates from which singular point should a residue curve calculation be started to give a distillation boundary. This methodology is applied to the automatic computation of several residue curve maps shown on figure 8 with thermodynamic parameters displayed in tables 1 to 3. VLLE envelope is computed externally through VLLE flash calculation using BibPhyAddIn in an Excel spreadsheet (Prosim, 2000) using the same thermodynamic model and set of binary parameters. VLLE envelope results fully agree with the information on the VLLE envelope found during the integration: compositions at the phase pattern change events and compositions of the two coexisting liquid phases in the VLLE region.

In figure 8a, the well known cyclohexane -water -ethanol diagram is displayed with an unstable ternary heteroazeotrope. Figure 8b shows the chloroform -acetone -methanol diagram which saddle ternary homogeneous azeotrope if found solving the topology equation. Its composition and temperature is the mean of the distillation boundaries end composition computed with the dichotomy procedure explained above. The saddle point temperature and composition maximum deviation is less than 0.1 °C and 0.03 in molar fraction. Such accuracy is far better than for values obtained from flash calculations and boiling temperature maps.

CONCLUSION

A general and polyvalent model for the dynamic simulation of a vapor, liquid, liquid-liquid, vapor-liquid or vapor-liquid-liquid tray is proposed: the positive flash model. It handles any phase pattern changes that may occur in the dynamic simulation of equilibrium distillation columns or during the computation of heterogeneous residue curve maps. Whatever the number of phases, the set of equation and the number of variables is kept unchanged, thanks to pseudo phase variables τ that take a positive value when a phase is not physically existent. Phase number changes are considered as events described by switch equations involving the variables τ. Nevertheless, practical implementation requires using a modification of the τ -method introduced as a minimization problem by Han & Rangaiah (2000) for steady-state simulation. It enables to find initial values of the pseudo phase variables τ suitable for the positive flash model. Validation of the model is showed through testing of its capacity to overcome phase number changes and to compute heterogeneous ternary residue curves maps. Further application will concern the building of a multipurpose distillation column dynamic simulator. [START_REF] Gmehling | Azeotropic Data[END_REF] Table 3. Thermodynamic model, binary parameters and calculated azeotropic data for the acrylonitrile -acetonitrile -water system. V (mol/s)

NOMENCLATURE

FIGURE CAPTION

L 1 L 2 V L 1 L 2 L 1 L 2 L 1 L 1 L 1 V

Figure 1 .FigureFigure 1 .

 11 Figure 1.Adiabatic vapor -liquid -liquid equilibrium stage for the τ-method

Figure 5 .Figure 7 .

 57 Figure 5. Evolution of τ variables ( τ I ▲ τ II τ V ) and holdup, vapor flow rate and temperature during simulation of scenario 1 ( U I , ○ U II ,

TABLE CAPTION Table 1 .

 CAPTION1 Thermodynamic model, binary parameters and calculated azeotropic data for the cyclohexane -water -ethanol system. dependent NRTL thermodynamic model: A ij,0 + A ij,T .T ; A ji,0 + A ji,T .T ; α ij,0 + α ij,T .T

	Cyclohexane (Cy) -Water (W) -Ethanol (E) at 1 atm.
	Binary interaction parameters A ij,0	A ji,0	α ij,0	A ij,T .	A ji,T	α ij,T
	Cy-W		2835	3125	0.274	0	0	0
	Cy-E		1390.42	876.793	0.4485	0	0	0
	Table 2. W-E	Thermodynamic model, binary parameters and calculated azeotropic data for the 1616.8 635.56 0.1448 2.0177 0.9907 0
		chloroform -acetone -methanol system. molar fraction		T (°C)	Relative
	Table 3. azeotrope node type Thermodynamic model, binary parameters and calculated azeotropic data for the calculated experimental * error vs. calc. exp.* exp
	Cy-W	acrylonitrile -acetonitrile -water system. saddle heterogeneous 0.701 / 0.299 / 0.000 0.700 / 0.300 / 0.000 69.4	69.5	0.1%
	Cy-E	saddle homogeneous	0.548 / 0.000 / 0.452 0.545 / 0.000 / 0.455 64.7	64.8	0.4%
	W-E	saddle homogeneous	0.000 / 0.116 / 0.884 0.000 / 0.107 / 0.893 78.1	78.2	1.2%
	Cy-W-E	unstable heterogeneous	0.158 / 0.536 / 0.306 0.167 / 0.527 / 0.306 62.9	62.6	1.3%
	* (Gmehling et al., 1994)

temperature

Table 1 .

 1 Thermodynamic model, binary parameters and calculated azeotropic data for the cyclohexane -water -ethanol system.

	Chloroform (C) -Acetone (A) -Methanol (M) at 1 atm.
	NRTL thermodynamic model: A ij ; A ji ; α ij
	Binary interaction parameters	A ij,0	A ji,0	α ij,0
	C-A		228.457	643.277	0.3043
	C-M		2736.860	-1244.030	0.0950
	A-M		184.701	222.645	0.3084
				molar fraction	T (°C)	Relative
	azeotrope node type	calculated	experimental *	calc. exp.*	error vs. exp
	C-A	stable homogeneous	0.646 / 0.354 / 0.000 0.659 / 0.341 / 0.000 65.1	64.4	1.8%
	C-M	unstable homogeneous	0.654 / 0.000 / 0.346 0.650 / 0.000 / 0.350 53.3	53.4	0.6%
	A-M	unstable homogeneous	0.000 / 0.787 / 0.213 0.000 / 0.794 / 0.206 55.4	55.5	1.0%
	C-A-M	saddle homogeneous	0.207 / 0.373 / 0.420 0.215 / 0.356 / 0.429 57.0	57.5	2.1%
	* (Gmehling et al., 1994)

Table 2 .

 2 Thermodynamic model, binary parameters and calculated azeotropic data for the chloroform -acetone -methanol system.

	Acrylonitrile (Acr) -Acetonitrile (Ace) -Water (W) at 1 atm.
	NRTL thermodynamic model: A ij ; A ji ; α ij
	Binary interaction parameters	A ij,0	A ji,0	α ij,0
	Acr-Ace			-336.251	475.621	0.3042
	Acr-W			584.640	2122.310	0.2960
	Ace-W			364.836	1321.730	0.2858
				molar fraction	T (°C)	Relative
	azeotrope node type	calculated	experimental *	calc. exp.*	error vs. exp
	Acr-W	unstable heterogeneous	0.684 / 0.000 / 0.316 0.698 / 0.000 / 0.302 71.2	71.0	2.0%
	Ace-W	saddle homogeneous	0.000 / 0.675 / 0.325 0.000 / 0.690 / 0.310 76.8	76.5	2.1%
	* (


Pierre-Olivier Mazagol
email: pierre.olivier.mazagol@univ-st-etienne.fr

P Niogret

Jérémie Riquier

Michel Depeyre

R Ratajczak

Carlos F Crispim-Junior

L Tougne

P.-O Mazagol

Tools Against Oblivion: 3D Visualization of Sunken Landscapes and Cultural Heritages Applied to a Dam Reservoir in the Gorges de la Loire (France)

Keywords: 3D, geovisualization, Geographic Information System, sunken landscape, cultural heritage, dam reservoirs

published or not. The documents may come

Introduction

Dams and associated reservoirs have multiple socioeconomic functions, such as generating renewable energy production, supplying water to people and crops, and facilitating tourism and leisure activities [START_REF] Marnezy | Les barrages alpins. De l'énergie hydraulique à la neige de culture[END_REF][START_REF] Loloum | La vie touristique des grands barrages hydroélectriques[END_REF]. Dams can also have positive health impacts [START_REF] Lerer | Health impacts of large dams[END_REF]. However, dams are usually associated with numerous adverse ecological and biophysical effects on the environment, including threats to connectivity for freshwater fish (e.g. [START_REF] Barbarossa | Impacts of current and future large dams on the geographic range connectivity of freshwater fish worldwide[END_REF], alteration of flows, sediment transfer and transport leading to vegetation encroachment, riverbed incision and narrowing, and in-stream habitat degradation (e.g. [START_REF] Petts | Dams and geomorphology: Research progress and future directions[END_REF].

Since the 1950s, the number of large hydroelectric dams, that is, those with heights of 15 meters or more from lowest foundation to crest or a dams between 5 meters and 15 meters impounding more than 3 million cubic meters (ICOLD 2011), has rapidly increased and can be seen as a symbol of the Anthropocene [START_REF] Steffen | The trajectory of the Anthropocene: The Great Acceleration[END_REF]. Because large dams and associated reservoirs cover larger land surface areas, they are more likely to include historical and cultural heritages. Several case studies dealing with reservoir impoundment losses have been reported in the literature. In 1970, Egypt's manmade Lake Nasser, created by the Aswan Dam on the Nile River, flooded archeological sites including the Buhen fort and the cemetery of Fadrus, while several other temples that were physically moved were preserved [START_REF] Harris | Orthodontics' contribution to UNESCO's campaign to save the monuments of Nubia: A 1970 field report[END_REF]. In 2003, China's Three Gorges Dam, on the Yangtze River, sank two large cities, 11 towns, 114 townships, numerous villages and farms [START_REF] Jackson | Resettlement for China's Three Gorges Dam: socio-economic impact and institutional tensions[END_REF] and over 1280 archaeological and cultural heritage sites [START_REF] Ponseti | The Three Gorges Dam Project in China: history and consequences[END_REF]. In 2020, the Turkish Ilısu Dam on the Tigris River will sink the ancient town of Hasankeyf and its 12 000 years of history, as well as 180 other villages [START_REF] Drazewska | Hasankeyf, the Ilisu Dam, and the Existence of "Common European Standards" on Cultural Heritage Protection[END_REF]. These dams are among the most emblematic examples worldwide.

In France, after the Second World War, the state embarked on a major policy of building large dams for the stated purpose of electricity production (Barjot, 2013). In his book "La France des villages engloutis" (The France of sunken villages), G. Guérit (2019) listed 44 valleys with sunken villages due to reservoir impoundment between 1914 and 1992. Dam construction inevitably entails different degrees of trauma. Difficulties experienced by the local population include the expropriation of land, resettlement of the population, cemetery relocation, police interventions, and residence destruction. Regardless of the location, the trauma of these events is constant and often generates social mobilization [START_REF] Schapper | Social mobilization against large hydroelectric dams: A comparison of Ethiopia, Brazil, and Panama[END_REF][START_REF] Flaminio | Representations of dams: infrastructures, concessions and controversies[END_REF]. This phenomenon was observed in, for example, the hydroelectric schemes of Tignes (in 1953) and Serre-Ponçon (in 1960), where a classic opposition between the "best interest of the nation" and "individual interests" was observed [START_REF] Bodon | La Modernité au village -Tignes, Savines, Ubaye… La submersion de communes rurales au nom de l'intérêt général. 1920-1970[END_REF]. In other cases, projects were more easily accepted by the local population and the situation remained calm, as with the Grangent dam in the Gorges de la Loire (1957).

The construction of such dams is often associated with the disappearance of anthropogenic infrastructures (e.g. bridges, houses, and railroads), socioeconomic activities, and cultural and landscape heritages under these flooded valleys. These landscape features remain mostly invisible, except in instances in which the water level decreases significantly. This invisibility has both a temporal and a spatial dimension. The invisibility is temporal because all of these underwater infrastructures and landscapes belong to the past and do not exist in the present, and it is spatial because these areas are flooded in the present. However, these heritages remain alive in the memories of local inhabitants, for example, through visits during reservoir draining periods or during dives. These legacies are now of interest to a growing number of people. Such renewed interest is also sustained by legends regarding sunken cities, where Ys (Brittany), Kitezh (Russia), Atlantis, Lemuria, and the Great Flood of Gun-Yu (China) are among the best known. This interest raises the question of how to keep these lost heritages "afloat" and accessible to all. To meet such an objective, several digital methods are generally used in heritage studies involving the preservation or management of specific sites, including remote sensing and GIS [START_REF] Summerby-Murray | Analysing Heritage Landscapes with Historical GIS: Contributions from problembased inquiry and constructivist pedagogy[END_REF][START_REF] Agapiou | Cultural heritage management and monitoring using remote sensing data and GIS: The case study of Paphos area, Cyprus[END_REF][START_REF] Domínguez-Ruíz | Contribution to the Knowledge of Cultural Heritage via a Heritage Information System (HIS). The Case of "La Cultura del Agua" in Valverde de Burguillos[END_REF][START_REF] Elfadaly | Cultural Heritage Management Using Remote Sensing Data and GIS Techniques around the Archaeological Area of Ancient Jeddah in Jeddah City, Saudi Arabia[END_REF][START_REF] Griffiths | Mapping spatial cultures: contributions of space syntax to research in the urban history of the nineteenth-century city[END_REF], 3D GIS [START_REF] Campanaro | 3D GIS for cultural heritage restoration: A 'white box' workflow[END_REF][START_REF] Dell'unto | Experiencing Ancient Buildings from a 3D GIS Perspective: a Case Drawn from the Swedish Pompeii Project[END_REF][START_REF] Deguy | Ricostruire e comunicare il paesaggio storico[END_REF][START_REF] Ferdani | 3D reconstruction and validation of historical background for immersive VR applications and games: The case study of the Forum of Augustus in Rome[END_REF], (H)BIM [START_REF] Colucci | HBIM-GIS Integration: From IFC to CityGML Standard for Damaged Cultural Heritage in a Multiscale 3D GIS[END_REF][START_REF] Raner | Regularized Reconstruction of HBIM for Built Heritage: Case Study with Chinese Ancient Architecture[END_REF]), 3D modeling and virtual landscaping [START_REF] Jacquinod | Production, pratique et usages des géovisualisations 3D dans l'aménagement du territoire[END_REF][START_REF] Berti | Disegni d'acqua. Acque e trasformazioni del territorio[END_REF], and more generally, 3D digital technologies [START_REF] Jacquinod | Production, pratique et usages des géovisualisations 3D dans l'aménagement du territoire[END_REF][START_REF] López | An efficient GPU approach for designing 3D cultural heritage information systems[END_REF]). These tools and technologies are not mutually exclusive; they can be used on different types of objects and applied to a wide range of spatial scales (e.g., vessels and such buildings as tombs, temples, castles, cities, and landscapes). Tools should be chosen according to the type of heritage being documented, the scale, and research or management objectives. Each tool involves a relevant methodology. In this paper, we propose a 3D GIS methodology associated with 3D geovisualization to recreate the sunken landscape of the Gorges de la Loire (Loire gorges).

Methodology

Study area

The study site under investigation is located in France, near the city of Saint-Etienne, in the Gorges de la Loire (Figure 1). Following a construction period from 1955 to 1957, the impoundment of the Grangent reservoir profoundly and abruptly changed the socioecological landscape of the valley between the cities of Aurec-sur-Loire and Saint-Just-Saint-Rambert (Figure 1).

Figure 1. Location of the study area in France

This valley has been used for centuries and remains in the collective memory. Before the establishment of the Grangent dam, there were villages, farms, factories, roads, and a canal, as well as a railroad and related infrastructure, such as bridges and tunnels. All these landscape elements symbolized the history of an industrial region through the 1950s. Furthermore, the river had an important social role for inhabitants of the surrounding towns, including Saint-Etienne. These residents appreciated relaxing on the banks of the Loire, dancing and drinking wine in the open-air dance halls, swimming, or attending water jousting events (Peyraud 2006). During the postwar economic boom and in the context of reconstruction and economic growth, the construction of a hydroelectric dam began a radical transformation of this area. The landscape underwent a radical change as a 365 ha and 21 km long lake with a total reservoir capacity of 57 million cubic meters, which submerged the Loire Valley. The change symbolized modern life in the second half of the 20th century abruptly replacing the former landscape of industrialization. This new landscape is dedicated to electricity production but also to leisure activities and tourism with a marina, a nautical base and a beach. However, during partial decennial emptying for maintenance purposes, the level of the reservoir drops from approximately 420 m asl to 405 m asl, and indicators of the past temporarily reappear as ruins and other remains resurface. During these events, the people of the Saint-Étienne region visit the remnants of the past, use bridges, and enter former tunnels. Such an infatuation reflects their interest in this area and the need to reappropriate these historically significant locations. The methodology proposed in this study constitutes an alternative method of providing access to sunken heritages to inhabitants and local authorities.

Data collection

The first step of our work is based on relevant archival research on valley memory, including socioeconomic activities, villages, hamlets, water mills, isolated houses, ownership information, bridges, and railways (Peyrot, 2006). Other archival records have been collated to complete this dataset [START_REF] Niogret | Géovisualisation 3D du patrimoine englouti des Gorges de la Loire[END_REF], including land registers, historical iconography, maps, newspaper articles, and interviews with associations and inhabitants.

Georeferencing

We created a spatial database with a 3D component and then constructed a 3D geovisualization using the projected coordinate system for France: RGF93/Lambert 93 -EPSG: 2154. To this end, we used ArcGIS Pro (ESRI), mainly because it provides an all-in-one solution. We used different sources of referential spatial data. Georeferenced and orthorectified aerial photographs of France from 1950 to 1965, provided by the French National Institute of Geographic and Forest Information (IGN), provide national coverage that is heterogeneous in terms of available dates and photograph quality. In our area of interest, the pre-dam orthophotographs date from 1953 but are unfortunately too blurred to be used for geovisualization. The IGN also provides several non-georeferenced aerial photograph series from different periods (https://remonterletemps.ign.fr/). We selected a total of 68 images dating from 1955 to cover the study area which encompasses approximately 120 km². These images were georeferenced from 1953 images. We chose approximately fifty ground control points per image (scale: 1/28978) evenly distributed over the scene within the area of interest and used the "adjust" transformation. The RMS error is near 0. We obtain 68 warped images that need to be cropped to obtain a mosaic of 75 regular raster images. We created a color-balanced mosaic dataset and ensured that transitions from one image to the adjoining image appear seamless (Figure 2).

Automatic colorization of grayscale aerial images

As we needed to obtain the best possible visual reproduction of the pre-dammed landscape, we chose to colorize these grayscale images. This step is the only one performed without ArcGIS Pro. We adapted the method developed by [START_REF] Ratajczak | Toward An Unsupervised Colorization Framework for Historical Land Use Classification[END_REF] to automatically colorize historical aerial images (Figure 3). This method is based on a relatively small cyclic generative adversarial neural network, inspired by [START_REF] Zhu | Unpaired Image-to-Image Translation using Cycle-Consistent Adversarial Networks[END_REF], that the authors designed to handle relatively large images (1024x1024 image pixels), considering hardware constraints. The strength of this method is in its ability to learn colorization in an unsupervised way, eliminating the need for panchromatic-color image pairs that could be difficult to obtain for historical data. In our experiments, we reused the model trained by the authors after 60 epochs, whose colorization results were evaluated through a Mean Opinion Score study. However, while [START_REF] Ratajczak | Toward An Unsupervised Colorization Framework for Historical Land Use Classification[END_REF] proposed working patchwise to keep the images at full resolution, we noticed that the color bands usually have a lower resolution than the panchromatic band. Inspired by classical pansharpening methods, and in order to make the colorization process faster while avoiding undesired mosaicking effects, we experimented with applying the colorization method on downscaled panchromatic images (from original size to 1024x1024 image pixels) instead of applying it patchwise. Generated downscaled color bands were further upscaled before being merged with the original panchromatic image. Next, we created another mosaic dataset from colorized photographs and added it to the geodatabase.

Digitization

We digitized different elements of the past landscape into several thematic vector layers and edited attribute tables from the archive work for roads, streets, paths, ferries, railways, houses, factories, mills, castles, canals, bridges, and tunnels. We chose to digitize only elements within a 300-m buffer zone along the centerline of the river, which already represents a significant amount of work. Buildings have required special attention. Two approaches were implemented depending on whether the buildings were (1) in the buffer zone or (2) outside the buffer zone.

1. The buildings in the buffer zone, which are now submerged, destroyed or still lakeside, have cultural, historical and heritage significance. The attribute data linked to these features are as exhaustive as possible, and include construction date, demolition date, current condition, type of building (e.g. farm, factory, and house), owner, cadastral parcel number, and height (sometimes estimated from period photographs).

2. For buildings outside of the buffer zone, we chose to downgrade a pre-existing layer provided by the IGN (BD TOPO®). This approach implies correction of the height, removal of buildings constructed since 1955 and addition of those that have been demolished.

DEM creation

The digital elevation model (DEM) is a crucial element of our work. The IGN currently provides a 5-m resolution DEM (RGE ALTI®, version 2.0). However, the IGN represents the surface of the lake and not the sunken relief.

To complete these data, we used a 1-m-resolution bathymetric DEM derived from multibeam echosounder measurements performed in the lake in 2010 by Électricité de France (EDF), the company managing the dam. The EDF DEM, representing the permanently lands, was resampled to match the 5-m resolution of the IGN DEM. Next, using a raster calculator, we replaced IGN DEM values with bathymetric ones to create a new DEM representative of the relief in 1955. At this stage, some artifacts remained, primarily due to the bridge ruins (Figure 4). An analytical hillshading calculation applied on the new DEM, combined with the bridge vector layer created, enabled us to easily locate such artifacts. Then, we removed artifacts from the DEM and applied a nearest neighbor interpolation using the Nibble function of ArcGIS Pro to replace cells. It was then possible to transform 2D vector layers to 3D vector layers from this final DEM.

3D Modeling

This step involved integrating the colorized aerial photographs draped over the DEM representation of the 1955 relief and the 3D vector layers into a 3D scene. To obtain a visual appearance that was as authentic as possible, we added three principal types of 3D objects: multipatch features, 3D objects from non-GIS software and extruded polygons. The Esri multipatch data format is a geometry used as a boundary representation for 3D objects. This format provides levels of detail ranging from a simple cube to complex 3D models. We used multipatch features to model complex details of buildings, such as roofs, chimneys, windows and doors. We also added wall textures derived from colorized period postcards or photographs when available (Figure 5). Other 3D objects were created and imported from Sketchup software, such as bridges or important buildings, such as the "castle of Essalois" (Figure 6), and even a steam train. For other buildings, we opted for a simple representation with extruded polygons (Figure 6). For comparison purposes, we also created a basic 3D scene for 2016 based on available IGN databases (Figure 9) (buildings from the BD TOPO®, DEM from the BD ALTI® and aerial photographs from the BD ORTHO®).

Results and Discussion

Interactive 3D geovisualization

Our methodology provides a realistic 3D geovisualization of the past Gorges de la Loire landscape (Figure 7). The 3D scene (https://arcg.is/0LOeTH) and the 2D map (https://arcg.is/1jqam4) are available on the web. This visualization is dynamic and interactive. Indeed, users can navigate easily into this 3D scene, interact with all objects by clicking on them (information included in attribute tables is displayed in a pop-up window, Figure 8), or choose what layers they want to activate and see in the 3D scene or the 2D map. Using an ArcGIS online account, it is possible to share 3D scenes and 2D maps and to provide access to anyone through a web browser.

The last way to share the modeling results is video exports (https://mediacenter3.univ-stetienne.fr/videos/?video=MEDIA200930111345487) or screenshot.

Memories

This interactive 3D geovisualization can be considered a snapshot, at a time when what can be referred to as a bifurcation (sensu Grossetti et al., 2010) occurs, i.e. a radical transformation. In fact, in our case, we observe a quadruple bifurcation. The first marked the turning point in the French energy policy after the Second World War.

The second induced a drastic upheaval in the trajectory of a territory whose socioeconomic activity is overturned. The third transfigured a landscape considered relatively stable over time. Finally, the fourth bifurcation is those of memories. We identified two types of memories in the present study. These types of memories are an individual and collective memory of material and immaterial heritages that must be preserved, even when it has already disappeared, and memory as data storage. Regarding this latter type of memory, the creation of a spatial database can be considered a "backup copy" tool for heritages. It seems important to implement such an approach to create a comprehensive territorial 3D GIS knowledge database before the bifurcation event through, for example, the geolocalization of each component of the landscapes or taking photographs. After the bifurcation, the 3D GIS allows safeguarding of the original landscapes and cultural heritages in memory (Figure 9). Without anticipation, the only way to proceed is to conduct archive work as performed in this study, which does not seem to be the best way. Indeed, when the work is done a posteriori, information retrieval can be more difficult, often partial, and sometimes impossible. These two types of memories previously mentioned are linked. The sunken memory (i.e. material and immaterial heritages) could "resurfaced" through the "backup memory".

Different purposes of use

This type of 3D geovisualization turns out to be a relevant backup tool for landscape and cultural heritage, not only in the context of dams and their reservoirs, but more generally in anticipation of natural hazards, such as earthquakes or armed conflicts. It could be an interesting tool for the management of these heritage sites, to publicize them and to promote the reappropriation of these lost places. In this regard, our modeling has been requested by two local authorities to be presented to the visitors in two heritage interpretation centers: that of the Gorges de la Loire and that of the Saint-Étienne metropolis. Finally, such tools could be of great interest to social scientists, as historians or sociologists, for example, who are not always accustomed to using digital technologies. These approaches allow further analysis of the lost components of a territory, enabling a new point of view that does not always come out through simple analysis of archival records. Indeed, our methodology requires significant research in archival records and allows the creation of important knowledge databases. This analysis makes such knowledge visible, as it was previously hidden in documents or in inhabitants' memories. Second, modeling is a way to gain a better understanding of past territories' organization and to detect certain relations that were not visible simply by reading archival records. We can also present the past landscape to the inhabitants, regardless of whether they are familiar with it. For example, providing the public with a comparison of the predam (1955) and the post-dam (2016) 3D modeling could serve as baseline data to address their perception and values of the pre-damming landscape among other potential perception studies.

Limitations and improvement opportunities

The 3D geovisualization presented in this study is the first step in a longer-term process. The primary limitations of this work have been identified. First, a large part of the proposed workflow is based on a proprietary solution whose cost could be an obstacle to reproducibility and generalization. To date, we have made ArcGIS Pro our first choice, given the relative simplicity of implementation and the quality of results obtained with regard to our original needs, objectives and constraints. Open-source GIS software, such as QGIS, will clearly become a good alternative solution in the coming years. Currently, none of them is able to obtain similar rendering (e.g. in terms of color balanced mosaic datasets or 3D modeling and rendering) or to facilitate sharing of 3D scenes on the web. Another important aspect that should not be overlooked is the amount of time needed to create complete 3D modeling on such a large scale.

There are many ways to improve the work performed. In the more global context of our project, 3D geovisualization is a first step towards more elaborate technologies, such as augmented reality or virtual reality, which we will test soon. Indeed, virtual reality interaction with this landscape would be a way to give it back its humanity. The integration of sounds (e.g. active factories, bathing areas, and birds) may better reflect this lost landscape and enhance the user experience, as experimented in the Bretez project (Pardoen 2019), whose objective is patrimonial valorization through 3D digital restoration and spatial sound restoration.

Conclusions

In this study, we proposed an innovative and interdisciplinary application of 3D mapping. More specifically, we developed an original approach for the 3D reconstruction of a landscape engulfed by a dam reservoir, as well as reconstructing its cultural heritage in the context of memory preservation. To this end, we notably employed an original automatic colorization method of old grayscale aerial photographs, which considerably improves the visual rendering of the modeled scene. The combination of archive work and geomatic tools enabled us to make visible again a territory that became doubly invisible as a result of time and human action. Our approach can also be viewed as a means of making visible a bifurcation of history in a small territory. This approach is a relevant solution that gives back to the population what it no longer has, but it still feels like a phantom limb and cannot be given back any other way. Crossing disciplinary perspectives makes it possible to take a different look at the same objects and to enhance the user experience.

Figure 2 .

 2 Figure 2. Georeferenced images from 1955 (a) before and (b) after balancing the mosaicked image

Figure 3 .

 3 Figure 3. Comparison between grayscale (top) and colorized (bottom) aerial photographs of 1955

Figure 4 .

 4 Figure 4. Bathymetric artifact removal from the DEM: (a) bridge ruins remain and (b) the final DEM

Figure 5 .Figure 6 .

 56 Figure 5. Multipatch feature with texture from colorized period postcard (Hôtel des Bains, Destroyed and sunk in 1957)

Declarations

This work was supported by the LABEX IMU (ANR-10-LABX-0088) of Université de Lyon within the program "Investissements d'Avenir" (ANR-11-IDEX-0007) operated by the French National Research Agency (ANR).


Sylvain Tourancheau
email: sylvain.tourancheau@univ-nantes.frkjell.brunnstrom@acreo.se

Patrick Le Callet

Kjell Brunnström

Börje Andrén

Psychophysical study of LCD motion-blur perception

Motion-blur is still an important issue on liquid crystal displays (LCD). In the last years, efforts have been done in the characterization and the measurement of this artifact. These methods permit to picture the blurred profile of a moving edge, according to the scrolling speed and to the gray-to-gray transition considered. However, other aspects should be taken in account in order to understand the way LCD motion-blur is perceived.

In the last years, a couple of works have adressed the problem of LCD motion-blur perception, but only few speeds and transitions have been tested. In this paper, we have explored motion-blur perception over 20 gray-to-gray transitions and several scrolling speeds. Moreover, we have used three different displays, to explore the influence of the luminance range as well as the blur shape on the motion-blur perception.

A blur matching experiment has been set up to obtain the relation between objective measurements and perception. In this experiment, observers must adjust a stationary test blur (simulated from measurements) until it matches their perception of the blur occuring on a moving edge. Result shows that the adjusted perceived blur is always lower than the objective measured blur. This effect is greater for low contrast edges than for high contrast edges. This could be related to the motion sharpening phenomenon.

INTRODUCTION

Though LCD (liquid crystal display) has become the market leader of monitors and televisions in the last years, motion-blur is still an important issue on these type of displays. This artefact is due to the sample-and-hold behaviour of LCD: the light intensity is held on the screen for the duration of the frame (LC displays are so called hold-type displays). This is particularly annoying when observers are tracking a moving object on LCD: the object is held stationary on the screen during the frame period while observers'eyes are still moving slightly (smooth pursuit). The edges of the object are then displaced on the retina resulting in a blur. [START_REF] Kurita | Moving picture quality improvement for hold-type AM-LCDs[END_REF] In the last years, efforts have been done in the characterization and the measurement of this artifact. One common measurement system used to characterize LCD motion-blur is referred as MPRT measurement system (for Moving Picture Response Time). [START_REF] Oka | Moving picture response time (MPRT) measurement system[END_REF] It consists in a high-speed camera tracking a moving edge in order to simulate the smooth pursuit perform by the eyes. Several configurations exists (moving camera on a rail, rotating mirror, etc.) as well as alternative procedures such as the usage of a stationary high-speed camera combined with motion compensation. [START_REF] Feng | Comparison of motion blur measurement in LCD[END_REF][START_REF] Tourancheau | LCD motion-blur estimation using different measurement methods[END_REF] Moreover, analysis of LCD motion-blur formation blur has shown that the blurred profile of a moving edge can be obtained from the temporal step response of the LCD pixels. [START_REF] Pan | LCD motion blur modeling and analysis[END_REF][START_REF] Watson | The Spatial Standard Observer: A human vision model for display inspection[END_REF] As temporal step response measurements are easy to carry out and can be done rapidly for a lot of gray-to-gray transitions, such a method is very useful. A comparison between this analytic method based on temporal response and some camera measurements has been performed for a large number of gray-to-gray transitions, on five recent monitors. [START_REF] Tourancheau | Motion blur estimation on LCDs[END_REF] The study has shown that both set of results were very close, even for displays with motion-blur reduction systems such as backlight flashing.

These objective physical measurements permits to know the edge profile that reaches the retina (if the effects of the optics of the eye are neglected). However, some process occur in the retina and further in the visual cortex, and one can assume that the perception of the edge will be different from physical measurements. Only few works have been done previously to explore the perception of LCD motion-blur. Someya 8 performed some blur matching experiment between a moving edge on LCD and a moving edge on CRT. Forty-two transitions have been tested but with only one velocity. A certain accordance between subjective and measured data has been demonstrated, but the complexity of the experimental setup, and particularly the physical differences between LCD and CRT, does not permit to draw some straightforward conclusions. Teunissen et al. [START_REF] Teunissen | Method for predicting motion artifacts in matrix displays[END_REF] have also developed a blur matching experiment in which the speed of a moving edge was adjusted by observers until the motionblur induced by the speed matches the blur of stationary blurred edge simulated from measurements. Eight transitions and four velocities have been tested, and one display. A very good match has been reported between speeds adjusted by the observers and speeds used to simulate the stationary blurred edge. Finally, Feng 10 has performed the same type of experiments with one transition and four velocities. Observers were asked to adjust a simulated stationary blurred edge until it matches their perception of a moving edge. A good correlation between subjective results and measurements has been reported, and moving edge was perceived sharper than predicted by measurements.

Those previous results only cover a limited number of transitions, few velocities and one display each. Moreover, some contradictions exists between results of Teunissen et al. and those of Feng, maybe due to the difference of procedures (speed adjustment versus blur adjustment). The goal of this paper is to explore the perception of blur on a moving edge for different gray levels and scrolling speeds with a blur matching procedure, and to repeat experiment on a couple of displays. The relation between psychophysical results and measurements would permit to determine the actual perception of motion-blur by human observers.

PHYSICAL MEASUREMENTS

In these study, the physical measurement method is based on the temporal step-response, because it presents a higher sample-rate than the camera method, and is easier to carry out and less time-consuming. Moreover, it is easily reproducible from one lab to another and this was an important condition for us since our experiments have been led in two different lab.

Since this analytic method has been widely described in some previous works, [START_REF] Tourancheau | LCD motion-blur estimation using different measurement methods[END_REF][START_REF] Pan | LCD motion blur modeling and analysis[END_REF][START_REF] Watson | The Spatial Standard Observer: A human vision model for display inspection[END_REF] details are not provided here. The blurred profile of a moving edge is obtained by a convolution of the temporal step-response with a unit one-frame-period wide window. The blurred edge profile is obtained in the temporal domain and the blurred edge time (BET ) is measured on it between 10% and 90% of the edge dynamic (cf. Figure 1). As a result, a value of BET (in temporal units) is obtained for each gray-to-gray transition.

BET is independent of the speed but the actual blur width (in spatial units) will depend on it. It has been shown from measurements as well as from analysis that, for a given gray-to-gray transition N i → N j , the blurred edge width (BEW ) varies linearly with the velocity of the moving edge V with a factor equal to BET . Here, BEW m denotes the measured blurred edge width:

BEW m = BET • V (1)
In the following, V is expressed in pixels per frame, BEW m in pixels and BET in frames.

PSYCHOPHYSICAL EXPERIMENT

Description

Psychophysical experiment has been conducted according to a procedure that is very similar to that used by Feng: 10 observers have to adjust the blur of a stationary edge until it matches their perception of the motion-blur occurring on a moving edge. The test field is split into two parts as illustrated in Figure 2a. The stationary

Id

Display type Size Native resolution and frame-rate Peak luminance (L max ) 1 monitor display 26" 1920×1200@60Hz 300 cd/m² 2 monitor display 30" 2560×1600@60Hz 160 cd/m² 3 raw TV panel 37" 1920×1080@60Hz 440 cd/m² blurred edge is displayed in the upper part while a sharp edge is moving from left to right on the lower part. The method of adjustment has been used: observers have to track the moving edge in order to perceive the motion-blur and then they adjust the blur of the stationary edge until it matches their perception of the moving one. The moving edge was scrolling across the screen until observers decided to validate their adjustment, with a few-seconds break between each crossing.

The moving edge is defined by a start gray level N s and a final gray level N f (here the edge is moving from left to right, so N s is on the right part of the edge and N f on the left part). As described previously, the spatial profile of a moving edge can be obtained from the temporal step-response of the corresponding gray-to-gray transition N s → N f . The blurred edge profiles obtained from the measurements described in the previous part have been used to set the stationary edge.

Observers have to vary the blur on the stationary edge (following the method of adjustment) until the point of subjective equality between both edges was reached. The adjustment was done by up-or down-scaling the stationary blurred edge obtained from measurements. The stationary edge could vary from a sharp edge to a very blurred edge as shown in Figure 2b.

Stimuli were generated with Matlab on a PC micro-computer using the PsychToolbox extension. [START_REF] Brainard | The Psychophysics Toolbox[END_REF] Twenty gray-to-gray transitions have been tested:

N s → N f | {N s , N f } ∈ [0; 63; 127; 191; 255] 2 , N s = N f

Displays under test

Three liquid-crystal displays have been tested in this study, they are described in Table 1. Brightness and contrast parameters have been set up according to manufacturers recommendations, and a color calibration have been performed. It should be noticed that display 2 has a backlight flashing system in order to reduce motion-blur.

Viewing conditions

Experiments were performed in a psychophysics test room, with a surrounding luminance around 10 cd/m 2 and a D65 chromaticity on the walls. The screen is filled with a mean luminance L mean corresponding to the geometric mean of the maximum and minimum luminance of the display, L mean = √ L max L min . The viewing distance has been set in such a way that one pixel subtend 1.5 visual minute, for each display. This corresponds to a viewing distance of 1.87H for display 1, 1.37H for display 2, and 2.2H for display 3. Whatever the display, the area where stimuli were displayed was 25-degree wide and 7.5-degree high.

Observers

Two labs have been involved in these experiments: IRCCyN-IVC lab (University of Nantes, France) and Acreo AB, Video and Display Quality lab (Sweden). Displays 1 and 3 have been tested in IRCCyN-IVC. Three subjects participated to the experiments: one of the authors and two naive observers. Display 2 has been tested in Acreo. Three subjects participated to the experiments: two authors and one naive observers. All subjects possessed normal or corrected-to-normal vision (visual acuity of 0.9 or better on both eyes). They were familiar with the procedure after several training sessions, the results of which have not been considered.

RESULTS

Figure 3 illustrates the results of the experiment, for the three displays under test, for transitions 0 → 63 and 63 → 191. On display 1 (first row), nine speeds V have been tested: 4, 6, 8, 10, 12, 14, 16, 18, and 20 pixels per frame. At the viewing distance used in the experiment, that corresponds to 6, 9, 12, 15, 18, 21, 24, 27, and 30 degrees per second. The dashed line represents the measured blurred edge width BEW m = BET • V (cf. Equation 1) and dots represents the perceived blurred edge width that we write BEW p , with the standard deviation of the subjective data represented by vertical bars. Regarding the relative stability of results on display 1 (the standard deviation of the subjective data is quite low despite of the low number of observers) we decided to test only four speeds on displays 2 and 3: 10, 15, 20, and 25 pixels per frame (i.e. 15, 22.5, 30, and 37.5 degrees per second).

It has been observed that, for each transition and for each display, the perceived blur varies linearly with the scrolling speed:

BEW p = k • V (2)
This model is represented by the red line in Figure 3. The linear correlation coefficient between the model and the subjective data was always higher than 0.982 on the 20 transitions and 3 displays.

From Equations 1 and 2, it comes that the relationship between the perceived blurred edge width BEW p and the measured blurred edge width BEW m is linear:

BEW p = A • BEW m (3)
A is the ratio between perceived motion-blur and measured motion-blur. For the majority of the tested transitions and display, this ratio is lower than 1. This indicates that the blurred edge is actually perceived sharper than the edge obtained from objective measurements. This result may be related to motion sharpening phenomenon: [START_REF] Bex | Sharpening of drifting, blurred images[END_REF][START_REF] Hammett | Motion sharpening: evidence for the addition of high spatial frequencies to the effective neural image[END_REF][START_REF] Takeuchi | Sharpening image motion based on the spatio-temporal characteristics of human vision[END_REF] blurred edges look sharper when they are moving on the retina than when stationary. [START_REF] Georgeson | Seeing blur: 'motion sharpening' without motion[END_REF] Since, LCD motion-blur occurs because of the retina move on a stationary edge, one can assume that motion sharpening is involved here.

DISCUSSION 5.1 Comparison with other studies

The main result of our experiments is that the ratio A between perceived and measured motion-blur is almost always lower than 1 (from 0.68 to 1.03 over all transitions and displays, with an average value of 0.91). This result is in accordance with that of Feng, 10 but his study did not explore a lot of transitions.

Teunissen et al. have obtained some quite different results. [START_REF] Teunissen | Method for predicting motion artifacts in matrix displays[END_REF] They tested 8 transitions with 4 motion speeds and 17 observers on one display. Their blur matching procedure consisted, as ours, in the equalization of a static blurred edge obtained from measurements with a moving edge on which motion-blur occured. The adjustment was done by varying the speed of the moving edge: this is main difference with our experiment. For each of the 8 transitions and for the 4 tested speeds, their results showed that the adjusted speed was very close to the speed used to simulate the still blurred edge. In other terms, the ratio between perceived and measured motion-blur was equal to 1 for all transitions. Moreover, observers were asked to rate the match between both blurred edges with a five-grade quality scale and the average score were above 4.4 in all cases, indicating a very good match.

The difference of tasks between our two experiments may explain the different results, as well as the fact that, in their experiment Teunissen et al. used a moving block instead of a moving edge in our experiment. As a consequence, for a block with a gray level N i moving on a background with a gray level N j , the blur matching procedure was done on both transitions N i → N j and N j → N i , whereas in our experiment only one transition was tested at a time. The equalization of both rising and falling edges at the same time may have introduced a bias.

Evolution of ratio A as a function of transitions

It has been observed that the ratio A between perceived motion-blur and measured motion-blur is lower for edges with low dynamic and tends towards 1 as the dynamic of the edge is increasing. For example, Figure 4 presents the perceived motion-blur as a function of measured motion-blur for transitions 0 → 63 and 0 → 255 on display 1. The red line is fitted to subjective data according to Equation 2. The ratio A is the slope of the line: 0.69 for transition 0 → 63 and 0.98 for transition 0 → 255.

On Figure 5, the ratio A between perceived and measured motion-blur is plotted as a function of the edge dynamic, for displays 1 and 2. The edge dynamic is defined as the absolute difference between the start luminance L s and the final luminance L f . It can be observed that A is increasing as the edge dynamic increases. The relationship is slightly different according to the type of transitions. For rising transitions (L s < L f ), linear correlation coefficients are quite reliable (0.868 and 0.836), but for falling transitions (L f < L s ), the growing of the relationship is less obvious and the linear correlation coefficients are not so relevant (less than 0.7).

Figure 6 presents the same results as Figure 5.a but here the transitions are gathered regarding the value of the start gray level N s . One can see that for a given start gray level N s , the evolution of the ratio A is increasing rather linearly as N f increases. Of course, more results are needed to draw sound conclusions, but it is another evidence for the fact that ratio A is decreasing with contrast.

The study of the ratio A on display 3 did not permit to bring out a particular tendency. This may be due to the fact that display 3 presented some important overshoots on its temporal step-responses, and so on the blurred edge profile, due to overdriving. As a consequence, the matching between simulated and perceived edges was more difficult to obtain.

CONCLUSION

The psychophysical study of LCD motion-blur described in this work has shown that objective motion-blur obtained from physical measurements and perceived motion-blur are highly correlated. However, perceived blur is lower than physical blur for almost all cases. This result may be related to the motion sharpening phenomenon. Results show that the ratio between perceived blur and measured blur is close to 1 for high contrasted transitions and is decreasing as the dynamic of the edge descreases. This indicates that sharpening is stronger for low contrasted edges. This effect is particularly visible for rising transitions, less for falling transitions.

Current VESA recommendation about motion-blur measurements [START_REF] Vesa | Flat Panel Display Measurements (FPDM2)[END_REF] suggest to take the average of motionblur measurements over transitions as a motion-blur value for the display. Results of this work tend to prove that perception of LCD motion-blur is not constant with transitions and that measurements should be weighted as a function of contrast to obtain a "perceptual" motion-blur score.

  (a) Temporal step-response. (b) Blurred edge profile.

Figure 1 .Figure 2 .

 12 Figure 1. Blurred edge profile and blurred edge time BET (b) obtained from the temporal step-response of the considered transition (a) (see text).

  per frame) Perceived Measured (e) 0 → 63 on display 3.

  per frame) Perceived Measured (f) 63 → 191 on display 3.

Figure 3 .

 3 Figure 3. Measured motion-blur (dashed line) and perceived motion-blur (dots) as a function of speed. Transitions 0 → 63 (left) and 63 → 191 (right), on display 1 (first row), display 2 (second row) and display 3 (third row). Vertical bars represent the standard deviation of the subjective data.

Figure 4 .

 4 Perceived motion-blur as a function of measured motion-blur, for transitions 0 → 63 (left) and 0 → 255 (right), for display 1. Vertical bars represent the standard deviation of the subjective data. The red line represents the model BEWp = A • BEWm, with A = 0.69 for 0 → 63 (left) and A = 0.98 for 0 → 255 (right).

Figure 6 .

 6 Figure 6. Ratio A between perceived and measured motion-blur as a function of the edge dynamic, for rising transitions on display 1. Idem as Figure5.a but transitions are linked regarding the start gray level Ns.

Table 1 .

 1 Displays under test.