T Heitz

C Giry

B Richard

F Ragueneau
email: frederic.ragueneau@ens-paris-saclay.fr

Identification of an equivalent viscous damping function depending on engineering demand parameters

Keywords: Dissipations, Reinforced concrete, Viscous damping, Experiments, Model identification

When performing a nonlinear time-history analysis of a reinforced concrete structure, it is necessary for the used structural model to dissipate the correct amount of energy. For the sake of computational efficiency, viscous damping models are still commonly used to account, partially or not, for non-viscous dissipations (e.g. friction between the crack surfaces, bond slip at the steel-concrete interface). In order to improve the physical relevance of such a substitution, an evolving equivalent viscous damping ratio estimated for a simply supported reinforced concrete beam is proposed in this paper. This work takes place in the scope of a moderate seismicity context for which steel yielding is not expected. The results are not directly identified from experimental results but rather from numerical simulations carried out thanks to an equivalent single-degree-of-freedom model, itself calibrated by means of quasi-static experiments. To begin with, the experimental setup used to calibrate the single-degree-of-freedom model and the equivalent viscous damping ratio assessment method are presented. Then, the single-degree-of-freedom model and the identification procedure are exposed. The resulting outputs are presented and commented. Finally, numerical experiments are performed in order to obtain equivalent viscous damping ratio values corresponding to a given maximum time-history curvature and a curvature demand.

Introduction

Despite the increasing accuracy of the models dedicated to the nonlinear behavior of reinforced concrete (RC) structures, their combination with complex finite element (FE) meshes still leads to high computational cost. In practice, an additional viscous damping is often used to account for dissipations not taken into account by the structural model (e.g. [START_REF] Ragueneau | Damage model for concrete including residual hysteretic loops: application to seismic and dynamic loading[END_REF][START_REF] Crambuer | Experimental characterization and modeling of energy dissipation in reinforced concrete beams subjected to cyclic loading[END_REF][START_REF] Charney | Unintended Consequences of Modeling Damping in Structures[END_REF][START_REF] Correia | Seismic energy dissipation in inelastic frames: Understanding state-of-the-practice damping models[END_REF]). However, it is important to keep in mind that civil engineering structures often do not show purely viscous damping. Some authors [START_REF] Frizzarin | Damage detection based on damping analysis of ambient vibration data[END_REF] state that cracks within damaged elements dissipate a significant amount of energy and this dissipation is better described by a Coulomb friction model rather than a linear viscous one, due to the bond slipping between steel and concrete. This paper focuses on the case of earthquakes that do not induce steel yielding in structures (generally, weak to moderate intensity earthquakes). This choice is based on the French regulatory seismic hazard map [START_REF]Plan Séisme[END_REF] in figure 1 where the vast majority of the country presents at most a moderate seismic hazard (scale 3 over 5). In particular, it is considered that the main dissipations occur in the concrete (e.g. crack initiation, propagation and friction). This choice is also supported by the French regulation regarding the design of nuclear power plants where steel yielding is generally considered as non-acceptable. Combining hysteretic and viscous dissipations can compromise the validity of a study [START_REF] Charney | Unintended Consequences of Modeling Damping in Structures[END_REF][START_REF] Hall | Problems encountered from the use (or misuse) of Rayleigh damping[END_REF] and could require a reduction of the viscous damping in the nonlinear range [START_REF] Correia | Seismic energy dissipation in inelastic frames: Understanding state-of-the-practice damping models[END_REF]. In order to describe this dependency, several evolving Rayleigh-type viscous damping models have been proposed [START_REF] Jehel | Initial versus tangent stiffness-based Rayleigh damping in inelastic time history seismic analyses[END_REF], but experimental evidences on slender buildings tend to show that an adequate modal viscous damping could be an intermediary between constant viscous damping for all modes and stiffness-proportional damping [START_REF] Cruz | Evaluation of Damping Ratios for the Seismic Analysis of Tall Buildings[END_REF][START_REF] Satake | Damping evaluation using full-scale data of buildings in Japan[END_REF]. Before performing such a fine analysis, simplified methodologies should be used in order to obtain a first design of a structure. Alternatively, displacement-based methods [START_REF] Priestley | Displacementbased seismic design of structures[END_REF] have became popular for the past decades, but they generally require the knowledge of a capacity curve and the value of an equivalent viscous damping ratio (EVDR). It is worth mentioning alternative approaches to model nonlinearities, such as the ones relying on the concept of ductility [START_REF] Fajfar | A Nonlinear Analysis Method for Performance-Based Seismic Design[END_REF][START_REF] Lagaros | Evaluation of ASCE-41, ATC-40 and N2 static pushover methods based on optimally designed buildings[END_REF]. In this analysis, the damping ratio has a key role when assessing maximum structural responses (i.e. reduction of the demand spectra to represent the nonlinearities of the structure). However, some studies show it is the second source of uncertainties after the ground motions [START_REF] Celik | Seismic fragilities for non-ductile reinforced concrete frames Role of aleatoric and epistemic uncertainties[END_REF][START_REF] Lee | Evaluation of equivalent damping ratio of a structure with added dampers[END_REF], leading so to uncertainty and wrong quantification of the structural response. Furthermore, this assumed equivalence with the energy dissipated hysteretically may be subjected to discussions. Indeed, several studies have emphasized the fact that the viscous damping depends on the response amplitude for RC buildings [START_REF] Li | Field measurements of amplitude-dependent damping in a 79-storey tall building and its effects on the structural dynamic responses[END_REF][START_REF] Liu | Observed natural frequencies, damping ratios, and mode shapes of vibration of a 30-story building excited by a major earthquake and Typhoon[END_REF][START_REF] Satake | Damping evaluation using full-scale data of buildings in Japan[END_REF] and on their degradation state often described by a displacement ductility level [START_REF] Rosenblueth | On a kind of hysteretic damping[END_REF][START_REF] Gulkan | Inelastic responses of reinforced concrete structure to earthquake motions[END_REF][START_REF] Iwan | The effective period and damping of a class of hysteretic structures[END_REF][START_REF] Kowalsky | Displacement based design: a methodology for seismic design applied to RC bridge columns[END_REF][START_REF] Priestley | Myths and fallacies in earthquake engineering, revisited[END_REF][START_REF] Elmenshawi | Hysteretic energy and damping capacity of flexural elements constructed with different concrete strengths[END_REF][START_REF] Rodrigues | A comparative analysis of energy dissipation and equivalent viscous dampung of RC columns subjected to uniaxial and biaxial loading[END_REF]. In fact, there is actually no reason to keep the EVDR constant throughout the nonlinear time-history analysis (NLTHA). For instance, [START_REF] Charney | Unintended Consequences of Modeling Damping in Structures[END_REF][START_REF] Correia | Seismic energy dissipation in inelastic frames: Understanding state-of-the-practice damping models[END_REF] have shown that EVDR is not constant in that case due to the continuous change in the modal properties of the structure. To sum-up, the aforementioned papers suggest that the EVDR depends on the degradation state of the structure and on the displacement demand. A correct quantification of energy dissipation according to the degradation state is of main importance for design purposes as well as for the estimation of margin. With a method like the linear equivalent approach, a small change of the EVDR leads to non-negligible differences in the pseudo-acceleration and in the displacement observed by the structure. For instance, according to the formula provided in [START_REF]Eurocode 8 Design of structures for earthquake resistance Part 1: General rules, seismic actions and rules for buildings[END_REF] (1), changing the EVDR from 1% to 2% (ξ=ξ 0 + ξ eq = 5 + ξ eq ) while considering the nonlinearities leads to a decrease of the design load comprised between 5% and 9%.

η = 10 5 + ξ (1)
Furthermore, it has been observed and commented that even in the case of a more complex nonlinear computation (see for instance [START_REF] Richard | SMART 2008: Overview, synthesis and lessons learned from the International Benchmark[END_REF][START_REF] Richard | SMART 2013: Lessons learned from the international benchmark about the seismic margin assessment of nuclear RC buildings[END_REF]), a small change in the value of the initial viscous damping chosen for the analysis leads to large differences in the response (displacement, floor response spectra,).

The RC structural elements show dissipation capacity in their pre-yielding regime. Indeed, nonlinear phenomena also occur in concrete and at its interfaces with steel reinforcements. The quantification of these dissipations is of main concern in the energy production domain. Indeed, the capacity of the structural components to dissipate energy influences directly the loading seen by the sensitive equipments. In the context of accidental loading such as earthquake, one should ensure the functioning of these equipments to maintain the global safety.

The influence of the degradation state of the structure and the displacement demand is investigated in this paper. In the first section, the experimental setup used to calibrate the beam equivalent single-degree-of-freedom (SDOF) model and the EVDR assessment method are presented. Then, the SDOF model formulation and its identification procedure are detailed. Finally, numerical experiments are performed on this model in order to obtain EVDR values corresponding to a maximum relative history displacement and a relative displacement demand.

Experimental campaign

Context and motivations

An experimental campaign has been set up on RC beams by means of the AZAL ÉE shaking table, as part of the TAMARIS experimental facility operated by the French Alternative Energies and Atomic Energy Commission (CEA). The main objective is to provide the scientific and engineering communities reference data through both dynamic and quasistatic tests in order to evaluate the dissipations depending on structural, material and signal characteristics. The measure of damping forces is challenging because their levels are unknown and generally more difficult to investigate than acceleration or restoring forces [START_REF] Jehel | A critical look into Rayleigh damping forces for seismic performance assessment of inelastic structures[END_REF]. Furthermore, they originate from several sources and evolve along the time-history analysis. In order to quantify and to analyze the dissipation and its evolution during quasi-static and dynamic loadings, a dense network of various sensor types has been mounted on the experimental setup. These indirect measurements associated to the hypothesis that all dissipations are related to internal dissipations (i.e. concrete nonlinearities and rebars/concrete sliding) allow for the assessment of the hysteretically dissipated energy. The latter hypothesis is fulfilled by an important design effort on the technological choices for the samples and the boundary conditions as described comprehensively in the Ph.D. thesis of the first author [START_REF] Heitz | Nonlinear local behaviors and numerical modeling of damping in civil engineering structures in dynamics[END_REF] and more concisely in a dedicated paper [START_REF] Heitz | Dissipations in reinforced concrete components: static and dynamic experimental identification strategy[END_REF]. A focus is made in this paper on the quasi-static full cyclic loading of one RC beam. This section describes briefly this part of the experimental campaign setup which is necessary to provide the data used in this work.

RC Beams

For the first beam (referenced as HA12-C1B-1), the reinforcement is composed of ten 12mm-diameter reinforcing steel bars in the section of the beam. The standard concrete C25/30 (in the sense of Eurocode 2 [START_REF]Eurocode 2: Design of concrete structures -Part 1-1: General rules and rules for buildings[END_REF]) has a 29.7 MPa strength and 28.1 GPa Young's modulus. The reinforcement drawing is given in figures 2a and 2b. For the second beam studied in section 5 (referenced as HA20-C1A-1), the reinforcement is composed of four 20mm-diameter reinforcing steel bars in the section of the beam. The concrete is a standard C25/30 with a 35.0 MPa strength and 26.2 GPa Youngs modulus. The steel reinforcement drawing is given in figures 2c and 2d. More details regarding the geometrical and material characteristics of the beams can be found in [START_REF] Heitz | Nonlinear local behaviors and numerical modeling of damping in civil engineering structures in dynamics[END_REF][START_REF] Heitz | Dissipations in reinforced concrete components: static and dynamic experimental identification strategy[END_REF]. In this study, the global methodology and the identification of the EVDR is developed for the beam HA12-C1B-1. The second beam (HA20-C1A-1) is used in the last part of the paper to validate the methodology and the influence of the degradation state of the structure and the displacement demand on the EVDR.

Setup

To better understand and visualize the items described in the following sections, a general view of the final experimental setup is presented in figure 3a and schematized in figure 4. Only the quasi-static part of the experimental campaign is considered in this paper. Elastic hinges based upon high performance steel blades are used at both ends to link the beam with the shaking table or the strong floor. The blades have been designed to behave elastically all along the tests. This technical solution avoids by this way the dissipations due to friction or assembly clearance. The weight of the beam is supported by air cushions that considerably reduce the friction forces between the beam and the ground during bending. Indeed, supporting the vertical forces through the beam-end supports would have caused bending cracks due to the 5.9 meters span.

Equivalent viscous damping identification procedure

A viscous damping model is a practical -but not entirely satisfactory -way to describe the dissipations occurring in a RC structure to which a dynamic loading is applied. Jacobsen [START_REF] Jacobsen | Steady forced vibration as influenced by damping[END_REF] proposed a method to assess an optimal equivalent viscous damping ratio from an energetic point of view. This approach is briefly described below.

3.1. Assessment of the viscous damping ratio for a linear massspring-damper system Let us consider the classical linear mass-spring-damper problem pictured in figure 5, where the different notations used are presented. For an harmonic ground acceleration üg (t) = -U g • ω 2 • cos(ω • t), the steady-state displacement and velocity responses are given by equations ( 2) and (3), where U is the response displacement amplitude, ω is the excitation angular frequency and Φ is the phase angle:

u(t) = U • cos(ω • t -Φ) (2)
u(t) = -U • ω • sin(ω • t -Φ) (3)
The energy E d dissipated by the viscous damper during one displacement half-cycle Γ u starting at t = 0 and ending at t = π ω is expressed in equation ( 5) with ω 0 the natural angular frequency of the system:

E d = Γ u F D • du = π ω 0 c • u(t) 2 dt = π ω 0 c • U 2 • ω 2 • sin(ω • t -Φ) 2 • dt (4)
Hence:

E d = π 2 • c • ω • U 2 (5)
For a purely linear spring, the maximum stored elastic energy is reached at displacement U:

E s = U 0 k • u • du = 1 2 • k • U 2 (6)
From equations ( 5) and ( 6), and

provided that c = 2•m•ξ •ω 0 and k = m • ω 2
0 , the viscous damping ratio writes:

ξ = 1 2 • π • ω 0 ω • E d E s (7)
At this point, it is clear that the viscous damping ratio of a linear oscillator not only depends on the ratio of the dissipated energy during one half-cycle over the maximum stored energy during this same half-cycle, but also on the ratio of the natural angular frequency of the oscillator over the excitation angular frequency.

Jacobsen's areas method

The expression in equation ( 7) is exact for the aforementioned ideal case. Jacobsen [START_REF] Jacobsen | Steady forced vibration as influenced by damping[END_REF][START_REF] Jacobsen | Damping in composite structures[END_REF] proposed to extend its use to non-viscous damping. The method was originally used to evaluate the equivalent viscous damping ratio (EVDR) for nonlinear frictional system. Since the publication in 1960, many authors have used the here so-called Jacobsen's areas method, with recurrent simplification hypotheses:

• For a linear viscous damper, the dissipated energy over a cycle is linearly proportional to the excitation frequency (see equation ( 5)). This compensates the inverse function of ω in equation ( 7) and gives a constant viscous damping ratio ξ. For non-viscous damping, the dissipated energy does not depend on the excitation frequency. Hence, the EVDR actually depends on the excitation angular frequency ω. Except if the excitation frequency is close to the eigenfrequency of the oscillator (i.e. ω 0 /ω = 1), this dependence should be taken into account, otherwise, a significant error on the EVDR is made as illustrated in figure 6. For the sake of simplicity, the maximum stored energy E s is often evaluated as the area of the right triangle linking the origin and the maximum displacement point. However, in its original publication, Jacobsen [START_REF] Jacobsen | Damping in composite structures[END_REF] rather suggests to use the area under the actual restoring force curve which he calls "skeleton", as depicted in figure 7b to consider nonlinear elasticity.

• When the loops are not symmetric, there is no actual reason to pick up the maximum relative displacement rather than the minimum one when assessing the stored elastic energy. For this reason, Kumar et al [START_REF] Kumar | Cyclic Response of Sand Using Stress Controlled Cyclic Triaxial Tests[END_REF] proposed an approach adapted to asymmetric hysteretic behaviors in soils as depicted in figure 7c. A method inspired from this work and more suitable for nonlinear restoring forces is proposed in figure 7d. The skeleton is computed as the mean force for a given displacement (reached during both loading and unloading). This prevents problems related to softening behaviors.

These observations make the full-cycle skeleton-based method depicted in figure 7d a reasonable choice to assess the EVDR from the full cyclic tests results further exposed. This last method is used after for the computation of the EVDR. Furthermore, the EVDR is identified for the eigenfrequency of the oscillator.

(a) Linearly stored energy evaluated on a half-cycle (b) Nonlinearly stored energy evaluated on a half-cycle (c) Linearly stored energy evaluated on a full-cycle proposed by [START_REF] Kumar | Cyclic Response of Sand Using Stress Controlled Cyclic Triaxial Tests[END_REF] (d) Nonlinearly stored energy evaluated on a full-cycle inspired by [START_REF] Kumar | Cyclic Response of Sand Using Stress Controlled Cyclic Triaxial Tests[END_REF] Figure 7: Different ways to apply Jacbosen's areas method [START_REF] Jacobsen | Damping in composite structures[END_REF][START_REF] Kumar | Cyclic Response of Sand Using Stress Controlled Cyclic Triaxial Tests[END_REF] 4. Identification of a hysteretic model on experimental data 4.1. Quasi-static full cyclic test A quasi-static full cyclic loading, labeled QSC1, is carried out thanks to two actuators mounted at positions L 4 and 3•L 4 with an increasing cycle amplitude. Each cycle of this displacementdriven loading is repeated three times to stabilize the damage state of the beam. The corresponding time-displacement evolution is given in figure 8. To interpret the whole displacement response of the beam as that of a simple degree of freedom oscillator (SDOF), a full-field measurement based upon a commercial digital image correlation method has been used (Videometric [START_REF] Videometric | Videometric[END_REF]). This method provides displacement data throughout the beam, allowing for the projection over the eigenmode shapes [START_REF] Heitz | Dissipations in reinforced concrete components: static and dynamic experimental identification strategy[END_REF][START_REF] Heitz | Damping identification and quantification: experimental evidences and first numerical results[END_REF]. The validation of the assumption of similar displacement fields between 4-point bending test and dynamic loading is exposed in [START_REF] Heitz | Dissipations in reinforced concrete components: static and dynamic experimental identification strategy[END_REF]. In the present case, the presented displacement data result from the projection of the displacement field of the beam on the theoretical first mode shape computed thanks to a finite element model of the experimental setup. Thus, they can be considered as an estimation of the first mode displacement: v(x, t) = q 1 (t)φ 1 (x), with q 1 (t) the first mode displacement and φ 1 (x) the first mode shape. The corresponding forcedisplacement graph is represented in figure 9. In order to describe the behavior of the beam, a nonlinear oscillator is associated to the first eigenmode. The nonlinear behavior of the first mode has been observed experimentally with a frequency drop from 6.79 Hz to 2.83 Hz (beam HA12-C1B-1) and from 7.16 Hz to 2.97 Hz (beam HA20-C1A-1). This observation is explained by the nonlinear phenomena described and modeled in section 4.2.2. This model is not intended to be general but rather adapted to the QSC1 test (i.e. bending of the beam following the first mode shape). It provides a useful basis for numerical experiments that will be carried out on virtual quasi-static tests in section 5. Moreover, in paper [START_REF] Heitz | Dissipations in reinforced concrete components: static and dynamic experimental identification strategy[END_REF], a spacetime projection error map is presented, showing no significant error growth during the degradation process. These numerical experiments investigate the uncoupled influence of both the loading amplitude (e.g. displacement, curvature, force, bending moment) and degradation state.

Modeled phenomena

The model presented in this section strongly relies on the work of Moutoussamy [START_REF] Moutoussamy | Essais hybrides en temps réel sur structures de Génie Civil[END_REF]. For small displacements applied at the initial state, the beam in bending has a linear elastic behavior, i.e. has a constant stiffness K which can be identified. Beyond a limit displacement δ d (i.e. displacement at which first crack occurred in concrete) , the beam starts to exhibit nonlinearities. A stiffness decrease is observed during the unloading phase, this indicates the creation and propagation of cracking and motivates the use of the damage mechanics framework [START_REF] Lemaître | A course on damage mechanics, 2nd Edition[END_REF].

Steel yielding is not observed during the tests and is consequently not taken into account. Nevertheless, hysteresis loops are observed during loading-unloading cycles. It is explained by the existence of friction within the reinforced concrete, e.g. between cracks surfaces or between the steel reinforcements and the surrounding cracked concrete. Hence, a sliding displacement variable u π is defined [START_REF] Ragueneau | Fonctionnement dynamique des structures en béton -Influence des comportements hystérétiques locaux[END_REF][START_REF] Richard | Continuum damage mechanics based model for quasi brittle materials subjected to cyclic loadings: Formulation, numerical implementation and applications[END_REF].

The unilateral effect is taken into account by splitting each internal variable related to damage in two parts. These parts correspond to two independent families of cracks on upper or lower half cross-sections of the beam from either side of the neutral axis depending on the direction of deflection. Supposing a linear behavior of concrete in compression, the beam behavior is mainly driven by the tension area (and the associated damage state). The index i will be used to stand for both direction indexes "+" and "-". A kinematic hardening variable α π is associated to the friction displacement u π .

Finally, the last observed phenomenon is the so-called "pinching" effect which consists in a stiffness reduction in the neighborhood of the zero-displacement point that explains the reversed S-shape of the load displacement curve. The origin of this effect is not clearly understood but different hypotheses exist. Several researchers state that this effect indicates a failure driven by shear [START_REF] Brown | Reinforced Concrete Beams Under Load Reversals[END_REF][START_REF] Penzien | Experimental Investigation into the Seismic Behavior of Critical Regions of Reinforced Concrete Components as Influenced by Moment and Shear[END_REF][START_REF] Bertero | Seismic Behavior of Ductile Moment-Resisting Reinforced Concrete Frames[END_REF][START_REF] Otani | Nonlinear dynamic analysis of reinforced concrete building structures[END_REF][START_REF] D'ambrisi | Modeling of Cyclic Shear Behavior in RC Members[END_REF][START_REF] Vintzileou | Experimental validation of seismic code provisions for RC columns[END_REF][START_REF] Archundia-Aranda | Behavior of reinforced concrete haunched beams subjected to cyclic shear loading[END_REF]. However, a shear failure is unlikely to occur with such a flexible beam in four-point bending test. Hence, two explanations are here proposed:

Mechanism

Observable • in the case of a crossing crack, the crack surfaces get into contact when the displacement is high enough (as illustrated in figure 10) thus enlarging the virtual cross section of the beam at the crack location;

• the steel-concrete bond failure induces slipping until the adherence is found again at the interface (inter-locking).

State potential

The thermodynamical framework [START_REF] Lemaître | A course on damage mechanics, 2nd Edition[END_REF] of the model is defined in this part. Helmoltz's free energy is chosen as state potential, which is a function of the internal and observable variables summarized in table 1:

Ψ(u, d i , u π , α π , η) = Ψ e (u, d i , η) + Ψ π (u π , α π , η, d i ) (8)
For a linear elastic behavior:

Ψ e = 1 2 • K • u 2 (9)
To take into account the stiffness degradation, the damage variable d i is introduced:

Ψ e = 1 2 • K • (1 -d i ) • u 2 (10)
The pinching effect only exists along with damage. This is taken into account thanks to a crack closure variable η that evolves from 0 (cracks completely closed) to 1 (cracks opened):

η = d i • exp - u U c (11)
U c corresponds to the crack closure displacement. Hence, the altered stiffness K p is defined as follows:

K p (η) = K • (1 -η) (12)
which means that the pinching appears progressively with the damage. Examples of "pinching curves" for a fully damaged beam (i.e. d i = 1) are given in figure 11. Then:

Ψ e = 1 2 • K p (η) • (1 -d i ) • u 2 (13)
If the friction state potential is now considered:

Ψ π = 1 2 • K p (η) • d i • (u -u π ) 2 + 1 2 • b π • (α ) 2 (14)
To summarize, the overall state potential writes:

Ψ = 1 2 •K p (η)•(1-d i )•u 2 + 1 2 •K p (η)•d i •(u-u π ) 2 + 1 2 •b π •(α π ) 2 (15)

State laws

The next step is to formulate the state laws by differentiating the state potential in equation ( 15), namely: 16)

F e = ∂Ψ ∂u = K p (η) • (1 -d i ) • u + K p (η) • d i • (u -u π ) (
F π = - ∂Ψ ∂u π = K p (η) • d i • (u -u π ) (17)
X π = - ∂Ψ ∂α π = b π • α π (18) (19)

Flow rules

The energy rate Y d i that drives damage is defined by:

Y d i = 1 2 • K p (η) • u i 2 ( 20
Y 0 = 1 2 • K • δ d 2 (21)
The damage evolution law is chosen similarly to the one of Moutoussamy [START_REF] Moutoussamy | Essais hybrides en temps réel sur structures de Génie Civil[END_REF]:

d i = d ∞ •       1 -       Y 0 Y d i       q       (22)
with d ∞ ≤ 1 the maximum potential damage (i.e. giving the secant stiffness

K 0 • (1 -d ∞ )
for an infinite displacement) and q a coefficient driving the slope right after the damage initiation in the force-displacement curve.

For the friction phenomena and the kinematic hardening, the threshold surface is:

f π = |F π -X π | (23)
However, to take into account hysteresis loops that are nonlinear by nature, a non-associative potential of dissipation is used, as proposed by Armstrong-Frederick [START_REF] Frederick | A mathematical representation of the multiaxial Bauschinger effect[END_REF]:

Φ π = |F π -X π | + 1 2 • a π • (X π ) 2 (24)
Then, the flow rules are:

uπ = λπ • ∂Φ π ∂F π = λπ • sign(F π -X π ) (25)
and:

απ = -λπ • ∂Φ π ∂X π = λπ • sign(F π -X π ) ( 26
)
λπ is the plastic multiplier.

Model identification strategy

The model identification is performed thanks to a minimization algorithm applied on the global relative force error ε F given in equation [START_REF] Kumar | Cyclic Response of Sand Using Stress Controlled Cyclic Triaxial Tests[END_REF].

ε F = N i=1 F i -F exp i 2 F exp i 2 ( 27
)
F i is the total force. There are 7 parameters to identify, as summarized in table 2. To facilitate the identification process, a sequential procedure is followed: The effects of the different parameters are graphically described in figures 11 and 12.

Identification results

The procedure described section 4.2.6 has been carried out on the different beams tested with QSC1. The results obtained for the beam HA12-C1B-1 are summarized in table 2 and the corresponding error indicator of equation 27 equals ε F = 0.0376. The behavior of the identified model is illustrated in figure 13a in comparison with QSC1 experimental data. The quality of the identification validates the choices made for the macro-element (RC beam) model. The total force (F e + F π ) is plotted in figure 13a.

Numerical study of the EVDR dependence on engineering demand parameters

Engineering demand parameters

The dependence of the EDVR on the degradation and loading levels are investigated. The choice is made in this paper to select the mid-span curvature γ (see figure 14) as the quantity of interest to quantify the loading level. Also, a degradation index Γ is defined as the maximum historic curvature measured at mid-span γ m over the theoretical first steel yielding curvature γ y as expressed in equation [START_REF] Li | Field measurements of amplitude-dependent damping in a 79-storey tall building and its effects on the structural dynamic responses[END_REF].

The curvature has the advantage over the displacement to be bending-type independent (e.g. 3-point or 4-point bending) since it is defined at the most critical section. Moreover, γ y can be computed in a predictive way on the basis of design rules (such as Eurocode 2 [START_REF]Eurocode 2: Design of concrete structures -Part 1-1: General rules and rules for buildings[END_REF]) knowing the material and geometrical properties of the RC section. By using the assumption of nonevolution of the mode shape, one can compute directly by using the analytical mode shape of a dynamically bending beam, the coefficient k r of the curvature mid-span displacement law in equation [START_REF] Lagaros | Evaluation of ASCE-41, ATC-40 and N2 static pushover methods based on optimally designed buildings[END_REF].

γ = k r • δ (28)
with

k r = π 2 L 2 ≈ 0.28 m -2 (29)
For the computation of steel-yielding curvature, a Eurocode 2-type constitutive law for the concrete is considered:

σ c =                0 ε c < 0 f c • 1 -1 -ε c ε c2 2 0 ≤ ε c ≤ ε c2 f c ε c2 < ε c (30)
where ε c2 and f c are the (nonlinear) elastic limits in strain and stress depicted in figure 15. The value of steel-yielding curvature obtained the beam is given in equation [START_REF] Liu | Observed natural frequencies, damping ratios, and mode shapes of vibration of a 30-story building excited by a major earthquake and Typhoon[END_REF]. It is recalled that the dissipations associated to steel yielding are not in the scope of this study. Hence, Γ is expected to remain within the range [0, 1].

Γ = γ m γ y (31)
with

γ y = 0.0298 m -1 (32)
The study of the influence of the degradation level and the loading on the EVDR is challenging for pragmatic reasons. A full cyclic quasi-static test as QSC1 takes about 4 hours to perform (not counting the setup phase). Then, the testing for many degradation states could be very time-consuming. When practical constraints (for example related to the schedule, the cost, or the available equipment) make impossible to investigate experimentally the influence of different parameters on a quantity of interest, numerical experimentation represents an elegant way to address the problem, under the condition that the numerical model used has been validated experimentally. Considering the ability of the model identified in section 4.3 to describe the phenomena taking place during the bending, it can be used as support of virtual quasi-static testing.

δ m δ 1 δ 2 δ 3 . . . δ N (cm) (cm) (cm) (cm) (cm) 1 0.2 0.

Numerical loading

As a side positive consequence for the use of a virtual experimental, there is no spurious dissipation due to external sources. Hence, the validity of the representation of hysteretic dissipation by an equivalent viscous damper is assessed. To study the influence of the degradation level Γ and the prescribed curvature amplitude of the cycle γ, an increasing full cyclic quasi-static displacement loading has been designed and is illustrated in figure 16. Each cycle of displacement amplitude δ is repeated twice so the hysteretic behavior in the 2 nd cycle can be considered as stabilized and the energy dissipation due to damage initiation is not taken into account. It is important to note that the first two cycles are equal to the last two in order to remain at the same degradation index Γ of the beam all along the loading.

Observations resulting from the numerical simulations

It is clear from the results given in figures 17a,17b,17c that the degradation index Γ and the displacement amplitude play a major role in the EVDR. Moreover, these two parameters are model-independent. The model includes damage mechanisms, hence, the maximum elastic energy storage decreases with the degradation index Γ (figure 17b). The other important observation is that the energy dissipated does not seem to depend on Γ for the present nonlinear model. However, since the equivalent viscous damping ratio depends on the ratio of the dissipated energy over the stored one, the EVDR seems to increase with respect to the degradation index (see figure 17c). One can guess a relatively smooth surface fitting the EVDR points. Then, the highest EVDR point is chosen as the "identification" point of coordinates (Γ id , γ id , ξ id ) and is used for the next steps:

• γ = γ id is fixed at this point, and the best-fitting function f (Γ) is identified (thanks to a nonlinear least squares algorithm implemented in the curve fitting toolbox of Matlab) as in figure 18a:

f (Γ) = 3.679 • Γ Γ + 0.2806 (33)
then

f (Γ) = 0 % if Γ = 0 3.679 % if Γ → +∞ (34)
• Γ = Γ id is fixed at this point, and the best-fitting function g(γ) is identified as in figure 18b:

g(γ) = 0.03759 • γ γ 2 -0.004994 • γ + 8.397 • 10 -5 (35) then g(γ) = 0 % if γ = 0 0 % if γ → +∞ (36)
• the surface given by equation [START_REF] Priestley | Myths and fallacies in earthquake engineering, revisited[END_REF] where α = 1 ξ id is a normalization coefficient:

ξ = α • f (Γ) • g(γ) (37)
The search for an equation of surface in the form of a product of rational functions is motivated by the form of the expression of the EVDR evaluated by Jacobsen's areas method as a ratio between two energies. The degrees of these polynomial functions are deduced from the shape of the curved guessed from the red dots in figures 18a and 18b.

The values obtained for the EVDR depending on both Γ and displacement amplitude γ are represented with the red dots in figure 19a and19b. In figure 19a, the blue arrow illustrates the loading path actually performed experimentally (curvature amplitude and degradation index increasing together). This is also the loading path used for the identification of the nonlinear oscillator model (figure 8). Each point in the figure 19a (red and black points) corresponds to a different and independent quasi-static numerical computation and not to an extrapolation of the black points. Finally, the black points are, among all the tests performed (red dots), those used to identify the functions f (Γ) and g(γ) (identification illustrated in figures 18a and 18b).

In figure 19b, the surface obtained with the ξ(Γ, γ) (equation [START_REF] Priestley | Myths and fallacies in earthquake engineering, revisited[END_REF]) is compared to the numerical experiments performed (red and black dots). The so-called "damping surface" fits very well the numerical experiments as seen in figure 20. Indeed, only the black dots in figures 19a, 19b have been used for the calibration. The set of red points corresponds to additional numerical experiments considering different loadings for the same beam.

Even though this paper focuses on one particular beam (HA12-C1B-1), several beam designs have been tested as explained in section 2.1. The model proposed in section 4.2 is able to reproduce the behavior of these different designs because the

Conclusions

In order to improve the relevancy of simplified nonlinear time-history analysis by enhancing the additional viscous damping term, an experimental campaign has been carried out on reinforced concrete beams in both quasi-static and dynamic tests. In particular, a RC beam has been subject to a quasi-static full cyclic four-point bending test. Then, a single-degree-of-freedom model has been proposed and identified to reproduce the experimental restoring force for a prescribed modal displacement measured thanks to a dedicated image correlation based technique. This model has been used to perform additional numerical tests in order to assess both the influence of a degradation index Γ and the amplitude of the curvature at mid-span γ due to the displacement prescribed by the actuators. It is observed that both parameters have a major influence on the equivalent viscous damping. More specifically, the equivalent viscous damping ratio (EVDR) increases with degradation index Γ in an asymptotically way according to the model identified in figure 19b. However, the evolution of the EVDR as a function of the curvature amplitude is not monotonous according to the same model identified in figure 19b: a maximum is reached before the EVDR decreases asymptotically to zero. A function resulting from the product of the two previously identified ones sums-up these observations. It corresponds to the equation of a smooth surface that can be used as a virtual abacus for a more predictive choice of EVDR. It has been shown with two different RC beams (differences regarding concrete and geometry of the reinforcing steel bars) that the evolution of EVDR stays similar according to degradation state of the structure and the displacement demand. Nevertheless, these local quantities (quality of concrete, steel rebars diameter, etc.) affect the coefficients involved in the EVDR surface equation. A perspective of this work is the definition of this link by considering additional experimental and numerical tests.

At this point of this study, it seems important to take into account the variations of EVDR during the nonlinear time-history analysis since this value can evolve from almost 0 % to 3 % for the studied RC beam. For instance, this evolution of the EVDR can lead to a reduction of the local pseudo-acceleration spectrum and as a consequence, of the loading observed by the structure up to 10%. It is worthy noting that the steel yielding is not taken into account in this model, while it is a highly dissipating phenomenon. Many studies focusing on the dissipations related to steel yielding during earthquakes are found in the literature, often proposing an EVDR evolving with respect to the displacement ductility. Hence, the work carried out in this paper is intended to be complementary to these models. The next step of this study will be to use the damping surface as a basis for a fast and physical-wise updating strategy for the EVDR during nonlinear time-history analyses.

Figure 1 :

 1 Figure 1: Regulatory seismic hazard map in France (1. very weak [PGA = 0.4 m.s -2 ], 2. weak [PGA = 0.7 m.s -2 ], 3. moderate [PGA = 1.1 m.s -2 ], 4. average [PGA = 1.6 m.s -2 ], 5. strong [PGA = 3 m.s -2 ])[START_REF]Plan Séisme[END_REF]

3 (

 3 a) Longitudinal reinforcements for HA12-C1B-1 beam (b) Cross-section reinforcements for HA12-C1B-1 beam (c) Longitudinal reinforcements for HA20-C1A-1 beam (d) Cross-section reinforcements for HA20-C1A-1 beam

Figure 2 :Figure 3 :Figure 4 :

 234 Figure 2: RC beams reinforcements (lengths in centimeters)

Figure 5 :

 5 Figure 5: Notations used for the considered mass-spring-damper oscillating system

Figure 6 :

 6 Figure 6: Dependence between EVDR and excitation frequency, example for an EVDR estimated with an excitation frequency of 3 Hz

Figure 8 :

 8 Figure 8: Imposed actuators' displacement and resulting modal displacement of the experimental QSC1 loading

Figure 9 :

 9 Figure 9: Force-displacement (first mode displacement) measurement for the considered QSC1 test

Figure 10 :

 10 Figure 10: Illustration of the increase of moment of inertia due to crack closure under bending

Figure 11 :

 11 Figure 11: Plot of the pinching function over displacement for d i = 1 (0 value corresponds to completely closed cracks)

Figure 12 :

 12 Figure 12: Description of some parameters influence on the model behavior

Figure 13 :Figure 14 :

 1314 Figure 13: Comparison between model and experiments

Figure 15 :

 15 Figure 15: Concrete constitutive law from Eurocode 2 [6]

Figure 16 : 5 . 2 .

 1652 Figure 16: Loading procedure for the sensitivity study (example for Γ = 0.091)

Figure 17 :

 17 Figure 17: Influence of the degradation index γ over energies and equivalent viscous damping ratio ξ eq for different cycle amplitudes obtained by Jacobsen's areas method

1 0Figure 18 :

 118 Figure 18: Uncoupled identification of dependence between EVDR, degradation index and curvature amplitude

  EVDR values from numerical simulations (b) Same as figure 19a with superposed fitted surface

Figure 19 : 35 Fitted

 1935 Figure 19: EVDR values for an increasing degradation index Γ and curvature amplitude γξ id represents the identification point for uncoupled functions plotted in figures 18a and 18b; the black dots correspond to the data used for the functions identification of equations 33 and 35

Figure 20 :

 20 Figure 20: Comparison of the EVDR estimated from the numerical simulations and the corresponding values on the best-fitting surface (interpolated on the identification domain)

Figure 21 :Figure 22 :

 2122 Figure 21: Uncoupled identification of dependence between EVDR, degradation index and curvature amplitude for beam HA20-C1A-1

Table 1 :

 1 Internal variables of the model

		variable Internal variable
	Displacement	u
	Damage	d i
	Friction	u π
	Kinematic hardening	α π
	Crack closure	η

Table 2 :

 2 ) Model parameters and identified values Y 0 is an initial energy threshold:

	Par.	Description	Value	Unit
	δ d K 0 p	Displacement at first crack in concrete Initial stiffness Stiffness loss coefficient (=1-d ∞ )	6.4 1.88 0.228	mm N•mm -1 -
	q	Fragility coefficient	0.426	-
	a π b π	Hysteresis loops width Initial stiffness of the hysteresis loops	1480 74.7	N N•mm -1
	U c	Crack closure displacement	54.2	mm

Table 3 :

 3 Loading procedure with increasing degradation indexδ m : maximum displacement; δ i : amplitude of i th block

			018	0.2
	2	0.3	0.027	0.2	0.3
	3	0.4	0.036	0.2	0.3	0.4
	. . .	. . .	. . .	. . .	. . .	. . .	. . .
	N	10.0	0.91	0.2	0.3	0.4	. . . 10.0

Acknowledgments

The authors wish to express their most grateful thanks to CEA/DEN for its financial support. The work carried out under the SINAPS@ project has benefited from French funding managed by the National Research Agency under the program Future Investments (SINAPS@ reference No. ANR-11-RSNR-0022). The work reported in this paper has also been supported by the SEISM Institute (http://www.institut-seism.fr).


Huan Du

Thi-Bich-Thuy Tran

Patrick Perré
email: patrick.perre@centralesupelec.fr

A 3-variable PDE model for predicting fungal growth derived from microscopic mechanisms

Keywords: Reaction-diffusion, continuum modeling, upscaling, fungus

In this work, we present a new PDE model of the growth of Postia placenta, a species of brown rot fungus. The formulation was derived mainly from the biological mechanisms embedded in our discrete model, validated against experimental data. In order to mimic the growth mechanisms, we propose a new reaction-diffusion formulation, based on three variables: the concentration of tips, the branch density and the total hyphal density. The evolution of tips obeys a reaction-diffusion model, with constant diffusivity, while the evolution of the two other variables results from time integrals.

The numerical solution is in excellent agreement with the averaged radial tip/hyphal densities of the mycelial network obtained by the discrete model. Thanks to the efficient exponential Euler method with Krylov subspace approximation, the solution needs only 3.5 s of CPU time to simulate 104-day of mycelium growth, in comparison with 8 hours for the discrete model. The great reduction of the RAM memory and computing time gives the possibility to upscale the simulation. The novelty of the PDE system is that the spatial colonization is formulated as a diffusion mechanism, which is self-standing, contrary to models based on an advection term. The continuous model can also reproduce the radial densities when the growth parameters in the discrete model are varied to adapt to different growth conditions. The correlation constructed between the two models provides us a tool for mutual insights between local biological mechanisms to the global biomass distribution, especially when analyzing experimental data.

Introduction

Wood and wood-based products have been popular constructive materials due to their excellent structural properties for many years. The good thermal performance and design of new manufacturing technology enable the reduction of energy consumption and emission, which addresses the present needs in building construction. This is why these materials became even more promoted in construction. However, fungal decay in bio-based materials is of major concern in relation to the service time of buildings and even to human safety, since it weakens the structural support of wood enough to cause mechanical failure. Every year an enormous amount of wood and woodbased products is destroyed by fungal decay, among which brown rot fungi are the most common and destructive within buildings throughout Europe and North America due to their rapid decay mechanisms [START_REF] Schmidt | Wood and Tree Fungi[END_REF]. Many researches have studied the growth and decay patterns of brown rot fungi as well as the impact of environmental factors to their growth [START_REF] Jin | Structural characterization of brown-rotted lignin[END_REF][START_REF] Eriksson | Microbial and Enzymatic Degradation of Wood and Wood Components[END_REF][START_REF] Goodell | Low molecular weight chelators and phenolic compounds isolated from wood decay fungi and their role in the fungal biodegradation of wood[END_REF][START_REF] Yelle | Multidimensional NMR analysis reveals truncated lignin structures in wood decayed by the brown rot basidiomycete Postia placenta[END_REF][START_REF] Xie | Microcalorimetric characterization of the recovery of a brownrot fungus after exposures to high and low temperature, oxygen depletion, and drying[END_REF][START_REF] Gonzalez | Effects of environmental factors on decay rates of selected white-and brown-rot fungi[END_REF][START_REF] Wadsö | The activity of rot fungi (Postia placenta) during drying and rewetting cycles measured by isothermal calorimetry[END_REF]. Nevertheless, the laboratory observations are constrained by the scales of study and their cost. As a supplementary tool, mathematical modeling in combination with laboratory experiments can realize a deeper and wider insight in a more efficient way.

The modeling techniques to simulate the spatial distribution of fungi are classified into two categories: discrete models and continuous models. The previous ones adapt well to model the mycelial growth at a small scale in identifying individual hypha to form mycelial networks in different environments. For example, [START_REF] Boswell | The development of fungal networks in complex environments[END_REF] and [START_REF] Boswell | Modelling mycelial networks in structured environments[END_REF] respectively developed a two-dimensional (2D) and a three-dimensional (3D) lattice-based model to simulate mycelial growth in a soil-like environment [START_REF] Boswell | The development of fungal networks in complex environments[END_REF][START_REF] Boswell | Modelling mycelial networks in structured environments[END_REF]. [START_REF] Fuhr | Modelling the hyphal growth of the wooddecay fungus Physisporinus vitreus[END_REF] constructed a 3D lattice-free model of hyphal growth in the heartwood of Norway spruce in simplifying the structure of wood [START_REF] Fuhr | Modelling the hyphal growth of the wooddecay fungus Physisporinus vitreus[END_REF]. These models are generally constrained in a small scale by the necessary RAM memory and the computing time.

In contrast, the continuous models provide a potential of modeling the fungal growth from the colony-scale to the macroscale, in which the density of fungal matter was viewed as averaged quantities. The early models were extended from reaction-advection systems which included nutrients and geometries of mycelial spread [START_REF] Edelstein | The propagation of fungal colonies: a model for tissue growth[END_REF][START_REF] Edelstein | A model for Fungal Colony Growth Applied to Sclerotium rolfsii[END_REF][START_REF] Edelstein | Growth and metabolism in mycelial fungi[END_REF][START_REF] Edelstein-Keshet | Models for branching networks in two dimensions[END_REF]. However, the advection process requires the development direction to be known at each point, which limits the prediction ability of the code and/or augments the mathematical complexity especially in heterogeneous media. Thus, in the recent models, the diffusion process has replaced the advection to model the movement of hyphal biomass [START_REF] Davidson | Modelling the qualitative response of fungal mycelia to heterogeneous environments[END_REF][START_REF] Davidson | A mathematical model for fungal development in heterogeneous environments[END_REF][START_REF] Davidson | Translocation induced outgrowth of fungi in nutrient-free environments[END_REF][START_REF] Falconer | Biomass recycling and the origin of phenotype in fungal mycelia[END_REF][START_REF] Falconer | Modelling interactions in fungi[END_REF][START_REF] Cunniffe | Scaling from mycelial growth to infection dynamics: a reaction diffusion approach[END_REF] or combined with the advection process to model the movement of tips and hyphae [START_REF] Boswell | Functional consequences of nutrient translocation in mycelial fungi[END_REF][START_REF] Boswell | Growth and function of fungal mycelia in heterogeneous environments[END_REF]. However, in these models, some of the mechanisms were based on assumptions and no precise validation process has been carried out.

The strategy of our work is to enlarge the spatio-temporal scale of the study of fungal growth step by step in an accurate and efficient way. Thus, we started from the micro-scale experiment [START_REF] Du | Morphological Characterization and Quantification of the Mycelial Growth of the Brown-Rot Fungus Postia placenta for Modeling Purposes[END_REF] to observe the hyphal growth behaviors of a species of the brown rot fungi, Postia placenta. Then, a discrete model [START_REF] Du | A lattice-based system for modeling fungal mycelial growth in complex environments[END_REF], derived from the biological mechanisms of hyphal development, was precisely validated with experimental data. Compared to the experiment, this model allows us to predict the mycelial growth at a larger spatial-scale in much shorter time. However, to reach the macroscale with an even a higher computational efficiency, a continuous model is needed.

In this work, we propose a PDE formulation mainly derived from the mechanisms of hyphal growth at the microscopic level. A reaction-diffusion equation mimics the movement and proliferation of tips. Then, the evolution of hyphal biomass is related to the tip density following the tip elongation mechanism. Most of parameters of the continuous model were directly obtained from the discrete model, while three remaining parameters were determined by an inverse procedure: the model parameters were adjusted to minimize the objective function, i.e. the gap between the PDE profiles and the radial profiles of hyphal density calculated by the validated discrete model. Thanks to the exponential Euler method with Krylov subspace approximation, the model was accurately solved in an efficient way. Then, the growth parameters in the discrete model have been varied to investigate their influences to the distributions of tips and hyphae. The good agreements between the numerical solutions and the profiles obtained in the discrete model verified the excellent prediction capacity of the continuous model and its facility to adapt to different growth conditions.

Additionally, thanks to the good correspondence of the two models, the local mechanisms of the mycelial network can also be estimated via the change of biomass distribution. Furthermore, due to the pure diffusion process, this model is simple enough to extend to study the fungal growth in multi-dimensional spaces and in porous media, and even provides a base to the multi-scale modeling of the fungal development in wood.

Method and materials

Discrete mechanisms and radial profiles

In our discrete model [START_REF] Du | A lattice-based system for modeling fungal mycelial growth in complex environments[END_REF], three main microscopic mechanisms (i.e., tip elongation, branching and anastomosis) of hyphal growth were implemented in a regular square-lattice to simulate the development of a mycelial network. The two fundamental features of mycelial growth are tips and hyphae. Tips belong to two classes: active tips which extend at a certain, non-null, rate and dormant tips which do not elongate. A hypha elongates by the extension of an active tip, following its dominant direction. The elongation rate of each active tip follows a corrected gamma distribution. Each site of the lattice can contain numerous hyphae and/or tips. Lateral and apical branching are considered as way to produce new tips. The rate of lateral branching depends on the lateral branching probability and the total number of tips and the lateral branches are uniformly distributed throughout the mycelium. The apical branching occurs to the active tips and its rate is determine by the apical branching probability. The branching is inhibited in three special cases: (i) no lateral branches emerge in the proximity of an existed branch; (ii) the formation of a lateral branch around a hyphal tip is inhibited by apical dominance; (iii) the apical branching does not occur if its trailing hypha has not a certain length. Apart from the process of tip proliferation, tips can also disappear through anastomosis. The fusion of a tip with a hypha can occur when they contact each other to yield an interconnected mycelial network. This model has been validated by a precise process using the experimental data of the growth of P. Placenta during 17 days. With the determined parameters, the spatio-temporal development of a mycelial network was generated from a spore during a longer timescale of 9 × 10 6 s (∼ 104 days).

We divided the mycelial network into a series of rings with a gap of ∆r = 1.5×10 -4 m , of which the radial densities of tips and hyphae were respectively calculated by averaging the discrete network over radius increments. The tip density is defined as the number of tips per unit volume, while the hyphal density denotes the length of hyphae per unit volume. Due to the numerous stochastic processes in the discrete model, a single realization contains fluctuations. Therefore, we repeated the simulations 20 times with the same parameters and averaged the two resulting densities for each ring: (i) the averaged tip density denoted by ρ tip (r, t) in the ring around the radius r and at time t and (ii) the averaged hyphal density denoted by ρ hypha (r, t). In the following, the temporal evolution of the tip and hyphal distributions, the two averaged densities, were calculated at nine instants from t = 2.25 × 10 5 to t = 9 × 10 6 s.

Continuous model

A continuous model was directly derived from the three microscopic mechanisms of mycelial growth of the above-presented discrete model. The fungal mycelium is represented by two components: tip density and hyphal density, which represent respectively the proliferation capacity of the mycelium and the total biomass quantity. The tip density, denoted by C(x, t), is the number of tips in a unit volume at the spatial position x and time t, while B(x, t) is the hyphal density defined as the length of hyphae per unit volume.

The three main mechanisms of mycelial growth, tip elongation, branching and anastomosis, lead to respectively the migration, the proliferation and the elimination of tips. Thus, as shown in Eq. (1a), a reaction-diffusion equation was applied to describe the variation of the tip density.

The diffusion term denotes the migration of tips with a diffusion coefficient D. The two reaction terms, S and P , represent the proliferation and the elimination of tips. As the hyphal extension is generated only by tip elongation, the amount of hyphae created in a unit volume over a short time interval is linearly related to the actual quantity of tips by two factors, the averaged tip elongation rate R and the proportion of active tips P active (Eq. (1b)).

∂C ∂t = ∇(D∇C) + S(B, C) -P (B, C), (1a)
∂B ∂t = RP active C. (1b)
As mentioned above, the term S denotes the production of tips by lateral and apical branching. According to the branching mechanisms in our discrete model, the number of lateral branches emerged per unit of time is the product of lateral branching probability P lat br by the total quantity of tips in the mycelium. The latter is obtained as the integral of the tip density over all the entire volume (i.e., Cdx). The lateral branches should be uniformly distributed in the mycelium. Thus, the amount of branches distributed into each unit volume depends on the proportion of hyphae containing in this volume over the whole mycelium (i.e., B/ Bdx). In addition, the lateral branching is inhibited in two special cases related to the number of branches and tips per unit hyphal length.

As shown in Eq. ( 2), we applied an exponential function to account for this inhibition, where β is defined as the inhibition coefficient and Br(x, t) denotes the number of branches per unit volume at position x and time t. Since the branches do not move, Br is the accumulation of the term S over time, calculated by Eq. (3).

Y (x, t) =      e -β(C(x,t)+Br(x,t))/B(x,t) B(x, t) > 0, 0 B(x, t) = 0. ( 2
)
∂Br ∂t = S. (3)
For apical branching, each active tip possesses a probability P api br to form a new tip. The increase of tips due to the apical branching is hence given by P api br P active C. Considering that the apical branching occurs when the hypha attains a certain length, we modified C(x, t) with a timedelayed term C(x, t -τ ) where τ is the necessary time for a tip to reach that given length with the average elongation rate. Together, the proliferation term S is modeled as:

S(B, C) = P lat br Cdx B Bdx Y + P api br P active C(x, t -τ ). ( 4
)
The disappearance of tips is caused by anastomosis. Such a fusion occurs when a hyphal tip contacts another hypha, which implies that the disappearance rate depends on the contact probability between tips and hyphae. Since tips elongate at a low rate (∼ 10 -9 m • s -1 ), they move a very short distance during a unit time interval. We focused on one tip and we took the volume formed by a unit surface and the short distance the tip passes through during a unit time interval.

All the hyphae contained in this volume were projected on the bottom surface. The possibility of this tip contacting another hypha is equal to the proportion of the projected area, which is related to the hyphal length. In all, the elimination rate of tips is proportional to the tip density C, the tip elongation rate R, and the hyphal density B and can be expressed as follows:

P (B, C) = γRBC. (5)
The coefficient γ is named as the anastomosis coefficient for correlating the projected area with the hyphal density.

To summarize, the reaction-diffusion system takes the following form Tip density

∂C ∂t = ∇(D∇C) + S -P, (6)
Hyphal density

∂B ∂t = RP active C, (7)
Branch density

∂Br ∂t = S. (8)
where

S =      P lat br
Cdx B Bdx e -β(C+Br)/B + P api br P active C(x, t -τ ) B(x, t) > 0,

P api br P active C(x, t -τ ) B(x, t) = 0, (9a)
P = γRBC. (9b)

Numerical solution

The model was solved in a one-dimensional (1D) axisymmetric configuration to be compared with the radial profiles obtained with the discrete model. The 1D axisymmetric formulation reads as follows:

∂C(r, t) ∂t = 1 r ∂ ∂r (Dr ∂C(r, t) ∂r ) + S -P, ∂B(r, t) ∂t = RP active C, ∂Br(r, t) ∂t = S. (10)
The term C(x, t -τ ) in Eq. ( 9a) is approximately transformed to an term able to be solved:

C(t -τ ) ≈ C(t) -τ ∂C ∂t (t). (11)
The initial condition for this model was defined as the averaged density profiles of tips, hyphae and branches obtained with the discrete model at t = 2.25 × 10 5 s.

The reaction-diffusion system was solved in cylindrical coordinates with a computational domain large enough to assume no-flux boundary conditions. A control-volume finite-element spatial discretization turned the PDEs system (Eq. ( 10)) into an initial value problem of the form

du dt = g(u), u(0) = u 0 , (12)
where u = (C, B, Br) ∈ R 3N , g: R 3N ⊂ D → R 3N is a nonlinear function of u, and N is the size of the 1D mesh. The calculation domain was meshed uniformly with a spacing of 1.5 × 10 -4 m which is equal to the spatial resolution of the profiles obtained in the discrete model.

We performed the time integration process to get a linearized version of Eq. ( 12),

du dt = g(u n ) + J g (u n )(u -u n ), t > t n , (13)
whose exact solution is

u(t) = u n + (t -t n )ϕ((t -t n )J g (u n ))g(u n ), t > t n , (14)
where J g ∈ R 3N ×3N is the Jacobian matrix of g and function ϕ is defined as ϕ(z) = e z -1 z .

Applying the exponential Euler method to Eq. ( 14) at time t n+1 , we get the approximate solution of Eq. ( 12)

u n+1 = u n + dt n ϕ(dt n J g (u n ))g(u n ), (15)
with time step dt n = t n+1 -t n .

The exponential integrators method was first presented in 1960 by Certaine [START_REF] Certaine | The solution of ordinary differential equations with large time constants[END_REF] and has been well-known since late of 1990s by works of M. Hochbruck and others [START_REF] Hochbruck | Exponential integrators for large systems of differential equations[END_REF][START_REF] Hochbruck | Exponential integrators[END_REF]. Thanks to the approximation of the Krylov subspace to the matrix-function vector product ϕ(dt n J g (u n ))g(u n ), the exponential integrators method became very useful for large system of stiff equations [START_REF] Carr | Efficient simulation of unsaturated flow using exponential time integration[END_REF][START_REF] Carr | A variable-stepsize Jacobian-free exponential integrator for simulating transport in heterogeneous porous media : application to wood drying[END_REF][START_REF] Michels | Exponential integrators for stiff elastodynamic problems[END_REF].

In this article, our work was inspired from the work of E. Carr, I.Turner and P. Perré [START_REF] Carr | A variable-stepsize Jacobian-free exponential integrator for simulating transport in heterogeneous porous media : application to wood drying[END_REF] using the "variable-stepsize exponential Euler method" to solve the system of non-linear transport equations with great success.

The Krylov subspace methods compute a small m×m matrix H m (the projection matrix) instead of the full Jacobian matrix J g . The approximation of the matrix-vector product is as follows:

ϕ(dt n J g (u n ))g(u n ) ≈ β 0 V m ϕ(dt n H m )e 1 . (16)
Thanks to this, we can have a large time step dt n without preconditioning techniques, which is very efficient in terms of computational time and RAM memory.

Parameter determination

The continuous model was calibrated for P. Placenta via the validated discrete model mentioned in Part 2.1. All parameters are summarized in table 1 using the international unit standard. Most of these parameters were defined directly from the discrete model [START_REF] Du | A lattice-based system for modeling fungal mycelial growth in complex environments[END_REF]. For example, since the tip elongation rates follow a corrected gamma distribution, its mean was calculated as the averaged tip elongation rate R. With this averaged rate, we obtained τ ≈ 1048 s, which is the average time for a tip to extend by 5 µm. The proportion of active tips, the lateral and the apical branching probability are simply those of the discrete model.

The determination of the three remaining, macroscopic parameters, the diffusion coefficient D, the inhibition coefficient β and the anastomosis coefficient γ, were determined by inverse analysis.

The objective function used for optimization compares, in the sense of mean squared difference, the radial profiles obtained with the discrete and continuous models. Since the tip and the hyphal density are tightly related and our final objective is to fit the latter one, we defined the relative error E(B) of the hyphal density profiles:

E(B) = t∈{t1,...,t9} N k=1 (ρ hypha [r(k), t] -B[r(k), t]) 2 t∈{t1,...,t9} N k=1 ρ 2 hypha [r(k), t] , (17)
where the total number of radius increments is denoted by N (here N = 299) and r(k) represents the radius of the k th ring. This function was minimized using the Nelder-Mead method. * Parameter values come from [START_REF] Du | A lattice-based system for modeling fungal mycelial growth in complex environments[END_REF].

Results

Radial profiles calculated from the discrete model

The spatial distribution of tips and the mycelial network were simulated during a long period (9 × 10 6 s, ca. 100 days) by the discrete model. Fig. 1 shows the evolution of the tip distribution and the mycelial network at three elapsed times: t = 1.8 × 10 6 , 4.5 × 10 6 and 9 × 10 6 s. These profiles show very clearly how the tips proliferate and colonize the free space starting from the inoculum zone. After a long-term growth, the distribution of tips seems uniform throughout the mycelium except at the edge. Clearly, the area of the colony increases with time, but we can also see the increase of the hyphal density inside the colony from Fig. 1 (c) which shows the zoom of a tiny part of the mycelial network.

Then, nine profiles of the averaged radial density at different times were obtained from the discrete model (Fig. 2). The profiles confirm that, after a long-term growth, the tip density is almost uniformly distributed, even though it decreases slightly towards the center of the colony.

The higher distribution of tips at the edge is a strategy that favors the exploitation of new space.

The uniform distribution implies the existence of an equilibrium between the emergence of branches and the anastomosis. In spite of the wave-like shape of the tip density profiles, its higher value at the edge maintains a local gradient. It seems therefore reasonable to propose a diffusion process to simulate the tip movement in a continuum approach. One may note also that the tip density in the center increases first and then decreases at longer times. Indeed, when the mycelial network becomes dense, the anastomosis and the inhibition occur more frequently, leading to more elimina-tion of tips. From Fig. 2 (b), we observed that the hyphal density profiles continuously increased in time keeping a negative gradient towards the edge. Note however that the hyphal density increases more and more slowly in the center due to the reduction of tip density, revealing the existence of an upper limit to the biomass even for the free growth of mycelia on a two-dimensional substrate.

Validation of the PDE formulation

The profiles averaged over 20 realizations of the stochastic model (Fig. 2) were used to optimize the three macroscopic parameters of the continuous model. The optimized solution of Eq. ( 10) has a very small relative error (E(B) = 2.70 × 10 -2 ), which proves the ability of the PDE formulation to represent the discrete mechanisms. The profiles, plotted at different times in Fig. 3, give another piece of information: even though the objective function considers only the hyphal density, both the tip and the hyphal profiles of the PDE solution are in excellent agreement with the discrete model. This confirms that the two variables, tip density and hyphal density, are tightly related to each other. Moreover, the wave-shape of the tip profiles was well captured by the diffusion process, even though it resembles more a front-like progression than a diffusion phenomenon. Compared to the advection process, which requires the propagation direction to be known a priori, the diffusion process greatly facilitates multi-dimensional modeling. Indeed, as the propagation results are from a diffusion process, the propagation direction is defined locally from the actual gradient of the tip concentration. The other discrete mechanisms are also well described by the continuous terms, even the implicit ones. This is the case in particular for the anastomosis rate at the macroscale, assumed to be linearly related to the tip density, the tip elongation rate and the hyphal density, and for the branching inhibition expressed as an exponential function depending on the branching and tip density per unit length of hyphae. The value of the three parameters (i.e., D, β and γ) determined by inverse procedure are listed in Table 1.

Ability of the PDE model to simulate different growth conditions

In the previous section, we proved that the continuous model perfectly reproduces, at the macroscopic scale, the growth of P. Placenta for the reference parameter set of the discrete model ( [START_REF] Du | Morphological Characterization and Quantification of the Mycelial Growth of the Brown-Rot Fungus Postia placenta for Modeling Purposes[END_REF][START_REF] Du | A lattice-based system for modeling fungal mycelial growth in complex environments[END_REF]). However, the change of the growth conditions, such as temperature, humidity or nutrient concentration, influences the behavior of hyphal growth, particularly branching and elongation rate [START_REF] Meyer | Reconstruction of signaling networks regulating fungal morphogenesis by transcriptomics[END_REF][START_REF] Rousk | Growth of saprotrophic fungi and bacteria in soil[END_REF][START_REF] Jacobs | Solubilization of calcium phosphate as a consequence of carbon translocation by Rhizoctonia solani[END_REF][START_REF] Dighton | Fungi in Ecosystem Processes[END_REF]. It is therefore very important to check whether the continuous model is able to predict the growth behaviour for different parameters of the discrete model. To that purpose, we altered the three key parameters of hyphal growth in the discrete model, (i.e., the apical branching probability P api br , the lateral branching probability P lat br and the tip extension rate R). The values of the three parameters listed in Table 1 were considered as reference. Then, for each parameter, we performed simulations with 0.5, 1.5 and 2 times the reference value. These factors were arbitrarily chosen to represent a quite significant change of growth conditions. This means that nine series of simulations (20 realizations per case study) have been carried out by changing one parameter of the discrete model at a time. The corresponding profiles were obtained in the same way as described in Part 2.1.

Obviously, these changes were reported when these discrete parameters are involved explicitly in the continuous formulation. The three macroscopic parameters identified in the previous section were used as initial guess. Then, for each series of profiles, we optimized four times the objective function E(B) (Eq. ( 17)) by using different set of free parameters (D, γ, β), (D, γ), (D, β) and D.

The relative errors of the hyphal density between the discrete and continuous solutions are presented in Fig. 4 before and after optimization. According to this bar chart, a good continuous-discrete match requires the value of D to be modified when the discrete parameter P api br changed, while both D and γ should be adjusted when parameters P lat br and R changed. Accordingly, parameter β can be kept unchanged in any case. The relative errors of the optimal solutions are all around 4.00 × 10 -2 , which is low enough to ensure the capacity of the continuous model to replicate the mycelial growth under different conditions.

To better understand the impact of the three discrete parameters on the mycelial growth, we compared the profiles of the tip and the hyphal density obtained with different values of P api br , P lat br and R (respectively Fig. 5, 6 and 7) after optimization. In accordance with the results of Fig. 4,D was the unique adjusted parameter to match the changes of P api br while two parameters, D and γ, were adjusted when changing the value of P lat br or R. The numerical solutions agree well with the discrete results for all of the nine cases. With the increase of the apical and lateral branching probability, the local tip and hyphal density increase more rapidly since the proliferation rate of tips is higher. As a subtler effect, we may notice that the extension rate of the colony also augments with the increase of local tips. This result is from a statistical effect: with the increase of tips at the edge of the mycelium, the chance to have tips whose extension direction is close to the radial direction also increases, resulting in a more efficient spatial extension of the network. With the increase of the tip elongation rate, the extension rate of the colony rises but the local hyphal density reduces significantly. Since the tip proliferation rate does not change, the amount of tips produce by branching is similar, but, as the colony occupies a larger space, less tips will be distributed per unit volume. To summarize, even though the impact of P lat br and R is greater than that of P api br , the three parameters influence the morphology of the mycelial network by intricate effects. This necessarily Consistently, the tip elongation rate R has a direct impact on D. The effect of R is almost linear with a slope equal to 1.5. The fact that D varies more rapidly than R, can be explained by the indirect effect of tip elongation rate on tip density.

The effect of branching probabilities on D reveals more subtle mechanisms. Let's focus first on the negative correlation between apical branching probability P api br and D. We have to keep in mind that a higher value of P api br slightly increases the colony expansion rate, because the probability to have active tips right along the macroscopic radial direction increases with the tip density. Besides, the tip density at the colony edge is higher, giving rise to a larger tip gradient. Consequently, the diffusivity D must be smaller to obtain the right expansion rate. This coupled effect of tip gradient and diffusivity to build up the macroscopic extension rate is also likely to explain the positive correlation between lateral branching probability P lat br and D. Due to apical inhibition, a higher lateral branching probability increases the tip density everywhere except right at the colony edge. This non-uniform tip source term tends to decrease the tip gradient at the colony edge and should be compensated by a slightly larger diffusivity D.

Regarding the effect of growth condition on γ, equation 9b reminds us that this factor corrects the linear relation between biomass density B and anastomosis. Indeed, due to the statistical partition of biomass B in a volume, the projected area along tip extension is not simply proportional to B. In a certain volume, the proportion of overlapping hyphae increases with density. When P lat br increases or R decreases, the hyphal density increases and γ should be reduced to obtain the correct anastomosis term.

Discussion

As summary of this work, it is worth to come back to our global upscaling strategy of mycelial growth modeling. We first proposed a discrete model which mimics the mycelium development at the microscopic scale [START_REF] Du | A lattice-based system for modeling fungal mycelial growth in complex environments[END_REF]. This model accounts for the fundamental mechanisms of network development (tip extension, lateral and apical branching, anastomosis). It has been calibrated and validated from experimental data by a rigorous procedure. For that purpose, the network development has been observed during 17 days by confocal microscopy, and the image were subsequently processed to extract relevant information [START_REF] Du | Morphological Characterization and Quantification of the Mycelial Growth of the Brown-Rot Fungus Postia placenta for Modeling Purposes[END_REF]. In the discrete model, some parameters such as the lateral and apical branching probabilities, were determined by an inverse analysis. Without the discrete model, these parameters could not have been accurately determined from experimental data. Moreover, via the validation process, some underlying mechanisms have been discovered, such as the inhibition of branching which is almost impossible to observe directly from the mycelial morphology.

Acting as a bridge between the microscopic behaviors of hyphal growth and the mycelial development at the macroscale, the validated discrete model was used to simulate the mycelial growth during ∼ 104 days. This is quite computationally intensive as one realization requires about 50 Go of RAM memory and a computational time of ca. 8 hours on a single core of a HPC cluster with the processor Intel Xeon CPU E5-2670 v3 @ 2.30 GHz. By this way, the spatial dimension was extended from 0.9 cm × 0.9 cm (experiment) to 9 cm × 9 cm (simulation). Therefore, the discrete model already allowed the spatio-temporal scale to be efficiently extended. 20 realizations of the discrete model, ca. 160 hours of CPU time, were averaged to obtain tip/hyphal density profiles. This piece of information served as input data for the continuous model proposed in the present work. Based on the hyphal growth mechanisms at the discrete due to the decrease of the tip density in the center at increasing hyphal density, it seems that there exists a limit for the biomass density on the 2D substrate. This last part told us that only two key parameters (the diffusion coefficient D and the anastomosis coefficient γ) out of the three macroscopic parameters, need to be adapted when the growth conditions change. It is certain that efforts should be made to figure out the quantitative change of the growth parameters depending on environmental factors, and this work is in progress in our team. However, the good quality of the continuous-discrete match proves the potential of the continuous model to describe the fungal growth under different conditions. It also provides a tool to realize the bidirectional transfer, from the local biological mechanisms to the global biomass distribution and vice-versa. For example, if the local mechanisms of hyphal growth change (different growing conditions or another species of interest), only their local growth characteristics are needed to be measured to obtain both the mycelial morphology and the biomass distribution. On the contrary, the variation of the local mechanism and the mycelial network can also be estimated via the experimental observation of the biomass distribution together with an inverse procedure. In this sense, the numerical efficiency of the continuous model is of particular interest.

Conclusion

This work proposes a new three-variable reaction-diffusion PDE model capable of simulating and predicting the mycelial growth under different environmental conditions. This model was calibrated using variable profiles predicted by the average over several realizations of a stochastic discrete model. For source and sink terms of tips, biomass production and branching density, the formulation rigorously follows the discrete mechanisms of mycelial development, which implies that all parameters have a physical meaning. Thanks to the discrete model, the ability of this new model to account for various growth conditions was successfully tested. Moreover, the correlation between the discrete and continuous model has been constructed, which provides us a tool to realize the bidirectional transfer, from the local biological mechanisms to the global biomass distribution and vice-versa.

As the mycelium network propagation is formulated as a diffusion process, this model has the potential to be easily extended in multi-dimensions and to heterogeneous media such as a Low Density Fiberboard. In the future, additional features such as translocation and substrate-degradation mechanisms will be incorporated in the model to simulate the mycelial growth in the real morphology of bio-based materials as obtained by nano-tomography. Once combined with decomposition mechanisms, this model can also be applied to solid-state fermentation.

Fig. 1 .

 1 Fig. 1. Simulation of the discrete model: temporal evolution of (a) the distribution of tips, (b) the mycelial network and (c) the zoom of the red squared part in (b).

Fig. 2 .

 2 Fig. 2. Radial densities of tips and hyphae obtained by averaging respectively the tip distributions and the mycelial networks over radius increments at nine times. Profiles on the right are the average of 20 realizations of the stochastic discrete model.

Fig. 3 .

 3 Fig. 3. Comparison of the numerical solutions (solid line) after optimization with the averaged radial profiles (dashed line) of the discrete model at eight times. The black solid lines depict the initial condition for the continuous model, which are the averaged profiles obtained with the discrete model at t = 2.25 × 10 5 s.

Fig. 4 .

 4 Fig. 4. Relative errors of the hyphal density (E(B)) between the discrete profiles and the numerical solutions in nine cases (in each case we varied one growth parameter in the discrete model). The five colors represent the numerical solutions before the optimization (blue) and after the optimization executed in specifying D, β and γ (orange), D and β (grey), D and γ (green), and D (yellow).

Fig. 5 .

 5 Fig. 5. Comparison of the numerical solutions after the optimization in identifying D and the radial profiles in the three cases of the variation of the apical branching probability P api br in the discrete model.

Fig. 6 .

 6 Fig. 6. Comparison of the numerical solutions after the optimization in identifying (D, γ) and the radial profiles in the three cases of the variation of the lateral branching probability P lat br in the discrete model.

Fig. 7 .

 7 Fig. 7. Comparison of the numerical solutions after the optimization in identifying (D, γ) and the radial profiles in the three cases of the variation of tip elongation rate R in the discrete model.

Fig. 8 .

 8 Fig. 8. Relationship between D or γ and the three growth parameters. X ref represents the reference value of each parameter listed in Table 1. Dots are the optimal values of D or γ obtained for different growth parameters. Dash lines are the smoothing functions.

Table 1 .

 1 Description of variables and parameter values in the continuous model.

	Symbol	Description	Unit	Value
	C	Tip density	m -3	--
	B	Hyphal density	m • m -3	--
	Br	Branch density	m -3	--
	t	Time	s	--
	R *	Tip elongation rate	m • s -1	4, 77 × 10 -9
	P * active	Proportion of active tips	--	0, 7
	P lat * br	Lateral branching probability	s -1	1.11 × 10 -5
	P api * br	Apical branching probability	s -1	3.22 × 10 -6
	τ *	Lag time for apical branching	s	1048
	D	Diffusion coefficient	m 2 • s -1	1.13 × 10 -12
	β	Inhibition coefficient	--	7.78 × 10 -5
	γ	Anastomosis coefficient	m	5.40 × 10 -2

Acknowledgments

This work has benefited from the financial support of the LabeX LaSIPS (ANR-10-LABX-0040-LaSIPS) managed by the French National Research Agency under the "Investissements d'avenir" program (n°ANR-11-IDEX-0003-02). The calculus presented in this work were performed on the Fusion HPC Center, CentraleSupélec, Université Paris-Saclay and on the ROMEO HPC Center, Université de Reims Champagne-Ardenne. Financial supports from Grand Reims, the Marne department and the Grand Est region are gratefully acknowledged.

scale, a new reaction-diffusion system, involving three spatio-temporal variables (tip, mycelium and branch densities), has been derived. While most discrete parameters are involved explicitly in the continuous formulation, three macroscopic parameters had to be defined by inverse analysis.

Once done, the discrete profiles are perfectly reproduced, in space and in time, by the PDE model.

Thanks to the exponential Euler method with Krylov subspace approximation, the computational time to solve the continuous model was again greatly compressed to less than 4 s on a single core of the same cluster. In brief, by following the upscaling strategy, the continuous model realizes the prediction of the mycelial growth during 104 days in less than 4 s. In addition, since the memory required by the continuous model is much less than the discrete model, it is possible to further enlarge, by several orders of magnitude, the spatial and time scales.

Our continuous model is based on a reaction-diffusion formulation. Unlike advection models, it is not necessary for the diffusion process to know the development direction at each point a priori. This feature opens the door to multi-dimensional modeling and growth simulation in porous media. For the reaction part, we rigorously followed the mechanisms of the discrete model, which implies that all parameters have a physical meaning. In addition, thanks to the balance between diffusion and production (a global effect of branching), an almost constant gradient establishes at the edge of the network. This can be seen at the edge of the colony in Fig. 3(a): a tip density peak develops towards the free space. Therefore, an almost constant propagation rate arises from this diffusion-reaction mechanism. Consequently, a constant diffusivity is likely to simulate a given average tip elongation rate. The reaction-diffusion process is applied to the tip density, while the hyphal density is the accumulation of the tip density over time. Similarly, the branch density at any location is obtained as the time integration of the source terms of branches at that location.

The 3-variable PDE model involves coupling between variables. For example, the tip density is connected to the hyphal density as the latter influences the lateral branching distribution and the anastomosis. Even so, the relation between the two main densities (tips and biomass) is nicely reproduced by the continuous model.

Apart from the reference simulation, additional series of profiles were produced by the discrete model with various parameter values, which are artificially set just to give us insight into the link between the possible microscopic mechanisms and the correponding macroscopic behaviours. These configurations have been reproduced by the continuous model. Whether the branching rate or the elongation rate has great impact on the morphology of the mycelium by changing the colony extension rate and the local biomass distributions. They can also alter the maximal value of the tip density, revealing that they change the balance between branching and anastomosis. In addition,