Didier Bazin

R Lasseux

M Lenormand

Quintard

F Golfier A'b

B Bazin

D Lasseux

R Lenormand

M Quintard

A discussion on Darcy-scale modeling of porous media dissolution in homogeneous and heterogeneous systems

A discussion on Darcy-scale modeling of porous media dissolution in homogeneous and heterogeneous systems

INTRODUCTION

Dissolution mechanism in porous media is frequently encountered in several domains: stimulation of petroleum wells by acid injection to increase the rock permeabili@, NAPL 2 transport in hydrogeology, or raisin[~ the ground level of the Dutch coast by sulfuric acid injection into subsurface limestones ~. This dissolution process is coupled with the fluid momentum equation in an unstable way: flow velocity is higher in the largest pores, which in generally produces faster dissolution processes. These processes increase locally the pore diameter and this may in turn facilitate the acid transport to these large pores. These physical mechanisms can lead to the formation of highly conductive flow channels called wormholes. As a consequence, prediction of wormhole propagation and the understanding of dissolution regimes present a major interest to optimize acidizing treatments. However, the description of dissolution macroscopic patterns is a very complex problem, determined by the flow microscopic characteristics, and with many factors in play such as the injection rate, the acid volume, and rock permeability. At the pore-scale, the dissolution mechanism involves acid transport by diffusion and advection to the solid surface, chemical reaction at the surface and product transport away from the surface. If the chemical reaction characteristic time is very short compared to the mass-transfer kinetics, the reaction is called mass-transfer-limited, and this is generally the case for limestone dissolution with HC1. On the other side, if mass-transfer kinetics is slow, then the reaction is reaction-rate limited, and this is the case for dolomite dissolution at the room temperature. Several experiments of dissolution have been performed in a variety of fluid-mineral systems, for mass transfer limited as well as reaction-rate limited processes. Hoefner et al. 4 and Fredd et al. 5 used limestone core samples and HC1 injection. The wormhole structure is visualized by injecting a low melting point alloy and then dissolving the porous medium or by neutron radiography. In addition to HCl-limestone systems, Several researchers 6'7 investigated the effect of temperature, acid concentration, rock mineralogy and injection rate for some Indiana or Glenn Rose limestone and dolomite cores. All the experiments led to the observation of similar dissolution regimes: face dissolution, conical wormhole, dominant wormhole, ramified wormholes and uniform dissolution. Another result is the existence of an optimum injection rate 4'6. It corresponds to the maximum penetration of the wormhole for a given volume of acid injected. It is reported that the optimum conditions are related to the formation of a dominant wormhole with little branching through the core. The optimum flow rate depends on several parameters, including the rock mineralogy (calcite or dolomite for carbonate formations), the temperature and the acid concentration. In order to predict this optimum injection rate and the wormhole development at the Darcyscale, several numerical models have been developed. Among the different approaches found in the literature, we can distinguish 9 the network models 4'5 which consider the dissolution of the solid grains in a two or three dimensional structure at the pore-scale, or the capillary tube model 6, based on a pre-existing shape of the wormhole. We can also cite the dimensionless apl~roach 5'8 based on correlations at the scale of the core sample itself, or the model of Liu et al." which was a first attempt to develop a macroscopic dissolution model. The reader can refer to the paper by Fredd et al. ~~ for a detailed presentation of these different models. Nevertheless, since none of these models could describe all the features of the dissolution physics, in particular they failed to describe the coupled nature of flow and reaction without assuming a wormhole geometry, a macroscopic acid-transport model has been implemented based on a macro-scale description involving a non-equilibrium mass balance model ~'~2. Concerning the momentum equation, assuming that the velocity of the pore-scale interfaces is small enough, we have adopted a Darcy-Brinkman formulation which allows to simulate correctly the flow in the fluid or porous zones in a continuous manner 13. The effect of dissolution history on the effective properties, such as the permeability, is approximated by direct relationships between the macro-scale parameters. For instance, a direct relationship is adopted between the permeability and the porosity, which is a classical assumption made in geochemistry. We will present in this paper original results obtained by this model for both homogeneous and heterogeneous systems and for 2D and 3D configurations. Finally, we will use the 2D results to discuss the possible features of a core-scale dissolution model.

PRESENTATION OF THE DARCY-SCALE MODEL

The dissolution macroscopic equations can be obtained from the coupled equations of flow and species transport at the pore-scale, and this lead to a local non-equilibrium dissolution problem. Most of the theoretical questions associated with this development have been discussed in Quintard et al. We just recall below the Darcy-scale model we have used and some of the various results that have been obtained. The solid phase and the fluid phase are identified as the o-phase and the I3-phase respectively. In order to be clear about the equations we have used, we will follow the notations associated with the volume averaging theory presented in Quintard et al. 14. For instance the superficial velocity is given by

1 ~vadV v, V p-vhase (I)
where Vfl-plu~se represents the volume of the 13-phase contained within the averaging volume, V, represented Figure 1. The intrinsic mass concentration is defined by

C A fl --( C A fl ) fl -" ~l I C A fl d V (2)
Vfl-phase Vi3_;~s e and we have a similar definition for the intrinsic pressure. We consider that the reaction is mass-transfer limited, and the acid is immediately consumed as soon as it reacts with the solid surface. Additional convective terms which appear during the upscaling are neglected. By this way, we obtain a complete system of equations which combines a Darcy-Brinkman model for the flow description in both zones (fluid and porous region) with a local non-equilibrium dissolution model. It can be written under dimensionless form as follows Flow equations:

/z AV# -VP, -/zK~.V# = 0 (3) c#
where Ix is the viscosity, K and @ are the permeability tensor and the porosity respectively. This equation is coupled with the overall mass balance equation for the fluid phase which may be approximated by 6). This will eventually lead to a new porosity field which will be used in the next iteration. In practice, the permeability field is directly obtained from the porosity field, by the use of a classical Kozeny-Carman relation K -el3. This particular choice is not a constraint of the model, and other correlations can be used. Various calculations have been performed on two-dimensional and three-dimensional domains, and the different results are discussed in the next section.

RESULTS AND DISCUSSION

Homogeneous system

In a recent paper ~ we have shown that the simulations allow to capture the different dissolution regimes for "Darcy-scale homogeneous systems". We can see in Figure 2 the different 2D dissolution figures obtained numerically for different acid injection rates. All these figures represent the porosity field. The obtained dissolution patterns are remarkably equivalent to the experimental dissolution patterns obtained for calcite dissolution 5'~~ or for water/salt system l~. These results allowed us to study in some details the transitions between the different regimes in terms of behavior diagrams. In addition, the performed simulations

Heterogeneous systems

Such a model can also be applied easily to heterogeneous systems. We consider a stratified system made of two different porous media of permeability K~ = 10 ll m 2 and K2 = 10 s m 2 respectively. We will limit here our study to two simple limit cases" a domain in which the acid solution is injected perpendicularly to the strata and a domain in which the acid solution is injected parallel to the strata. The numerical data for both simulations are given in Table 1. The results in Figure 5 shows how complex is the evolution of the unstable dissolution front due to the heterogeneity effects: early stage in one stratum similar to the homogeneous case, then interaction of the dominant wormhole with the second medium, and subsequent complicated evolution, different from the homogeneous case. On the contrary, for the flow parallel to the strata, Figure 6, dissolution occurs essentially in the more permeable region like in the homogeneous case. This illustrates the interest of the proposed model. However, numerical limitations exists, and domain of very large extent are beyond our possibilities. A larger-scale averaged model, or core-scale model (Figure 1), if available, would open new possibilities, and this is the subject of the next section.

CORE-SCALE MODELING

In this section, we discuss the introduction of the macroscopic equations which may be used to describe the problem at the core-scalenE. For the flow description at this scale, a classical Darcy's law is obtained by averaging the Darcy-Brinkman formulation. With regard to the transport and dissolution part, it is not obvious to apply an upscaling method which leads to some core-scale equations valid in a general way. In fact, the value of the mass transfer coefficient can strongly modify the acid transport behavior and the corresponding Darcy-scale dissolution pattern.

Figure 1 .

 1 Figure 1. Different scales of the problem.

V.Vt~ =0 ( 4 )

 4 Species transport equations :OC Afl ep /)t + Vp "VCap = V.(D* .VCap)-ct Cap the stoechiometric coefficient of the chemical reaction. The mass transfer coefficient ct and the dispersion tensor D* are calculated from closure problems 14. A numerical model has been developed to solve these equations I~. Given the different mechanisms involved, the numerical model uses a multiple-step approach. First, the stationary Darcy-Brinkman problem is solved for a given porosity field by a predictor-corrector method coupled with an Uzawa algorithm. The effective coefficients in the transport equations are obtained previous to the simulation as a function of the porosity and velocity field by numerically solving on a representative periodic unit cell the pore-scale closure problems presented in Golfier et al. 1~. The acid concentration profile resulting from the transport equation can be used to solve for the dissolution equation, Eq. (

Figure 2 .

 2 Figure 2. Porosity fields representative of the dissolution patterns (From Golfier et al.~l).It must be emphasized here the importance of the choice of the Darcy-Brinkman formulation for the flow model. In fact, we could have used a classical Darcy / Darcy model, which is simpler compared to the Darcy-Brinkman formulation. However, it requires the introduction of an appropriate value for the permeability in the wormhole, i.e., in the fluid zone (noted Knuia and called fluid permeability). This concept of fluid permeability is artificial and it must be chosen by a posteriori arguments to recover the correct physics. For the Darcy-Brinkman model, the fluid permeability needs only to be large enough so that the model leads asymptotically towards the Stokes equation, i.e., the Darcy term vanishes as the porosity is close to one and the fluid permeability does not play any role. On the contrary, with the Darcy-Darcy formulation, the obtained velocity varies with the fluid permeability value. Therefore, an appropriate value of the fluid permeability must be chosen to model correctly the flow into the wormhole. If this value may be estimated in some cases, flow between two parallel plates, for example, this estimation is practically impossible in a 3D flow without assuming a tube geometry for the wormhole.

Figure 3 .

 3 Figure 3.3D porosity field in conical regime.

Figure 4 .

 4 Figure 4. 3D porosity field in wormholing regime.

Figure 5 .

 5 Figure 5. Wormhole length as a function of the pore volume injected and associated dissolution pattern for strata in series.

Figure 6 .

 6 Figure 6. Wormhole length as a function of the pore volume injected and associated dissolution pattern for strata in parallel.

Figure 7 .

 7 Figure 7. K*-~wcorrelation in the wormholing regime.

Figure 8 .

 8 Figure 8. Wormhole length as a function of time for different models. 5 CONCLUSION A three-dimensional Darcy-scale dissolution model has been presented in this paper. Its originality lies in the fact that it is based on a combination of a Darcy-Brinkman model coupled with a local non-equilibrium dissolution model. Simulations have been performed for 2D and 3D configurations for both homogeneous and heterogeneous systems and present a good agreement with the literature. Based on the 2D results, different approaches have been explored to model dissolution phenomena at the core-scale. A double-porosity model has shown a relatively good agreement with the direct simulations, at least under local equilibrium conditions, compared to a simple one-equation model. This indicates that complicated dynamic effects are present which may require the use of a more complex oneequation model. Developments are made in this direction.

One-equation model. First, a simple one-equation model can be proposed similar to the Darcy-scale model presented in the previous sections.

. +c:, . . (D" " ) " "

~'fl ()t "]" Vfl "VCAfl -" V.

.VCAfl --a CAr l [START_REF] Bazin | SPE 30085[END_REF] 04 / dci The dissolution patterns obtained numerically at the Darcy-scale are used to determine the correlations for the core-scale effective coefficients (macroscopic permeability K*, core-scale mass transfer coefficient a*, fluid fraction #w and macroscopic porosity e'p) by spatial integration. The obtained results show an acceptable independence of these coefficients with respect to the dissolution history at least for certain conditions (wormholing regime in particular, see an example of correlation in Figure 7).

From the system of equations obtained at the core-scale and the mass transfer and permeability correlations provided by the numerical simulations at the Darcy-scale, a 1D core-scale model can be developed for both approaches to verify the validity of such a representation. To check the validity, we compared predictions provided by these core-scale models with direct Darcy-scale simulations. Figure 8 presents the comparison of these simulations for Q = 50 cm3.h 1 and Cap = 545 kg.m 3. If the results obtained by the oneequation model are not satisfactory, the two-equation model correctly predicts the wormhole propagation and presents a relatively good agreement with the Darcy-scale model. This suggests that a simple, fickian, dispersive model with chemical reaction is not well suited to model the dissolution of a porous medium at the core-scale in the presence of wormholes. If one is looking at a one-equation model, it must be more complicated. Work is currently undertaken in this direction.


Grégoire Allaire
email: gregoire.allaire@polytechnique.fr

Samuel Kokh
email: samuel.kokh@cea.fr

Test-case number 19: Shock-Bubble Interaction (PN)

Introduction

We describe here a test proposed by [START_REF] Quirk | On the dynamics of a shock-bubble interaction[END_REF] based on the experiments of [START_REF] Haas | Interaction of weak shock waves with cylindrical and spherical gas inhomogeneities[END_REF]. Its goal is to simulate the propagation of a shock through a helium bubble in air. Haas and Sturtevant initial purpose was to get a better understanding of the Richtmeyer-Meshkov instabilities. More generally they wanted to understand how pressure waves in heterogeneous media can generate turbulent phenomena which tend to mix the fluids. From a numerical point of view the goal of this test is to validate compressible multifluid flows models as well as numerical methods used for solving these models. This test has been performed at least in the following studies by [START_REF] Abgrall | How to prevent pressure oscillations in multicomponent flow calculations: a quasi conservative approach[END_REF], [START_REF] Fedkiw | A non-oscillatory eulerian approach to interfaces in multimaterial flows (the ghost fluid method)[END_REF], [START_REF] Karni | Hybrid multifluid algorithms[END_REF], [START_REF] Kokh | Numerical simulation of 2-D two-phase flows with interface[END_REF] and [START_REF] Saurel | A simple method for compressible multifluid flows[END_REF].

Description

Geometry. The test is two-dimensional. The computational domain is a rectangular box which is 890 mm long (horizontal axis) and 89 mm high (vertical axis). At time t = 0, the bubble has a 50 mm diameter and its center is located at (x c , y c ), x c = 420 mm, y c = 44.5 mm (the origin being the low left corner of the domain). The initial location of the shock is a vertical line which is 222.5 mm away from the right side of the domain. Physical Model. behavior of each fluid is governed by the gas dynamics compressible Euler equations (without any diffusion term, neither surface tension, nor gravity). Each fluid is assumed to obey the perfect gas equation of state. Thus the fluid i = 1, 2 is modeled by the following equations,

Gas

γ R (kJ.kg -1 .K -1 ) C v (kJ.kg -1 .K -1
           ∂ t ρ i + div (ρ i u i ) = 0 ∂ t (ρ i u i ) + div (ρ i u i ⊗ u i + p i I) = 0 ∂ t (ρ i e i ) + div [(ρ i e i + p i )u i ] = 0 (1)
where ρ i is the density, u i the velocity (it is a two-component vector), e i the specific total energy such that e i = ε i + |u i | 2 /2 with ε i the specific internal energy, the pressure p i being provided by,

p i = (γ i -1)ρ i ε i , (2)
where γ i is the ratio of the heat capacities of the i th gas. The interface modeling and its numerical treatment are free choices (most of the simulation references use an isothermalisobaric mixture law in order to thicken the interface which is then captured on an Eulerian mesh).

Numerical Data. The equation of state parameters for air and the helium bubble are provided in table 1 (data for R and C v are, a priori, not required). Let us note that the parameters for the helium bubble describe indeed a 28% mass mixture between helium and air. The shock travels from the right side of the domain to the left side with a 1.22 Mach velocity. This means its velocity is 1.22 times higher than the sound velocity in the pre-shock air at rest (atmospheric pressure and density equal to 1 kg.m -3 ). Let us recall that for a perfect gas, the sound velocity is given by c = γp/ρ and that the horizontal shock velocity is -1.443523 × 10 3 m.s -1 . The helium bubble is supposed to be initially at mechanical equilibrium with the surrounding air. Thanks to the Rankine-Hugoniot conditions it is possible to find the initial values for air after the shock. Using compatible units with those used in the equation of state, the initial state is defined in table 2 (let us recall that: 1 bar = 10 5 Pa, and Pa = 1 J m -3 ).

Boundary conditions. The horizontal boundaries of the domain are solid walls where "mirror" boundary conditions are to be applied (i.e. non-penetration conditions). The right vertical boundary is set to be an "inflow" boundary condition equal to the initial data for air to the right of the shock. The left vertical boundary is treated as a free "outflow", which means a zero order extrapolation of the variables has to be performed out of the computational domain.

Measures and Comparisons. It is required to compare the shape of the bubble with Haas and Sturtevant experimental results at the following time steps: 32, 52, 62, 72, 82, 102, 245, 427, and 674 µs. We shall also plot the pressure evolution in time downstream from the position of the bubble.

Remark. Data provided in [START_REF] Fedkiw | A non-oscillatory eulerian approach to interfaces in multimaterial flows (the ghost fluid method)[END_REF] are different from those in [START_REF] Quirk | On the dynamics of a shock-bubble interaction[END_REF] which are being used here. There exists a similar test where the helium is replaced with a gas heavier than air (refrigerant fluid R22, see [START_REF] Quirk | On the dynamics of a shock-bubble interaction[END_REF])

Remark. Let us note that similar tests are available (see e.g. [START_REF] Allaire | A five-equation model for the simulation of interfaces between compressible fluids[END_REF][START_REF] Shyue | A fluid-mixture type algorithm for compressible multicomponent flow with van der Waals equation of state[END_REF], however no direct comparison is possible as these tests deal with different equations of state than those which are proposed here.

Figure 1 :

 1 Figure 1: Computational domain.

Table 1 :

 1 Equation of state parameters.

	)

Table 2 :

 2 Initial state.


Jean-Baptiste Camps
email: jean-baptiste.camps@chartes.psl.eu

Julien Randon-Furling
email: julien.randon-furling@cantab.net

Lost Manuscripts and Extinct Texts: A Dynamic Model of Cultural Transmission

Keywords: Agent-based models, Stochastic models, Loss of cultural artefacts, Text transmission, Stemmatology

How did written works evolve, disappear or survive down through the ages? In this paper, we propose a unified, formal framework for two fundamental questions in the study of the transmission of texts: how much was lost or preserved from all works of the past, and why do their genealogies (their "phylogenetic trees") present the very peculiar shapes that we observe or, more precisely, reconstruct? We argue here that these questions share similarities to those encountered in evolutionary biology, and can be described in terms of "genetic" drift and "natural" selection. Through agent-based models, we show that such properties as have been observed by philologists since the 1800s can be simulated, and confronted to data gathered for ancient and medieval texts across Europe, in order to obtain plausible estimations of the number of works and manuscripts that existed and were lost.

Introduction

How much do we preserve of the written knowledge, science [START_REF] Cisne | How science survived: Medieval manuscripts' "demography" and classic texts' extinction[END_REF] and culture of the past? And how representative is what we know compared to what existed? Such fundamental questions depend on the process through which texts were distributed materially.

Before the advent of the printing press, written texts were circulated in manuscript form. In order to make the text available, the author would dictate it to a secretary, or write a draft on wax tablets, papyrus, parchment or, eventually, paper, and this original, authorial, manuscript would then have to be copied manually by a scribe in the form of a new manuscript, and then circulated.

Copies could then be used to create more manuscripts, again by manual copying, perhaps by other scribes in other regions at a later date. During this process, successive modifications were introduced in the text, either by error or intentionally, to make the text more suited to its intended audience. These alterations in the written sequence that forms the text could then be transmitted to its "descendants" by a given manuscript. Wear and tear, accidents, fashions caused the destruction of some manuscripts, while others enjoyed the long life of library preservation. In the end, knowledge was lost and some texts went extinct, while other texts gained traction or were eventually preserved for future generations.

In textual studies and philology, since the development of the "common errors" methods from the 19th century onward [START_REF] Timpanaro | La genesi del metodo del Lachmann[END_REF][START_REF] Haugen | 2 The genealogical method[END_REF], the analysis of the alterations to the text allows the philologists to reconstruct the relations between the surviving copies (witnesses) of a given work, and to represent them as a tree-like graph, called a stemma codicum (fig. 1). This reconstructive process is not entirely different from the methods used by biologists to reconstruct the links between existing or extinct species, based on their shared characteristics that are supposed to derive from a common ancestry, and to represent it as a phylogenetic tree. Methods from this biological subfield, known as cladistics, have even sometimes been directly applied to texts, with controversial results [START_REF] Barbrook | The phylogeny of the Canterbury Tales[END_REF][START_REF] Howe | Manuscript evolution[END_REF][START_REF] Spencer | Phylogenetics of artificial manuscripts[END_REF][START_REF] Macé | The Evolution of Texts: confronting stemmatological and genetical methods[END_REF][START_REF] Hoenen | History of computer-assisted stemmatology[END_REF]. These trees represent the result of an enquiry into the relationships between the surviving witnesses, together with the hypothetical lost nodes that can be deduced from them. To arrive at it, researchers have to study the variants -more precisely, common innovations or errors (i.e., mutations) -observed in the text of the surviving witnesses (fig. 1,A). The reconstructed tree will show only what can be deduced from surviving witnesses: the witnesses themselves, and as many hypothetical nodes as are needed to explain their relationships (fig. 1,B). While it may be the case with more recent works that all nodes are known and the graph represents the full transmission of the text -for instance with the genealogy of printed editions (fig. 1,C) -, most of the time the tree represents only a (potentially very small) subset of what existed (fig. 1

, D).

A long standing observation, not yet fully understood, about the structure of the reconstructed trees was made by the French philologist Joseph Bédier, almost a century ago, in 1928 [START_REF] Bédier | La tradition manuscrite du lai de l'ombre. réflexions sur l'art d'éditer les anciens textes (premier article)[END_REF]. He observed that most trees reconstructed by philologists show a root bifurcation (a root with outdegree 2): in most reconstructions, the original (or archetype) had two, and only two, direct descendants, preventing an accurate reconstruction of the original text by majority principle (e.g., two witnesses vs one). Instead of searching an explanation for such data in the dynamics of text transmission, he interpreted this "forest of bifid trees" as the result of a methodological flaw or an unconscious bias, spurring a century long debate and causing a long lasting methodological schism inside textual scholarship that remains to be fully resolved [START_REF]L'ombre de Joseph Bédier: théorie et pratiques éditoriales au XXe siècle[END_REF][START_REF] Duval | La 'Tradition manuscrite du Lai de l'Ombre' de Joseph Bédier ou la critique textuelle en question[END_REF].

Bedier's initial observation has been replicated, with estimates of the proportion of root bifurcation varying, from Bédier's 95.5%, to somewhat lower estimates ranging from 70% to 83% [START_REF] Shepard | Recent theories of textual criticism[END_REF][START_REF] Castellani | Bédier avait-il raison?: La méthode de Lachmann dans les éditions de textes du Moyen Age[END_REF][START_REF] Haugen | The silva portentosa of stemmatology: Bifurcation in the recension of Old Norse manuscripts[END_REF]. Yet, some have argued that the prevalence of root bifurcation could be an explainable feature of manuscript transmission of texts as it reached us. Tentative explanations include combinatorial estimations of the proportion of root bifurcation for a given number of witnesses, under the assumption that all configurations are equally likely [START_REF] Maas | Leitfehler und stemmatische Typen[END_REF][START_REF] Castellani | Bédier avait-il raison?: La méthode de Lachmann dans les éditions de textes du Moyen Age[END_REF][START_REF] Hoenen | How many stemmata with root degree k?[END_REF], or consider the effects of decimation (i.e., manuscript loss) [START_REF] Greg | Recent theories of textual criticism[END_REF], for instance by applying a uniform loss probability to static preexisting trees [START_REF] Guidi | Sugli stemmi bipartiti. decimazione, asimmetria e calcolo delle probabilità[END_REF] or by calculating a node specific loss probability to simulated trees [START_REF] Hoenen | Silva Portentosissima -Computer-Assisted Reflections on Bifurcativity in Stemmas[END_REF]. Even if it generated little follow-ups, there have also been rare attempts of using birth and death process for exploring the dynamics of manuscript transmission [START_REF] Weitzman | Computer simulation of the development of manuscript traditions[END_REF][START_REF] Weitzman | The Evolution of Manuscript Traditions[END_REF].

The abundance of roots -and more generally nodes [START_REF] Haugen | The silva portentosa of stemmatology: Bifurcation in the recension of Old Norse manuscripts[END_REF] -with out-degree 2, is not the only property that can be observed in many stemmas. The asymmetry between branches (cf. fig. 1, C andD), the presence or not of lateral transmission (generally called "contamination"; fig. 1,E) are other properties worthy of investigations, as well as those that can indicate that Figure 1: From the real tree to the reconstruction. A an artificial example of the transmission of a very short text and its plausible reconstruction, where circles depict extant manuscripts, rectangles lost ones), and gray, nodes that are recoverable when reconstructing the tree; B another artificial example, where the distribution of lost manuscripts causes the loss of a whole branch, and where, consequently, the last common ancestor is not the original (root), but a lower witness (archetype); C the observed phylogeny of the successive editions of Fortunio, Regole grammaticali della volgar lingua (1516) [START_REF] Richardson | Giovan Francesco Fortunio: Regole grammaticali della volgar lingua[END_REF][START_REF] Guidi | Sugli stemmi bipartiti. decimazione, asimmetria e calcolo delle probabilità[END_REF]; D the reconstructed phylogeny (stemma) of the Old French Song of Roland [START_REF] Segre | La chanson de Roland, number 16 in Documenti di filologia[END_REF]; E the reconstructed phylogeny of the Anglo-French Guy de Warewic (Ewert, 1932), showing many cases of lateral transmission ("contamination") in dashed lines, and even one instance of lateral transmission from outside the tree; F the same tree, but including (box shape) the outline of a lost branch (cf. [START_REF] Hoenen | The stemma as a computational model[END_REF]).

the tree made from extant witnesses represents only a small portion of the original tradition (i.e., lateral transmission from outside the tree; root identifiable not with the original but with a later manuscript; fig. 1,F). It is reasonable to assume that some of these properties reflect the dynamics of manuscript transmission, while others keep trace of important destruction, decimating manuscripts and removing even full branches. For the texts, this can be seen as an evolutionary process, where two antagonist tendencies are at work: the apparition of textual variants in individuals, causing the increase of diversity in the tradition, and the extinction of full branches, causing some variants to prevail upon others and so reducing diversity.

Here again, these observations can be put in perspective with problems occurring in evolutionary biology, where, too, two antagonists tendencies are at work, mutation and fixation, either by drift or natural selection, in a context where processes of speciation and extinction are strongly linked and where extant species represent only a small subset of the species that have existed [START_REF] Yessoufou | Reconsidering the Loss of Evolutionary History: How Does Non-random Extinction Prune the Tree-of-Life?[END_REF]. In both cases, survival might be the exception and extinction the rule, be it by "bad genes or bad luck", a process that can be seen in terms of gambler's paradox [START_REF] Raup | Extinction: bad genes or bad luck?[END_REF][START_REF] Canettieri | Philology and Information Theory[END_REF].

Basic processes of reproduction and destruction create complex shapes in the trees, from which one might want to deduce whether they can be fully explained by random process akin to genetic drift, or if differences in selective values are to be suspected. In other terms, going back to textual traditions, if cultural context, through literary taste, canon or fashion for instance, creates a form of evolutionary pressure on textual traditions. Any insights gained on this question would have an applicability beyond the question of the transmission of antique and medieval texts, because it seems that similar dynamics are at work in the diffusion of content in other medium, including print [START_REF] Guidi | Sugli stemmi bipartiti. decimazione, asimmetria e calcolo delle probabilità[END_REF] or even the web [START_REF] Barabási | Emergence of Scaling in Random Networks[END_REF], and have been observed in areas such as the cognitive evolution of scientific fields and the dynamics of scientific memes [START_REF] Bentley | Random Drift versus Selection in Academic Vocabulary: An Evolutionary Analysis of Published Keywords[END_REF][START_REF] Chavalarias | Phylomemetic Patterns in Science Evolution-The Rise and Fall of Scientific Fields[END_REF].

Data on cases as different as the songs of the Medieval Occitan troubadours from southern France or the incunabula editions printed in Renaissance Italy outline the same Pareto-like world, where a large number of texts are kept only in a single or handful of documents, while a limited number of "successful" texts are kept in a large number of copies (fig. 2, A andB), where most authors are known only for one or two texts, while a very limited number of writers can have dozens of texts preserved (fig. 2, C andD). . . Such a process is also apparent in the constitution of a literary canon of a limited number of authors and texts. This 'canonialisation' can be seen has a progressive loss of diversity, where an ever shrinking number of authors and texts take on an ever growing share of the circulated documents (fig. 2, E andF). But is this due to chance or to a selective process ?

In fact, some properties as were just mentioned for textual traditions have some pendants in evolutionary models, concerning for instance the very unequal distribution of descendants [START_REF] Cosette | Single cell dynamics causes Pareto-like effect in stimulated T cell populations[END_REF], the varying patterns of biodiversity varying in time and space, studied in macroecology and biogeography [START_REF] Cabral | Mechanistic simulation models in macroecology and biogeography: state-of-art and prospects[END_REF][START_REF] Rangel | Modeling the ecology and evolution of biodiversity: Biogeographical cradles, museums, and graves[END_REF] and the dynamics of speciation and extinction that manifest themselves in the shape of phylogenies and the loss of branches from the tree of life [START_REF] Mace | Preserving the Tree of Life[END_REF][START_REF] Yessoufou | Reconsidering the Loss of Evolutionary History: How Does Non-random Extinction Prune the Tree-of-Life?[END_REF]. For this reason, there are inspiration and resources to be found in the study of mathematical properties of evolutionary trees, regarding the establishment of a null-model [START_REF] Bienvenu | The split-and-drift random graph, a null model for speciation[END_REF].

Can we gain some insight on what existed, what was lost, and the driving forces between extinction or survival of texts -by drift or selection? In order to do so, we need a better understanding of the dynamical processes of manuscript transmission.

An agent-based model of a stochastic process

In this paper we present a selection of the first results obtained with a stochastic model for the transmission of manuscripts in the Middle Ages. Following Weitzman [START_REF] Weitzman | Computer simulation of the development of manuscript traditions[END_REF][START_REF] Weitzman | The Evolution of Manuscript Traditions[END_REF], we use so-called birth-and-death processes. These are random processes introduced in probability theory to describe, among other things, simple population dynamics and genealogies. For the simplest versions of theses processes, it is possible to derive analytically (i.e. with mathematical formulae) certain quantities of interest: the expected number of individuals (here for us, manuscripts) still present at a time 𝑡, the extinction probability, the survival probability,... But also, quantities of particular interest in the context of manuscript genealogies: the probability that the latest common ancestor (lca) be an archetype rather than the original, or the probability of root bifidity for the reconstructed stemma. Here we favour a numerical approach, that allows us to explore more complicated variants of a birth-and-death process through an agent-based computer simulation. The agents correspond to manuscripts, and during each time step of the simulation each agent has a probability 𝜆 of being copied (once) and a probability 𝜇 of disappearing. Following Cisne, we limit the increase rate of the population by letting the birth rate 𝜆 depend on the size of the population at the previous step 𝑘 𝑡-1 :

𝜆 𝑡 = 𝜆 𝐾 𝑘 𝑡-1 , (1)
where 𝐾 is a theoretical upper bound on the maximum possible size of the manuscript population extant at any given time during the period under consideration.

Phase diagrams obtained through computer simulations

We ran agent-based simulations of a Cisne-type tradition with a total time frame of 500 pseudoyears -that is, a time-step in the simulation corresponds roughly to 3 months. This was derived using an estimate for the time taken to produce an average 200 page manuscript (though the speed of scribes is known to vary a lot, from 1 to 10 leaves a day) [START_REF] Overgaauw | Fast or slow, professional or monastic. The writing speed of some latemedieval scribes[END_REF]. We chose 𝐾 = 100 000 = 10 5 , as an order of magnitude (rather than 10 4 or 10 6 ) for the total number of manuscripts in a given tradition that could be extant simultaneously at any one time 1 .

As for the remaining free parameters, namely the base "birth" or copy rate, 𝜆, and the "death" rate 𝜇, we explored all possible pairs of values of these parameters within their range (from 0 to 1 in theory, but reduced here to 10 -4 to 10 -3 by field expertise, i.e. rough estimates from philological knowledge) 2 . The space of all possible values is usually called the phase space in physics and other mathematical sciences, and a representation of the value taken by a given 1 This order of magnitude is a very rough estimate based on human population and our own assumptions: the medieval population of countries such as France or Italy was in the 10 7 range; we estimate that the maximum saturation of this market for a given book would be reached if around 1% of the population were to own a copy. 2 Historical and philological knowledge of loss rates is very scarce and elusive, but can still be approached from various angles, such as the collection of data from ancient library catalogues, inventories, wills, as well as allusions and intertextuality [START_REF] Wilson | The Lost Literature Of Medieval England[END_REF][START_REF] Bardon | La littérature latine inconnue[END_REF][START_REF] Buringh | Medieval Manuscript Production in the Latin West[END_REF]. Buringh [START_REF] Buringh | Medieval Manuscript Production in the Latin West[END_REF] provides estimates for the Latin West, with a geometric mean of loss around -25% per century, with variations from -11% in the 9th to -32% in the 14th and 15th centuries (with local variations between medieval institutions from -3% to -71% per century). The global loss rate for non-illustrated manuscripts of several well known collections have been estimated around 93-97% [START_REF] Kestemont | Estimating the Loss of Medieval Literature with an Unseen Species Model from Ecodiversity[END_REF][START_REF] Neddermeyer | Von der Handschrift zum gedruckten Buch: Schriftlichkeit und Leseinteresse im Mittelalter und in der frühen Neuzeit quantitative und qualitative Aspekte[END_REF][START_REF] Wijsman | Luxury Bound: Illustrated Manuscript Production and Noble and Princely Book Ownership in the Burgundian Netherlands (1400-1550[END_REF][START_REF] Oostrom | Stemmen op schrift: geschiedenis van de Nederlandse literatuur vanaf het begin tot 1300[END_REF]. But there is a potential bias in accounting only for well known institutional collections, from which some manuscripts are known to have survived: trying to account for fully lost libraries, Buringh [START_REF] Buringh | Medieval Manuscript Production in the Latin West[END_REF] is compelled to revise his estimates higher, to -25% by century until the 12th, up to -43% in the 15th. For incunabula, using editions whose original number of print made is known, it is possible to gather loss estimates by counting known surviving exemplars in public or private collections: doing so for Venetian incunabula, Trovato [START_REF] Trovato | Everything you always wanted to know about Lachmann's method: a nonstandard handbook of genealogical textual criticism in the age of post-structuralism, cladistics, and copy-text[END_REF] finds very variable loss rates according to textual and material typology, from 73% for the Decretales printed on parchment to 99.3% for more popular chivalrous literature (Orlando furioso for instance). This shows the importance both of variation in time and space, and of textual contents and material typology. In some extreme cases, loss can be very close to 100%, for reasons that may combine the fragility of the document form, lack of consideration for the documents or large scale historical events such as political instability, invasions or major cultural changes; examples are provided by cases as different as the Merowingian royal diplomas on papyrus or the Lombard royal charters [START_REF] Ganz | Charters Earlier than 800 from French Collections[END_REF], the Mayan (pre-colombian) manuscripts or medieval notarial acts [START_REF] Holtz | Überlieferungs-und Verlustquoten spätmittelalterlicher Herrscherurkunden[END_REF]. Production estimates have also been attempted on the basis of the observable quantity (eg the extinction probability) when parameters are varied across the phase space is called the phase diagram of this quantity. We thus produced phase diagrams for a number of relevant observables (Fig. 3). Note that, since there is a stochastic component in the model, each phase diagram is computed by averaging over the results of a relatively large number of simulations, here 100 -this means that we produced 100 artificial manuscript traditions for each pair (𝜆, 𝜇) ∈ [10 -4 , 10 -3 ], varying values by increments of 10 -4 ; hence each phase diagram required 10 000 simulations.

The approach then consists in identifying, within the phase diagrams (also called heat maps), regions in which the values for the observables are consistent either with measured quantities (as is done in the natural sciences) or with estimates for these quantities coming from other, independent and altogether different, models. We have circled in red such regions on the phase diagrams in Fig. 3.

These heat maps show that the results obtained through these simulations are internally consistent in terms of not only population size and survival rates, but also in terms of structural properties of the resulting trees. In particular, the results obtained for a ratio 𝜇 𝜆 between 5 8 and 6 7 are surprisingly consistent with the observed properties of some medieval traditions, in particular those from chivalric narratives in Old French. In particular, values of 0.55 for the survival of works and 0.05 for the survival of manuscripts (fig. 3), A and B, red squared area, bottom-right tile) are identical to those provided by Kestemont et al. for Old French chivalric romances, using unrelated methods from ecodiversity [START_REF] Kestemont | Forgotten books: The application of unseen species models to the survival of culture[END_REF]. Yet, for what regards specifically Old French epics, known as chanson de geste -a genre predating the later form of the roman, and whose circulation and reception considerably differs for a long time -, the median final population of 2 and the third quartile of LCA outdegree of 2 (though, median LCA Shannon index is 0.69) 3 . This would lead us to revise Kestemont et al. estimates to 0.22 (instead of 0.55) for the specific survival of Old French epics (chansons de geste) and 0.01 (instead of 0.5) for the survival of epic manuscripts (a figure closer to that observed by Trovato [START_REF] Trovato | Everything you always wanted to know about Lachmann's method: a nonstandard handbook of genealogical textual criticism in the age of post-structuralism, cladistics, and copy-text[END_REF] for their later Italian successors). On the other hand, specific values, this time, for (Arthurian and Antique matters) romans yield a third quartile of LCA outdegree of 3 more coherent with Kestemont quantity of sealing wax acquired by a given producer (a chancellery for instance [START_REF] Bautier | Introduction[END_REF]). More founded loss estimates have also been gathered by counting how many of the acts mentioned in imperial or royal registers are kept in original or consigned in the archives of the recipients: this gives a loss rate of originals varying from 80% (acts from the emperor Charles IV in 1360-1361) to 90% for the acts from Louis X of France, increasing to 99% for the judgements rendered by his Parliament, suggesting here as well a massive effect of typological variation [START_REF] Holtz | Überlieferungs-und Verlustquoten spätmittelalterlicher Herrscherurkunden[END_REF][START_REF] Canteaut | Quantifier l'activité des chancelleries à l'aune de la tradition des actes : l'exemple de la chancellerie des derniers Capétiens (1314-1328)[END_REF], resulting in very strong biases in the body of documents available to us. For our simulation needs, if we start from Trovato's estimates (potentially more reliable, because based on editions whose original number of copies is known), we get a survival rate whose order of magnitude is between 0.1 and 0.01 (between 10 -2 and 10 -3 ) in 500 years (2000 steps in our model), in similar ranges as Buringh's 0.75 in 100 years and those reported by Holtz and Canteaut; from this we can deduce a step loss rate for a given total survival rate. For instance, for 1% survival rate, (1 -𝜇) 2000 = 0.01, which simplifies to 𝜇 = 0.002. So we retain values of 𝜇 between 10 -4 to 10 -3 . Given that books could not have been produced order of magnitudes higher or slower than they were destroyed (or we would be either drown in medieval manuscripts or keep none), we explore the same range for 𝜆. Of course, fixed rates are a limitation, and do not yet account for substantial variations in time (such as massive extinction events, like, e.g., the fall of the Roman empire, the fire in Alexandria library, the shift from volumen to codex or from caroline to gothic script, etc.). et al. general estimate (median tradition size, according to Martina [START_REF] Martina | La produzione manoscritta del romanzo francese in versi : modelli materiali e modelli di cultura, These de doctorat[END_REF], is 2 -and mean 4.8for the sole romans en vers, similarly to chansons de geste, but is expected to be higher for later romans en prose). ; B survival rate of manuscripts (i.e., nodes); C the median final population of surviving witnesses for traditions with at least one; second row contains data about the age of manuscripts in the simulation, for D all manuscripts ever created; E extant ones; as well as F the date of the oldest extant manuscript; the third row concerns structural properties of the trees themselves, F the median distance between the lowest common ancestor of the surviving manuscripts and the actual root of the original tree; G the third quartile of the LCA (archetype) out-degree and H and the median Shannon index for the families (the main branches stemming from the LCA). For each pair of parameter values between 0.0001 and 0.001, 100 simulations were run for 1000 active and 1000 inactive steps. Grey areas correspond to irrelevant values of the parameters, or unstable regions. The red square shows parameter regions where the observables computed on the simulated manuscript populations are consistent with observations made for Medieval French epics or with plausible estimates made using different methods.

The situation is, for the moment (and until further data is acquired) consistent yet deserving of further inquiry concerning the distribution of age of surviving manuscripts: in the simulation's red-squared area, the median date of surviving manuscripts would be in the 800's step (around 200 years after the original) and the median date of the oldest for each tradition in the 150-200 years range. For chansons de geste, the median date of surviving manuscripts would be between 1250 and 1300 [START_REF] Vitale-Brovarone | La diffusion manuscrite des chansons de geste: une vue d'ensemble[END_REF][START_REF] Camps | La 'Chanson d'Otinel': édition complète du corpus manuscrit et prolégomènes à l'édition critique[END_REF] -150 to 200 years after the first documents of the genre itself (the end of the 11th century for the composition of the oldest version of the Roland). For romans en vers, it is, like the genre itself, slightly offset in time, with a peak between 1275 and 1325 [START_REF] Martina | La produzione manoscritta del romanzo francese in versi : modelli materiali e modelli di cultura, These de doctorat[END_REF].

Discussion

Combining previous inquiries by Weitzman and Cisne with the power of computer simulations and the methodology of statistical physics, we are able to reproduce the evolutionary process that underlies the observable data for, at least, some textual traditions such as those from medieval French epics and romances. The results obtained can even corroborate or refine results obtained by unrelated methodologies, such as those recently published by Kestemont et al. [START_REF] Kestemont | Forgotten books: The application of unseen species models to the survival of culture[END_REF], indicating that these relatively simple birth-and-death process have relevancy in philology as well as they have in Evolutionary Biology for instance. This method then provides us a way to account both for population dynamics in time, loss or production estimates, as well as the shape of the stemmata (the phylogenies) of manuscripts, answering the century long Bédier observation [START_REF] Bédier | La tradition manuscrite du lai de l'ombre. réflexions sur l'art d'éditer les anciens textes (premier article)[END_REF], whose lack of solution until now has been at the core of a lasting schism in philological studies.

The range of further investigations opened by this research is considerably large. Models using individual variable rates of 𝜆 and 𝜇 could be used to account for phenomenon such as efforts of preservation of old and venerable artefacts, or higher selective values of some copies, or accelerated destruction due to small scale (e.g., burnt libraries), larger scale (e.g., the Dissolution of English monasteries, French Wars of Religion,. . . ) or global events (e.g., shift in book types such as from volumen to codex or caroline to gothic scripts, major cultural changes like the Renaissance, . . . ). Modelling should also include the actual variation of the texts, to reflect the introduction of variants (mutations) in some families, and processes of transmission of inherited variants, as well as lateral transmission.

More generally, once having established this 'null model', deviations due to different factors should be explored, such as higher selective values for some mutations (variants) or individuals, fluctuations in time and space and the existence of different 'ecological niches' (e.g., the Anglo-Norman public versus the readers of Franco-Italian epics), typological variation in books or texts, chocks and bottlenecks, etc. The question of the age of surviving manuscripts should also be explored and accounted for, especially in the light of potential variations of 𝜆 and 𝜇 in time. For instance, the demand and rate of copy for a given text could be expected to be highest shortly after its initial release, when it is most fitted to the taste and fashion of the time, perhaps reinforce itself if the text gets a quick breakthrough, and then decrease with the passing of years. Similarly, the rate of destruction could vary at a global or local level, as some shocks lead to peaks of destruction or canonicalisation and conservation efforts lead to lower rates. Last but not least, the model should account for non standard transmission, in particular lateral transmission (contamination), a process not uncommon in textual transmission but that is also encountered in the natural world (e.g., lateral gene transfer).

The generality of the models considered here makes them applicable not only to medieval texts, but to any type of cultural transmission, at least in written form, from manuscript circulation to the elaboration of a canon of works. Further investigations should try to verify it on the broadest possible range of cases, starting with Western Medieval and Antique texts, but preferably encompassing cultural productions from very different time periods and continents.

Figure 2 :

 2 Figure 2: A Pareto-like world with diminishing diversity. A distribution of witnesses per troubadour text; B distribution of the extant copies per work for incunabula printed in Italy; C distribution of texts per author for troubadours and D incunabula; E Shannon diversity, Generations as sites, texts as species, witnesses as individuals; F number of texts and witnesses per author.

Figure 3 :

 3 Figure3: Heat maps (phase diagrams) for the simulation of Cisne-like models; first row contains population size and survival properties, namely A Survival rate of traditions (i.e., trees); B survival rate of manuscripts (i.e., nodes); C the median final population of surviving witnesses for traditions with at least one; second row contains data about the age of manuscripts in the simulation, for D all manuscripts ever created; E extant ones; as well as F the date of the oldest extant manuscript; the third row concerns structural properties of the trees themselves, F the median distance between the lowest common ancestor of the surviving manuscripts and the actual root of the original tree; G the third quartile of the LCA (archetype) out-degree and H and the median Shannon index for the families (the main branches stemming from the LCA). For each pair of parameter values between 0.0001 and 0.001, 100 simulations were run for 1000 active and 1000 inactive steps. Grey areas correspond to irrelevant values of the parameters, or unstable regions. The red square shows parameter regions where the observables computed on the simulated manuscript populations are consistent with observations made for Medieval French epics or with plausible estimates made using different methods.

Data about the traditions of the chansons de geste follow Vitale-Brovarone's[START_REF] Vitale-Brovarone | La diffusion manuscrite des chansons de geste: une vue d'ensemble[END_REF] and Camps'[START_REF] Camps | La 'Chanson d'Otinel': édition complète du corpus manuscrit et prolégomènes à l'édition critique[END_REF]. Information on the shape of stemmata have been computed based on a restriction to chansons de geste and deduplicating of the collection provided by OpenStemmata[START_REF] Camps | Open Stemmata: A Digital Collection of Textual Genealogies[END_REF][START_REF] Camps | Open Stemmata: Database[END_REF].