domain: info spi

text:
Audrey Ahoukeng Donwoung
email: dreyswoung@gmail.com

Bachir Nazim Agoulmine
email: nazim.agoulmine@univ-evry.fr

Bachir Djafri
email: bachir.djafri@univ-evry.fr

Tooled

Tooled Methods for the Simulation of Deployment of Safe Services in the Internet of Vehicles

Keywords: Internet of Vehicles, Edge Computing, Cloud Computing, Transportation Services, Formal Modelling, SysML, SUMO Simulator, Distributed Systems

Introduction

Nowadays vehicles have become essential items in daily life of people and in the management of big cities. The growing number of these vehicles on the roads increases the need to significantly reduce the traffics, pollution emission, risk of personal injury, material damage and so on. According to Transport and Environment organisation, 1 , the almost a third of trucks tested emit dangerous levels of air pollution. Similarly, road safety situation report carried out by the World Health Organization (WHO), road accidents are responsible for 1.35 million2 deaths every year worldwide. Many researchers have looked at the Internet of Vehicles (IoV) as one potential solution to help improving the situation. This concept originated to connect multiple users with vehicles to ensure safe and efficient road conditions through one communication from vehicle to any entity (V2X) which must be reliable and in real time.

Since vehicles are systems equipped with computing and communication capabilities and more and more edge computing infrastructure is being deployed (Edge Computing), vehicles can communicate directly together if close enough (V2V) and/or with the infrastructure (V2I) at any time. With their more and more increased computing capacity, they are also able to collect and process a large amount of information in real time (weather, road conditions, traffic conditions, etc.) coming from their embedded sensors or from the infrastructure [START_REF] Bertolino | A tour of secure software engineering solutions for connected vehicles[END_REF]. However, since vehicles are mobile and move at variable speeds, the conditions in which they evolve are not fully predictable including their connectivity [START_REF] Hajlaoui | Résolution à base d'heuristiques du problème de routage dans les réseaux ad hoc de vehicules[END_REF]. Hence some services can not be executed only onboard and therefore the vehicles need to interact with far-end applications, such as intelligent driving assistance, and high resolution maps creation services. Cloud Computing and Edge Computing are those technologies that will help to execute these applications and reduce the response times distributing the processing between the vehicles, the edges and the cloud computing. Thus, this system formed by connected vehicles, the edges and the cloud computing can be seen as a distributed system of fixed and mobile components, connected through an underlying communication infrastructure such as 5G [START_REF] Rachedi | End-to-end network slicing for its-g5 vehicular communications[END_REF].

In this paper, the objective is to model, analyze and validate a representative since that is a route optimization service for the collection, transportation and deposit of waste in a city. The aim is to formally design the waste management planning service and verify its operation and behaviour by simulation to increase confidence in the services before its real deployment. The used simulator is called SUMO (Simulation of Urban MObility) [START_REF] Lopez | Microscopic traffic simulation using sumo[END_REF]. It is an open-source and highly portable simulator tool that allows modelling of intermodal traffic systems including road vehicles, public transport and pedestrians. This paper is structured in 5 sections. Section 1 is a general introduction to the context. Section 2 presents the objectives of the paper. Section 3 presents the problem statement where we describe the system model and system properties formal modelling. Section 4 presents an implementation and an evaluation in SUMO. Section 5 is a performance evaluation after doing simulation tests in SUMO.

Related Works

Internet of Vehicles has been the topic of numerous contributions in the literature. Several surveys have been published [START_REF] Xiaolong | Artificial intelligence for edge service optimization in internet of vehicles: A survey[END_REF] highlighting the challenges in the area in terms of infrastructures and mechanisms to support efficient services deployments and connectivity among vehicles and infrastructures and between vehicles [START_REF] Jerbi | Protocoles pour les communications dans les réseaux de véhicules en environnement urbain : routage et geocast basés sur les intersections[END_REF]. More recently, more focused papers have been published addressing the interactions between the vehicles, the Edge and the Cloud Computing infrastructures (that corresponds to our considered context). For instance, in [START_REF] Liu | Edge computing for autonomous driving: Opportunities and challenges[END_REF] authors present a framework to build advanced computing systems for autonomous vehicles. The paper addresses three categories of problems: cutting-edge computer system designs (real-time data processing as autonomous vehicles are mobile and often have power restrictions in term of energy), V2X applications (how vehicles cooperate with each other and with the infrastructure) and autonomous vehicle security (protecting autonomous driving on-board computer systems from attack). In [START_REF] Kwon | Web-based car control and monitoring for safe driving of autonomous vehicles[END_REF] authors focus on the design of a web-based car monitoring solution that allows network vehicles to exchange their mobility information to get quick road awareness using the specification of the service vehicle information and vehicle signal specification. In [START_REF] Chen | Optimization of classified municipal waste collection based on the internet of connected vehicles[END_REF], authors propose a solution to optimize waste collection and optimize vehicles service routes in real time based on signals sent by intelligent waste bins. In [START_REF] Nowakowski | Vehicle route planning in e-waste mobile collection on demand supported by artificial intelligence algorithms[END_REF], authors investigate using artificial intelligence algorithms for solving the vehicle routing problem with time windows for a heterogeneous fleet of waste collection vehicles. Authors of [START_REF] Edigleison Da | A platform for cloudification of network and applications in the internet of vehicles[END_REF] provide a contribution that aim to evaluate the performance and mobility of vehicular networks using the SUMO simulator used in this work. In [START_REF] Atiquzzaman | Dfcv: A novel approach for message dissemination in connected vehicles using dynamic fog[END_REF], authors aim to reduce the delay and guarantee the delivery of messages in city environment where the density of vehicles is high. The authors introduce the concept of Dynamic Fog (DFCV) that consists on using the Fog Computing and Cloud Computing for the dissemination of messages. All these works have been inspiring this contribution. However our focus in this contribution is different as it aims to address the specification of a vehicular service evaluate its behavior in term of adaptation in a dynamic context that is the collection of waste in a city.

The research questions we aim to address in this paper are the following: How to model a vehicular service using a formal modelling language? How to model constraints in terms of mobility starting points and ending points (i.e. trucks parks)? And how to model constraints in terms of collecting point of waste containers?

The used method consists on specifying the different representative scenarios of the waste containers collection and identify the functional and non-functional goals (i.e. requirements) of such a service in the frame of these scenario. In addition, the method aims also to identify the constraints on the connected vehicles in the context of each scenario i.e. the context in which the waste collection planning service will be deployed (i.e. of paths of connected trucks in a city to collect waste containers) and verify whether the designed service is compliant with the identified requirements at run-time.

4 Proposed Solution

System Model

The used high level defined scenario is the following : A company is responsible of waste collection in a city. It has deployed a set of waste collection points in the city where a set of connected containers are installed. Citizens use these containers to deposit their waste over time. Containers could only be unloaded at specific time intervals by the company defined by the city regulation. For that, the company has a set of electrical and connected trucks that are used for collecting the waste. Trucks are parked in different sites in the city. Containers are also installed in different sites in the city. When they are filled above a certain predefined level (threshold), a message is sent by the containers embedded sensors to the back-end waste collection management service deployed in the cloud (via appropriate wireless network and protocol). At periodic time, this cloud service processes the messages and identifies the set of containers that need to be unloaded and their respective GPS positions. Next, this cloud service selects the set of trucks to use for the collection of waste and plan accordingly all their trips.

The proposed system model introduces three layers namely the connected vehicles layer (i.e trucks), the Edge layer and the Cloud layer:

• The three layers are connected together using network infrastructures. The Cloud Computing and the Edge layers are connected via fixed backbone network.

• The Edge layer consists on a set of road side computing and networking infrastructure to which vehicles (in the vehicular layer) can connect using wireless communications technologies (e.g. 802.11p).

• The terminal layer is composed of all the connected vehicles, connected waste containers, traffic lights, etc. These objects are responsible for the collection of data environment and transmit it to the edge layer so that data can be processed and stored.

The Cloud layer contains all services that are accessible by the connected vehicles. It is responsible for processing and storing the information received from each Edge node and sometimes returns a response following a request.

System Properties Model

The services offered by the distributed system (based on mobile and fixed components) in connection possibly with the safety of vehicles and the safety of passengers and pedestrians, may have requirements and constraints in terms of reliability and response time (so-called functional requirements) and so-called non functional requirements that result from the cooperation between distributed components of the system, fixed (Edge and Cloud) and mobile (vehicles).

Non functional requirements are not directly related to the services offered by the system, but can be nevertheless critical for the proper functioning of the system. These requirements may be related with time and/or space: V2I and V2V communication constraints (e.g. delay, speed), availability of processing resources (e.g. CPU, memory), in relation to arrivals and departures of vehicles (e.g. mobility model, breakdowns). These requirements may also be neither temporal nor spatial, but simply resulting from a failure or unavailability of resources in the infrastructure. Therefore, these non-functional requirements may no longer be satisfied (temporarily, partially or totally) over time.

In order to model the system properties, we will base on specific use-cases to identify them. These properties include the functional and the non functional requirements. Once these requirements are identified, we will derive sub-goals from these goals which will identify sub properties of the system. The process will stop when we'll know that these properties are either intrinsic properties of the vehicles, network, edge computing or cloud computing infrastructure or these are properties that should be provided by the cloud services or the connected vehicles controllers.

In the following, we present the different considered scenario to identify these system properties and model them formally using SysML modeling language [START_REF]Information technology -Object management group systems modeling language (OMG SysML). Standard, International Organization for Standardization[END_REF]. SysML supports all the cycle of software engineering from the specification, analysis, design, to the verification and validation of systems and systems-of-systems design. This language is chosen in this work for its high capacity to model any engineered system and also because it is a full ISO standard since 03/2017 [START_REF]Information technology -Object management group systems modeling language (OMG SysML). Standard, International Organization for Standardization[END_REF].

The requirements of the system have been captured in Papyrus SysML and specified as shown in figure 1. We mention in the diagram to verify that certain requirements are indeed satisfied by blocks from the block diagram. To also make a verification using the activity diagram.

The block diagram describes the architecture of a system as a block which can be decomposed recursively into other blocks to the basic elements of the structure. The structure of system which resulting is therefore a hierarchy of blocks. The blocks are linked together by relation which can mainly be either a composition or an aggregation relation. Our system is essential made up of vehicles, wastes containers, edges nodes and Cloud, hence the use of composition link to connect our different blocks. The wastes containers are equipped with an electronic card and sensor for their operation. However we will not enter into the composition of the vehicle.

In each block we define the different variables, operations to be performed and if necessary the constraints.

Preliminary Implementation and Evaluation

We have implemented these scenarios and corresponding use cases in the SUMO simulator (Simulation of Urban MObility). 1 In this tool, each vehicle is individually simulated and has its own route. We have then used SUMO version 1.14.1 for Windows and OpenStreetMap2 to create a map of the City of Evry (in France), convert it in a format understood by the simulator using netconvert tool and finally imported it in SUMO. Next, We have added additional objects (polygons) in the map representing mobility constraints (e.g building, rivers, etc.)3 .

Conclusion and future work

In this paper, we presented a work which objective was to formally model a innovative service in the context of connected vehicles. The chosen services is the waste collecting management service in a city and the chosen modeling language is SysML.

We have conducted all the phase of the proposed tooled method from the formal specification of the service towards it simulation in a realistic urban mobility simulator that is SUMO. The learned lesson is that this global process is very important for the specification, development and test of innovative services in the context of Internet of Vehicles to understand their functioning and make them more optimized and safe.

Our future work is to advance the modelling part by introducing more artefact in the formal models related to this ecosystem and also enhance the capaV2V communications that have not been considered in this work. These V2V communication can be of paramount benefit to this future ecosystem to leverage the provided services.

Figure 1 :

 1 Figure 1: SysML Requirements Diagram Specification

https://www.transportenvironment.org

1.35 million: number from https://apps.who.int/iris/bitstream/handle/10665/277372/WHO-NMH-NVI-18.20-fre.pdf

https://www.eclipse.org/sumo/

https://www.openstreetmap.org

https://sumo.dlr.de/docs/Networks/Import/OpenStreetMap.html


𝑎 A Utrobicic

Y Angelis

𝑏 J Bortfeldt

𝑐 F Brunbauer

𝑑 E Chatzianagnostou

𝑏 K Dehmelt

𝑒 G Fanourakis

𝑓 K J Floethner 𝑑

𝑔 M Gallinaro

F Garcia

𝑖 P Garg

𝑒 I Giomataris

𝑗 K Gnanvo

𝑘 T Gustavsson

𝑗 𝑙 F J Iguaz

𝑎 D Janssens

𝑛 A 𝑚

𝑗 M Kallitsopoulou

𝑜 P Kovacic

𝑗 M Legou

𝑝 J Lisowska 𝑑

𝑞 M Liu

𝑔 𝑟 S Lupberger

𝑘 I Malace

𝑑,𝑏 Maniatis

Y Meng

𝑞 H Muller 𝑑

𝑟 E Oliveri

𝑑 G Orlandini 𝑑

𝑠 T Papaevangelou

𝑗 M Pomorski

𝑡 L Ropelewski

𝑏 𝑑 D Sampsonidis

𝑑 𝑢 L Scharenberg

𝑟 T Schneider

𝑗 𝑑 L Sohl

M Van Stenis

𝑑 Y Tsipolitis

𝑣 S E Tzamarias 𝑏

𝑢 R Veenhof 𝑑

𝑤 X Wang

𝑞 S White 𝑑

𝑥 Z Zhang

𝑞 Y Zhou

A large area 100 channel Picosec Micromegas detector with sub 20 ps time resolution

𝑠 Friedrich-Alexander-Universität Erlangen-Nürnberg, Schloßplatz 4, 91054 Erlangen, Germany 𝑡 CEA-LIST, Diamond Sensors Laboratory, CEA Saclay, F-91191 Gif-sur-Yvette, France 𝑢 Center for Interdisciplinary Research and Innovation (CIRI-AUTH), Thessaloniki 57001, Greece 𝑣 National Technical University of Athens, Athens, Greece 𝑤 Bursa Uludaǧ University, Görükle Kampusu, 16059 Niufer/Bursa, Turkey 𝑥 University of Virginia, USA E-mail: antonija.utrobicic@irb.hr A : The PICOSEC Micromegas precise timing detector is based on a Cherenkov radiator coupled to a semi-transparent photocathode and a Micromegas amplification structure. The first proof of concept single-channel small area prototype was able to achieve time resolution below 25 ps. One of the crucial aspects in the development of the precise timing gaseous detectors applicable in high-energy physics experiments is a modular design that enables large area coverage. The first 19-channel multi-pad prototype with an active area of approximately 10 cm 2 suffered from degraded timing resolution due to the non-uniformity of the preamplification gap. A new 100 cm 2 detector module with 100 channels based on a rigid hybrid ceramic/FR4 Micromegas board for improved drift gap uniformity was developed. Initial measurements with 80 GeV/c muons showed improvements in timing response over measured pads and a time resolution below 25 ps. More recent measurements with a new thinner drift gap detector module and newly developed RF pulse amplifiers show that the resolution can be enhanced to a level of 17 ps. This work will present the development of the detector from structural simulations, design, and beam test commissioning with a focus on the timing performance of a thinner drift gap detector module in combination with new electronics using an automated timing scan method.

K

: Micropattern gaseous detectors (MSGC, GEM, THGEM, RETHGEM, MHSP, MI-CROPIC, MICROMEGAS, InGrid, etc), Timing detectors, Cherenkov detectors 1

Introduction

One of the many challenges for future physics experiments is the development of precise timing detectors (O 100 ps) with a large area coverage and capabilities of withstanding high radiation fluence [START_REF] Colaleo | The 2021 ECFA Detector Research and Development Roadmap[END_REF]. PICOSEC Micromegas (MM) is a gaseous detector for precise timing applications with a time resolution of tens of ps [START_REF] Bortfeldt | for the PICOSEC Micromegas Collaboration PICOSEC: Charged particle timing at sub-25 picosecond precision with a Micromegas based detector[END_REF]. It consists of a MM amplification structure that is coupled to a Cherenkov radiator with a CsI semi-transparent photo-cathode. The detector is operating as a two-stage amplification structure with a working gas mixture of Ne:C 2 H 6 :CF 4 (80:10:10). Operating principles are in detail described in [START_REF] Bortfeldt | for the PICOSEC Micromegas Collaboration PICOSEC: Charged particle timing at sub-25 picosecond precision with a Micromegas based detector[END_REF][START_REF] Aune | for the PICOSEC Micromegas Collaboration, Timing performance of a multi-pad PICOSEC-Micromegas detector prototype[END_REF]. A single channel prototype achieved excellent time resolution below 25 ps for measurements with Minimum Ionizing Particles (MIPs) [START_REF] Bortfeldt | for the PICOSEC Micromegas Collaboration PICOSEC: Charged particle timing at sub-25 picosecond precision with a Micromegas based detector[END_REF] and developments continued towards a multi-channel detector that can cover larger areas [START_REF] Aune | for the PICOSEC Micromegas Collaboration, Timing performance of a multi-pad PICOSEC-Micromegas detector prototype[END_REF]. An improved detector prototype with 19 hexagonal pads was built and tested. It was observed that this prototype suffered from the drift field non-uniformity. The effect was investigated and flatness corrections were made based on a particle hit position to obtain a uniform time response. After applying the necessary corrections, the detector yielded a time resolution of 25 ps in the best case [START_REF] Aune | for the PICOSEC Micromegas Collaboration, Timing performance of a multi-pad PICOSEC-Micromegas detector prototype[END_REF]. It was found that the MM Printed Circuit Board (PCB) had initial deformations, most likely caused by the production process and the stretching of the mesh, which can be exacerbated by the tightening of the PCB to the detector housing. In order to produce an even larger multipad PICOSEC, new design procedures had to be introduced.

Mechanical aspects and timing performance of the 100 channel prototype

The design of the new 100-ch prototype was driven with improvements in the flatness and uniformity of the drift gap as the main objectives. Structural mechanics simulations of the mesh tension interface and coupling to the detector casing were performed to select the substrate that could ensure planarity better than 10 𝜇m over 100 cm 2 active area, Figure 1. The Finite Element Analysis suggested the use of a more rigid and thicker PCB. The MM board was fabricated as a combination of 4 x 1 mm thick ceramic core with a 0.3 mm thick FR4 outer layers [START_REF] Utrobicic | for the PICOSEC Micromegas Collaboration Mechanical aspects of CERN GDD 10cm x 10cm new PICOSEC detector[END_REF]. The ceramic core was used to achieve rigidity, while the use of the thin FR4 layers enabled the classic PCB manufacturing process, precise flattening, and polishing. The readout structure was arranged in a 10 by 10 matrix of square pads.

To prevent bending of the MM board or the flange with a Cherenkov radiator when pressed with screws, it was crucial to mechanically detach all the active detector components, including the MM board and the Cherenkov radiator, from the housing. The design choice was to use an Outer Board (OB) on the side of the chamber where the MM board is located, and a flange with a quartz window on the other side where the Cerenkov radiator is located, Figure 2 (left). The signal from each pad was routed to the SMB connectors on the top of the OB via spring-loaded pins. In addition to the electrical connection, the compressive force of the spring-loaded pins ensures that the radiator is pressed firmly and evenly into the positioning groove in the chamber. The time response of the 100-channel Detector Under Test (DUT) to 80 GeV muons was measured at the CERN SPS H4 secondary beamline. The setup consisted of a triple GEM detector telescope for a particle hit position information and MCP-PMTs (later in the text MCP) used as a timing reference or trigger [START_REF] Bortfeldt | for the PICOSEC Micromegas Collaboration PICOSEC: Charged particle timing at sub-25 picosecond precision with a Micromegas based detector[END_REF][START_REF] Aune | for the PICOSEC Micromegas Collaboration, Timing performance of a multi-pad PICOSEC-Micromegas detector prototype[END_REF]. In a standard readout chain configuration, an induced signal from Hamamatsu MCP-PMT R3809U-50 https://www.hamamatsu.com/jp/en/product/type/R3809U-50/index.html a detector is amplified using high bandwidth CIVIDEC C2 amplifier and digitized at a sampling rate of 10 GS/s by LeCroy WR8104 oscilloscope. The same oscilloscope is used for the acquisition of the reference time detector signal and serial bitstream with event ID obtained from the tracker. The recorded waveforms are analyzed using the methods described in previous publications [START_REF] Bortfeldt | for the PICOSEC Micromegas Collaboration PICOSEC: Charged particle timing at sub-25 picosecond precision with a Micromegas based detector[END_REF][START_REF] Aune | for the PICOSEC Micromegas Collaboration, Timing performance of a multi-pad PICOSEC-Micromegas detector prototype[END_REF].

The first measurements were made with a prototype detector with a "standard" drift gap thickness of 220 𝜇m and a CsI photocathode on a thin Cr layer. Two MCP detectors were aligned with each other in the x-y plane and the DUT was positioned behind them. To be able to take data from particles crossing an entire pad area, one MCP was used as a trigger detector. It was mounted parallel to the DUT and placed in front of it. In addition, another MCP, aligned with the trigger MCP in the x-y plane, was placed in front of it and used as a timing reference. For the timing measurements, the DUT was fixed so that the measured pad (or detector area) was aligned with the trigger and timing MCP. Measurements of more than 10 pads showed that a time resolution below 25 ps can be obtained for each pad in its central square area of 5 mm x 5 mm, thus reproducing timing properties of the single channel prototype. The timing of the detector can be further improved if the dependence of the signal arrival time (SAT) on the electron peak (e-peak) charge is taken into account as a time-walk correction. If the drift field is uniform over the entire surface of the pad, the dependence of the timing characteristics on the e-peak charge should be the same for all pads. One of the most important requirements for a large area PICOSEC MM detector is that all pads have a similar dependence of the mean SAT and time resolution on the e-peak charge. This would suggest that it is possible to use only a single global time walk correction function for all the pads, making it easier to calibrate the detector. To investigate the timing performance over a larger area, timing measurements were taken from seven different pads in the same row, Figure 2 (right). The timing analysis shows that all pads closely follow the time-walk correction function, indicating a consistent global behaviour that is important for many-channel detectors, Figure 3 [4].

Method for timing characterization of large area detectors

Measuring the timing performance of large-area detectors presents a challenge for the selection of the reference detector, which must have better timing characteristics and spatial uniformity than the DUT. Although large-area MCPs are available, they are more expensive and the timing performance over the entire area must be known to ensure that the DUT response measurements are not biased. In the previous work [START_REF] Bortfeldt | for the PICOSEC Micromegas Collaboration, Timing performance of a micro-channel-plate photomultiplier tube[END_REF][START_REF] Sohl | for the PICOSEC Micromegas Collaboration Spatial time resolution of MCP-PMTs[END_REF], the timing of the reference detector MCP was investigated as a function of the particle hit position at different radii from the centre of the detector. It was observed that the MCP can be used as a timing reference over the effective photocathode diameter (11 mm), with a time resolution below 6 ps. In the outer regions, the time resolution deteriorates rapidly. To conduct a measurement with minimal bias caused by the non-uniformity of the reference detector over larger DUT areas (> 11 mm), a modification to the measurement method had to be made.

The analysis of the measurements at the reference MCP taken with MIPs shows a doublepeak amplitude spectrum, Figure 4 (left). The correlation of the amplitude spectrum with the hit positions obtained using the tracker shows that the events in the lower signal amplitude peak belong to the outer region of MCP (Figure 4, middle), while those in the larger amplitude peak correspond to the hits in the center of the MCP, Figure 4 (left). The method is based on using the central ("good") part of the MCP to perform timing measurements. This was realized by using the same MCP as a trigger and as a time reference with a sufficiently high threshold to reject the events from the outer region. The MCP was mounted to the movable stage and an automated scanning method was employed to cover the large DUT area with the central part of MCP. The movement was synchronised with the number of events recorded to ensure an even amount of data over the entire area. The MCP was moved on a grid of 2.5 mm x 2.5 mm every 1000 recorded events. A dense grid of movement ensured an overlap of about 75 %, which is beneficial to average possible imperfections of the MCP. Although the method is time-consuming due to the geometric scan, it benefits from the reduced load on the tracking system and oscilloscope recording. This approach can be employed for measurements on a single pad, signal sharing, or entire detector scan studies.

Beam test measurements with custom electronics and 180 𝜇m drift gap thickness

Earlier measurements and simulations indicated that the reduction of the drift gap improves the time resolution of the PICOSEC MM [START_REF] Tzamarias | for the PICOSEC Micromegas Collaboration, Modeling the timing characteristics of the PICOSEC Micromegas detector[END_REF][START_REF] Sohl | for the PICOSEC Micromegas Collaboration, Single photoelectron time resolution studies of the PICOSEC-Micromegas detector[END_REF]. After the initial tests described in section 2 showed the preservation of the timing performance of the single-channel detector and confirmed the uniformity of the detector, the drift gap was reduced to 180 𝜇m. In addition, 50 channels of the detector were fitted with custom-made 10-channel, 38 dB, 650 MHz amplifier cards with built-in discharge protection. The RF amplifier from [START_REF] Hoarau | RF pulse amplifier for CVD-diamond particle detectors[END_REF] served as the basis for the design. The signals were digitised and recorded with oscilloscopes. Measurements on this detector and electronics configuration were performed during two RD51 beam test campaigns in Jul. and Oct. 2022. In Jul., two sets of measurements were taken over a four-pad area, while in Oct. a six-pad measurement was performed. The motivation was to study the single-pad time response and signal sharing over multiple pads. The RMS of SAT values for tracks passing through the 4 mm x 4 mm pad central area during the scans are summarised in tables 1 and 2. SAT distributions of all measured pads show RMS values of 17.1 ps on average. To check whether the uniformity of the drift field is maintained, the dependence of the SAT and the time resolution on the e-peak charge was also investigated. A similar behaviour to the prototype with the standard drift gap detector is found, as shown in figure 5.

Conclusions

The aim of developing a new 100-channel PICOSEC MM detector with an area of 100 cm 2 was to achieve a flat and uniform drift gap that would ensure uniform timing over a large area. A prototype based on an embedded ceramic MM board confirmed the uniformity and preserved the time resolution of the single-channel detector. The drift gap thickness was reduced to 180 𝜇m and timing tests with a newly developed custom RF amplifier were conducted. A timing method for characterising larger areas was implemented and used for measurements. Scans of 10 different pads with a readout chain based on custom amplifiers and oscilloscopes show excellent time resolution averaging 17 ps for all measured pads and the first results show uniformity of the detector with shorter drift gap. New developments are underway to improve stability and robustness and to develop a full readout chain with a multi-channel SAMPIC digitizer [START_REF] Lisowska | Towards robust PICOSEC Micromegas precise timing detectors[END_REF].

aspects and timing performance of the 100 channel prototype 1 3for timing characterization of large area detectors 4 4test measurements with custom electronics and 180 𝜇m drift gap thickness 4 5

 145 Method Beam Conclusions 5

Figure 1 .

 1 Figure 1. Left: Simulation results of the mesh tension influence on 4 mm thick ceramic. A shear force of 15 N/cm was applied on the coverlay faces directed inward the active area. Right: Simulation result of the spring-loaded pin pressure. A Force of 23 g was applied to 100 pressing points. Ceramics Young's modulus was 320 GPa and the Poisson ratio was 0.23.

Figure 2 .

 2 Figure 2. Left: Exploded view of the detector layout. Right: MM PCB pad layout. Each pad is 9.82 mm x 9.82 mm with 200 𝜇m spacing between the edges of adjacent pads. The blue squares represent the pads used for evaluating the timing performance over the detector area.

Figure 3 .

 3 Figure 3. Mean SAT (left) and time resolution (right) values for different electron peak charge bins [4].

Figure 4 .

 4 Figure 4. Left: Spectrum of the MCP amplitudes (below 1 V in yellow and above 1 V in green). Middle: Hit map of the lower amplitudes (yellow). Right: Hit map of the higher amplitudes (green).

Figure 5 .

 5 Figure 5. Scan measurement of 6 pads: SAT (left) and time resolution (right) for different e-peak charge bins. A similar behaviour is measured for all pads.

Table 1 .

 1 RMS of the SAT distributions for 2 x 4 pad scan

	measurements. Cathode/anode voltage: -465 V/+275 V.
	PAD	23 24 33 34 27 28 37 38
	RMS [ps] 17.4 17.1 17.0 17.9 17.2 16.4 17.1 16.6

Table 2 .

 2 RMS of the SAT distributions for 6 pad scan

	measurements. Cathode/anode voltage: -460 V/+275 V.
	PAD	27 37 28 38 29 39
	RMS [ps] 17.7 17.3 17.1 17.1 16.4 16.9

C2-HV Broadband amplifier (2 GHz, 40 dB), https://cividec.at/electronics-C2-HV.html

Acknowledgments

We acknowledge the financial support of the EP R & D, CERN Strategic Programme on Technologies 168 for Future Experiments; the RD51 collaboration, in the framework of RD51 common projects; the 169 Cross-Disciplinary Program on Instrumentation and Detection of CEA, the French Alternative 170 Energies and Atomic Energy Commission; the PHENIICS Doctoral School Program of Université 171 Paris-Saclay, France; the Fundamental Research Funds for the Central Universities of China; 172 the Program of National Natural Science Foundation of China (grant number 11935014); the 173 COFUND-FP-CERN-2014 program (grant number 665779); the Fundação para a Ciência e a 174 Tecnologia (FCT), Portugal (CERN/FIS-PAR/0005/2021); the Enhanced 175 Eurotalents program (PCOFUND-GA-2013-600382); the US CMS program under DOE contract 176 No. DE-AC02-07CH11359.


Alexander Mamontov

Liang Chang

Héloïse Dossmann

Benoît Bertrand

Luc Dechoux

Serge Thorimbert

Iron Catalyzed Dearomatization of Pyridines into Annelated Azepine Derivatives in a One-Step, Three-Component Reaction

come

Iron Catalyzed Dearomatization of Pyridines into Annelated Azepine Derivatives in a One-step, Three-component Reaction

Alexander Mamontov, ‡ Liang Chang, ‡ † Héloïse Dossmann, Benoît Bertrand, Luc Dechoux, and Serge Thorimbert* Sorbonne Université, CNRS, Institut Parisien de Chimie Moléculaire (IPCM), F-75005 Paris, France. Supporting Information Placeholder ABSTRACT: Commercially available Fe(TTP)Cl catalyzes three-component dearomative formal cycloaddition reactions between pyridines, diazo compounds and coumalates. Diversely substituted annelated seven-membered N-heterocycles could be generated in less than ten minutes, in one step at room temperature. The reaction is compatible to gram scale. The extension to benzimidazoles in place of pyridines has been successfully demonstrated. The mechanism of this reaction has been carefully examined by computational studies that corroborate the observed regioselectivities.

Organic chemists are continuously discovering and improving methodologies to build complex molecules. Some colleagues from the pharmaceutical industry recently stated: "Drug molecules typically contain amines and N-heterocycles as well as unprotected polar groups, which present significant synthesis challenges. More synthetic methods are needed which tolerate these functionalities". [START_REF] Blakemore | Organic Synthesis Provides Opportunities to Transform Drug Discovery[END_REF] In the expending world of azaheterocyclic compounds, 2 (polycyclic) seven-membered azepan and azepine derivatives are key structural fragments found in natural products isolated from fungi, plants and marine sponges (Figure 1). [START_REF] Ryan | Chapter 7 -Seven-Membered Rings[END_REF] These motifs possess a wide range of biological activities leading to various applications in medicine and pharmacology. 4 As a result, many approaches have been developed to construct the 7-membered core, 5 among which nucleophilic substitution, lactamisation, ring reorganization, ring closing metathesis, cycloadditions, including the use of azomethine ylides. 6

Figure 1. Natural Products Containing Fused Azepine Moiety

In comparison to other 7-membered N-heterocycles, dihydropyrido [1,2-α]azepines synthesis remains underexplored and efficient and general approaches are still to be discovered. [START_REF] Fischer | The Chemistry of Pyrido[1,2-a]Azepines and Their Hydro Derivatives[END_REF] The first but unsuccessful approach toward pyrido [1,2-α]azepines was reported by Jones. 8 In a more recent approach, Xia prepared a pyridinium [1,2-α]azepine ruthenium complex that only led to the azepino-pyridinium salt after deprotection with tetrafluoroboric acid. [START_REF] Zhou | A Missing Member of Conjugated N-Heterocycles: Realizing Pyrido[1,2-α]Azepine by Reacting Ruthenium Alkenylcarbene Complex with Alkyne[END_REF] Eberbach and co-workers reported and studied the 8-electrocyclization of conjugated butadienylpyridinium ylides into dihydropyrido [1,2-α]azepines (Figure 2a). Altogether, the existing methods require multiple steps to access the desired compounds, with moderate yields or difficulties to isolate intermediates. Therefore, searching for an efficient, easily adaptable method to synthesize dihydropyrido [1,2-α]azepines is still highly desirable.

Azomethine ylides have been used as C-N-C partners to build N-heterocycles. [START_REF] Wei | Catalytic Asymmetric Reactions with N-Metallated Azomethine Ylides[END_REF]11 Among others, dihydroazepines have also been synthesized from iminoesters and pyrones (Figure 2 b). [START_REF] Liu | Et3N-Catalyzed Tandem Formal [4 + 3] Annulation/Decarboxylation/Isomerization of Methyl Coumalate with Imine Esters: Access to Functionalized Azepine Derivatives[END_REF] Pyridinium ylides which are stable dipoles are known to be versatile partners in [4+3] cycloadditions. 13 The possible salt-free preparation of pyridinium ylides by in situ reaction of pyridines with diazo compounds in the presence of metal catalyst is even more interesting in terms of sustainability. Methyl coumalate, a 4-carbomethoxy pyrone, has demonstrated a wide range of partner-dependent reactivities. In addition to photoinduced intramolecular [2+2] electrocyclization leading to cyclobutenyl synthons, 15 it could deliver 2 π-or 4 π-electrons in cycloadditions with dienes, dienophiles, or dipoles. 16 It is an excellent electrophile that participates in cascade reactions 17 or can be converted into a nucleophilic species via the transient addition of organocatalysts. 18 As a cherry on the cake, coumalic acid is produced from the dimerization of malic acid, 19 providing a basis for sustainable synthesis of more complex molecules.

In this context, we developed a simple and general method to access diversely functionalized annelated azepines in a one-step multicomponent reaction, starting from commercially available compounds (Figure 2c). To initiate the study, pyridine 1a, ethyl diazoacetate 2a and methyl coumalate 3a were mixed in DCM in the presence of 1 mol% of 5,10,15,20-Tetraphenyl-21H,23H-porphyrin iron (III) chloride (Fe(TTP)Cl) as catalyst. To our delight, the desired dihydropyridoazepine 4a, was isolated in 73% yield after 10 minutes (Table 1, entry 1). No other product was detected. A concentration-dependance study was carried out and 0.33 M seems the optimal concentration. We also divided per two or doubled the c) This study : multicomponents approach to dihydropyrido [1,2]azepines

N R EWG C4 synthon O O CO 2 R' + EWG N 2 + N R CO 2 R' N O O CO 2 Me + 3a N EtO 2 C EtO 2 C N 2 catalyst + CO 2 Me DCM / Time Temperature 2a 1a 4a
amount of catalyst without improvement (see SI, Table S1). Substitution of DCM by other solvents did not provide better results. Toluene or THF afforded compound 4a in yield reaching 70% whereas EtOH led to 63% of 4a.

Different metal salts were investigated and selected results are summarized in Table 1 (see SI, Table S1). Various iron-based catalysts such as 5,10,15,20-Tetrakis(4-methoxyphenyl)-21H,23H-porphyrin iron (III) chloride (Fe(TTPOMe)Cl) and 5,10,15,20-Tetrakis(pentafluorophenyl)-21H,23H-porphyrin iron (III) chloride (Fe(TTPF)Cl) or FeCl 3 were tested. While similar results were observed with Fe(TTP)Cl and Fe(TTPOMe)Cl, the use of Fe(TTPF)Cl reduced the yield in 4a down to 58% (Table 1, entries 1-4). Unfortunately, FeCl 3 did not lead to the formation of 4a even with 1 h reaction time or increasing its quantity up to 10 mol% (see SI). Initial results with copper diacetate (II) were encouraging, giving 4a in a moderate 23% yield, however all attempts to improve the conversion were unsuccessful (Table 1, entries 5-7). The use of catalysts based on palladium, rhodium, silver or chromium did not allow the formation of the desired compound even with longer reaction time (See SI, Table S1).

In summary, we chose the following conditions: pyridine, diazoacetate 2a (1 equiv), coumalate 3a (1.2 equiv.) and Fe(TTP)Cl (1 mol %) in DCM and tested at the gram scale. It delivered 2.1 g (61% yield) of the desired dihydropyridoazepine 4a in less than 10 min (Table 1, entry 2).

Starting from commercially available benzyl 2-diazoacetate 2b and tert-butyl 2-diazoacetate 2c (Table 2), the desired dihydropyridoazepines 4b and 4c were obtained in moderate yield, whereas diazo amides 2d-f led to the expected fused azepines 4d, 4e and 4f in yield up to 67%.

Pyridines 1g-k with substituents in para position were first tested for their symmetry and low steric hindrance. Isopropyl or tertbutyl groups had no significant influence and the desired compounds 4g-h were obtained in yields around 65 %. The use of 4phenylpyridine 1i as substrate led to the corresponding dihydropyridoazepine 4i in 63%. However, the presence of the strongly electron-withdrawing trifluoromethyl group significantly decreased the reactivity of 1j. The same trend was observed with 4-cyano group as no reaction occurred. Gratifyingly, the dihydropyridoazepine 4k bearing a 4-methoxy group could be isolated in 83% yield after precipitation. Altogether, these results suggest that electronic density on the nitrogen atom of the nucleophilic pyridine has a strong influence on the overall reactivity and electron-neutral or donating substituents are preferred. Although the conversion of the reactants was sometimes almost quantitative, unfortunately the stability of the products should be mentioned as a drawback of this efficient dearomative reaction. For example, 4k degraded over silica or alumina column chromatography.

Table 2. Scope and Limitations

R n R R'OC CO 2 R"
Subsequently, we tested pyridines 1l-q substituted in meta position. Interestingly, products 4l-q were obtained as single regioisomers. Indeed, although positions 2 or 6 of pyridines 1l-q could be engaged to fuse the rings, the cyclization took place exclusively at position 2. Yields for products presenting a methyl or a methoxy group (4l-m) are comparable to the one of 4q bearing a fluorine atom (60 to 71%). Electronic effects seem to be less pronounced than steric ones. Indeed, from iodide to fluoride (compare yields from 4n to 4q), the smaller the atomic radius of the halogen is, the higher the isolated yield. These results could be attributed to steric hindrance but may also be explained by the possible insertion of iron into the carbon-halogen bond. When 3-(trifluoromethyl)pyridine was used, no bicyclic product could be isolated. To complete the study, we checked ortho-substituted pyridines. Unfortunately, 2-methyl-and 2-fluoro pyridine did not react at all, proving the strong steric influence appearing during the cyclisation.

We turned our attention toward polysubstituted pyridines. We first investigated symmetrically disubstituted ones such as 3,5dimethyl-and 3,5-difluoropyridine. The corresponding dihydropyridoazepines 4r and 4s were isolated in low yields. These results were attributed to steric hindrance in the case of 1r and to a strong decrease of nucleophilicity for 1s. Positively, the 3-bromo-4methylpyridine 1t led to the formation of dihydropyridoazepine 4t in 65 % yield with a total regioselectivity. Finally, we varied both pyridines and diazo compounds, the desired dihydropyridoazepines 4bt to 4fm were isolated as expected.

The third reactant of this three-component reaction, namely pyrone 3, could also be varied by changing the nature of the ester, introducing various alkyls or aryls instead of methyl group (see Table 2). As representatives, coumalates bearing allylic, propargylic or phenolic ester function were prepared from coumalic acid and the corresponding alcohols using EDC. The coumalates 3u-x reacted smoothly giving the expected dihydroazepines 4u-x in yields from 67 to 76%. The reaction supported the double functionalization on both pyridine 1 and coumalate 3 leading to 4uq and 4vm in 61 and 70% yields respectively.

We also investigated the replacement of pyridine 1 by benzimidazole 5. Using the previously optimized conditions, compound 6a resulting from both N-H insertion and cyclization reaction was isolated with a modest 23% yield together with unreacted 3a. In order to fully consume starting materials, the quantity of ethyl diazoacetate 2a was increased to 2.5 equivalents. In these conditions, the tricyclic fused azepine 6a was isolated in 65% yield (Table 3).

Subsequently, we explored the influence of the substituents on the regioselectivity and efficiency of the reaction. The 5-bromo-, and 5-methoxybenzimidazole 5b-c were tested as representatives. The reaction proceeded well leading to products 6b-c as an almost equimolar mixture of two regioisomers. Brominated tricyclic benzimidazolo azepine 6b was obtained with a slightly lower yield compared to the methoxy substituted product 6c (52 vs 76 % yields respectively). The N-methylbenzimidazole 5d required only one equivalent of 2a and reacted cleanly to provide the expected tricyclic fused azepine in 63% isolated yield (Table 3).

Table 3. Scope and Limitations in Benzimidazole Series a,b

a Reaction conditions : benzimidazoles 5 (1 mmol), 2a (2.5 equiv.), 3a (1 equiv.), Fe(TTP)Cl (1 mol%), r.t., DCM (3 mL). b Isolated yields. c mixture of regioisomers -6b 2/3 ratio : 47/53 -6c 2/3 ratio : 43/57. d 1 equiv. of 2a.

Concerning the mechanism of the reaction, a pyridinium ylide is stated to be catalytically produced from the reaction between 1a and 2a in the presence of Fe(TPP)Cl. [43][44][45][46][47][48] Then, two plausible pathways could be envisaged for the observed cyclization between 3a and the pyridinium ylide, i.e. a concerted [4+3] cycloaddition or a two-step asynchronous addition-cyclization.

Based on previous reactivity observed for 3a, we privileged a two-step process. On the other hand, the total regioselectivity observed with 3-substituted pyridines preferentially supported a concerted mechanism, or at least a pre-organized association between the two reactants. To discriminate these two pathways, we carried out DFT calculation (Figure 3, see also SI, Table S2).

A preassociative dipole-dipole complex A linking the substrates was found to align the two separate polar molecules (+ 6.8 kcal/mol to access A). The first C-C bond is created through a transition state lying at 20.7 kcal/mol above reactants. This 1,6addition leads to an enolate B in an endothermic reaction (8 kcal/mol). It rapidly undergoes an intramolecular endothermic C-C bond formation leading to the lactone-bridged nitrogen tricyclic compound C that presents a strain of 28.3 kcal/mol. Surprisingly, the concerted [4+2] retro Diels-Alder reaction that expulses CO 2 , requires the highest energy to convert C into azepine D but releases 22.4 kcal/mol by decongestion of the structure. Intermediate D stands at 0.9 kcal/mol lower than complex A. The last step is the formation of the highly conjugated dihydropyridoazepine 4a via a 1-5 hydrogen shift. This step which maximizes -electrons conjugation is the driving force of the observed reaction which releases 9.5 kcal/mol. The total regioselectivity observed with 3-substituted pyridines, leading to products 4m-r, could be justified by the development of steric hindrance that appears during the initial approach of the two reacting molecules (see Figure S1 andS2). The absence of reactivity observed for substituted 2-pyridines could similarly be explained by the stronger steric interactions that appear for such substrates.

N N EtO 2 C CO 2 Me N R R 1 6d, R = H ; R 1 = Me ; 63% d 6a, R = H ; R 1 = CH 2 CO 2 Et ; 65% 6b, R = Br ; R 1 = CH 2 CO 2 Et ; 52% c 5d, R = H ; R 1 = Me 5a, R = H ; R 1 = H 5b, R = Br ; R 1 = H 6c, R = OMe ; R 1 = CH 2 CO 2 Et ; 76% c 5c, R = OMe ; R 1 = H R 1 2 3 R

Figure 3. Computed Pathway for the Formation of 4a

In summary, we have developed a multicomponent domino reaction to quickly access fused bicyclic [6+7] dihydropyridoazepines or tricyclic [6+5+7] fused benzimidazolo azepines starting from commercially available substrates, renewable feedstock and Fe(TTP)Cl as catalyst. A wide range of pyridines, benzimidazoles, diazo compounds and coumalates were used and more than 35 compounds were successfully synthesized. A total regioselectivity was observed for 3-substituted pyridines. The proposed mechanism was rationalized by DFT calculations. This dearomatization process leads to the formation of molecules previously reported as difficult to synthesize. The reaction appeared scalable up to the gram scale and releases only CO 2 and N 2 as by products.

ASSOCIATED CONTENT ACKNOWLEDGMENT

We thank Sorbonne Université and CNRS for funding, the Fédération de Recherche de Chimie Moléculaire de Paris Centre (FR2769) for providing technical access to analytical instruments. The calculation was granted access to the HPC resources of the SACADO MeSU platform at Sorbonne-Université. L.C. thanks the Chinese Scholarship Council (CSC) for a PhD grant.

  14

Figure 2 .

 2 Figure 2. Pathways to Pyrido[1,2-α]azepines

  b Yield after purification by flash-chromatography. c Fe(TTP)Cl : 5,10,15,20-Tetraphenyl-21H,23H-porphyrin iron (III) chloride d Reaction on 1 gram scale e Fe(TTPOMe)Cl : 5,10,15,20-Tetrakis(4-methoxyphenyl)-21H,23H-porphyrin iron (III) chloride f Fe(TTPF)Cl : 5,10,15,20-Tetrakis(pentafluorophenyl)-21H,23H-porphyrin iron (III) chloride. g Reaction in dichloroethane

a

  Reaction conditions : pyridines 1 (1 mmol), diazo compounds 2 (1.2 mmol), coumalates 3 (1 mmol), Fe(TTP)Cl (1 mol%), DCM (3 mL), r.t., 10 min. b Isolated yields

  10

Table 1 . Optimization of the Reaction Conditions a

 1

	Entry	Catalyst (mol%)	Temp.	Time	4a (%) b
	1	Fe(TTP)Cl (1) c	r.t	10 min	73
	2	Fe(TTP)Cl (1)	r.t	10 min	61 d
	3	Fe(TTPOMe)Cl (1) e	r.t	10 min	73
	4	Fe(TTPF)Cl (1) f	r.t	10 min	58
	5	Cu(OAc) 2 (10)	r.t	4 h	23
	6	Cu(OAc) 2 (10)	r.t	12h	27
	7	Cu(OAc) 2 (10)	50 °C	4 h	22 g

a All reactions were carried out with pyridine 1a (1 mmol), 2a (1.2 mmol), 3a (1 mmol), Fe(TTP)Cl (1 mol%), DCM (3 mL), r.t. 10 min.

Data Availability Statement

The data underlying this study are available in the published article and its Supporting Information

Supporting Information

The Supporting Information is available free of charge on the ACS Publications website. It contains information regarding experimental procedures, optimization table, characterization of organic molecules, and mechanistic studies (PDF).

AUTHOR INFORMATION

Corresponding Author * Serge Thorimbert . serge.thorimbert@sorbonne-universite.fr Orcid : 0000-0002-0290-2820

Authors

Alexander Mamontov : Orcid : 0000-0002-5444-7778 Benoit Bertrand Orcid : 0000-0003-0542-4454 Luc Dechoux Orcid : 0000-0003-3806-1848 Héloise Dossmann Orcid : 0000-0003-4844-229X Present Addresses † Liang Chang : 138 Xianlin Rd, Nanjing, Jiangsu, China, 210023. Orcid 0000-0001-9552-519X

Author Contributions

The manuscript was written through contributions of all authors. / All authors have given approval to the final version of the manuscript. / ‡These authors contributed equally.


P Horgue
email: phorgue@imft.fr

F Renard

G S Gerlero

R Guibert

G Debenest

porousMultiphaseFoam v2107: An open-source tool for modeling saturated/unsaturated water flows and solute transfers at watershed scale

Keywords: Porous media, Unsaturated flow, Richards' equation, Newton's algorithm, OpenFOAM

The first releases of porousMultiphaseFoam proposed an open-source software suite to solve the equations for multiphase flow (generalized Darcy's law) in porous media or groundwater flows (Richards' equation) by taking advantage of OpenFOAM, a finite volume platform with automatic discretization on three-dimensional unstructured grids and good parallel efficiency. Recently, the porousMultiphaseFoam toolbox has been confronted with complex cases of fast water flows and solute transfers in realistic hydrological configurations with variable forcing conditions (heterogeneous infiltration and local tracer injection). Several developments have been carried out to make it possible to simulate those cases, which extend the toolbox with: (i) a set of solvers dedicated to groundwater flows, including coupled water flow and solute transport and simplified 2D approaches, (ii) improved numerical techniques for problems with strong non-linearities, (iii) libraries/executables for pre-processing of input data (geographical information and time-variable forcing terms) and (iv) passive or coupled scalar transport (tracer) with groundwater solvers that support any number of species. New solvers are validated on several (un-)saturated configurations by a direct comparison with a well validated finite element code.

Introduction

The modeling of saturated/unsaturated flow in soils is of importance in a wide range of scientific fields such as environmental monitoring, waste management or hydrogeology. Multiphase flow through porous media can be modeled by solving mass conservation for each phase expressing velocities using the Generalized Darcy's law [START_REF] Muskat | Physical principles of oil production[END_REF]. However, a classical approach for the air-water flow modeling (for hydrology for example) consists in neglecting the air pressure gradients to reduce the problem to the resolution of the water mass conservation, via the so-called Richardson [START_REF] Richardson | Weather prediction by numerical process[END_REF] or Richards' equation [START_REF] Richards | Capillary conduction of liquids through porous mediums[END_REF].

Despite its apparent simplicity, the Richards' equation remains challenging to be accurately solved due to its highly non-linear terms which require iterative procedures such as Picard's or Newton's methods. Moreover, the modeling of hydro-geological events may involve a wide range of time and space scales, inducing strong requirements on space grids or time steps. To reduce the computational cost, it is then necessary to use advanced numerical techniques for time step management and/or adaptive mesh refinement. Main difficulties and recent works related to the numerical solution of the Richards' equation, as well as an inventory of the techniques that have been developed, can be found in two recent reviews [START_REF] Farthing | Numerical Solution of Richards' Equation: A Review of Advances and Challenges[END_REF] and [START_REF] Zha | Review of numerical solution of Richardson-Richards equation for variably saturated flow in soils[END_REF]. Over the last two decades, several open-source tools dedicated to the 3D modeling of Richards' equation have been developed, such as ParFlow [START_REF] Ashby | A parallel multigrid preconditioned conjugate gradient algorithm for groundwater flow simulations[END_REF], PFLOTRAN [START_REF] Hammond | Evaluating the performance of parallel subsurface simulators: An illustrative example with PFLOTRAN[END_REF], Dumux [START_REF] Flemisch | DuMux: DUNE for multi-{phase,component,scale,physics,...} flow and transport in porous media[END_REF] and RichardsFOAM [START_REF] Orgogozo | An open source massively parallel solver for Richards equation: Mechanistic modelling of water fluxes at the watershed scale[END_REF], with special attention brought to the parallel efficiency of the different solvers since the three-dimensional modeling of groundwater at the watershed scale can be particularly expensive. The open-source tools listed previously exhibit parallel linear behavior up to at least one thousand processors, such as the porousMultiphaseFoam (PMF) used in this work [START_REF] Horgue | An open-source toolbox for multiphase flow in porous media[END_REF][START_REF] Horgue | An extension of the open-source porousMultiphaseFoam toolbox dedicated to groundwater flows solving the Richards' equation[END_REF].

The groundwater modeling is often simplified by neglecting the unsaturated upper part and using the Dupuit-Forchheimer assumption to consider that the vertical component of the velocity is zero [START_REF] Bear | Dynamics of fluids in porous media[END_REF]. This groundwater free surface approximation allows then to simulate groundwater flow as a two-dimensional flow by solving the water height without capillary effects. The strong assumptions required by the 2D approach are counterbalanced by the drastic reduction of the computational costs which is particularly effective under specific conditions. A more accurate simplified approach consists in coupling 1D (for the unsaturated part) and 2D (for the saturated part) models to simulate unsaturated and saturated groundwater flow at the watershed scale while keeping low computational costs [START_REF] Twarakavi | Evaluating interactions between groundwater and vadose zone using the hydrus-based flow package for modflow[END_REF][START_REF] Kuznetsov | Quasi 3d modeling of water flow in vadose zone and groundwater[END_REF][START_REF] Renard | A new quasi-3D unsaturated-saturated hydrogeologic model of the Plateau de Saclay (France)[END_REF].

Simultaneously to the fluid flow, the solute transport through porous media has also been widely studied with some specific effects such as the dispersion [START_REF] Golfier | Comparison of theory and experiment for solute transport in highly heterogeneous porous medium[END_REF][START_REF] Urteaga | Transverse solute dispersion in microfluidic paper-based analytical devices (µPADs)[END_REF] (velocity fluctuations at the micro-scale which increase the spreading of the tracer) or the retardation factor [START_REF] Deng | Upscaling retardation factor in hierarchical porous media with multimodal reactive mineral facies[END_REF] (linked to the sorption coefficients). Where the precise characterization of the parameters remains difficult to establish [START_REF] Gan | A comprehensive evaluation of various sensitivity analysis methods: A case study with a hydrological model[END_REF], some usual formulations are commonly accepted in the community and generally implemented in the available numerical tools.

A general open-source toolbox based on OpenFOAM was initially developed for generic two-phase flow [START_REF] Horgue | An open-source toolbox for multiphase flow in porous media[END_REF], including the commonly used soil water retention models (for relative permeability and capillary pressure). The PMF toolbox was first extended to Richards' equation [START_REF] Horgue | An extension of the open-source porousMultiphaseFoam toolbox dedicated to groundwater flows solving the Richards' equation[END_REF], using existing libraries and functions, and providing an efficient tool dedicated to air-water flows. The toolbox has also served as the basis for other codes such as the hybridPorousInterFoam [START_REF] Carrillo | Multiphase flow modeling in multiscale porous media: An open-source micro-continuum approach[END_REF] or for two-phase flow with adaptative mesh refinement [START_REF] Sugumar | Grid Adaptation of Multiphase Fluid Flow Solver in Porous Medium by OpenFOAM[END_REF]. However, it has become apparent that multiple features were missing in order to offer a complete modeling suite such as: (i) A fast two-dimensional free-surface modeling tool, (ii) Passive/coupled tracer solvers for 2D-saturated/unsaturated flows, (iii) Pre-processing tools for realistic configuration with Geographical Information System (GIS), (iv) Improved numerical techniques to allow the modeling of unsaturated complex cases.

The work presented in this paper intends to validate the multiple developments which have been made during the last years. Due to the good parallel efficiency observed for OpenFOAM executables, one will be able to simulate flows and transfers in three-dimensions at the watershed scale.

Mathematical model and governing equations

Mathematical model for saturated flow (2D modeling)

The groundwater flow modeling of the saturated zone can be simplified to a two-dimensional problem by considering a constant head-over-depth assumption (called Dupuit-Forchheimer's assumption). In this configuration, the momentum equation reduces to the Darcy's law and the continuity equation reads :

ε ∂ψ ∂t -∇ • H water Kρ∥g∥ 2 µ ∇ψ = V source , (1)
where

ψ [m] is the groundwater head, ε [-] the porosity, K [m 2 ] the permeability, ρ [kg.m -3
] the density of the fluid, µ [P a.s] the dynamic viscosity of the fluid and g [m.s -2 ] the gravity. The water depth H water [m] is defined as:

H water = ψ -z 0 , (2)
where z 0 [m] is the height of the impermeable bottom of the aquifer. The flow velocity V source [m.s -1 ] is used for rainwater infiltrations or other inflow sources.

The boundary conditions available for 2D-saturated flows are:

-Dirichlet condition to impose potential :

ψ boundary = ψ f ixed ,
-Neumann condition on the potential gradient to impose flow velocity u f ixed . Fixed gradient value is computed as :

∇ψ boundary • - → n = µ H water Kρ∥g∥ 2 u f ixed • - → n (3)

Seepage condition

During medium and high intensity infiltrations, it happens that the local capacity of the aquifer is reached and that any additional infiltration results in a surface runoff which does not further affect the groundwater flow. This physical feature needs to be handled to correctly model the behavior of a watershed over a long period. Each computation cell i gets a maximal value z max,i relative to the local thickness of the aquifer. This additional geometric information is obtained using bilinear interpolation from Digital Elevation Model (DEM) files and represents the local maximum value for the potential ψ i . At any linear iteration k, the potential ψ k i and the flow balance ϕ k i for each cell i can be used to determine if the cell is overflowing and therefore if a fixed potential value should be applied:

if ψ k i ⩾ z max,i and ϕ k i ⩾ 0 then ψ k+1 i = z max,i (4)
The dynamic forcing is applied after each time-step (transient simulation) or solver iteration (steady simulation).

Mathematical model for unsaturated flow (vadose zone)

Modeling variably saturated flow in the vadose zone generally requires solving Richards' equation, a physics-based equation where air pressure gradient has been neglected to keep only the liquid phase conservation equation. Several forms of the Richards' equation can be found in the literature such as the pressure-head form, the saturation-based form and the mixed form proposed by Celia et al. [START_REF] Celia | A general mass-conservative numerical solution for the unsaturated flow equation[END_REF], the latter of which has been implemented in the software. This modified version consists of using the water content θ for the accumulation term and pressure head h for the other terms, which reads:

∂θ ∂t + θS s ε ∂h ∂t -∇ • Kρ∥g∥ 2 µ k r,h ∇h + ∇ • Kρ µ k r,h g = Q source , (5)
where θSs ε ∂h ∂t is the aquifer specific storage term (S s [-] is the specific storage coefficient) and k r (h) [-] is the relative permeability and Q source [m 3 .s -1 ] is a volumic source term. The mixed-formulation has the advantages of being mass conservative (like the θ-based form) as well as being able to handle locally saturated areas (as the h-based form).

This two dependent variables (θ, h) equation cannot be directly linearized, and water content θ has to be expressed as a function of head pressure h by considering the capillary capacity C h defined as:

C h = ∂θ ∂h (6)
The capillary capacity C h [m -1 ] and the relative permeability k r,h [-] can be computed using usual relationships such as Brooks and Corey [START_REF] Brooks | Hydraulic Properties of Porous Media and Their Relation to Drainage Design[END_REF], Van Genuchten [START_REF] Van Genuchten | A Closed-form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils[END_REF] or Ippisch et al. [START_REF] Ippisch | Validity limits for the van Genuchten-Mualem model and implications for parameter estimation and numerical simulation[END_REF].

Boundary conditions available for unsaturated flows are :

-Dirichlet condition to impose pressure head value:

h boundary = h f ixed
-Neumann condition on the pressure-head gradient to impose flow velocity u f ixed . Fixed value is computed as :

∇h boundary • - → n = µ Kρ∥g∥ 2 k r,h u f ixed • - → n (7)
Note that the Richards' equation ( 5) must be solved using an iterative procedure to handle the (possibly strong) non-linearities of the relative permeability k r,h and θ -h relationships (Picard or Newton's algorithm for example, see Section 3.2).

Mathematical model for tracer transport

The general conservation equation for the transport of a tracer with concentration C [-] reads [START_REF] Bear | Modeling groundwater flow and contaminant transport[END_REF][START_REF] Bear | Dynamics of fluids in porous media[END_REF] :

R ∂ΘC ∂t + ∇ • U w C -Θ∇ • (D ef f ∇C) = -Q C -ΘRλC, (8)
where λ [s -1 ] is the decay coefficient, U w the water flux equals to:

U w = U [m.s -1 ]
for unsaturated solvers

U × H water [m 2 .s -1 ] for 2D-saturated solvers , (9)
Θ equals to:

Θ = θ [-] for unsaturated solvers ε H water [m] for 2D-saturated solvers , (10)
and R[-] the retardation factor :

R = 1 + (1-ε total )ρsK d θ [-] for unsaturated solvers 1 + (1-ε total )ρsK d ε [-] for 2D-saturated solvers , (11)
where ε total [-] is the total porosity (larger than kinematic porosity ε), ρ s [kg.m -3 ] is the solid matrix density and [START_REF] Fried | Dispersion in porous media[END_REF][START_REF] Bear | Dynamics of fluids in porous media[END_REF]:

K d [m 3 .kg -1 ] the volumetric partitioning coefficient. The effective dispersion coefficient D ef f [m 2 .s -1 ] is computed using local values of velocity U [m.s -1 ] field as
D ef f = D m τ + α|U | ε (12)
where D m [m 2 .s -1 ] is the molecular diffusion coefficient, τ [-] the tortuosity and α [m] the dispersion tensor constructed using velocity direction and α L /α T [m] longitudinal and transversal dispersion coefficients. An equation taking the form of (8) has to be solved for each tracer species transported by the fluid. Q c stands for source term whose dimension depends on the solver considered: [s -1 ] for unsaturated and [m.s -1 ] for 2D-saturated solvers.

Numerical methods

Equation discretization

The toolbox has been developed using the OpenFOAM environment which is a development framework for solving equations using the finite volume method. The OpenFOAM native functions allow one to automatically discretize classical operators ( ∂ ∂t , ∇, ∆) on unstructured three-dimensional meshes using a wide range of numerical schemes. In the validation cases studied in this work, a classical second order centered scheme is used for spatial discretization of convective and diffusive terms. A second order backward scheme is used for time discretization of tracer transport and saturated flow equations. The time discretization of Richards' equation (unsaturated flow) is carried out using a backward Euler scheme to simplify the implementation of non-linear Picard and Newton's methods (cf Section 3.2).

Non-linear methods for Richards' equation

To linearize the Richards' equation [START_REF] Carrillo | Multiphase flow modeling in multiscale porous media: An open-source micro-continuum approach[END_REF], it is necessary to express the implicit water content θ at iteration k + 1 as a function of pressure head h using capillary capacity C h (cf eq. 6) as:

θ k+1 = θ k + C k h h k+1 -h k , (13)
where θ k and h k are the current iteration values (C k h computed using h k ). The derivative ∂θ ∂t can then be discretized using a backward Euler scheme:

∂θ ∂t ≈ θ k+1 -θ ∆t = 1 ∆t C k h h k+1 -h k + θ k -θ ( 14
)
where θ is the water content at the previous time step. This makes possible to express the linearized form of conservative equation 5 as a function of solved variable h k+1 and current iteration values h k and θ k :

1 ∆t C k h h k+1 -h k + θ k -θ -∇ • Kρ∥g∥ 2 µ k r,h k ∇h k+1 + ∇ • Kρ µ k r,h k g = Q source ( 15
)
This linearized mixed form of the Richards' equation ( 15) can be rewritten as:

C k h k+1 -h ∆t -M k h k+1 + L k = S k (16)
with the diagonal matrix C k containing the accumulative terms, h the old-time pressure head vector, M k head-pressure gradient-related terms, L k gravity-related terms and S k the source-term vector. Picard's algorithm consists of solving the equation 16 iteratively while updating the matrix coefficients using the pressure-head vector h k at iteration k. Letting R k be the equation residual at iteration k:

R k = C k h k -h ∆t -M k h k -L k -S k (17)
This method converges linearly to the solution when:

|R k | < r tol ( 18
)
with r tol the user-defined residual tolerance. This method was the one initially implemented in PMF, yet it may converge slowly for highly non-linear configurations and/or low tolerance values. It is particularly unsuitable for the modeling of high flow rates in the vadose zone (see test case in section 5.1.4) where the time-steps required for the convergence of this algorithm imply a prohibitive computational cost. However, the method is robust and has a wide range of convergence.

To improve the convergence rate and reduce the computational cost, one can use Newton's method, which consists of solving the equation:

J k h k+1 -h k = -R k (19)
where J k is the Jacobian matrix of the residual function R which is equal (row i, column j):

J ij = ∂R i ∂h j (20)
This method converges to the solution quadratically and thus allows solving highly non-linear problems, as shown in section 5.1.4.

Both methods are used sequentially in unsaturated solvers (groundwaterFoam and groundwaterTrans-portFoam), by solving the equation first using Picard's method until |R k | < r tol,P icard and then refining the solution with Newton's method until |R k | < r tol,N ewton . This approach used by several numerical codes, including METIS, takes the advantage of the wide convergence range of Picard's algorithm to approach the solution roughly and then converges quadratically with Newton's method.

Time step management

An efficient time step management is required for modeling transient saturated/unsaturated water flow because the characteristic time scales can vary by several orders of magnitude. A short-time and high amplitude water injection rate into an unsaturated column for example would require small time-steps at the beginning (∆t ≈ 10 seconds for several hours of physical time) while the following capillary drainage should be simulated over a long period of several months (∆t ≈ 1 day or more). Because of the non-linearity of Richards' equation, a general CFL condition cannot be directly derived. Two methods for time step management are generally used:

(i) the heuristic one is based on the number of linear iterations (Picard or Newton) required at each time step, (ii) the adaptive one is based on the time truncation error estimates (related to the numerical scheme used for time discretization).

D'Haese et al. [START_REF] D'haese | Assessment of adaptive and heuristic time stepping for variably saturated flow[END_REF] have proposed a comparison study of these two methods and their variants by analyzing the advantages and limitations of both approaches which are case-dependent. In this work, time stepping is based on a rigorous mathematical analysis performed on the numerical scheme used for time discretization which can be applied to all implemented solvers (2D-saturated, unsaturated, passive transport or coupled). The time step is computed to maintain a relative time truncation error (estimated using Taylor series) below a user-defined value ϵ error . The time-step computations for the available time schemes in the presented code are :

1. backward Euler: ∆t = 2ϵerrorxmax ( ∂ 2 x /∂t 2 ) max 2. 2 nd order backward scheme:∆t = 3 4ϵerrorxmax ( ∂ 3 x /∂t 3 ) max 3. Crank-Nicolson: ∆t = 3 12ϵerrorxmax ( ∂ 3 x /∂t 3 ) max
where x max is the value of the variable at the cell where the maximal derivative has been computed. A third order time scheme had to be implemented in the toolbox to allow adaptive time step for available time schemes.

Accessory numerical tools

Toolbox file structure

The file structure of the toolbox is depicted in Figure 1. One can see old libraries or solvers developed for PMF in 2015 and the features added in the new version.

eventFile structure

To help the set up of realistic cases with time-varying conditions, a generic OpenFOAM class eventFile has been developed with three major concrete subclasses available :

-eventFileInfiltration for homogeneous or heterogeneous infiltration (2D-saturated solvers), -patchEventFile for solute/water injection on boundaries (unsaturated solvers), -sourceEventFile for local water or solute injection (all solvers).

For each time iteration, forcing terms are automatically computed from linear interpolation between the two closest dates. An example of sourceEventFile is given in Figure 2.

An additional concrete class outputEventFile allows specifying output writing times while avoiding any constraints on the time step, with the solutions being constructed by time interpolation. Extension .evt has been chosen to identify the files.

Pre-processing tools

OpenFOAM executables have been developed to convert GIS data into informations for OpenFOAM simulations:

-setFieldsFromDEM: a variant of the official OpenFOAM setFields tool which performs bilinear interpolation from DEM file in ASCII format (.xyz). This tool is generally used to set up the local height of the aquifer (z min and z max ).

-setFieldsFromXY: similar to setFieldsFromDEM with an interpolation performed using the three closest points of any list of points. This tool is used in this work to set up a permeability field with non-uniform values.

-setBoundaryHeadPressure: utility to set up boundary conditions for unsaturated (h) or 2D-saturated (ψ) solvers using a uniform water height or using a DEM or STL input file.

A simple Darcy utility has also been added (darcyFoam) to obtain an initial single-phase velocity field which can be used for initialization and/or passive scalar transport.

Multi-species transport

The transport properties of species (dispersion coefficients, retardation, decay and other parameters) have been developed inside a new class multiScalarMixture derived from the OpenFOAM class basicMixture. This makes easy to manage the transport of a N-component mixture with homogeneous or heterogeneous properties.

Validation of new developments

To validate the development presented in the previous sections, each solver is tested independently in a configuration adapted to the studied mathematical model. The validation of the code is performed by comparing results with those obtained using METIS [START_REF] Goblet | Simulation d'écoulement et de transport miscible en milieu poreux et fracturé[END_REF], a finite-element solver for saturated/unsaturated flow that can handle multiple features such as solute transport or reactions. This code, developed during the last 30 years and still in development, has been used in several studies [START_REF] Renard | A new quasi-3D unsaturated-saturated hydrogeologic model of the Plateau de Saclay (France)[END_REF][START_REF] Robineau | A double medium approach to simulate groundwater level variations in a fissured karst aquifer[END_REF][START_REF] Mattei | Can soil hydraulic parameters be estimated from the stable isotope composition of pore water from a single soil profile?[END_REF] and is considered here as our reference solution. Numerical schemes for the finite-element validation code METIS are similar to those used by the PMF toolbox (i.e., second order in time and space).

The validation cases presented in the following are inspired by real cases and the physical values (time, flux, concentration) are representative of realistic configurations. The differences in orders of magnitude between output variables can thus be different and to improve the readability of this section, the choice is made (i) to present the real quantities in the configuration description of each test and (ii) to plot numerical comparisons with dimensionless variables. Quantities are systematically made dimensionless using the largest value from the reference solution, i.e. the METIS simulation.

Unsaturated flow solvers

The development of new solvers in OpenFOAM is independent of dimensionality of the domain (1D/2D/3D) thanks to the automatic discretization of the equations. Nevertheless, this is not the case for the reference code (METIS) and therefore validation is performed on a one-dimensional case even though the solvers in PMF can be used in two or three dimensions.

The physical configuration for the validation of unsaturated solvers consists of a 1D column 34.50 m in height discretized with 5009 computational cells (each with height ∆z = 0.0249 m). We consider the bottom of the column to be in contact with the saturated area by imposing a Dirichlet boundary condition (H bottom = 0 m) while an infiltration velocity (steady or transient) is imposed on the top. The model and physical properties are summarized in Table 1.

Dispersivity coefficient α L 1.0 [m] α T 0.2 [m] Volume partitioning coefficient K d 5.10 -5 [m 3 .kg -1 ] Radioactive decay coefficient λ 1.10 -9 [s -1 ] (b) Transport properties
Table 1: Properties for unsaturated validation cases

Steady flow (steadyGroundwaterFoam)

The first case consists in solving a steady flow with a constant infiltration velocity on the top of the column (V top = 4.753 • 10 -9 m/s). The convergence residual for the pressure head equation is set to 10 -10 and the numerical results are presented in Figure 3 (H is scaled using the height of the column). Although the agreement is excellent between the two codes, a slight difference can be seen close to the bottom where the height tends to 0. Due to the relatively large α value in the Van Genuchten model, the saturation variation occurs over a small height ( H H column < 0.01) with the presence of a strong pressure head gradient on few cells. Even if the size of the elements is similar between the two codes, METIS solves mass conservation at the nodes while OpenFOAM (with its finite volume method) solves conservation at the cell centers, so that boundary conditions are imposed slightly differently. This fact explains the small differences in the gradient computation near the boundaries.

Transient flow (groundwaterFoam)

Starting from the steady solution of the previous case, the infiltration velocity at the top takes the form of a piecewise linear function, with values varying between 0 and 1000 times the reference injection value (V top = 4.753 • 10 -9 m/s) over 2.10 8 seconds (∼ 2315 days). This variation of infiltration was chosen to verify robustness and stability in the case of a high amplitude and sudden injection in an unsaturated zone. The results over time are compared in terms of water flow rate out of the column (cf Fig 4) and show good agreement. Some small temporary differences are observed during the drainage phase (second half of the simulation), when flow dynamic is driven by the local pressure head gradients which may slightly differ between the discretization methods (as for the steady configuration).

Passive transport (porousScalarTransportFoam)

Using the pre-computed saturation and velocity field of the steady case, a passive tracer is injected from the top (Q = 1.374.10 -8 kg.m -2 .s -1 ) for a period of 6622 seconds. The physical time of the simulation is 1.262.10 9 seconds (∼ 40 years) and solute flux going out the column for both codes are presented in Figure 5. The good agreement observed allows to validate the good implementation of multiple elements influencing the solute transport, namely: the heterogeneous saturation, the dispersion, the decay and the retardation factor.

Coupled water/tracer simulation (groundwaterTransportFoam)

Starting from the steady configuration, water is injected inside the column (at z = 24 m) at fixed flow rate Q in = 3.82 • 10 -5 m 3 /s (i.e. 0.135 m 3 /h) with an imposed solute concentration C in = 6.12 • 10 -11 kg.m -3 from t start = 0 to t stop = 91 days; the injection is then stopped and the simulation continues until t end = 730 days. This configuration is challenging because it implies successively: (i) an abrupt variation of local water content, (ii) the displacement of a water front close to the saturated point (h ≈ -0.0033 and θ θmax ≈ 0.997) followed by a quasi-steady case, (iii) an abrupt change related to the end of injection and (iv) a capillary drainage over a long time period. Water and solute fluxes going out the column are presented in Figure 6 This configuration highlights the efficiency of the time step management method based on the time-step truncation error. This case presents various time scales which can be seen by observing the simulation time-steps reported in Figure 7. At the early beginning (t < 1 h), small time steps (0.1 < ∆t < 20 s) are required to correctly capture the local and rapid increase of saturation at the injection point. Then, the front displacement inside the column implies sharp local saturation and pressure head changes from (t ≈ 1 to t ≈ 21 h) with ∆t ≈ 9 s. Finally, the quasi-steady flow is established when the front reaches the bottom, which allows time step management to greatly increase the time step size to reach the maximal value ∆t max = 90 hours (user-defined value corresponding to the solver outputs). The abrupt changes of saturation can be visualized at t = 2160 h (90 days) when the injection is stopped. However, the capillary drainage that follows this event is slow and allows the solver to quickly restore the large time step.

Moreover, this configuration highlights the convergence issues of the Picard's algorithm compared to Newton's algorithm. Due to the strong non-linearities of the Van Genuchten function, particularly close to the saturation point (the derivative of the capillary function tends to infinity when saturation tends to 1), Picard's algorithm does not converge easily and requires both many linear iterations and small time steps. This results in a prohibitive computational cost which makes it impossible to simulate this configuration in a reasonable time (actual execution times range from fewer than 4 minutes when using the newly implemented Newton's algorithm to several hours with Picard's algorithm on a standard CPU). The direct comparison between the two codes on a single-core simulation highlights computation times of the same order of magnitudes with an increased efficiency of PMF, more or less important depending on the case.

2D-saturated flow solver (2D watershed modeling)

In the same way as for the validation of unsaturated solvers, the 2D-saturated solvers are tested on a realistic watershed. The same unstructured mesh presented in Figure 8a is used for both codes. However, it must be noted that due to the different numerical methods used, METIS (Finite Element) computes the solution on the mesh nodes while OpenFOAM (Finite Volume) computes the solution on mesh cells. Figure 8 shows the permeability field (figure 8b) and Digital Elevation Models (DEM) for the upper (figure 8c) and the lower (figure 8d) limits of the aquifer. As for the unsaturated simulations, METIS is considered our reference solution for the validation of PMF and comparisons are carried out on probes placed in the domain and visible in Figure 8d. All probes are used for the validation of the steady case but only the marked probes (probe A, shown in green , and B, marked in orange) are used for transient comparison.

Steady flow (steadyGroundwater2DFoam)

The first case consists of solving a steady flow with a constant and homogeneous infiltration velocity on the whole watershed (V inf iltration = 5.59 • 10 -9 m/s). The convergence residual for the potential equation is set to 10 -10 . Figure 9a shows the potentials and velocity fields in the steady configuration with the watershed outlet located in the lower right corner. Figure 9b shows the good agreement in terms of potentials between METIS and PMF at each probe.

Transient flow (groundwater2DFoam)

Starting from the steady solution of the previous case, the infiltration velocity on the top varies in time linearly and disorderly between 0 and 4 times the reference injection value (V inf,ref = 5.59 • 10 -9 m/s, the amplitude of variation corresponds to realistic values at the level of the saturated zone) over 4.2.10 8 seconds (∼ 4861 days). Temporal variation of dimensionless infiltration value ( 10 as well as the potential variations of probes A and B. Comparisons show good agreement with a slight difference between codes that occur during the event. This difference may be due to several factors related to numerical methods used such as the seepage condition which occurs in cells for PMF and in nodes for METIS or the post-processing operations (PMF results in mesh cells are interpolated to get probes value placed on mesh nodes).

V inf V inf,ref ) are visible in Figure

Passive transport (porousScalarTransport2DFoam)

Using the pre-computed saturation and velocity field, a passive tracer is injected (Q = 5.10 -6 kg.m -2 .s -1 ) for a period of 4.10 8 seconds while the physical final time of the simulation is t end = 4.73.10 9 seconds (∼ 40 years). Concentrations fields at t = 0.2 × t end and t = t end are presented in Figure 11. To validate both transversal and longitudinal dispersion, probes A and B have been placed respectively in the direction of and perpendicular to the main flow. The good agreement between the observed probes (c.f. Figure 12) allows to validate the correct implementation of the multiple elements influencing the solute transport, namely: the heterogeneous potential and water height fields, the dispersion related to the velocity field, the decay and the retention coefficient.

Coupled simulation (groundwaterTransport2DFoam)

Starting from the steady 2D simulation, a passive tracer is injected in a similar manner to the previous case (Q = 5.10 -6 kg.m -2 .s -1 ) for a period of 4.10 8 seconds with the final physical time of the simulation equal to t end = 4.73.10 9 seconds (∼ 40 years). A transient and variable infiltration is applied during the simulation inducing slight concentration variations related to water height changes as can be seen in Figure 13. As for previous 2D validation cases, the specificities of numerical codes induce slight differences with smaller transient variations in the case of METIS but with a globally similar dynamic behavior.

Parallel Efficiency

To ensure that the good parallel efficiency of PMF is preserved within the new developments, we performed a weak scaling test consisting in solving the same problem with the discretized domain splitted among an increasing number of processors. Test is performed on the two-dimensional steady case (cf 5.2.1) using a much finer mesh with a mesh resolution of ∆x = ∆y = 3 m and a total of 2.05 M cells. Simulations are performed using between 18 and 576 cores on the CALMIP Olympe cluster (Intel® Skylake 6140). Results are reported in Table 2 with CPU time t CP U computed as:

t CP U = t SIM × n cores , (21)
where t SIM is the total physical time of the simulation and n cores the number of cores used for computation. The CPU time per linear iteration t CP U,iter is obtained as:

t CP U,iter = t CP U n iter (22)
with n iter the number of linear iterations which may vary according to the domain decomposition. The results are plotted in Figure 14 and compared to linear behavior expected from the reference, i.e. one computation node composed of 36 cores (the 18-cores simulation occurs on a half node and can be shared with other tasks). The results highlight a super-linear behavior until 144 cores, as was previously observed with the OpenFOAM platform [START_REF] Horgue | An open-source toolbox for multiphase flow in porous media[END_REF][START_REF] Horgue | An extension of the open-source porousMultiphaseFoam toolbox dedicated to groundwater flows solving the Richards' equation[END_REF]. The parallel efficiency measured as CPU time per linear iteration remains lower than the reference (t CP U,iter = 86 ms) up to 288 cores before showing a deterioration. This deterioration has been observed previously with an efficiency decrease occurring when the number of cells per processor becomes lower than 5000 (around 3500 in this case). Poor results observed for 576 cores (with a large number of linear iterations inducing large physical time), are mainly due to the increase of the under-relaxation (0.1 vs. 0.5 for other configurations). This modification is mandated by the decomposition in small subdomains that makes the solving process unstable in the standard configuration.

Conclusion

The porousMultiphaseFoam toolbox, initially designed as a generic tool for the efficient parallel solution of multiphase flow through porous media, had shortcomings in hydrological modeling, with some missing preprocessing and solver features. The tool has now been improved and offers a very complete hydrogeological modeling suite with advanced numerical methods (such as automatic time step management and nonlinear algorithms), fast modeling methods (including simplified 2D approaches), passive/coupled transport modeling, and preprocessing tools essential to the setup of realistic cases. Every part of the new developments has been validated by a dedicated test case and a comparison with an existing reference code.

Main perspective of this work concerns the 3D modeling of water and solute transfers of a real configuration. The cases studied in this work are 1D or 2D to reduce computational cost and allow direct comparison with an well 1D/2D validated code. However, the developments made in the OpenFOAM environment can be directly used for 3D simulations and its good parallel efficiency make possible to compare simplified approach (1D+2D) with full modeling (3D) on real cases. Moreover, some improvements are planned in term of features and ergonomics. The development of a simple reactions module is in the implementation stage, as well as a coupling with a chemical modeling code as has been done by analogous packages based on OpenFOAM. In terms of ergonomics, an interface between the simulation code and a GIS management tool will be set up to allow use of the toolbox by a larger community.

Figure 1

 1 Figure 1: File structure

Figure 2 :

 2 Figure 2: sourceEventFile example with N dates and M injections points

Figure 3 :

 3 Figure 3: Saturation distribution with PMF and METIS in the steady case.

Figure 4 :

 4 Figure 4: Comparison of water outflow with PMF and METIS, transient flow case

Figure 5 :

 5 Figure 5: Tracer flux when solving for passive transport with PMF and METIS

Figure 6 :

 6 Figure 6: Fluxes on the bottom of the column as a function of time: (a) water flux, (b) solute flux

Figure 7 :

 7 Figure 7: Evolution of (automatically adjusted) time steps during the simulation

Figure 8 :

 8 Figure 8: Watershed for the 2D validation cases: (a) mesh (b) permeability field (c) DEM for upper limit (d) DEM for lower limit. Probes used for METIS/PMF comparison are visible on (a)

Figure 9 :

 9 Figure 9: Steady flow: (a) visualization of velocity and potential fields, (b) comparison of potential values between METIS and PMF at each probe.

  ncores

Figure 10 :

 10 Figure 10: Comparison of transient potential variations for probes A and B

Figure 11 :Figure 12 :Figure 13 :Figure 14 :

 11121314 Figure 11: Concentration visualization at (a) t = 0.2 × t end (b) t = t end

Table 2 :

 2 Simulation costs for parallel efficiency evaluation

	279.9	5037	111

n cells/core n iter t SIM [s] t CP U [s] t CP U,iter [ms] 18 113.6 × 10 3 45.2 × 10 3

Funding information

This work was performed using HPC resources from CALMIP (Grant P21009).