Michael Mercier

Adrien Faure

Olivier Richard

Considering the Development Workflow to Achieve Reproducibility with Variation

The ability to reproduce an experiment is fundamental in computer science. Existing approaches focus on repeatability, but this is only the first step to reproducibility: Continuing a scientific work from a previous experiment requires to be able to modify it. This ability is called reproducibility with Variation.

In this contribution, we show that capturing the environment of execution is necessary but not sufficient; we also need the environment of development. The variation also implies that those environments are subject to evolution, so the whole software development lifecycle needs to be considered. To take into account these evolutions, software environments need to be clearly defined, reconstructible with variation, and easy to share. We propose to leverage functional package managers to achieve this goal.

Motivation

Reproducibility is a wide notion that needs to be specified. Feitelson [START_REF] Dror G Feitelson | From repeatability to reproducibility and corroboration[END_REF] has defined a taxonomy of the different way to reproduce scientific results. In this taxonomy, The first level of reproduction is the Repetition, i.e. do exactly the same experimental process to obtain the same results. The second level, Replication is similar but the experience's input is changed.

Currently, most of the reproducibility tools only support these two levels by capturing the software environment. Indeed, software environment is hard to reproduce, and without it, it is impossible to run the experiment. Also, experiment software environment tightly depends on the Operating System (OS) distribution it was built on. It is sometimes impossible to install it on an other distribution because of inter dependency issues. One approach to solve this problem is to snapshot the software environment into an image. But, even if an image of the experiment runtime environment is provided by the original author, continuing his work requires more than just repeating the experiment; to be able to corroborate someone's approach, we not only need to rerun experiments, but we also need to modify them: test new variations, add more parameters, and develop new features. This is the next level of reproducibility in the aforementioned taxonomy, called Variation.

Enabling scientists to reproduce an experiment with variation requires that the reproducer is able to rebuild experiment software with some modifications (even if the software is unmaintained and the tools necessary for building it, are long gone). It means that the reproducibility with variation can be achieved if the reproducer is able to reproduce not only the experiment "production" environment, but also the "development" environment which is necessary to modify the software. Moreover, when a variation of previous experiment produces new results, this experiment should also be reproducible.

In this context, we are proposing a new way of seeing reproducibility through the scientific software development lifecycle. Each step in this lifecycle requires a software environment. We define a software environment by a set of applications and libraries, with all their dependencies, and their configurations, required to achieve a step in a scientific workflow. In computer science, a scientific workflow contains a software development lifecycle that starts by setting up a development environment with build tools and dependencies. Then, this environment is used to build a production environment that will, in turn, be used to run the actual experiment. But, software development is an iterative process: one can produce different versions of the production environment, or even modify the development environment to update or to add tools. This process is in the middle of the scientific workflow and all the software environments produced, for development and production should be captured to enable reproducibility. The Figure 1, exhibits that the first two levels of reproducibility can be achieved with only one environment, but reproducibility with variation requires taking into account the whole development process.

Repetition

To contribute to the scientific workflow, it has to be reproducible by itself: Thus, internal (e.g. colleague, intern) and external (e.g. scientist from other laboratory, reviewer) contributors have the capability to reproduce this workflow.

So, to provide this capability the development environment should be defined entirely, and every changes should be tracked. This statement also holds for others software environments involved in the authors' workflow, like the ones use for input data generation, or output data analysis.

Reproducible Software Environments with Nix

Courtès et al. [START_REF] Courtès | Reproducible and User-Controlled Software Environments in HPC with Guix[END_REF] emphasize that functional package managers (FPM), like Guix and Nix, are good candidates to share complex and upgradable environment.

In the following of this paper, we will focus on Nix, but most of the assertions also hold for Guix. The FPM are applying the concept of mathematical function to software packaging. Each software building process is described through a function. The dependencies are also functions that are given as inputs. This function, or package definition, allows to precisely describe a package: where and how to gather the source code, which commit to use, the dependencies and their versions, and finally how to build the package. When a package is built, the dependency graph is resolved by a lazy evaluation of the function parameters, and all the necessary piece of software are also built. The result of the evaluation of a package definition is called a derivation. A derivation is concretely a set of files that contains the results of the building process of the software, which placed on the special place that contains all the derivations: the store. Finally installing a software is simply exposing a derivation from the store through symbolic links.

Nix packages are written in a functional Domain Spe-cific Language (DSL). This ensures that each build is pure, i.e. it only depends on its inputs, and the same inputs give the exact same package even on a different machine.

To implement the workflow described in Section 2 with Nix, the critical feature is the capacity to create a software environment without virtualization. This feature is used to create an isolated environment for the package building process. An environment can be seen as a set of derivations and relies on the fact that an FPM can infer all the dependencies of a derivation, and only expose these dependencies on the specified environment. Installing the environment will expose to the user the packages described in the environment. Archiving an environment will extract the whole dependency tree of the environment and create a self-contained archive. The resulting tarball, also called a closure, contains every binary and file necessary to run the packaged applications. Thus, the environment can be installed on another machine without any external download or building process.

To achieve each level of reproducibility with Nix, the first requirement is to create a package definition for each application and its dependencies. Thankfully, Nix is a very active community and more than 40000 applications are already packaged in the main repository called "nixpkgs". It is also possible to maintain a private set of packages that import dependencies from "nixpkgs".

The first level of reproducibility, the Repetition, can be achieved by providing to an external scientist the production closure, that can be installed to repeat the experiment. The environment could also contain all scripts to deploy and run the experiment.

The second level of reproducibility, the Replication, that consists of replaying an experiment while changing its input, have the same requirement as the repetition: Only the production environment is needed.

The third level of reproducibility, the Variation, is where using Nix is the most advantageous. Nix provides the interesting feature, called the "nix-shell", permitting to enter the package build environment. Hence, packaging a software with Nix have the side effect to provide also the build environment for the users. Nix capacity to define software environment and software package in a unified way gives the scientist the ability to share a reproducible production environment, and the associated development environment, with a single definition.

Related works

The Popper method, proposed by Jimenez et al. [START_REF] Jimenez | The popper convention: Making reproducible systems evaluation practical[END_REF], describes a structural framework for dependencies and artifacts. They identified a generic workflow describing an experimental methodology, from source code to the final manuscript of a contribution. Our contribution is compatible with this approach, the popper method proposes a structural organization of the experiment, whereas our approach proposes to implements a part of this workflow using Nix.

Repeatability has been the focus of previous works on reproducibility. The platform presented in [START_REF] Ricci | Apt: A platform for repeatable research in computer science[END_REF], has the ability to instantiate an experiment environment in their infrastructure from a previously captured environment. The approach is interesting as it provides a way for a scientist to repeat experiments that requires specific hardware. Our approaches could be complementary to cover both hardware and software to provide a higher level of reproducibility. Boettiger et al. [START_REF] Boettiger | An introduction to Docker for reproducible research[END_REF] survey how to use docker to do reproducibility, and also introduce the development environment. From our implementation with Nix, the docker approach shows similarities. However, Nix closure is more adapted than the Docker images for application packaging because Docker provides an inappropriate level of abstraction: Docker is about constructing and configuring a complete OS, instead of declaring application dependencies.

Constructing reproducible experiment and workflow with an FPM has already been explored. In [START_REF] Wurmus | Reproducible genomics analysis pipelines with GNU Guix[END_REF], they build a toolset upon an FPM (GUIX), to facilitate the usage of bioinformatics common pipelines. They argue that using an FPM is a good foundation for reproducible computational experiment workflow.

The Blue Brain project [START_REF] Devresse | Nix based fully automated workflows and ecosystem to guarantee scientific result reproducibility across software environments and systems[END_REF], is a big project, with a complex software stack, that aims to build a mammalian brain with a computer. In addition to a structured development workflow (using git, agile methodologies, code reviewing), they decided to package their workflow with Nix. They identified nine properties that are facilitated with Nix, from reproducibility to deployment and cross compilation. Their approach is based on internal needs and specific use case, whereas our contribution focuses on the role of the development lifecycle in the reproducibility, with Nix as a possible implementation.

Discussion

The proposed workflow is not considering input and output data. One approach is to package the experiment data inside the production environment. It is a viable solution for small amount of data, but most of the experiments have managed data separately with other tools.

The data related to the experiment is not only input and output data, the software environments that are used to develop and run the experiment also have external inputs, like the source code of the experiment, the dependencies, the build tools binaries, and the configuration files, i.e. all the other artefacts necessary to build the environments. With Nix, it is easy to extract these artefacts for the experiment production environment. Nix is capable to export all the dependencies of an environment in a closure that can be imported on any machine where Nix is installed. This feature of Nix is very hard to achieve with other kind of tools because of the lack of clear dependency definition. However, even if it is straightforward for the production environment, how to extract this closure for the development environment is unclear for now.

Capturing and archiving the environments closures is necessary for the Variation reproducibility, because the lifespan of Internet links is only a few months [START_REF] Lawrence | Persistence of web references in scientific research[END_REF]. Even if Nix is capable to rebuild everything from source, the source code repository can be unavailable, breaking the environment reproducibility. The problem that emerges is that closures also have to be safely archived, versioned, and accessible for a long time period. A mid-term solution would be to store those closures (or even replace the closure content) using trusted centralized archives like Internet Archive 1 and Software Heritage 2 .

In the case where the experiment results depend on the OS Kernel, (e.g. performance evaluation) this approach is not sufficient. Indeed, Nix packaging handles the whole application software but not the OS Kernel. So, the proposed workflow needs to be supplemented with a building process definition of the entire OS -and not just the application layerto be able to reconstruct a complete OS image with variation. A Linux distribution based on Nix, called NixOS, is a good candidate.

When a specific hardware is necessary to achieve reproducibility, an additional layer of control is needed. Tesbeds like Grid'5000, Chameleon, and Emulab, are giving this level of control with the capability to create and deploy OS image on the fly on different hardware.

Conclusion

The reproducibility with variation is the next level of reproducibility that the Computer Science community should aim at. The Variation requires to take into account the software development workflow, including the capability to modify and rebuild environments. The use of functional package managers is a promising approach. This kind of tool permits to achieve this last mile to the reproducibility with variation, with a unified way to describe environments and packages, and a simple method to backup and to restore them.

Figure 1 :

 1 Figure 1: The different level of reproducibility in regard to the development lifecycle: Variation requires to enclose the development environment and to provide a way to modify it while keeping reproducibility.


• Parametrization and tests based on higher-order statistics (curvature, energy flux through scales, bispectrum, …)

• Mimic barotropization

• Girsanov theorem for MLE and Bayesian estimations with satellite images

• Learning on SWOT data

• Filtering / smoothing with (reduced) models under location uncertainty 14 V. Resseguier -valentin.resseguier@gmail.com

  Θ @ t ⇥ + w ? • r⇥ + Ḃ • r⇥ = r • t ⇥ + w ? • r⇥ + Ḃ • r⇥ = r • Θ @ t ⇥ + w ? • r⇥ + Ḃ • r⇥ = r • t ⇥ + w ? • r⇥ + Ḃ • r⇥ = r • t ⇥ + w ? • r⇥ + Ḃ • r⇥ = r • t ⇥ + w ? • r⇥ + Ḃ • r⇥ = r • t ⇥ + w ? • r⇥ + Ḃ • r⇥ = r • t ⇥ + w ? • r⇥ + Ḃ • r⇥ = r • -valentin.resseguier@gmail.comRelated works and outlooks • Bifurcations (SQG) and attractor (Lorenz 63) exploration • Stabilization / destabilization of Reduced Order Model • Comparisons with data-driven and Stochastic Lie Derivative approaches (Holm and coauthors)


C.-Y Hung

P Aussillous

H Capart

Granular surface avalanching induced by drainage from a narrow silo

Using theory and experiments, we investigate granular surface avalanching due to material outflow from a narrow silo. The assumed silo geometry is a deep rectangular box, of moderate spanwise width and small gap thickness between smooth front and back walls. A small orifice deep below the free surface lets grains drain out at a constant rate. The resulting granular flows can therefore be assumed quasi-two dimensional and quasi-steady over most of the surface descent history. To model these flows, we couple a kinematic model of deep granular flow with a dynamic model of shallow surface avalanching. We then compare the calculated flow fields with detailed particle tracking measurements, letting the silo descend relative to the high-speed camera to increase spatial resolution. The results show that the avalanching surface shape and near-surface flow are controlled by the spanwise gradient in subsidence velocity, and how this gradient is in turn controlled by the height above orifice and the gap thickness. Whereas the deep flow pattern is rateindependent, shallow avalanching is paced by the granular rheology.

Introduction

In various environmental and industrial contexts, the removal of material beneath the surface of a granular substrate can cause the surface above the withdrawal zone to subside and steepen. Examples include mining-induced ground subsidence [START_REF] Vivanco | The effect of rock decompaction on the interaction of movement zones in underground mining[END_REF], and material drainage from granular hoppers or pebble bed reactors [START_REF] Cleary | DEM modelling of industrial granular flows: 3D case studies and the effect of particle shape on hopper discharge[END_REF][START_REF] Rycroft | Analysis of granular flow in a pebble-bed nuclear reactor[END_REF]. If steepening causes the surface inclination to exceed the angle of repose, avalanching will occur. One common example occurs in the upper chamber of an hour glass, where a funnel-shaped avalanching surface forms due to drainage through the underlying narrow throat. In the present work, we examine these phenomena in a simplified geometry: rectangular silos holding grains between two closely spaced smooth walls, making the flows quasi-two-dimensional. Such flows have earlier been investigated using laboratory experiments [START_REF] Gray | Pattern formation in granular avalanches[END_REF][START_REF] Samadani | Shocks in sand flowing in a silo[END_REF][START_REF] Benyamine | Discharge flow of a bidisperse granular media from a silo[END_REF][START_REF] Maiti | Experiments on eccentric granular discharge from a quasi-two-dimensional silo[END_REF], continuum models [START_REF] Staron | Continuum simulation of the discharge of the granular silo. A validation test for the µ(I) visco-plastic flow law[END_REF][START_REF] Dunatunga | Continuum modelling and simulation of granular flows through their many phases[END_REF][START_REF] Daviet | Nonsmooth simulation of dense granular flows with pressure-dependent yield stress[END_REF], and discrete element simulations [START_REF] Cleary | DEM modelling of industrial granular flows: 3D case studies and the effect of particle shape on hopper discharge[END_REF][START_REF] Staron | Continuum simulation of the discharge of the granular silo. A validation test for the µ(I) visco-plastic flow law[END_REF]Zhou et al. 2015).

To simplify the problem further, we consider silos of small width to depth ratios, restricting the span of the avalanching layer (Fig. 1). Experimental observations then suggest that the following simplifications can be made. When drainage at depth through a bottom or side orifice occurs at a constant rate, first, the resulting granular surface flow can be considered quasi-steady. After a starting transient, during which the granular free surface deforms from its initial shape to its avalanching shape, the surface flow becomes nearly steady in a frame of reference descending at constant speed with the free surface. Avalanching thus takes the form of a travelling wave. Because it is driven by a deep subsidence flow that varies spatially with height, the avalanching flow is not perfectly steady. Except for low submergence above the orifice, however, the spatial evolution of this deep flow is gradual, hence the descending avalanching flow experiences forcing changes that are slow with respect to its adaptation time. At any given time over most of its descent history, therefore, the avalanching flow can be approximated as a steadily flowing layer that continually exchanges grains with the substrate through entrainment and detrainment. This makes avalanching flows in narrow silos analogous to other steady granular layer flows subject to erosion and deposition, like flows in rotating drums [START_REF] Gray | Granular flow in partially filled slowly rotating drums[END_REF][START_REF] Hung | Granular flow regimes in rotating drums from depth-integrated theory[END_REF]) and bounded heaps [START_REF] Khakhar | Surface flow of granular materials: model and experiments in heap formation[END_REF][START_REF] Fan | Kinematics of monodisperse and bidisperse granular flows in quasi-two-dimensional bounded heaps[END_REF]).

In the present paper, we investigate quasi-steady, quasi-two-dimensional silo avalanching flows using a combination of theory and experiment. Using either continuum or discrete models, one could attempt to model together both the deep flow and surface avalanching regions of the phenomenon (see e.g. [START_REF] Staron | Continuum simulation of the discharge of the granular silo. A validation test for the µ(I) visco-plastic flow law[END_REF]). This is however made difficult by the highly different scales that govern the two regions. Whereas the deep flow is steady, slow, and gradually varying in space, the avalanching flow is confined to a thin surface layer featuring faster velocities and steeper velocity gradients. Analogous to the boundary layer approach for viscous and turbulent flows, therefore, we adopt a modeling strategy that couples different equations for different zones of the flow. For the deep flow, we adopt the kinematic model proposed by [START_REF] Nedderman | A kinematic model for the flow of granular materials[END_REF], which reduces the problem to the solution of a diffusion equation.

To describe surface avalanching, on the other hand, we adopt the model proposed by [START_REF] Capart | Depth-integrated equations for entraining granular flows in narrow channels[END_REF] and recently applied by [START_REF] Hung | Granular flow regimes in rotating drums from depth-integrated theory[END_REF] to granular flows in rotating drums. Instead of solving local partial differential equations [START_REF] Jop | Initiation of granular surface flows in a narrow channel[END_REF][START_REF] Lusso | A free interface model for static/flowing dynamics in thin-layer flows of granular materials with yield: simple shear simulations and comparison with experiments[END_REF][START_REF] Fernández-Nieto | 2D granular flows with the µ(I) rheology and side walls friction: A well-balanced multilayer discretization[END_REF], the model uses depth-integrated equations for the balance of mass, momentum and kinetic energy of the shallow avalanching layer. This model applies to dry, dense, two-dimensional granular flows in narrow channels, over deep erodible deposits, when the flow and bed are composed of the same monodisperse grains. The model further assumes that the flow is subject to friction along the walls and governed by a linearized form of the viscoplastic µ(I) rheology [START_REF] Jop | Crucial role of sidewalls in granular surface flows: consequences for the rheology[END_REF][START_REF] Da Cruz | Rheophysics of dense granular materials: Discrete simulation of plane shear flows[END_REF], such that the base of the flow coincides with the yield locus. As a result, it cannot describe flows over non-erodible beds [START_REF] Parez | Unsteady granular flows down an inclined plane[END_REF][START_REF] Sarno | Experimental investigation on the effects of the fixed boundaries in channelized dry granular flows[END_REF] or the transitions between flow and rest observed for shallow granular layers in wide channels [START_REF] Pouliquen | Friction law for dense granular flows: application to the motion of a mass down a rough inclined plane[END_REF][START_REF] Edwards | Erosion-deposition waves in shallow granular free-surface flows[END_REF].

For the experiments, we adopt a configuration aimed at maximizing the resolution at which the deep and surface flows can be recorded by a high speed camera. A halfsilo geometry is adopted, and the silo is mounted on a motorized traverse to make it rise steadily relative to the camera, keeping the avalanching flow stationary relative to the observation window. Particle tracking velocimetry is then applied to the footage. Averaging of the deep flow and avalanching layer velocity fields is performed in the separate frames of reference in which they can be considered steady.

The paper is structured as follows. In section 2, we describe the theories adopted to model deep flow and surface avalanching, and how they are coupled together. In section 3, we present the experimental set-up, conditions, and imaging analysis. In section 4, we compare and discuss results, before drawing conclusions in section 5.

Theory

The geometrical conditions considered in the present work are defined in Fig. 1. We assume silos of rectangular shape, with smooth front and back walls separated by a small gap thickness B equal to some multiple (of the order of 10) of the grain diameter D. We therefore consider three-dimensional grain configurations, not monolayers. Because the gap is narrow and the walls smooth, however, we can assume that the flow field will be effectively two-dimensional. For simplicity, we consider either the symmetric geometry of Fig. 1(a), of spanwise width 2L on either side of a centered botttom orifice, or the asymmetric geometry of Fig. 1(b) having width L and a side orifice at X = 0. For the near-field around the orifice, differences between bottom and side outlets are expected, leading to different variants of the Beverloo relation expressing the rate of outflow in terms of grain diameter, gap thickness and orifice opening [START_REF] Zhou | Experiments on[END_REF]. Regardless of the case, however, a constant flow rate is obtained when the orifice opening is large enough to prevent jamming, and when the silo free surface is sufficiently high above the orifice. Both conditions will be assumed, hence we consider the outflow rate known and constant. Throughout the silo, only a single type of grains is considered.

Upon starting the silo discharge, grains will descend and converge towards the orifice, drawing down the free surface. Avalanching will then occur once uneven subsidence has steepened the granular free surface beyond its angle of repose. After this inital transient, we consider separately the shallow avalanching layer and the deep flow region comprised between this layer and the orifice near-field. Since the outflow rate is constant, flow in the deep region will be assumed steady in the silo frame of reference defined by axes (X, Z). For sufficiently small silo widths L, the span of the avalanching flow will be bounded by the side walls. We therefore expect the avalanching layer to attain and maintain a quasisteady flow state in a frame of reference descending with the surface at the constant mean rate Ḣ, where H(t) is the time-varying mean height of the free surface above the orifice. After formulating the theory, we will be able to define more clearly the range of conditions for which this assumption is reasonable. In the following two sub-sections, we describe the distinct models adopted for the deep flow and avalanching layers, respectively. In both zones, variations in solid fraction are neglected, hence incompressible flow is assumed.

In earlier work [START_REF] Da Cruz | Rheophysics of dense granular materials: Discrete simulation of plane shear flows[END_REF][START_REF] Jop | Crucial role of sidewalls in granular surface flows: consequences for the rheology[END_REF], incompressibility was found to be a good approximation for dense, slow flows characterized by low inertia numbers. These are the conditions encountered in the experiments described below. Our objective is to predict the gradually evolving granular velocity field (U, W )(X, Z, t) (averaged over the gap width) and shape of the free surface Z(X, t).

Deep flow kinematics

To describe the steady deep flow (U 0 , W 0 )(X, Z), we adopt the kinematic model of [START_REF] Nedderman | A kinematic model for the flow of granular materials[END_REF]. In this model, vertical shear is assumed to cause grains to drift laterally at velocity

U 0 = -K ∂W 0 ∂X , (2.1)
where U 0 (X, Z), W 0 (X, Z) are respectively the horizontal and vertical components of the mean granular velocity in the deep flow region, and K is a coefficient with dimensions of length called the kinematic constant by [START_REF] Nedderman | A kinematic model for the flow of granular materials[END_REF]. Substitution of (2.1) into the continuity equation

∂U 0 ∂X + ∂W 0 ∂Z = 0 , (2.2)
then yields the diffusion equation with diffusivity

K ∂W 0 ∂Z -K ∂ 2 W 0 ∂X 2 = 0 , (2.3)
where diffusion proceeds upwards along the timelike Z axis to even out horizontal variations in W 0 . The micro-mechanical basis of (2.1) and (2.3) remains to be clarified, but a simple physical explanation is as follows [START_REF] Litwinisyn | On sand displacements in neighborhood of an immobile vertical plane[END_REF][START_REF] Caram | Random-walk approach to granular flows[END_REF]. At the discrete level, the downward motion of individual grains is associated with the upward migration of void spaces, or holes, into which individual grains fall. At each fall, the hole can move either right or left with equal probability. The resulting upward random walk of individual void spaces yields a diffusion process in the continuum limit. The opposite granular flux W 0 is then governed by the diffusion equation (2.3). Although the equation lacks a rigorous mechanical basis, it describes a physically plausible deep flow velocity field that can be used to drive surface avalanching. Upon calibrating the diffusivity K, equation (2.3) has been found to yield good agreement with silo experiments [START_REF] Tüzün | The flow of granular materials-II. velocity distributions in slow flow[END_REF][START_REF] Choi | Velocity profile of granular flows inside silos and hoppers[END_REF] and discrete element simulations [START_REF] Rycroft | Analysis of granular flow in a pebble-bed nuclear reactor[END_REF][START_REF] Balevičius | Analysis and DEM simulation of granular material flow patterns in hopper models of different shapes[END_REF].

For the silos we consider, the diffusion equation must be solved subject to the lateral boundary conditions ∂W 0 ∂X = 0 , at X = -L and X = 0 , (2.4) imposing no flux across the left wall and symmetry plane (for the symmetric silo), or across the left and right walls (for the asymmetric silo). At the outlet level Z = 0, on the other hand, a bottom boundary condition can be written

W 0 = -2Qδ(X) , (2.5)
where Q is either half the rate of outflow (for the symmetric case), or the full outflow (for the asymmetric case), divided by the gap width B (see Fig. 1). Disregarding the details of the flow near the orifice, a delta function is used to represent the outflow as a point sink. Since the diffusion equation features a first derivative in Z and a second derivative in X (see e.g. [START_REF] Kevorkian | Partial Differential Equations[END_REF]), no condition is required or can be applied along the granular free surface, nor can additional conditions restraining slip be applied along the side walls. The deep flow is driven by drainage from the orifice, and is not affected by avalanching at the free surface. In addition to the velocity field, we will also be interested in the streamfunction Ψ 0 (X, Z) defined by

∂Ψ 0 ∂Z = U 0 , ∂Ψ 0 ∂X = -W 0 , (2.6)
the contours of which give the granular streamlines. For definiteness, we associate Ψ 0 = 0 with the vertical streamline through the origin X = 0. The above mathematical problem can be made dimensionless by defining variables X = X/L, Ẑ = ZK/L 2 , Ŵ0 = W 0 L/Q, Û0 = U 0 L 2 /(KQ) and Ψ0 = Ψ 0 /Q. As detailed in [START_REF] Kevorkian | Partial Differential Equations[END_REF], the solution can be obtained using either Fourier series, advisable for Ẑ moderate to large, or mirrored Green functions, advisable for Ẑ small to moderate. Since the latter case better fits the conditions of our experiments, we will calculate the solution using the Green function

G( X, Ẑ) = 1 (π Ẑ) 1/2 exp - X2 4 Ẑ , (2.7)
which solves (2.3) subject to bottom boundary condition (2.5) on a laterally unbounded domain. To enforce lateral boundary conditions (2.4), mirrored green functions must be added to obtain

Ŵ0 ( X, Ẑ) = -G( X, Ẑ) - n k=1 G( X -2k, Ẑ) + G( X + 2k, Ẑ) , (2.8)
which is an infinite sum truncated at n terms that converges to the solution for n large [START_REF] Kevorkian | Partial Differential Equations[END_REF]). The solution for the streamfunction, likewise, can be obtained by integration of (2.8) along X as

Ψ0 ( X, Ẑ) = Φ( X, Ẑ) + n k=1 Φ( X -2k, Ẑ) + Φ( X + 2k, Ẑ) , (2.9)
where

Φ( X, Ẑ) = erf X 2 Ẑ1/2 .
(2.10)

The resulting dimensionless solutions are illustrated in Fig. 2. Figure 2(a) shows spanwise profiles of horizontal velocity at different heights above the orifice. Figure 2(b) shows the corresponding profiles of vertical velocity. Finally Fig. 2(c) shows the streamlines Ψ = constant for equally spaced values of the streamfunction. Near the orifice, flow is apparent only in the central region, but there is no sharp boundary between flowing and stationary regions. The diffusion solutions feature small non-zero motions even away from the central zone. Using the kinematic model, similar deep flow fields in silos have earlier been calculated by [START_REF] Nedderman | A kinematic model for the flow of granular materials[END_REF], using Fourier series, and by [START_REF] Choi | Velocity profile of granular flows inside silos and hoppers[END_REF] using numerical computations.

Shallow layer governing equations

To model surface avalanching (Fig. 3), we consider a distinct coordinate system (x, z) tilted at the angle of repose α, with origin at (X, Z) = (-L/2, H(t)) and descending with the granular free surface at speed Ḣ = -Q/L. The x-axis of the coordinate system therefore coincides with the inclined profile (Fig. 3(a))

Z(X, t) = H(t) -(X + L/2) tan α .

(2.11)

Throughout this section, we assume the restricted domain -L X 0, with mirror symmetry implied in case of a symmetric silo. In the inclined coordinate system, we define z(x, t) and z(x, t) to be the free surface and basal boundary of the avalanching layer. The depth of the avalanching layer is then given by h(x, t) = z -z. We assume shallow flow, h L, and small relief z L relative to profile z = 0 inclined at the angle of repose. Consistent with these approximations, flow acceleration in the z-direction will be neglected relative to acceleration in the downslope x-direction. However both velocity components, u(x, z, t), in the x-direction and w(x, z, t) in the z-direction, will intervene in the flow kinematics. Inspired from [START_REF] Liggett | Fluid Mechanics[END_REF], and following [START_REF] Capart | Depth-integrated equations for entraining granular flows in narrow channels[END_REF] and [START_REF] Hung | Granular flow regimes in rotating drums from depth-integrated theory[END_REF], we use different lines above and below variables to denote values sampled at different locations: • denotes a variable sampled along the basal boundary, • a variable sampled along the free surface, and • a variable sampled along the angle of repose profile (2.11). The overbar • is reserved for quantities averaged over the depth h.

Following [START_REF] Capart | Depth-integrated equations for entraining granular flows in narrow channels[END_REF], we start from local equations and integrate them from z to z to obtain depth-integrated equations governing the avalanching layer. As for the deep flow, we assume incompressibility of the granular medium, hence the local continuity equation

∂u ∂x + ∂w ∂z = 0 .
(2.12)

Variations of u and w in the transverse y-direction are neglected, hence u and w coincide with the corresponding velocity components averaged over the gap thickness. Local balance of momentum in the x-direction, on the other hand, yields the width-averaged equation of motion [START_REF] Jop | Initiation of granular surface flows in a narrow channel[END_REF])

ρ ∂u ∂t + u ∂u ∂x + w ∂u ∂z = ρg - ∂p ∂x + ∂τ ∂z - 2τ W B , (2.13)
where ρ = c S ρ S is the mean density of the bulk phase, and p = ρg ⊥ ( z -z) is the granular pressure, assumed lithostatic. Following [START_REF] Jop | Crucial role of sidewalls in granular surface flows: consequences for the rheology[END_REF], the normal stresses are assumed to reduce to an isotropic pressure. For shallow flow [START_REF] Liggett | Fluid Mechanics[END_REF], granular accelerations in the z direction normal to the slope can be neglected, hence a lithostatic pressure balancing the weight of the grains above. In the same equation, g = g sin α and g ⊥ = g cos α are the parallel and perpendicular components of the gravitational acceleration, τ W = µ W p is the shear stress along the sidewalls (with µ W the wall friction coefficient), and τ is the internal shear stress. For the latter, we assume the viscoplastic, dense granular flow rheology [START_REF] Da Cruz | Rheophysics of dense granular materials: Discrete simulation of plane shear flows[END_REF])

τ = µ 0 p + χD(ρp) 1/2 γ , (2.14)
expressed as the sum of two components. The first, rate-independent, is a plastic yield stress proportional to p, with coefficient µ 0 = tan α. The second, rate-dependent, is a viscous stress that varies linearly with the shear rate γ = ∂u/∂z, with effective viscosity χD(ρp) 1/2 dependent on grain diameter D, granular pressure p, and dimensionless rheological coefficient χ. This is a linearized version of the µ(I) rheology proposed by [START_REF] Jop | Crucial role of sidewalls in granular surface flows: consequences for the rheology[END_REF] and applied earlier to liquid-granular flows by [START_REF] Berzi | A theoretical analysis of free-surface flows of saturated granular-liquid mixtures[END_REF]. The non-linear µ(I) rheology has earlier been applied to shallow granular flows over rigid beds by [START_REF] Gray | A depth-averaged µ(I)-rheology for shallow granular free-surface flows[END_REF]. The two rheological coefficient µ 0 and χ will be determined in the next section by linearizing the non-linear µ(I) rheology calibrated earlier by [START_REF] Jop | Crucial role of sidewalls in granular surface flows: consequences for the rheology[END_REF] for the same granular material. Note that the constitutive relations for wall friction and internal shear stress are written assuming u > u W , and γ > 0. For steady uniform flow over a loose deposit of velocity ( u, w), these equations are satisfied by the velocity profile (Fig. 3(b))

u( η) = u + (u -u) 7 3 -35 6 η 3/2 + 7 2 η 5/2 , (2.15)
where η = ( z -z)/h and u is the depth-averaged velocity (for a derivation, see [START_REF] Berzi | A theoretical analysis of free-surface flows of saturated granular-liquid mixtures[END_REF]). As discussed in [START_REF] Capart | Depth-integrated equations for entraining granular flows in narrow channels[END_REF], experiments exhibit some deviations from the profile (2.15) even for equilibrium flows. For unsteady, non-uniform flows, moreover, the velocity profile is expected to deviate from its equilibrium shape. As long as such deviations are small, however, integrals over depth based on the profile (2.15) should furnish a good approximation, and provide a simple way to close depth-integrated equations. Because the corresponding shear rate at the base γ = ∂u/∂z( z) is equal to zero, note that this profile is appropriate only for flows over erodible beds, when the basal interface can be assumed to coincide with the yield locus. For flows over rigid beds, say over rock, a finite shear rate at the base would be expected, and the base elevation prescribed, instead of the present case with zero basal shear rate and a basal elevation set by the flow itself. Likewise, the description would not apply for granular flows over deposits composed of a different type of grains. When assumptions are met, the velocity profile u(x, z, t) at a given coordinate x and time t can be calculated using (2.15) from z(x, t), h(x, t), u(x, t) and u(x, t). The corresponding streamfunction, likewise, can be calculated from

ψ(η) = ψ -h uη -h(u -u)( 7 3 η -7 3 η5/2 + η7/2 ) , (2.16)
where ψ satisfies ∂ψ/∂z = u, ∂ψ/∂x = -w, and such that

ψ -ψ = hu .
(2.17)

To obtain unsteady, non-uniform governing equations for the avalanching layer, we first integrate (2.12) and apply the Leibniz integral rule to get

z z ∂u ∂x + ∂w ∂z dz = ∂(hu) ∂x -u ∂ z ∂x + u ∂ z ∂x + w -w = 0 , (2.18)
where q = hu = z z udz. This can be further worked out using the kinematic boundary condition along the free surface,

∂ z ∂t + u ∂ z ∂x = w , (2.19)
and the definition h = z -z to obtain the depth-integrated continuity equation

∂h ∂t + ∂(hu) ∂x = e = w -u ∂ z ∂x - ∂ z ∂t , (2.20)
where e(x, t) is the erosion or entrainment rate, or rate of granular volume transfer across the basal interface z(x, t) and into the avalanching layer. This rate is defined as the deviation from volume conservation of the layer, expressed by the left-hand-side. Likewise, the local equation of motion (2.13) can be integrated over depth to obtain the depth-integrated momentum equation

∂ ∂t (hu) + ∂ ∂x h u 2 = e u + g h -g ⊥ h ∂ z ∂x - 1 ρ τ - 2 ρB hτ W .
(2.21)

Assuming u(x, t), w(x, t) to be otherwise given, we must solve for three evolving profiles h(x, t), z(x, t) and u(x, t). We therefore supplement (2.20) and (2.21) with an additional equation expressing the balance of depth-integrated kinetic energy [START_REF] Capart | Depth-integrated equations for entraining granular flows in narrow channels[END_REF][START_REF] Hung | Granular flow regimes in rotating drums from depth-integrated theory[END_REF]. This is obtained by integrating over depth the product of (2.13) with velocity u(x, z, t), yielding

∂ ∂t h u 2 2 + ∂ ∂x h u 3 2 = e u 2 2 + g hu -g ⊥ hu ∂ z ∂x -τ u - 1 ρ h τ γ - 2 ρB h τ W u . (2.22)
In this equation, the next-to-last and last terms represent dissipation by granular contacts inside the avalanching layer and dissipation by friction along the walls, respectively. Equations (2.20)-(2.22) represent a generalization of the equations derived in [START_REF] Capart | Depth-integrated equations for entraining granular flows in narrow channels[END_REF] for the case of stationary deposits ( u = w = 0). We now specialize them to the problem of avalanching flows over loose deposits undergoing uneven subsidence due to silo discharge. Ideally, we would wish to match the velocity at the base of the avalanching layer with the deep flow velocity along the same curve, i.e. let ( u, w) = ( u 0 , w 0 ). The precise shape of the basal interface z(x, z, t), however, is not known in advance. As in linearized groundwater seepage problems (see e.g. [START_REF] Ni | Groundwater drainage and recharge by networks of irregular channels[END_REF]), therefore, we choose instead to apply kinematic matching to the approximate basal interface given by Z(X, t). This corresponds to the linear profile (2.11), inclined at the angle of repose. We therefore approximate

u ≈ u 0 = U 0 cos α -(W 0 -Ḣ) sin α , (2.23) w ≈ w 0 = U 0 sin α + (W 0 -Ḣ) cos α , (2.24)
where velocities in the upright, stationary silo frame of reference are converted to the tilted frame of reference descending at constant speed Ḣ. Likewise we approximate

w + u ∂ z ∂x ≈ w 0 + u 0 ∂z ∂x = w 0 . (2.25)
In terms of streamfunction, equivalently, we let

ψ ≈ ψ 0 = Ψ 0 + ḢX , (2.26)
where Ψ 0 is the streamfunction (2.9) sampled along (X, Z)(x, t) = (-L/2+x cos α, H(t)x sin α). The solutions presented below are obtained by applying the matching condition at this known, but approximate interface position Z(X, t). We could instead proceed by iterations, applying the matching condition again at each step to the interface location calculated from the previous step. Because the deep flow varies gradually with height above orifice, however, the error of the approximate solution with respect to the converged iterative solution is small, hence we prefer the simpler approach. With this approximation, the depth-integrated continuity equation becomes

∂ z ∂t + ∂(hu) ∂x = w 0 . (2.27)
For quasi-steady flow, the dominant balance is between the second and third term, hence w 0 ∼ hu/L. For shallow flow the ratio h/L = is small, thus the deep flow velocity components u 0 ∼ w 0 ∼ u are one order of magnitude smaller than the depth-averaged avalanche velocity u.

= -g ⊥ h ∂ z ∂x - µ W B g ⊥ h 2 , (2.28) ∂ ∂t 77 96 hu 2 + ∂ ∂x κhu 3 = -g ⊥ hu ∂ z ∂x -35 9 χDg 1/2 ⊥ u 2 h 1/2 -5 9 µ W B g ⊥ h 2 u , (2.29)
where κ = 342853 233376 ≈ 1.469, and where the depth-integrated terms u 2 and u 3 were integrated using profile (2.15) with u = 0. The last two terms of (2.29), representing energy dissipation, are obtained by integration of the corresponding terms τ γ and τ W u in (2.22), using profile (2.15) and the rheology (2.14). The next-to-last term of (2.29) is thus controlled by the rheology coefficient χ that determines the effective viscosity. After simplification, equations (2.28) and (2.29) are identical to the momentum and kinetic energy equations derived earlier in [START_REF] Capart | Depth-integrated equations for entraining granular flows in narrow channels[END_REF]. Together with (2.27) they form a closed set of governing equations for profiles z(x, t), h(x, t), and u(x, t), driven by the deep flow velocities U 0 (x, t) and W 0 (x, t). These in turn can be estimated using the kinematic model described in the previous section. While these unsteady equations could in principle be solved numerically to describe the transient response, here we consider only the quasi-steady response obtained by neglecting the time derivatives in (2.27)-(2.29). The semi-analytical approach used to construct solutions for this case is described in the next sub-section.

Quasi-steady avalanching solutions

For quasi-steady flow, the continuity equation (2.27) reduces to

∂q ∂x = w 0 = - ∂ψ 0 ∂x , (2.30)
where q = hu is the depth-integrated granular discharge of the avalanching layer.

Integration with boundary conditions q = 0, ψ 0 = 0 at x = L/(2 cos α) then yields

q(x, t) = -ψ 0 (x, t) , (2.31)
such that q = 0 is also automatically satisfied at x = -L/(2 cos α). As expected, it follows from (2.17) and (2.31) that ψ = 0, i.e. the quasi-steady free surface coincides with the zero streamline. For quasi-steady flow, the downslope distribution of the avalanching discharge q(x, t) is therefore fully known, including its maximum value q max (t) at any given time t during the descent. An example is illustrated in Fig. 4(a). Hereafter, we drop the explicit time-dependence, but note that the evolving height H(t) of the avalanching layer during its descent will affect its discharge profile q(x).

To solve the momentum and kinetic energy equations (2.28) and (2.29) subject to the known discharge profile q(x), we substitute u = q/h and introduce the dimensionless variables x = x/(L/ cos α), q = q/Q, ĥ = h/h a , Ŝ = (-∂ z/∂x)/S a , where where Q is the dimensionless silo discharge

h a = χDQ g 1/2 ⊥ µ W /B 2/
Q = Q (µ W /B) 1/8 g 1/2 ⊥ (χD) 3/4 (L/ cos α) 7/8 . (2.35)
An equivalent dimensionless number, found to control avalanching behavior in rotating drums, was given the name entrainment number by [START_REF] Hung | Granular flow regimes in rotating drums from depth-integrated theory[END_REF]. It can be seen from (2.33) and (2.34) that this number controls the strength of the left-hand-side flux divergence terms associated with convective inertia. For low Q, the avalanching layer is locally in equilibrium and its behavior is dominated by the dense granular flow rheology.

In the inertialess limit Q → 0, in particular, the flux divergence terms disappear, and the solution is simply ĥ(x) = 35 4 q(x) 2/7 , Ŝ(x) = ĥ(x) .

(2.36)

For high Q, by contrast, inertia effects dominate and take the flow far from local equilibrium [START_REF] Hung | Granular flow regimes in rotating drums from depth-integrated theory[END_REF]. For intermediate Q, both inertia and rheology intervene.

The equations in that case can be solved by eliminating Ŝ between (2.33) and (2.34), yielding for ĥ the ordinary differential equation of the first order d ĥ dx = 35 9

ĥ1/2 q -4 9 ĥ4 + Q8/7 (3κ -77 24 ) ĥq dq/dx Q8/7 (2κ -77 48 )q 2 .

(2.37)

Starting from ĥ(-1/2) = 0 at the upstream boundary, this non-linear equation can be integrated numerically in the phase plane (x, ĥ) to obtain the solution curve ĥ(x) (Fig. 4(b)). Because the solution has vertical tangents at end points x = ±1/2, we integrate using discrete steps of constant distance d ξ = (dx 2 + d ĥ2 ) 1/2 rather than constant intervals dx. Next, we deduce the excess inclination profile Ŝ by substituting the solution to (2.37) back into (2.33). Finally, the free surface profile ˆ z(x) = z/(S a L/ cos α) is determined by integrating d ˆ z/dx = -Ŝ numerically subject to the integral condition

1/2 -1/2 ˆ zdx = 0.
Results are plotted in Fig. 4(c), together with the streamlines calculated from (2.16). For the experiments below, the dimensionless silo discharge Q takes relatively low values, hence the inertialess solution (2.36) provides a good approximation to the full solution of (2.37). Because the basal forcing term and discharge distribution q(x) gradually evolve as the avalanching layer descends down the silo, note that solutions must be calculated anew for each time at which a snapshot of the flow is desired.

Experiments

The experiments were performed at the granular flow laboratory of IUSTI, Aix-Marseille Univ., using the set-up shown in Fig. 5. The rectangular silo is composed of parallel glass walls separated by wooden strips, planed to uniform thickness and sandwiched between the glass plates along the perimeter. By switching between different sets of wooden strips, three gap thicknesses B 1 = 3.5 mm, B 2 = 5 mm and B 3 = 7 mm were investigated. To one side, wood pieces with mitered edges were used to control the orifice opening in the range 6 mm to 11 mm, yielding mass flowrates in the range ṁ = 3 to 12 g/s, as measured using a scale collecting the silo efflux. The lower edge of the orifice was set at a fixed height of 20 mm above the silo floor, and taken as origin (X, Z) = (0, 0) of the coordinate system. The internal silo chamber had spanwise width L = 60 mm and maximum depth H max = 600 mm.

Slightly polydisperse glass spheres were used as granular material. Two different sphere sizes were tested, having diameters D 1 = 0.5 mm and D 2 = 0.76 mm. The D 1 spheres were drawn from a batch used in previous experiments, for which the following properties were determined in earlier work [START_REF] Jop | Crucial role of sidewalls in granular surface flows: consequences for the rheology[END_REF][START_REF] Jop | Initiation of granular surface flows in a narrow channel[END_REF]: density ρ S = 2450 kg/m 3 , sphere-wall friction coefficient µ W = tan(13.1 • ), and parameters µ s = tan(20.9 • ), µ 2 = tan(32.76 • ) and I 0 = 0.279 in the non-linear µ(I) rheology

µ = τ σ = µ s + µ 2 -µ s I 0 /I + 1 , (3.1)
where I = γD/(p/ρ S ) 1/2 is the inertia number. Provided that I is not too large, this nonlinear relation between µ an I can be linearized by Taylor expansion around I = 0. This yields for the coefficients of our assumed linear rheology (2.14) the values µ 0 = µ s = 0.38 and χ = (µ 2 -µ s )/(I 0 c

1/2 0 ) = 1.2, where c 0 = 0.6 is the estimated granular volume fraction. These values of µ 0 and χ are used without adjustment for all the comparisons between theory and experiment presented below. Based on the theory described in the previous section, we can estimate as follows the range of inertia numbers encountered in the present experiments. Using the theoretical velocity profile (2.15), first, the inertia number averaged over the depth of the avalanching layer is given by

I = 35 8 uD (c 0 g ⊥ h 3 ) 1/2 . (3.2)
Using the approximation (2.36) and adopting the silo discharge Q as an upper bound on the avalanching discharge q, next, a maximum depth-averaged inertia number for each experiment can be obtained from

I max = 1 2 35 4 2/7 QD (c 0 g ⊥ h 5 a ) 1/2 . (3.3)
The different conditions tested are listed in Table 1. Measurements were acquired using a Photron high-speed camera of resolution 640 by 640 pixels operated at frame rates in the range 250 to 500 frames per second. To maximize both the resolution and coverage of the footage acquired for each experiment, the silo was mounted on a motorized track allowing it to rise at constant speed V ≈ Q/L relative to the camera. In this fashion, the camera could zoom in on an observation window of width L, within which the descending avalanching layer was nearly stationary. The lower part of the window, moreover, could scan the undisturbed deep flow, travelling past different regions of the deep flow as the silo moved relative to the camera. In one single experiment, therefore, the camera can record the deep flow over a height greater than the observation window, and capture the avalanching layer at different stages of its descent relative to the silo. Lighting was provided by a high-intensity LED, giving a point-like source of illumination that produced sharp highlights near the centers of the spheres, facilitating their capture and tracking.

Experiments were started with the silo filled up almost to maximum depth, with a level free surface. An initial transient was therefore needed to deform the free surface from its horizontal initial state to its tilted avalanching state. The time needed for a discharge q ∼ Q to resorb the corresponding triangular area A = tan αL 2 /2 is approximately T = A/Q ∼ L 2 /Q. Over this time, the surface has descended by ∆H = -ḢT ∼ QT /L = L.

To avoid this initial transient, measurements were therefore acquired only after an initial drawdown ∆H of at least 8 times the silo width L. For this reason, the silo was designed to be twice as high as the target observation domain.

[mm] D [mm] b [mm] ṁ [g/s] Q [mm 2 /s] Q [-] Imax [-] 1
Particle tracking velocimetry was performed using the methods of [START_REF] Capart | Voronoï imaging methods for the measurement of granular flows[END_REF], upgraded to take into account path regularity when matching particles from one frame to the next. For each experiment, of the order of 10 8 velocity vectors were acquired. Averaging over multiple frames was conducted by binning vectors into non-overlapping Cartesian mesh cells, using distinct frames of reference for the deep flow and for the avalanching layer. For the deep flow, stationarity in the silo frame of reference was assumed to average vectors from the lower part of the window using a mesh of coordinates (X, Z). For the avalanching layer, averaging was performed in the descending frame of reference (X, Z -Ḣt), over time intervals of shorter duration to capture different stages of the layer evolution during its descent. During each experiment, an electronic scale was also used to monitor the mass outflow from the silo. Flow fields acquired in this fashion are illustrated in Fig. 6, and compared to long exposure images acquired with the silo moving or held stationary with respect to the camera.

As a check on the PTV measurements and some of the assumptions of the theory, the mass flow rates ṁ estimated by weighting the collected efflux can be compared with mass flow rates c 0 ρ S BQ estimated from particle tracking velocimetry, where Q = -LW 0 is the volume flow rate per unit width and W 0 is the mean downward velocity of the grains averaged from 12 cross sections at different heights of the measured deep flow field. Values of both ṁ and Q are listed in Table 1. The PTV values understimate the collected efflux by up to 20 %, for the smaller grains, and by up to 10 % for the larger grains. This discrepancy indicates that the granular velocities measured by PTV along the tank front wall are somewhat slower than the granular velocities averaged over the gap thickness and associated with the bulk efflux. Since the discrepancy is not large, however, it appears reasonable to assume for simplicity that the flow is two-dimensional in character. Reflecting the crucial role of wall friction [START_REF] Jop | Crucial role of sidewalls in granular surface flows: consequences for the rheology[END_REF], the silo gap thickness was found to exert a significant influence on the observed avalanching flow patterns. For small to medium gap thickness (B = 3.5 and 5 mm), shallow avalanching layers are obtained, with lenticular shapes that span the entire silo width L (Fig. 6). The resulting flow pattern (Fig. 6(b)) is qualitatively similar to the flow field expected from the theory (see Fig. 4(c)). The largest gap thickness tested, however (B = 7 mm), produces deep avalanching layers with slow velocity dead zones upstream and downstream (Fig. 7). This pattern no longer matches that expected from our shallow avalanching theory. Judging from the comparison between outflow rates measured by the scale and by PTV, respectively, the flows in that case remain roughly two-dimensional, hence the change of character does not seem due to a breakdown of this assumption. Instead, we suspect that for these flows there is no longer a sufficient separation of scales between the deep flow and surface avalanching. For the avalanching layer, quantitative comparisons with the dynamic theory will therefore be restricted to experimental runs conducted with small to medium gap thickness. For the deep flow, by contrast, the kinematic theory will be compared with data from all runs. These comparisons are described in the next section.

Comparison and discussion

Deep flow comparison

Prior to comparing theory with experiments for the deep flow region, the kinematic constant or diffusivity K must be calibrated. This is done for each experiment using the following approach. Consider the rate of material transfer (per unit interface area) across an imaginary horizontal interface Z(X, t) = H(t) descending at speed Ḣ = -Q/L down the silo, and given by

W 0 (X, Z) = W 0 (X, Z) -Ḣ , (4.1)
where only the deep flow is considered. Assuming that the flow above this interface is steady in the descending frame of reference, this influx must be balanced by a horizontal granular current J(X, Z) (per unit gap width) given by

J(X, Z) = - X 0 W 0 (X , Z)dX . (4.2)
The control volumes associated with this balance are illustrated in Fig. 6(a). Let J max (Z) denote the maximum value of the granular current J(X, Z). Taken along horizontal interface Z = H instead of the tilted interface Z(X), it provides an approximation to the maximum avalanching discharge q max . Here we use measurements of J max for different elevations Z to calibrate the diffusivity K for each experiment. Using the deep flow theory of sub-section 2.1, the current J is given by

J(X, Z) = Ψ 0 (X, Z) + ḢX , (4.3)
or in dimensionless form

Ĵ = J Q = Ψ0 ( X, Ẑ) -X . (4.4)
where as before X = X/L and Ẑ = ZK/L 2 . Let Ĵmax = F ( Ẑ) denote the maximum dimensionless current for a given dimensionless elevation Ẑ. This relation can be inverted (numerically) to obtain for Z as a function of J max the relation

Z L = L K F -1 J max Q , (4.5)
where we have reverted to dimensional variables. This provides a simple linear relation in coefficient L/K that can be fitted to measurements of J max /Q for different relative elevations Z/L. Results obtained via this calibration procedure are presented in Fig. 8 for one of the experimental runs. In Fig. 8(a), relation (4.5) is compared with the measured data for the K value yielding the best fit. In Figure 8(b), we then compare the theoretical profiles W 0 (X, Z) obtained with this K value with the measured vertical velocity profiles. Finally, in Fig. 8(c), we compare the theoretical profiles J(X, Z) with the corresponding experimental profiles. In addition to matching well the values of maximum current J max (Z), on which the fitting was performed (Fig. 8(a)), the results show good agreement for other features of the deep flow, including the detailed velocity profiles (Fig. 8(b)) and the locus X max (Z) where the maximum current J max is attained (Fig. 8(c)).

Considering that the diffusivity K is the only parameter requiring calibration, agreement is good overall. Nevertheless, there are some significant differences between the predicted and measured velocity profiles. Most conspicuously (Fig. 8(b)), the experiments feature a zone of low velocity near the sidewall above the orifice, the thickness of which grows with height. This departure from our assumed condition of perfect slip may be partly due to the use of wood strips for the sides. The same behavior, however, was also observed by [START_REF] Maiti | Experiments on eccentric granular discharge from a quasi-two-dimensional silo[END_REF] in their experiments conducted in an eccentric silo made entirely from acrylic (polymethyl methacrylate) plates. Regardless of sidewall material, it is likely that grain mobility near the sides is reduced due to corner effects (grains in contact with both the front and sidewalls). This vertical boundary layer is not modeled by the theory, which instead assumes perfect slip along the side walls. Fortunately (Fig. 8(c)), this local discrepancy exerts a limited influence on the integrated granular current profiles J(X, Z), hence its effect on the dynamics of the avalanching layer is expected to be small.

Using the larger set of experiments, in which we varied the outflow discharge, silo gap thickness, and grain diameter (see Table 1), we can further investigate the influence of these parameters on the calibrated diffusivity K. The results, plotted in Fig. 9, are unexpectedly simple. For the range of conditions examined, K is simply proportional to the gap thickness B, satisfying to a close approximation the linear relation K = λB, where λ = 0.26. When the gap thickness is held constant, K is independent of the discharge Q, consistent with the rate-independent character of the kinematic model. Likewise, K is found to be independent of the grain diameter D.

As a consequence of this simple dependence, the Ĵmax and Xmax data for all the experimental runs should collapse together when plotted against dimensionless number Ĥ = HB/L 2 . This is checked in Fig. 10. Upon substituting Ẑ = λ Ĥ in (4.4), the data collapse and agree closely with the theoretical curve Ĵmax ( Ĥ) (Fig. 10(a)). Agreement for the locus Xmax ( Ĥ) is also good (Fig. 10(b)). The two-dimensional deep flow pattern above the orifice, therefore, is effectively the same for all cases when plotted in the dimensionless variables ( X, Ĥ). To test whether this holds more generally, one would need to vary parameters over a broader range. For instance, it is unclear if the same reduction would hold for grains and walls of materials different from glass, which would affect their frictional properties. For the present experimental conditions, however, this reduction to a single dimensionless flow pattern accounts well for the deep flow observations.

Surface flow comparison

Once the basal forcing is known from the kinematic model, the dynamic model can be used to predict the surface avalanching flow. For this purpose, we estimate the basal normal velocity w 0 or equivalently the basal streamfunction ψ 0 from the analytical solutions (2.8) and (2.9). For the diffusivity, we adopt the values K = λB estimated from the gap thickness B using the linear relation calibrated in the previous section. For the rheological coefficients µ 0 and χ, we adopt the values determined in section 3 from previous experiments by [START_REF] Jop | Crucial role of sidewalls in granular surface flows: consequences for the rheology[END_REF] using the same granular material. Detailed comparisons will be presented for two runs, corresponding respectively to gap thicknesses B = 3.5 and 5 mm. In Fig. 11, we first compare predicted and measured streamfunctions for the entire flow field, at three successive times. The streamlines shown represent equally spaced streamfunction contours, hence the distance between contours is inversely proportional to the discharge through the corresponding streamtube. For the theory, the streamfunction is obtained from

Ψ = Ψ 0 + ψ -ψ , (4.6)
where the first term represents the deep flow, the second term surface avalanching, and the third their common part that must be subtracted in order not to be counted twice. This approach is used in boundary layer problems to match inner and outer solutions [START_REF] Kevorkian | Partial Differential Equations[END_REF]. By analogy, the same approach can be used here to bridge between the avalanching and deep flow models. In our case, however, a more formal boundary layer analysis is not feasible because we do not know the more general equations that would reduce to the avalanching and deep flow models in the appropriate limits.

For the experiments, the streamfunction is obtained from the measured velocity field (U, W )(X, Z) by integrating U column-wise from the bottom and W row-wise from the far sidewall, then averaging the two results. Side by side comparisons between the theoretical and experimental streamlines obtained in this fashion are shown in Fig. 11. Going from top to bottom, the streamlines are first evenly spaced across the silo width, due the quasisteady drawdown of the free surface. They then get focused by the avalanching layer, where the streamlines become very closely spaced, to connect with the uneven deep flow streamlines at their base. Thereafter, the streamlines gradually converge toward the orifice where they again become closely spaced. Theory and experiment are seen to be in good agreement regarding this overall pattern. In Fig. 12, we compare predictions and measurements for the granular velocity fields. At five equally spaced times during the descent, we compare calculated results for the downslope velocity u(x, z, t) (Fig. 12(a,c)) with measured results for the magnitude of the velocity norm ((U -U 0 ) 2 + (W -W 0 ) 2 ) 1/2 , from which the deep flow was subtracted (Fig. 12(b,d)). On all panels, the last, lowermost snapshot shows the magnitude of the total velocity when the avalanching layer approaches the orifice level. In both theory and experiment, the avalanching layers exhibit lenticular depth profiles, with vanishing depths at both sidewalls and maximum depths h max that occur towards the orifice side of the silo. This asymmetry and the steepness of the avalanching layer become stronger as the surface descends closer to the orifice, becoming more strongly influenced by its localized pull. Surface velocities likewise accelerate as the avalanching layer drops.

Comparison of the left and right panels of Fig. 12 further shows the influence of gap thickness B on the avalanching flows. The maps of Fig. 12(a,b) correspond to gap thickness B = 3.5 mm, whereas those of Fig. 12(c,d) correspond to gap thickness B = 5 mm. For the larger gap thickness (Fig. 12(c,d)), kinematic diffusion is stronger, hence a less localized surface drawdown by the deep flow. The maximum discharge q max that the avalanching layer must transfer to even out the drawdown is therefore reduced, hence weaker avalanching flows. The larger gap thickness, moreover, reduces the relative magnitude of wall friction, leading to deeper and slower avalanching layers for the same discharge. Note that this does not imply any change in the magnitude of the wall friction coefficient, assumed constant in the theory. For increasing gap width, the magnitude of the wall friction force relative to the other terms decreases because, in the widthaveraged equation (2.13) and the ensuing depth-integrated equations, the wall friction force per unit width -2τ W /B is the only term that gets divided by the width B. For the smaller gap thickness (Fig. 12(a,b)), kinematic diffusion is weaker hence the deep flow drawdown is more concentrated above the orifice. The induced maximum avalanching discharge q max is therefore greater. Due to the smaller gap thickness, wall friction causes the avalanching layers to reduce their depth and steepen their inclination. As a result of these different effects, surface velocities are significantly increased compared to Fig. 12(c,d), even though the discharge Q through the orifice is nearly the same.

Although the free surface is nearly linear for the wider gap (Fig. 12(c,d)), for the thinner gap the surface curvature is more pronounced, transitioning from concave downward to concave upward going from upstream to downstream (Fig. 12(a,b)). As it nears the orifice wall, the avalanching layer cusps upwards as it thins to a sharp tip. Similar behavior was earlier described for symmetric silos, with avalanching layers that cusp upwards as they reach the silo axis of symmetry [START_REF] Samadani | Shocks in sand flowing in a silo[END_REF]. Although this behavior was interpreted as resulting from shock formation, our shallow flow theory captures this phenomenon without a shock. Similar behavior can be observed for flows in rotating drums, where the granular free surface goes from straight to curved as the rotation rate increases [START_REF] Rajchenbach | Flow in powders: from discrete avalanches to continuous regime[END_REF].

To make comparisons more precise, Fig. 13 plots together theoretical and experimental profiles for two depth-integrated quantities associated with the avalanching layers: the depth-integrated discharge q = udz, and the depth-integrated kinetic energy, k = 1 2 u 2 dz, both varying with height above orifice H and with the downslope coordinate x. Since the discharge profile is controlled by the deep flow, the comparisons for q represent a test of the kinematic model (Fig. 13(a,c)). For a given discharge profile, on the other hand, the kinetic energy profile is controlled by the avalanche dynamics, hence comparisons for k represent a test of the dynamic model (Fig. 13(b,d)). In both the theory and experiments, peaks in q and k occur further downstream as the granular surface drops and approaches the orifice. Also, profiles for the smaller gap thickness (Fig. 13(a,b)) are more asymmetric than profiles for the larger gap thickness (Fig. 13(c,d)). Finally, the profiles for k are more strongly peaked than the profiles for q. For all these features, fairly good agreement is registered between theory and experiment. Although the measured signal for the kinetic energy is more noisy, the level of agreement obtained for q and k indicates that both the kinematic and dynamic models perform reasonably well.

Conclusion

In the present work, theory and experiment were used to investigate granular surface avalanching in narrow silos. Provided that the deep and shallow flows exhibit a clear separation of scales, it was shown that narrow silo avalanching can be described well by combining kinematic and dynamic models, matched along an approximate basal interface. The resulting theory captures both the rate-independent character of the deep flow and the rate-dependent behavior of the avalanching layer, each in the frame of reference where they can be considered steady. For the deep flow, the measured velocities and the induced forcing at the base of the avalanching layer are described fairly well by the kinematic model, subject to a linear relation between the kinematic diffusivity and the silo gap thickness. For the surface flow, on the other hand, the shape of the avalanching layer and longitudinal profile of depth-integrated kinetic energy are well predicted by the dynamic model. Combined together, the two models also produce velocity fields and streamlines that agree fairly well with experiments over the entire domain.

More general continuum formulations have been proposed [START_REF] Staron | Continuum simulation of the discharge of the granular silo. A validation test for the µ(I) visco-plastic flow law[END_REF][START_REF] Dunatunga | Continuum modelling and simulation of granular flows through their many phases[END_REF][START_REF] Lee | A three-dimensional continuum model incorporating static and kinetic effects for granular flows with applications to collapse of a two-dimensional granular column[END_REF][START_REF] Daviet | Nonsmooth simulation of dense granular flows with pressure-dependent yield stress[END_REF], and may be able to model subsidence-driven granular avalanching under less restrictive assumptions. Nevertheless, the simple approach proposed in the present work presents various advantages. First, it applies to precisely the conditions likely to tax the capabilities of modeling approaches like discrete element simulations or Navier-Stokes solvers. The number of grains contained in the present experimental silos, for instance, would require very large computational resources to be modeled by discrete particles. Resolving the steep gradients and large convective inertia of the present avalanching layers, likewise, would require very small mesh cells and small time steps for Navier-Stokes-type computations. Another advantage of our semi-analytical approach is to clarify the physical role played by different parameters.

Nevertheless, the approach of the present paper is currently subject to various limitations. In particular, it is restricted to two-dimensional, quasi-steady flows of monodisperse grains over deep erodible beds of the same composition, assumed governed by wall friction and by the linearized µ(I) rheology. Further research will be needed to overcome these limitations and possibly address other applications like flows partly controlled by non-erodible beds [START_REF] Parez | Unsteady granular flows down an inclined plane[END_REF][START_REF] Fernández-Nieto | 2D granular flows with the µ(I) rheology and side walls friction: A well-balanced multilayer discretization[END_REF][START_REF] Sarno | Experimental investigation on the effects of the fixed boundaries in channelized dry granular flows[END_REF], and unsteady non-uniform flows like granular column collapse [START_REF] Chou | Discrete element modelling and experimental validation for the falling process of dry granular steps[END_REF][START_REF] Ionescu | Viscoplastic modelling of granular column collapse with pressure-dependent rheology[END_REF][START_REF] Lee | A three-dimensional continuum model incorporating static and kinetic effects for granular flows with applications to collapse of a two-dimensional granular column[END_REF].

A post-doctoral stay by C.-Y. Hung and a sabbatical visit by H. Capart at IUSTI, Aix-Marseille Univ., provided the opportunity for the present research. In Marseille, the hospitality of Olivier Pouliquen and Elisabeth Guazzelli is gratefully acknowledged. We also thank Yoël Forterre for providing materials and advice. Some of the theoretical ideas were earlier explored with Colin Stark. This work was undertaken under the auspices of the Labex MEC (ANR-10-LABX-0092) and of the A*MIDEX project (ANR-11-IDEX-0001-02), funded by the Investissements d'Avenir program managed by the French National Research Agency (ANR). Additional financial support was provided by the Ministry of Science and Technology, Taiwan, and the program for research excellence of National Taiwan University.

and discrete and continuous simulations of, the discharge of granular media from silos with a lateral orifice. J. Fluid Mech 829, 459-485. Zhou, Y., Ruyer, P. & Aussillous, P. 2015 Discharge flow of a bidisperse granular media from a silo: Discrete particle simulations. Phys. Rev. E 92, 062204.

Figure 1 .

 1 Figure 1. Sketch of the assumed silo configurations: (a) symmetric silo with a centered botttom orifice; (b) asymmetric silo with a side orifice.

Figure 2 .

 2 Figure 2. Solutions to the kinematic deep flow model: (a) spanwise profiles of horizontal velocity; (b) spanwise profiles of vertical velocity; (c) streamlines.

Figure 3 .

 3 Figure 3. Definition sketch of surface flow model: (a) overview; (b) local cutaway. The surface flow is depicted in a frame of reference descending with the free surface, in which the flow is approximately steady. In this frame of reference, the erosion flux across the basal interface exhibits a change of sign, from positive upstream to negative downstream, hence the apparent kink in panel (a).

Figure 4 .

 4 Figure4. Dynamic model solution for the avalanching layer: (a) imposed spanwise profile of normal to slope velocity ŵ0 = dq/dx (dashed line) and the resulting avalanching discharge q(x) (continuous line); (b) dimensionless phase diagram (x, ĥ) with the solution curve ĥ(x) in black; (c) basal interface (dashed line) and streamlines (continuous lines) in frame of reference descending with the free surface.

Figure 5 .

 5 Figure 5. Experimental setup: (a) schematic; (b) overview; (c) close-up showing the side orifice of adjustable opening.

Figure 6 .

 6 Figure 6. PTV process applied to run 3: (a) long exposure image for silo moving relative to the camera; (b) color map of the measured magnitude of the granular velocity, averaged over 250 images; (c) reconstructed color map of the deep granular velocity field, in the silo frame of reference; (d) composite long exposure image of the deep silo flow, with the silo held stationary with respect to the camera. In (a), the colored lines illustrate the control volumes used to derive (4.2). In (b), the region below the dashed line was assumed unaffected by surface avalanching, yielding the measurements used to produce the deep flow map (c).

Figure 7 .

 7 Figure 7. Qualitatively different avalanching flow pattern observed in silo of large gap thickness (B = 7 mm , run 7): (a) long exposure image; (b) color map of velocity magnitude.

Figure 8 .

 8 Figure 8. Comparison of the kinematic theory (blue curves) with experimental measurements (black symbols) for the deep flow region (run 2): (a) dependence of maximum horizontal current Jmax on elevation above orifice Z, used to calibrate the diffusivity K; (b) vertical velocity profiles W (X, Z); (c) horizontal current profiles J(X, Z), with locus of maximum values indicated by dashed line (theory) and circles (experiment). The error bars indicate ± one standard deviation, estimated from measurements of three repeat experiments for the same conditions, and multiplied for clarity by factors 20, 40 and 3 for panels (a),(b) and (c) respectively.

Figure 9 .

 9 Figure 9. Variation of the dimensionless diffusivity with the dimensionless gap thickness for all experimental runs.

Figure 10 .

 10 Figure 10. Comparison of the deep flow theory (blue line) with all experimental measurements (black dots): (a) maximum current Jmax; (b) location Xmax where the maximum current occurs.

Figure 11 .

 11 Figure 11. Comparison of predicted and measured streamlines for two different gap thicknesses, at three successive times: (a) prediction, and (b) measurements, for gap thickness B = 3.5 mm (run 3); (c) prediction, and (d) measurements for gap thickness B = 5 mm (run 5). Streamlines are plotted at equal streamfunction increments, each corresponding to 1/8 of the total outflow discharge.

Figure 12 .

 12 Figure 12. Comparison of predicted and measured granular velocity fields for two different gap thicknesses, at six successive times: (a) prediction, and (b) measurements, for gap thickness B = 3.5 mm (run 3); (c) prediction, and (d) measurements for gap thickness B = 5 mm (run 5). Colors from blue to red indicate increasing velocity magnitude.

Figure 13 .

 13 Figure 13. Comparison of theory (blue lines) and experiment (black lines) for profiles of discharge and kinetic energy, at six successive times during the descent of the avalanching layers, for two different gap thicknesses:(a) discharge profiles q(x), and (b) kinetic energy profiles k(x), for gap thickness B = 3.5 mm (run 3); (c) discharge profiles q(x), and (d) kinetic energy profiles k(x), for gap thickness B = 5 mm (run 5). Circles indicate peaks in the predicted profiles. The error bars in panels (a) and (b) indicate ± one standard deviation, estimated from three repeat experiments for the same conditions.

Table 1 .

 1 Experimental conditions for the different runs

	Run	gap	grain	orifice mass flow	silo	dimensionless	maximum
	no. thickness	size	opening	rate	discharge silo discharge inertia number
	B

† Email address for correspondence: hcapart@yahoo.com


Part II Location uncertainty (LU)

Usual random CFD

• Random parameters, boundary conditions, forcing

• Random initial conditions

• Arbitrary Gaussian forcing

• Averaging, homogenization w

Ḃ a = a(x, x) = E{ dB ( dB) T } dt Advection of tracer Θ @ t ⇥ + w ? • r⇥ + Ḃ • r⇥ = r • ✓ 1 2 ar⇥ ◆ 8
Large scales:

Small scales:

Variance tensor:

w

Ḃ a = a(x, x) = E{ dB ( dB) T } dt Advection Advection of tracer Θ @ t ⇥ + w ? • r⇥ + Ḃ • r⇥ = r • ✓ 1 2 ar⇥ ◆ 8
Large scales:

Small scales:

Variance tensor:

w

Ḃ a = a(x, x) = E{ dB ( dB) T } dt Advection Diffusion Advection of tracer Θ @ t ⇥ + w ? • r⇥ + Ḃ • r⇥ = r • ✓ 1 2 ar⇥ ◆ 8
Large scales:

Small scales:

@ t ⇥ + w ? • r⇥ + Ḃ • r⇥ = r • ✓ 1 2 ar⇥ ◆ 8
Large scales:

Small scales:

Variance tensor:

w

Ḃ a = a(x, x) = E{ dB ( dB) T } dt
Reduced models under location uncertainty • Instabilities triggered, possibly followed by extreme events

• Uncertainty quantification to address filter divergence Related works, outlooks and application

• Bifurcations (SQG) and attractor (Lorenz 63) exploration

• Comparisons with data-driven parametrisation and SALT (Stochastic Advection by Lie Transport) (Holm and coauthors)

• Parametrization and tests based on higher-order statistics (curvature, energy flux through scales, bispectrum, …)

• (Surface gravity) wave / turbulence interaction

• Data assimilation (DA) :

Filtering / smoothing

• EnKF with LU model as a R&D tool (for e.g. airplanes, drones)

• PF with reduced LU model for real-time monitoring and flow control (for e.g. polluant dispersion monitoring, drag and damage reduction in e.g. wind turbines)

  Drift correction@ t ⇥ + w ? • r⇥ + Ḃ • r⇥ = r • + w ? • r⇥ + Ḃ • r⇥ = r • + w ? • r⇥ + Ḃ • r⇥ = r • t ⇥ + w ? • r⇥ + Ḃ • r⇥ = r •

•

  Rigorous and low-cost estimators • Stabilization of the unstable modes • Maintain variability of stable modes • Uncertainty quantification A word about reduced order models 9 • Very fast simulation of very complex system (e.g. for industrial application) • Physical model (PDE) simplified using observations State of art • Possible parametrization with eddy viscosity • Impossible to parametrize with additional dissipation Need ad hoc closure like MQG (Sapsis and Majda, 2013a,b,c)

V. Resseguier -valentin.resseguier@scalian.com


Qm 1 φm 1 (xT 1 , S2), •

Maintenance optimization

In order to optimize some criterion minimize a cost: functioning, maintenance, . . . maximize a reward: availability, . . .

Our approach

propose a general model for the evolution of the equipment state based on PDMPs formalize the maintenance problem as an impulse control problem for PDMPs derive a numerical scheme to approximate the value function (with error bounds) compute the approximate optimal maintenance cost Piecewise deterministic Markov processes [Davis 93] General class of non-diffusion dynamic stochastic hybrid models: deterministic motion punctuated by random jumps.

Hybrid process X t = (m t , x t ) discrete mode m t ∈ {1, 2, . . . , p} Euclidean state variable x t ∈ R n

1 = S 1 X t = m, φ m (x , t) , P (m,x ) (S 1 > t) = e -t 0 λm φm(x ,s) ds Em x T1
Q m φ m (x , T 1 ), • E m1 Em x T1 Q m φ m (x, T 1 ), • xT 1 XIV Colloque Franco-Roumain de Mathématiques Appliquées Bordeaux August 2018 6/31
Dynamics X t follows the flow until the next jump time

T 2 = T 1 + S 2 X T 1 +t = m 1 , φ m 1 (x T 1 , t) , t < S 2 Em 1 E m x T 1 Qm φm(x, T1), • x T1

S2

XIV Colloque Franco-Roumain de Mathématiques Appliquées Bordeaux August 2018 6/31

Dynamics

Post-jump location (m 2 , x T 2 ) selected by Markov kernel

Q m 1 φ m 1 (x T 1 , S 2 ), • . . . Em 1 E m x T 1 Q m φ m (x, T 1 ), • x T1 S 2
Embedded Markov chain Mathematical definition

Strategy S = (τ n , R n ) n≥1 τ n intervention times R n new positions after intervention
Value function Value function Dynamic programming Costa, Davis, 1988 For any function g ≥ cost of the no-impulse strategy

J S (x ) = E S x ∞ 0 e -αs f (Y s )ds + ∞ i=1 e -ατ i c(Y τ i , Y τ + i ) V(x ) = inf S∈S J S (x ) f , c cost functions, α discount factor Y t controlled
J S (x ) = E S x ∞ 0 e -αs f (Y s )ds + ∞ i=1 e -ατ i c(Y τ i , Y τ + i ) V(x ) = inf S∈S J S (x )
v 0 = g v n = L(v n-1 ) v n (x ) ---→ n→∞ V(x )
XIV Colloque Franco-Roumain de Mathématiques Appliquées Bordeaux August 2018 15/31

Dynamic programming operator Example: N (0, I 2 ): Example: N (0, I 2 ):

L(w )(x ) = L(Mw , w )(x ) = inf t≤t * (x ) E x F (x , t) + e -αS 1 w (Z 1 )1 {S 1 <t∧t * (x )} +e -αt∧t * (x ) Mw (φ(x , t ∧ t * (x ))1 {S 1 ≥t∧t * (x )} ∧E x F (x , t * (x )) + e -αS 1 w (Z 1 ) with F (x , t) = t∧t * (x ) 0 e -αs-Λ(x ,s) f φ(x , s) ds Mw (x ) = inf y ∈U c(x , y ) + w (y )
-3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3 XIV Colloque Franco-
-3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3 XIV Colloque Franco-Roumain de Mathématiques Appliquées Bordeaux August 2018 19/31
Horizon and control set

finite set U of new starting points select horizon N such that v N (x ) -V(x ) small enough -→ numerical approximation of v N (x )

Main idea

Replace the dynamic programming iteration of functions by an iteration of random variables Backward dynamic programming

Change of notation

For a well chosen function g and large enough N

v N = g v n = L(v n+1 ) v 0 (x ) V(x )
XIV Colloque Franco-Roumain de Mathématiques Appliquées Bordeaux August 2018 21/31

Dynamic programming

Markov property Recursion on random variables

v n (Z n ) = L(Mv n+1 , v n+1 )(Z n ) = inf t≤t * (Zn) E F (Z n , t) + e -αS n+1 v n+1 (Z n+1 )1 {S n+1 <t∧t * (Zn)} +e -αt∧t * (Zn) Mv n+1 φ(Z n , t ∧ t * (Z n )) 1 {S n+1 ≥t∧t * (Zn)} | Z n ∧E F (Z n , t * (Z n )) + e -αS n+1 v n+1 (Z n+1 ) | Z n with F (x , t) = t∧t * (x ) 0 e -αs-Λ(x ,s) f φ(x , s) ds Mv n+1 (x ) = inf y ∈U c(x , y ) + v n+1 (y )
v n (Z n ) expression of v n+1 (Z n+1 ), Z n , S n+1 + v n+1 (y ) for all y in U Numerical scheme first compute recursively v n (y ) approximation of v n (y ) for all y in U then compute recursively v n ( Z n ) approximation of v n (Z n )

Discretization

In the expression of operator L replace inf by min over a discretized grid Step 1: Exact simulation of the PDMP Step 4: Calibrating N the number of allowed jumps + interventions

Z n , Z n+1 ,
Horizon N (number of iterations)

  Local characteristics for each mode m E m open subset of R d Flow φ m : R d × R → R d deterministic motion between jumps, one-parameter group of homeomorphisms Intensity λ m : E m → R + intensity of random jumps Markov kernel Q m on (E m , B(E m )) selects the post-jump location Dynamics X t follows the deterministic flow until the first jump time T

f

  unavailability cost proportional to time spend in failed state c fixed cost for going to the workshop + repair < change costs α = 0 (finite horizon)

  S n+1 by their quantized approximation starting fromZ 0 ∈ U n (Z n ) expression of v n+1 (Z n+1 ), Z n , S n+1 + v n+1 (y ) for all y in U Numerical scheme first compute recursively v n (y ) approximation of v n (y ) for all y in U then compute recursively v n ( Z n ) approximation of v n (Z n )

	Recursion on random variables Properties of the numerical scheme
	the quantized process has no Markov property → a different approximation of L for each time step and each starting point
	Under Lipschitz regularity assumption, convergence of the
	scheme with errors bounds depending on
	XIV Colloque Franco-Roumain de Mathématiques Appliquées Discretization In the expression of operator L replace inf by min over a discretized grid Z n , Z n+1 , S n+1 by their quantized approximation starting Bordeaux August 2018 from Z 0 = x XIV Colloque Franco-Roumain de Mathématiques Appliquées Bordeaux August 2018 the time discretization step inf → min the quantization error (Z XIV Colloque Franco-Roumain de Mathématiques Appliquées Bordeaux August 2018	23/31 23/31 24/31

v n , S n ) → ( Z n , S n )

Dynamics

Starting point

Step 2 : Discretisation of the control st U Step 5: Approximation of the value function