R Aboul Hosn

L Sibille

N Benahmed

B Chareyre

A discrete numerical model involving a partially fluid-solid coupling to describe suffusion effects in soils

Keywords: Suffusion, Constitutive behavior, Discrete numerical model, Pore finite volume

Suffusion, a particular case of internal erosion in soils, is a particle-scale mechanism involving the selective erosion of fine particles within the matrix of coarse particles under the effect of water seepage. Modifications in the micro-structure as a consequence of the loss of a fine fraction may affect both the hydraulic and the mechanical properties of the soil. Thus leading, in some cases, to drastic consequences. In this study, a discrete numerical model methodology is introduced to investigate the initiation and development of suffusion as well as to analyze its effects on the soil mechanical properties. For that purpose, an original numerical extraction procedure was developed allowing us to mimic the suffusion process by taking into account both the micro-structure of the granular packing and the hydraulic loading in the suffusion development. Such a procedure is based on a one-way fluid-solid coupling where the interstitial flow is solved with a finite volume approach defined at the pore scale. Numerical soil samples subjected to different hydraulic gradients show that depending on the amount and the role of eroded particles in the granular assembly, the eroded medium either shows negligible deformations during erosion but then collapses suddenly once sheared or it deforms significantly during erosion and recovers part of its strength once subjected to shear forces. A non-linear non-monotonic

Introduction

Internal erosion is a major cause of the failure of hydraulic earthen structures [START_REF] Stewart | The wac bennett dam sinkhole incident[END_REF][START_REF] Fell | The state of the art of assessing the likelihood of internal erosion of embankment dams, water retaining structures and their foundations[END_REF][START_REF] Zhang | Seepage failure mechanism of the gouhou rockfill dam during reservoir water infiltration[END_REF]. It can be initiated under different forms: concentrated leak erosion, backward erosion, soil contact erosion, and suffusion [START_REF] Fell | The state of the art of assessing the likelihood of internal erosion of embankment dams, water retaining structures and their foundations[END_REF]. This study focuses on the suffusion process involving the selective erosion of fine particles within the matrix of coarse particles. Such a particle-scale mechanism is characterized by the detachment and the migration of fine grains by the interstitial flow leaving behind the granular skeleton. These modifications in the soil microstructure may lead to deformations at the macroscopic scale and may influence significantly the mechanical behavior of the soil. Soils susceptible to suffusion are described as internally unstable. Different factors may influence the internal stability of soils such as geometric conditions (i.e., grain sizes, pore sizes and constriction sizes) and hydraulic conditions (i.e., hydraulic gradient, flow velocity and flow direction). Several geometric methods have been defined in the literature to evaluate the likelihood of internal instability. They are classically based on the particle size distribution (see for instance [START_REF] Kenney | Internal stability of granular filters[END_REF][START_REF] Burenkova | Assessment of suffusion in noncohesive and graded soils[END_REF]). However, such criteria are usually found to be conservative [START_REF] Li | Comparison of two criteria for internal stability of granular soil[END_REF][START_REF] Wan | Assessing the potential of internal instability and suffusion in embankment dams and their foundations[END_REF].

According to [START_REF] Kenney | Controlling constriction sizes of granular filters[END_REF], the controlling void entities in the filtration processes are the constrictions, defined as the narrow sections between pores. Thus, detached particles should overcome such key obstacles to pass through the filter. Therefore, recent research works based on geometric criteria take into account the constriction size distribution, CSD, instead of using directly the particle size distribution, PSD [START_REF] Vincens | Approaches to determine the constriction size distribution for understanding filtration phenomena in granular materials[END_REF][START_REF] Silveira | An analysis of the problem of washing through in protective filters[END_REF][START_REF] Indraratna | Analytical modeling and experimental verification of granular filter behaviour[END_REF][START_REF] Locke | Time-dependent particle transport through granular filters[END_REF]. Other recent research works study the degree of interlocking of fine particles with neighboring particles with respect to their participation in the stress transfer in the granular skeleton [START_REF] Shire | Micromechanical assessment of an internal stability criterion[END_REF].

2

Even if soils satisfy a geometric criterion, sufficiently strong hydraulic loading is necessary to detach particles [START_REF] Kovacs | Seepage hydraulic[END_REF][START_REF] Terzaghi | Soil mechanics: a new chapter in engineering science[END_REF][START_REF] Skempton | Experiments on piping in sandy gravels[END_REF]. Moreover, suffusion may be accompanied by the filtration of some detached particles causing a modification in the pore spaces. This may induce a clogging process and may result in a decrease in the hydraulic conductivity [START_REF] Marot | Internal flow effects on isotropic confined sand-clay mixtures[END_REF].

During the suffusion process, the soil structure is modified and the post-suffusion behavior may be strongly affected. Intensive researches have been devoted to study the role of fine particles in the soil structure. Studies showed that the fines content may influence significantly the micro-structure of the soil. The participation of the particles of different sizes in the inter-particle contact force network changes with the fines content [START_REF] Thevanayagam | Effects of fines and confining stress on undrained shear strength of silty sands[END_REF][START_REF] Thevanayagam | Liquefaction potential and undrained fragility of silty soils[END_REF][START_REF] Shire | Fabric and effective stress distribution in internally unstable soils[END_REF][START_REF] Shire | The influence of finer fraction and size-ratio on the micro-scale properties of dense bimodal materials[END_REF]. Moreover, fines content has a significant influence on the shear strength and the onset of mechanical instabilities in soils as well as on the position of the critical state line of silty sands [START_REF] Benahmed | An experimental investigation into the effects of low plastic fines content on the behaviour of sand/silt mixtures[END_REF][START_REF] Yin | A micromechanics-based model for sandsilt mixtures[END_REF]. Therefore, the soil behavior seems to be dependent on the range of fines content, which may explain the changes in the soil strength after suffusion. Nevertheless, volumetric deformation of soil (subjected to its own weight and/or an external loading) may also be induced by suffusion development and changes in soil mechanical properties may not be deduced from considerations of the fines content only.

Besides, several laboratory investigations have been carried out to study the development of suffusion and its induced effects on the soil properties. Erosion tests, by developing suitable apparatus for that purpose, showed the influence of soil gradation and hydraulic gradient on erosion rates. Changes in the void ratio and the hydraulic conductivity were also reported. It was found in [START_REF] Sibille | A description of internal erosion by suffusion and induced settlements on cohesionless granular matter[END_REF] that the initial fines content and the history of the hydraulic loading play a major role in the suffusion development and an expression of the erosion rate was proposed to predict the cumulative eroded mass. Concerning the mechanical behavior, contradictory conclusions were found. An increase or a decrease in the soil strength may occur [START_REF] Ke | Strength reduction of cohesionless soil due to internal erosion induced by one-dimensional upward seepage flow[END_REF][START_REF] Ke | Experimental investigations on suffusion characteristics and its mechanical consequences on saturated cohesionless soil[END_REF][START_REF] Ke | Triaxial erosion test for evaluation of mechanical consequences of internal erosion[END_REF][START_REF] Xiao | Experimental investigation of the effects of suffusion on physical and geomechanic characteristics of sandy soils[END_REF][START_REF] Chang | A stress-controlled erosion apparatus for studying internal erosion in soils[END_REF][START_REF] Chang | Laboratory investigation of initiation and development of internal erosion in soils under complex stress states[END_REF]. Therefore, further investigations are necessary in order to describe and model how a suffusive erosion affects the constitutive behavior of soils. One of the objectives of this paper is thus to analyze the relation between the degree of development of a suffusion process and the changes of the soil mechanical properties from a direct description of the degradation of the micro-structure induced by the water seepage.

In addition to experimental studies, few numerical studies are available to describe the relation between the removal of soil particles and its consequences at the macroscopic scale. Some include direct coupling between solid and liquid phases [START_REF] Sari | Investigation of internal erosion process using a coupled dem-fluid method[END_REF] to analyze the migration and filtration processses during internal erosion. More recently, [START_REF] Tejada | Role of blockages in particle transport through homogeneous granular assemblies[END_REF] investigated the effects of the intermittent blockage of constrictions and its consequences on the particle transport showing an exponential decay between consecutive trapping events. Other numerical studies immitate erosion by defining a procedure allowing the removal of particles without considering a fluid phase [START_REF] Scholtès | Multiscale approaches to describe mechanical responses induced by particle removal in granular materials[END_REF][START_REF] Wood | Changing grading of soil: effect on critical states[END_REF][START_REF] Wood | Modelling mechanical consequences of erosion[END_REF]. This latter is an efficient method from a numerical point of view being simple and requiring less computational resources. However, the removal of particles is usually based on the particle's size and the stress that a particle holds. So this method lacks other criteria that play a major role in the erosion of fine particles, such as, the constriction sizes to allow particles to pass through or the driving fluid forces to cause particles to move. Hence, a second objective in this work is to develop the existing extraction criteria by taking into account complex geometric and hydraulic processes which govern the suffusion phenomenon.

With the aim of studying the mechanical behavior of soils subjected to suffusion, this paper introduces a numerical model based on the discrete element method (DEM). In the first part of the article, an extraction process involving both a detachment criterion and a transport criterion of potentially erodible particles is defined to mimic the suffusion development. This process introduces a one-way coupling between the fluid and the solid phases by associating the DEM with the pore scale finite volume method (PFV) as developed in [START_REF] Chareyre | Pore-scale modeling of viscous flow and induced forces in dense sphere packings[END_REF][START_REF] Catalano | Pore-scale modeling of fluidparticles interaction and emerging poromechanical effects[END_REF] to solve the interstitial fluid flow and to compute the fluid forces applied on each soil particle. In a second part, suffusion like simulations are presented for a granular assembly subjected to a water seepage under controlled hydraulic gradients. On one hand, the modification of the granular micro-structure and its macroscopic response during suffusion development are analyzed. On the other hand, eroded samples are sheared under triaxial compression conditions to characterize their mechanical properties and the evolution of the latter with the degree of development of the suffusion. Finally, limitations and further development of the numerical model are discussed.

Discrete numerical model

Solid granular phase

The solid phase is represented as a granular assembly with the discrete element method. Periodic boundary conditions are adopted and gravity is neglected in all simulations, which may result in having rattlers. The periodic cell is formed as a parallelepipedic block initially filled with a cloud of spherical particles. Then, an isotropic compression is applied until reaching the required confining pressure and initial density. From this state, either the erosion process (as described in section 3) is carried out, or triaxial compressions are simulated to characterize the macroscopic mechanical properties of the granular assembly (as discussed in section 5).

The inter-particle interaction is modeled by a linear elastic relationship between contact forces and relative displacements of particles, with a slip Coulomb model. Thus, it is based on a very classical interaction law characterized by constant normal and tangential (or shear) stiffnesses, K n and K s , respectively, and a contact friction angle, ϕ c , such that the normal and tangential contact forces, F n and F s , are defined by:

F n = K n δ n n (1)
∆

F s = -K s ∆ U s with || F s || ≤ || F n || tan ϕ c ( 2
)
where n is the normal to the contact plane, δ n , the overlapping distance between spheres, and U s , the relative tangential (or shear) displacement at the contact point. Only compressive normal forces are modeled, and the contact is lost as soon as the overlap, δ n , vanishes.

Since spherical shapes lead to excessive rolling of particles, resulting in small internal friction angles at the macroscopic scale too far to be representative of the actual soil properties [START_REF] Plassiard | A spherical discrete element model: calibration procedure and incremental response[END_REF][START_REF] Widuliński | Numerical simulations of triaxial test with sand using dem[END_REF], contact rolling resistance is taken into account in the local constitutive law [START_REF] Oda | Experimental micromechanical evaluation of strength of granular materials: Effects of particle rolling[END_REF][START_REF] Iwashita | Micro-deformation mechanism of shear banding process based on modified distinct element method[END_REF][START_REF] Etsrada | Comparison of the effects of rolling resistance and angularity in sheared granular media[END_REF][START_REF] Aboul Hosn | Discrete numerical modeling of loose soil with spherical particles and interparticle rolling friction[END_REF]. The rolling resistance at the contact is defined by the rolling stiffness, K r , and the coefficient of rolling friction, η r (commonly called friction because of the mathematical form it takes although it does not directly involve the friction between two surfaces). Thus, the rolling moment at a contact point, M r , acting against the relative rolling rotation of particles, θ r , is expressed as:

∆ M r = -K r ∆ θ r with || M r || ≤ || F n || η r min(R 1 , R 2 ) (3)
where R 1 and R 2 are the radii of the two spheres in contact respectively.

115

Note that all details of the discrete numerical model implemented in this paper are discussed in [START_REF] Aboul Hosn | Discrete numerical modeling of loose soil with spherical particles and interparticle rolling friction[END_REF].

Interstitial fluid flow

As mentioned previously, the hydraulic loading governs the onset and development of internal erosion. Therefore, the interstitial fluid dynamics is solved to take into consideration the effect of the fluid phase in the erosion process.

The interstitial flow is described with a numerical method based on a pore-scale discretization in finite volumes (PFV), and assuming incompressible pore fluids. A detailed presentation of this approach can be found in [START_REF] Chareyre | Pore-scale modeling of viscous flow and induced forces in dense sphere packings[END_REF][START_REF] Catalano | Pore-scale modeling of fluidparticles interaction and emerging poromechanical effects[END_REF] and only a short description is given here.

In the PFV method, pore bodies are defined locally through a regular trian-gulation of the granular packing (i.e. the triangulation of the particle centers weighted by the particle radii). The void space being discretized by a tetrahedral mesh, each tetrahedron representing a pore, a finite volume formulation of the Stokes flow between adjacent pores is made possible. Dual to the triangulation, the Laguerre tessellation represents the connections/paths of the fluid between adjacent pores. Thereafter, viscous flow equations are upscaled at the pore level and approximated with a finite volume scheme. The flux q ij between two connected pores, i and j characterized by pressures p i and p j respectively, is related to the local pressure gradient and the pore-space geometry:

q ij = g ij p i -p j l ij (4)
where g ij is the hydraulic conductance of the throat between pores i and j, and l ij the length of this throat.

The pressure field is obtained by solving a linear system composed of the continuity equations applied to the fluid phase in each pore. The latter implies that the rate of volume change of one pore is equal to the sum of the fluxes through the four facets of the tetrahedral pore element. Considering the pore i, the continuity equation writes:

V f i = j=4 j=1 q ij = g ij p i -p j l ij (5)
where V f i is the volumetric strain rate of pore i deduced from the solid particle velocities at the four apexes of the tetrahedron i.

Finally, the forces exerted by the fluid on each particle can be derived from the pressure field and expressed as the contour integrals of the pressure p and of the viscous shear stress τ . For a particle k with a boundary surface Γ k , the fluid force F k reads:

F k = ∂Γ k p n ds + ∂Γ k τ n ds ( 6
)
where n is the unit normal to the surface Γ k .

The way to take into account this force induced by the water seepage on the solid particles to describe the suffusion phenomenon is presented in the next section. Numerical simulations were performed using YADE software [START_REF] Šmilauer | The Yade Project[END_REF] in which both the DEM and the PFV are implemented.

Numerical erosion procedure

A numerical extraction procedure is defined to mimic the suffusion process. Describing accurately this phenomenon, a fully coupled hydro-mechanical model is required. However, the typical time scale related to the transport of solid particles in the interstitial space under the action of a water seepage is quite large and it would require to simulate a relatively long physical time resulting in an expensive computational cost. Therefore, a simplified extraction procedure, involving a one-way fluid-solid coupling [START_REF] Aboul Hosn | A discrete numerical description of the mechanical response of soils subjected to degradation by suffusion[END_REF] and keeping a reasonable computational cost, is defined. It is based on the association of the DEM and PFV methods presented in the previous section. As illustrated by the numerical scheme shown in Figure 1, the proposed method consists in splitting the suffusion process in two steps: (i) the detachment of particles under the action of the water seepage and (ii) their transport through the interstitial space. These two steps are considered sequentially from an equilibrium state of the granular assembly reached under a prescribed stress state, i.e. an isotropic compression up to a given confining pressure or a triaxial compression up to a given stress deviator. Equilibrium is assessed from a dimensenionless unbalanced force. It is computed as the ratio of the mean resultant particle force to the mean contact force and tends to zero for a perfect static equilibrium. A value of 10 -3 was considered here as representative of a state sufficiently close to this limit.

Particle detachment induced by water seepage

Prior to the detachment step, the seepage flow through the granular medium for a fixed global hydraulic gradient is resolved from a single iteration of the PFV method. Then, fluid forces acting on each solid particle are computed and stored. As for the solid phase, periodic limit conditions are considered for the interstitial water flow [START_REF] Šmilauer | Reference manual[END_REF]. The detachment step consists in checking the equilibrium of solid particles under the combined action of inter-particle contact forces and fluid forces. For this purpose, DEM cycles are iterated by applying permanently on solid particles the fluid forces previously stored. At this stage seepage flow and fluid forces are not updated by taking into account new particle positions and velocities, this is why we call it a one-way coupling. Nevertheless, this simplification is quite reasonable since the idea is not to give enough time for particles to be transported (i.e. particle displacements and velocities stay low) but to check for each particle if the fluid force can be balanced by the contact forces transmitted by neighboring particles in a configuration close to its initial equilibrium state. Formally, this is done by comparing the magnitude of the unbalanced force on the particle p subjected to a fluid force F p f and contact forces F α c (where α spans over all the contacts involving the particle p):

F p unb = α∈p F α c + F p f , (7)
with the mean magnitude of the contact forces transmitted by neighboring particles on p, F α c , and considered as stabilizing forces. Finally, a particle is considered as detached when:

F p unb > λ F α c , for α ∈ p (8)
where λ represents the threshold of the unbalanced force ratio below which a particle is considered at equilibrium.

Perfect equilibrium of the granular assembly is never reached with the DEM 150 (it is reached only asymptotically) even without seepage forces. As one would not expect particle detachment when fluid forces are zero, λ in Equation 8 is chosen high enough in order not to detect detached particles in a granular assembly in a quasi-static state not subjected to a water seepage. However, it should be small enough to detect any particles out of equilibrium when a quasistatic state is sought under a hydraulic loading. Finally, there is still the question of how many DEM cycles should be run

Particle transport through interstitial space

The possibility of detached particles (as identified in the previous section) to be transported through the pore network is now checked. The migration of fine particles is controlled by the topology of the interstitial space. In other words, it depends on how the pores of different sizes are connected via constrictions.

Thus, a detached particle should be able to cross a sufficient number of connected pores and constrictions to leave the soil matrix and be finally eroded.

Simulating such a process with the fully coupled DEM-PFV method is possible [START_REF] Tejada | Role of blockages in particle transport through homogeneous granular assemblies[END_REF] but requires a relatively high computational time. Moreover, this study does not aim to describe thoroughly the suffusion process in itself, but to investigate the effect of erosion, by the suffusion of fine particles, on the mechanical behavior of the soil. Hence, a simplified transport criterion is followed.

This criterion is based on the constriction size distribution, CSD, from the set of particles including the coarse ones and the non-detached finer particles as identified in the previous detachment step (i.e. all particles not detached by the water seepage). The computation of the CSD is the sub-product of the PFV method since it results from the regular triangulation of the granular packing and associating each tetrahedron with a pore. The constriction radius is identified at the interface between two tetrahedra as the radius of the inscribed circle between the three spheres intersecting the facet. A controlling constriction size is extracted from the CSD and the transport criterion consists in comparing the size of each detached particle to this controlling constriction size. This latter is chosen based on the suggestion of [START_REF] Raut | Further advancement in filter criteria through constriction-based techniques[END_REF] in their filter criterion which corresponds to D c35 (the constriction size for which 35 % of the constrictions of the granular assembly are finer than this size). Then, a detached particle is considered as being transported through the interstitial space if its diameter is smaller than

D c35 .
Finally, a detached particle being transported is identified as an eroded particle and is erased from the granular assembly. DEM cycles are iterated again with the new particle configuration (i.e. without the eroded particles) and without any fluid forces. Due to the particles removal, the granular assembly may deform under the constant applied stress state. The simulation is run until a new equilibrium state (if any) is reached. Such a process can be repeated to pursue the erosion development, or a drained triaxial compression test can be simulated to evaluate the new mechanical properties of the eroded soil.

Suffusion development and soil deformation

Suffusion-like simulations were performed from an isotropic stress state (under a confining pressure of 100 kPa), following the numerical extraction procedure defined previously. The relatively narrow particle size distribution displayed in Figure 4 was used in this first approach to limit the number of discrete elements necessary to form a representative elementary volume (REV) and to keep a quite low computational cost. Note that, such a particle size distribution represents an internally stable soil as considered in previous studies such as [START_REF] Wood | Modelling mechanical consequences of erosion[END_REF][START_REF] Scholtès | Multiscale approaches to describe mechanical responses induced by particle removal in granular materials[END_REF]. Therefore, to be able to apply the previous erosion procedure, some simplifications are taken in what follows. These include: (1) Considering only the CSD of the coarse fraction of the granular assembly for the calculation of the controlling constriction size; (2)Assuming the fine fraction is that correspond- ing to particles whose diameter is smaller than D 50 . The extraction algorithm in Figure 1 has been applied repeatedly for increasing values of the prescribed hydraulic gradient, i. For i < 0.2, the hydraulic gradient is increased by steps of 0.05. For 0.2 ≤ i ≤3, an increment of 0.2 is applied. As i exceeds 3, the increment step increases to 0.5. For a given hydraulic gradient, the extraction is repeated until no more particles are eroded. Figure 5 shows, for a given hydraulic gradient (i = 1.2), the relative detached and eroded masses, Md and Me respectively (with Me = Me M0 and M e the mass of eroded particles), as long as the extraction procedure is repeated at this value of the hydraulic gradient. The hydraulic loading process followed in the simulations immitates the experimental multi-stage procedure [START_REF] Chang | A stress-controlled erosion apparatus for studying internal erosion in soils[END_REF] where the hydraulic gradient is increased by steps as well as to take into account the effect of hydraulic loading history on the suffusion development [START_REF] Sibille | A description of internal erosion by suffusion and induced settlements on cohesionless granular matter[END_REF]. Suffusion like simulations are resumed in Figure 6 showing the cumulative eroded mass ratio for the prescribed hydraulic gradients. Here, the total mass of eroded particles M e is distinguished from the eroded mass of active particles M ea (i.e, particles participating in the inter-granular force transfer and involved in at least two contacts), with M ea ∈ M e . Then, eroded mass ratios Me = M e /M 0 and Mea = M ea /M 0 are respectively defined.

Under the initial isotropic stress state, rattlers (or non-active particles) of the fine fraction represent 38 % (by mass) of the total mass of fine particles. By definition, these rattlers are considered as detached for any non zero hydraulic gradient and their ability to be transported depends only on the transport criterion. That is why Me is much higher than Mea and a sufficiently strong hydraulic gradient (here about i = 1) is necessary to destabilize active particles involved in the contact force network (stabilizing them with respect to the fluid action). Besides, for increasing hydraulic gradient, new particles can become inactive (due to particle re-arrangement during the first step in Figure 1) and constriction sizes may globally increase (due to the previous erosion step) giving the opportunity for new inactive particles to be eroded. Finally, this results in a modified post-suffusion particle size distribution and constriction size distribution due to the selective erosion of particles as displayed in Figure 7.

Mechanical response of the granular assembly during suffusion can be characterized in terms of volumetric strain. Under the external applied isotropic stress the sample may compact only when active particles are eroded (since inactive particles do not participate in the load bearing capacity of the assembly). Nevertheless, Figure 8 shows that the soil compaction is not concomitant with the erosion of active particles. For Mea > 0.2 %, the sample compresses significantly whereas below this threshold it presents negligible volumetric deformations, evidencing the negligible role in the contact force network of the active particles firstly eroded. However, as noticed from Figure 8, the porosity is always larger after erosion even if the soil shows a large contractive behavior. The creation of a more open micro-structure, induced by particle removal, is predominant over the sample compaction during suffusion, as shown by [START_REF] Wood | Discrete element modelling of soil erosion[END_REF][START_REF] Wood | Changing grading of soil: effect on critical states[END_REF][START_REF] Scholtès | Multiscale approaches to describe mechanical responses induced by particle removal in granular materials[END_REF]. Although the eroded sample is always looser than the virgin one, porosity change with erosion development is not monotonous here. For a high degree of erosion development (i.e. for Mea > 1 %) the eroded medium is slightly denser than that at erosion initiation ( Mea < 0.2 %) as a result of the competition between compaction due to particle rearrangement and pore space increase due to particle removal. However, even if the soil becomes denser and gains some strength, the soil had already suffered from large settlements which may lead to serious drawbacks at the scale of hydraulic structures.

Laboratory experiments, aiming to characterize soil mechanical properties of soil samples eroded by suffusion, face a major issue related to the heterogeneous state of the sample after erosion. The fine concentration may present a gradient in the main seepage direction [START_REF] Ke | Experimental investigations on suffusion characteristics and its mechanical consequences on saturated cohesionless soil[END_REF][START_REF] Chang | A stress-controlled erosion apparatus for studying internal erosion in soils[END_REF] (with a more important fine concentration in the downstream part of the sample) and the density itself may also be affected with possibly counterintuitive distributions as observed by [START_REF] Sibille | A description of internal erosion by suffusion and induced settlements on cohesionless granular matter[END_REF] showing a denser state in the upstream part of the sample where fine fraction is the lowest. However, interpretation of the triaxial compression tests performed to evaluate the post-erosion shear strength properties usually assumes a homogeneous state of the sample. The numerical erosion process defined here implicitly induced a homogeneous removal of eroded particles. For instance, the homogeneity of the granular assembly is checked in Figures 9 and10

Effects of suffusion on the macroscopic constitutive behavior

Drained triaxial compressions have been simulated to highlight how such induced modifications of the micro-structure, due to the removal of a fine fraction, affect the mechanical behavior of the soil. Results in terms of deviatoric stress ratio, q/p, and volumetric strains, ε v , are displayed in Figure 11 for the virgin sample (i.e. not eroded (NE)) and for the eroded ones having different degrees of erosion. The eroded samples are characterized by the highest hydraulic gradient, i, reached before performing the triaxial compression; the corresponding eroded masses can be deduced from Figure 6.

Globally, the eroded samples present a reduced shear strength with respect to the virgin sample together with a significant change of the volumetric strains tendancy. Two distinct groups are clearly identifiable. For sufficiently low hydraulic gradients (i ≤ 1), group 1, the eroded samples show the most important shear strength reduction with a well marked contractive behavior, whereas for higher hydraulic gradients, group 2 (i > 1), a part of the initial shear strength is recovered and samples present a first phase of contractant behavior followed by a slight dilatant behavior. Moreover, for the group 1, an atypical softeninghardening behavior (developing after an initial peak of q/p, from an axial stress of about 3%), and volumetric reduction happening in two steps (a strong contractancy developing from ε a ≈ 5% after an initial moderate volume reduction), are also observable.

Despite this particular point that will be discussed further in the paper, changes observed in the shear strength and volumetric deformation are consistent with the sample porosity reached after the erosion process. The porosity of samples of group 1 is more important (n ≈ 0.46) than for the samples of group 2 (n ≈ 0.43) resulting, for the latter, from the strong compaction experienced during the extraction process (cf. Fig. 8) as more active particles were eroded at higher hydraulic gradients.

Change of the peak friction angle with erosion development, represented via the eroded mass ratio of active particles Mea , is presented in Figure 12. A non-linear and non-monotonic relation is found, in good agreement with the post-erosion porosity following the same trend (Fig. 8). In other words, for the case studied here, the shear strength is greatly affected even for a low degree of development of erosion where the friction angle decreases by about 19 % (from 37.5 o to 30.2 o ). A low degree of development of erosion is slightly more damaging to the shear strength than an advanced degree of erosion development.

The latter should be interpreted carefully in the sense that, for a large erosion development, the granular assembly has already shown a very high contractancy during erosion (up to ε v ≈ 14% to 15%). Regarding a structure made of soil, as a dike or an embankment, if such a high volume reduction concerns a sufficiently large volume of soil, it could lead to an important settlement of the structure, possibly as damaging as a shear strength reduction.

Another interesting point highlighted by these numerical results concerns the mechanical behavior of a soil when only inactive particles are eroded by suffusion. In this case and for the current mechanical state, as inactive particles do not participate in the contact force network and thus in the load bearing of the granular assembly, no deformation is induced by the suffusion and the original contact force network remains unchanged. However, does it affect the mechanical properties of the sample? Or, in other words, may the inactive particles in the initial state have a role in the load bearing capacity when the medium is sheared and deformed? Eroded samples belonging to group 1 correspond to this case: the eroded mass of active particles is almost negligible ( Mea < 0.2%) whereas the inactive particles have been largely removed from the granular assembly ( Me ≥ 25%). This is confirmed by the vanishing volumetric strain occuring during erosion for this group. Nevertheless, as shown by Figures 11 and12, these eroded samples present mechanical properties largely degraded with respect to the virgin sample, in particular with a reduction of the peak internal friction angle from 37.5 o to about 30.2 o . Therefore, although some particles are inactive in a given mechanical state, they may play a key role in subsequent states with respect to the macroscopic behavior of the granular assembly (even if this role, for each particle, may be very transient). It is not beacause no deformation (or settlement) is observed during the suffusion occurence (if the mechanical loading applied to the soil is unchanged) that the soil mechanical properties have not been degraded by suffusion.

Finally, properties of eroded granular assemblies at the critical state are discussed in the last part of this section. Changes in the critical state may be of importance for many phenomenological constitutive relations relying on the definition of such a state. It can be deduced from the Figure 11 that the mobilized friction angle at large strains is unchanged (with q/p ≈ 1.2) whatever the hydraulic gradient is and the degree of erosion development. For the numerical model used, and for a given particle grading, shear strength at large deformations depends almost only on the inter-particle rolling friction (reflecting the angularity of the particle shapes) as shown for instance in [START_REF] Aboul Hosn | Discrete numerical modeling of loose soil with spherical particles and interparticle rolling friction[END_REF]. As the rolling friction η r (Eq. 3) is here identical for all the particles (the eroded and noneroded ones), the removal of a part of the particles does not affect the overall strength due to this resistance to rolling introduced at the contact level. Besides, the change in the gradation of the particle sizes due to the selective erosion of particles (Fig. 7) is apparently not sufficiently marked in this case to have an influence on the shear strength at the critical state. However, even if it is not very significant, this change in the gradation clearly affects the porosity reached at the critical state as shown in Figure 13. As the finest particles, possibly filling the voids between coarser particles, are removed, the porosity increased at the critical state. Unlike the post-erosion internal friction angle (Fig. 12) depending on the initial (post-erosion) density, the relation between the porosity at critical state and the eroded mass is monotonous since it relies on the change of the grading (and obviously not on the initial configuration of the granular packing) also changing monotonously with the erosion development (Figure 7).

Micro-structure changes induced by suffusion

In order to understand how the microstructure of the soil is modified during erosion and to explain the particular behavior of eroded soils subjected to low hydraulic gradients (group 1), the probability density function, P , of the normal contact forces is calculated for each eroded soil specimen and compared to that of the non-eroded soil. The probability density, P , is displayed in Figure 14 for the different specimens in the isotropic stress state. Furthermore the mechanical coordination number Z m is also computed for each specimen. According to [START_REF] Throton | [END_REF] Z m is defined as:

Z m = 2N c -N 1 N p -N 0 -N 1 (9)
where N c is the total number of contacts, N 1 is the number of particles with one contact, N 0 is the number of particles with zero contact, and N p is the 415 total number of particles. Probability density of contact forces and mechanical coordination number are displayed in Figure 14. An additional granular assembly denoted LS is also considered in this figure. LS is a granular assembly reconstituted in a loose state by using the PSD of the specimens of group 1.

The objective is to try to constitute a reference loose sample, possibly as loose as those of group 1, but without subjecting it to the erosion process (contrary to those of group 1). The methodology proposed in [START_REF] Aboul Hosn | Discrete numerical modeling of loose soil with spherical particles and interparticle rolling friction[END_REF] to create initially very loose samples by numerically imitating the moist tamping technique by adding inter-particle cohesion forces during the compaction phase has been followed.

However, the highest porosity we manage to reach for LS at the isotropic state is n=0.444. Note that this porosity is lower than the porosity of specimens of group 1 (slightly above than 0.460) but higher than the porosity of group 2 (ranging from 0.435 to 0.439). Hence, if specimens of group 1 are discarded in a first time, specimens NE, specimens of group 2, and LS are respectively in order from the denser to the looser. As expected this order appears also clearly in Figure 14 with respect to the force distributions (distributions are different for NE and LS and distributions of group 2 are in between) and the coordination number (Z m is decreasing as the porosity is increasing). However, it is interesting to note that specimens of group 1 do not follow this expected trend and constitute a very particular case. Although they are the loosest, the erosion process has lead to unique granular assemblies, very loose compared to the initial non-eroded one, NE, but characterized by very similar coordination numbers and quasi identical contact force distributions. Specimen LS shows it is not possible to reconstruct at low initial density, without using the erosion process, such a micro-structure typical of a dense sample.

In addition, Figure 15 shows the response to a drained triaxial compression of the reconstituted loose sample compared with the other soil specimens. Here again, as expected from their intermediary densities, specimens of group 2 present volumetric changes which are between the clearly dilatant behavior of specimen NE and the globally contractant response of LS. With similar probability density of contact forces and Z m to the non-eroded granular assembly NE at isotropic state, specimens of group 1 show a kind of transitional volumetric trend. It is characterized by, apparently slightly dilatant behavior with peak stress after the initial contractancy, and as the shearing increases, the hidden loose character appears again by exhibiting a large contractant behavior (more important than for LS specimen) in accordance with the very high initial porosity. Finally, as the LS specimen is looser than the NE specimen it is characterized by a lower initial stiffness as showed in Figure 15. This is not the case for specimens subjected to the erosion process (both groups 1 and 2) where the initial stiffness stay quasi-unchanged with respect to the non-eroded one.

All these observations illustrate that suffusion may result in atypical soil microstructures and macroscopic behaviors, particularly when erosion concerns only inactive particles without accompanying effects on the coarse matrix integrity.

Conclusion

A numerical approach has been presented in this paper to describe the erosion of soil by suffusion. A simplified extraction procedure is defined based on the discrete element method with a one-way fluid-solid coupling. By taking into account the hydraulic loading, to determine the detachment of the solid particles, and the constrictions of the interstitial space, the definition of this extraction procedure constitutes new advances while limiting the computational cost. The simulations of the erosion tests show that depending on the role and the amount of the eroded particles, two different post erosion mechanical behaviors may be concluded. On one hand, if only inactive particles are removed, the soil structure presents almost no deformation during erosion, maintaining the same probability density of contact forces and mechanical coordination number as the non-eroded soil. However, if sheared, this soil may suffer from a sudden collapse revealing its hidden loose character induced by the voids left at the place of the eroded particles. On the other hand, if more active particles are eroded, the granular assembly loses its equilibrium during the erosion process and the soil deforms with a global volume reduction to reach another equilibrium state.

Once sheared, the soil still presents a slightly dilatant behaviour (not as marked as in the non eroded state) and recovers part of its shear strength with respect to the previous case. In all cases, the erosion process leads to the creation of more porous materials (even if it may be partially compensated by a compaction during the erosion). However, a comparison with a granular assembly generated in an initial loose state showed that the post-erosion mechanical behaviour cannot be only described by an increase in porosity. It results in atypical microstructures and macroscopic behaviors, particularly when only inactive particles have been eroded.

It would be interesting in future work to consider a wider PSD similar to those used experimentally and which are internally unstable. For such PSD, the numerical granular assembly should include at least several hundreds of thousands of particles to constitute a representative elementary volume. Moreover, a limitation of the proposed model is the simplified transport criterion considering a single characteristic constriction size. The ability of the detached particles to be transported inside the porous network is determinant in the development of suffusion. This feature will be evaluated more thoroughly since the whole constriction size distribution and pore connectivity is available from the numerical model.

Acknowledgement

The environment and support provided by the French research group GDR MeGe 3176 is gratefully acknowledged. Laboratory 3SR is part of the LabEx Tec 21 (Investissements dAvenir -grant agreement n • ANR-11-LABX-0030)

Figure 1 :

 1 Figure 1: Dissociation of the detachment step and transportation step in the simplified numerical description of suffusion.

Figure 2

 2 presents the mass ratio Md = M d /M 0 (where M d is the mass of detached particles and M 0 the initial mass of fine particles) identified as detached in the absence of water seepage for the discrete numerical model (as defined in Section 4) in a quasi-equilibrium state. For λ ≥ 0.1, the mass of detached particles is negligible, therefore, the lower boundary λ = 0.1 is chosen as the threshold unbalanced force ratio for the next erosion simulations.

Figure 2 :

 2 Figure 2: Mass ratio of particles identified as detached for different values of the threshold λ applied on the unbalanced force ratio for a nil hydraulic gradient (i.e. without water seepage simulated in the granular assembly).

Figure 3 :

 3 Figure 3: Number of detached particles for different DEM cycles.

Figure 4 :

 4 Figure 4: Initial particle size distribution of the granular assembly.

Figure 5 :

 5 Figure 5: Relative masses of detached and eroded particles at each step of repetition of the extraction procedure until no more particles are eroded under the prescribed hydraulic gradient, i = 1.2.

Figure 6 :

 6 Figure 6: Simulated eroded mass for an increased by step of the prescribed hydraulic gradient i.

Figure 7 :

 7 Figure 7: Post-erosion Particle Size Distribution (PSD) and Constriction Size Distribution (CSD).

Figure 8 :

 8 Figure 8: Contractive volumetric deformation of granular assembly and the evolution of porosity with the erosion of active particles.

  after it has been eroded for a hydraulic gradient up to i = 3. Porosity and grain size distribution were computed around different cross sections of the sample (bottom/downstream, middle, and top/upstream). Both porosity and grading are almost identical in all sections. Consequently, the modification of the micro-structure of the numerical sample, being kept homogeneous, is indeed not directly representative of relatively large samples eroded in laboratory, but presents the advantage to get around this problem of interpretation of mechanical tests on heterogeneous samples. This constitutes an important basis for the next section of the paper about the post-erosion mechanical properties.

Figure 9 :

 9 Figure 9: Local porosities at different sections of the eroded sample under a hydraulic gradient, i = 3.

20

 20

Figure 10 :

 10 Figure 10: Local particle size distributions around different sections of the eroded sample under a hydraulic gradient, i = 3.

Figure 11 :

 11 Figure 11: Simulated responses to drained triaxial compression of granular assemblies eroded at different hydraulic gradients.

Figure 12 :

 12 Figure 12: Maximum friction angle ϕmax in terms of eroded mass ratio of active particles, determined from simulations of triaxial compressions on eroded samples.

Figure 13 :

 13 Figure 13: Porosity at critical state ncs in terms of eroded mass ratio of active particles, determined from simulations of triaxial compressions on eroded samples.

Figure 14 :

 14 Figure 14: Probability density functions of normal contact forces at the isotropic state of the eroded specimens up to a hydraulic gradient i and non-eroded initially dense NE and loose LS specimens.

Figure 15 :

 15 Figure 15: Comparison of the response of an initially loose sample (LS) to a drained triaxial compression to other granular assemblies eroded at different hydraulic gradients.


Lorenzo Brandolese

Fernando Cortez

BLOWUP FOR THE NONLINEAR HEAT EQUATION WITH SMALL INITIAL DATA IN SCALE-INVARIANT BESOV NORMS

Keywords: Mathematics Subject Classification. 35K55 (primary), 30H25 (secondary) Blowup, Nonlinear Heat Equation, Besov. 1

Blowup for the nonlinear heat equation with small initial data in scale-invariant Besov norms

Lorenzo Brandolese, Fernando Cortez

Introduction

In this paper we study the Cauchy problem for the nonlinear heat equation (1.1)

∂ t u = ∆u + |u| α u, x ∈ R n , t ∈ [0, T ] u(0, x) = u 0 (x),
where α > 0, 0 < T ≤ ∞, and u : R + × R n → R is a real function. This problem attracted a considerable interest and we refer to, e.g., [1, 6-9, 12, 15, 16, 19, 25, 28] for a small sample of the huge existing literature. Several well-posedness results are available for the Cauchy problem (1.1). For example, if u 0 ∈ C 0 (R n ), then there is T = T (u 0 ) > 0 and a unique u ∈ C([0, T ), C 0 (R n )) which is a classical solution to (1.1) on (0, T ) × R n . For more singular data, say u 0 ∈ L p (R n ), we know the following, see [START_REF] Brezis | A nonlinear heat equation with singular initial data[END_REF][START_REF] Weissler | Local existence and nonexistence for semilinear parabolic equations in L p[END_REF][START_REF]Existence and non-existence of global solutions for a semilinear heat equation[END_REF].

-When p > nα 2 and p ≥ 1, there exists a constant T = T (u 0 ) > 0 and a unique function u(t) ∈ C([0, T ], L p (R n )) that is a classical solution to (1.1) on (0, T ) × R n . -When p < nα 2 , there is no general theory of existence. Besides, A. Haraux and F. Weissler [START_REF] Haraux | Non-uniqueness for a semilinear initial value problem[END_REF] established the non-uniqueness, by showing that there is a positive solution in C([0, T ],

L p (R n )) ∩ L ∞ loc ((0, T ), L ∞ (R n ))
, arising from zero initial data. -When p = nα 2 , see Theorem 2.1 below. We will be interested in the issues of the blowup in finite time v.s. the global existence of the solutions. The first works in this direction are due to H. Fujita. Fujita proved that for the positive solutions of (1.1), if the initial data u 0 is of class C 2 (R n ) with derivatives up to the second order bounded on R n , then a necessary condition for u to be unique in

C(R n × [0, T )) is that ∀ x ∈ R n , |u 0 (x)| ≤ M e |x| β ,
for some constants M > 0 and 0 < β < 2. See [START_REF] Fujita | On the blowing up of solutions of the Cauchy problem for ut = ∆u + u 1+α[END_REF][START_REF]On some nonexistence and nonuniqueness theorems for nonlinear parabolic equations[END_REF].

About the problem of the existence of regular global solutions, there are two possible scenarios: if nα/2 < 1, then no nontrivial positive solution of this problem can be global (a situation now referred as Fujita's phenomenon), while for nα/2 > 1, there are global non-trivial positive solutions under small initial data assumptions. K. Hayakawa [START_REF] Hayakawa | On nonexistence of global solutions of some semilinear parabolic differential equations[END_REF] and F. Weissler [START_REF] Weissler | Local existence and nonexistence for semilinear parabolic equations in L p[END_REF][START_REF]Existence and non-existence of global solutions for a semilinear heat equation[END_REF] later proved that Fujita's phenomenon occurs in the case of the critical exponent nα/2 = 1.

Motivations and overview of the main result

To motivate our results, we introduce the concept of a scale-invariant space. For λ > 0, let us set

u λ (t, x) = λ 2 α u(λ 2 t, λx) and u 0,λ (x) = λ 2 α u 0 (λx). (2.1)
For every solution u(t, x) of (1.1), u λ (t, x) is also a solution of (1.1) with initial data u 0,λ (x). In this setting, we say that a Banach space E is scale-invariant, if

u(t, •) E = u λ (t, •) E . (2.2)
Scale-invariant space are known to play an essential role in issues like well-posedness, global existence or blow-up of the solution.

The purpose of the present paper is to study the borderline cases of explosion and global existence for solutions of (1.1), in a scale-invariant Banach space. In the case of problem (1.1), the only L p (R n )-space invariant under the above scaling (2.1) is obtained for p = nα/2. Notice that p ≥ 1 if and only if α is larger or equal to the Fujita critical exponent. Therefore, we will be especially interested in solutions in L nα/2 (R n ).

Our starting point is the following theorem, where we collect some of the results of Brezis, Cazenave and Weissler, in this scaling invariant setting.

Theorem 2.1 (See [START_REF] Weissler | Local existence and nonexistence for semilinear parabolic equations in L p[END_REF]. See also [START_REF] Brezis | A nonlinear heat equation with singular initial data[END_REF] for the uniqueness). Let u 0 ∈ L nα/2 (R n ), and assume that nα/2 > 1. There exists a time

T = T (u 0 ) > 0 and a function u ∈ C([0, T ], L nα/2 (R n ))∩ L ∞ loc ((0, T ], L ∞ (R n )) such that u is a classical solution of (1.1) on (0, T ) × R n . Moreover, (i) sup 0<t<T t σ/2 u(•, t) p < +∞, (ii) lim t→0 t σ/2 u(•, t) p = 0, for any nα 2 < p < nα(α+1) 2 and σ = 2 α -n p . The uniqueness of classical solutions to (1.1) holds in the class C([0, T ], L nα/2 (R n )).
Moreover, there exists δ = δ(α, n) such that if u 0 nα/2 < δ then such solution is global, i.e., one can take T arbitrarily large.

The solution of Theorem 2.1 satisfies the integral equation

(2.3) u(t) = e t∆ u 0 (x) + t 0 e (t-s)∆ |u| α u(s) ds,
where e t∆ f = G t * f and

G t (x) = (4πt) -n/2 e -|x| 2 /(4t)
is the standard Gaussian. The uniqueness of weak solutions to the integral equation (2.3) have been also addressed in [START_REF] Brezis | A nonlinear heat equation with singular initial data[END_REF]. The authors show that there is at most one weak solution to (3.1) in the class C([0, T ],

L nα/2 (R n )) ∩ L ∞ loc ((0, T ], L ∞ (R n )), provided n 2 α ≥ α + 1 > 1.
Under the additional restriction n 2 α > α + 1 > 1 the uniqueness of weak solutions to (3.1) holds in the larger class C([0, T ], L nα/2 (R n )). We refer to E. Terraneo's paper [START_REF] Terraneo | Non-uniqueness for a critical non-linear heat equation[END_REF] for further uniqueness/non uniqueness results of weak solutions.

The problem of obtaining global-in-time solutions by relaxing the stringent smallness assumption u 0 nα/2 < < 1 was also addressed. New ideas in this direction were brought by M. Cannone and Y. Meyer's works on the Navier-Stokes equations [START_REF] Cannone | Harmonic analysis tools for solving the incompressible Navier-Stokes equations[END_REF][START_REF] Meyer | Oscillating Patterns in Some Nonlinear Evolution Equations[END_REF].

In the model case α = 2 and n = 3, i.e. for the cubic heat equation in R 3 ,

∂ t u = ∆u + u 3 , Y.
Meyer observed in his lecture notes [START_REF] Meyer | Oscillating Patterns in Some Nonlinear Evolution Equations[END_REF] that if u 0 ∈ L 3 (R 3 ), with

u 0 Ḃ-1/2 6,∞ < < 1,
(this condition is considerably weaker than requiring u 0 3 < < 1) then the maximal time T * of the solution is T * = +∞. In fact, the method described therein would go through provided u 0 Ḃ-1+3/p p,∞ < < 1 and 3 < p < 9. In [START_REF] Meyer | Oscillating Patterns in Some Nonlinear Evolution Equations[END_REF], he also raised the question whether or not, for u 0 ∈ L 3 (R 3 ), the even weaker smallness condition u 0 Ḃ-1 ∞,∞ < < 1 would still imply T * = +∞. See next section for the definition of Besov spaces. Notice that these Besov spaces enjoy the same scaling invariance properties as L 3 (R 3 ) and we have the continuous injections

L 3 (R 3 ) ⊂ Ḃ-1+3/p p,∞ (R 3 ) ⊂ Ḃ-1 ∞,∞ (R 3 ) (3 < p ≤ +∞). Moreover, Ḃ-1 ∞,∞
is known to be the largest function space invariant under translation and satisfying such scaling property. In this sense, a smallness condition on the Ḃ-1 ∞,∞ -norm would be the least demanding restriction that one could put in a scale-invariant setting.

In the same spirit, but for the general case of problem (1.1), the best result for the global-in-time existence are due to Miao, Yuan, and Zhang [START_REF] Miao | Strong solutions to the nonlinear heat equation in homogeneous Besov spaces[END_REF]. They proved (among other things) that the solution of Theorem 2.1 is global, provided u 0 ∈ L nα/2 (R n ), with nα/2 > 1, under the smallness condition

u 0 Ḃ-2/α+n/p p,q < < 1, for some 1 < nα 2 < p < nα(α+1) 2 , 1 ≤ q ≤ ∞.
The restriction nα 2 < p, together with the condition q ≥ p, ensure the embedding of

L nα/2 (R n ) into Ḃ-2/α+n/p p,q (R n ).
On the other hand the authors of [START_REF] Miao | Strong solutions to the nonlinear heat equation in homogeneous Besov spaces[END_REF] left open the limit case p = nα(α + 1)/2. In other words, they left open the question whether or not initial data u 0 ∈ L nα/2 (R n ), small in the Ḃ-2/(α+1) nα(α+1)/2,q -norm, give rise to global-in-time solutions. Our main result below provides a negative answer to the above problem, thus settling the borderline problem of the global solvability of (2.4), at least in the case of integer nonlinearity exponents.

More specifically, for b ∈ N, we consider the Cauchy problem for the non-linear heat equation

(2.4) ∂ t u = ∆u + u b u(0, x) = u 0 (x) x ∈ R n , t ∈ (0, T ).
where 0 < T ≤ ∞. The results recalled for the problem (1.1) -in particular Theorem 2.1remain valid for (2.4), with b = α+1. These two Cauchy problems in fact agree for positive solutions, or for real solutions of any sign, when b is an odd integer.

Theorem 2.2. For any δ > 0 and b < q ≤ +∞, with b ∈ N and n(b -1)/2 > 1, there exists

u 0 ∈ S(R n ) (the Schwartz class) such that (2.5) u 0 Ḃ-2/b nb(b-1)/2,q ≤ δ,
and such that the maximal time

T * of the solution u ∈ C([0, T * ), L n(b-1)/2 (R n )) to (2.4)
arising from u 0 is finite.

As mentioned in the introduction, for n(b -1)/2 ≤ 1, because of Fujita's phenomenon, finite time blow up occurs for positive solutions, no matter which norms of the initial data are assumed to be small.

In the case b = n = 3, Theorem 2.2 negatively answers Y. Meyer's question [18, Conjecture 1].

There are several blowup results for (1.1) based on the maximum principle, energy functionals, concavity methods, or the spectral properties of the Laplacian, etc. See e.g. [START_REF] Bandle | Blowup in diffusion equations: a survey[END_REF] for a review of these classical methods. But none of them seems to be effective to establish Theorem 2.2, as the smallness condition (2.5) represents a severe obstruction for their applicability.

The proof of Theorem 2.2 is constructive: suitable initial data are given by (4.5) below with N = N (δ) large enough. Our approach, inspired by Palais [START_REF] Palais | Blowup for nonlinear equations using a comparison principle in Fourier space[END_REF], rather uses the positivity of the Fourier transform inherited from its initial condition u 0 . Even though conceptually similar to [START_REF] Palais | Blowup for nonlinear equations using a comparison principle in Fourier space[END_REF], our paper is technically completely different (for example, we are able to remove the restriction b ≤ 1 + 2/n that appears therein). From the technical point of view, our paper is somehow closer to [START_REF] Lemarié-Rieusset | Sobolev multipliers, maximal functions and parabolic equations with a quadratic nonlinearity[END_REF][START_REF] Montgomery-Smith | Finite time blow up for a Navier-Stokes like equation[END_REF], where the authors studied the blowup for different equations, namely diffusion problems with nonlocal quadratic nonlinearity. Our method bears also some relation with that of [START_REF] Ru | The blow-up solutions of the heat equations in FL 1 (R n )[END_REF]. However, the blowup result in F(L 1 ) of [START_REF] Ru | The blow-up solutions of the heat equations in FL 1 (R n )[END_REF] is not put in relation with the size of the data in scale-invariant norms. As such, our blowup result looks more precise, and its proof shorter.

Since the work of Cannone [START_REF] Cannone | Harmonic analysis tools for solving the incompressible Navier-Stokes equations[END_REF] we know that fast enough oscillations of the initial give rise to global-in-time smooth solutions for a large class of semilinear dissipative system, and that size conditions on Besov norms with negative regularity represent an effective way to measure such oscillations. The main interest of our result is to illustrate a limitation of this principle, by showing that there are scale invariant Besov norms that turn out to be too weak to be used for this purpose.

In our blowup result, the maximal time T * can be taken arbitrarily small, as one easily checks applying Theorem 1 to rescaled data u 0,λ , that have the same Besov norm as in (2.5), and existence time T * λ = λ -2 T * . In the more difficult case of the Navier-Stokes equations, a similar problem was addressed by Bourgain and Pavlović [START_REF] Bourgain | Ill-posedness of the Navier-Stokes equations in a critical space in 3D[END_REF] (see also [START_REF] Yoneda | Ill-posedness of the 3D-Navier-Stokes equations in a generalized Besov space near[END_REF]). These authors considered the Cauchy problem for Navier-Stokes with small data in Ḃ-1 ∞,∞ . While they left open the hard problem of the blowup, they succeeded in constructing a solution featuring a "norm-inflation" phenomenon in such Besov space, after an arbitrarily short time. But it was later realized by O. Sawada [START_REF] Sawada | A remark on the Navier-Stokes flow with bounded initial data having a special structure[END_REF] that Bourgain's-Pavlović solution, in fact, does not blow up in finite time.

Preliminaries

Let us recall the definition of the Besov norms and the Littlewood-Paley decomposition:

let ψ ∈ S(R n ) such that supp ψ ⊂ {3/4 ≤ |ξ| ≤ 8/3} and 1 = ∞ j=-∞ ψ j (ξ) (ξ ∈ R n , ξ = 0),
where ψ j (x) = 2 nj ψ(2 j x), j ∈ Z. Here and throughout, f denotes the Fourier transform of f . The homogeneous Besov spaces Ḃs.p q can be defined as follows, at least for s < n/p and 1 ≤ p, q ≤ ∞, which will be our case (in this paper we will only deal with the case s < 0):

Ḃs p,q = {f ∈ S ′ (R n ) : f = j∈Z (2 js ψ j * f ) in the S ′ (R n )-sense, and f Ḃs,p q < ∞},
where, for 1 ≤ q < +∞,

f Ḃs p,q =   ∞ j=-∞ 2 js ψ j * f q p   1/q
and f Ḃs p,∞ = sup q∈Z 2 js ψ j * f p . As mentioned before, one can obtain in Theorem 2.1 the global existence of the solution, dropping the smallness assumption on u 0 n(b-1)/2 , and putting instead a smallness assumption on the Ḃ-2/(b-1)+n/p p,q -norm of the data, which is weaker than the L n(b-1)/2norm. Let us sketch a proof of this fact, following the arguments of [START_REF] Meyer | Oscillating Patterns in Some Nonlinear Evolution Equations[END_REF][START_REF] Miao | Strong solutions to the nonlinear heat equation in homogeneous Besov spaces[END_REF], putting in evidence the admissible range for p, which is n(b -1)/2 < p < nb(b -1)/2.

One rewrites Equation (2.4) in the equivalent Duhamel formulation

(3.1) u(t, x) = e t∆ u 0 (x) + t 0 e (t-τ )∆ u b (τ, x) dτ =: Φ(u)(t, x). If u 0 ∈ L n(b-1)/2 (R n ), then the solution u ∈ C([0, T ], L n(b-1)/2 (R n )) of Theorem 2.1 (recall that α = b -1
) is obtained through the contraction mapping theorem, as the limit u = lim u l of approximate solutions (where u 1 = e t∆ u 0 and u l+1 = Φ(u l ), for l = 1, 2, . . .) in the X-norm, where

u X := sup 0<t<T u(t) n(b-1)/2 + sup 0<t<T t 1/(b-1)-n/(2p) u(t) p =: u Y + u Z .
Indeed, first notice that e t∆ u 0 ∈ X by standard heat kernel estimates. Next, the key estimates for the nonlinear term are the following:

t 0 e (t-s)∆ u b (s) ds p ≤ C t 0 (t -s) -n 2 ( b p -1 p ) u b (s) p/b ds ≤ C u b Z t 0 (t -s) -n 2 ( b p -1 p ) s -b/(b-1)+nb/(2p) ds ≤ C u b Z t -1/(b-1)+n/(2p) , (3.2)
and

t 0 e (t-s)∆ u b (s) ds n(b-1)/2 ≤ C t 0 (t -s) -n 2 ( b p -2 n(b-1) ) u b (s) p/b ds ≤ C u b Z . (3.3)
These estimates are valid when 1 < n(b -1)/2 < p < nb(b -1)/2 (one also needs here 1 < b ≤ p, but the restriction b ≤ p can be dropped after the solution is constructed, by interpolation). These estimate ensure that

Φ(u) X ≤ C u b Z .
The Lipschitz estimates

Φ(u) -Φ(v) X ≤ C( u b-1 Z + v b-1 Z ) u -v Z ,
is established in a similar way. But u 0 Ḃ-2/(b-1)+n/p p,∞ ≃ u 1 Z owing to the heat kernel characterization of Besov spaces (see [START_REF] Cannone | Harmonic analysis tools for solving the incompressible Navier-Stokes equations[END_REF]). Hence, starting with u 0 small enough in the Ḃ-2/(b-1)+n/p p,∞ -norm allow to construct a solution with maximal lifetime T * = +∞. Without any smallness assumption, a well known variant [START_REF] Brezis | A nonlinear heat equation with singular initial data[END_REF][START_REF] Weissler | Local existence and nonexistence for semilinear parabolic equations in L p[END_REF] of the above argument still allows to construct a solution u = lim u l in the X-norm, at least when T > 0 is small enough. This relies on the observation that the approximate solutions (and hence the solution u itself) satisfy the additional condition lim t→0 t 1/(b-1)-n/(2p) u l (t) p = 0, for all l.

Proof of main theorem

We start with a simple general remark about the properties of solutions of Theorem 2.1, arising from initial data in the Schwartz class. In this case, or more in general when u 0 ∈ L 1 ∩L n(b-1)/2 , the corresponding solution obtained in Theorem 2.1 remains in L 1 (R n ) during the whole lifetime of the solution, and u ∈ C([0, T ], L 1 (R n )). This could be seen applying Gronwall-type estimates, or otherwise with the following argument: our claim is immediate if b > n(b -1)/2. Indeed, in this case we may take p = b in Theorem 2.1, and we have sup 0<t<T t 1/(b-1)-n/(2b) u(t) b < ∞. So,

u(t) 1 = Φ(u)(t) 1 ≤ u 0 1 + C(u 0 ) t 0 s -b/(b-1)+n/2 ds ≤ C(u 0 , T ),
because our condition n(b -1)/2 > 1 ensures that the above integral is finite for all finite T > 0. Moreover, the continuity with respect to t is obvious. On the other hand, if n(b -1)/2 ≥ b then we first observe that u b ∈ C([0, T ], L n(b-1)/(2b) ), next that e t∆ u 0 ∈ L n(b-1)/(2b) by interpolation, and from the integral equation u(t) = Φ(u)(t) we deduce u ∈ C([0, T ], L n(b-1)/(2b) ∩ L n(b-1)/2 ). If n(b -1)/(2b 2 ) ≤ 1, then using again the equation u = Φ(u) we get by interpolation that

u b ∈ C([0, T ], L 1 (R n )) and so u ∈ C([0, T ], L 1 (R n )).
Otherwise we iterate this argument, until we find m ∈ N such that n(b -1)/(2b m ) ≤ 1 and we conclude as before.

In the same way, going back to the sequence (u l ) of approximate solutions introduced in the previous section, one can prove that when u 0 ∈ L 1 ∩L n(b-1)/2 not only the convergence u l → u holds in the X-norm, but also u l → u in the C([0, T ], L 1 )-norm, as l → ∞.

Later on we will choose a specific u 0 ∈ S(R n ) such that u 0 ≥ 0 and u 0 is even (in a such way that u is real-valued). All the approximate solutions u l constructed from such datum u 0 satisfy u l (t, •) ≥ 0. The convergence of (u l ) in the C([0, T ], L 1 (R n ))-norm implies that u(t, •) ≥ 0 during the whole lifetime of the solution.

We introduce the following notation: for b ∈ N, and a non-negative measurable function f , we denote Let us now state a useful lemma.

Lemma 4.1. Let δ > 0, b ∈ N (b ≥ 2) and w ∈ S(R n ), such that w ≥ 0. Let c δ = 1 -e -δ 2 (b 2 -1)
. Also assume that the support of w is contained in the ball B(0, 1). Let w k , α k and t k be defined by the recursive relations (k ≥ 1):

w 0 = w w k = w b k-1 , α 0 = 1 α k = α b k-1 b -2k c δ , t 0 = 0 t k = t k-1 + b -2k δ 2 (b 2 -1).
Then, if u is the solution of (3.1) with initial condition u 0 (x) ∈ L n(b-1)/2 (R n ), and if u 0 (ξ) ≥ A w(ξ) with A > 0, then, for any k ∈ N,

(4.2) u(t, ξ) ≥ A b k α k e -b k t 1 t≥t k w k (ξ),
where

1 t≥t k is the indicator function of the interval [t k , +∞).
Proof. Using Fourier transform, we have that (3.1) becomes

u(ξ, t) = e -t|ξ| 2 u 0 (ξ) + t 0 e (s-t)|ξ| 2 [ u(s, ξ)] * b ds. (4.3)
We start with the case k = 0. We have u(t, •) ≥ 0, because u 0 (ξ) ≥ A w(ξ) ≥ 0, as observed at the beginning of this section. Then, using that supp w ⊂ {|ξ| ≤ 1}, we get u(ξ, t) ≥ e -t|ξ| 2 u 0 (ξ) ≥ e -t|ξ| 2 A w(ξ) ≥ A e -t w(ξ)

= A α 0 e -t w 0 (ξ), ∀t ≥ 0. (4.4)
This agrees with (4.2) for k = 0. Suppose now that inequality (4.2) holds for k -1. Then we get, for all t ≥ t k :

u(ξ, t) ≥ t 0 e (s-t)|ξ| 2 [ u(s, ξ)] * b ds ≥ t t k-1 e (s-t)|ξ| 2 (A b k-1 α k-1 ) b e -b k s [ w k-1 (ξ)] * b ds. ≥ A b k α b k-1 w k (ξ) t t k-1 e (s-t)|ξ| 2 e -b k s ds ≥ A b k α b k-1 w k (ξ)e -b k t t t k-1 e (s-t)b 2k ds,
where in the last inequality we used that supp

w k ⊂ {|ξ| ≤ b k }. But, for t ≥ t k , we have t t k-1 e (s-t)b 2k ds = b -2k (1 -e -b 2k (t-t k-1 ) ) ≥ b -2k c δ , because t k -t k-1 = b -2k δ 2 (b 2 -1)
, and so 1 -e -b 2k (t k -t k-1 ) = c δ . Hence we get,

u(ξ, t) ≥ A b k α b k-1 b -2k c δ e -b k t 1 t≥t k w k (ξ) ≥ A b k α k e -b k t 1 t≥t k w k (ξ),
by the recursive relation defining α k . Our claim now follows by induction.

For later use, let us observe that closed form for the sequences introduced in the previous lemma w k , α k and t k , are

w k = w b k (k ≥ 0), next α k = b -2b (b-1) 2 b k + 2 b-1 k+ 2b (b-1) 2 c b k -1 b-1 δ (k ≥ 0), and
t k = δ 2 (b 2 -1) k j=1 b -2j (k ≥ 1)
as it is easily checked.

Next lemma provides a first blowup result for equation (3.1).

Lemma 4.2. Let δ > 0 and w ∈ S(R n ) (w = 0) be a Schwartz function such that w ≥ 0. Also assume that the support of w is contained in the ball B(0, 1). Let u 0 ≥ Aw, with

A ≥ b 2b/(b-1) 2 c -1/(b-1) δ e δ/2 w -1 1 . If u is the classical solution of (3.1) arising from u 0 and belonging to C([0, T * ], L n(b-1)/2 (R n )), then T * ≤ δ 2 .
Proof. Assume, by contradiction, T * > δ 2 . Applying Lemma 4.1, and using that t k ↑ δ 2 as k → +∞, we get, for t = δ/2 and all k ∈ N,

u(δ/2, •) 1 ≥ A b k α k e -b k δ/2 w k 1 = A b k α k e -b k δ/2 w b k 1 ,
by Tonelli's theorem and the non-negativity of w. The size condition on A ensures that, taking sup k∈N in the right-hand side, one gets u(δ/2, •) 1 = +∞. But by the positivity of u(t, •) and Fourier inversion formula,

u(δ/2, •) L ∞ ≥ (2π) -n u(δ/2, •) L 1 = +∞.
This contradicts the fact that the lifetime of Weissler solution satisfies T * > δ/2.

Proof of Theorem 2.2. Let δ > 0 fixed and w ∈ S(R n ) such that w = 0 and w ≥ 0. We also assume that w is an even function and its support is contained in the ball B(0, 1 2b ). Let also u 0,N ∈ S(R n ) be defined as

(4.5) u 0,N (x) = ǫ N N k=0 2 2k/b η k cos( 3 2 2 k x 1 ) w(x) (N ∈ N),
where the sequences (η k ) and (ǫ N ) are chosen in the following way:

(4.6) η k = 1/(1 + k) 1/b , ǫ N = 1/ log(log(3 + N )).
In fact, the only thing that does matter in what follows are the following properties of (η k ) and (ǫ N ): they must be nonnegative and such that (η k ) ∈ ℓ q if and only if q > b, and Thus, if q ≥ 1 we get

ǫ N → 0,
u 0,N Ḃ-2/b n(b-1)b/2,q = ǫ N   j∈Z 2 -2 b qj ∆ j u 0,N q n(b-1)b/2   1 q = ǫ N   N j=0 η q j w(x) cos( 3 2 2 j x 1 ) q n(b-1)b/2   1 q ≤ ǫ N   N j=0 η q j w q n(b-1)b/2   1 q = ǫ N   N j=0 η q j   1/q w n(b-1)b/2 .
Thus, for any fixed δ > 0 and q > b, we can find N 0 ∈ N such that

u 0,N 0 Ḃ-2/b n(b-1)b/2,q < δ. (4.7)
Now, let T * N be the maximal time of the solution obtained in Theorem 2.1, arising from u 0,N (x). We denote u N this solution. We are going to prove that T * N < +∞, and more precisely that T * N < δ. If N ≥ N 0 and T * N < δ then there is nothing to prove. We thus pick N ≥ N 0 and assume T * N ≥ δ. By the remark at the beginning of Section 4, we have u N (t, ξ) ≥ 0 for all t ∈ [0, T * N ). Thus, if 0 < t < T * N , we get

u N (t, ξ) = e -t|ξ| 2 u 0,N (ξ) + t 0 e -(t-s)|ξ| 2 [ u N (s, ξ)] * b ds ≥ e -t|ξ| 2 u 0,N (ξ)
We have

u 0,N (ξ) = ǫ N N k=0 2 2k/b η k 1 2 [ w(ξ + 3 2 2 k e 1 ) + w(ξ -3 2 2 k e 1 )].
Hence, using that the support of w is contained in {|ξ| ≤ 1 4 },

u N (t, ξ) ≥ ǫ N N k=0 e -t 2 2k+2 2 2k/b η k 1 2 w(ξ + 3 2 2 k e 1 ) =A k (t,ξ) + e -t 2 2k+2 2 2k/b η k 1 2 w(ξ -3 2 2 k e 1 ) =B k (t,ξ)
.

This implies that

( u N ) * b (t, ξ) ≥ ǫ b N N k 1 =0 • • • N k b =0 (A k 1 + B k 1 ) * • • • * (A k b + B k b )(t, ξ) . (4.8)
It is now convenient to distinguish two cases.

The case b even. We bound from below (4.8) retaining just a few terms of the above summation:

( u N ) * b (t, ξ) ≥ ǫ b N N k=0 (A k * B k ) * b/2 (t, ξ). But, w(• + 3 2 2 k e 1 ) * w(• -3 2 2 k e 1 ) = w * w, hence, ( u N ) * b (t, ξ) ≥ ǫ b N N k=0 e -b t 2 2k+2 2 2k-b η b k ( w) * b (ξ).
Using that supp( w * b ) ⊂ B(0, 1), we deduce that

u N (t, ξ) ≥ t 0 e (s-t)|ξ| 2 u N (s, •) * b ds ≥ ǫ b N N k=0 2 2k-b η b k t 0 e (s-t)-b s 2 2k+2 ds ( w) * b (ξ) ≥ ǫ b N N k=0 2 -b η b k 4b 1 -e -t(b 2 2k+2 -1) e -t ( w) * b (ξ) ≥ ǫ b N N k=0 η b k 2 -b 4b
1 -e -t(4b-1) e -t ( w) * b (ξ).

Our choice of (η k ) and (ǫ k ) ensure that ǫ b N ( N k=0 η b k ) → +∞ as N → +∞. Now let us take t = δ/2 and N ≥ N 0 large enough in a such way that

ǫ b N N k=0 η b k 2 -b 4b 1 -e -δ(4b-1)/2 e -δ/2 ≥ b 2b/(b-1) 2 e δ/2 c 1/(b-1) δ w b 1 .
Hence Lemma 4.2 applies and implies that the lifetime of the solution of u t = ∆u + u b arising from the initial datum u N (δ/2, •) must blow up before the time δ/2. By the uniqueness result of Theorem 2.1, this implies that T * N < δ. The case b odd. In this case we can write b = 2m + 3, with m ∈ N. Going back to (4.8), we bound this expression from below in the following way: 1 -e -t(8(m+3)-1) e -t w * b (ξ).

( u N ) * b (t, ξ) ≥ ǫ b N N -1 k=0 (A k * B k ) * m * A k+1 * B k * B k (t, ξ).

By the invariance of convolution products under translation, supp(A

k * B k ) is contained in {|ξ| ≤ 1/b} and supp(A k+1 * B k * B k ) is contained in {|ξ| ≤ 3/(2b)}.
But ǫ b N N -1 k=0 η b-1 k η k+1 → +∞ as N → +∞ and therefore we can conclude taking t = δ/2 and applying Lemma 4.2, exactly as in the case b even.

Conclusions

The global-in-time solvability of the Cauchy problem for the nonlinear heat equation (2.4) in R n is usually obtained putting a smallness assumption on a suitable scale invariant norm of the initial data. However, in the present paper we proved that the scaleinvariant norm of the Besov space Ḃ-2/b nb(b-1)/2,q is not suitable for this purpose: in fact, arbitrarily small initial data in this space (or in any larger scale invariant space) can give rise to solutions that blow up in finite time. While our method provides a few quantitative estimates on the solution, it gives little information on the nature of the blowup. The issue of the type of blowup has been thoroughly investigated, e.g., in [START_REF] Collot | Dynamics near the ground state for the energy critical nonlinear heat equation in large dimensions[END_REF][START_REF] Herrero | Explosion de solutions d'équations paraboliques semilinéaires supercritiques[END_REF][START_REF] Matano | Classification of type I and type II behaviors for a supercritical nonlinear heat equation[END_REF].

Our result is sharp, in the sense that, for all s > -2/b, a smallness condition on the Ḃs p,q -norm of u 0 (with s -n/p = -2/b-2/(b(b-1)), to respect the scale invariance), which is slightly more stringent, does ensure that the solution is globally defined. On the other hand, the precise role of the third index q (that does not affect the scaling of the Besov norm) on this blowup issue is less clear: as the proof of Theorem 2.2 requires q > b, and breaks down when q = b, the following open problem naturally arises.

(4. 1 )

 1 f * b = f * . . . * f b times .

  k → +∞ as N → +∞. In the case b ≥ 3 is odd we will additionally needǫ b N N -1 k=0 η b-1k η k+1 → +∞ The choice (4.6) does satisfy these requirements. Observe that the Fourier transform of cos( 3 2 2 k x 1 )w(x) is 1 2 [ w(ξ + 3 2 2 k e 1 )+ w(ξ -3 2 2 k e 1 )], which is contained in the union of two balls centered at ± 3 2 2 k e 1 and radius 1/(2b) (and hence in a single dyadic annulus). Let us consider the homogeneous Littlewood-Paley decomposition, ∈Z ψ(2 -j ξ) = 1, for ξ = 0, obtained using a radial functionψ ∈ C ∞ 0 (R n ) which is supported in { 3 4 ≤ |ξ| ≤ 8 3 } and constant equal to 1 in { 5 4 ≤ |ξ| ≤ 7 4 }.If we denote by ∆ j f = ψ j * f the Littlewood-Paley dyadic blocks, then we see that∆ j cos( 3 2 2 k x 1 ) w(x) = 0, if k = j,and ∆ j cos( 3 2 2 j x 1 ) w(x) = cos( 3 2 2 j x 1 ) w(x).

e

  Hence, supp (A k * B k ) * m * A k+1 * B k * B k (t, •) ⊂ {|ξ| ≤ 1/2}. Hence, ( u N ) * b (t, ξ) ≥ ǫ b N N -1 k=0 e -(m+3) t 2 2k+3 2 2k-b η 2m+2 k η k+1 w * b (ξ).Arguing as before, we obtain,u N (t, ξ) ≥ t 0 e (s-t)|ξ| 2 ( u N ) * b (s, •) ds (s-t)-(m+3) s 2 2k+3 ds w * b (ξ)

  (Here b does not need to be an integer): Let n(b -1)/2 > 1 and u 0 a smooth and well decaying initial data as |x| → +∞; does the smallness assumption u 0 Ḃ-2/b nb(b-1)/2,b < < 1 imply that the solution of the Cauchy problem for u t = ∆u + |u| b-1 u in R n is global in time ?

Acknowledgement

The authors are grateful to Paul Acevedo and to the Referee for their useful comments and suggestions.

Fernando Cortez was supported by Escuela Politcnica Nacional, Proyecto PII-DM-08-2016.