Alexandre Malley-Ernewein

Sylvie Lorente
email: sylvie.lorente@villanova.edu

The shape of the elemental system in a porous medium designed for Thermochemical Energy Storage

Keywords: thermochemical energy storage, constructal, porous reactor

Here we document how the design of the elemental component of a designed reactor for thermochemical energy storage can be improved for better energy performance and less pumping power. The pore channels are made of cylinders with walls are covered by salt reacting with the water vapor molecules of the blown humid air. The design of such open systems is changed to dendritic configurations while keeping the salt and fluid volumes constant. Morphing the configuration from I-shaped to Y-shaped architectures allows to increase the heat exchanges.

Additional improvements are obtained when the geometric features of the dendritic configurations are chosen to follow constructal design.

Nomenclature

Introduction

Tree-shaped networks are an icon of constructal design [START_REF] Bejan | Constructal tree network for fluid flow between a finite-size volume and one source or sink[END_REF][START_REF] Bejan | Constructal tree-shaped flow structures[END_REF][START_REF] Wechsatol | Tree-shaped insulated designs for the uniform distribution of hot water over an area[END_REF][START_REF] Zimparov | Constructal tree-shaped parallel flow heat exchangers[END_REF][START_REF] Ledezma | Constructal tree networks for heat transfer[END_REF]. Dendritic shapes are ubiquitous both in nature (river basins, blood vessels networks) and in engineering. When the objective is to connect one point to an infinite volume of points at constant fluid volume and with the smallest effort possible, it was demonstrated that the most performant architecture corresponds to a treeshaped network [START_REF] Wechsatol | Optimal tree-shaped networks for fluid flow in a discshaped body[END_REF]. The diameter ratio between parent tube and children tubes depends on the type of flow (laminar or turbulent) and on the number of children tubes [START_REF] Bejan | Design with Contructal Theory[END_REF]. This ratio can be predicted, but its value changes with the driving force. For example when the latter is a concentration difference, the best configuration in a Y-shaped assembly made of cylinders is D2/D1 = 2 -1/2 [START_REF] Wattez | From pore network prediction based on the Constructal law to macroscopic properties of porous media[END_REF], whereas it is D2/D1 = 2 -1/3 [START_REF] Bejan | Thermodynamic optimization of geometry: T-and Yshaped constructs of fluid streams[END_REF], when a pressure difference drives the laminar movement of fluid. Yet, if the flow network is embedded within a solid which emits heat uniformly in order to extract that heat, we showed that the combination of 2 objectivesminimum pumping power and maximum heat transfer -leads to architectures that may not remain dendritic, depending on the flow characteristic parameters (Be number, Re number, level of bifurcation) [START_REF] Wechsatol | Dendritic heat convection on a disc[END_REF]. An exhaustive effort was implemented recently in Refs. [START_REF] Miguel | Constructal branching design for fluid flow and heat transfer[END_REF][START_REF] Cetkin | Constructal branched micromixers with enhanced mixing efficiency: Slender design, sphere mixing chamber and obstacles[END_REF] for Newtonian and non-Newtonian fluids, and for constant and pulsatile flow systems. In this work, the fluid flows through cylindrical ducts which walls are receiving a constant and uniform heat flux. The analysis accounts for the pressure losses due to friction but also for the acceleration of the fluid due to the temperature increase associated to the heat flux along each wall. Therefore the diameter ratio leading to smaller pressure losses and hence smaller pumping power depends on the strength of the heat flux.

Constructal design finds numerous applications in the domain of energy storage. The essence of heating processes and their link to the constructal law was covered recently in [START_REF] Bejan | Thermodynamics of heating[END_REF][START_REF] Lorente | Current trends in constructal law and evolutionary design[END_REF]. Latent energy storage designs were presented in [START_REF] Lorente | Constructal design of latent thermal energy storage with vertical spiral heaters[END_REF][START_REF] Ziaei | Constructal design for convection melting of a phase change body[END_REF][START_REF] Joshi | Constructal enhancement of thermal transport in metal foam-PCM composite-assisted latent heat thermal energy storage system[END_REF][START_REF] Nazzi Ehms | Theoretical and numerical analysis on phase change materials (PCM): A case study of the solidification process of erythritol in spheres[END_REF] with the common objective of showing how the energy storage performances can be improved, while Refs. [START_REF] Solé | Geometry optimization of a heat storage system for concentrated solar power plants (CSP)[END_REF][START_REF] Xie | Optimization of Pin-Fins for a Heat Exchanger by Entropy Generation Minimization and Constructal Law[END_REF] proved the superiority of constructal designs in the design of fins for heat exchangers. Thermochemical Energy Storage which consists in storing and releasing heat by means of the chemical reaction occurring between a reactant and a gas [START_REF] Lougou | Heat transfer and fluid flow analysis of porous medium solar thermochemical reactor with quartz glass cover[END_REF][START_REF] Lehmann | Modelling sorption equilibria and kinetics in numerical simulations of dynamic sorption experiments in packed beds of salt/zeolite composites for thermochemical energy storage[END_REF] is receiving increasing attention in building applications mainly because heat leakages are absent from such reactions. Our previous work [START_REF] Malley-Ernewein | Constructal design of thermochemical energy storage[END_REF][START_REF] Malley-Ernewein | Analysis of thermochemical energy storage in an elemental configuration[END_REF] demonstrated the relevance of constructal design for Thermochemical Energy Storage.

In this work, we document the merit of a dendritic network, in its elemental configuration of one pairing level, vis-à-vis a straight channel for Thermochemical Energy Storage. The networks are submitted to several driving forces: a pressure difference drives the flow of humid air, a concentration difference generates a flow of mass (represented here by water vapor), while the flow of heat comes from the non-uniform and unsteady heat released or absorbed by the salt chemical reaction on the walls of the channels. The reacting module is a designed porous material which pores allow the heat transfer fluid (humid air) to cross from one side to the opposite side (Fig. 1). Because the pore walls are lined with reacting salt, the water molecules transported by the fluid flow can diffuse through the very thin salt layer (grain scale) and react with the salt. Cold humid air entering the material allows the hydration of the salt, and less humid and warmer air leaves the elemental channel through the opposite face of the porous material. In a dehydration case, hot and dry air enters the duct, the concentration difference in water molecules (under the appropriate conditions of pressure and temperature) leads to the salt dehydration and to the storage of the incoming heat.

The shape of the pores is a degree of freedom. Known are the constant mass flow rate crossing the module, the constant salt volume available and the constant fluid volume, which means that ultimately the module porosity is fixed. In this work we consider the pore as an elemental two-dimensional system. The pore could be a straight slot connecting directly one face of the module to the opposite face, or it could have a dendritic shape.

Configuration and model

Assume a rectilinear channel of hydraulic diameter Dh = 2D which walls are covered with a thin salt layer of thickness e (Fig. 2). The length of the channel is L. The heat transfer fluid is blown through the channel while water vapor reacts with the salt thanks to a molecular diffusion process.

While the fluid volume and the salt volume are kept constant, the shape of the element is changed from an I to a Y. The fluid enters at a mass flow rate m in the mother channel, then splits into 2 daughter tubes at a mass flow rate m /2. The mother channel has a hydraulic diameter D = 2D and a length L1, while the daughter channels have a hydraulic diameter D = 2D2 and a length L2 (Fig. 2). We have

DL = D L + 2D L (1)
eL = e L + 2e L

where e, e1 or e2 are the salt thicknesses.

Such elemental configuration is aimed at being a component of a designed vascularized reactor.

At first assumption the outlets (and inlets) are equidistant as in [START_REF] Lorente | Heterogeneous porous media as multiscale structures for maximum flow access[END_REF], and the distance between the two outlets is termed d. Hence, the Y-shaped configurations exhibit two degrees of freedom: the hydraulic diameter ratio and the length ratio (or angle of bifurcation α). and its corresponding flow resistance is given by:

∆p m ~ L D (3)
The Y-shaped channel flow resistance is:

∆p m ~ L D + L 2D (4)
The latter is minimum for a constant fluid volume when the Di ratio is D /D = 2 / . This result is obtained by invoking the method of the Lagrange multipliers. The aggregate function Φ is constructed from the flow resistance and the constant channels volume (assuming a unitary thickness).

Φ = L D + L 2D + λ D L + 2D L (5)
D /D = 2 / is the value for which ∂ϕ / ∂D = ∂ϕ / ∂D = 0.

Such result is meaningful provided the Svelteness number is greater than 10. The Svelteness number is a non-dimensional macroscopic parameter that allows at a glance to determine if and when the local pressure losses can be neglected when compared to the friction losses [START_REF] Lorente | Svelteness, freedom to morph, and constructal multi-scale flow structures[END_REF]. Sv = # 10 is the minimum value for which such assumption is acceptable. Here, we define Sv as L + L / D L + 2D L / . In every configuration considered, Sv > 10.

Finally the pressure drop becomes:

∆p m ~ 1 D %L + 2 / L & (6)
or

∆p m ~ L + L D Sv (7)
The problem is described thanks to mass conservation for the humid air in the channel (Eq. ( 8)), mass conservation for the water vapor in the channel and in the salt layer (Eqs. ( 9) and ( 10)) and momentum conservation for laminar flow in the channel (Eq. ( 11)). The first law of thermodynamics is written along the tube and through the salt (Eqs. ( 12) and ( 13)). Finally the reaction kinetics equation gives the reaction advancement rate for the hydration (Eq. ( 14)) and for the dehydration (Eq. ( 15)).

Dρ ( Dt = 0 (8) 10 / 29
where ρ ( is the humid air density. It is assumed that the changes in water vapor content do not impact the humid air properties [START_REF] Michel | Experimental investigation of an open thermochemical process operating with a hydrate salt for thermal storage of solar energy: Local reactive bed evolution[END_REF].

Dc Dt + +%,D -,( +c & = 0 ( 9
)
where c is the water vapor concentration, and D -,( is the diffusion coefficient of vapor in the fluid.

∂c ∂t + +%,D -, +c & = , γ n ∂a ∂t (10)
where D -, is the diffusion coefficient of vapor in the salt, γ is the reaction stochiometric coefficient, ns is the salt bulk molar concentration and a is the reaction advancement.

D2 Dt = , 1 ρ ( +p + ν+ 2 (11)
where 2 is the velocity vector, p is the fluid pressure and ν is the kinematic viscosity.

ρ ( c ( DT Dt + + • k ( +T = 0 (12)
where c ( is the fluid heat capacity, T is the temperature and kf is the fluid thermal conductivity.

c ∂T ∂t + + • k +T = n ∂a ∂t ∆h ( 13
)
where c is the salt heat capacity, ks is the salt thermal conductivity and Δh is the reaction heat.

∂a ∂t = k 89: ;1 , p < p - 1 , a (14)
∂a ∂t = k 89: ;1 , p < p - a (15)
where k 89: is the reaction kinetic coefficient, peq is the equilibrium pressure, and pv is the vapor partial pressure.

The boundary and initial conditions are the following: at the entrance of the fluid channel, the water vapor concentration, c 9: , and the temperature, T 9: , are fixed. An entrance length greater than 0.05 Re > ? × D is fixed to ensure that the flow is fully developed [START_REF] Langhaar | Steady flow in the transition length of a straight tube[END_REF]. The inlet mass flow rate m is constant. The outlet boundary conditions are atmospheric pressure and thermal outflow.

The walls of the channels are thermally insulated without any vapor flux. A no-slip condition is imposed at the interface between fluid and salt.

The initial temperature of the system is Tinit (Tinit < Tin), the fluid is at rest and at atmospheric pressure, the initial vapor concentration is cinit corresponding to the equilibrium vapor pressure at temperature Tinit (cinit < cin). Table 1 gives the set of data. The salt characteristics are the ones of the SrBr2 [START_REF] Michel | Experimental investigation of an open thermochemical process operating with a hydrate salt for thermal storage of solar energy: Local reactive bed evolution[END_REF].

Table 1: Data used in the numerical simulations, adapted from [START_REF] Malley-Ernewein | Analysis of thermochemical energy storage in an elemental configuration[END_REF][START_REF] Michel | Experimental investigation of an open thermochemical process operating with a hydrate salt for thermal storage of solar energy: Local reactive bed evolution[END_REF]. The set of equations is solved by means of a FEM code [30]. Validation was made by comparing the numerical results to experimental ones in a different configuration. The details were carefully presented in a previous paper [START_REF] Malley-Ernewein | Constructal design of thermochemical energy storage[END_REF]. Rectangular elements of constant thickness are placed along the salt thickness, the elements shape in the fluid is rectangular except in the vicinity of the bifurcation (see Fig. 3). The strategy chosen was based on the verification of both mass and energy conservation. Mass conservation was found much less sensitive to the mess fineness than energy conservation. The mesh was refined regularly from one test to the other one, until the relative difference in the result of energy conservation was below 1%. Based on this assessment,

Initial and boundary conditions

Y or I-shaped architecture

First, the thickness of the module was allowed to decrease when moving from the straight pore to the dendritic one (see Fig. 4). We also chose D = D while the fluid and salt volumes remained fixed. Under such assumptions, the friction losses corresponding to the Y-shaped configuration become

∆p m ~ 1 D %L + 2 / L & (16)
And because the amount of fluid is constant,

∆p m I = ∆p m (17)
In other words, we work first at identical fluid mechanics performance and look for thermal improvements. The aim of the Y-shaped design is to increase the size of the interface between the heat transfer fluid and the reacting salt. The salt thickness is assumed to remain the same along the mother and the 2 daughter tubes. Its value is obtained from the definition of the constant salt volume (Eq. ( 2)). Table 2 and Fig. 4 summarize the configurations obtained. The last column is the relative difference between the heat transfer exchange surface (fluid/salt interface) in the Y-shaped arrangement and the I-shaped configuration. Remember that the mass flow rate, fluid volume and salt volume are constant. The diameter ratio is always D2/D1 = 2 -1/2 in accord with the constructal law. For example, when the bifurcation angle α = 20°, the dendritic configuration corresponds to a length ratio of 0.48. Morphing the geometry from a straight channel to a dendritic one, leads to an overall porous material thickness -from one side to the opposite side of the corresponding porous reactor -of 0.88 L (Fig. 4).

We show in Fig. 5 the reaction advancement as a function of time, obtained in the different configurations of Fig. 4. Time is plotted on the abscissa in a non-dimensional way. The nondimensional values are related to the moment when the reaction advancement is 1 for the slowest case. The I-shaped configuration is the one needing more time to reach full salt hydration. At the start of the process, the reaction advancements are identical, whatever the channels shape; this is true until the non-dimensional time reaches a value of 0.02. This result is expected as the entrance dimension of each configuration is constant. Then, in time, the salt hydration is faster in the case of the Y-shaped architectures which in the case α = 20° reaches full hydration for a nondimensional time of 0.79. In this configuration, the change in geometry comes with an increase of the heat transfer surface of 17% as shown in Table 2, facilitating the diffusion of water molecules through the salt layer and therefore the chemical reaction. The results are confirmed in Fig. 6 where the heat released is plotted as a function of the reaction advancement. The nondimensional values are based on the heat released in the I-shaped configuration. Even though the salt volume is the same, the change in design allows a maximum increase in heat released of 20% when the bifurcation angle α = 10°. Note that the improvement in performance lasts all along the salt hydration process, as the curve corresponding to the Y-shaped design is always 20% above the straight channel configuration in Fig. 6.

Decreasing the bifurcation angle allows to increase the interface between fluid and salt, which in turns leads to a much faster reaction advancement. According to Fig. 5, if the Y-shaped structure has an angle of bifurcation greater than 50°, then the salt hydration advancement is mainly identical to a I-shaped design. The same happens in terms of heat released by the chemical reaction. A small angle between the two daughter channels guaranties a higher heat transfer rate during the entire hydration process (see Fig. 6) as the relative impact of the daughter ducts increases. To get more insight, the full hydration-dehydration process was modeled. The energy performance of the configurations is assessed by dividing the heat transfer rate due to the chemical reaction by the power necessary to blow the fluid through the channels. This ratio is termed β. Next the ratio β I / Kβ LMN K was calculated. The non-dimensional time is based on the time necessary to complete a full hydration-dehydration cycle which corresponds to the straight channel case.

The evolution of the energy performance indicator β is presented in Fig. 7 for I-shaped and Y-shaped (α = 20°) configurations in the case of the full process. The results show that the full cycle for the Y-shaped configuration takes 20% less time than the I-shaped configuration. This feature is favorable for short term storage applications, composed by several hydrationdehydration cycles in a short period. In addition, the absolute value of the β indicator in the Yshaped configuration is 10% higher than for the I-shaped configuration. Note that the difference in the hydration and dehydration durations is due to the inlet conditions used here (summarized in Table 1). They are taken close to standard reactor operation conditions [START_REF] Michel | Experimental investigation of an innovative thermochemical process operating with a hydrate salt and moist air for thermal storage of solar energy: Global performance[END_REF]. 4 shows the configurations for which the entrance width remains identical to the I-shaped case while again the heat exchange surface is kept constant. Here D2/D1 varies and so does α as the spacing d between the two outlets does not change. Figure 8 is a summary of the results obtained with Tables 3 and4. The values presented correspond to a non-dimensional time of 0.005. The design performances can be improved in two ways under the imposed constraints (salt and fluid volume, fluid-salt heat exchange surface):

either by increasing the inlet thickness of the Y-shaped networks, or by increasing the daughtermother channels D /D ratio. Either solution is a path to less flow resistances as the fixed heat transfer performance is maintained constant with the heat exchange surface.

More constraining is the problem where the overall material dimensions are fixed. In such case, while the fluid and salt volumes remain unchanged, opportunities to improve both the heat transfer and fluid mechanics aspects exist. This is typically a case when the constructal law provides the directions on how the configurations should evolve. In Table 5, the bifurcation angle is constant because the exchange surface between fluid and salt and the spacing d between outlets are fixed when the ratio D /D changes. The constructal ratio D /D = 2 / leads to minimum pumping power, while the heat released is insensitive to the morphing of the channels. Therefore the corresponding energy performance indicator is also the highest in this case.

The search for continuous improvement requires to change the bifurcation angle which, based on the constraints, means that α can vary between 8° and 80°. The change in bifurcation angle allows to vary the heat transfer exchange surface as shown in Table 6. Figure 9 provides an overview of the results for a non-dimensional time of 0.005. For a fixed angle of bifurcation, the performance map reaches a maximum which corresponds to the daughter-mother channels ratio D /D obeying the constructal design methodology. Increasing the energy performance indicator requires to liberate the design configuration and to allow the bifurcation angle to increase which means shorter daughter channels and higher heat transfer exchange surfaces.

Conclusion

The objective of this work was to discover if a drastic change in flow architecture is a way to improve the thermal performance of an elemental configuration for thermochemical energy storage. From the reference case of a straight channel as basic configuration of a pore network part of a reactor, the shape of the pore was moved to a Y-shaped channel made of a main duct and 2 identical daughter channels. The fluid volume was kept constant. The fundamental result of the study is that the Y-shaped configuration leads to a faster reaction advancement and heat released than the I-shaped one.

When the overall thickness of the reactor is allowed to change, the improvements are noteworthy provided the dendritic architecture exhibits an angle of bifurcation smaller than α = 50°. In such cases, the morphing Y-shaped configurations permits to increase the heat transfer surface and therefore the thermal performance of the element.

The dramatic change in shape comes with a more compact pore network as the best performances are obtained when the angle of bifurcation decreased to α = 10°. This configuration corresponds to a reduction of 25% of the module thickness. A more compact structure leads to increases in the overall thermal and fluid mechanics results.

For a given set of constraints, we can predict how the shape of the pores should evolve in order to benefit the most of a given amount of salt. When maintaining the reactor thickness constant, the main findings are that:

-changing the diameter ratio to the one predicted by the constructal law is the only way to improve the energy performance indicator when the heat transfer exchange surface is fixed.

-the road to more improvements is to be found in allowing the bifurcation angle to increase which leads to more heat exchanges while maintaining the pressure as low as possible.

Finally, this work contributed to highlight the interdisciplinary and versatile nature of the constructal law approach to engineering problems.

Figure 1 :

 1 Figure 1: Two different pore configurations in a designed porous material: assemblies of straight channels or Y-shaped channels (not to scale).

Figure 2 :

 2 Figure 2: I-shaped and Y-shaped channels (not to scale). The total fluid volume is fixed, together with the salt volume and the mass flow rate.

  10 elements were located on the transversal direction of the salt. From the fluid/salt interface to the mid distance D/2, a sensitivity analysis proved that 15 elements are sufficient away from the bifurcation. They are very thin at the fluid/salt interface and increase in size when moving toward the center of the channel. The maximum thickness ratio is 5. The elements in the fluid flow direction are 4 000, equidistantly positioned. In the case of Y-shaped channels, the elements size both in the fluid and in the salt are even more refined. Finally, the number of elements is 116 000 for the I-shaped configuration, and from 100 000 to 250 000 depending on the value of the bifurcation angle in the Y-shaped designs. It takes about 2000 hours to obtain full salt hydration.For example full hydration happened after 1850 hours of exposure to humid air in the case of the Y-shaped channel with a bifurcation angle of 50°. A carefully study was implemented, reducing the time step until no change in results was noticed from one simulation to the other one. A time step of 1 hour appeared to fulfill this requirement.

Figure 3 :

 3 Figure 3: Example of mesh for a Y-shaped configuration. Global view of the bifurcation (left)and detail of the mesh at the junction between the mother duct and bifurcation.

Figure 4 :

 4 Figure 4: Overview of the different designs.

Figure 5 :

 5 Figure 5: The reaction advancement as a function of the non-dimensional time. Y-shaped configurations with different angles of bifurcation (α = 10°, 20°, 50° and 120°) are compared to I-shaped design.

Figure 6 :

 6 Figure 6: The non-dimensional heat released by the salt chemical reaction as a function of the advancement; Y-shaped configurations with different angles of bifurcation (α = 10°, 20°, 50° and 120°) are compared to a I-shaped design.

Figure 7 :

 7 Figure 7: Evolution of the energy performance during a hydration-dehydration cycle for a Ishaped and a Y-shaped configuration. In the latter, the bifurcation angle is α = 20°.

Figure 8 :

 8 Figure 8: Evolution of the energy performance indicator as a function of the ratio S /S for fixed entrance dimension S (blue plane) and as a function of the entrance dimension S for fixed diameter ratio S /S (green plane), when the reactor overall thickness is free to vary.

Figure 9 :

 9 Figure 9: Evolution of the energy performance indicator as a function of the D /D ratio for fixed bifurcation angle α (blue plane), and as a function of the bifurcation angle α for fixed ratio D /D (red plane), when the overall reactor thickness is fixed.

Table 2 :

 2 Y-shaped configurations for constant salt and fluid volumes, and constant inlet width D.

			Relative increase in heat
	α (degree)	L2 / L1	transfer exchange surface compared to the I-shaped
			design (%)
	120	0.06	3.3
	50	0.14	6.9
	20	0.48	16.8
	10	3.33	33.6

Table 3 :

 3 Y-shaped configurations with different entrance dimensions D for fixed salt-fluid exchange surface and fixed ratio D /D . The bifurcation angle α and length ratio LY/LI-shaped L I = L + L cos R are determined from the previous constraints.

	S/SI-shaped	D2/D1	D1/D	α	LY/LI-shaped
	1	2 -1/2	0.94	13.2°	0.62

Table 4 :

 4 Y-shaped configurations with different ratios D /D for fixed salt-fluid exchange surface and fixed entrance dimension. The variation in D /D allows to calculate the corresponding bifurcation angle α and length ratio LY/LI-shaped L I = L + L cos R .

	S/SI-shaped	D1/D	D2/D1	α	LY/LI-shaped
			0.83	11.3	0.56
			0.77	15.7	0.68
			2 -1/2	20.0	0.75
	1	1	0.67	22.8	0.78
			0.50	34.5	0.85
			0.40	41.7	0.87
			0.33	46.6	0.88
	Table

Table 5 :

 5 Y-shaped configurations with different ratios D /D for fixed salt-fluid exchange surface and fixed reactor overall thickness. The bifurcation angle α and the entrance dimensions D are determined from the previous constraints.

	S/SI-shaped	α	D2/D1	D1/D
			1.5	1.03
			1.0	1
	1	20
			2 -1/2	0.86
			0.5	0.57

Table 6 :

 6 Y-shaped configurations with different bifurcation angle α for fixed reactor overall thickness and fixed ratio D /D . The entrance dimensions D and salt-fluid exchange surfaces are determined from the previous constraints.

	D2/D1	α	D1/D	S/SI-shaped
		8	0.63	1.36
		10	0.69	1.24
		20	0.86	1
	2 -1/2
		40	0.97	0.88
		60	1.01	0.85
		80	1.03	0.83

Acknowledgements

Alexandre Malley-Ernewein is funded by the French ANR research project DECARTH, ANR-16-CE22-0006-03.


Pascal Dominique
email: p.dominique@brgm.fr

Hideo Aochi

Jacques Morel

Triggered Seismicity in a Flooded Former Coal Mining Basin (Gardanne Area, France)

Keywords: Post-mining seismicity, Mine aquifer, Seismic migration, Rock burst

We studied the correlation between seismicity and the water table level in an abandoned coal mine (Bouches-du-Rhône, France), closed in 2003, where groundwater has been pumped out since 2010 to prevent underground flooding. Microseismicity was first felt by the population in 2010 and a strongly felt seismic swarm occurred in November 2012. The origin of the seismicity was therefore questioned, in relation to both the potential instability of old, shallow galleries that might generate damage at the surface and a local seismic hazard assessment. A temporary dense seismic network in the area allowed us to analyse the spatial distribution of the seismicity in detail. Most of the seismicity was clearly located under the mine workings, highlighting that an existing fault system crossing the mining operation was being hydraulically activated, in accordance with the known tectonic extension regime. Our analysis clearly shows a spatiotemporal relationship between seismic migration and the level of the mine aquifer between 2013 and 2017. Thus, seismicity will persist with oscillations of the mining aquifer, depending on the pumping capacities and effective rainfall. Continuous hydraulic and seismic monitoring is necessary to better understand these phenomena and assess the associated risks.

Introduction

The Gardanne coalfield is located in Provence, southern France, between the towns of Marseille and Aix-en-Provence. Coal was exploited from the beginning of the seventeenth century until 2003. The first exploitation was performed along shallow seams directly from the surface. Then, room-and-pillar techniques were used down to a depth of 700 m, while the deepest exploitation, down to a depth of 1300 m, used longwall mining. Approximately 130 Mt of coal was extracted from the Gardanne coal mine.

The long historical experience of the mining industry has shown that underground works exploited with room-and-pillar techniques can lead to mechanical instabilities, especially when groundwater fluctuations occur. With time, pillars left in place deteriorate, and the process that can be accelerated in the presence of water, sometimes leading to the failure of one or several pillars, with potential damage to the surface [START_REF] Goldbach | What is the seismic risk of mine flooding? CSIR[END_REF][START_REF] Senfaute | Induced micro-seismicity and mechanical response during the experimental flooding of an iron ore mine[END_REF].

As coal exploitation ceased in 2003, the halt of pumping operations led to progressive flooding of the underground mine workings until the mine groundwater level quasistabilized in 2010, though it still fluctuates. To survey the risk of mechanical instabilities associated with room-andpillar sectors below inhabited areas, a permanent microseismic network was installed at the end of 2007, which has since been monitored by Ineris (French National Institute for Industrial Environment and Risks) on behalf of BRGM (French Geological Survey). This network consists of five microseismic stations located in the immediate vicinity (a few hundred meters) of these room-and-pillar sectors and seismic sensors positioned in boreholes at different depths. This network detects any microseismic event within a radius of 400 m around the stations and monitors potential microseismic activity induced by rock fracturing associated with underground mine workings deterioration. The intent is to identify any risks with potential consequences to the surface, and thus protect the population.

3

Since 2010 and especially November 2012, several seismic swarms lasting a few days have occurred outside these ground instability surveillance zones. The seismic vibrations regularly felt by the population in the Fuveau-Greasque (FG) area have led the administration and local authorities to request the deployment of temporary complementary surface seismic network to determine their origin, including the possibility that pre-existing faults were being reactivated.

Therefore, the first purpose of this study was to verify whether these seismic ground motions are the result of the collapse of old mining galleries (roof collapse) or related to ruptures along pre-existing tectonic faults. In both cases, the role of water and fluctuations of the mine aquifer level are important. The second purpose was to correlate the spatiotemporal evolution of the seismicity with the water level in the mine workings.

Geological and Hydrogeological Settings

Geological Setting

Western Provence is part of the Pyrenean-Provençal domain in southern France (Fig. 1). The tectonic formation of western Provence occurred at the end of the Eocene. It has four structural units: the Bandol, Beausset, Arc, and Durance units [START_REF] Rouire | Carte géologique de la France à 1/250 000[END_REF].

The Arc unit lies between the North Provençal and South Provençal overlaps. In the west, it disappears beneath the Rhodanian Quaternary to the west of the Salon-Cavaillon fault (SCF). In the east, it ends at the level of the Trias of Barjols, beyond which the tectonic style and the Mesozoic characteristics are different. The E-W faults of La Nerthe, L'Etoile, Sainte-Baume, and Allauch belong to the Arc unit. Therefore, the Arc unit extends to the south under the southern Provençal overlaps of Etoile, Nerthe, and Sainte-Baume, while to the north, the overlaps of Fare and Aix-Eguilles cover the unit.

Locally, the basin is limited to the south by the overlapping east-west running Diote fault (DF), which extends further to the Jean-Louis fault (JLF) in the east. The submeridian Meyreuil fault (MF) is an extension of the Durance fault from the north to south. It crosses the Gardanne Basin and affects both the indigenous lands of the Arc Basin and the overlapping Diote fault (Fig. 2).

Stress State

Based on seismicity data and focal mechanisms, large-scale studies near the outer boundaries of the western Alpine Belt highlight the strain/stress state of this western Provence area [START_REF] Baroux | Analyses of the stress field in southeastern France from earthquake focal mechanisms[END_REF][START_REF] Delacou | Present-day geodynamics in the bend of the western and central Alps as constrained by earthquake analysis[END_REF], suggesting an extensional stress (normal faulting) regime represented by a NNW-SSE trend for the horizontal principal stress axis, σ 1 . Geodetic data confirm this trend [START_REF] Calais | Current strain regime in the Western Alps from continuous global positioning system measurements, 1996-2001[END_REF], even if the strain rate is very low (less than 0.01 µstrain/yr), as well as the horizontal and vertical GPS velocities (less than 1 mm/yr) [START_REF] Devoti | A combined velocity field of the Mediterranean region[END_REF].

At the Gardanne mine site, [START_REF] Gaviglio | Consequences of structural change on the physical characteristics of fuvelian limestones in the Arc Basin (Bouchesdu-Rhône, France): fracturing, mechanical properties, stress states[END_REF] and [START_REF] Gaviglio | Measurements of natural stresses in a Provence mine (Southern France)[END_REF] obtained stresses at a very local scale from in situ measurements (flat jack method or hydraulic fracturing). They showed that two zones of stress state can be distinguished (Fig. 2): one isotropic northern extension zone and another anisotropic southern compression sector. The northern zone presents mainly as an extension to N45°E, particularly in the inflexion zone (with a vertical stress very close to the N135°E maximal horizontal stress), whereas the southern zone briefly corresponds with an anisotropic compression of E-W to NE-SW. A sub-meridian extensional stress in the northeastern part of the area has been mentioned [START_REF] Chalumeau | Hydrogeological functioning of lignite operations in the Fuveau basin (Bouches-du-Rhône)[END_REF][START_REF] Ch | Coal basin of Provence. Study of the effects on surface terrains of underground mining works and their foreseeable evolution over time[END_REF].

Hydrological and Hydrogeological Contexts

The above stress state studies were confirmed in [START_REF] Chalumeau | Hydrogeological functioning of lignite operations in the Fuveau basin (Bouches-du-Rhône)[END_REF] thesis, using the geomorphology method. Moreover, the author described this disconnection and its separate hydrogeological behaviour (Fig. 2). Indeed, due to the very strong inflow of water from the surrounding massifs

Diote Fault

Fig. 2 Gardanne mine area with stress state (after [START_REF] Gaviglio | Consequences of structural change on the physical characteristics of fuvelian limestones in the Arc Basin (Bouchesdu-Rhône, France): fracturing, mechanical properties, stress states[END_REF]Gaviglio et al. 1996, modified) separating the southwestern dry zone in compression and the northeastern aquifer zone in extension (after Chalumeau 2000, modified). Arrows represent magnitudes and orientations of the principal stresses determined for each site (Gaviglio et al. 2016). Large hollow arrows give the stress state (extension or compression) and its mean azimuth [START_REF] Gaviglio | Consequences of structural change on the physical characteristics of fuvelian limestones in the Arc Basin (Bouchesdu-Rhône, France): fracturing, mechanical properties, stress states[END_REF][START_REF] Chalumeau | Hydrogeological functioning of lignite operations in the Fuveau basin (Bouches-du-Rhône)[END_REF][START_REF] Ch | Coal basin of Provence. Study of the effects on surface terrains of underground mining works and their foreseeable evolution over time[END_REF] (Sainte-Victoire and Regagnas), the development of underground operations was particularly difficult in the northeastern part of the Gardanne mine area (Fuveau-Greasque (FG) area). Thus, many drainage galleries were excavated to limit the costs of pumping at the level of the shafts. Likewise, even the construction of several dams to block the galleries could not prevent the abandonment of some mining districts, such as the Rocher Bleu district. These significant inflows of water are linked to several factors:

• The presence of a vast natural impluvium corresponding to the Arc Basin collecting all rainwater from the limestone massifs of Sainte-Victoire to the north and northeast, from the Regagnas (southeast) and from l'Etoile (south). [START_REF] Dheilly | Management of the Gardanne mining reservoir, Bouches du Rhône: Environmental and hydraulic approaches to pumping and its discharge into the Mediterranean Sea. Past and present mining operations: environmental and societal impacts[END_REF] credited the origin of the main water inflows to the Gardanne mine to surface infiltration from the Fuvelian and Begudian outcrops and deep-water inflows (mainly Jurassic) from the sources of the Sainte Victoire at -220 m ASL (above sea level), Source 90 at -50 m ASL, and Massif de l'Etoile. [START_REF] Chalumeau | Hydrogeological functioning of lignite operations in the Fuveau basin (Bouches-du-Rhône)[END_REF].

The main level of mining, called the Grande Mine, is ≈ 2 m thick. Outcropping over the Greasque sector, the Grande Mine layer gradually deepens towards the northwest, reaching a depth of ≈ 330 m below the FG area. This deepening continues to the west until it reaches 1300 m west of Gardanne. The residual voids of the Grande Mine layer correspond to a privileged drain of all these mining aquifersmainly the Jurassic Aquifer-facilitating the circulation of water (volume, flow, transfer time) from East to West.

In contrast, compressed zones of lower permeability are found in the southwestern sector close to certain faults that sometimes act as barriers, such as the Jean-Louis fault to the south of the Rocher Bleu [START_REF] Chalumeau | Hydrogeological functioning of lignite operations in the Fuveau basin (Bouches-du-Rhône)[END_REF]. Likewise, the clayey Valdonian (upper Cretaceous; 72-84 My) and marl limestone Neocomian (lower Cretaceous; 129-145 My) impede the circulation of water between aquifers [START_REF] Gonzalez | Public information project on the Gardanne mining basin (Bouches-du-Rhône)[END_REF].

Management of Mine Water in the Gardanne Coal Mine

After mine closure, pumping of the Jurassic aquifer at a flow rate of 1500 to 2000 m 3 /h at the Arc shaft was shut down on May 3, 2003. This led to the natural rise of the groundwater level, with a progressive flooding of the mine workings starting from the deeper parts on the western side and progressing towards the east. To avoid overflow of mineralized water rich in iron and sulfate (due to sulfides in the oxidized lignite) contaminating the port of Marseille and flooding hazards at the surface, water flow has been again controlled by pumps set in a former Gerard mine shaft (see Fig. 2). As pumping resumed (August 4, 2010), the water level in the mine reservoir rose from -1100 m ASL to -14 m ASL. Between 15 and 45 Mm 3 of mine voids were saturated. Hydrological models [START_REF] Dheilly | Management of the Gardanne mining reservoir, Bouches du Rhône: Environmental and hydraulic approaches to pumping and its discharge into the Mediterranean Sea. Past and present mining operations: environmental and societal impacts[END_REF] integrating water level pumping capacities and effective rainfall show that a constant mine aquifer (as described above) of 800 m 3 /h supplies the mine reservoir. An additional volume of 200 to 1000 m 3 /h depends on the annual effective rainfall. All these volumes are pumped at the Gerard shaft in Gardanne (+ 18 m ASL) and discharged by a 14 km pipeline outside the Marseille harbour, at a depth of 30 m, in the Mediterranean Sea. Thanks to the pumping system, the mine groundwater is stabilized between -30 m ASL and + 10 m ASL, ≈ 250 to 350 m below the ground surface in the FG area. To prevent pollution and flooding, the mine water level must be maintained between + 18 m ASL (flooding level in the Galerie de la Mer) and -35 m ASL (pump level at -40 m ASL).

The mine water level is monitored remotely. A high-flow pumping test (3000 m 3 /h) was performed between mid-July and the end of September 2014 in the Arc shaft, located 12 km east of the Gerard shaft. The water in the Jurassic aquifer comes from rainfall at the heights of the Sainte-Victoire and Regagnas massifs. This pumping and the slope rupture of the water level recorded in the Gerard shaft established the connectivity of the aquifer between the two shafts with a delay of 1.5-2 days [START_REF] Dewandel | Research project ARK. Improved knowledge of the aquifer in the Aix-Gardanne basin[END_REF], and the strong influence of inflows from the Jurassic aquifer in the Arc Basin [START_REF] Dheilly | Management of the Gardanne mining reservoir, Bouches du Rhône: Environmental and hydraulic approaches to pumping and its discharge into the Mediterranean Sea. Past and present mining operations: environmental and societal impacts[END_REF].

Seismic Context

Seismicity Near the Gardanne Coal Mine

Over the past five centuries, several earthquakes have strongly shaken the Provence region. Historical data show a relatively sparse epicentral distribution of earthquakes (Sis-France 2016; [START_REF] Jomard | The SISFRANCE database of historical seismicity. State of the art and perspectives[END_REF]. Some earthquakes highlight regional faults, such as the Moyenne Durance, Nîmes, or Trevaresse Faults (Figs. 1,3a).

The Lambesc earthquake on June 11, 1909 [START_REF] Baroux | Source parameters of the 11 June 1909, Lambesc (southern France) earthquake: a reappraisal based on macroseismic, seismological and geodetic observations[END_REF] was the most destructive of the twentieth century in France. With an intensity of VIII-IX and a moment magnitude of Mw5.7 [START_REF] Baroux | Source parameters of the 11 June 1909, Lambesc (southern France) earthquake: a reappraisal based on macroseismic, seismological and geodetic observations[END_REF]Manchuel et al. 2017), it caused 46 deaths, 250 injuries, and the partial destruction 1 3 of some villages near the epicentre (Rognes, Saint-Cannat, and Lambesc) and part of the town of Salon-de-Provence (SisFrance 2016). Other notable tremors occurred in the Bouches-du-Rhône (Fig. 3) on May 28, 1985, (intensity V and Mw2.6), in Vernègues on December 23, 1934, in Salon-de-Provence (IV-V and Mw3.3) on December 12, 1846, in Gemenos (VI and Mw3.9), on March 25, 1783, and in Mallemort (V and Mw3.2) and Arles between 1725 and 1734. The oldest recorded earthquake dates from May 26, 1397 (V-VI and Mw3.5) in Arles. We also note the events of the Middle Durance (1509, 1708, 1812, 1863, and 1913), some of which reached epicentral intensity VIII.

Specifically, in the Gardanne coal mine area, the February 19, 1984 (VI and Mw3.6) Mimet earthquake has been studied regarding its origin: artificial or natural [START_REF] Haessler | The seismic event of February 19, 1984 in western Provence: tectonic earthquake or mining blow[END_REF]. Based on a spectral analysis of aftershocks, the authors showed that the nature of the signal was very close to that of a tectonic earthquake. Hence, regardless of the contradiction between the focal mechanism, the regional stress state, and the hypocentral location that favours an event associated with mine exploitation, they concluded that this February 19, 1984 shock was not caused by a mine collapse, but was likely due to a local modification of the regional stress field induced by mine exploitation.

Great effort has been made by the SI-Hex Working Group [START_REF] Cara | SI-Hex: a new catalogue of instrumental seismicity for metropolitan France[END_REF] to discriminate between natural and artificial earthquakes. Nevertheless, Fig. 3b shows a cluster of seismic events centred on the Gardanne coal mine. It is possible that this seismic activity is linked to natural activity or to artificial activity during the exploitation of the coal mine. Based on the very well-constrained location for recent seismic events and on the work performed on the slight amelioration over time of the RMS precision of the location between 1962 and 2009 in this Provence area [START_REF] Cara | SI-Hex: a new catalogue of instrumental seismicity for metropolitan France[END_REF], a location error of approximately 10 km is considered, especially for lower magnitudes, i.e., less than Mw2.5. Then, if we take into account a radius of almost 10 km around the Gardanne coal mine, and even less, there is a clear cluster of seismicity in this area, particularly during exploitation until February 2003. This seismicity seems to have continued in recent years but at a lower level. These observations obtained from French national networks (gathered within the research infrastructure of the French seismological and geodetic network RESIF-EPOS-www. resif. fr) provide a global view in recent decades. After the closure of the coal mine in February 2003, a local seismic network at the scale of the Gardanne mining basin with the specific objectives detailed above has been operated by Ineris since the end of 2007. Following the occurrence of events strongly felt by the population, seismic monitoring has been reinforced by the BRGM (since mid-2013). Local networks reveal local-scale seismic activity in the former coal mine area.

As stated above, this paper focuses on the recent seismic activity obtained by the deployment of a temporary seismic network during the 2013-2017 period in a specific area of the former Gardanne coal mine (FG area).

Fuveau-Greasque Seismic Activity

The inhabitants of the Le Jas de Bassas and Galoubet districts (commune of Fuveau) first felt "significant vibrations" during the month of June 2010. These observations coincided with the seismic activity recorded by the permanent Ineris network [START_REF] Dünner | annual monitoring report for areas overlying mining buildings at risk of sinking of a brittle nature -coal basin of Provence[END_REF]. No recorded seismic events or evidence of vibrations was known in this sector before 2010, following the flooding of the mine workings.

On November 1, 2012, at 3:58 am at local time (2h58 UT), a seismic event followed by a seismic swarm was widely felt by the population of the communes of Fuveau (west) and Greasque (north) for about three days. The lack of automatic reports by the French Alternative Atomic Energies and Atomic Energy Commission (CEA) and SISMALP Observatory and the local testimonies of populations indicate a low energy event (magnitude probably less than 2.5) and a superficial origin (less than 5 km). This seismic event was manually determined by different networks/organizations, such as Ineris, CEA, OCA (Côte d'Azur Observatory, Nice) and RéNaSS (National Network of Seismic Surveillance, Strasbourg) (Table 1). As mentioned above, epicentral determination uncertainties for national seismic observatories (CEA, RéNaSS) can reach 10 and even 30 km.

Methods

Local Temporary Seismic Network of the Fuveau-Greasque Area

Following this seismic swarm, the State decided to strengthen the microseismic monitoring system in the eastern Provence coal basin by supplementing the system set up in 2013. A temporary network of four stations, growing to five in 2015 and then nine in 2018, was installed by BRGM in the area (Fig. 4). Three-component accelerometers are located on the ground in private homes. Although the background noise level is not optimal, this network was deployed quickly at a low cost and has allowed us to follow the continued seismic activity, including the seismic swarms that occurred in 2014, the end of 2016, and in 2017. To achieve good sampling for signal processing, the acquisition was set to 1000 Hz. The triaxial broadband accelerometric sensors Implantation des réseaux permanent et temporaires thus allow the recording of regional seismicity as well as small local events without risk of signal saturation in the case of strong motion. This complementary seismic network acquires data continuously or at a triggering threshold. Figure 4a shows the deepening of mine workings from east to west. In the FG area (Fig. 4b), the deepening angle is very low, with a dip of 3° to 7° towards the northwest.

Data Processing

We use an adapted version of the HYPO71 location program [START_REF] Lee | Hypo71 (revised): a computer program for determining hypocenter, magnitude, and first motion of local earthquakes[END_REF] to locate seismicity. The velocity model is based on borehole log profiles and on knowledge of the geology. A velocity model over the Gardanne area provided values for compressional P-wave velocities (Vp) and shear S-wave velocity (Vs) for different geological layers [START_REF] Bertil | Study of the seismicity of the Gardanne basin, French Alternative Energies and Atomic[END_REF]. We have refined this model, including the ratio of compressional wave velocity to shear wave velocity, Vp/ Vs (Table 2) by minimizing the residues on the arrival times of the waves and the root mean square error (RMS) obtained at the locations. Some tests on quarry shots allowed us to check the velocity model and Vp/Vs ratio, and to validate the location results. The significant variations in water saturation of the locally saturated environment at ≈ 0 m ASL depth and oscillations of several meters to tens of meters could also affect this ratio. The RMS on all residuals of the travel times to the stations obtained by location by HYPO71 is less than 0.01 s for all events, which indicates a correct quality of the locations. We determined that the location uncertainty of local seismic events is almost 50-100 m inside the temporary network and 100-200 m outside.

Another point is the local response of the sites, which is a key element of any seismological and seismic hazard analysis. The use of background noise recordings generally allows the experimental determination of lithological site effects that may significantly modify the soil response [START_REF] Nakamura | A method for dynamic characteristics estimation of subsurface using microtremor on the groundsurface[END_REF]. This is also ideally suited to regions with moderate seismic activity, where there are few records of ground motion, and to urban or industrial contexts, where the noise level is potentially high [START_REF] Bour | On the use of microtremor recordings in seismic microzonation[END_REF]. The characterization of the soil responses permits us to evaluate any overestimation of the signal that could have an impact on the calculation of the magnitude.

We checked the noise environment of seismic stations that record the events. For all seismic stations, the soil responses were flat, and the horizontal to vertical spectral ratio H/V was less than 2, i.e., there was no observation of apparent site effect. In contrast, the systematic calculation of local magnitude based on amplitude and distance shows that the farthest northeastern station, "1466", could underestimate the amplitude of the seismic signal regardless of the location of the event inside the network. This could be due to the poor condition of the sensor installation, which is not firmly anchored to a hard floor.

Finally, a seismic catalogue based on manual pickings of P and S arrivals of 2688 events was established between mid-2013 and December 2017. Then, for some large magnitude events of each sequence, we calculated the focal mechanism using full waveform inversion in a limited frequency band [START_REF] Aochi | Mechanism of the ML4.0 25th April 2016 earthquake in southwest of France in the vicinity of the Lacq gas field[END_REF]. The elastic response to a unit double-couple force (Green's functions) are calculated using the same 1D model in Table 2. We fixed the epicentre position and varied the focal depth to minimize the misfit, which was calculated for the L-2 norm (least squares). We chose to directly invert the focal mechanism of the doublecouple seismic moment (three parameters: strike, dip and rake) rather than the six independent moment tensors. Additionally, the moment magnitude (Mw) and origin time (t 0 ) are unknown. The parameters were searched using a genetic algorithm [START_REF] Goldberg | Genetic Algorithms in Search, Optimization, and Machine Learning[END_REF]) that has previously been adopted in similar inversions [START_REF] Ulrich | Rapidness and robustness of finite source inversion from elliptical patches method using continuous GPS and acceleration data: 2011 Mw9.0 Tohoku earthquake[END_REF][START_REF] Aochi | Mechanism of the ML4.0 25th April 2016 earthquake in southwest of France in the vicinity of the Lacq gas field[END_REF].

Results

Seismicity and Ground Motion Record

The BRGM temporary network successfully revealed the unexpected number of local seismic events in the FG area (Fig. 5, Supplemental Fig. S-1). We could detect three times more events than those detected using the Ineris permanent network. Correlatively, the threshold of magnitude detection was improved from 0 down to -0.5 or even less.

Nearly 2688 seismic events were detected and located between October 2013 and December 2017, with some swarms generating up to 250 events per day. Typically, these swarms for the most part were felt by the population and lasted a few days: November Figure S-2 in the supplemental information file shows a plot of the peak ground acceleration (PGA), defined as the maximum values of the vector sum of the three component acceleration time histories, with respect to hypocentral distance. Based on the relationship between PGA and intensity [START_REF] Worden | Probabilistic relationships between ground-motion parameters and modified Mercalli intensity in California[END_REF], almost 15% of the recorded PGA values were susceptible to be felt by the population, and the highest PGA values reached the threshold of very slight damage to buildings from mid-2013 to late 2017. We also observed that the hypocentral distances of the recorded events were mostly between 0.4 km and 2 km. As most focal depths were from 400 to 800 m, the epicentral distance was of the same order as, or even less than the depth.

Focal Mechanism and Implications

Table 3 and Fig. 6 summarize the focal mechanisms we obtained. We chose a different bandpass filter at ≈ 1 Hz to clearly extract the pulse and varied the focal depth every 100 m with respect to the location catalogue. We do not always use the closest station for each event to avoid any uncertainty of the seismic event location. Globally, we could fit the waveforms well. For seismic events of magnitudes less than 1.3, the signals were not strong in this frequency range, which made it difficult to extract the main phase of the signals from the background noise level. We adopted the focal depth that produced the minimum misfit, and the result was consistent with the catalogue. All the analysed seismic events were found to be deeper than the mining gallery, indicating reactivation of the natural fault system.

All obtained mechanisms show normal faulting, mostly NW-SE strikes. The strike trend inferred from the fault mechanism should be discussed statistically, as the quality of the result varies one after another due to the number and coverage of the stations, the S/N ratio and a complex source mechanism. Under these limitations, the solution may not always be unique. For the example presented in Fig. 6, the waveform fitting is sufficiently good for the largest event on January 31, 2017. Although the best solution indicates a fault strike in the NW-SE direction, the second preferable solution still has a strike of N210°E (NNE-SSW). We cannot simply average them, as each solution might have a local minimum. However, among the different global searches of the parameters indicating the same strike direction, the variation remained small for this event (for example, strike and dip less than a few degrees and rake less than 10°). On the other hand, the records of the second largest event on December 1, 2014, showed two bursts of energy in the seismograms, indicating two events. Thus, the obtained mechanism should be regarded as a mean feature of the total process. From Table 3, the NW-SE striking mechanism (four among the six events analysed here) had an average of N301°E (N59°W), indicating the representative tectonic system shown in Fig. 2. The tectonic interpretation also indicates different directions of distensive stress and open structures (partens) in the surroundings; therefore, the diversity of the focal mechanism may have mechanical causality (Fig. 2).

Seismicity and Hydrological Considerations

As described above, the mine water level often fluctuates. These variations are due to the pumping capacity mainly at (a) 1 3

the Gerard shaft and effective rainfall. This term "effective" is used because it alone contributes to the reconstitution of aquifers and depends regionally on the type of soil. The data are recorded in Aix-en-Provence.

In 2010 (February, May, and especially June), some seismic activity was recorded on the permanent Ineris seismic network in the FG area, although no events had been detected since its installation on November 1, 2007[START_REF] Dünner | annual monitoring report for areas overlying mining buildings at risk of sinking of a brittle nature -coal basin of Provence[END_REF]. Moreover, no event had been felt since the mine closure in 2003. As shown in Fig. 7, this time coincided with the reaching of the -30 m ASL of the mine water table at the pumping Gerard shaft (Fig. 4). We can also distinguish three periods of mine aquifer evolution, depending on pumping capacity and rainfall after 2010 (green curve in Fig. 7):

• In the first period, from mid-2010 to mid-2016, the seasonal fluctuation in the level of the mine aquifer was on the order of 20 m (approximately between -10 m ASL and + 10 m ASL). • The second period, from mid-2016 to the end of 2017, when four pumps were replaced at the same time as a very dry period, which produced a decrease in the water level of ≈ 20 m, reaching the water levels in early 2010 (from -10 m ASL to -30 m ASL) by the end of 2016. • Since the end of 2016, these levels have stabilized, with smaller fluctuations of about -30 m ASL.

We focus now on these three periods, with an emphasis on the locally instrumented period from the end of 2013 to the end of 2017.

First Period (mid-2010 to mid-2016)

In the first period, a previous study [START_REF] Dominique | Microseismicity in a mining context-The Gardanne coal basin (Bouches-du-Rhône)[END_REF]) distinguished two types of seismicity between mid-2013 and the end of 2015 triggered by oscillations of the water table :   • One, rather sparse type, on a kilometric or plurikilometric scale, with migration for several months, which followed the level of the mine water table. This microseismicity, often not felt by the population, was of low magnitude and can be considered microseismic background noise triggered by water table oscillations, particularly after rainy episodes. This seismic background noise extended over a few square kilometers and mostly below the former mine workings; and • The other, which was more concentrated, occurred in swarms, often with higher magnitudes at the start of the swarm. The rainfall peaks again strongly affected the triggering of the very local microseismic swarms in the FG area. These swarms are most likely linked to the activity of faults crossing mine voids. This triggering occurred ≈ 10 days after the effective rainfall peak (and especially sharp rise of the mine water level) during the swarms in November 2012 and December 2014. While the rainfall peaks did not trigger a microseismic swarm, these swarms were systematically preceded by intense rainy episodes sharply raising the level of the mine water table.

The December 2014 swarm was imaged as a clear cloud with a NW-SE trend and a dip towards the NE of 60° (Fig. 8a). This dipping direction does not correspond to the focal mechanism of the 2014/12/01 Mw1.76 event (Table 3 and Fig. 8a), which seems to have a complex source mechanism, as explained above.

Figure 8c also shows the lateral extension to the northwest and southeast of the seismic zone in 2015 compared with that of the seismic swarm in December 2014 (Fig. 8b). This lateral extension is highlighted by the very shallow (less than 300 m) seismic events (in orange in Fig. 8c) towards the south-eastern part during April 2015 when the mine water table was at its highest level (+ 11 m ASL). This observation seems to correlate perfectly with a rise in the water table to levels never seen before. These events were mainly above the mine workings, which are intersected by the seismogenic zone. At this location, the mine workings are approximately 300 m deep (+ 20 m ASL). This indicates that these seismic events were mostly above the mine workings and could have affected them. We also observed a very slight trend in the seismicity deepening towards the northwest in parallel with the mine workings.

The estimate of the magnitude considering activation of the entire active surface in 2015-in the event of rupture of the entire seismogenic area-can be made by making the following assumptions (Fig. 8c). Geometrically, we can consider a fault width of 500 m and a fault length of 1400 m for an equivalent source radius of 440 m. A stress drop value Δσ of 1 MPa can be considered acceptable in relation to the type of movement for these low magnitude and shallow depth events, which seem to be factors favourable to low values of stress drop [START_REF] Hough | Shaking from injection-induced earthquakes in the central and eastern United States[END_REF][START_REF] Hough | Shaking intensity from injection-induced versus tectonic earthquakes in the central-eastern United States[END_REF][START_REF] Pacor | Spectral models for ground motion prediction in the L'Aquila region (central Italy): evidence for stress-drop dependence on magnitude and depth[END_REF]. Using a simple model [START_REF] Brune | Tectonic stress and the spectra of seismic shear waves from earthquakes[END_REF], we obtained a moment magnitude (Mw) close to 3.5. Detailed analysis of sparse seismicity makes it possible to identify concentrations or alignments of events, highlighting new segments that could potentially play a role in the swarms. There are thus three parallel zones trending NW-SE and another smaller swarm trending WNW-ESE to the south (Fig. 9). The microseismicity seems to follow the orientation of known or supposed tectonic faults in the area. This agrees with the identified clusters obtained by automatic location by [START_REF] Namjesnik | Automatic detection and location of microseismic events from sparse network and its application to post-mining monitoring[END_REF].

Second and Third Periods (mid-2016 to the end of 2017)

The second and third periods from mid-2016 to end-August 2017 were marked by frequently repeated seismic swarms stronger in number and energy than those recorded since 2014 (Figs. 5 and7, supplemental S-1). Nine swarms appeared during these periods compared with two swarms during the first period. The effective rainfall peak and especially sharp acceleration of the rise of the mine water level can explain the first period triggering. This is not the case for these two periods. Therefore, another mechanism for triggering seismicity is proposed.

The beginning of the second period was defined by the dewatering of the flooded mine workings after six years (Fig. 7). Thus, this suggests that over six years, hydrostatic (and seismic) stability was linked to the loading of a water column oscillating between 20 and 40 m (-10 m to + 10 m ASL) above the -30 m ASL. Since mid-2016, the dewatering and depressurization caused by a decrease of 20 m of this mine water level (at -30 m ASL) broke a certain hydrostatic equilibrium (green curve in Fig. 10b). It should be noted that this level of the mine aquifer of -30 m ASL corresponds to the first local seismicity in the FG area in 2011.

Since December 2016, the mine water level has been relatively stable at -30 m ASL. Seismic activity increased progressively until March 2017 and then decreased slightly. The last peak of activity began in August 2017, showing a response time of up to ≈ 8 months in reaction to the modification of the state of stress. Since the end of 2017, the mine water table has not experienced significant fluctuations, and seismic activity has remained moderate, with about a dozen events per month.

Seismicity Migration with Mine Water Level

We attempted to statistically correlate the spatial migration of seismicity to the water table evolution of the mine aquifer. For this purpose, the FG area is divided into ten bands 280 m wide (200 m latitude and 200 m longitude) oriented in a NE-SW direction (Fig. 10a).

For each band, the monthly number of seismic events was counted from October 2013 to December 2017. The horizontal histograms in Fig. 10b represent proportion of monthly events in each of the ten bands. For each month, the colour of the band (or geographic portion) varies according to the number of seismic events. The red bar corresponds to the band having the greatest number of seismic events. The black bars correspond to the bands having a median number of seismic events. The grey bars correspond to the bands containing the fewest seismic events. The months whose histograms are completely grey correspond to a monthly seismicity less than or equal to 10 events in the entire FG area (extent of the 10 bands). We consider this low monthly seismicity to be "seismic background noise", and not representative of the analysis. In other words, the brown curve giving the number of seismic events per month in the FG area indicates the representativeness of the data provided in the histograms. If the number of events is less than or equal to 10, the histogram is not very representative. If the number of events is large (e.g. > 50 events during a seismic swarm), the histogram can be considered very representative for seismicity analysis.

We found a good correlation between the level of the mine aquifer measured at the Gerard shaft (green curve) and the maximum monthly number of seismic events by band area (red bar of the histograms). Fig. 10b clearly shows that over the 2013 to 2017 period, the seismicity was spatially well correlated with the mine aquifer.The low dip of the mine workings (3 to 7° to the northwest) facilitates the lateral drainage of water and the horizontal migration of seismicity. We observe the migration of seismicity to the southeast when the level rises and to the northwest when the level falls. It is interesting to note that this migration is almost synchronous -at the monthly time scale -with the observed water levels. There is probably a slight delay of a few days because the water level measurement is at the Gerard shaft, 6 km west of the FG area (Fig. 4). However, the karstic environment in the FG area greatly facilitates the circulation of water and could explain the small time lag between variations in the water level and the generated seismicity. This rapid circulation was also confirmed by the delay -relatively short -of ≈ 10 days between effective rains arriving in the basin and the development of seismicity in the FG area [START_REF] Dominique | Microseismicity in a mining context-The Gardanne coal basin (Bouches-du-Rhône)[END_REF], as also observed in other contexts [START_REF] Bollinger | Hydrological triggering of the seismicity around a salt diapir in Castellane. France[END_REF][START_REF] Dominique | Instrumentation of a small seismic crisis in Salsigne (Aude)[END_REF][START_REF] Hainzl | Evidence for rainfall-triggered earthquake activity[END_REF][START_REF] Kraft | Meteorological triggering of earthquake swarms at Mt Hochstaufen[END_REF][START_REF] Rigo | Seismicity rate and wave-velocity variations as consequences of rainfall: the case of the catastrophic storm of September 2002 in the Nîmes Fault region (Gard, France)[END_REF].

One hypothesis related to this migration could be the arrival of water, which would locally modify the velocity The horizontal monthly histograms represent the proportion of seismic events normalized in each band. The red bar shows the band containing the greatest number of monthly seismic events. The black bar shows the band containing a median number of monthly seismic events. In grey, the histograms for which the monthly number of seismic events is the lowest number of monthly seismic events and considered not to be representative (i.e., less than or equal to 10 seismic events in the FG area). The bold number on the right side of the figure is the band number 1 3

model sufficiently to influence the location of the seismic events and create apparent migration. We did not truly consider this because the "water" interface in the mine voids (= galleries) is too thin (1 to 2 m) to play a role in seismic wave propagation. The waves propagate mostly vertically. For comparison, the wavelength of the P waves remains at a few hundred meters. Thus, such a fluctuation in the structure is always possible, but this does not influence the location analysis, especially if one considers other sources of error (e.g. picking of arrival times and velocity model). For focal mechanism analysis, we use lower frequencies and uncertainty in location is less important. We can estimate that for a water height variation of 1 m in the mine aquifer, the seismicity migrates on the order of 20 to 40 m towards the SE (if the water level increases) or towards the NW (if the water level decreases). Of course, this estimate is only a trend and does not prevent the development of sparser seismicity.

We also observe that all seismic swarms after the end of 2016 (water level at ≈ -30 m ASL) are found further northwest than those in December 2014 (water level at ≈ -5 m ASL), that is, in bands with a number greater than 5. Likewise, the last major swarm in August 2017 was the most northwestern swarm (band 8) (Figs. 9 and10b).

The Case of Pillar or Roof Fall/Rock Burst

From the beginning of significant seismicity in the FG area with the seismic swarm in November 2012, collapses of the roof of the galleries or pillar bursts were suspected. This is one of the main reasons why the local temporary seismic network was deployed in the FG area. However, the frequency, nature, epicentre location, and focal depth of the seismic events were not compatible with such a phenomenon. Finally, out of thousands of seismic events analysed, only the event on February 1, 2017, 4h45 UT, shows the characteristics of an underground mining seismic event (blast or roof collapse; asterisk in Table 4).

Several arguments support this hypothesis of a mining event. First, the shape of the signal can be easily differentiated from the events that occurred just before. Two events of equivalent magnitude at 1.5 and 1.2 s earlier (Table 4 and Supplemental Fig. S-3a) are compared. The gradual monochromatic rise "in the funnel" of the wave differs from that observed just before, representative of an earthquake. The Fourier spectrum shows large differences between both types of waves. The spectrum of the mining event shows a higher frequency content (between 65 and 80 Hz) than that of natural seismic events (between 15 and 30 Hz) .

Then, this event occurred after a series of three seismically felt events: on January 18, January 31, and February 1. The swarm of January 31-February 2 generated 38 events on January 31. The strongest event on January 31 at 2h54 UT had the greatest acceleration (nearly 70 cm/s 2 ) recorded on a sensor (ROSS station) since installation of the local network, producing the maximum moment magnitude of Mw1.94 of the catalogue (Table 3). It also occurred just 1 s after two nearby significant events of ≈ magnitude 1. We can therefore assume that a series of these significant events may have destabilized the roof of the galleries or certain pillars, causing their collapse. We also note that the depth of this event was more superficial than most of the seismicity in the sector. Finally, it was isolated without any direct association with swarms already identified in this area (e.g. Figure 9 and supplemental Fig. S-4). All these elements suggest that a local destabilization of mine workings occurred on Feb. 1, 2017.

Conclusion

This paper explores the seismicity triggered by the rise and oscillations of the water table in former mine workings in Gardanne, southern France, after mine closure in 2003. We focused on the period between 2013 and 2017 using a new seismic network installed after the seismic swarm of December 2012.

Our observations showed that triggered seismicity closely depends on the period and amplitude of mine water table variations. The low dip of the mine workings used as a privileged drain favoured a seismic migration in a NW-SE direction according to the water table oscillations. Most seismicity was found below the mine workings between 400 and 800 m below the ground surface (100 to 500 m below the mine workings). In addition, it was observed that the shallowest events (≈ 200 m) occurred above the mine workings and in the southeastern part of the FG area (Fig. 8). This also seems consistent with a rise in mining levels towards the southeast.

The focal mechanisms achieved in this study, the northeastern dip of the main faults, the NW-SE direction of the events of the 2014 swarm, and the extension regime described by several authors in the FG area [START_REF] Chalumeau | Hydrogeological functioning of lignite operations in the Fuveau basin (Bouches-du-Rhône)[END_REF][START_REF] Gaviglio | Consequences of structural change on the physical characteristics of fuvelian limestones in the Arc Basin (Bouchesdu-Rhône, France): fracturing, mechanical properties, stress states[END_REF][START_REF] Gaviglio | Measurements of natural stresses in a Provence mine (Southern France)[END_REF]) agree with normal faults in a NW-SE direction (Figs. 2,6,and 9). The two main clusters in the NW-SE direction are on the eastern part of the exploited mine workings. The upward extension of the fault surface crossing the mine workings at a depth of ≈ 300 m (≈ 0 m ASL) arrives in the middle of the exploited panels, which can be considered as reservoirs. This configuration would correspond to the destabilization scheme proposed by [START_REF] Talwani | On the nature of reservoir-induced seismicity[END_REF], taking up an analysis by [START_REF] Roeloffs | Fault stability changes induced beneath a reservoir with cyclic variations in water level[END_REF].

As described by other authors for reservoir-induced seismicity (RIS), this seismicity can be qualified as continuously triggered seismicity or protracted seismicity [START_REF] Talwani | On the nature of reservoir-induced seismicity[END_REF]) due to its persistence over time. In addition to the highly activated seismic periods, nearly ten events are still recorded every month. This implies that this seismicity will exist as long as the water level remains and fluctuates above -30 m ASL. A clear migration of seismicity with the oscillations of the mine aquifer depending on the pumping capacities and effective rains is highlighted. The seismically activated area corresponds to a fault size of Mw3.5 earthquake. All these points lead to the necessity of continuous hydraulic and seismic monitoring to better understand these phenomena and to assess the related risks.

Fig. 1

 1 Fig.1Geological units in the Pyrenean-Provençal domain (afterRouire 1979, modified)

Fig. 3

 3 Fig. 3 Seismicity around the Gardanne Mine. a Historical seismicity from the macroseismic database (SisFrance 2016). The black rectangle is the zoomed area shown on the right. b Instrumental seismic-

Fig. 4

 4 Fig. 4 Evolution of the temporary seismic BRGM network. a Five temporary BRGM stations (purple triangles) are superimposed on a map of the Gardanne coal basin. Ineris permanent network (green triangles). Main faults (bold black lines): DF Diote Fault, MF Meyreuil Fault, and JLF Jean-Louis Fault). Location of the seismic swarm in 2012 (red star). The pumping Gerard shaft (orange dot) and type of mine workings: room and pillars (vertical lines), long walls (dots),

Fig. 5

 5 Fig. 5 Histogram of the FG area seismicity showing the number of seismic events per month from October 2013 to December 2017. We can distinguish the seismic swarms of December 2014 and August

LPFig. 7

 7 Fig. 7 Rise of the mine water table recorded at the Gerard shaft (orange dot in Fig. 4) after the 2003 mine closure of Gardanne with the three periods of evolution after 2010. Mine aquifer level at the Gerard shaft (green) and cumulative number of seismic events in the FG area (blue) recorded by the Gardanne permanent seismic network

Fig. 8 a

 8 Fig. 8 a 3D distribution of seismic events in the December 2014 swarm. View from the southeast looking to the northwest perpendicular to the fault plane. The size is proportional to the magnitude (from M -1.3 to M 1.4). The colour corresponds to the depth. In black is the mine workings. In yellow is the surface digital elevation model (DEM) in metres (0 corresponds to + 330 m ASL). The four seismic stations that recorded the seismic swarms are positioned on the surface DEM. b 3D distribution of seismic events in the December

Fig. 9

 9 Fig. 9 Seismicity 2013-2017 in the FG area. In the background, the chamber-and-pillar type mine workings. The lines represent ASL. The seismicity is developed mainly along NE-SW structures parallel to the main trend of the mine workings. Large red stars symbol-

Fig. 10 a

 10 Fig. 10 a Ten geographical bands with their number where the monthly seismicity is counted from October 2013 to December 2017 in the FG area. In the background, the chamber-and-pillar type mine workings and seismicity from October 2013 to December 2017. The lines represent ASL. b Mine aquifer level at the Gerard shaft (green curve) and monthly seismicity in the entire FG area (brown curve). The horizontal monthly histograms represent the proportion of seismic events normalized in each band. The red bar shows the band

Table 1

 1 Characteristics of the seismic event on November 1, 2012. Md (duration magnitude), MLv and Ml (local magnitude), and Mw (moment magnitude)

	Organization	Origin Time (TU)	Latitude	Longitude	Magnitude	Depth (km)	Suspected type
	Ineris	02h58 31	43°26.520'N	5°31.678'E	Md 2.9	0.55	Mine collapse
	OCA	02h58 30	43°29.100'N	5°36.900'E	MLv 2.57 ± 0.38	-	-
			± 42.5 km	± 19.4 km	Mw 3.15
	CEA	02h58 33	43°24.852'N	5°32.826'E	Ml 2.2 ± 0.38	2	Induced event
			± 3.1 km	± 1.4 km	Md 2.2
	RéNaSS	02h58 28	43°09.900'N	5°28.200'E	MLv 2.2	0 (fixed)	Quarry blast
			± 10 km	± 5 km

  At the end of 2021, only one swarm has occurred on June 12 and 13, 2021 since 2017. The main swarms occurred during the winter season, except for that

	1 to 3, 2012; December 1 to
	4, 2014; December 8, 2016; January 21, 2017; January 31
	to February 2, 2017; February 6 to 7, 2017; February 17 to
	18, 2017; March 9 to 10, 2017; August 22 to 24, 2017; and
	September 8, 2017.

Table 2

 2

	Seismic velocity model adopted for this study with Vp/ Vs = 1.9	Depth (km)	Vp (km/s)
		0-1	3.85
		1-4	5.0
		4-15	5.7
		15-25	6.7
		> 25	8.0

Table 3

 3 Parameters of the obtained focal mechanism analysis for large magnitude events. Mw: moment magnitude; M BRGM : local magnitude; and M Ineris : duration magnitude. Although this paper focuses

	(b)	EW			NS			UD
	A	A=7.2e-04			A=5.5e-04			A=9.0e-04
	0									1466
	-A	A=3.6e-04			A=6.6e-04			A=1.6e-03
										1418
		A=7.0e-04			A=7.8e-04			A=1.4e-03
										ROSS
	Ground velocity (m/s)	A=2.4e-03 A=7.5e-04 A=6.0e-04 A=5.6e-04			A=2.3e-03 A=5.1e-04 A=4.3e-04 A=6.4e-04			A=3.2e-03 A=1.2e-03 A=1.1e-03 A=1.1e-03	BULL SAVA VILO RAMP
		A=2.3e-04			A=9.2e-04			A=1.0e-03
										VERW
		A=6.0e-04			A=7.6e-04			A=1.4e-03
										BARL
	0		4	8	0	4	8	0	4	8
			Time (s)			Time (s)			Time (s)

Fig. 6 a Focal mechanisms for the seismic events listed in Table 3. The solid and open triangles represent the BRGM stations and the other network, respectively. BRGM stations are only used for the focal mechanism analyses. The framework by the broken lines indicates the FG area shown in Fig. 4b. b The fitting of seismograms for the 2019/01/16 Mw1.34 event. The velocity waveforms are filtered between 1-2 Hz for both the synthetics (black lines) and the observations (grey lines). The amplitude is normalized for each panel with a constant A

21/10/16 8/12/16 31/1/17 6/2/17 21/1/17 17/02/17 22/8/17 1/12/14 9/3/17 08/09/17 1 km Rock burst Depth

Table 4

 4 Parameter characteristics of the main seismic events in January and early February 2017. The event marked by an asterisk is suggested as a mining event. The seismograms of 01/02 16h45 are shown in Fig.13. The map of seismic events is shown in Fig.14

	Date	18/01	31/01	01/02 16h45			01/02	02/02
	Hour, Minute	23h02	02h54				16h48	05h01
	Second (TU)	5.3 s	28.1 s	33.9 s	34.2 s	35.4 s (*) 49.3 s	23.9 s
	Latitude	43.441°N	43.443°N	43.444°N	43.444°N	43.442°N	43.444°N	43.444°N
	(° decimal)	± 50 m	± 50 m	± 50 m	± 50 m	± 50 m	± 75 m	± 50 m
	Longitude	5.539°E	5.534°E	5.534°E	5.534°E	5.533°E	5.533°E	5.525°E
	(° decimal)	± 50 m	± 50 m	± 50 m	± 50 m	± 50 m	± 75 m	± 50 m
	Depth (m)	500	640	640	640	240	860	610
		± 50	± 50	± 50	± 50	± 150	± 75	± 50
	Local Magnitude 1.4	1.6	1.1	1.1	0.8	1.6	1.2

Supplementary Information

The online version contains supplementary material available at https:// doi. org/ 10. 1007/ s10230-022-00860-z. The BRGM seismic catalogue is freely downloadable at through the BRGM https:// doi. org/ 10. 18144/ 0__ f4f64 ebf-b3b5-4462-b374-61955 a40a0 fe. We used MATLAB ver R2018a for some figures and the Software SeisGram2K Seismogram Viewer v7.0.0X10 (29Mar2016) (BETA)-Copyright (C) 2000-2016 Anthony Lomax (www. alomax. net) for data viewing and processing. Finally, we are very grateful to the three anonymous reviewers whose comments greatly improved the quality of the manuscript.

Acknowledgements

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.