A Pouya

M N Vu

& D Seyedi

Modelling effective permeability of fracture networks in permeable rock formation by singular integral equations method

Keywords: facture network, porous rock, steady-state flow, effective permeability

In this paper, theoretical and numerical formulations of plane steady-state fluid flow in a fractured porous rock are used to investigate its effective permeability. If the far field inflow is uniform, the theoretical solution shows that the pressure field in the matrix is a function of the discharge in the fracture network. A numerical resolution based on singular integral equations is employed to derive the general problem of many intersected fractures in order to obtain the pressure field in anisotropic matrix. This solution allows computing the flux in the fractures which is the key issue for upscalling the equivalent permeability. This paper presents in detail the method for deriving the equivalent permeability from this solution. This method is applied to two real cases: an Excavation Damage Zone (EDZ) around a deep underground gallery and a geological rock formation presenting several families of fractures. The results of the both cases show that the developed method provides an easy and efficient way to determine the equivalent permeability of the fractured porous rock medium. This equivalent permeability can be implemented in analytical and numerical tools for continuous media towards estimating the flow characteristics in the rock formation.

Introduction

In the context of important industrial applications such as underground nuclear waste disposal in claystone, geological CO 2 storage or hydrocabure reservoir, the numerical modelling of fluid flow in fractured porous materials is intensively investigated. More specifically, in damaged zones, cracks and geological faults generally have a great effect on the effective permeability of porous media or rock masses. The interest of this question also extends to other fields, such as hydrogeology for water resources management in aquifers, etc. Different approaches are used to study the problem of flow in porous media containing fractures with mass exchange between matrix and fracture, a review of which can be found in Sanchez-Vila et al. [START_REF] Sánchez-Vila | A synthesis of approaches to upscaling of hydraulic conductivities[END_REF], Goméz-Hernandez and Wen [START_REF] Goméz-Hernández | Upscaling hydraulic conductivities in heterogeneous media: an overview[END_REF] and Renard and de Marsily [START_REF] Renard | Calculating equivalent permeability: a review[END_REF]. However, these works are essentially based on some numerical-empirical modelling. Pouya and Fouché [START_REF] Pouya | Permeability of 3D discontinuity networks: new tensors from boundary-conditioned homogenization[END_REF] presented some rigorous theoretical-numerical results based on homogenisation theories for the equivalent permeability of heterogeneous or fractured media. The self-consistent scheme is investigated by Dormieux and Kondo [START_REF] Dormieux | Approche micromécanique du couplage perméabilité-endommagement[END_REF] and Barthélémy [START_REF] Barthélémy | Effective permeability of media with a dense network of long and micro fractures[END_REF] for upscalling the effective permeability of cracked porous media. In their method, cracks are assimilated to ellipsoidal inclusion obeying to a Darcy's law of flow. Following a different way, based on theory of Cauchy integral, Liolios and Exadaktylos [START_REF] Liolios | A solution of steady-state fluid flow in multiply fractured isotropic porous media[END_REF] studied mathematically the problem of plane steady-state flow in an infinite isotropic porous media with multiple non-intersecting fractures. In this approach, the fracture is considered as zero thickness discontinuity line and complex numbers are used. Recently, Pouya and Ghabezloo [START_REF] Pouya | Flow around a crack in a porous matrix and related problems[END_REF] used a direct approach to obtain a general potential pressure solution function of flux in the fracture network for anisotropic matrix containing intersecting curvilinear fractures. By continuing this work, Pouya and Vu [START_REF] Pouya | Modeling steady-state flow in fractured porous rocks by singular integral equations method[END_REF] employed the singular integral equation method to establish a fast numerical calculation that can be used to derive the general solution and to show the fluid flow around many intersecting fractures. Using the recent theoretical and numerical advances, this paper presents a method to evaluate the effective permeability of the permeable rock formation containing a dense fracture network. The numerical solution of pressure and discharge in the fracture, and then, the equivalent permeability using the formulas of average velocity and of average pressure gradient. The numerical calculation is applied to two cases. The first one is an EDZ around a deep underground gallery where the fractures are curved surfaces with a well defined shape and are almost regularly distributed in this zone. The second is a fractured geological formation with several families of fractures in which the fractures are plane surfaces and their characteristic parameters (orientation, length and position) are defined by a stochastic law. This approach provides an easy and efficient tool to determine the effective permeability of fractured porous rock formation which includes large applied areas.

An infinite homogeneous body Ω containing a fracture network is considered (fig. 1). Fracture number m is denoted Γ m . It is generally presented by a smooth curve z m of the curvilinear abscise s. The intersection points of two or more fractures and the extremities of fractures ending in the matrix constitute a set of singular point S. Fluid velocity v(x) in the matrix is given by Darcy's law:

Γ - Ω ∈ ∀x ( ) ( ). ( ) v x x p x = - ∇ k (1)
where k is the matrix permeability and p(x) the pressure field. The flow through the discontinuity is usually expressed in terms of the Poiseuille's law:

Γ ∈ ∀s p s c s q s ∂ - = ). ( ) ( (2)
where, s is the abscise along the fracture and c is the fracture's conductivity.

p ∞ (x) Ω ∂Ω Γ k(x), p(x), v(x)
Γ - Ω ∈ ∀x ∇.v(x) = 0 (3)
The fracture-matrix mass exchange law on the fracture excluding singular points is obtained by considering masse balance in a portion ds of the fractures (fig. 2a) At the intersection points of fractures (fig. 2b), Pouya and Ghabezloo [START_REF] Pouya | Flow around a crack in a porous matrix and related problems[END_REF] analyzed the masse balance in a small disc centred on the singular point and deduced the following relationship:

Γ ∈ ∀s [ ] [ ] 0 ) ( ) ( . ) ( = ∂ + s p s n z v s (4)
S z ∈ ∀ ( ) 0 ) ( . 0 = -       + ∇ ∑ z x q x v b b δ (5)
with b q 0 is the outgoing flow on the branch of fracture b and δ is the Dirac distribution. Nevertheless, a deeper mathematical analysis allows us to go further than eqn. ( 5) (Pouya and Vu [START_REF] Pouya | Modeling steady-state flow in fractured porous rocks by singular integral equations method[END_REF]) and shows that:

0 0 = ∑ b b q (6)
In order to calculate the effective permeability, it is sufficient to impose a constant pressure gradient at infinity:

x A x p . ) ( = ∞
. The previous problem was derived rigorously by Pouya and Ghabezloo [START_REF] Pouya | Flow around a crack in a porous matrix and related problems[END_REF] and a general solution of potential was given:

m m m m 2 m m 1 ( ) ( ) ( ) ( ) . . ( ) 2 . ( ) -1 -1 k
k x z s p x p x q s t s d s x z s πκ ∞ Γ - = +   -   ∑ ∫ (7)
where, x is a current point in the material, z m (s) is the point the Γ m at the curvilinear abscise s, ) (s q m and ) (s t m are respectively the discharge along the fracture and the unit vector tangent to the fracture at this point, and κ is the square root of determinant of k: κ = k .

Numerical calculation

When field point x is located in the fractures; we obtain a singular integral equation of unknown pressure field in the fracture networks by replacing eqn. (2) into eqn. [START_REF] Goméz-Hernández | Upscaling hydraulic conductivities in heterogeneous media: an overview[END_REF]. This equation was resolved numerically by collocation method (Bonnet [START_REF] Bonnet | Equations intégrales et éléments de frontière[END_REF]) which consists of enforcing exactly eqn. [START_REF] Goméz-Hernández | Upscaling hydraulic conductivities in heterogeneous media: an overview[END_REF] at finite number of points called collocation point. Curvilinear fractures are approximated by a series of small linear segments and are discretized by E elements numbered n and denoted E n (1≤n≤E). Thus, eqn.( 7) is written as follows:

p(x) ∑ = ∞ + = N n n I x p 1 ) ( (8) with: n n n n 2 n ( ) ( ) ( ) . . . ( )
- =   -   ∫ n E x z s I x q s t ds x z s -1 -1 k k (9)
Two types of elements are distinguished: extremity elements and current elements. Linear interpolation of pressure is used for current elements. General theoretical results of velocity field singularity around a fracture tip lead us to choose an interpolation function for the flux q(s) corresponding to a variation as s 1/2 where s is the distance on the fracture line to the extremity point. Therefore, interpolation function for pressure on the left extremity elements (tip at s=0), current element, and right extremity elements (tip at s=L) are respectively expressed as following:

3 2 1 1 3 2 1 1 3 2 1 2 3 ( ) ( ) ( ) p p p s p s L p p p s p s L p p p s p L s L - = + - = + - = - - (10)
At first, a collocation point x is selected per element which verifies I n (x) = 0. And then, the set of extremity points of fracture network is added as supplementary collocation points. It should be noted that no matter how the fractures intersect and how the discretization is made, the number of nodes N is always smaller than number of collocation points M. Otherwise, this method of selection of collocation points always leads to the number of equations M larger than the number of nodal pressure unknown. It is interesting that all elementary integrals I n (x) in eqn. ( 9) can be calculated analytically by using the variable interpolations in eqn. [START_REF] Maleki | Modélisation numérique du couplage entre l'endommagement et la perméabilité des argiles raides[END_REF], which provides a very fast calculation method. Enforcing eqn. [START_REF] Goméz-Hernández | Upscaling hydraulic conductivities in heterogeneous media: an overview[END_REF] in M collocation points leads to following matrix equation: H.P=Y [START_REF] Pouya | Permeability of 3D discontinuity networks: new tensors from boundary-conditioned homogenization[END_REF] where, P(Nx1) is the column of nodal pressure unknown, Y(Mx1) is the column of infinite pressure field at collocation points and H(MxN) is computed from the assembly operation after computing all elementary integrals. The approximated solution can be found by the least squares fitting method. Pouya and Vu [START_REF] Pouya | Modeling steady-state flow in fractured porous rocks by singular integral equations method[END_REF] have well shown the validity of this method by comparing the numerical solution of a single straight fracture in an infinite homogeneous media with the known closed-form solution of the same problem. Once the eqn. ( 11) is solved, the nodal pressure P{p n } (1≤n≤N) are known, using eqn. (2), we can deduce the flux in each element. In addition, in the postprocessing stage, eqn. [START_REF] Pouya | Permeability of 3D discontinuity networks: new tensors from boundary-conditioned homogenization[END_REF] shall be employed again to compute directly the field pressure in porous matrix. For a homogeneous rock mass Ω, containing a fractures network Γ m that is considered as the limit case of thin permeable layers, the average velocity and pressure gradient in Ω is given, respectively, by:

        + Ω Ω ≡ ∑ ∫ ∫ Γ Ω m m ds t q d v V 1 ; Ω ∇ Ω ≡ ∫ Ω d p G 1 (12)
Let us to suppose that a linear pressure x A x p . ) ( = is applied at the infinity. The pressure field in matrix is known, especially on the boundary Ω ∂ , hence, the vector G can be deduced such as:

1 . G p nds ∂Ω ≡ Ω ∫ (13)
where, n is the unit outward normal on Ω ∂ . Besides, Pouya and Fouché [11] demonstrated that 1 .

k vd G Ω Ω = - Ω ∫
. Hence, the linearity of all velocities and flux with respect to G implies that there exists a tensor k f that satisfies the following relationship:

1 . k m f m qtds G Γ = - Ω ∑ ∫ (14)
Thus, under boundary condition of linear pressure, the global permeability tensor of Ω domain is given by k p =k+k f . In the case of weak fracture density, the Mori-Tanaka estimation of k p is determined by neglecting the interaction between fractures i.e. assuming that the flux q in the left-hand side of eqn. ( 14) is deduced from the results obtained for a single fracture in an infinite body which is given in Pouya and Ghabezloo [START_REF] Pouya | Flow around a crack in a porous matrix and related problems[END_REF]. However, in the case of high fracture density (or dense fracture networks), the integral of flux on the fracture network, i.e. the term described in left-hand side of eqn. ( 13) could be computed numerically for a fracture network discretized by E element.

In the next section, we shall take up some real examples for calculating the effective permeability in using the recent advance of theoretical development and numerical tools which are described above.

Application

Effective permeability of EDZ

Excavation of underground galleries generates cracks and fractures in a zone around the gallery called Excavation Damaged Zone (EDZ). The effect of fractures on the hydraulic properties of EDZ can be estimated through a numerical method taking into account the complex geometry and density of discontinuities. Three zones around the gallery (fig. 3) can be distinguished: the zone crossed by fractures; damaged or micro-fractured zone beyond fractured zone and the zone of intact rock beyond EDZ. The first zone is considered here to calculate its equivalent permeability. The fractures appearing in this zone has more or less regular geometry with the same shape, spacing, extensions that are repeated along the axis of gallery.

X Y a b Z λ + = (15)
The parameters a, b, λ depend on the orientation of the gallery. And then, the values: a=b=4, λ=0.57 are determined by comparing the trace of theoretical surface to observed fracture in field. In addition, the last parameter d=60 cm of spacing between fractures is fixed by the field measurement.

The above geometrical model is used to evaluate the effective permeability of EDZ (fig. 5) with the isotropic matrix permeability k = 10 -12 m/s and the hydraulic conductivity of fracture c=10 -11 m 2 /s. The chevron fracture form is discretized by series of small linear segments. The linear pressure x A x p . ) ( = is imposed at infinity. By solving the problem for two vectors of pressure gradient A 1 =(1,0) and A 2 =(0,1), two average velocity vectors 1 V , 2 V and two average pressure gradient vectors 1 G , 2 G are obtained. The equivalent permeability k p is then deduced from the equation V = -k p .G. Launching the numerical program for the EDZ zone, the following effective permeability matrix is found:

k p - -   =     x 10 -12 (m/s)
This result shows that the effective permeability matrix is almost symmetric according to the results announced by Pouya and Fouché [START_REF] Pouya | Permeability of 3D discontinuity networks: new tensors from boundary-conditioned homogenization[END_REF]. The principal directions of this permeability tensor can be calculated easily. They are found to be approximately parallel and orthogonal to the fractures surfaces. Otherwise, permeability predicted in this way is smaller than that given by the Mori-Tanaka estimation based on theoretical results of Pouya and Ghabezloo [START_REF] Pouya | Flow around a crack in a porous matrix and related problems[END_REF]. The difference can be explained by the fact that the Mori-Tanaka estimation does not take into account the fractures interaction.

Fracture network in a permeable rock

The study of fluid flow in fractured permeable rocks needs, first, a geometrical model of fractures that may be described the deterministic models (Kolditz and Clauser [START_REF] Kolditz | Numerical Simulation of Flow and Heat Transfer in Fractured Crystalline Rocks: Application to The Hot Dry Rock Site in Rosemanowes (U.K.)[END_REF]) or by stochastic approach (Cacas et al. [START_REF] Cacas | Modeling Fracture Flow with A Stochastic Discrete Fracture Network: Calibration and Validation: 1. The Flow Model[END_REF]; Billaux [START_REF] Billaux | Hydrogéologie des milieux fracturés. Géométrie, connectivité et comportement hydraulique[END_REF]; Gervais [START_REF] Gervais | Modélisation géométrique d'un réseau de fractures dans un massif rocheux stratifié. Application aux carrières marbrières de Comblanchien (Côte d'Or[END_REF] and Maleki and Pouya [START_REF] Maleki | Modélisation numérique du couplage entre l'endommagement et la perméabilité des argiles raides[END_REF]). The stochastic simulation of joint network is based on a hierarchical probabilistic model that realistically reproduces fracture connectivity using minimal data such as the number of fracture sets, the fracture length, spacing and density. In practice, the stochastic distribution is often used for modelling fluid flow in hydrogeology or for natural geological reservoir since their large flexibility. In this work, the study introduced by Maleki and Pouya [START_REF] Maleki | Modélisation numérique du couplage entre l'endommagement et la perméabilité des argiles raides[END_REF] is adopted. A square domain is characterized by two corners (X min ,Y min ) and (X max ,Y max ). Three fracture families are defined each one by a density ρ (number of fracture per area unit), the half of fracture length r = L/2 and the orientation θ (angle between fracture line and x-axis, 0 ≤ θ < π). The fracture line is completely determined by these parameters and the coordinates (x 0 ,y 0 ) of its center. A uniform distribution of the fracture center points and of the angle θ; and an exponential law of length distribution are supposed.

In order to study the effect of fracture density on the effective permeability of fractured porous media, three families are generated in a domain 10x10m 2 of rock formation. Several configurations are considered for which all the fractures parameters are constant except for the density that increases (fig. 6). Herein, we introduce the dimensionless parameter The equivalent permeability increases first linearly with the fracture density ρ, for small values of the density. A sharp increase is observed for a critical density ρ≈1.5 that represents the percolation threshold. The dimensionless density R ρ = π ρ , where R is the mean half-length of the fractures, is found to be around 1.2. After percolation threshold, the effective permeability function of density increases with a different pace as observed also by Maleki and Pouya [START_REF] Maleki | Modélisation numérique du couplage entre l'endommagement et la perméabilité des argiles raides[END_REF].

Conclusion

The effective permeability of fractured porous rocks is studied numerically using analytical and numerical developments based on the singular integral equations. This method that is clearly presented in this paper allows taking into account the real geometry of the fractures and of their interactions effects. The geometrical shape of individual fractures may be introduced precisely in the model as it has been done for the case of the EDZ studied here above. In the case of great number of fractures, the geometry of the fractures network can be described and generated numerically in the model by stochastic laws. The parameters of these stochastic laws are deduced from observation data. The numerical program that has been developed using this method provides an efficient tool for quick evaluation of the effective permeability of fractured porous rocks. The results presented in this paper to illustrate this method show also the important effect of percolation which modifies significantly the magnitude of effective permeability.

Figure 1 :

 1 Figure 1: Rock mass containing fractures Mass conservation in the matrix reads:

Figure 2 :

 2 Figure 2: a. Mass exchange between matrix and fracture. b Mass balance in a disc around an intersection point (Pouya and Ghabezloo [14])

Figure 3 :Figure 4 :

 34 Figure 3: Representation of the EDZ around the gallery and the shape of the fracture (ANDRA 2008 [18])

Figure 5 :

 5 Figure 5: Representative segment of EDZ in calculating the effective permeability

Figure 6 :

 6 and then we use the numerical procedure to compute the equivalent permeability of fractured porous rock with homogeneous permeability k for matrix formation. Fracture network in porous rock formation: (a) non-percolated network, (b) percolated network.

Figure 7 Figure 7 :

 77 Figure 7 depicts the evolution of the determinant of the effective permeability tensor eff p K = k with the fracture density for different value of λ .

Acknowledgement

This work is co-funded by IFSTTAR and BRGM.


Piotr Biler

AND Cyril Imbert
email: cyril.imbert@u-pec.fr

Grzegorz Karch
email: grzegorz.karch@math.uni.wroc.pl

NONLOCAL POROUS MEDIUM EQUATION: BARENBLATT PROFILES AND OTHER WEAK SOLUTIONS

Keywords: . 2000 Mathematics Subject Classification. 35K55, 35B45, 35C06 porous medium equation, nonlocal equation, hypercontractivity, self-similar solutions

A degenerate nonlinear nonlocal evolution equation is considered; it can be understood as a porous medium equation whose pressure law is nonlinear and nonlocal. We show the existence of sign changing weak solutions to the corresponding Cauchy problem. Moreover, we construct explicit compactly supported selfsimilar solutions which generalize Barenblatt profiles -the well-known solutions of the classical porous medium equation.

Introduction

In this work, we study the following degenerate nonlinear nonlocal evolution equation (1.1)

∂ t u = ∇ • |u|∇ α-1 (|u| m-2 u) , x ∈ R d , t > 0,
where m > 1 and ∇ α-1 denotes the integro-differential operator ∇(-∆) for merely integrable initial data, and we prove that these solutions satisfy sharp hypercontractivity L 1 → L p estimates.

A nonlocal operator. Equation (1.1) involves a nonlocal operator denoted by ∇ α-1 which can be defined as the Fourier multiplier whose symbol is iξ|ξ| α-2 . This notation emphasizes that it is a (pseudo-differential)

operator of order α -1. Recalling the definition of the fractional Laplace operator (-∆)

α 2 (v) = F -1 (|ξ| α F v)
and the Riesz potential I β = (-∆) -β 2 , i.e. Fourier multipliers whose symbols are |ξ| α and |ξ| -β respectively (see for instance [23, Ch. V]), the fractional gradient ∇ α-1 can also be written as ∇I 2-α . Finally, let us emphasize that the definition of ∇ α-1 is consistent with the usual gradient: ∇ 1 = ∇; the components of ∇ 0 are the Riesz transforms; moreover we have

∇ • ∇ α-1 = ∇ α 2 • ∇ α 2 = -(-∆) α 2 .
It is also possible, following the reasoning from [START_REF] Droniou | Fractal first order partial differential equations[END_REF]Th. 1] to define the fractional gradient via the singular integral formula for smooth and bounded functions

v : R d → R (1.3) ∇ α-1 v(x) = C d,α v(x) -v(x + z) z |z| d+α dz with a suitable constant C d,α > 0.
Related equations and results. Our preliminary results for the problem (1.1)-(1.2) have been announced in [START_REF] Biler | Barenblatt profiles for a nonlocal porous medium equation[END_REF].

First, we would like to shed light on the link between (1.1) and other partial differential equations. Notice that when α = 2 equation (1.1) coincides with the classical (nonlinear parabolic) porous medium equation (1.4)

∂ t u = ∇ • |u|∇(|u| m-2 u) = ∇ • (m -1)|u| m-1 ∇u .
For the theory of porous media equations, the interested reader is referred to [START_REF] Vázquez | Smoothing and decay estimates for nonlinear diffusion equations[END_REF][START_REF]The porous medium equation[END_REF] and references therein. Of course, for m = 2, the Boussinesq equation is recovered.

The following nonlinear and nonlocal equation

(1.5) ∂ t v + |v x | - ∂ 2 ∂x 2 α 2 v = 0
in the one-dimensional case x ∈ R was studied by the first, the third authors and R. Monneau [START_REF] Biler | Nonlinear diffusion of dislocation density and self-similar solutions[END_REF]. Such an equation was derived as a model for the dynamics of dislocations in crystals. In [START_REF] Biler | Nonlinear diffusion of dislocation density and self-similar solutions[END_REF], the existence, uniqueness and comparison properties of (viscosity) solutions have been proved, and explicit self-similar solutions have been constructed. Notice that the function u = v x , where v is a solution to (1.5), solves the one-dimensional case of (1.1) with m = 2. Thus, equation (1.1) is a multidimensional generalization of the one in (1.5).

Recently, Caffarelli and Vázquez [START_REF] Caffarelli | Nonlinear porous medium flow with fractional potential pressure[END_REF][START_REF] Caffarelli | Asymptotic behaviour of a porous medium equation with fractional diffusion[END_REF] studied nonnegative weak solutions of (1.1) in the case m = 2 in the multidimensional case. Precisely, they studied the following (nonlocal) porous medium equation in R d

(1.6)

∂ t u = ∇ • (u∇p),
with the nonlocal pressure law p = (-∆) -s u, 0 < s < 1, obtained from the density u ≥ 0. Notice that, for α = 2 -2s ∈ (0, 2), equation (1.6) reads ∂ t u = ∇ • (u∇ α-1 u). For sign changing u's, our equation (1.1) is a (formally parabolic) extension of equation (1.6) of the structure of (1.5). In [START_REF] Caffarelli | Nonlinear porous medium flow with fractional potential pressure[END_REF], Caffarelli and Vázquez constructed nonnegative weak solutions for (1.6), i.e. for (1.1) with m = 2, with initial data satisfying:

u 0 ∈ L 1 (R d ) ∩ L ∞ (R d
) and such that 0 ≤ u 0 (x) ≤ Ae -a|x| for some A, a > 0. Besides the positivity and the mass preservation, the properties of solutions, listed in the next paper [7, p. 4], include the finite speed of propagation proved using the comparison with suitable supersolutions. Further regularity properties of solutions of (1.1) with m = 2 and α ∈ (0, 1) are studied in [START_REF] Caffarelli | Regularity of solutions of the fractional porous medium flow[END_REF].

Another nonlocal porous medium equation has been proposed in [START_REF] De Pablo | A fractional porous medium equation[END_REF][START_REF]A general fractional porous medium equation[END_REF][START_REF] Vázquez | Barenblatt solutions and asymptotic behaviour for a nonlinear fractional heat equation of porous medium type[END_REF]]

∂ t u + (-∆) α 2 (|u| m-1 u) = 0
for α = 1 and α ∈ (0, 2), respectively. Among several other properties like smoothing effects and decay estimates, solutions of this generalization of the porous medium equation enjoy the L 1 -contraction property, so, they are unique. But self-similar solutions are not compactly supported [START_REF] Vázquez | Barenblatt solutions and asymptotic behaviour for a nonlinear fractional heat equation of porous medium type[END_REF]Th. 1.1].

Finally, we recall that the following nonlocal higher order equation, appearing in the modeling of propagation of fractures in rocks,

∂ t u = ∇ • (u n ∇(-∆) 1 2 u)
(with u ≥ 0 and n > 1), has been studied in [START_REF] Imbert | Existence of solutions for a higher order non-local equation appearing in crack dynamics[END_REF] in a one-dimensional bounded domain. At least formally, this equation with n = 1 corresponds to (1.1) with α = 3 and m = 2.

Notation. In this work,

Q T denotes (0, T ) × R d . The usual norm of the Lebesgue space L p (R d ) is denoted by • p for any p ∈ [1, ∞],
and H s,p (R d ) with the norm • H s,p is the fractional order Sobolev space, see Section 3. The Fourier transform F and its inverse transform

F -1 of a function v ∈ L 1 (R d ) are defined by F v(ξ) = (2π) -d 2 v(x) e -ix•ξ dx, F -1 v(x) = (2π) -d 2 v(x) e ix•ξ dξ.
Here, all integrals with no integration limits are over the whole space R d if one integrates with respect to x and over the whole half-line R + = [0, ∞) if the integration is with respect to t. As usual, w + = max{0, w}, w -= max{0, -w}, so w = w + -w -. Constants (always independent of x and t) will be denoted by the same letter C, even if they may vary from line to line. Sometimes we write, e.g., C = C(p, q, r) when we want to emphasize the dependence of C on particular parameters p, q, r, for instance.

Main results

In this work, we show two main results: we construct explicit self-similar solutions of equation (1.1), as well as we prove that the initial value problem (1.1)-(1.2) has a global-in-time weak solution which satisfies certain optimal decay estimates.

We first "recall" the appropriate notion of weak solutions for Equation (1.1), see for instance [START_REF] Vázquez | Smoothing and decay estimates for nonlinear diffusion equations[END_REF][START_REF]The porous medium equation[END_REF].

Definition 2.1 (Weak solutions). A function u : Q T → R is a weak solution of the problem (1.1)-(1.2) in Q T if u ∈ L 1 (Q T ), ∇ α-1 (|u| m-2 u) ∈ L 1 loc (Q T ) and |u|∇ α-1 (|u| m-2 u) ∈ L 1 loc (Q T ), and
u∂ t ϕ -|u|∇ α-1 (|u| m-2 u) • ∇ϕ dt dx + u 0 (x)ϕ(0, x) dx = 0 for all test functions ϕ ∈ C ∞ (Q T ) ∩ C(Q T )
such that ϕ has a compact support in the space variable x and vanishes near t = T .

The first main result of this work says that there is a family of nonnegative explicit compactly supported self-similar solutions of (1.1), i.e. nonnegative solutions that are invariant under a suitable scaling. Observe that if u(t, x) is a solution of (1.1), then so is L dλ u(Lt, L λ x) for each L > 0, where λ = (d(m -1) + α) -1 . Thus, the scale invariant solutions should be of the following form

(2.1) u(t, x) = 1 t dλ Φ x t λ with λ = 1 d(m -1) + α ,
for some function Φ : R d → R satisfying the following nonlocal "elliptic type" equation

(2.2) -λ∇ • (yΦ) = ∇ • (|Φ|∇ α-1 (|Φ| m-2 Φ)) where y = x t λ .
Theorem 2.2 (Self-similar solutions). Let α ∈ (0, 2], m > 1. Consider the function Φ α,m : R d → R defined as

(2.3) Φ α,m (y) = k α,d (1 -|y| 2 ) α 2 + 1 m-1 with the constant k α,d = d Γ d 2 (d(m -1) + α)2 α Γ 1 + α 2 Γ d+α 2 .
Then, the function u Zel'dovich solutions of the porous medium equation (1.4), see for instance [START_REF]The porous medium equation[END_REF][START_REF] Vázquez | Smoothing and decay estimates for nonlinear diffusion equations[END_REF].

: (0, ∞) × R d → R + defined by (2.1) with Φ = Φ α,m is a weak solution of (1.1) in the sense of Definition 2.1 in Q η,T ≡ (η, T ) × R d for every 0 < η < T < ∞. Moreover, u(t,
Remark 2.4. For each M ∈ (0, ∞) we can find a nonnegative self-similar solution u with prescribed mass M ≡ u(t, x) dx (which is conserved in time) by a suitable scaling of the profile Φ α,m . Indeed, this self-similar solution is given by the formula

(2.4) u(t, x) = t - d d(m-1)+α k α,d R 2 -xt - 1 d(m-1)+α 2 α 2 + 1 m-1
, where, for each M > 0, there exists a unique R > 0 such that u(t, x) dx = M . Remark 2.5. Self-similar solutions of equation (1.6) (which is a particular case of equation (1.1)) have been proved to exist in [START_REF] Caffarelli | Asymptotic behaviour of a porous medium equation with fractional diffusion[END_REF] by studying the following obstacle problem for the fractional Laplacian. For α ∈ (0, 2) and Ψ(y) = C -a|y| 2 where a = a(d, α) and C > 0, one looks for a function P = P (y) with the following properties:

P ≥ Ψ, ( -∆)
α 2 P ≥ 0, and either P = Ψ or (-∆)

α 2 P = 0.
The novelty of our approach is that we exhibit the explicit self-similar profile Φ α,2 defined in (2.3) and, consequently, the explicit solution of this obstacle problem: P (y) = I α (Φ α,2 ) y R , where I α = (-∆) -α 2 is the Riesz potential and R > 0 is a suitable constant.

Next, we prove the existence of weak solutions to the initial value problem (1.1)-(1.2).

Theorem 2.6 (Existence and decay of L p -norms). Let α ∈ (0, 2) and Remark 2.8. After proving those hypercontractivity estimates in [START_REF] Biler | Barenblatt profiles for a nonlocal porous medium equation[END_REF], we learned that a similar result is obtained in [START_REF] Caffarelli | Regularity of solutions of the fractional porous medium flow[END_REF], however, for a less general model: α ∈ (0, 1), m = 2, and nonnegative u 0 . Moreover, analogous decay estimates for another fractional porous medium equation of the form ∂ t u + (-∆) α 2 (|u| m-1 u) = 0 were proved recently in [START_REF]A general fractional porous medium equation[END_REF].

(2.5) m > 1 + 1-α d if α ∈ (0, 1], m > 3 -2 α if α ∈ (1, 2). Given u 0 ∈ L 1 (R d ),
Compared with the methods used in [START_REF] Caffarelli | Nonlinear porous medium flow with fractional potential pressure[END_REF], we propose an alternative strategy of the proof of the existence of solutions. In this paper, we consider approximating solutions u = u δ,ε of the equation

∂ t u = δ∆u + ∇ • (|u|∇ α-1 G ε (u)),
considered in the whole space R d , where G ε (u) is a sufficiently smooth approximation of u|u| m-2 , and then we pass to the limit with the parameters ε ց 0, and δ ց 0. Solutions of the approximating equation exist because the parabolic regularization term δ∆u is strong enough to regularize equation (1.1) when 0 < α < 2, but of course not for α = 2. Our approach resembles the approach to the one-dimensional model achieved in [4, Sec.

4 and 5] via viscosity solutions.

Preliminaries

In this section, we collect known results that we will used in proofs of the main theorems.

Bessel and hypergeometric functions. Bessel functions of order ν are denoted by J ν (z), and they behave for small and large values of the (complex) variable z like

J ν (z) ∼ 1 Γ(ν+1) ( z 2 ) ν as z → 0, J ν (z) ∼ ( 2 π ) 1 2 cos(z -νπ 2 -π 4 )z -1 2 as |z| → ∞,
where for functions f , g, the relation f ∼ g means that f g → 1. For the proofs of those properties of J ν , the reader is referred to, e.g., [START_REF] Watson | A treatise on the theory of Bessel functions[END_REF].

The hypergeometric function, denoted by 2 F 1 (a, b; c; z), is defined for complex numbers a, b, c and z as the sum of the series

2 F 1 (a, b; c; z) = ∞ n=0 (a) n (b) n (c) n n! z n for |z| < 1,
where (a) n ≡ Γ(a+n) Γ(a) , and Γ denotes the Euler Gamma function. This series is absolutely convergent in the open unit disc and also on the circle

|z| = 1 if ℜ(a + b -c) < 0.
It is known [19, p. 39] that when b = -n is a negative integer, 2 F 1 (a, -n; c; z) is a polynomial function of degree n. In particular, we have

(3.1) 2 F 1 (a, -1; c; z) = 1 - a c z.
We will also use the following differentiation formula [19, p. 41] (3.2)

d dz 2 F 1 (a, b; c; z) = ab c 2 F 1 (a + 1, b + 1; c + 1; z).
The Weber-Schafheitlin integral. If 0 < b < a and if integral (3.3) below is convergent, then the following identity holds true

(3.3) ∞ 0 t -λ J µ (at)J ν (bt) dt = b ν 2 -λ a λ-ν-1 Γ( ν+µ-λ+1 2 ) Γ( -ν+µ+λ+1 2 )Γ(1 + ν) 2 F 1 ν + µ -λ + 1 2 , ν -µ -λ + 1 2 ; ν + 1; b 2 a 2 .
According to Watson, [30, pp. 401-403], this result was obtained by Sonine and Schafheitlin. However, it is usually referred to as the Weber-Schafheitlin discontinuous integral since there occurs a discontinuity for a = b.

The Stroock-Varopoulos inequality. We next recall the the Stroock-Varopoulos inequality, see [18, Theorem 2.1 and Condition (1.7)] for a proof.

Proposition 3.1. For α ∈ (0, 2], w ∈ C ∞ c (R d
) and q > 1, the following inequality holds true

(3.4) sgn w |w| q-1 (-∆) α 2 w dx ≥ 4(q -1) q 2 ∇ α 2 sgn w |w| q 2 2 dx ≥ 4(q -1) q 2 ∇ α 2 |w| q 2 2 dx.
Fractional order Sobolev spaces. The fractional order Sobolev spaces are defined as

H s,p (R d ) = {v ∈ L p (R d ) : ∇ s v ∈ L p (R d )} = {v ∈ L p (R d ) : (I -∆) s 2 v ∈ L p (R d )},
here with p ∈ (1, ∞), supplemented with the usual norm denoted by • H s,p , and we refer the reader to the books [START_REF] Taylor | Tools for PDE[END_REF][START_REF]Partial differential equations[END_REF] for properties of those spaces. In particular for s = α -1 with α ∈ (1, 2), the following well-known continuous embedding will be used repeatedly

(3.5) H α-1,p (R d ) ⊂ L ∞ (R d ) provided p > d α -1 (> 1).
We also recall the fractional integration theorem [23, Ch. V, §1.2]: the Riesz potential

I s = (-∆) -s 2 satisfies (3.6) I s u q ≤ C(p, q, s) u p for all s ∈ (0, d) and p, q ∈ (1, ∞) satisfying 1 q = 1 p -s d .
Some functional inequalities. We will use the following Nash inequality Moreover, we will use the following Gagliardo-Nirenberg type inequality Lemma 3.2. For p > 1 and p ≥ m -1, the inequality

(3.7) v 2 1+ α d 2 ≤ C N ∇ α 2 v 2 2 v 2α d 1 valid for all functions v ∈ L 1 (R d ), such that ∇ α 2 v ∈ L 2 (R d ),
(3.8) u a p ≤ C N ∇ α 2 |u| r 2 2 2 u b 1 holds with (3.9) a = p p -1 d(r -1) + α d , b = a -r = d(m -1) + pα d(p -1) , r = p + m -1.
Proof. This inequality is a consequence of the Nash inequality (3.7) written for v = |u| r 2 , i.e.

(3.10) u r 1+ α d r ≤ C N ∇ α 2 |u| r 2 2 2 u rα d r 2
, and two Hölder inequalities

u p ≤ u γ r u 1-γ 1 with γ = p -1 r -1 r p , and
u r 2 ≤ u δ p u 1-δ 1 with δ = r -2 p -1 p r .
Combining the above three inequalities, we get (3.8).

Proof of Theorem 2.2

This section is devoted to the study of nonnegative self-similar solutions for (1.1) with m > 1. As explained above, this problem reduces to a study of the elliptic-like equation (2.2) which for nonnegative Φ takes the form

-λy Φ = Φ∇ α-1 (Φ m-1 ).
Moreover, since we want to construct compactly supported solutions, we are interested in solutions Φ vanishing outside the unit ball B 1 . This is the reason why we consider the Dirichlet problem

-λy = ∇ α-1 (Φ m-1 ) in B 1 , Φ = 0 in R d \ B 1 . (4.1)
It is well known that, in the case of nonlocal operators (such as ∇ α-1 ), the homogeneous Dirichlet condition should be understood in the form Φ ≡ 0 outside the domain B 1 , and not only Φ = 0 on the boundary ∂B 1 .

The reader is referred to, e.g., [START_REF] Barles | On the Dirichlet problem for second-order elliptic integro-differential equations[END_REF] for more explanations.

We claim that the proof of Theorem 2.2 reduces to the following key computation. Here, 2 F 1 denotes the classical hypergeometric function defined in Section 3.

Lemma 4.1. For all β ∈ (0, 2), β < d, and γ > 0, we have

(4.2) I β (1 -|y| 2 ) γ 2 + =    C γ,β,d × 2 F 1 d-β 2 , -γ+β 2 ; d 2 ; |y| 2 for |y| ≤ 1, Cγ,β,d |y| β-d × 2 F 1 d-β 2 , 2-β 2 ; d+γ 2 ; 1 |y| 2 for |y| > 1, with C γ,β,d = 2 -β Γ( γ 2 +1)Γ( d-β 2 ) Γ( d 2 )Γ( β+γ 2 +1) and Cγ,β,d = 2 -β Γ( γ 2 +1)Γ( d-β 2 ) Γ( d 4 )Γ( d+γ 2 +1) .
The proof of this lemma is postponed to the end of this section.

The following corollary is an immediate consequence of Lemma 4.1 with 2 -β = α = γ, of the property of 2 F 1 formulated in (3.1), and of the identity (-∆)

α 2 = (-∆)I 2-α . It has an important probabilistic interpretation,
and recently, related results and generalizations have been proved in [START_REF] Dyda | Fractional calculus for power functions and eigenvalues of the fractional laplacian[END_REF].

Corollary 4.2 (Getoor [14, Th. 5.2]). For all α ∈ (0, 2], the identity

K α,d (-∆) α 2 (1 -|y| 2 ) α 2 + = 1 in B 1 holds true with the constant K α,d = Γ( d 2 ) 2 α Γ(1+ α 2 )Γ( d+α 2 )
.

Before proving Lemma 4.1, we first use it to derive Theorem 2.2.

Proof of Theorem 2.2. We check that u

(t, x) = t -dλ Φ α,m (t -λ x) is a weak solution of (1.1) in the sense of Definition 2.1. First, u ∈ L 1 (Q T ) if and only if Φ α,m ∈ L 1 (R d ), which is obviously true. For later use, it is convenient to introduce the function Φ α (u) = 1 -|y| 2 α 2 + , so that k α,d Φ α = Φ m-1 α,m .
The fact that, for all η, T such that 0 < η < T , u∇ α-1 (u m-1 ) and ∇ α-1 (u m-1 ) are locally integrable in (η, T ) × R d follows from

I 2-α (Φ α ) ∈ H 1,1 loc (R d ),
which we prove by computing I 2-α (Φ α ). In order to do so, we first assume that α > 2 -d, and we apply Lemma 4.1 with γ = α ∈ (0, 2) and β = 2 -α, we use equation (3.1), and we get

(4.3) I 2-α (Φ α )(y) = C α,2-α,d 1 -d+α-2 d |y| 2 if |y| ≤ 1, Cα,2-α,d |y| 2-(d+α) 2 F 1 d+α 2 -1, α 2 ; d+α 2 ; 1 |y| 2 if |y| > 1.
The right-hand side of equation ( 4.3) defines a locally integrable function because

a + b -c = ( d + α 2 -1) + α 2 - d + α 2 < 0.
We then deduce that

∇ α-1 (Φ α )(y) = ∇I 2-α (Φ α )(y) = - λ k α,d y for y ∈ B 1 .
Note also that ∇ α-1 (Φ α ) can be computed outside B 1 thanks to the differentiation formula (3.2).

We now remark that

Φ α,m (y) ∇ α-1 Φ m-1 α,m (y) = -λy Φ α,m (y) for all y ∈ R d , which is in L 1 (R d ).
Moreover, the following equalities hold true in the sense of distributions in Q T ,

∂ t u(t, x) = -λt -dλ-1 ∇ y • (yΦ α,m )(t -λ x), ∇ x • (u∇ α-1 (|u| m-1 ))(t, x) = t -dλ-1 ∇ y • (Φ α,m ∇ α-1 Φ m-1 α,m )(t -λ x).
This allows us to conclude that u is indeed a weak solution of (1.1) in (η,

T ) × R d for all 0 < η < T < ∞ if α > 2 -d. Assume now that 0 < α ≤ 2 -d,
F (Φ γ )(ξ) = 2 γ 2 Γ γ 2 + 1 1 |ξ| d+γ 2 J d+γ 2 (|ξ|).
Since I β is the Fourier multiplier of symbol |ξ| -β , I β (Φ γ ) is the (inverse) Fourier transform of the following radially symmetric function

2 γ 2 Γ γ 2 + 1 1 |ξ| d+γ 2 +β J d+γ 2 (|ξ|).
We recall that by properties of Bessel functions collected in Section 3, we have J ν (r) = O (r ν ) as r → 0 and

J ν (r) = O r -1 2
as r → ∞. We see that the previous function is integrable since β < d and d < γ + 2β + 1.

Thanks to [24, Th. 3.3], we get

I β (Φ γ )(y) = 2 γ 2 Γ γ 2 + 1 |y| 1-d 2 ∞ 0 1 t d+γ 2 +β J d+γ 2 (t) t d 2 J d 2 -1 (t|y|) dt = 2 γ 2 Γ γ 2 + 1 |y| 1-d 2 ∞ 0 t -( γ 2 +β) J d+γ 2 (t)J d 2 -1 (t|y|) dt.
We obtain (4.2) applying (3.3) with the following choice of parameters:

• if |y| ≤ 1, we put λ = γ 2 + β, µ = d+γ 2 , ν = d 2 -1, a = 1 and b = |y|, • if |y| > 1, we put λ = γ 2 + β, µ = d 2 -1, ν = d+γ 2 , a = |y| and b = 1.

A regularized problem

In order to construct weak solutions of (1.1) for general initial data, we first consider the following regularized problem (5.1)

∂ t u = δ∆u + ∇ • |u|∇ α-1 (G(u)) , u(0, x) = u 0 (x),
where G : R → R satisfies (5.2)

   G differentiable and increasing, G(0) = G ′ (0) = 0, G ′ locally Lipschitz continuous.
Remark that for m ≥ 3 or m = 2, the function G(u) = |u| m-2 u satisfies (5.2). For m ∈ (1, 3), we consider the

following approximation G = G ε of |u| m-2 u G ε (u) = sgn u u 2 + ε 2 m-1 2 -ε m-1
with ε > 0. The following theorem holds true for a general function G satisfying (5.2).

Theorem 5.1 (Existence of solutions to the regularized problem). Let δ > 0 and assume that G is an arbitrary function satisfying (5.2). Moreover, assume

(5.3) u 0 ∈ L 1 (R d ) ∩ L ∞ (R d ) if α ∈ (0, 1], L 1 (R d ) ∩ ∩ p>pα H α-1,p (R d ) if α ∈ (1, 2), with p α = d α-1 > 1.
There exists a unique function u in the space

(5.4) u ∈ C [0, ∞), L 1 (R d ) ∩ L ∞ (R d ) if α ∈ (0, 1], ∩ p>pα C [0, ∞), L 1 (R d ) ∩ H α-1,p (R d ) if α ∈ (1, 2),
satisfying problem (5.1) in the usual weak sense

(5.5) u∂ t ϕ -|u|∇ α-1 (G(u)) • ∇ϕ -δ∇u • ∇ϕ dt dx = 0 for all ϕ ∈ C ∞ c (Q T ).
Moreover, u(t, x) is nonnegative if the initial condition u 0 is so, and for all t > 0 and q ∈ [1, ∞] we have

(5.6) u(t, x) dx = u 0 (x) dx and u(t) q ≤ u 0 q .
Local-in-time existence of mild solutions.

Proposition 5.2. Let p > p α = d α-1 . There exists T > 0 depending only on u 0 , and a function u in the space

(5.7) u ∈ C [0, T ), L 1 (R d ) ∩ L ∞ (R d ) if α ∈ (0, 1], C [0, T ), L 1 (R d ) ∩ H α-1,p (R d ) if α ∈ (1, 2)
such that

(5.8) u(t) = e δt∆ u 0 + t 0 ∇e δ(t-s)∆ • Ψ(u(s)) ds with Ψ(u) = |u|∇ α-1 G(u), in C[0, T ], L 1 (R d ) ∩ L ∞ (R d ))
where e t∆ denotes the heat semigroup.

Remark 5.3. We identify the heat semigroup e t∆ and its kernel (4πt) -d 2 exp -|x| 2 4t . We will use the following classical fact (5.9)

∇ β e δt∆ v p ≤ C(p, r, β, δ)t -d 2 ( 1 r -1 p )-β 2 v r with 1 ≤ r ≤ p ≤ ∞, and β ∈ [1, 2).
Now we turn to the proof of Proposition 5.2.

Proof of Proposition 5.2. We look for a solution u ∈ C([0, T ], X) as a fixed point of the map

T : u → e δt∆ u 0 + t 0 ∇e δ(t-s)∆ • Ψ(u(s)) ds,
where X is chosen as follows (5.10)

X = L 1 (R d ) ∩ L ∞ (R d ) if α ∈ (0, 1], L 1 (R d ) ∩ H α-1,p (R d ) if α ∈ (1, 2).
The associated norms are

u 1 + u Y with Y = L ∞ (R d ) and Y = H α-1,p (R d )
, respectively. We show that T has a fixed point by the Banach contraction principle as soon as T = T ( u 0 X ) > 0 is sufficiently small.

In both cases, it is enough to prove the following lemma.

Lemma 5.4. For all T ∈ (0, 1), the operator T maps C([0, T ], X) into itself. Moreover, there exist C > 0 and

γ > 0 such that for all u, v ∈ B(0, R) ⊂ C([0, T ], X), (5.11) T (u) -T (v) C([0,T ],X) ≤ C 1 (R)T γ u -v C([0,T ],X) ,
where C 1 (R) is a constant which also depends on α, d, ε, m, δ (and on p if α ∈ (1, 2)).

Indeed, once this lemma is proved, we first derive (5.12)

T (u) C([0,T ],X) ≤ u 0 X + RC 1 (R)T γ
by choosing v = 0 in (5.11) and using estimate (5.9). Now it is enough to choose R = 2 u 0 X and T > 0 such that C 1 (R)T γ ≤ 1 2 in order to ensure that T maps B(0, R) into itself, and is a contraction.

The case α ∈ (0, 1]. In order to get estimate (5.11), we first write

(5.13) T (u)(t) -T (v)(t) = t 0 ∇e δ(t-s)∆ • (Ψ(u) -Ψ(v)) (s) ds,
and the difference of Ψ's is represented as

(5.14) Ψ(u) -Ψ(v) = (|u| -|v|)∇ α-1 G(u) + |v|∇ α-1 (G(u) -G(v)).
Lemma 5.5. For every α ∈ (0, 1] and p ∈ (1, ∞) there exists a constant C(p, α) > 0 such that for all u ∈

L ∞ (R d ) ∩ L p (R d ) the following inequality (5.15) ∇ α-1 (G(u) -G(v)) q ≤ C(p, α) sup |z|≤ u ∞+ v ∞ G ′ (z) u -v p holds true with 1 q = 1 p -1-α d .
Proof. Since, for α ∈ (0, 1], we have ∇ α-1 = ∇ 0 I 1-α , where the components of ∇ 0 are the Riesz transforms (which are bounded operators on L p (R d ) for each p ∈ (1, ∞)), we obtain (5.15) from estimate (3.6) as follows

∇ α-1 (G(u) -G(v)) q ≤ I 1-α (G(u) -G(v)) q ≤ C(p, α) G(u) -G(v) p ≤ C u -v p .
Now, we come back to the proof of (5.11) with

X = L 1 (R d ) ∩ L ∞ (R d ),
First, for all u, v ∈ B(0, R) ⊂ X and some q ∈ (1, ∞), ,d,q, RG ′ (2R)) and 1 q + 1 q * = 1 and 1 q = 1 p -1-α d . We used estimate (5.15) twice to get the second line in (5.16), and the inequality u r ≤ u X which is valid for all r ∈ [1, ∞] to obtain the last one.

Ψ(u) -Ψ(v) 1 ≤ C u -v q * ∇ α-1 G(u) q + C v q * ∇ α-1 (G(u) -G(v)) q ≤ C u -v q * G ′ ( u ∞ ) u p + C v q * G ′ ( u ∞ + v ∞ ) u -v p (5.16) ≤ C 0 (R) u -v X with C 0 (R) = C(α
The estimate of the second norm in X is obtained similarly: for all u, v ∈ B(0, R) ⊂ X and some q ∈ (1, ∞),

(5.17) Ψ(u) -Ψ(v) q ≤ u -v ∞ ∇ α-1 G(u) q + v ∞ ∇ α-1 (G(u) -G(v)) q ≤ C 0 (R) u -v X .
Now, we apply inequalities (5.9) with β = 1, p, r = 1, 1 and p, r = ∞, q , respectively, and the estimates (5.16), (5.17) yield

T (u) -T (v) C([0,T ],L 1 (R d )) ≤ CC 0 (R)T 1 2 u -v C[0,T ],X) , (5.18) T (u) -T (v) C([0,T ],L ∞ (R d )) ≤ CC 0 (R)T 1 2 -d 2 1 q * u -v C[0,T ],X) . (5.19)
Combining (5.13), (5.14), (5.18) and (5.19), we thus get (5.11)

for α ∈ (0, 1], with γ = 1 2 -d 2 1
q * , now with a new constant C 1 (R) = C(α, d, q * , RG ′ (2R)), where we have chosen q * > d to ensure d 2 1 q * < 1 2 . As far as the continuity of T (u) with respect to time is concerned, it is enough to study

S(u) = t 0 ∇e δ(t-s)∆ • Ψ(u(s)) ds.
We fix t ∈ [0, T ] and write for h small enough (and positive if t = 0, negative if t = T ),

S(u(t + h)) -S(u(t)) = t+h t ∇e δ(t+h-s)∆ • Ψ(u(s)) ds + t 0 ∇e δ(t-s)∆ • (e δh∆ Ψ(u(s)) -Ψ(u(s))) ds.
As above, use two key estimates (5.16) and (5.16) together with (5.9) (and the dominated convergence theorem) to conclude the proof of Proposition 5.2 for α ∈ (0, 1].

The case α ∈ (1, 2). We argue as before, using (5.13) and (5.14). We need now to state and to prove the corresponding key technical lemma. Lemma 5.6. For α ∈ (1, 2), p > d α-1 , and u ∈ H α-1,p (R d ),

(5.20)

G(u) -G(v) H α-1,p ≤ C 2 ( u H α-1,p + v H α-1,p ) u -v H α-1,p , where C 2 ( u H α-1,p + v H α-1,p ) depends on u H α-1,p , v H α-1,p and on the W 1,∞ -norm of G ′ in the interval [0, u ∞ + v ∞ ].
Proof. We use the classical identity

(5.21) G(u) -G(v) = K(u, v)(u -v) with K(u, v) = 1 0 G ′ (τ u + (1 -τ )v) dτ,
and the Moser estimate for the product of two functions in H α-1,p (R d ), see, e.g., [25, Ch. 2, ineq. (0.22)], to obtain the following inequality

(5.22) G(u) -G(v) H α-1,p ≤ u -v ∞ K(u, v) H α-1,p + u -v H α-1,p K(u, v) ∞ .
Moreover, we recall [25, Ch. 2, Prop. 4.1] that for every increasing locally Lipschitz function H we have

H(u) H α-1,p ≤ C|H ′ |( u ∞ ) u H α-1,p
for every p ∈ (1, ∞) and α -1 ∈ (0, 1). Choosing H = G ′ , we deduce that

K(u, v) H α-1,p ≤ 1 0 G ′ (τ u + (1 -τ )v) H α-1,p dτ ≤ C|G ′′ |( u ∞ + v ∞ )( u H α-1,p + v H α-1,p ). (5.23)
Moreover, we have the trivial estimate (5.24)

K(u, v) ∞ ≤ |G ′ |( u ∞ + v ∞ ).
Combining (5.22), (5.23), (5.24) and (3.5), we finally complete the proof of inequality (5.20). Now, we are in a position to obtain the estimate of the first component of the norm of X. From (5.14), we get

(5.25) Ψ(u) -Ψ(v) 1 ≤ C 2 (R)R u -v p ′ + C 2 (2R) v p ′ u -v X ≤ C 3 (R) u -v X with C 3 (R) = CR(C 2 (R) + C 2 (2R)
). We used inequality (5.20) twice, as well as the fact that u r ≤ C u X for all r ∈ [1, ∞] and u ∈ X. From (5.25) and (5.9) with β = 1 and p, r = 1, 1 , we get inequality (5.18) where

C 0 (R) is replaced with C0 (R) = C 3 (R).
The estimate of the H α-1,p -norm is obtained analogously. First, we have

(5.26) Ψ(u) -Ψ(v) p ≤ u -v ∞ ∇ α-1 G(u) p + v ∞ ∇ α-1 (G(u) -G(v)) p ≤ C 3 (R) u -v X .
From inequalities (5.26) and (5.9) with β = 1 and β = α and p, r = p, p , we get (5.27)

T (u) -T (v) C([0,T ],H α-1,p (R d )) ≤ CC 3 (R)(T 1 2 + T 1 2 -α 2 ) u -v C[0,T ],X) .
Finally, combining (5.18) and (5.27), we complete the proof of (5.11) with γ = 1 2 -α 2 and with some C1 (R). The time continuity of T (u) is proved as in the case α ∈ (0, 1], and this achieves the proof of Proposition 5.2.

Regularity of the solutions.

Corollary 5.7 (Regularity of the solutions). Consider u

0 ∈ L 1 (R d ) ∩ L ∞ (R d ) if α ∈ (0,
1] (as before), and

u 0 ∈ L 1 (R d ) ∩ ∩ p>pα H α-1,p (R d ) if α ∈ (1,
2) (a strengthened assumption). Then the solution constructed in Proposition 5.2 enjoys the following regularity

(5.28) u ∈ C 1 (0, T ), L p (R d ) ∩ C (0, T ), H 1,p (R d ) for every p ∈ (p α , ∞) with pα = d d-(1-α) if α ∈ (0, 1], d α-1 if α ∈ (1, 2).
In particular, u = u(t, x) is a weak solution of equation (5.1), i.e.

∂ t u = δ∆u + ∇ • Ψ(u) with Ψ(u) = |u|∇ α-1 G(u)
in (0, T ) × R d in the sense of distributions (cf. equation (5.5)), and

(5.29) u(t, x) dx = u 0 (x) dx for all t ∈ (0, T ).

Proof. If α ∈ (0, 1], for every u ∈ C([0, T ], X), where the space X is defined in (5.10), we obviously have

u ∈ L ∞ ((0, T ), L ∞ (R d ))
. We derive from inequality (3.6) that for all p ∈ (p α , ∞), we have

∇ α-1 G(u) ∈ L ∞ ((0, T ), L p (R d )) with pα = d d-(1-α) > 1. This implies ∇ • Ψ(u) ∈ L q (0, T ), H -1,p (R d )
for all q ∈ (1, ∞) and p ∈ (p α , ∞). Thus, the maximal regularity of mild solutions for the nonhomogeneous heat equation [START_REF] Ladyzhenskaya | tseva, Linear and quasi-linear equations of parabolic type[END_REF] gives us ∇u ∈ L q ((0, T ), L p (R d )) for every q ∈ (1, ∞) and p ∈ (p α , ∞). Consequently,

∇ • |u|∇ α-1 G(u) = f 1 + f 2 with
(5.30)

f 1 = ∇|u| • ∇ α-1 G(u), f 2 = -|u|(-∆) α 2 G(u).
First, we remark that f 1 ∈ L q ((0, T ), L p (R d )) for every q ∈ (1, ∞) and p ∈ (p α , ∞). Second, we notice that

G(u) ∈ L q ((0, T ), H 1,p (R d )), hence (-∆) α 2 G(u) ∈ L q ((0, T ), H 1-α,p (R d )) ⊂ L q ((0, T ), L p (R d ))
. Hence, we also have f 2 ∈ L q ((0, T ), L p (R d )) for every q ∈ (1, ∞) and p ∈ (p α , ∞). Using again the maximal regularity result, we obtain that

∂ t u ∈ L q ((0, T ), L p (R d ))
for every q ∈ (1, ∞) and p ∈ (p α , ∞). Thus, using the following representation in L p (R d ) for all 0 < s < t

u(t) -u(s) = t s ∂ t u(τ ) dτ,
and the Hölder inequality we obtain for every p ∈ (p α , ∞)

u(t) -u(s) p ≤ t s ∂ t u(τ ) p dτ ≤ t s ∂ t u(τ ) q p dτ 1 q (t -s) 1 q * ,
i.e. u ∈ C 0,β ((0, T ), L p (R d )) for all p ∈ (p α , ∞) and β ∈ (0, 1). Using inequality (5.15), this estimate implies that Ψ(u) ∈ C 0,β ((0, T ), L p (R d )) for every p ∈ (p α , ∞). Now, the classical theory of linear parabolic equations, see, e.g., [START_REF] Pazy | Semigroups of linear operators and applications to partial differential equations[END_REF]Ch. 4,Theorem 3.5], implies that ∂ t u ∈ C 0,β ((0, T ), L p (R d )) and u ∈ C 0,β ((0, T ), H 1,p (R d )) which is the desired regularity result.

If α ∈ (1, 2), we have u ∈ L q ((0, T ), H α-1,p (R d )) for all q ∈ (1, ∞) and p ∈ (p α , ∞), where p α = pα = d/(α -1). Equations (3.5) and (5.20) imply that (5.31)

∇ α-1 G(u) ∈ L q ((0, T ), L p (R d ))
for all p ∈ (p α , ∞) and q ∈ (1, ∞). In particular, Ψ(u

) = |u|∇ α-1 G(u) ∈ L q ((0, T ), L p (R d ))
. Hence, the maximal regularity gives us

(5.32) u ∈ L q ((0, T ), H 1,p (R d ))
for all p ∈ (p α , ∞) and q ∈ (1, ∞). We now write once again

∇ • Ψ(u) = f 1 + f 2 with f 1 = sgn u ∇u • ∇ α-1 G(u), f 2 = |u|∇ α-1 (G ′ (u)∇u).
In view of (5.31) and (5.32), we have

f 1 ∈ L q ((0, T ), L p (R d )).
We now claim that

f 2 ∈ L q ((0, T ), H -1+(2-α),p (R d )).
Indeed,

G ′ (u)∇u ∈ L q ((0, T ), L p (R d )).
This implies that

∇ α-1 (G ′ (u)∇u) ∈ L q ((0, T ), H -1+(2-α),p (R d ))
which, in turn, implies the claim. Hence,

∇ • Ψ(u) ∈ L q ((0, T ), H -1+(2-α),p (R d ))
for all p ∈ (p α , ∞) and q ∈ (1, ∞). Then, the maximal regularity implies that

u ∈ L q ((0, T ), H 1+(2-α),p (R d ))
for all for all p ∈ (p α , ∞) and q ∈ (1, ∞), thus we see that the space regularity of u is improved. More generally, the same argument shows that if

u ∈ L q ((0, T ), H β,p (R d )), with β ≤ α, then u ∈ L q ((0, T ), H β+(2-α),p (R d )).
Now choose the least integer k ≥ 1 such that β k = 1 + k(2 -α) > α, and notice that β k < 2. Then

f 2 ∈ L q ((0, T ), H β k -α,p (R d )) ⊂ L q ((0, T ), L p (R d )).
Then the maximal regularity implies that ∂ t u ∈ L q ((0, T ), L p (R d )), which implies, as was in the case α ∈ (0, 1],

that u ∈ C 0,β ((0, T ), L p (R d )) for all p ∈ (1, ∞) and β ∈ (0, 1). The previous reasoning in spaces of the form L q ((0, T ), Y ) extends readily to spaces of the form C 0,β ((0, T ), Y ). This yields the desired regularity result in the case α ∈ (1, 2).

Moreover, it is known by [20, Ch. 4, Theorem 3.2] that mild solutions of the equation

∂ t u = δ∆u + ∇ • f with f ∈ C 0,β ((0, T ), L p (R d ))
are in fact weak solutions, i.e. they satisfy the equation in the sense of distributions.

Under these regularity properties, the proof of the mass conservation property (5.29) is completely standard.

Convexity inequalities. First, we show a simple but useful technical result involving monotone functions and the fractional Laplacian.

Lemma 5.8. Let α ∈ (0, 2]. Assume that g, h ∈ C 1 [0, ∞) are strictly increasing functions. Then, for every

nonnegative v ∈ C ∞ c (R d ) we have h(v)(-∆) α 2 g(v) dx ≥ 0.
Proof. Notice that for α = 2 this lemma is obviously true, which one checks integrating by parts. Now, let α ∈ (0, 2). Since (-∆) α 2 C = 0 for every constant C ∈ R, we can assume that g(0) = 0. In the same way, we can assume that h(0) = 0, because (-∆)

α 2 w dx = 0 for every w ∈ C ∞ c (R d ). Defining w = g(v), it suffices to show that h(g -1 (w))(-∆) α 2 w dx ≥ 0 for all w ∈ C ∞ c (R d
) such that w ≥ 0. To do it, notice that f ∈ C 2 [0, ∞) defined via the relation f (s) = s 0 h(g -1 (τ )) dτ for s ≥ 0 is convex (it suffices to check that f ′′ (s) ≥ 0). Hence, using the pointwise inequality

(5.33) (-∆) α 2 g(v) ≤ g ′ (v)(-∆) α 2 v,
(see, e.g., [START_REF] Córdoba | A maximum principle applied to quasi-geostrophic equations[END_REF], [12, Lemma 1]) we obtain h(g -1 (w))(-∆)

α 2 w dx ≥ (-∆) α 2 f (w) dx = 0.
Next, we formulate a crucial technical tool used in the derivation of various integral estimates for solutions of the regularized problem (5.1).

Proposition 5.9 (Convexity inequalities). Consider a C 2 function ϕ : R → R + such that, for all r ∈ R, r = 0, ϕ ′′ (r) > 0, and

(5.34) ϕ(r) + |ϕ ′ (r)| + ϕ ′′ (r) ≤ C(|r| M1 + |r| M2 ) for some constant C > 0 and M 1 , M 2 ∈ [1, ∞).
Then for all 0 ≤ s < t ≤ T , the function u given by Proposition 5.2 satisfies

(5.35) ϕ(u(t, x)) dx+ t s ψ(u(τ, x))(-∆) α 2 G(u(τ, x)) dx dτ +δ t s ϕ ′′ (u)|∇u| 2 dx dτ ≤ ϕ(u(s, x)) dx,
where ψ(r) = |r|ϕ ′ (r) -ϕ(r) sgn r.

The proof of Proposition 5.9 is more or less classical, and we recall it in Appendix for the sake of completeness.

Remark 5.10. Remark that ψ ′ (r) = |r|ϕ ′′ (r) > 0, hence the function ψ is increasing. Since G is also increasing, the result stated in Lemma 5.8 can be applied to show that the first dissipation term

t s ψ(u(τ, x))(-∆) α 2 G(u(τ, x)) dx dτ
is nonnegative. The fact that this quantity is finite is a part of the result stated in Proposition 5.9. Moreover,

[18, Theorem 2.2] implies that g ϕ (G(u)) ∈ L 2 ((0, T ), H α 2 ,2 (R d ))
for a function g ϕ constructed from ϕ, see [START_REF] Liskevich | Some problems on Markov semigroups[END_REF] for the detailed presentation. The special case ϕ(r) = |r| p is treated below.

Remark 5.11. The convexity of ϕ also implies that the second dissipative term in (5.35) is nonnegative. Hence, Proposition 5.9 implies that ϕ(u(t, x)) dx decreases along the flow of the regularized equation (5.1).

Corollary 5.12 (Estimates of the L p -norms). For all p ∈ (1, ∞) and 0 < s < t,

(5.36) |u(t)| p dx + (p -1) t s |u| p-1 u(-∆) α 2 G(u) dx dτ + δp(p -1) t s |u| p-2 |∇u| 2 dτ ≤ |u(s)| p dx.
In particular, for p ≥ pα (see Corollary 5.7),

d dt |u| p dx ≤ -(p -1) |u| p-1 u(-∆) α 2 G(u) dx -δp(p -1) |u| p-2 |∇u| 2 . (5.37)
Thus, for all p ∈ [1, ∞], the norm u(t) p decreases as t increases.

Proof. If p ≥ 3, we can apply Proposition 5.9 with 0 ≤ s < t ≤ T , ϕ(r) = |r| p ; indeed, in this case, ϕ is a C 2 -function and satisfies the growth assumption (5.34) with M 1 , M 2 = p, p -2 . Next, since u ∈ C 1 ((0, T ), L p (R d )), we obtain (5.37) from the inequality in Proposition 5.9 by a direct computation. We leave the details to the reader.

If p ∈ (1, 2), we consider, for each η > 0, the function ϕ η such that ϕ η (0) = ϕ ′ η (0) = 0 and

ϕ ′′ η (r) = p(p -1)((r 2 + η 2 ) p 2 -1 -η p-2 ).
In particular, the function ϕ η satisfies the assumptions of Proposition 5.9, hence, we have

ϕ η (u(t, x)) dx + t s ψ η (u(s, x))(-∆) α 2 G(u(s, x)) dx ds + δ t s ϕ ′′ η (u(s, x))|∇u(s, x)| 2 dx ds ≤ ϕ η (u(s, x)) dx, with ψ ′ η (r) = |r|ϕ ′′ η (r), ψ η (0) = 0.
Letting now η → 0 and using the Fatou lemma yields the integral formulation of inequality (5.37). Thus, we just proved that u(t) p ≤ u 0 p for all t > 0 and p ∈ (1, ∞). By computing the limits as p → ∞ and p → 1, the bounds u(t) 1 ≤ u 0 1 and u(t) ∞ ≤ u 0 ∞ are also obtained.

We can now complete the proof of Theorem 5.1.

Proof of Theorem 5.1. In view of Proposition 5.2 and Corollary 5.7, it remains to prove that solutions are nonnegative if initial data are so, and that solutions are global in time.

The positivity property is derived immediately in a usual way from the conservation of mass property (5.6) and the monotonicity of the L 1 -norm. Indeed, (5.6) yields

u(T, x) dx = u + (T, x) dx -u -(T, x) dx = ((u 0 ) + -(u 0 ) -)(x) dx and |u(T, x)| dx = u + (T, x) dx + u -(T, x) dx ≤ ((u 0 ) + + (u 0 ) -)(x) dx
(here, as usual, u + = max{0, u} and u -= max{0, -u}). These inequalities imply u -(T, x) dx ≤ (u 0 ) -(x) dx, and, in particular, the assumption (u 0 ) -= 0 gives us u ≥ 0 a.e.

As far as the global existence of solutions is concerned, we argue as follows.

For α ∈ (0, 1], the time interval, where a solution is constructed via the Banach fixed point theorem, depends only on u 0 1 + u 0 ∞ , and this norm of the solution does not increase. Hence, we extend u = u(t, •) to the whole half-line [0, ∞), step-by-step.

For α ∈ (1, 2), the Duhamel formula (5.8) and inequality (5.20)

with v = 0 yield u(t) H α-1,p ≤ u 0 H α-1,p + C 1 ( u 0 p * ) t 0 (t -s) -α 2 u(s) H α-1,p ds.
Due to the singular Gronwall lemma, see, e.g., [20, Ch. 5, Lemma 6.7], we deduce that the norm u(t) H α-1,p cannot explode in finite time. This shows that local-in-time solutions of the regularized equation (5.1) can be also continued to global-in-time ones.

Hypercontractivity and compactness estimates

Hypercontractivity estimates. We now turn to prove certain L 1 → L p estimates for solutions of problem Proof. We first remark that it is enough to prove the decay estimate (6.1) for large p's, since the general result follows by the interpolation of the L p -norms combined with the estimate u(t) 1 ≤ u 0 1 from Corollary 5.12. This is the reason why we will now prove (6.1) for p ≥ max{m -1, 1, pα } = p m (see Corollary 5.7 for a definition of pα ).

We also remark that we can assume that M = u 0 1 = 1 by rescaling the solution u in the following way. First, we consider the function u(t) = 1 M u t M m-1 which satisfies equation (5.1) with suitably rescaled parameters: δ = δ M m-1 and (if applicable) ε = ε M . Scaling back, we recover the desired inequality (6.1).

We first prove the result when Corollary 5.12 holds true with G(r) = |r| m-2 r and in the differential form (5.37). This is the case when m = 2 or m ≥ 3. In the case m ∈ (1, 3), Corollary 5.12 holds true only in the integral sense (5.36) and for a regularized function G. We will see below how to pass to the limit as the regularization parameter ε goes to 0 and get (5.36) with G(r) = |r| m-2 r. For expository reasons, we prefer to present the proof when we indeed have a differential inequality, and then to explain how to adapt it if only an integral version of it is available.

The proof on Theorem 6.1 in the cases m ≥ 3 and m = 2 is split into two steps: first, we show inequalities (6.1) with non-optimal constants C which blow up for p = ∞; then, we improve those constants by an iteration method.

Decay estimates with optimal exponents and nonoptimal constants. Our computation consists in getting the following differential inequality for p ∈ (p m , ∞). Lemma 6.2. There exists a constant K = K(p, m) > 0 independent of ε > 0 and δ > 0 such that K and K -1 are bounded as p → ∞ and

(6.2) d dt u p p ≤ -K u a p ,
with a defined in (3.9).

Proof. We get from Corollary 5.12

d dt |u| p dx ≤ -(p -1) u|u| p-1 (-∆) α 2 (u|u| m-2 ) dx (6.3) ≤ - 4p(p -1)(m -1) (p + m -1) 2 ∇ α 2 |u| p+m-1 2 2 2
, after applying the Stroock-Varopoulos inequality (3.4) with w = u|u| m-2 and q = p m-1 + 1. We use next the inequality from Lemma 3.2 combined with u(t) 1 ≤ u 0 1 = 1 to estimate the right-hand side of the above inequality. Thus, we get the differential inequality (6.2) with the constant (6.4)

K = 4(m -1)p(p -1) C N (p + m -1) 2 .
With the differential inequality (6.2) in hand, a direct computation shows that every nonnegative solution of the inequality d dt f (t) ≤ -Kf (t) a p has to satisfy the algebraic decay

f (t) ≤ K a p -1 t -1 a p -1
.

We recognize (6.1) for p > p m with the constant C p = K a p -1 -1 a-p . Here, let us notice that the constants C p = C(d, α, m, p) which are obtained at this stage of the proof blow up as p → ∞. Thus, we cannot get the L ∞ -bound directly in this way.

Lemma 6.4. Consider functions f : [0, T ] → (0, ∞) nonincreasing, and g : [0, T ] → (0, ∞) increasing, smooth, and g(0) = 0. Assume that for a.e. t ∈ (0, T ), s < t,

f (t) + K t s f (τ ) γ+1 g ′ (τ ) dτ ≤ f (s).
Then for a.e. t ∈ [0, T ], we have

f (t) ≤ (Kγg(t)) -1 γ .
The proof of this lemma is given in Appendix C for the readers' convenience. Now, the proof of Theorem 6.1 is complete.

Compactness estimates. Now, we prove estimates which will allow us to pass to the limit as ε → 0 and δ → 0, successively, in the regularized problem (5.1). We are going to use the fact that the approximating functions G = G ε are γ-Hölder continuous with γ = min{m -1, 1} on the interval 0, u 0 ∞ uniformly in ε ∈ [0, 1]. These estimates are formulated in the following technical lemmas.

Lemma 6.5. Assume that G is γ-Hölder continuous with γ = min{m -1, 1}. For every α ∈ (0, 1] and

m > 1 + 1-α d ∈ (1, 2), we have (6.8) ∇ α-1 G(u) q ≤ C u γ γp with q ≥ (m -1 -(1 -α)/d) -1 and p = (1/q + (1 -α)/d) -1 .
Proof. Here, it suffices to combine the estimate on the Riesz potential with G(u) p ≤ C u γ γp . The choice of q (or equivalently, the restriction on m) ensures that γp ≥ 1.

The compactness of a sequence of solutions to the regularized problem (5.1) in the case α ∈ (1, 2) is a consequence of the following estimate. Lemma 6.6. Consider a γ-Hölder continuous function G : R → R. Then, for every α ∈ (1, 2) such that

(2 -γ)α < 2, we have ∇ α-1 G(u) p ≤ C u H α 2 ,2 with α -1 γ - d p = α 2 - d 2 .
Remark 6.7. Note that the assumptions of Lemma 6.6 ensure that p > 2.

Proof of Lemma 6.6. We use successively the fact that H

∇ α-1 G(u) p ≤ C G(u) F α-1 p,2 ≤ C[G] γ u F α-1 γ p,2γ ≤ C u H α 2 ,2 ,
where [G] γ := sup r =s |G(r) -G(s)|/|r -s| γ and p > 2 is chosen so that the so-called differential dimension is constant; this yields the condition appearing in the statement of the lemma. The proof is now complete. Lemma 6.8. For α ∈ (1, 2) and (3 -m)α < 2, there exists p ∈ (1, ∞) such that (6.9)

∇ α-1 G(u) p ≤ C u H 1,2 (R d )
where C depends on the γ-Hölder seminorm of G in (0, u ∞ ) (with γ = min{m -1, 1}).

Proof. Apply Lemma 6.6 with γ = m -1.

Lemma 6.9. If α ∈ (1, 2) and m > α, there exists p ∈ (1, ∞) and r ∈ (1, ∞) (and r > m -1) such that (6.10)

∇ α-1 |u| m-1 sgn u p ≤ C |u| r+m-1 2 sgn u H α 2 ,2 (R d ) . Proof. Consider v = |u| r+m-1 2
sgn u. Then estimate (6.10) is equivalent to the following one

∇ α-1 |v| 2(m-1) r+m-1 sgn v p ≤ C v H α 2 ,2 (R d ) .
We then apply Lemma 6.6 with γ = 2(m-1) r+m-1 and get the desired result if γ > 2(1 -1 α ), or equivalently,

r < m -1 α -1 .
Here, in order to find r ∈ 1, m-1 α-1 , the condition m > α is needed. The proof is now complete.

Proof of Theorem 2.6

This section is devoted to the proof of our main result on the existence of solutions to problem (1.1)-(1.2) satisfying decay estimates (2.6).

Proof of Theorem 2.6. We first consider very regular initial data, i.e. we assume that u 0 satisfies (5.3). Then, this condition is relaxed by considering initial data that are merely integrable.

The proof proceeds in three steps: passage to the limit with the parameter of the regularization of the nonlinearity, then with the parameter of the parabolic regularization, and finally -stability with respect to initial data.

Passage to the limit as ε → 0. Consider u 0 satisfying (5.3). From Theorem 5.1, we have a sequence of solutions u ε of (5.1) for G ε defined in such a way that there exists C > 0 such that for all ε ∈ (0, 1],

[G ε ] γ := sup r =s |G ε (r) -G ε (s)| |r -s| γ ≤ C
where γ = min{m -1, 1}. Thanks to Corollary 5.12, there exists also a constant C > 0 (depending on δ > 0)

such that for all ε ∈ (0, 1],

(7.1) u ε L ∞ ((0,T )×R d ) ≤ C, u ε L 2 ((0,T ),H 1,2 (R d )) ≤ C.
Hence, we can construct a sequence (ε n ) n such that

u n ⇀ u in L 2 ((0, T ), H 1,2 (R d ))
where u n denotes u εn . Moreover, for all R > 0, the embedding H Hence, we infer from [START_REF] Rakotoson | An optimal compactness theorem and application to elliptic-parabolic systems[END_REF] (which contains an optimal result on the compactness for Hilbert space valued vector functions) that, up to a subsequence, for every R > 0,

u n → u in L 2 ((0, T ) × B R ).
Moreover, passing again to a subsequence if necessary, we can assume that for every R > 0, (7.2)

∇u n ⇀ ∇u in L 2 ((0, T ) × B R ), u n → u for a.e. (t, x) ∈ Q T , which imply for all ϕ ∈ C ∞ c ( QT ), QT u n ∂ t ϕ dt dx → QT u∂ t ϕ dt dx, QT ∇u n • ∇ϕ dt dx → QT ∇u • ∇ϕ dt dx.
As far as the nonlinear term in equation (5.1) is concerned, we have

|u n | → |u| for a.e. (t, x) ∈ Q T .
Hence, the dominated convergence theorem combined with (7.1) imply that

|u n |∇ϕ → |u|∇ϕ in L 2 ((0, T ), L q ′ (R d ))
for all q ′ ∈ (1, ∞). 6.5 and 6.8 imply that

(7.3) ∇ α-1 G εn (u n ) is bounded in L q (Q T )
for some q ∈ (1, ∞).

We deduce from (7.3) that we can extract a weakly converging subsequence of (

∇ α-1 G εn (u n )) n ⊂ L q (Q T ) which limit is denoted by L. Since G εn (u n ) → |u| m-1 sgn u a.e. in Q T , we have ∇ α-1 G εn (u n ) ⇀ L in L q (Q T ), where L = ∇ α-1 (|u| m-1 sgn u) in D ′ (Q T ).
In particular, ∇ α-1 (|u| m-1 sgn u) ∈ L q (Q T ) holds, and

QT ∇ α-1 G εn (u n ) • |u n |∇ϕ dt dx → QT ∇ α-1 (|u| m-1 sgn u) • |u|∇ϕ dt dx
as n → ∞. Hence, we obtain a weak solution of (5.1) with G(r) = |r| m-1 sgn r for u 0 satisfying (5.3).

By the Fatou lemma, we derive from Corollary 5.12 and Proposition 3.1 the following estimates: for a.e.

t > s > 0, q ∈ [1, ∞], p ∈ (1, ∞), uniformly in δ ∈ (0, 1], u(t) q ≤ u 0 q , (7.4) u(t) Remark that now we can derive the hypercontractivity estimates in the case m ∈ (1, 3) from inequality (7.5).

We also know that the solution u is nonnegative if u 0 is so.

Passage to the limit as δ → 0. Let u = u δ denote the solution constructed above. We consider v δ = |u| m sgn u. Using (7.5) with s = 0 and p = m + 1 > 1, we get that, for all T > 0,

v δ is bounded in L 2 ((0, T ), H α 2 ,2 (R d )).
Hence, there exists a subsequence δ n → 0 such that v n ⇀ v in L 2 ((0, T ), H Hence, we can use [START_REF] Rakotoson | An optimal compactness theorem and application to elliptic-parabolic systems[END_REF] once more, and conclude that for all R > 0

v n → v in L 2 ((0, T ) × B R ).
Moreover, passing to a subsequence if necessary, we can assume that u n → u for a.e. (t, x) ∈ Q T , where u n denotes u δn . Now Lemmas 6.5 and 6.9 imply that (7.8)

∇ α-1 (|u n | m-1 sgn u n ) is bounded in L q (Q T )
for some q > 1. Hence, we can pass to the limit in the nonlinear term in the weak formulation of equation (5.1).

To complete this proof, notice that inequality (7.7) implies that

δ n ∇u n → 0 in L 2 (Q T ).
In particular,

δ n QT ∇u n • ∇ϕ dt dx → 0
as n → ∞. We thus conclude that u is a weak solution of (1.1)-(1.2).

with 1 q = 1-2 p . As far as Θ η is concerned, we choose it such that Θ η (s) = 0 and Θ ′ η (τ ) = ρ η (τ -s-η)-ρ η (τ -t+η)

where ρ η is an even mollifier supported in Hence, l n does not blow up, and neither does κ n .

Appendix C. Proof of Lemma 6.4

Proof of Lemma 6.4. First, we remark that we can reduce to the case g(t) = t and K = 1 through a change of variables.

The proof is simple if f is smooth. If f is not, extend f by 0 to R and consider a mollifier ρ ε . Then write for t 1 < t 2 ,

f (t 1 -s) + K t1-s -∞ f γ+1 (τ ) dτ ≤ f (t 2 -s).
Now integrate against ρ ε (s) and use the Jensen inequality to get

f ε (t 1 ) + K t1 -∞ f γ+1 ε (τ ) dτ ≤ f ε (t 2 ).
We are now reduced to the case f = f ε , which is smooth. Passing to the limit, the proof is now complete.

α 2 - 1 ,

 21 α ∈ (0, 2). The equation is supplemented with an initial condition (1.2) u(0, x) = u 0 (x). First, we construct nonnegative self-similar solutions of equation (1.1) which are explicit and compactly supported. They generalize the classical Barenblatt-Kompaneets-Pattle-Zel'dovich solutions of the porous medium equation, see (1.4) below. Second, we prove the existence of sign changing weak solutions to problem (1.1)-(1.2)

  x) satisfies the equation in the pointwise sense for |x| = t 1 d(m-1)+α , and is min α 2(m-1) , 1 -Hölder continuous at the interface |x| = t 1 d(m-1)+α . Remark 2.3. When α = 2 in expression (2.4) below, we recover the classical Barenblatt-Kompaneets-Pattle-

  and with a constant C N = C(d, α) > 0. The proof of (3.7) for d = 1 can be found in, e.g.,[START_REF] Karch | On convergence of solutions of fractal Burgers equation toward rarefaction waves[END_REF] Lemma 2.2], and this extends easily to the general case d ≥ 1.

2 +

 2 which means that d = 1 and α ≤ 1. The critical case α = 1 can be obtained by passing to the limit as α ց 1; indeed, the constants C α,2-α,d (d + α -2) and Cα,2-α,d (d + α -2) appearing in (4.3) simplify thanks to the relation zΓ(z) = Γ(z + 1). If now α < 1, we can argue as above by analytic continuation. The proof is now complete. Now we turn to the proof of the main technical lemma. Proof of Lemma 4.1. We first assume that β ∈ (0, d) and that γ > max{0, d -2β -1}, and we then argue by the analytic continuation with a choice of parameters corresponding each time to 2 F 1 defined and bounded for all |y| ≤ 1. The Fourier transform of Φ γ (y) = 1 -|y| 2 γ is expressed in terms of Bessel functions, see, e.g., [23, Ch. IV, Sec. 3]

( 5 .

 5 1). Theorem 6.1 (L p -decay of solutions to the regularized problem). Let u = u(t, x) be a solution to the regularized problem (5.1) constructed in Theorem 5.1. There exists a constant C = C(d, α, m) > 0 such that for all ε > 0, δ > 0 and p ∈ [1, ∞],

α 2 , 2 (

 22 R d )),where v n denotes v δn . Moreover, inequality (7.4) with q = ∞ implies that for R > 0 lim meas(E)→0,E⊂[0,T ] E BR |v n (t, x)| 2 dt dx = 0.

(

  [-η, η]. Now we can write(∂ t u)ϕ ′ (u)Θ η dτ dx = ∂ t (ϕ(u))Θ η dτ dx = ϕ(u)(ρ η (τ -t + η) -ρ η (τ -s -η)) dτ dx.Hence,(∂ t u)ϕ ′ (u)Θ η dτ dx → ϕ(u(t, x)) dx -ϕ(u(s, x)) dx.Moreover,|u|∇φ(u) = |u|ϕ ′′ (u)∇u = ∇(ψ(u)) in L p (R d )thanks to the Stampacchia theorem (u ∈ L ∞ (R d ) hence ϕ ′ and ψ locally Lipschitz is enough). Hencet s |u|∇ α-1 G(u) • ∇φ(u) dτ dx = t s ∇ α-1 G(u) • ∇(ψ(u)) dx dτ. It remains to prove that 1 G(u) • ∇(ψ(u)) dx dτ < ∞.The last inequality comes from the computations we made above. As far as the second inequality is concerned, we use (1.3) to writet s ∇ α-1 G(u) • ∇(ψ(u)) dx dτ = lim η→0 t s (G(u)(τ, x) -G(u)(τ, y))F η (y -x) • ∇(ψ(u(τ, x))) dx dy dτwhere F η is defined as followsF η (z) = C d,α z η d+α + |z| d+α .Through an integration by parts, we now gett s (G(u)(τ, x) -G(u)(τ, y))F η (y -x) • ∇(ψ(u(τ, x))) dx dy dτ = -t s ∇G(u)(τ, x) • F η (y -x)ψ(u(τ, x)) dx dy dτ + t s (u(τ, x) -u(τ, y))(∇ • F η )(y -x)ψ(u(τ, x)) dx dy dτ.The first term on the right-hand side equals 0 since F η is odd. Moreover, we have∇ • F η (z) = C d,α η d+α + (d + α + 1)|z| d+α (η d+α + |z| d+α ) 2 .Hence, the Fatou lemma yields (A.1). The proof of the proposition is now complete.Appendix B. From (6.6) to the boundedness of κ n Consider l n = log κ n and write (6.6) as follows l n+1 ≤ a n + b n l n with

  there exists a global-in-time weak solution u of the Cauchy problem (1.1)-(1.2). Moreover,

			u(t, x) dx = u 0 (x) dx
	and
	(2.6)	u(t) p ≤ C(d, α, m) u 0	d(m-1)/p+α d(m-1)+α 1	t -	d d(m-1)+α 1-1 p	for all t > 0,
	holds with the constant C(d, α, m) independent of p and u 0 .
	The solution u is nonnegative if the initial condition u 0 is so. If u 0 ∈ L p (R d ) for some p ∈ [1, ∞], then
			u(t) p ≤ u 0 p
	holds for all t > 0.
	Remark 2.7. Estimates (2.6) are sharp since the decay in Theorem 2.6 corresponds exactly to that for self-similar
	solutions constructed in Theorem 2.2. Moreover, for α = 2, they are similar to those for degenerate partial
	differential equations like the porous medium equation (showing the regularization effect on the L p -norms of
	solutions); see, e.g., [28], [8, Ch. 2].

  1,2 (B R ) ⊂ L 2 (B R) is dense and compact.

	Using (7.1) we also obtain
	lim meas(E)→0,E⊂[0,T ] E BR	|u n (t, x)| 2 dt dx = 0.

The authors wish to thank Jean Dolbeault and Régis Monneau for the fruitful discussions they had together. The authors were supported by an EGIDE project (PHC POLONIUM 20078TL, 0185, 2009-2010). The second author was supported by an ANR project (EVOL). The first and the third authors were supported by the MNSzW grant N201 418839, and the Foundation for Polish Science operated within the Innovative Economy Operational Programme 2007-2013 funded by European Regional Development Fund (Ph.D. Programme: Mathematical Methods in Natural Sciences).

Recurrence step. To improve the constant C p and to handle the limit case p = ∞, we apply a variation on the Moser-Alikakos method of estimating the L p -norms with p = 2 n recursively, see, e.g., [START_REF] Alikakos | An application of the invariance principle to reaction-diffusion equations[END_REF] and [START_REF] Karch | On convergence of solutions of fractal Burgers equation toward rarefaction waves[END_REF]Lemma 3.1].

The starting point is the already obtained estimate for p = 2 k ≥ m -1 with the least integer k. Lemma 6.3. For each n ≥ k, the following estimate (6.5) u(t) 2 n ≤ κ n t -µn for all t > 0, holds true with µ n = 1-2 -n α d and with a positive κ n satisfying the recursive estimate

where K n is given by formula (6.4) with p = 2 n+1 .

Thus, having this estimate we see that lim sup n→∞ κ n < ∞ (irrespective of the value of κ k at the beginning of the recurrence), essentially since ∞ n=k n2 -n < ∞, see Appendix B for details. Recall that the constants K = K n in the preliminary estimates have been such that K n and K -1 n were bounded uniformly when p = 2 n → ∞.

Proof of Lemma 6.3. We combine (6.3) with the Nash inequality (3.10) and two Hölder inequalities with

As the result, we have

with some K as in (6.2) corresponding to p = 2 n+1 . Next we estimate the L 2 n -norm of u(t) using (6.5), and arrive to the differential inequality for

.

Finally, we integrate this inequality on [0, t] and take a suitable negative power to arrive to (6.6).

Interpolation between p = 2 n and p = 2 n+1 and the passage to the limit p → ∞ finish the proof of the hypercontractivity estimates in the cases m ≥ 3 and m = 2.

To deal with the case m ∈ (1, 3), instead of the differential inequality from Corollary 5.12, we use its integral counterpart (5.36). In particular, inequalities (6.2) and (6.7) have their integral counterparts as well (see inequality (7.5), below). Hence, the previous proof works in this case by applying the following lemma with g(t) = t and g(t) = Ct ν for some well chosen positive constants ν, C, successively.

The conservation of mass, the positivity property, the monotonicity of L p -norms, and the hypercontractivity estimates follow from (7.4)-(7.7) in a standard way.

Stability with respect initial conditions. Assume now that

Consider an approximating sequence u n 0 satisfying (5.3). Then, the sequences of L p -norms are bounded and we can pass to the limit as we did already when letting δ → 0. We also recover the expected properties of the weak solution.

Assume finally that u 0 is merely integrable. Then

T n (r) = min{max{r, -n}, n}.

Hence, there exists a weak solution u n of (1.1). In view of the hypercontractivity estimates, we can extract a converging subsequence in (η, T ) × R d , with arbitrary η > 0, to a function u η in the following sense

Moreover, the L 1 -norm of u n is bounded in (0, T ) × R d . Then we have

as n → ∞ for each η > 0. We can now conclude the proof of Theorem 2.6 through a diagonal procedure.

Appendix A. Proof of Proposition 5.9

Proof of Proposition 5.9. We fix an arbitrary p ≥ p α = d α-1 . First, we remark that

for all φ compactly supported in time and in

To justify the previous equality, it is enough to mollify the function φ in time and space, and remark that the mollified function φ η satisfies (for η small enough) φ η (0, x) = φ η (T, x) = 0, uφ η ∈ C 1 ((0, T ), L 1 (R d )) and u(∂ t φ η ) dt dx = -(∂ t u)φ η dt dx.

Letting η → 0 and using the regularity of Ψ(u) yields the desired result.

Next, consider φ(τ, x) = ϕ ′ (u(τ, x))Θ η (τ ), where Θ η is truncation function in time of [s, t]. Remark that φ ∈ L 1 ((0, T ), H 1,p ′ (R d )); indeed, ϕ ′ (u) ∈ L ∞ ((0, T ), L p ′ (R d )), and we also have for all w ∈ (H


Jesús González Feliu
email: jesus.gonzales-feliu@let.ish-lyon.cnrs.fr

Frédéric Henriot

Jesús González Feliu

Estimating traffic flows and environmental effects of urban commercial supply in global city logistics decision support

We present an integrated transport and land-use model for prediction of traffic flows for the entire urban logistics system. The model shows the relation between several aspects of commercial strategies (store location, commercial area's supply, and urban planning policies) and both the upstream distribution of goods and the downstream usage of private vehicles for shopping. Several case studies based on representative scenarios are proposed.

Introduction

City logistics studies the urban part of the logistic chain and more precisely freight movements in order to reduce costs and externalities related to urban freight distribution, mainly congestion and pollution.

Two main movements can be identified for freight in urban areas: logistic movements are related to freight distribution between the different activities, i.e. the exchange of goods between establishments; end-consumer movements can be defined as the trips where the purchased goods are moved by the consumer [START_REF] Russo | Demand models for city logistics : a state of the art and a proposed integrated system[END_REF]. The first group of urban freight movements has been one of the most studied subjects in city logistics research. The second group, commonly studied in other fields [START_REF] Thill | Analyse des décisions spatiales et du processus de choix des consommateurs : théorie, méthodes et exemple d'applications[END_REF] but often not modelled as a part of the global urban freight distribution chain, is however an important component in terms of traffic volumes. In France, shopping trips represent from 10% to 25% of the total person trips, depending on the day of the week [START_REF] Routhier | Mesurer l'impact du transport de marchandises en villele modèle Freturb (version 1)[END_REF]. Moreover, end-consumer movements represent among 50% of the total km.PCUs (Passenger Car Units) of urban freight transportation. In terms of pollution, the urban transport of goods, including end-consumer movements, produces about 25% of the total CO 2 emissions for transport, 35% of the NO

x emissions and 40 to 50% of the solid particles [START_REF] Segalou | The environmental assessment of urban goods movement[END_REF].

The aim of this paper is to propose an integrated road traffic flow prediction model for urban logistics systems, then to apply it, using a set of scenarios to represent several future states, to the metropolitan area of Lyon in France.

Freight movements estimation: a global approach

The model is built considering the interactions between urban planning for commercial areas, land-use strategies, distribution solutions and shopping trip flows, in order to evaluate the impact of the different city logistics approaches on the overall system, defined by the different types of goods movements. The model is composed of two modules: logistic movements simulation (FRETURB) and end-consumer movements generation (Shopping Trip Generation, STG).

FRETURB is a prevision and diagnosis model for decision support in urban freight distribution that integrates several modules [START_REF] Routhier | Mesurer l'impact du transport de marchandises en villele modèle Freturb (version 1)[END_REF]. In our integrated approach, we will use the following two:

• Pick-ups and deliveries demand generation in each urban area. This module's variables are related to the establishments: type of activity, number of jobs and number of subsidiaries. • Road occupancy by running vehicles module. This module estimates the traffic flow of freight distribution vehicles in each area, according to the number of stops in the rounds, the vehicle type, the operating mode and the density of activities. The distance to the town centre is considered in case of direct trips and the first leg in case of rounds.

The model has been calibrated from logistic behaviours observation, for several surveys in the late 1990's . Organisation and management of the logistic chains have been analysed according to a large typology of establishments (116 types according to the nature of their activity and their size).

We introduce a new Shopping Trip Generation model built from the following definition of private car shopping trips: for each purchasing activity, two related trips are defined, one ending at the shopping activity and the other starting at it. First of all, the demand is estimated considering that the shopping destination is chosen according to the commercial zone attractiveness, which depends on the commercial supply, the type of purchase and the available modal choices. Finally, the O/D matrix is defined, analogously to FRETURB's "Road occupancy by running vehicles": the traffic flow of shopping trips by private car is estimated in each area from the trip generation data, according to the distance, to the type of urban area and to socio-economic variables.

Figure 1 A Flowchart of the integrated simulation model

The method works as follows. Consider a city, divided into representative zones. Each zone is defined by two data sets: one defines the economical activities of the zone (establishments data) and the other characterises its population (inhabitants data). From the first data set, the FRETURB module calculates the vehicle flow exchanges between zones for freight distribution and collection, and gives to the STG module the composition of each commercial area (in number of stores and employees, by type of commercial activity). The STG module calculates analogously private cars shopping trips. After that, these flows are integrated in order to calculate the total number of km.PCU for freight-related vehicle trips, to calculate the congestion level. Finally, CO 2 emissions and local pollution rates can be also obtained by considering the characteristics of the different vehicles (professionals and individuals) in the considered city.

Proposed scenarios and computational experiments

The model is applied to the urban area of Lyon (near 2 000 000 inhabitants), using the information provided by the 2006 household survey (covering 454 towns and near 830 000 households), and the 2005 SIRENE directory, the national exhaustive firm register, providing establishments NACE among other information.

Different scenarios are proposed to simulate four contrasted extreme commercial land-use policies.

Starting from the current situation (year 2006), the simulation is conducted to know the effects, at the horizon 2025, of an exclusive development of:

• big commercial areas on the periphery, with specialised retailers and hypermarkets;

• medium supermarkets and specialised stores, in every district of the city;

• e-commerce with home delivery services, that should help to limit the shopping trips, but needs specialised distribution services on the upstream side of the chain; • e-commerce with reception points, where the consumer can go by car or by other modes, depending on the quality of the network and its sufficiency in every part of the city.

A fifth simulation should be a more realistic scenario with a mixture of the commercial policies presented above, in order to find a sustainable mix of policy choices for the coming decades.

The main results show that the decisions in urban planning and the choices of the stakeholders have significant effects on the overall traffic flows of the urban goods movements. Peri-urban location of commercial areas decreases the upstream supply-chain costs (concentration of commercial activities, no need to travel into the dense city area), but increases the downstream transportation costs (raise of the usage of private car for shopping), because of both unavailability of public transport and low density of the built urban space. The third and fourth scenarios appear to be better configurations in terms of congestion and environmental impacts of combined logistic and end-consumer movements.

Figure 2 A

 2 Figure 2 A schema of the four extreme scenarios