Nyassoke Titi

Gaston Clément

Sadefo Kamdem
email: jules.sadefo-kamdem@umontpellier.fr

Louis Aimé Fono

Optimal Renewable Resource Harvesting model using price and biomass stochastic variations: A Utility Based Approach

Keywords: Stochastic Gilpin-Ayala, CRRA utility, Viscosity solutions, Renewable Resources, Optimal Effort. 1

In this article, we provide a general framework for analyzing the optimal harvest of a renewable resource (i.e. fish, shrimp) assuming that the price and biomass evolve stochastically and harvesters have a constant relative risk aversion (CRRA) . In order to take into account the impact of a sudden change in the environment linked to the ecosystem, we assume that the biomass are governed by a stochastic differential equation of the 'Gilpin-Ayala' type, with regime change in the parameters of the drift and variance. Under the above assumptions, we find the optimal effort to be deployed by the collector (fishery for example) in order to maximize the expected utility of its profit function. To do this, we give the proof of the existence and uniqueness of the value function, which is derived from the Hamilton-Jacobi-Bellman equations associated with this problem, by resorting to a definition of the viscosity solution.

Introduction

Throughout the world, renewable resources have catapulted to the forefront of the environmental issues and economic developments. Many researchers have taken a particular interest in fishery fields. They used some stochastic biological production functions for resource growth for purposes of sustainable development. Some stochastic models that account for sudden, unforeseeable events are quite similar in most resource problems such as forestry, fishery, and so on. In general, these problems involve finding an optimal sustainable economic policy in continuous time with stochastic price or biomass growth.

Evolutionary ecology can be viewed as a branch of biology. She is particularly interested in describing and understanding the variability of ecological systems (e.g. fisheries or forestry), through biodemographic processes. It also allows the study or identification of the distribution of abundance of different types of organisms or species in an ecological system. As a corollary, it involves understanding the evolutionary properties of species and their interactions, depending on abiotic or biotic factors of their environment. The environment of a fish species, for example, includes both physical properties dependent on insolation, climate, geology, human activities exploiting ecosystem services, as well as other organisms that share its habitat. In particular, population ecology allows the study of how the size of populations of species living together in groups changes over time or space ( [START_REF] Begon | Ecology: From individuals to ecosystems[END_REF]).

Sustainable development comes down to considering how best to manage, on the basis of very long-term objectives, interactions between natural and social sources of variability. From an economic point of view, the behavior of an entrepreneur will be different depending on whether its investment cycle is shorter or longer than that of the ecosystem it operates. Economic or environmental decisions should be taken under constraint of maintaining the viability of ecosystems, in order to sustain the ecosystem services (forestry, fisheries, etc.) by economic agents.

In the early 1970s, economists were interested in renewable resources and the environment, over-exploitation of open access fisheries, pollution etc., long after pioneering work such as that of [START_REF] Schaefer | Some considerations of population dynamics and economics in relation to the management of the commercial marine fisheries[END_REF], [START_REF] Gordon | The economic theory of a common-property resource: the fishery[END_REF]. In economics, the literature considers resources as stocks to be exploited rationally. The problem is therefore the optimal rate of extraction. Since [START_REF] Hotelling | The economics of exhaustible resources[END_REF] considered more difficult by the existence of a renewal ( [START_REF] Pearce | Environment and Economic Development: The Sustainable Management of Natural Resources in the Deveoloping World[END_REF]). A fortiori, this living resource depends on others with which it interacts within an ecosystem. The evolution over time of their exploitation is a logic of both flow and variability. Authors such as Voltera in the 1930's, [START_REF] Clark | Mathematical optimization and the economics of natural resources[END_REF], among others, extended Hotelling's work in the case of resources renewable. Against all evidence, the assumptions of the model of [START_REF] Schaefer | Some aspects of the dynamics of populations important to the management of the commercial marine fisheries[END_REF], lead to hold the fishermen solely responsible for fluctuations in abundance. Such a model does not take into account the possible impact on the stock, of marine pollution, a tsunami, drastic climate change or an epidemic. Indeed, as indicated above, a living organism modifies the environment in which it lives, and adapts to exogenous modifications (i.e. natural disasters, extreme pollution, climate, etc.) impacting this environment. In the living world, the viability of a complex organism, a tree, a fish, relies on the combination of many and different rhythms (microorganisms, cells and certain functions). Organisms are alive when they have the faculties to reproduce, to modify the environment in which they evolve and to adapt to exogenous modifications of this environment. These organisms interact, exchange matter, energy, information. They confront each other, cooperate, coexist. They fluctuate in complex ways, at multiple spatial and temporal scales, as their interactions also fluctuate.

Based on the work of [START_REF] Schaefer | Some aspects of the dynamics of populations important to the management of the commercial marine fisheries[END_REF], scientific research has proposed several alternative models capable of taking into account exogenous environmental factors. Among these models, the stochastic version of the Gilpin-Ayala differential equation ( [START_REF] Jiang | A note on nonautonomous logistic equation with random perturbation[END_REF]) is increasingly used. However, this equation considers constant drift and variance. Like other authors, to characterize the stochastic dynamics of the stock, we consider in our article, a stochastic differential equation of the 'Gilpin-Ayala' type, with a jump component, and with a drift and a variance liable to change according to the regime. The model considered makes it possible to take into account the impact of possible extreme events (exogenous or not), liable to significantly reduce the stock of resources.

In general, the authors are concerned with the question of how to maximize profit, which is just the difference between revenues and costs. More recently, it is recognized that fishers are risk averse. We introduced in the commonly used model the economic notion of utility function to take into account the fisher's risk attitude. This has made economic model more complex, but richer.

Even though, economics has enhanced our understanding of the problems of resources management, the biological dynamics of the fishery resource that underline the optimization problem make the fisheries more difficult than forestry problems or others. There is a wide variety of population growth models in the literature. In this work, we focus on Gilpin-Ayala growth model for fishery resources, our results may be easily adapted to many other fields.

The model that has been traditionally considered by many authors is the logistic model described by the ordinary differential equation:

dX t = X t [r -aX t ] dt, (1)
where X t stands for the population size at time t, r > 0 represents the growth rate of the species while a > 0 is the environmental carrying capacity. Many authors ( [START_REF] Gilpin | Global models of growth and competition[END_REF], [START_REF] May | Complexity and stability in model ecosystems[END_REF], [START_REF] Gopalsamy | Stability and oscillations in delay differential equations of population dynamics[END_REF]) claimed that a little more complicated model was needed and proposed their following model:

dX t = X t r -aX λ t dt, (2)
where λ > 0 denotes the parameter to modify the classical deterministic logistic model, which is often called GA parameter.

The stochastic nature of renewable resource should be considered in the modeling process. To this end, stochastic versions Gilpin-Ayala model have been studied by many authors ( [START_REF] Jiang | Global stability and stochastic permanence of a non-autonomous logistic equation with random perturbation[END_REF], [START_REF] Jiang | A note on nonautonomous logistic equation with random perturbation[END_REF], [START_REF] Jovanović | Dynamics of non-autonomous stochastic gilpin-ayala competition model with time-varying delays[END_REF], [START_REF] Lian | Stochastic delay gilpin-ayala competition models[END_REF]):

dX t = X t r -aX λ t dt + σX t dW (t), (3)
where σ is volatility and W (t) is the standard Brownian motion.

The rest of the paper is organized as follows. Section 2 briefly describes the stochastic dynamics of the fishing population and prices, specifies the decision rule and the model is formulated. Having derived the stochastic dynamic programming related to our problem, we show in section 3 that the value function is viscosity solution, while section 4 concludes the paper.

The model setup

Throughout this paper we let (Ω, F, {F t } t≥0 , P) be a complete probability space with a filtration {F t } t≥0 satisfying the usual conditions (i.e. it is increasing and right continuous while F 0 contains all P-null sets).

Let W (t) and W Y (t), t ≥ 0, be two standard Brownian motions defined on this probability space which are supposed to be indepedent.

Let α(t) be a right-continuous-time Markov chain, F t -adapted with finite state space S = {1, 2} and generator

Q = (q ij ) ∈ R 2 × R 2 .
We consider in this paper that the Markov chain α(t) is irreducible, which is equivalent to the condition that q ij > 0 for i = j. We assume that the Markov chain α(.) is independent of the Brownian motions W Y (.) and W (.), t is time, t ∈ [0, T ] and T is finite-horizon of time.

Switching jump-diffusion Gilpin-Ayala population model

We set an SDE under regime switching of the form:

dX t = X t r α(t) -a α(t) X λ t -qE α(t) (t) dt + σ α(t) X t dW (t); α(t) = 1, 2 (4)
where r α(t) intrinsic rate of growth in regime α(t), E α(t) is the fishing effort exerted on the population at time which depends on the current regime α(t), q > 0 is the catchability coefficient, λ > 0 denotes the parameter to modify the classical deterministic logistic model, which is often called GA parameter,

a α(t) = r α(t) /K α(t) > 0, K α(t)
is the environmental carrying capacity and σ α(t) is volatility in regime α(t).

With initial value X 0 = x 0 ∈ (0, K), α(t) = i.

As said above, population equations may suffer abrupt environmental shocks. Introducing Levy jump into equations, we have extended Eq. ( 4) to consider the effect of bio-disaster in optimal exploitation of fishing resources. Therefore, the corresponding equations may be expressed as follows

dX t = X t-r i -a i X λ t--qE i (t) dt + σ i X t-dW (t) + R\{0} η(t, X t-, i, z) Ñi (dt, dz); i = 1, 2 (5)
where

X
|η(t, x, i, z)| < ρ(z)(1 + |x|). (7)
Readers are referred to [START_REF] Applebaum | Lévy processes and stochastic calculus[END_REF] for more precise properties of random measures.

When consider η(t, i, z) = 0 for all (t, i, z) ∈ [0, T ] × S × R, i.e there is no jump, many authors (see [START_REF] Li | The stationary distribution and ergodicity of a stochastic generalized logistic system[END_REF], [START_REF] Vasilova | Stochastic gilpin-ayala competition model with infinite delay[END_REF], [START_REF] Jiang | Global stability and stochastic permanence of a non-autonomous logistic equation with random perturbation[END_REF], [START_REF] Jiang | A note on nonautonomous logistic equation with random perturbation[END_REF] showed that Eq. ( 5) has a unique global positive solution for any given positive initial value and represented by

X λ t,i = exp λ t 0 r i -qE i (s) - 1 2 σ 2 i ds + t 0 σ i dW (s) 1/X λ 0,i + λa i t 0 exp λ s 0 r i -qE i (τ ) - 1 2 σ 2 i dτ + s 0 σ i dW (τ ) ds , t ≥ 0.

Mean -reverting price

The unit price of fish, Y t , will depend on the mean (or long-term) price Ȳi . Further, let the unit price evolve as the geometric mean-reverting stochastic process

dY t = θ i ( Ȳ0 -Y t )dt + σ Y Y t dW Y (t) (8)
with initial condition Y (0) = y 0 , where the parameters are positve constants, θ is the reversion speed, p0 is a maximum price, σ Y is the volatility of the price. W Y (t) is standardized Brownian motion as set before.

Model ( 8) is a generalization of the Orstein-Uhlenbeck process, where we use the level-dependent volatility instead of the constant volatilty. The mean price, Ȳ0 , attracts the prices in its direction. In others words, when Y t > Ȳ0 the trend term θ( Ȳ0 -Y t ) < 0 and therefore Y t decreases and when Y t < Ȳ0 a similar argument establishes that Y t grows.

Specifying the Decision Problem

The mostly used cost function in fishery management is the quadratic cost function (see [START_REF] Clark | Mathematical Bioeconomics: The Optimal Management Resources[END_REF], [START_REF] Clark | Mathematical Bioeconomics: The Optimal Management of Renewable Resources[END_REF], [START_REF] Brites | Fisheries management in random environments: Comparison of harvesting policies for the logistic model[END_REF], [START_REF] Titi | Fishery management in a regime switching environment: Utility theory approach[END_REF]).

It can be written as

c(E t ) = (c 1 + c 2 E t )E t
where c 1 , c 2 > 0 are constants.

The profit is the difference between revenue and cost. Here π(t) is the instantaneous profit from the harvest of the stock biomass and is given as:

π(t) = qE t X t Y t -(c 1 + c 2 E t )E t (9)
where, qE t X t denotes the volume of harvest and Y t the actual price of the harvest at the time of decision making.

We consider, as [START_REF] Titi | Fishery management in a regime switching environment: Utility theory approach[END_REF], that the fisher is risk-averse and we define the power utility function as

U (x) = (x 1-γ )/(1 -γ), ( 10
)
where γ is the constant relative risk aversion (CRRA) coefficient, γ = 0 denotes risk neutrality, γ < 0 and γ > 0 implies risk loving and risk aversion, respectively. When γ = 1, U (x) = ln(x).

For a time t in the horizon [0, T ], our problem is to maximize the present value V i (t, x, y) for each i ∈ S defined as:

V i (t, x, y) =E x,y,i T t e -β(s-t) π(X s , Y s , E s , i) 1-γ 1 -γ ds + e -β(T -t) V i (X T ) . (11)
Here E x,y,i is the conditional expectation given X(t) = x, Y (t) = y and α(t) = i under P, where T is the finite time horizon β > 0 is a discount factor.

We say that the control process E(t) is admissible if the following three conditions are satisfied:

1. the SDE (5) for the state process X(t) has a unique strong solution;

2. the SDE [START_REF] Jiang | A note on nonautonomous logistic equation with random perturbation[END_REF] for the state process Y (t) has a unique strong solution;

3. E x,y,i T t e -β(s-t) π(X x t , Y y t , E t , i) 1-γ 1 -γ dt + e -β(T -t) V (X T ) < ∞.
The effort is bounded then A, the set of admissible control, is bounded. The fisher strives to maximize total expected discounted utility of both profit and terminal biomass. Therefore, the value function can be written as:

v i (t, x, y) = sup E∈Ai V i (t, x, y). (12)
3 Dynamic Programming and Viscosity Solutions

In this section, we need to prove existence of the viscosity solution starting to prove some properties of the value function, which will be the viscosity solution.

We shall make the following assumptions: there exists ρ > 0 such that for all s, t ∈

[0, T ], x, x ∈ R + , y, y ∈ R + and E ∈ A |l(t, x, y, E) -l(s, x , y , E)| + |m(x, y) -m(x , y )| ≤ ρ [|t -s| + |x -x | + |y -y |] , (13)
and the global linear growth conditions:

|l(t, x, y, E)| + |m(x, y)| ≤ ρ [1 + |x| + |y|] . (14)
Here

l(t, x, y, E) = π(X x t , Y y t , E t , i) 1-γ 1 -γ , and m(x, y) = V (X x T ).
Let's consider the following two functions f and g defined by:

f (t, X t-, E i (t)) = X t-r i -a i X λ t--qE i (t) and g(t, X t-, E i (t)) = σ i X t-.
Let us define the operator I 1 , I 2 and I 3 of the value function

I 1 (t, x, y, D (x,y) v i , D 2 (x,y) v i ) =θ( Ȳ0 -y) ∂v i ∂y + f (t, x, E i (t)) ∂v i ∂x + 1 2 σ 2 y ∂ 2 v i ∂y 2 + 1 2 g 2 (t, x, E i (t)) ∂ 2 v i ∂x 2 . = f (t, x, E i (t)), θ( Ȳ0 -y) .D (x,y) v i + 1 2 T r       σ 2 y 0 0 σ 2 i x 2    .D 2 (x,y) v i    I 2 (t, x, y, v i ) =q ij (v j (t, x, y) -v i (t, x, y)). I 3 (t, x, y, D (x,y) v i ) = R\{0} v i (t, x + η(t, x, i, z), y) -v i (t, x, y) -η(t, x, i, z) ∂v i (t, x, y) ∂x ν i (dz) = R\{0} v i (t, x + η(t, x, i, z), y) -v i (t, x, y) -(η(t, x, i, z), 0).D (x,y) v i ν i (dz).
The Hamilton-Jacobi-Bellman equations associated with the problem ( 12) is:

∂v i ∂t + sup E∈Ai -βv i + π(x, y, E i (t)) 1-γ 1 -γ + I 1 (t, x, y, D (x,y) v i , D 2 (x,y) v i ) + I 2 (t, x, y, v i ) + I 3 (t, x, y, D (x,y) v i ) = 0, (15)
v i (T, x, y) = κ(x, y) x 1-γ 1 -γ f or i, j ∈ {0, 1} κ > 0. ( 16
)
To the best of our knowledge, there is not in general a smooth solution of Eq. ( 15). Consequently, these should be interpreted in a weaker sense, notably in the framework of viscosity solutions introduced by [START_REF] Crandall | Userâs guide to viscosity solutions of second order partial differential equations[END_REF].

Regularity of value functions

In this section, we study the growth and continuity properties of the value functions. First, we need some estimates on the moments of the price and population processes.

Lemma 3.1. For any k ∈ [0, 2] there exists C = C(k, K, T ) > 0 such that for all h, t ∈ [0, T ], x, y, x , y ∈ R + : E|X t,x h | k ≤ C(1 + |x| k ); E|Y t,y h | k ≤ C(1 + |y| k ). E|X t,x h -x| k ≤ C(1 + |x| k )h k 2 ; E|Y t,y h -y| k ≤ C(1 + |y| k )h k 2 . E|X t,x h -X t,x h | k ≤ C|x -x | 2 ; E|Y t,y h -Y t,y h | k ≤ C|y -y | 2 . E sup 0≤s≤h |X t,x h | k ≤ C(1 + |x| k )h k 2 ; E sup 0≤s≤h |Y t,y h | k ≤ C(1 + |y| k )h k 2 .
Proof of Lemma 3.1 see Appendix A.

Proposition 3.1. For any i ∈ S, the value function denoted by v i (t, x, y) satisfies a linear growth condition and is also Lipschitz in couple (x, y) uniformly in t. There exists a constant C > 0, such that

0 ≤ v i (s, x s , y s ) ≤ C(1 + |x s | + |y s |), ∀(s, x s , y s ) ∈ [0, T ] × R + × R + . |v i (s, x s , y s ) -v i (s, x s , y s )| ≤ C(|x s -x s | + |y s -y s |), ∀s ∈ [0, T ], x s , x s ∈ R + , y s , y s ∈ R + .
Proposition 3.2. Under assumptions ( 13) and ( 14)

the value function v ∈ C 0 ([0, T ] × R + × R + ). More precisely, there exists a constant C > 0 such that for all t, s ∈ [0, T ], x, x s ∈ R + , y, y s ∈ R + , |v i (t, x, y) -v i (s, x s , y s )| ≤ C (1 + |x| + |y|)|s -t| 1 2 + |x -x s | + |y -y s | .
Proofs of proposition 3.1 and proposition 3.2 are rejected in appendices B and C.

The Value Function, Viscosity Solution of Hamilton-Jacobi-Bellman equation

In this section, we will first define what we mean by viscosity solutions. Then we will prove that the value function is a viscosity solution.

We consider the Hamilton-Jacobi-Bellman equations as follows:

∂v i ∂t (t, x, y) + sup E∈Ai -βv i (t, x, y) + π 1-γ 1 -γ + θ(ȳ 0 -y) ∂v i ∂y (t, x, y) + x r i -a i x λ -qE i (t) ∂v i ∂x (t, x, y) + 1 2 σ 2 Y ∂ 2 v i ∂y 2 (t, x, y) + 1 2 σ 2 x 2 ∂ 2 v i ∂x 2 (t, x, y) + q ij [v j (t, x, y) -v i (t, x, y)] + R\{0} v i (t, x + η(t, x, i, z), y) -v i (t, x, y) -η(t, x, i, z) ∂v i (t, x, y) ∂x ν i (dz) = 0. ( 17
)
And the related systems:

       ∂u i ∂t + sup E∈Ai -βu i + l + I 1 (s, x s , y s , D (x,y) u i , D 2 (x,y) u i ) + I 2 (s, x s , y s , u i ) + I 3 (s, x s , y s , D (x,y) u i ) = 0, u i (T, x s , y s ) = κ γ x 1-γ s 1-γ f or i, j ∈ {0, 1} κ > 0 (i, s, x s , y s ) ∈ S × [0, T i ] × R + × R + . ( 18
)
Using a notion of viscosity solution introduced by [START_REF] Crandall | Userâs guide to viscosity solutions of second order partial differential equations[END_REF], we prove below the existence and uniqueness of a solution of [START_REF] Clark | Mathematical Bioeconomics: The Optimal Management Resources[END_REF]. To do this, we define the set of measurable functions on [0, T ] × R + × R + with polynomial growth of degree q ≥ 0 as, and( t, x, ȳ) such that u i -φ attains its local maximum at ( t, x, ȳ),

C q ([0, T ] × R + × R + ) = {φ : [0, T ] × R + × R + , measurable | ∃C > 0, |φ(t, x, y)| ≤ C(1 + |x| q + |y| q )}. Definition 3.1. We say that u i ∈ C 0 ([0, T ] × R + × R + ) is called i. a viscosity subsolution of (18) if for any i ∈ S, u i (T, x, y) ≤ κ γ x 1-γ 1-γ , for all x ∈ R + , p ∈ R + and for all functions φ ∈ C 1,2,2 ([0, T ] × R + × R + ) ∩ C 2 ([0, T ] × R + × R + )
∂φ ∂s + sup E∈Ai -βu i + l + I 1 ( t, x, ȳ, D (x,y) u i , D 2 (x,y) u i ) + I 2 ( t, x, ȳ, u i ) + I 3 ( t, x, ȳ, D (x,y) u i ) ≥ 0, (19)
ii. a viscosity supersolution of ( 18) if for any i ∈ S,

u i (T, x, y) ≥ κ γ x 1-γ 1-γ , for all x ∈ R + , y ∈ R + and if for all functions φ ∈ C 1,2,2 ([0, T ] × R + × R + ) ∩ C 2 ([0, T ] × R + × R + ) and (t, x, y) such that u i -φ attains its local minimum at (t, x, y), ∂φ ∂s + sup E∈Ai -βu i + l + I 1 (t, x, y, D (x,y) u i , D 2 (x,y) u i ) + I 2 (t, x, y, u i ) + I 3 (t, x, y, D (x,y) u i ) ≤ 0, (20)
iii. a viscosity solution of (18) if it is both a viscosity sub-and a supersolution of equation [START_REF] Clark | Mathematical Bioeconomics: The Optimal Management Resources[END_REF].

Theorem 3.1. Under assumption [START_REF] Jovanović | Dynamics of non-autonomous stochastic gilpin-ayala competition model with time-varying delays[END_REF], the value function v is a viscosity solution of Eq. ( 17).

Proof of Theorem 3.1. See Appendix D.

Uniqueness, Comparison Principle

In this section, we will use the notion of parabolic superjet and subjet defined by Crandall, Ishii and Lions [START_REF] Crandall | Userâs guide to viscosity solutions of second order partial differential equations[END_REF]. We state comparison principles, from which we obtain the uniqueness of the solution of the coupled system of partial differential equations.

Definition 3.2. Given v ∈ C o ([0, T ] × R × R × S) and (t, x, y, i) ∈ [0, T ) × R × R × S,
we define the parabolic superjet:

P 2,+ v(t, x, y, i) = (c, q, M ) ∈ R × R 2 × S 2 : v(s, x , y , i) ≤ v(t, x, y, i) + c(s -t) + q.((x -x), (y -y)) + 1 2 ((x -x), (y -y)).M ((x -x), (y -y)) + o(|((x -x), (y -y))| 2 ) as (s, x , y ) → (t, x, y)
and its closure:

P2,+ v(t, x, y, i) = (c, q, M ) = lim n→∞ (c n , q n , M n ) with (c n , q n , M n ) ∈ P 2,+ v(t n , x n , y n , i) and lim n→∞ (t n , x n , y n , v(t n , x n , y n , i)) = (t, x, y, v(t, x, y, i)) .
Similarly, we define the parabolic subjet P2,v(t, x, y, i) = -P2,+ (-v)(t, x, y, i) and its closure P2,v(t, x, y, i)

= -P2,+ (-v)(t, x, y, i).
It is proved in [START_REF] Linos | Optimal control of diffustion processes and hamilton-jacobi-bellman equations part i: the dynamic programming principle and application[END_REF] that

P 2,+(-) v(t, x, y, i) = φ ∂t (t, x, y, i), D (x,y) φ(t, x, y, i), D 2 (x,y) φ(t, x, y, i)
and v -φ has a global maximum (minimum) at (t, x, y, i) .

The previous notions lead to new definition of viscosity solutions.

Definition 3.3. u i ∈ C 0 ([0, T ] × R * + × R * + )
satisfying the polynomial growth condition is a viscosity solution of ( 18) if

(1) for any test-function φ ∈ C 1,2,2 ([0, T ]×R * + ×R * + ), if (t, x, y) is a local maximum point of u i (., ., .)-φ(., ., .)

and if (c, q, L 1 ) ∈ P2,+ u(t, x, y, i) with c = ∂φ(t, x, y)/∂t, q = D (x,y) φ(t, x, y) and L 1 ≤ D 2 (x,y) φ(t, x, y), then c + sup E∈Ai {-βv i + l(x, y, E i ) + I 1 (t, x, y, q, L 1 ) + I 2 (t, x, y, φ) + I 3 (t, x, y, q)} ≤ 0, in this case u is a viscosity subsolution;

(2) for any test-function φ ∈ C 1,2,2 ([0, T ]×R * + ×R * + ), if (t, x, y) is a local minimum point of u i (., ., .)-φ(., ., .)

and if (c, q, L 2 ) ∈ P2,u(t, x, y, i) with c = ∂φ(t, x, y)/∂t, q = D (x,y) φ(t, x, y) and L 2 ≥ D 2 (x,y) φ(t, x, y), then c + sup E∈Ai {-βv i + l(x, y, E i ) + I 1 (t, x, y, q, L 2 ) + I 2 (t, x, y, φ) + I 3 (t, x, y, q)} ≥ 0, in this case u is a viscosity supersolution.

The authors [START_REF] Barles | Second-order elliptic integro-differential equations: viscosity solutions' theory revisited[END_REF] proved that definitions 3.2 and 3.3 are equivalent. The second definition is particular suitable for the discussion of a maximum principle which is the backbone of the uniqueness problem for the viscosity solutions theory.

Before state next lemma, we first introduce the inf and sup-convolution operations we are going to use.

R α [U ](z, r) = sup |Z-z|≤1 U (Z) -r • (Z -z) - |Z -z| 2α R α [V ](z, r) = inf |Z-z|≤1 V (Z) + r • (Z -z) + |Z -z| 2α R α [U ](z, r) is called the modified sup-convolution and R α [V ](z, r) the modified inf-convolution. Notice that R α [V ](z, r) = -R α [-U ](z, r).
Lemma 3.2. (nonlocal Jensen-Ishii's lemma [START_REF] Barles | Second-order elliptic integro-differential equations: viscosity solutions' theory revisited[END_REF])

For any i ∈ S, let u i (., ., .) and v i (., ., .) be, respectively, a usc and lsc function defined on

[0, T ] × R + × R + and φ ∈ C 1,2,2 ([0, T ] × R 2 + × R 2 + ) ∩ C 2 ([0, T ] × R 2 + × R 2 + ) if ( t, ( b1 , p1 ), ( b2 , p2 )) ∈ [0, T ] × R 2 + × R 2 + is a zero global maximum point of u i (t, x, y) -v i (t, x , y ) -φ(t, (x, y), (x , y )) and if c -d := D t φ( t, (x 1 , ŷ1 ), (x 2 , ŷ2 )),
q := D (x,y) φ( t, (x 1 , ŷ1 ), (x 2 , ŷ2 )), r := -D (x ,y ) φ( t, (x 1 , ŷ1 ), (x 2 , ŷ2 )) , then for any K > 0, there exists α(K) > 0 such that, for any 0 < α < α(K), we have: there exist sequences

t k → t, (x k , y k ) → (x 1 , ŷ1 ), (x k , y k ) → (x 2 , ŷ2 ), q k → q, r k → r, matrices M k , N k and a sequence of functions φ k , converging to the function φ α := R α [φ](((x, y), (x , y )), (q, r)) uniformly in R 2 + × R 2 + and in C 2 (B(( t, (x 1 , ŷ1 ), (x 2 , ŷ2 )), K)), such that u i (t k , (x k , y k )) → u i ( t, (x 1 , ŷ1 )), v i (t k , (x k , y k )) → v i ( t, (x 2 , ŷ2 )) (t k , (x k , y k ), (x k , y k ))
is a global maximum of u i (., (., .)) -v i (., (., .)) -φ(., (., .), (., .))

(c k , q k , M k ) ∈ P2,+ u i (t k , (x k , y k )); (-d k , r k , N k ) ∈ P2,-v i (t k , (x k , y k )) - 1 α    I 0 0 I    ≤    M k 0 0 -N k    ≤ D 2 (x,y),(x ,y ) φ(t k , (x k , y k ), (x k , y k )).
Here

c k -d k = ∇ t φ(t k , (x k , y k ), (x k , y k )), q k = ∇ (x,y) φ(t k , (x k , y k ), (x k , y k )), r k = ∇ (x ,y ) φ(t k , (x k , y k ), (x k , y k ))
and φ α ( t, (x 1 , ŷ1 ), (x 2 , ŷ2 )) = φ( t, (x 1 , ŷ1 ), (x 2 , ŷ2 )), ∇φ α ( t, (x 1 , ŷ1 ), (x 2 , ŷ2 )) = ∇φ( t, (x 1 , ŷ1 ), (x 2 , ŷ2 )).

We refer the reader to the mentioned paper for a proof. Now we can state our comparison result.

Theorem 3.2. (comparison principle): If u i (t, x, y) and v i (t, x, y) are continuous in (t, x, y) and are, respectively, viscosity subsolution and supersolution of the HJB system [START_REF] Vasilova | Stochastic gilpin-ayala competition model with infinite delay[END_REF] with at most linear growth, then

u i (t, x, y) ≤ v i (t, x, y) f or all (t, x, y, i) ∈ [0, T ] × R + × R + × S.
Proof of Theorem 3.2. See E.

The following corollary follows from Theorems 3.1 and 3.2.

Corollary 3.1. The value function v is a uniques viscosity solution of Eq. ( 17) that has at most a linear growth.

Conclusion

In this paper we have from a theoretical point of view, analyzed the expected utility profit maximizing.

The main contribution of this paper is related to the risk aversion of the fisher with the assumption of a regime switching environment where the dynamic of stock resource is dependent and therefore affected by season. Given that these assumptions also are valid in other types of renewable resources, the model here may have wider applications. Possible examples include forestry resources. The numerical simulations will be examined in further research.

A Proof of Lemma 3.1

1. Let k ∈ [0, 2] and h = s -t.
According to Hölder inequality Eη k ≤ Eη 2 k/2 for ∀k ≥ 0,

E|X t,x h | k ≤ E|X t,x h | 2 k/2 . ( 21
)
Given that

dX t = f (t, X t-, E i (t), X t )dt + g(t, X t-, E i (t))dw(t) + R\{0} η(t, X t-, i, z) Ñi (dt, dz); i = 1, 2.
According to the elementary inequality

M i=1 a i 2 ≤ M M i=1 |a i | 2 , |X t,x h | 2 ≤ 4 |x| 2 + h 0 f (u + t, X t,x u , i)du 2 + h 0 g(u + t, X t,x u , i)dw u 2 + h 0 R\{0} η(t, X t-, i, z) Ñi (dt, dz) 2 .
Using Itô-isometry and Fubini's theorem

E|X t,x h | 2 ≤ 4 x| 2 + h 0 E f (u + t, X t,x u , i) 2 du + h 0 E g(u + t, X t,x u , i) 2 du + h 0 R\{0} η(t, X t-, i, z) 2 ν i (dz)du ,
Using the growth condition on f , g and η, there exists

C 1 ∈ R such that E|X t,x h | 2 ≤ C 1 1 + |x| 2 + h 0 E X t,x u 2 du .
Applying Gronwall's inequality we obtain

E|X t,x h | 2 ≤ C 1 e C1h 1 + |x| 2 i.e E|X t,x h | 2 ≤ C 1 + |x| 2 . ( 22
)
Using [START_REF] Titi | Fishery management in a regime switching environment: Utility theory approach[END_REF] and elementary inequalities (a

1 +a 2 ) k ≤ 2 k-1 (|a 1 | k +|a 2 | k ) and ( |a 1 + a 2 | ≤ |a 1 |+ |a 2 |)
we deduce

E|X t,x h | k ≤ C 1 + |x| k .
The same reasoning gives us

E|Y t,y h | k ≤ C 1 + |y| k .
2. We have

|X t,x h -x| 2 ≤ 3 h 0 f (u + t, P t,x u , i)du 2 + h 0 g(u + t, X t,x u , i)dw u 2 + h 0 R\{0} η(t, X t-, i, z) Ñi (dt, dz) 2 .
Similar arguments as above we obtain

E|X t,x h -x| 2 ≤ C 1 h 0 1 + E X t,x u 2 du, using (22)
we deduce

E|X t,x h -x| 2 ≤ C 1 + |x| 2 h. Hence E|X t,x h -x| k ≤ C 1 + |x| k h k/2 .
3. Let us define the process X t,x s -X t,x s . Put f (u + t, X t,x u , X t,x u , i) = f (u + t, X t,x u , i) -f (u + t, X t,x u , i), ḡ(u + t, X t,x u , X t,x u , i) = g(u + t, X t,x u , i) -g(u + t, X t,x u , i) and η(t, X t,x u , X t,x u , i, z) = η(t, X t,x u , i, z)η(t, X t,x u , i, z). Then Applying Itô's formula, we have

E|X t,x h -X t,x h | 2 ≤ 4 |x -x | 2 + E h 0 f (u + t, X t,x u , X t,x u , i)du 2 + E h 0 ḡ(u + t, X t,x u , X t,x u , i)dw u 2 + h 0 R\{0} η(t, X t,x u , X t,x u , i, z) Ñi (dt, dz) 2 . E|X t,x h -X t,x h | 2 ≤ C |x -x | 2 + h 0 E X t,x u -X t,x u 2 du . Hence E|X t,x h -X t,x h | 2 ≤ C|x -x | 2 .
Similar arguments as above we deduce

E|X t,x h -X t,x h | 2 ≤ C|x -x | 2 ; E|Y t,y h -Y t,y h | k ≤ C|y -y | 2 .
4. Using Doob's inequality for martingale. We get

E sup 0≤s≤h |X t,x h h| k ≤ C(1 + |x| k )h k 2 ; E sup 0≤s≤h |Y t,y h | k ≤ C(1 + |y| k )h k 2 .

B Proof of Proposition 3.1

We first show that v is Lipschitz in (x, y), uniformly in t and its linear growth condition.

v i (s, x s , y s ) = sup E∈Ai E T s e -β(u-s) l(i, u, X s,xs u , Y s,ys u , E u )du + e -β(T -s) m(X s,xs T , Y s,ys T
) .

1. Using elementary inequality | sup A -sup B| ≤ sup |A -B|, from Lipschitz condition (13) of the article on l, m, and from estimates (Lemma 3.1), with k=1, we have

|v i (s, x s , y s ) -v i (s, x s , y s )| ≤ sup E∈Ai E T s e -β(u-s) l(i, u, X s,xs u , Y s,ys u , E u ) -l(i, u, X s,x s u , Y s,y s u , E u ) du + e -β(T -s) m(X s,xs T , Y s,ys T ) -m(X s,x s T , Y s,y s T ) ≤ sup E∈Ai E T s l(i, u, X s,xs u , Y s,ys u , E u ) -l(i, u, X s,x s u , Y s,y s u , E u ) du + m(X s,xs T , Y s,ys T ) -m(X s,x s T , Y s,y s T ) ≤ sup E∈Ai E T s |X s,xs u -X s,x s u | + |Y s,ys u -Y s,y s u | du + |X s,bs T -X s,x s T | + |Y s,ys T -Y s,y s T | ≤ sup E∈Ai T s E |X s,xs u -X s,x s u | + E|Y s,ys u -Y s,y s u | du + E|X s,xs T -X s,x s T | + E|Y s,ys T -Y s,y s T | |v i (s, x s , y s ) -v i (s, x s , y s )| ≤ C |x s -x s | + |y s -y s | .
2. From linear growth condition ( 14) on l, m, and from estimates (Lemma 3.1), with k=1, we obtain

|v i (s, x s , y s )| ≤ sup E∈Ai E T s l(i, u, X s,xs u , Y s,ys u , E u ) du + m(X s,xs T , Y s,ys T ) |v i (s, x s , y s )| ≤ ρ sup E∈Ai E T s 1 + |X s,xs u | + |Y s,ys u | du + (1 + |X s,xs T | + |Y T |) |v i (s, x s , y s )| ≤ ρ sup E∈Ai T s 1 + E|X s,xs u | + E|Y s,ys u | du + (1 + E|X s,xs T | + E|Y s,ys T |) |v i (s, x s , y s )| ≤ C (1 + |x s | + |y s |) .

C Proof of Proposition 3.2

Let 0 ≤ t < s ≤ T . To prove continuity property in time t, we use the dynamic programming principle.

v i (t, x, y) = sup E∈Ai E T t e -β(u-t) l(i, u, X t,x u , Y t,y u , E u )du + e -β(T -t) m(X t,x T , Y t,y T ) = sup E∈Ai E s t e -β(u-t) l(i, u, X t,x u , Y t,y u , E u )du + e -β(s-t) v(s, X t,x s , Y t,y s , i) = sup E∈Ai E s-t o e -βu l(i, t + u, X t,x t+u , Y t,y t+u , E t+u )du + e -β(s-t) v(s, X t,x s-t , Y t,y s-t , i) . 0 ≤ v i (t, x, y) -v i (s, x s , y s ) = sup E∈Ai E s-t 0 e -β(u) l(i, u, X t,x u , Y t,pt u , E u )du + e -β(s-t) v(s, X t,x s-t , Y t,y s-t , i) -v(s, x s , y s , i) + e -β(s-t) -1 v(s, x s , y s , i) .
Applying linear growth condition (13) of the article on l, noting that 0 ≤ 1 -e -β(s-t) ≤ β(s -t) and v satisfies proposition (3.1) of the article, we deduce that:

|v i (t, x, y) -v i (s, x s , y s )| ≤ sup E∈Ai E s-t 0 l(i, u, X t,x u , Y t,y u , E u ) du + e -β(s-t) v(s, X t,x s-t , Y t,x s-t , i) -v(s, x s , y s , i) + e -β(s-t) -1 v(s, x s , y s , i) ≤ (s-t) 1 2 s-t 0 sup E∈Ai E l(i, u, X t,x u , Y t,y u , E u ) 2 du 1 2

+ sup

E∈Ai E e -β(s-t) v(s, X t,x s-t , Y t,y s-t , i) -v(s, x s , y s , i)

+ sup

E∈Ai E e -β(s-t) -1 v(s, x s , y s , i)

≤ (s -t) 1 2 s-t 0 ρ 2 sup E∈Ai 1 + E|X t,x u | + E|p t,y u | 2 du 1 2 + sup E∈Ai E v(s, X t,x s-t , Y t,y s-t , i) -v(s, x s , y s , i) + β|s -t| sup E∈Ai E v(s, x s , y s , i) ≤ |s -t| 1 2 s-t 0 ρ sup E∈Ai 1 + E|X t,x u | + E|Y t,y u | du + sup E∈Ai E v(s, X t,x s-t , Y t,y s-t , i) -v(s, x s , y s , i) + β|s -t| sup E∈Ai E v(s, x s , y s , i) ≤ C |s -t| 1 2 s-t 0 (1 + E|B t,x u | + E|P t,y u |)du + (|x -x s | + |y -y s |) + β(1 + |x s | + |y s |)|s -t| ≤ C (1 + |x| + |y|)|s -t| 1 2 + |x -x s | + |y -y s | .

D Proof of Theorem 3.1

We establish the viscosity super-and sub-solution properties, respectively in the following two steps.

Step 1. v i (t, x, y), i = 1, 2 is a viscosity super-solution of Eq. ( 17).

We already know that v ∈ C 0 ([0, T ] × R + × R + ). We first note that v i (T, x, y) = κ γ x 1-γ 1-γ so, the boundary condition at time

t = T is clearly satisfied. Let (s, x s , y s ) ∈ [0, T ] × R + × R + , i ∈ S and φ ∈ C 1,2,2 ([0, T ] × R + × R + ) ∩ C 2 ([0, T ] × R + × R + ) such that v i (., ., .
) -φ(., ., .) has a local minimum at (s, x s , y s ). Let N(x s , y s ) to be a neighborhood of (s, x s , y s ) where v i (., ., .) -φ(., ., .) take its minimum, we introduce a new test-function ψ as follows:

ψ(., ., ., j) =        φ(., ., .) + [v i (s, x s , y s ) -φ(s, x s , y s )], if j = i, v i (., ., .), if j = i. (23)
This helps us to suppose without any loss of generality that this minimum is equal to 0.

Let τ α be the first jump time of α(t) = α(t) xs,ys,i , i.e. τ α = min{t ≥ s : α(t) = i}. Then τ α > s, a.s. Let θ s ∈ (s, τ α ) be such that the state (X xs,i t , Y ys,i t ) starts at (x s , y s ) and stays in N(x s , y s )

for s ≤ t ≤ θ s . Applying the generalized Itô's formula to the switching process e -βt ψ(t, X t , Y t , α(t)), taking integral from t = s to t = θ s ∧ h, where h > 0 is a positive constant, and then taking expectation

we have

E xs,ys,i e -βθs∧h ψ(θ s ∧ h, X θs∧h , Y θs∧h , α(θ s ∧ h)) = e -βs ψ(s, X s , Y s , i) + E xs,ys,i θs∧h s e -βt -βψ(t, X t , Y t , α(t)) + ∂ψ(t, X t , Y t , α(t)) ∂t + X t r i -a i X λ t -qE i (t) ∂ψ(t, X t , Y t , α(t)) ∂x + θ( Ȳ0 -Y t ) ∂ψ(t, X t , Y t , α(t)) ∂y + 1 2 σ 2 X 2 t ∂ 2 ψ(t, X t , Y t , α(t)) ∂x 2 + 1 2 σ 2 Y ∂ 2 ψ(t, X t , Y t , α(t)) ∂y 2 + q α(t)j (ψ(t, X t , Y t , j) -ψ(t, X t , Y t , α(t))) + R\{0} ψ(t, X t + η(t, X t , Y t , α(t), z), Y t ) -ψ(t, X t , Y t , α(t)) -η(t, X t , Y t , α(t), z) ∂ψ(t, X t , Y t , α(t)) ∂x ν i (dz) dt , α(t) = j. (24)
From hypothesis, for any t ∈ [s,

θ s ∧ h] v i (t, X xs t , Y ys t ) ≥ φ(t, X xs t , Y ys t ) + v i (s, x s , y s ) -φ(s, x s , y s ) ≥ ψ(t, X xs t , Y ys t , i). (25)
Moreover, we have

q ij (v j (t, X t , Y t ) -v i (t, X t , Y t )) ≤ q ij (ψ(t, X t , Y t , j) -ψ(t, X t , Y t , i)), by definition of ψ.
Recalling that (X ys s , Y ys s ) = (x s , y s ) and using Eq. ( 23) and Eq. ( 25) we have Setting τ = E(θ s ∧ h) combining ( 26) and ( 27) and multiplying both sides by 1/(τ -s) > 0, we obtain

E xs,ys,i e -βθs∧h ψ(θ s ∧ h, X θs∧h , Y θs∧h , α(θ s ∧ h)) ≥ + e -βs v i (s, x s , y s ) + E xs,ys,i θs∧h s e -βt -βv i (t, X t , Y t ) + ∂ψ(t, X t , Y t , α(t)) ∂t + X t r i -a i X λ t -qE i (t) ∂ψ(t, X t , Y t , α(t)) ∂x + θ( Ȳ0 -Y t ) ∂ψ(t, X t , Y t , α(t)) ∂y + 1 2 σ 2 X 2 t ∂ 2 ψ(t, X t , Y t , α(t)) ∂x 2 + 1 2 σ 2 Y ∂ 2 ψ(t, X t , Y t , α(t)) ∂x 2 + q ij (v j (t, X t , Y t ) -v i (t, X t , Y t )) + R\{0} ψ(t, X t + η(t, X t , Y t , α(t), z), Y t ) -ψ(t, X t , Y t , α(t)) -η(t, X t , Y t , α(t), z) ∂ψ(t, X t , Y t , α(t)) ∂x ν i (dz) dt . (26) By Bellman's principle e -βs ψ(s, x s , y s , i) = e -βs v i (s, x s , y s ) = sup E∈Ai E xs, ys,i θs∧h s e -βt l(i, t, X s,xs t , Y s,ys t , E t )dt
+ e -β(θs∧h) v i (θ s ∧ h, X
sup E∈Ai E xs,ys,i 1 τ -s θs∧h s e -βt βv i (t, X t , Y t ) - ∂ψ(t, X t , Y t , α(t)) ∂t -X t r i -a i X λ t -qE i (t) ∂ψ(t, X t , Y t , α(t)) ∂x -θ( Ȳ0 -Y t ) ∂ψ(t, X t , Y t , α(t)) ∂y - 1 2 σ 2 X 2 t ∂ 2 ψ(t, X t , Y t , α(t)) ∂x 2 - 1 2 σ 2 Y ∂ 2 ψ(t, X t , Y t , α(t)) ∂x 2 -q ij [v j (t, X t , Y t ) -v i (t, X t , Y t )] - R\{0} ψ(t, X t + η(t, X t , Y t , α(t), z), Y t ) -ψ(t, X t , Y t , α(t)) -η(t, X t , Y t , α(t), z) ∂ψ(t, X t , Y t , α(t)) ∂x ν i (dz) -l(i, t, X t , Y t , E t ) dt ≥ 0. ( 28
)
Letting τ ↓ s and using the dominated convergence theorem, it turns out that e -βs -∂ψ(s, x s , y s , i) ∂t + inf

E∈Ai βv i (s, x s , y s )-

x s r i -a i x λ s -qE i (s) ∂ψ(s, x s , y s , i) ∂x -θ( Ȳ0 -y s ) ∂ψ(s, x s , y s , i) ∂y - 1 2 σ 2 x 2 s ∂ 2 ψ(s, x s , y s , i) ∂b 2 - 1 2 σ 2 Y ∂ 2 ψ(s, x s , y s , i) ∂y 2 -q ij [v j (s, x s , y s ) -v i (s, x s , y s )] - R\{0}
ψ(s, x s + η(s, x s , y s , i, z), y s ) -ψ(s, x s , y s , i)

-η(s, x s , y s , i, z) ∂ψ(s, x s , y s , i) ∂x ν i (dz) -l(i, s, x s , y s , E s ) ≥ 0. ( 29
)
This shows that the value function v i (t, x, y), i = 1, 2, satisfies the viscosity super-solution property (20).

Step 2. v i (t, x, y), i = 1, 2, is a viscosity sub-solution of [START_REF] Vasilova | Stochastic gilpin-ayala competition model with infinite delay[END_REF].

We argue by contradiction. Assume that there exist an i 0 ∈ S, a point (s,

x s , y s ) ∈ [0, T ] × R * + × R * + and a testing function φ i0 ∈ C 1,2,2 ([0, T ] × R * + × R * + ) ∩ C 2 ([0, T ] × R * + × R * + ) such that v i0 (., ., .) -φ i0 (., ., .)
has a local maximum at (s, x s , y s ) in a bounded neighborhood N(x s , y s ), v i0 (s, x s , y s ) = φ i0 (s, x s , y s ), and min -∂φ i0 (s, x s , y s ) ∂t + inf

E∈Ai 0 βv i0 (s, x s , y s )-

x s r i0 -a i x λ s -qE i0 (s) ∂φ i0 (s, x s , y s ) ∂x -θ( Ȳ0 -y s ) ∂φ i0 (s, x s , y s ) ∂y - 1 2 σ 2 x 2 s ∂ 2 φ i0 (s, x s , y s ) ∂x 2 - 1 2 σ 2 Y ∂ 2 φ i0 (s, x s , y s ) ∂y 2 -q i0j [v j (s, x s , y s ) -v i0 (s, x s , y s )] - R\{0} φ i0 (s, x s + η(s, x s , y s , i 0 , z), y s ) -φ i0 (s, x s , y s ) -η(s, x s , y s , i 0 , z) ∂φ i0 (s, x s , y s ) ∂x ν i0 (dz)-l( i0 , s, x s , y s , E s ) , v i0 (T, x s , y s )-κ γ x 1-γ s 1 -γ > 0, i 0 = j. (30)
By the continuity of all functions involved in (30

) (v i0 , ∂φi 0 ∂x , ∂ 2 φi 0
∂x 2 , q ij , l, ...), there exists a δ > 0 and an open ball B δ (x s , y s ) ⊂ N(x s , y s ) such that,

- ∂φ i0 (t, x, y) ∂t + inf E∈Ai 0 βv i0 (t, x, y)- x r i0 -a i x λ -qE i0 (t) ∂φ i0 (t, x, y) ∂x -θ( Ȳ0 -y) ∂φ i0 (t, x, y) ∂y - 1 2 σ 2 x 2 ∂ 2 φ i0 (t, x, y) ∂x 2 - 1 2 σ 2 Y ∂ 2 φ i0 (t, x, y) ∂y 2 -q i0j [v j (t, x, y) -v i0 (t, x, y)] - R\{0}
φ i0 (t, x + η(t, x, y, i 0 , z), y s ) -φ i0 (t, x, y)

-η(t, x, y, i 0 , z) ∂φ i0 (t, x, y) ∂x ν i0 (dz) -l( i0 , t, x, y, E t ) > δ, i 0 = j, (t, x, y) ∈ B δ (x s , y s ) (31) and

v i0 (T, x, y) -κ γ x 1-γ 1 -γ > δ (t, x, y) ∈ B δ (x s , y s ). Let θ δ = min{t ≥ s : (t, X t , Y t ) ∈ B δ (x s , y s )} be the first exit time of (t, X t , Y t ) (= (t, X s,xs t , Y s,ys t ))
from B δ (x s , y s ). Let θ s = θ δ ∧ τ α where τ α is the first stopping time of α(t) xs,ys,i0 . Then θ s > 0, a.s..

For 0 ≤ t ≤ θ s , we have

βv i0 (t, X t , Y t ) - ∂φ i0 (t, X t , Y t ) ∂t -X t r i0 -a i X λ t -qE i0 (t) ∂φ i0 (t, X t , Y t ) ∂x -θ( Ȳ0 -Y t ) ∂φ i0 (t, X t , Y t ) ∂y - 1 2 σ 2 x 2 s ∂ 2 φ i0 (t, X t , Y t ) ∂x 2 - 1 2 σ 2 Y ∂ 2 φ i0 (t, X t , Y t ) ∂y 2 -q i0j [v j (t, X t , Y t ) -v i0 (t, X t , Y t )] - R\{0} φ i0 (t, X t + η(t, X t , Y t , i 0 , z), y s ) -φ i0 (t, X t , Y t ) -η(t, X t , Y t , i 0 , z) ∂φ i0 (t, X t , Y t ) ∂x ν i0 (dz) -l( i0 , t, X t , Y t , E t ) > δ, i 0 = j, (t, X t , Y t ) ∈ B δ (x s , y s ) (32)
and

v i0 (T, x, y) -κ γ x 1-γ 1 -γ > δ (t, x, y) ∈ B δ (x s , y s ). (33)
As previously, we can replace φ i0 by a new test-function ψ defined as follows:

ψ(., ., ., j) =

       φ i0 (., ., .), if j = i 0 , v i0 (., ., .), if j = i 0 . (34)
+ X t r i -a i X λ t -qE i (t) ∂ψ(t, X t , Y t , α(t)) ∂x + θ( Ȳ0 -Y t ) ∂ψ(t, X t , Y t , α(t)) ∂y + 1 2 σ 2 X 2 t ∂ 2 ψ(t, X t , Y t , α(t)) ∂x 2 + 1 2 σ 2 Y ∂ 2 ψ(t, X t , Y t , α(t)) ∂y 2 + q α(t)j (ψ(t, X t , Y t , j) -ψ(t,
+ (θs∧τ ) s e -βt βv i0 (t, X t , Y t ) - ∂φ i0 (t, X t , Y t ) ∂t -X t r i0 -a i X λ t -qE i0 (t) ∂φ i0 (t, X t , Y t ) ∂x -θ( Ȳ0 -Y t ) ∂φ i0 (t, X t , Y t ) ∂y - 1 2 σ 2 X 2 t ∂ 2 φ i0 (t, X t , Y t ) ∂x 2 - 1 2 σ 2 Y ∂ 2 φ i0 (t, X t , Y t ) ∂y 2 -q i0j [v j (t, X t , Y t ) -v i0 (t, X t , Y t )] - R\{0} ψ(t, X t + η(t, X t , Y t , i 0 , z), Y t ) -ψ(t, X t , Y t , i 0 ) -η(t, X t , Y t , i 0 , z) ∂ψ(t, X t , Y t , i 0 ) ∂x ν i (dz) dt , i 0 = j (36)
i.e e -βs v i0 (s, x s , y s )

≥ E xs,ys,i0 e -βτ v i0 (τ, X τ , Y τ , α(τ ))1 {τ <θs} + e -βθs v i0 (θ, X θs , Y θs , α(θ s ))1 {τ ≥θs} + (θs∧τ ) s e -βt l( i0 , t, X t , Y t , E t ) + δ dt ≥ E xs,ys,i0 e -βτ [κ γ X 1-γ τ 1 -γ + δ])1 {τ <θs} + e -βθs v i0 (θ s , X θs , Y θs , α(θ s ))1 {τ ≥θs} + (θs∧τ ) s e -βt l( i0 , t, X t , Y t , E t ) + δ dt ≥ E xs,ys,i0 + (θs∧τ ) s e -βt l( i0 , t, X t , Y t , E t ) dt + e -βθ v i0 (θ s , X θs , Y θs , α(θ s ))1 {τ ≥θs} + e -βτ [κ γ X 1-γ τ 1 -γ ]1 {τ <θs} + δE xs,ys,i0 (θs∧τ )
+ e -βθ v i0 (θ, X θ , Y θ , α(θ))1 {τ ≥θ} + e -βτ [κ γ X 1-γ τ 1 -γ ]1 {τ <θs} + C 0 δ 1 -E xs,ys,i0 e -βτα (38)
which is a contradiction to the DP principle since E xs,ys,i0 e -βτα < 1. Therefore the value function v i (t, x, y), i = 1, 2, is a viscosity sub-solution of the system (18).

This completes the proof of Theorem 3.1.

E Proof of Theorem 3.2

For , , δ, λ > 0, we define the auxiliary functions φ :

(0, T ] × R 2 + × R 2 + → R and Ξ: [0; T ] × R 2 + × R 2 + × S by φ(t, (x, y), (x , y )) = t + 1 2 |(x, y) -(x , y )| 2 + δe λ(T -t) (|(x, y)| 2 + |(x , y )| 2 )
and Ξ(t, (x, y), (x , y ), i) = v i (t, x, y) -u i (t, x , y ) -φ(t, (x, y), (x , y )).

By using the linear growth of v i and u i , we have for each

i ∈ S lim |(x,y)|+|(x ,y )|→∞ Ξ(t, (x, y), (x , y ), i) = -∞.
Then, since v i and u i are uniformly continuous with respect to (t, x, y) on each compact subset of [0, T ]×R + × R + ×S and that S is a finite set, Ξ attains its global maximum at some finite point belonging to a compact K ⊂ [0, T ]×R 2 + ×R 2 + ×S say, (t δ , (x 1δ , y 1δ ), (x 2δ , y 2δ ), α δ ). Observing that 2Ξ t δ , (x 1δ , y 1δ ), (x 2δ , y 2δ ), α δ ≥ Ξ t δ , (x 1δ , y 1δ ), (x 2δ , y 2δ ), α δ + Ξ t δ , (x 1δ , y 1δ ), (x 2δ , y 2δ ), α δ and using the uniform continuity of v i and u i on K we have

1 |(x 1δ , y 1δ ) -(x 2δ , y 2δ )| 2 ≤ v i (t δ , (x 1δ , y 1δ )) -v i t δ , (x 2δ , y 2δ ) + u i (t δ , (x 1δ , y 1δ )) -u i (t δ , (x 2δ , y 2δ )) ≤ 2C|(x 1δ , y 1δ ) -(x 2δ , y 2δ )|. Thus, |(x 1δ , y 1δ ) -(x 2δ , y 2δ )| ≤ 2C (39)
where C is a positive constant independent of , , δ, λ. From the inequality, 2Ξ (T, (0, 0), (0, 0), α δ ) ≤ 2Ξ (t δ , (x 1δ , y 1δ ), (x 2δ , y 2δ ), α δ )

and the linear growth for v i and u i , we have:

δ (|(x 1δ , y 1δ )) | 2 + | x 2δ , y 2δ | 2 ≤ e -λ(T -t δ ) v i (t δ , x 1δ , y 1δ ) -v i (T, 0, 0) + u i T, 0, 0 -u i (t δ , x 2δ , y 2δ ) ≤ e -λ(T -t δ ) C 2 (1 + |(x 1δ , y 1δ )| + |(x 2δ , y 2δ )|) . (40) It follows that δ |(x 1δ , y 1δ ) | 2 + |(x 2δ , y 2δ | 2 (1 + |(x 1δ , y 1δ )| + |(x 2δ , y 2δ )| ≤ C 2 .
Consequently, there exists C δ > 0 such that,

|(x 1δ , y 1δ )| + |(x 2δ , y 2δ )| ≤ C δ . (41)
This inequality implies that for any fixed δ ∈ (0, 1), the sets {(x 1δ , y 1δ ), > 0} and {(x 2δ , y 2δ ), > 0} are bounded by C δ independent of . It follows from inequalities (40) and (41) that, possibly if necessary along a subsequence, named again t δ , (x 1δ , y 1δ ), (x 2δ , y 2δ ), α δ that there exists (x 1δ0 , y 1δ0 ) ∈ R 2 + , t δ 0 ∈ (0, T ] and α δ 0 ∈ S such that: lim

↓0 (x 1δ , y 1δ ) = (x 1δ0 , y 1δ0 ) = lim ↓0 (x 1δ , y 1δ ), lim ↓0 t δ = t δ0 , lim ↓0 α δ = α δ0 .
If t δ = T then writing that Ξ t, (x, y), (x, y), α δ ≤ Ξ T, (x 1δ , y 1δ ), (x 2δ , y 2δ ), α δ , we have

u i (t, x, y) -v i (t, x, y) - t -2δe λ(T -t) (|(x, y)| 2 ) ≤ u i (T, (x 1δ , y 1δ )) -v i (T, (x 2δ , y 2δ )) - T - 1 2 |(x 1δ , y 1δ ) -(x 2δ , y 2δ )| 2 -δ(|(x 1δ , y 1δ )| 2 + |(x 2δ , y 2δ )| 2 )
≤ u i (T, (x 1δ , y 1δ )) -v i (T, (x 2δ , y 2δ ))

= [u i (T, (x 1δ , y 1δ )) -v i (T, (x 1δ , y 1δ ))]

+ [v i (T, (x 1δ , y 1δ )) -v i (T, (x 2δ , y 2δ ))]

≤ C 1 |(x 1δ , y 1δ ) -(x 2δ , y 2δ )|
where the last inequality follows from the uniform continuity of v i and by assumption that u i (T, (x 1δ , y 1δ )) = κ γ x 1-γ 1δ 1-γ = v i (T, (x 1δ , y 1δ )). Sending , , δ ↓ 0 and using estimate (39), we have: u i (t, x, y) ≤ v i (t, x, y). Assume now that t δ < T . Applying Lemma 3.2 with u i , v i and φ(t, (x, y), (x , y )) at the point (t δ , (x 1δ , y 1δ ), (x 2δ , y 2δ ), α δ ) ∈ -q i0j [u j (t δ , x δ , y δ ) -u i0 (t δ , x δ , y δ )]

-l( i0 , t δ , x δ , y δ , E t δ ) , u i0 (T, x δ , y δ ) -κ γ x 1-γ δ 1 -γ ≤ 0, i 0 = j and min βv i0 (t δ , x δ , y δ ) -d + inf

E∈Ai 0 -x δ r i0 -a i x λ δ -qE i (s) 1 (x δ -x δ ) + 2δe λ(T -t) x δ -θ( Ȳ0 -y δ ) 1 (y δ -y δ ) + 2δe λ(T -t) y δ -1 2 (σx δ ; σ P )(N δ -2δe λ(T -t) I)

   σx δ σ Y   
-q i0j [v j (t δ , x δ , y δ ) -v i0 (t δ , x δ , y δ )]

-l( i0 , t δ , x δ , y δ , E t δ ) , v i0 (T, x δ , y δ ) -κ γ x 1-γ δ 1 -γ ≥ 0, i 0 = j.

Let us define operators A E (x, v, φ, X, Z) and B E (x, v).

A E (t, x, y, w, X, Y Z) = x r i0 -a i x λ -qE i (t) X + θ( Ȳ0 -y)Y + 1 2 wZw B E (t, x, y, v) = q i0j [v j (t, x, y) -v i0 (t, x, y)].

By subtracting these last two inequalities and remarking that min(x; y)-min(z; t) ≤ 0 implies either x-z ≤ 0 or y -t ≤ 0, we divide our consideration into two cases:

≤ sup

E∈Ai 0 l( i0 , t δ , x δ , y δ , E t δ ) -l( i0 , t δ , x δ , y δ , E t δ )

+ sup

E∈Ai 0 A E t δ , x δ , y δ , (σx δ ; σ P ), 1 (x δ -x δ ) + 2δe λ(T -t δ ) x δ , 1 (y δ -y δ ) + 2δe λ(T -t δ ) y δ , M δ + 2δe λ(T -t δ ) I

-A E t δ , x δ , y δ , (σb δ ; σ P ), 1 (x δ -x δ ) + 2δe λ(T -t) x δ , 1 (y δ -y δ ) + 2δe λ(T -t δ ) y δ , N δ -2δe λ(T -t δ ) I

+ sup

E∈Ai 0 B E (t δ , x δ , y δ u) -B E (t δ , x δ , y δ v)] ≡ I 1 + I 2 + I 3 . (42)

Definition 3 . 4 .

 34 For any usc (upper semi-continuous) function U : R m → R and any lsc (lower semicontinuous) function V : R m → R, we set

( 0 ,- t 2 δ- 2 ≤ 2 | 2 δ- 2 (

 022222 T ) × R 2 + × R 2+ × S, for any ζ ∈ (0, 1) there are d ∈ R, M δ , N δ ∈ S 2 such that:d λδe λ(T -t δ ) (|(x δ , y δ )| 2 + |(x δ , y δ )| 2 ), 1 ((x δ , y δ ) -(x δ , y δ )) + 2δe λ(T -t δ ) (x δ , y δ ), M δ + 2δe λ(T -t δ ) I ∈ P2,+ u(t δ , x δ , y δ , i) d, 1 ((x δ , y δ ) -(x δ , y δ )) -2δe λ(T -t δ ) (x δ , y δ ), N δ -2δe λ(T -t δ ) I ∈ P2,v(t δ , x δ , y δ , i) y),(x ,y ) φ(t δ , (x δ , y δ ), (x δ , y δ )) + ζ D 2 (x,y),(x ,y ) φ(t δ , (x δ , y δ ), (x δ , y δ )) + ζ(2 + 4δ e λ(T -t) ) 2δ + 4ζδ 2 e λ(T -t) )e λ(T -t) x δ , y δ ), (x δ , y δ )) x δ , y δ ), (x δ , y δ )) (x δ , y δ ) -(x δ , y δ )| 2 .Furthermore, the definition of the viscosity subsolution u i and supersolution v i implies thatmin βu i0 (t δ , x δ , y δ ) -dt λδe λ(T -t δ ) (|(x δ , y δ )| 2 + |(x δ , y δ )| 2 ) + inf E∈Ai 0 -x δ r i0 -a i x λ δ -qE i (s) 1 (x δ -x δ ) + 2δe λ(T -t) x δ-θ( Ȳ0 -y δ ) 1 (x δ -y δ ) + 2δe λ(T -t) y δ -1 σx δ , σ y )(M δ + 2δe λ(T -t) I)

Case 1 β 2 δ+

 12 u i0 (t δ , x δ , y δ ) -v i0 (t δ , x δ , y δ ) + t λδe λ(T -t δ ) (|(x δ , y δ )| 2 + |(x δ , y δ )| 2 )

  t-is the left limits of X t , Ñi (dt, dz) = N i (dt, dz) -dtν i (dz) denotes the compensated Poisson measure which is independent of W (t). N i (dt, dz) is a stationary Poisson random measure on [0, T ] × R\{0} with intensity measure dtν i (dz). ν i (dz) is the Lévy Measure for the jump size when the Markov chain α(t) is in state i, i.e. a σ-finite Borel measure on R\{0} with the property

	ρ 2 (z)ν i (dz) < +∞, such that for all t, s ∈ [0, T ],
	R\{0}
	x, y ∈ R and i ∈ S,
	|η(t, x, i, z) -η(t, y, i, z)| < ρ(z)|x -y|	(6)

R\{0}

min(1, z 2 )ν i (dz) < +∞.

And η : [0, T ] × R × S × R → R is a continuous function with respect to ν and is measurable satisfying the following assumption: there exist ρ : R → R + with

  For any first exit time τ ∈ [s, T ]. Applying Itô's formula to the switching process e -βt ψ(t, X t , Y t , α(t)), taking integral from t = s to t = (θ s ∧ τ )-and then taking expectation yield E xs,ys,i e -βθ∧τ ψ(θ s ∧ τ, X θs∧τ , Y θs∧τ , α(θ s ∧ τ )) = e -βs v i (s, x s , y s ) + E xs,ys,i

	s	(θs∧τ )-	e -βt -βψ(t, X t , Y t , α(t)) +	∂ψ(t, X t , Y t , α(t)) ∂t

  X t , Y t , α(t))) used E xs,ys,i e -βθs∧τ ψ(θ s ∧τ, X θs∧τ , Y θs∧τ , α(θ s ∧τ )) = E xs,ys,i e -βθs∧τ ψ(θ s ∧τ, X θs∧τ , Y θs∧τ , α(θ s ∧ τ )-) due to continuity. Noting that the integrand in the RHS of (35) is continuous in t. Using (32), (33) and that v i0 (t, X t , Y t ) ≤ φ i0 (t, X t , Y t ) in (35). Also noting that α(t) = i 0 for 0 ≤ t ≤ θ s , it follows e -βs v i0 (s, x s , y s ) ≥ E xs,ys,i0 e -βθs∧τ φ i0 (θ s ∧ τ, X θs∧τ , Y θs∧τ , α(θ s ∧ τ ))

	in which we
	+	ψ(t, X t + η(t, X t , Y t , α(t), z), Y t ) -ψ(t, X t , Y t , α(t))
	R\{0}
		-η(t, X t , Y t , α(t), z)	∂ψ(t, X t , Y t , α(t)) ∂x	ν i (dz) dt , α(t) = j. (35)

  -βt dt + e -βτ 1 {τ <θs} ≥ C 0 1 -E xs,ys,i0 e -βτα . -βt l( i0 , t, X t , Y t , E t ) dt

	constant C 0 such that,
		(θs∧τ )
	E xs,ys,i0
	It follows that
	v i0 (s, x s , y s )
				(θ∧τ )
	≥	sup	E xs,ys,i0 +
		τ ∈[s,T ],E∈A
				e -βt dt + e -βτ 1 {τ <θs} . (37)
				s
	Now considering the estimate of the term E xs,ys,i0	(θs∧τ ) s	e -βt dt+e -βτ 1 {τ <θs} , there exists a positive

s e s e

In view of condition [START_REF] Jovanović | Dynamics of non-autonomous stochastic gilpin-ayala competition model with time-varying delays[END_REF] on l and from estimate (??), we have the classical estimates of I 1 and I 2 :

Using the Lipschitz condition for u and v, we have

Writing that Ξ(t, (x, y), (x, y), i) ≤ Ξ(t δ , (x δ , y δ ), (x δ , y δ ), i) for i ∈ S and using the inequality (42),

this implies

Sending ↓ 0, with the above estimates of (I 1 ) -(I 2 ) -(I 3 ), we obtain:

Choose λ sufficiently large positive (λ ≥ 2C) and send , δ → 0 + to conclude that u i (t, x, y) ≤ v i (t, x, y).

Case 2 The second case occurs if u i0 (T, x δ , y δ ) -v i0 (T, x δ , y δ ) ≤ 0 and finally that u i (t, x, y) ≤ v i (t, x, y).

This completes the proof.


Olivier Christmann
email: olivier.christmann@ensam.eu

Sylvain Fleury

J Ér Ôme Migaud

Vincent Raimbault

Benjamin Poussard

Thibaut Guitter

Geoffrey Gorisse

Simon Richir

Visualizing the invisible: user-centered design of a system for the visualization of flows and concentrations of particles in the air

Keywords: Flow particles, data representation, visualization, virtual reality, user-centered design

This study presents 2 experiments addressing the representation of scientific data, in particular airflows, with a usercentered design approach. Our objective is to provide users feedback to data visualization designers to help them choose an air flow representation that is understandable and attractive for non-experts. The first study focuses on static markers allowing to visualize an airflow, with information characterizing the direction and the intensity. In a second study, carried out in an immersive virtual environment, two information were added, the temporal evolution and the concentration of pollutants in the air. To measure comprehension and attractiveness, participants were asked to answer items on Likert scales (experiment 1) and to answer User Experience Questionnaire (experiment 2). The results revealed that arrows seem to be a very common and understandable form to represent orientation and direction of flow, but that they should be improved to be more attractive by making them brighter and more transparent, as the representation could occlude the scene, especially in virtual reality. To solve this problem, we suggest giving the users the ability to define the specific area where they want to see the air flow, using a cross-sectional view. Vector fields and streamlines could therefore be applied in a virtual reality context.

Introduction

Many physical phenomena are not directly observable with the naked eye, even if some indirect effects can be understood. This is the case for atmospheric pollution (e.g. pollutants, pollen...), underwater pollution composed of microparticles, or even wind. However, the improvement of numerical computation software and the development of increasingly accurate models make possible to obtain precise simulations in 2 or 3 dimensions, with temporal evolution. Data visualization enables a graphical transcription of these numerical data to allow user a concept exploration which could not be fully apprehended without seeing them [START_REF] Donalek | Immersive and collaborative data visualization using virtual reality platforms[END_REF] . In addition to the "invisible" physical phenomena mentioned above, simulation visualization can be extended to many areas, such as the flow of people during the evacuation of buildings [START_REF] Guest | Interactive analysis and visualization of situationally aware building evacuations[END_REF] .

These simulations results needs to be communicated using representations which are suited to the level of users' understanding, who are still mostly experts in the specific domains [START_REF] Agrawal | Challenges and opportunities with big data visualization[END_REF] . Thus, the chosen representation metaphor will directly influence the users' understanding. To engage a user, expert or not, on a complex topic, a consistent and attractive visual representation of the data should be created [START_REF] Grainger | Environmental data visualisation for non-scientific contexts: Literature review and design framework[END_REF] . Few research works have studied the user experience when considering these representations [START_REF] Quispel | Graph and chart aesthetics for experts and laymen in design: The role of familiarity and perceived ease of use[END_REF] .

The medium can also play an important role in the understanding of visualizations. The advent of emerging technologies such as virtual reality and augmented reality bring new modalities that go far beyond the most mediated uses such as entertainment, education, healthcare, or industry. Recent works show that it is possible to combine computational fluid dynamics methods with real-time 3D engines in order to visualize the simulation results not on a screen or a display wall, but via a stereoscopic vision device, and thus to be immersed at scale 1 in the environment used for the simulation. This paradigm shift opens the possibility of communicating results to non-expert users for sharing, discussion, or sensitization, but introduces a new issue concerning the nature of the representation and the interaction with it. The intelligibility of scientific data representations in virtual reality has been little addressed in the literature.

We thus carried out 2 experimental studies aiming at assessing different representations of scientific data, in the particular case of an air flow. The first study focuses on the static marker allowing to visualize an air flow, with the information characterizing the direction and the force. It is the preliminary to our second study, performed in an immersive virtual environment, which allows two additional information to the user: the temporal evolution on the one hand, and the concentration of pollutants in the air on the other hand. In this work, we are not interested in the possibility of directly modifying simulation parameters in real time in 3D space (for example by dynamically modifying the position of an air purifier in a room), because we aim here at studying the intelligibility of the proposed representations.

Related work

Visualizing scientific data, which can be multidimensional and multivariate, allows for understanding, analyzing, assessing and sharing. These data come from various areas such as meteorology [START_REF] Helbig | Concept and workflow for 3d visualization of atmospheric data in a virtual reality environment for analytical approaches[END_REF][START_REF] Ware | Designing a better weather display[END_REF] , aerodynamics [START_REF] Harwood | A real-time modelling and simulation platform for virtual engineering design and analysis[END_REF] , public health related to pollen [START_REF] Desjardins | A space-time parallel framework for fine-scale visualization of pollen levels across the eastern united states[END_REF] , geosciences [START_REF] Hohl | Computationally enabled 4d visualizations facilitate the detection of rock fracture patterns from acoustic emissions[END_REF] or even archeology [START_REF] Knabb | Scientific visualization, 3d immersive virtual reality environments, and archaeology in jordan and the near east[END_REF] . The constant increase in computational power (CPU and GPU) gives today the possibility to modify some parameters "on-the-fly", like in flow simulations, based on CFD (computational fluid dynamics). CFD is a combination of physics and applied mathematics [START_REF] Kim | Augmented-reality visualization of aerodynamics simulation in sustainable cloud computing[END_REF] , that allow the calculation of gas or liquid flows, and much commercial software (e.g., ANSYS Fluent, Paraview) exist today to address issues such as ventilation inside a vehicle or the behavior of air flows in a room [START_REF] Yudin | Computational hydrodynamics in air flows modeling : Using the unreal engine based on the numerical solution of the navier-stokes equations[END_REF] . A CFD process has three steps 12 :

• Pre-processing: modeling of the geometries, initial and boundary conditions, solver; • Computation; • Post-processing: visualization of the results.

In this paper, we will focus on the last step regarding how to represent the results of a CFD computation, without mentioning the calculation process or the optimization of the whole chain as it is done in numerous research, like in [START_REF] Harwood | Interactive flow simulation using tegra-powered mobile devices[END_REF] .

Air flow simulation

Many works have proposed extensive literature reviews in the field of flow visualization [START_REF] Post | The state of the art in flow visualisation: Feature extraction and tracking[END_REF][START_REF] Mcloughlin | Over two decades of integration-based, geometric flow visualization[END_REF][START_REF] Salzbrunn | The state of the art in flow visualization: Partition-based techniques[END_REF] , resulting in a taxonomy of different types of representations [START_REF] Helbig | Concept and workflow for 3d visualization of atmospheric data in a virtual reality environment for analytical approaches[END_REF][START_REF] Giraldo | Perception of multisensory wind representation in virtual reality[END_REF] . These representations are most often "4D", as they allow the temporal evolution of the simulated data to be examined. Adding a temporal dimension to a 2D or 3D visualization changes it into an animation, which helps to understand temporal characteristics of a phenomenon as well as the correlation between different factors represented simultaneously [START_REF] Yeh | Animating streamlines with orthogonal advancing waves[END_REF] . Although professional software allows the visualization of results from numerical computations, they are hardly accessible to non-experts [START_REF] Stone | Exploiting gaming technologies to visualise dynamic thermal qualities of a domestic dwelling: Pilot study of an interactive virtual apartment[END_REF] . Emerging technologies such as virtual reality and augmented reality can reduce the difficulties of interpreting results, through better contextualization due to immersion in the virtual or real environment, as was suggested already in 1995 [START_REF] Zannetti | Is virtual reality the future of air pollution modelling?[END_REF] . If computational power was a limiting factor to link CFD with extended realities, many works demonstrate the feasibility of a complete chain outside computational software, by directly implementing the Navier-Stoke equations and using video game engines [START_REF] Yudin | Computational hydrodynamics in air flows modeling : Using the unreal engine based on the numerical solution of the navier-stokes equations[END_REF][START_REF] Stone | Exploiting gaming technologies to visualise dynamic thermal qualities of a domestic dwelling: Pilot study of an interactive virtual apartment[END_REF][START_REF] Berger | Cfd post-processing in unity3d[END_REF] , at the cost of lower accuracy, but which may be sufficient for exchanges with non-experts as it allows for rapid modification of scenarios and visualizations [START_REF] Harwood | A real-time modelling and simulation platform for virtual engineering design and analysis[END_REF] . Augmented reality and virtual reality bring complementary benefits: thanks to the former, it is possible to visualize insitu data, which is relevant in a room [START_REF] Kim | Augmented-reality visualization of aerodynamics simulation in sustainable cloud computing[END_REF] , in a house [START_REF] Fukuda | Integration of cfd, vr, ar and bim for design feedback in a design process -an experimental study[END_REF] or in architectural projects [START_REF] Ritterbusch | Augmented reality visualization of numerical simulations in urban environments[END_REF] , but it requires the user to move. Virtual reality, on the other hand, allows the integration of many scenarios, from the real world (via scanning or modelling) or not [START_REF] Giraldo | Perception of multisensory wind representation in virtual reality[END_REF] .

2D and 3D representations

Helbig et al. [START_REF] Helbig | Concept and workflow for 3d visualization of atmospheric data in a virtual reality environment for analytical approaches[END_REF] list the possible representations depending on whether the data is 2D, 3D scalar and 3D vector, emphasizing the interest of finding a representation adapted to each variable, with elements (mainly color) that allow to differentiate them. According to Post et al. [START_REF] Post | The state of the art in flow visualisation: Feature extraction and tracking[END_REF] , 2D representations refer to texture-based visualization, while 3D scalar and vector 3D can be represented by geometric visualization. We will list the most commonly encountered in the literature in the context of flows. When the data is 2D, it is possible to use textures, isolines (the value can then be color and/or scale coded for the last mentioned). Representations of 3D scalar data are slices, iso-surfaces, and iso-volumes, with a possibility of figuring the value as color or opacity. 3D vector data can be represented as streamlines or arrows with color or size matching. These representations are sometimes referred to as glyphs in some works, as they can combine different information by their size, shape, texture orientation and color [START_REF] Ware | Designing a better weather display[END_REF] . Finally, there is a last kind of representation, which can be described as indirect, which consists of representing the effects of a physical phenomenon on contextual elements, such as the movements of a tree and its foliage as a result of wind speed [START_REF] Giraldo | Perception of multisensory wind representation in virtual reality[END_REF] . In a similar approach, some works are centered on multisensory aspects [START_REF] Giraldo | Perception of multisensory wind representation in virtual reality[END_REF][START_REF] Ito | Alteredwind: Manipulating perceived direction of the wind by cross-modal presentation of visual, audio and wind stimuli[END_REF] .

Concerning the representation of flows in the literature, of air in particular, streamlines are the most common, giving the direction and allowing to associate an additional color information. For instance, we find this representation to illustrate the operation of an air purifier [START_REF] Kim | Augmented-reality visualization of aerodynamics simulation in sustainable cloud computing[END_REF] , to visualize the results of an aerodynamic simulation in a wind tunnel [START_REF] Harwood | A real-time modelling and simulation platform for virtual engineering design and analysis[END_REF] , to understand the behavior of the wind in a set of buildings [START_REF] Ritterbusch | Augmented reality visualization of numerical simulations in urban environments[END_REF] or to study the thermal comfort in a dwelling [START_REF] Berger | Cfd post-processing in unity3d[END_REF] . Vectors represented as 3D objects are frequently found in the literature to visualize the direction of an airflow, with additional velocity or force information given by the norm of the vector [START_REF] Helbig | Concept and workflow for 3d visualization of atmospheric data in a virtual reality environment for analytical approaches[END_REF][START_REF] Yudin | Computational hydrodynamics in air flows modeling : Using the unreal engine based on the numerical solution of the navier-stokes equations[END_REF] ; a color could add an informative element. There is most often a semantic link between the colors and the measurement of a variable, this link can also be related to the physical light spectrum [START_REF] Ware | Designing a better weather display[END_REF] .

To represent air quality, visualizations are frequently based on particles, as the concentration can be visualized by color or opacity. Colors typically range from blue to red to reflect quality degradation [START_REF] Kim | Augmented-reality visualization of aerodynamics simulation in sustainable cloud computing[END_REF] , and similarly higher opacity indicates high concentration [START_REF] Yudin | Computational hydrodynamics in air flows modeling : Using the unreal engine based on the numerical solution of the navier-stokes equations[END_REF] . In this perspective, we can mention the volumetric rendering approach based on voxels, used for example to visualize interpolated pollen values [START_REF] Desjardins | A space-time parallel framework for fine-scale visualization of pollen levels across the eastern united states[END_REF] , rock fracture patterns based on acoustic emissions [START_REF] Hohl | Computationally enabled 4d visualizations facilitate the detection of rock fracture patterns from acoustic emissions[END_REF] or geovisualization of COVID-19 data [START_REF] Lan | Geovisualization of covid-19: State of the art and opportunities[END_REF] and which can be applied to air quality visualization. These works use both color and transparency: low density voxels are high transparency to allow visualization of high density, opaque voxels; voxels of close value can be drawn as isovolumes.

The representation of a gas can also be illustrated by a set of "visible" particles, such as small spheres with the same semantic possibilities [START_REF] Kim | Augmented-reality visualization of aerodynamics simulation in sustainable cloud computing[END_REF][START_REF] Yudin | Computational hydrodynamics in air flows modeling : Using the unreal engine based on the numerical solution of the navier-stokes equations[END_REF][START_REF] Stone | Exploiting gaming technologies to visualise dynamic thermal qualities of a domestic dwelling: Pilot study of an interactive virtual apartment[END_REF] . Densities of fine particles are frequently represented as 2D textures (heatmap) where color is an indication of concentration, e.g. in cross-sectional planes [START_REF] Yudin | Computational hydrodynamics in air flows modeling : Using the unreal engine based on the numerical solution of the navier-stokes equations[END_REF][START_REF] Berger | Cfd post-processing in unity3d[END_REF][START_REF] Fukuda | Integration of cfd, vr, ar and bim for design feedback in a design process -an experimental study[END_REF] .

User-centered design

Most of the articles mentioned above are interested in the complete processing of computation and visualization, whether in virtual reality or augmented reality, and aim at technically assessing the use of emerging technologies for the visualization of data from numerical computations. Many works thus focus on optimizing computations to allow implementation on devices with low computational power such as smartphones, as they are interesting for low-cost AR [START_REF] Harwood | Interactive flow simulation using tegra-powered mobile devices[END_REF] . Others are technical demonstrators of an integration of physical equations up to visualization (e.g., Yudin et al. (2019) [START_REF] Yudin | Computational hydrodynamics in air flows modeling : Using the unreal engine based on the numerical solution of the navier-stokes equations[END_REF] ). The question of understanding the proposed representations, yet crucial when they are addressed to nonexperts, is slightly addressed. The development of new real-time 3D simulation tools requires a multidisciplinary approach, with knowledge in physics, mathematics and software development [START_REF] Harwood | A real-time modelling and simulation platform for virtual engineering design and analysis[END_REF] , but we need to take into account human factors and to consider how to make these data intelligible to experts and non-experts. Like Lan et al. (2021) [START_REF] Lan | Geovisualization of covid-19: State of the art and opportunities[END_REF] , we believe that visualizing complex data can be useful not only for decision makers, but also for the general public. The user experience of these augmented reality and virtual reality tools is rarely assessed, even though it determines their interest for the end user. Among these few works, we can mention the study of Giraldo et al. (2020) [START_REF] Giraldo | Perception of multisensory wind representation in virtual reality[END_REF] , which evaluates the effects of 3 different wind representations (including a tactile restitution and an external restitution) on the perception and the feeling of presence, or the study of Stone et al. (2014) [START_REF] Stone | Exploiting gaming technologies to visualise dynamic thermal qualities of a domestic dwelling: Pilot study of an interactive virtual apartment[END_REF] , which applies humancentered design principles in the design of the virtual reality tool dedicated to evaluate the thermal quality of a house. In the case of 2D representations, Quispel et al. (2015) [START_REF] Quispel | Graph and chart aesthetics for experts and laymen in design: The role of familiarity and perceived ease of use[END_REF] focused on ergonomic criteria to study the relationship between familiarity, perceived ease of use and attractiveness of graph designs. The need for a user-centered design approach is also evoked by Berger and Cristie (2015) [START_REF] Berger | Cfd post-processing in unity3d[END_REF] and Lan et al. (2021) [START_REF] Lan | Geovisualization of covid-19: State of the art and opportunities[END_REF] . Beyond intelligibility, a representation can also be unsuitable for virtual reality or augmented reality, such as particles with an opacity gradient that hinders the user from seeing enough of the environment, or 3D visual elements (streamlines, vectors) too close to the eyes of the user and therefore disturbing. It is therefore necessary to conduct ergonomic evaluations of these representations, and to take the user into account when developing tools for visualizing scientific results.

The aim of the user-centered approach is to design software that is adapted to its end-users [START_REF] Wilson | Costs and benefits of user involvement in design: Practitioners' views[END_REF] . The ISO 9241-210 (2019) [START_REF]Ergonomics of human-system interaction -Part 210: Human-centred design for interactive systems[END_REF] standard identifies six principles that characterize user-centered design: 1) the design is based upon an explicit understanding of users, tasks and environments; 2) users are involved throughout design and development; 3) the design process is driven and refined by user-centered evaluation; 4) the process is iterative: 5) the design address the whole user experience; 6) the design team includes multidisciplinary skills and perspectives. A human-centered approach requests four main tasks to be carried out. The first task consists in understanding and specifying the context of use. The second is the specification of user needs and the other stakeholder's requirements. The third is the production of design solutions that can be for instance mock-ups or prototypes. The fourth task is the assessment of the solutions at each stage of the project, from the early concept design to long-term usage, to make the best design choices. There are two main usercentered methods that can be used to carry out this last task, namely user testing and expert evaluation.

In the present work, we opted for two iterations of user testing on mock-ups and prototypes with an experimental comparison approach.

How to represent particle flows in an understandable and pleasant picture?

Hypotheses

In the present experiment, the two parameters that will be represented are the direction and the intensity of an air flow. Overall, it is possible to represent flows in a way that evokes a vector, i.e. with a line that goes directly in the direction of the flow. Another possibility is to use a metaphoric representation, such as a windsock with the position depending on the air flow, but this solution requires an effort of interpretation. Concerning the direction, our first hypothesis was that a vector-like representation (a line in the wind direction or an arrow) would be more understandable than a Windsock metaphorical representation, as it does not require interpretation.

The representation of flow intensity in a vector representation can be based on color, or on the shape of the representation (here, the length of the lines). The use of gradient color codes from green to red appears so common that it can be expected to be more easily interpreted and understood. Regarding intensity, our second hypothesis was that the use of a color code allows a better understanding than a difference in shape (i.e., line size).

Method

Forty-one participants (37 men and 4 women), all students in a French engineering school, carried out the experiment. The youngest were 22 years old and the oldest was 38 years old. Their mean age was 22 years old (SD=2.828).

Each participant was asked to observe on a computer screen a picture showing a 3D environment in which an air movement (i.e. a direction and a force) was symbolized. Three pictures were examined successively by the participants, each of them using different way of representing the air movement (see figure 1) color arrows pointing the direction of the wind, and color representing its intensity, 2) windsocks, like a real one, behave as if they were pushed by the wind, inclining in the direction of wind, more or less depending on its intensity, 3) tab are blue points which unfold under the effect of the wind, thus becoming blue lines in the direction of the wind, with a length depending on its intensity. The order of presentation of the pictures was counterbalanced between the participants. After each picture examination, the participants had to answer the four following item on Likert scales from 1 to 10:

• This representation allows me to visualize the direction of the air movements;

• This representation allows me to visualize the intensity of the air movements; • This representation has good aesthetic qualities; • This visualization prevents me from seeing the scenery.

Results

The Means and standard deviations for the scores concerning the perception of the three pictures are represented on figure 2.

Regarding the distribution of the scores for each item, a Levene test was used to verify homoscedasticity between the three conditions. It revealed acceptable homogeneity of variances for aesthetic qualities (F(2)=.662, p=.518) and the visualization prevention (F(2)=1.321, p=.271), but not for the perceived direction (F(2)=16.186, p<.001) and the perceived intensity (F(2)=3.554, p=.032).

ANOVA with repeated measures revealed a significant difference between the conditions for the evaluation of aesthetic qualities (F(2)=7.914, MS=34.96, η2=.227, p<.001). Pairwise comparison showed that the condition arrows differentiate with tab (p=.002) and windsock (p=.006), but there was no significant difference between tab and windsock for aesthetic qualities (p=.644). Concerning the visualization prevention, ANOVA with repeated measures revealed a significant difference between the conditions (F(2)=6.937, MS=18.179, η2=.221, p=.002). Pairwise comparison showed that the condition windsock differentiates with tab (p=.01), but there was no significant difference between arrows and windsock (p=.10) or tab (p=.34) for visualization prevention.

Friedman test revealed a significant difference between the conditions for the evaluation of the wind direction (chi2=35.796, df=2, W=.437, p<.001). Pairwise comparison showed that the condition arrows differentiate with tab (p<.001) and windsock (p<.001), and tab with windsock (p=.025). Concerning the evaluation of wind force, Friedman test revealed a significant difference between the conditions (chi2=28, df=2, W=.341, p<.001). Pairwise comparison showed that the condition windsock differentiates with tab (p<.001) and arrows (p<.001), and tab with arrows (p<.001).

Discussion

The objective of this first experiment was to evaluate the different representations' intelligibility on the two air flows parameters: direction and force. The first hypothesis was that a vector-like representation would be more understandable in terms of direction than a Windsock metaphorical representation, as it does not require any interpretation. The results are only partially consistent with this hypothesis, as participants rated the colored arrows as better at representing the direction of flow than the two other conditions, but Windsock and Tab are indistinguishable. It is possible that the low score for the Tab condition is due to the fact that the line does give an orientation, but in the absence of any arrow, it is not known in which direction the flow is going.

Regarding intensity, our second hypothesis was that the use of a color code allows a better understanding than a difference in shape (line size). The results support this assumption since, once again, the colored arrows condition is superior to the other two. The coloured arrows condition is thus judged to be the most comprehensible on the two parameters, while also being judged to be more attractive. Finally, we note that the Tab condition is the one that prevents the least from seeing the scenery. This last point will have to be the object of particular attention for the design of flow visualizations.

A limitation of the proposed representations is the necessary discretization of the space-time, which may lead to not accounting for some subtle variations of the flow. In a use case of decision support on the positioning of air purifiers in large volumes (subway station), this is not a blocking point. However, it would be relevant to be able to represent where these subtle changes occur and to provide a zoom in/out mechanism to inspect the corresponding volume.

How to represent understandable and pleasant dynamic particles in a 3D

environment?

Hypotheses

In this second experiment, a third parameter was added in addition to speed and direction: the particle concentration in the air. Overall, the air shown in the scene is cleaner when it comes out of a filter than before passing through it. The representations are no longer static images, but virtual reality animations. Thus, the direction and speed of the air flows are directly represented by the movements.

In virtual reality, it is tempting to try to represent the phenomena as realistically as possible. However, it is questionable whether imitation of reality is really the best way to achieve an intelligible representation. The first assumption was that a schematic representation in the form of 3-dimensional arrows is more comprehensible than a realistic representation of moving particles. In some cases, we are looking for beauty, while in others we are looking for intelligibility. Both criteria may therefore be important depending on the context. The second hypothesis was that a schematic representation appears unattractive, compared to a representation of bright particles with a trail like a comet.

Method

Eight participants (6 men and 2 women), students in computer science, have been asked to carry out the experiment. The youngest was 22 years old and the oldest was 27 years old, for a mean age of 24 years old (SD=1.941). All of them were familiar with using virtual reality. All the participants were asked to explore in virtual reality a modeled room in which dynamic air flows were represented by three specific types of animations (see figure 3). The objective of the three representations were to make understandable the impact of an air filter in the room, by visualizing the dirty air and the clean air flows. In the artistic representation, the flows were represented by blue light trails moving in the air. The schematic representation was flows of 3D arrows, dark for the dirty air and blue for the clean air that came out of the filter. Finally, in the realistic representation, the participants saw gray particles for the dirty air, and blue particles for the clean air.

All the participants were immersed in virtual reality with HTC Vive headset during 5 minutes with each of the three conditions. The order of presentation was counterbalanced. After each presentation, they were asked to answer the French version of the User Experience Questionnaire [START_REF] Laugwitz | Construction and evaluation of a user experience questionnaire[END_REF] .

Results

Regarding the distribution of the scores for each item, a Levene test was used to verify homoscedasticity between the three conditions. It revealed acceptable homogeneity The Means and standard deviations for the scores concerning the perception of the three pictures are represented on figure 4.

ANOVA with repeated measures revealed no significant difference between the three conditions for efficiency (F(2)=2.923, MS=1.255, η2=.227, p=.087), dependability (F(2)=3.51, MS=1.440, η2=.227, p=.006) and stimulation (F(2)=1.462, MS=.536, η2=.114, p=.265).

Concerning novelty, ANOVA with repeated measures revealed a significant difference between the three conditions (F(2)=5.672, MS=4.625, η2=.335, p=.002). Pairwise comparison showed that the condition artistic differentiates with schematic (p=.014) and realistic (p=.008), but there was no significant difference between schematic and realistic for novelty (p=.792).

ANOVA with repeated measures revealed a significant difference between the three conditions for attractiveness (F(2)=8.477, MS=2.272, η2=.279, p=.004). Pairwise comparison showed that the condition artistic differentiates with realistic (p=.010), but there was no significant difference between schematic and artistic (p=.121) or realistic (p=.234) for attractiveness.

Finally, ANOVA with repeated measures revealed a significant difference between the three conditions for perspicuity (F(2)=4.68, MS=3.456, η2=.341, p=.028). Pairwise comparison showed that the condition schematic differentiates with realistic (p=.003), but there was no significant difference between artistic and realistic (p=.151) or schematic (p=.086) for perspicuity.

Discussion

The first hypothesis was that a schematic representation in the form of 3-dimensional arrows is more comprehensible than a realistic representation of moving particles. The results are in line with this hypothesis, as the participants found the application more insightful in the schematic condition than in the realistic condition.

The second hypothesis was that a schematic representation appears unattractive, compared to a representation of bright particles with a trail like a comet. The results are only partially consistent with this hypothesis, because there is no difference between the three conditions concerning the stimulation, a superiority of artistic condition compared to realistic condition for attractiveness, and a superiority of artistic condition compared to schematic condition for novelty.

In virtual reality, we noted that the problem of not being able to see the environment can be detrimental when a shape passes near the user's eyes. This problem seems to be more important in realistic condition, which may explain why evaluations are so low in this condition.

General discussion and conclusion

The objective of the two experiments described in the present paper was to provide users feedback to support data visualization designers' choice for air flow representation in order to make it understandable and attractive. The results revealed that arrows seem to be a very common shape understandable to represent orientation and direction of flow. This kind of schematic shape could be improved to be more attractive by making it brighter and more transparent.

Moreover, a difficulty concerning the shapes preventing the participants from seeing the scene has been noted during the two experiments. To solve this problem, we suggest giving the user the possibility to delimit the area where he wants to see the air flow. We also develop a different concept of representation based on the idea of a specific need to decide where to position the filter in the room to maximize the efficiency. The user can visualize a sectional view of pollution densities in the room, which can be simpler and sufficient to make a relevant decision (see figure 5). This idea of a partial representation on a surface or a volume delimited by the user could also allow the use of vector fields, which does not seem relevant in virtual reality in the case of a representation of the whole data (figure 6). Streamlines are very common in the literature, and we consider that they could be adapted to virtual reality, although they are in fact different from arrows representations and are closer to an artistic visualization. In the same way, the discretization of the space-time must be able to be modulated with the representation (i.e, without hindering the user), according to the need for precision, in order to consider the subtle variations in the flows. These new developments have not yet been evaluated by user testing, which constitutes our research perspectives.

One important limitation of this research is that the relevancy of design choices is strongly dependent of the specificities on the use cases. It depends on the parameters to represent, the use that will be made of the representation (training, decision-making, etc.) and presumably the technical knowledge of users. Because of these differences, it is difficult to have a general discourse on this subject. Future work on this topic would be to categorize the types of user profiles, data types and uses in order to build a set of targeted recommendations.

Figure 1 .

 1 Figure 1. Three different representations of air movement (direction and force)

Figure 2 .

 2 Figure 2. Means and standard deviations for the scores concerning the perception of the three pictures.

Figure 3 .

 3 Figure 3. Three different representations of dynamic air particles

Figure 4 .

 4 Figure 4. Means and standard deviations for the dimensions of the User Experience Questionnaire 29 concerning the perception of the three visualizations.

Figure 5 .

 5 Figure 5. Partial views of the data

Arts et M étiers Institute of Technology, France

Mann+Hummel, France

AMValor, France