Capucine Legentil

Jeanne Pellerin

Paul Cupillard

Algiane Froehly

Guillaume Caumon

Testing scenarios on geological models: local interface insertion in a 2D mesh and its impact on seismic wave simulation

Keywords: Triangle meshing Gas-water contact Geomodeling Uncertainties Inverse problem Capucine Legentil: Conceptualization, Investigation, Methodology, Software, Writing (original draft) Conceptualization, Methodology, Supervision, Funding Acquisition, Resources, Validation, Writing Conceptualization, Supervision, Funding Acquisition, Validation, Conceptualization, Methodology

come

Testing scenarios on geological models: Local interface insertion in a 2D mesh and its impact on seismic wave simulation

Introduction

Geological models have been used for decades to make forecasts about the behavior of subsurface reservoirs (e.g., [START_REF] Ringrose | Reservoir Model Design[END_REF], seismic hazards [START_REF] Shaw | Unified Structural Representation of the southern California crust and upper mantle[END_REF] or underground waste disposal (e.g., [START_REF] Allen | Using simplified methods to explore the impact of parameter uncertainty on CO2 storage estimates with application to the Norwegian Continental Shelf[END_REF][START_REF] Mari | The application of high-resolution 3D seismic data to model the distribution of mechanical and hydrogeological properties of a potential host rock for the deep storage of radioactive waste in France[END_REF][START_REF] Mulrooney | Structural analysis of the Smeaheia fault block, a potential CO2 storage site, northern Horda Platform, North Sea[END_REF]. These multi-material models can be seen as three dimensional geological maps identifying the geometry of interfaces such as horizons, faults and intrusions. In general, these models are created from the observations at hand (e.g., borehole and seismic reflection data), then meshed and filled with petrophysical properties before solving the relevant physical equations [START_REF] Ringrose | Reservoir Model Design[END_REF][START_REF] Wellmann | 3-d structural geological models: Concepts, methods, and uncertainties[END_REF]. However, ambiguities and lack of observations often raise uncertainties about the location or even the existence of some interfaces (see [START_REF] Wellmann | 3-d structural geological models: Concepts, methods, and uncertainties[END_REF] and references therein). Additionally, temporal evolution at geological time scales or during subsurface engineering projects may involve geometric changes in the geomodel features. Assessing the impact of the uncertainties and/or the temporal changes on the model forecast is a complex computational challenge, which can be partly addressed by deforming an existing mesh to avoid model reconstruction [START_REF] Abrahamsen | Bayesian kriging for seismic depth conversion of a multi-layer reservoir[END_REF][START_REF] Tertois | Editing faults within tetrahedral volume models in real time[END_REF]. However, remeshing is needed when interfaces are displaced along faults [START_REF] Caumon | Constrained modifications of non-manifold B-reps[END_REF], or when new surfaces are introduced in the model [START_REF] Suter | A Novel Method for Locally Updating an Earth Model While Geosteering[END_REF].

However, the repeated updating of a geological model is a di cult task. There are several strategies to edit an interface or add a new one in an unstructured mesh. The most reliable is to rebuild the model from scratch, compute the connectivity of the layers and fault blocks, and generate a new mesh of the whole model [START_REF] Bidmon | Intuitive, Interactive, and Robust Modification and Optimization of Finite Element Models[END_REF].

Another strategy is to directly modify the meshed model [START_REF] Tertois | Editing faults within tetrahedral volume models in real time[END_REF], but the output mesh validity must be guaranteed.

The insertion of a polygonal line in a 2D triangulated surface is an example of such a mesh updating operation.

Although it seems simple at first sight, it can be challenging to implement in a robust way because of limited floating point precision. For example, sequential insertion of several tangential lines representing close stratigraphic surfaces or fractures in a fracture corridor may then raise invalid mesh features such as inverted triangles. A good test of robustness is, therefore, to check that results are consistent and stable when inserting the same interface several times. Exact geometric computations (Lévy, 2016;[START_REF] Li | Recent progress in exact geometric computation[END_REF][START_REF] Shewchuk | Adaptive precision floating-point arithmetic and fast robust geometric predicates[END_REF] do address these issues in combination with constrained remeshing algorithms (e.g., [START_REF] Si | TetGen, a delaunay-based quality tetrahedral mesh generator[END_REF]. However, they involve algorithmic complexity and a computational overhead in terms of memory and time, which makes the problem « solved but practically challenging » [START_REF] Li | Recent progress in exact geometric computation[END_REF]. Therefore, we consider in this work a subclass of mesh updating problems where the element to insert is represented implicitly by a level-set of a piecewise linear function defined on the mesh. This considerably simplifies the management of geometric accuracy, as the target interface and the mesh to update are represented on the same data structure. This implicit interface representation has been used in geosciences for structural modeling for several years (e.g., [START_REF] Chilès | Modelling the geometry of geological units and its uncertainty in 3d from structural data: the potential-field method[END_REF][START_REF] Houlding | 3D Geoscience Modeling[END_REF][START_REF] Mallet | Three-dimensional graphic display of disconnected bodies[END_REF], and it can be applied both on two-dimensional and three-dimensional models.

In the mesh updating examples mentioned above, the modification is often done locally in areas of particular interest (e.g., subdomains where the uncertainty is large, where sensitivity of the model features is important, or areas where new data have been acquired and the model needs to be edited). Therefore, local mesh editing approaches, where the vast majority of the mesh is kept intact, are essential for computational e ciency. However, they raise specific challenges to ensure the mesh conformity between the intact and the edited parts [START_REF] Caumon | Constrained modifications of non-manifold B-reps[END_REF][START_REF] Suter | A Novel Method for Locally Updating an Earth Model While Geosteering[END_REF].

If the insertion of lines in a 2D triangle mesh may be implemented e ciently, the corresponding 3D problem, i.e., inserting a triangulated surface or discretized lines in a tetrahedral mesh is more di cult. The steps are the same than in 2D: for each element to insert, compute the intersection with the mesh, then insert the intersected elements in the mesh.

In this paper, we propose a first step to address the above challenges by taking advantage of recent advances of mesh updating libraries developed for shape optimization. Because implementing robust and e cient mesh generation and mesh modification algorithms is extremely challenging, and since open-source implementations are available for triangle and tetrahedral meshes [START_REF] Dapogny | Three-dimensional adaptive domain remeshing, implicit domain meshing, and applications to free and moving boundary problems[END_REF][START_REF] Geuzaine | Gmsh: A 3-D finite element mesh generator with built-in pre-and post-processing facilities[END_REF][START_REF] Si | TetGen, a delaunay-based quality tetrahedral mesh generator[END_REF], we focus on the integration of one of these libraries in an automatic uncertainty evaluation workflow. Relying on the level-set formulation of the interface to be inserted in the model, we consider the problem of monitoring the evolution of a gas-water contact in a subsurface reservoir using seismic data. We start from an input meshed model conforming to the main geological structures that is locally updated by inserting the fluid contact in the reservoir layer. For simplicity, we keep these structures frozen and we only assume that the contact depth is unknown. Because fluid content directly a ects the e ective elastic properties of the porous medium, fluid substitution can, in principle, be detected using seismograms as done for instance in time-lapse seismic surveys (e.g., [START_REF] Arts | Monitoring of CO2 Injected at Sleipner Using Time Lapse Seismic Data[END_REF][START_REF] Landrø | Discrimination between pressure and fluid saturation changes from marine multicomponent time-lapse seismic data[END_REF]. We quantitatively evaluate the impact of contact depth on the simulation results. Under a simple assumption of horizontal fluid contact, we show that our approach can be integrated into an inversion loop to determine the depth of the fluid contact and reduce uncertainties in the reservoir model from only a few seismic sensors.

The contributions of this paper are the following: a practical approach to automatically update triangular meshes for physical simulations is proposed. A special care is taken to the mesh validity, so that the model is suitable for physical simulations. An assessment of the detectability of fluid contacts changes in reservoirs is presented. The algorithms are available on Github and distributed under a GPL license, allowing reproducibility.

The geological models we consider are meshed with triangles that are conforming to geological interfaces such as horizons and faults (Figure 1). We propose to modify a single layer (or region of the model) using the remeshing platform Mmg (Dapogny et al., 2014) (Section 3.2.1). This library, initially developed by [START_REF] Dobrzynski | Anisotropic delaunay mesh adaptation for unsteady simulations[END_REF], has been extended with level-set functionalities and is used to solve shape optimization problems (e.g., [START_REF] Allaire | Structural optimization using sensitivity analysis and a level-set method[END_REF][START_REF] Osher | Level Set Methods for Optimization Problems Involving Geometry and Constraints I. Frequencies of a Two-Density Inhomogeneous Drum[END_REF][START_REF] Sethian | Structural Boundary Design via Level Set and Immersed Interface Methods[END_REF]. A few adaptations enable to use this library to modify geological models (Section 3.2) and a careful choice of options enables to obtain meshes of good quality, for numerical simulations (Section 4). The approach is demonstrated on a simple synthetic model. Seismic wave propagation simulations are computed using the Discontinuous Galerkin approach implemented in Hou10ni2D [Hu-dee-nee] [START_REF] Barucq | E cient DG-like formulation equipped with curved boundary edges for solving elasto-acoustic scattering problems[END_REF][START_REF] Barucq | Activity report magique 3d 2016[END_REF] (Section 4.1). The simulation results for di erent gas-water contact depth are compared to a reference one in order to solve a Bayesian inverse problem, leading to a probabilistic estimation of the contact depth (Section 4.2). (1) From several level-sets defining the possible interface positions between water-saturated sandstone and gas-saturated sandstone, we use Mmg [START_REF] Dapogny | Three-dimensional adaptive domain remeshing, implicit domain meshing, and applications to free and moving boundary problems[END_REF] to locally remesh a selected layer (dark blue) and (2) the Hou10ni2D software [START_REF] Barucq | Activity report magique 3d 2016[END_REF] to simulate wave propagation and (3) compare the simulation results.

State of the art

Definitions Physical modeling in the subsurface generally involves the construction of a geological model representing the interfaces between rocks that have di erent physical properties (such as porosity, permeability, density, Young modulus, shear modulus) [START_REF] Ringrose | Reservoir Model Design[END_REF][START_REF] Wellmann | 3-d structural geological models: Concepts, methods, and uncertainties[END_REF]. These properties depend on multiple factors (e.g., facies, fluid saturation, compaction, temperature). Geological models, therefore, are multimaterial models, a term used more generally in mechanics. The interfaces between di erent materials correspond to horizons, faults, unconformities, intrusions, cavities and fluid contacts. Among the various possibilities to represent geological media, unstructured meshes o er the interesting ability to honor complex geometries and to spatially vary the level of detail. In previous works, 3D models have been considered, but for simplicity here we consider triangulated two-dimensional geological models similar to those defined in [START_REF] Anquez | Automatic correction and simplification of geological maps and cross-sections for numerical simulations[END_REF] or [START_REF] Pellerin | RINGMesh: A programming library for developing mesh-based geomodeling applications[END_REF] (Figure 2). Indeed 2D models are most e ective for proof of concept, as they allow to assess the impact of changes on physical simulation in acceptable time (Section 4.1) and they are commonly used for physical models in cylindrical domains.

Generalization to 3D will be discussed in Section 5. Triangulated models are adequate to run wave simulations with a Discontinuous Galerkin method provided that the mesh is valid and of su cient quality. By valid we mean that all intersections between two triangles of the mesh are either empty or a common edge/vertex to both triangles and that no element is empty. The quality of the mesh is related to the size and shape of the triangles. It has an impact on the accuracy and computation time, and its definition is specific to each solver (e.g., [START_REF] Loseille | Fully anisotropic goal-oriented mesh adaptation for 3D steady Euler equations[END_REF]. The interfaces between two materials are represented by discretized lines that separate two regions of the model. An alternative representation of interfaces in geological modeling, and more generally in multi-material modeling, is the level-set representation (also called implicit). In geological modeling this representation has gained popularity and is the base of the implicit stratigraphic modeling. In this context, the scalar field representing stratigraphy in a conformal layer can be seen as relative geological time [START_REF] Caumon | Three-dimensional implicit stratigraphic model building from remote sensing data on tetrahedral meshes: Theory and application to a regional model of la Popa Basin, NE Mexico[END_REF][START_REF] Collon | 3D geomodelling in structurally complex areas -Implicit vs. explicit representations[END_REF][START_REF] Chilès | Modelling the geometry of geological units and its uncertainty in 3d from structural data: the potential-field method[END_REF][START_REF] Houlding | 3D Geoscience Modeling[END_REF][START_REF] Ledez | Robustness and e ciency of geometric programs: The Predicate Construction Kit (PCK)[END_REF][START_REF] Mallet | Geomodeling[END_REF][START_REF] Tertois | Editing faults within tetrahedral volume models in real time[END_REF]. In some cases, level-set interfaces can be directly integrated in physical solvers using Heaviside basis functions as in XFEM [START_REF] Moës | Non-planar 3D crack growth by the extended finite element and level sets-Part I: Mechanical model[END_REF] or to define specific boundary conditions (e.g., [START_REF] Durand-Riard | Balanced restoration of geological volumes with relaxed meshing constraints[END_REF]. Level-set representation has also gained popularity to solve subsurface inverse problems (e.g., [START_REF] Cherpeau | Method for Stochastic Inverse Modeling of Fault Geometry and Connectivity Using Flow Data[END_REF][START_REF] Giraud | Generalization of level-set inversion to an arbitrary number of geological units in a regularized leastsquares framework[END_REF][START_REF] Guillén | Minimising the delta test for variable selection in regression problems[END_REF][START_REF] Zheglova | Multiple level-set joint inversion of traveltime and gravity data with application to ore delineation: A synthetic study[END_REF].

Editing models Di erent types of local modifications are described in the literature. Some authors propose to locally deform the mesh using geometric methods [START_REF] Laurent | Interactive editing of 3D geological structures and tectonic history sketching via a rigid element method[END_REF][START_REF] Tertois | Editing faults within tetrahedral volume models in real time[END_REF]. These constrained deformation methods are based on the definition of a smooth deformation field interpolated from a few anchor points.

One of the problems encountered when modifying the geometry of a set of vertices in a mesh is that the elements may become invalid, as some intersections may occur between the edges of the triangles or the faces of the tetrahedra [START_REF] Tertois | Editing faults within tetrahedral volume models in real time[END_REF].

However, general displacement, insertion and removal of boundaries lead to more substantial topological changes a ecting all model entities (nodes, edges, regions). Sword [START_REF] Sword | Building flexible interactive geologic models[END_REF] proposes methods to edit and update the object adjacencies and to verify the model topological validity at each modification (edge intersection verification, for example). In 3D, [START_REF] Euler | Editing and rapidly updating a 3d earth model[END_REF] propose a general way to edit boundary representations which relies on intersections and border projections of triangulated surfaces before globally reconstructing a boundary representation. Although this approach is very flexible in principle, the use of projection raises instabilities and the approach does not provide a strict guarantee that the model integrity is maintained throughout the editing. As part of the update of geomodels during drilling monitoring, [START_REF] Suter | A Novel Method for Locally Updating an Earth Model While Geosteering[END_REF] propose to modify the only geological block involved in the new data which does not correspond to the initial model. The region to modify is completely erased and meshed again taking into account the new data (Figure 3), but mesh conformity, required in most physical simulation codes, is not discussed. We follow a similar strategy in this work, but we also capitalize on recent advances in shape optimization that regroups a set of methods for finding the "best shape" of an object so that it fulfills its functions. This optimization is automated and the geometrical modeling also uses unstructured meshes. The management of heterogeneity as present in geological models, introduces specific di culties compared to the optimization of homogeneous objects [START_REF] Wang | A multi-material level set-based topology and shape optimization method[END_REF]. Level-set methods are very interesting because they allow tracking fronts and free boundaries. [START_REF] Sethian | Structural Boundary Design via Level Set and Immersed Interface Methods[END_REF] first introduced the level-set description into topology optimization. Level-set techniques were then developed in, e.g., [START_REF] Allaire | Structural optimization using sensitivity analysis and a level-set method[END_REF], [START_REF] Osher | Level Set Methods for Optimization Problems Involving Geometry and Constraints I. Frequencies of a Two-Density Inhomogeneous Drum[END_REF], [START_REF] Sethian | Structural Boundary Design via Level Set and Immersed Interface Methods[END_REF], [START_REF] Wang | A level set method for structural topology optimization[END_REF] to perform topological changes of the level-set components. An interface ⌦ is defined by an iso-value ' 0 of a level-set function (i.e, a scalar field) ' as

⌦ = x, '(x) = ' 0 .
(

) 1
This level-set representation is able to represent geometric interfaces between materials, such as geological interfaces, as shown in Figure 4. [START_REF] Allaire | Shape optimization with a level set based mesh evolution method[END_REF] propose to explicitly discretize the level-set at each iteration in a computational domain equipped with a simplicial mesh. In doing so, a sub-mesh of the domain is created at each step. mesh (composed of the white and grey elements). The iso-value 0 represents the interface between the white and grey regions.

In the shape optimization by a physical constraint, the accuracy of the border has a direct impact on the precision of the simulations. To build the obstacle boundaries, [START_REF] Santosa | A level-set approach for inverse problems involving obstacles[END_REF] proposes an inverse approach involving a set of level-set functions that minimizes a residual to fit data. To deal with multi-material models, [START_REF] Wang | Design of Multimaterial Compliant Mechanisms Using Level-Set Methods[END_REF] and [START_REF] Wang | A multi-material level set-based topology and shape optimization method[END_REF] introduce a "color" level-set method, where each material is defined by a level-set function. This approach is similar to a subsurface inverse problem considering, for example, the position of a horizon as a parameter set. To solve an inverse problem in full waveform case, Guo and de Hoop (2013) adapt the method of [START_REF] Santosa | A level-set approach for inverse problems involving obstacles[END_REF] and use concepts from shape optimization and image segmentation to determine the position of a salt body.

Implicit local interface insertion

In this section we propose to insert a level-set interface in a mesh by restricting the modifications to a prescribed region of the model. The chosen level-set iso-value is first discretized, and the modified region is remeshed conformably to prevent a steep decrease of the mesh quality that would slow down further computations or make them unstable.

When merging the edited region with the other parts of the model, specific care to preserve the mesh conformity is taken.

Input Data

Level-set function To perform the local insertion of an interface in the geomodels, the input data is composed of the global mesh and the interface implicitly defined (e.g, a fault, a gas-water contact...). In our case, the level-set is approximated by the iso-values of a regular piecewise linear scalar field represented on the nodes of the triangular mesh. Typically, the scalar field represents the signed distance to the line, which is defined by the 0-valued segment in each intersected triangle (Figure 5b, Figure 6b). The insertion of the new points and edges is based on this implicit line definition.

Model partitioning

In our method, the area to be modified (inner region) consists of one or several input surfaces representing geological units, but this can easily be changed to a list of triangles, e.g., to insert an internal fracture.

The implicit line will be inserted in the inner region while maintaining the conformity to the lines and corners of that region, e.g., the blue layers (# A1, A4, B2 and B3) in Figure 6a. The rest of the model (outer region), is not a ected by the modification, except a small set of triangles which needs to be modified to maintain the model conformity. In Figure 6a, all the triangles of the outer region that are in contact with the bold boundaries could be a ected by this conformity requirement; in practice, as will be further explained in Section 3.2.3, only the triangles that are in contact with the level-set need to be considered.

Local level-set discretization

Mmg software

We use the Mmg code for the level-set insertion in the inner region. Mmg is an open source 2D and 3D remeshing software with two major features: level-set discretization and mesh adaptation. Mmg is dedicated to the processing of triangle and tetrahedral meshes. It is a meshing software mainly used for computational fluid dynamics applications and for shape optimization applications [START_REF] Benard | Mesh adaptation for large-eddy simulations in complex geometries[END_REF][START_REF] Loseille | Fully anisotropic goal-oriented mesh adaptation for 3D steady Euler equations[END_REF]. It has also recently been e ectively integrated to adapt meshes in full waveform inversion problem to accelerate the computations [START_REF] Jacquet | FWI (Full Waveform Inversion) dans le domaine temporel utilisant des méthodes numériques hybrides pour la caractérisation de milieux élasto-acoustiques[END_REF].

The discretization of level-sets allows to insert a line or a surface that represents a scalar field iso-value in a mesh.

The scalar field is defined at all the nodes of the mesh, forming a piecewise linear function. The software takes into account multi-material models, so the discretization of the iso-value can be limited to determined materials and maintain existing materials.

The other main feature of Mmg is the mesh adaptation: the user can set several parameters to control the size of the mesh elements. For example the hmin and hmax parameters are used to define the minimum and maximum size of the elements; the parameter hgrad controls the ratio between the length of two adjacent edges. The user can set several command line parameters or provide an input size map to control the size of mesh elements.

Level-set discretization

The mesh modifications are made with the level-set discretization method. The algorithm implemented in Mmg [START_REF] Dapogny | Three-dimensional adaptive domain remeshing, implicit domain meshing, and applications to free and moving boundary problems[END_REF] consists of two main steps. First, the line is inserted in the defined region as simply as possible while keeping the mesh consistent: each triangle intersected by the level-set is cut into either two triangles (when the level-set goes exactly through one vertex) or one triangle on one side and two triangles on the other side (when the level-set cuts the interior of two edges) (Figure 7). However, most of the newly created triangles are of poor quality (Figure 8.a.), so remeshing is essential before running a seismic simulation. In the interface discretization, the major constraint is to guarantee the conformity of the entire mesh for both the inserted and pre-existing boundaries. To this end, only the triangles adjacent to the modified region and in contact with the ends of the implicit interface are modified. They are simply divided into two new triangles (Figure 9).

The inner region is then divided into two regions corresponding to the two sides of the inserted interface. The two parts of the mesh are remeshed separately to adapt the mesh to the constraints imposed by the user: the maximum or minimum element size and the gradation (i.e., the ratio of adjacent elements size). The operations performed on the mesh consist of insertion or displacement of nodes, edge swap or collapse. All decisions are made to adapt to the size map defined as input and to improve the mesh quality. To maintain the conformity of the overall mesh, the edges between sub-parts must not be modified (Figure 8b).

Merge

The last step consists in merging the modified and remeshed region and the preserved one. During this step some points are duplicated. The mesh cannot have points at the same position for our wave propagation simulation. Duplicated nodes are removed by using co-location detection from Gmsh 1 . This API function is based on a R-Tree [START_REF] Guttman | R-trees: a dynamic index structure for spatial searching[END_REF]).

1 https://gitlab.onelab.info/gmsh/gmsh

Mesh quality

To evaluate the mesh quality we propose to compute the minimal height of the worst element among the n t triangles:

Q = min ≈kÀ[[1;n t ]] h k , (2)
where h k is the lowest height of the triangle k. The quality Q has an important impact on elastic wave simulation time because of small time steps imposed by the Courant-Friedrichs-Lewy (CFL) condition (Section 4.2). To minimize the decrease of quality, if the distance between the level-set and a node located on another interface is less than a given tolerance ✏, the scalar field is locally modified: the iso-value is set at this existing node so that no new node is created during the line insertion in Mmg (Figure 10). In Figure 10, the di erence of the scalar field values of two vertices in a triangle ' 1 * ' 0  is normalized by the norm of gradient of the level-set function in the triangle  í ('. This value is then compared to the tolerance ✏ that can be assimilated to a distance. By adjusting the parameter ✏, one can limit the creation of very elongated triangles which otherwise appear at the intersection between the iso-value and the borders of the layer, at the expense of a small inaccuracy in the final geometry.

Figure 10: During the level-set discretization, if the distance between the iso-surface and an existing node is less than ✏ (2), the iso-value is approximated at this point and no other point is created when this approximation is not made, ill-shaped triangles are formed (in red).

Open-Source Implementation

Process automation We propose algorithms for the automation of the process. The code is available on Github2 and depends on Mmg3 [START_REF] Dapogny | Three-dimensional adaptive domain remeshing, implicit domain meshing, and applications to free and moving boundary problems[END_REF] and Gmsh4 [START_REF] Geuzaine | Gmsh: A 3-D finite element mesh generator with built-in pre-and post-processing facilities[END_REF]. It allows to create the two sub-meshes to restrict the discretization to the modified region, and then to merge them. Duplicated nodes are removed by co-location detection from Gmsh. It is then possible to create models that correspond to the input format of the simulator. To prevent remeshing of the entire mesh, the options -nomove, -noinsert, -noswap are used on the command line when calling Mmg for level-set discretization. The -nomove option blocks the movement of vertices and the -noinsert one prevents the creation of new vertices. The -noswap option prevents the edge swapping between two adjacent triangles.

Element quality improvement A large part of finite element schemes solving boundary values problems imposes the absence of triangles whose edges are on two di erent borders of the model (Figure 11). In our application, the edges are not domain boundaries but only material boundaries, hence they do not carry boundary conditions. Therefore, we activate the -nofem option in Mmg to improve the mesh quality. The options hmin and hmax can also be used to set the mesh size to adapt the mesh to the physical simulator. triangles with two constrained edges are allowed, yielding a higher quality mesh for our application.

However, there are very thin objects in geology, which inevitably yield poor quality triangles. During the physical simulations, such triangles show a very small value of the determinant of their Jacobian matrix, leading to potential inaccuracy in the simulation result. In dynamic problems, these triangles also imply a small time-step because of the CFL condition, which impacts the computational cost. Additionally, it could be combined with approaches to simplify models in order to reduce the simulation time. For example, during the level-set discretization, an option -rmc is available to remove small surface parts. The decision is taken by comparing the surface of the small component to the region surface.

For each new model created, the mesh quality (2) is computed. If the quality is too low (the choice of the quality threshold is discussed in Section 5), the process can be repeated after approximating the scalar field at several nodes by increasing ✏ to avoid problematic configurations (Figure 10, Figure 14). A new model is then created with less elongated triangles.

File format Mmg uses the Medit file format (.mesh). We use the tag that may be set on the elements to encode model entities. The scalar field and/or metric for remeshing values are defined on the vertices of the mesh and are stored in an other file format (.sol). The local parameters for multi-material models are defined in a 2 dimensional Mmg file (.mmg2d).

Gas-water contact uncertainty and elastic wave propagation

In this section, we apply our local level-set based model edition to evaluate the impact of a gas-water contact depth on seismic signals. We consider a synthetic model on which we quantitatively seek to evaluate the correlation between the contact depth and some seismic waveform misfits, demonstrating the usefulness of our approach to e ciently sample and reduce uncertainties on a material interface position.

Our synthetic geological model M ref is an anticline overlying horsts and grabens (Figure 12). The anticline is composed of folded layers of shale and sandstone. After the insertion of the gas-water contact in the reservoir, the gas zone is represented in grey (layer 6 in Figure 12). The domain size is 16.2 km horizontally and 5 km vertically.

The petrophysical values of the di erent layers of the model are summarized in Table 1. These values were chosen consistently with the case study of [START_REF] Hamada | Reservoir fluids identification using vp/vs ratio?[END_REF]. For simplicity, we use a constant density d = 2600 kg_m 3 . Lithological and petrophysical characterization of the model Figure 12. Density is constant d = 2600 kg_m 3 .

Figure 13: The initial mesh of our synthetic model generated with Gmsh and Mmg, contains 22,077 triangles.

Seismic simulation parameterization

In this study, we use the Hou10ni2D software [START_REF] Barucq | Activity report magique 3d 2016[END_REF] which implements a Discontinuous Galerkin (DG) method [START_REF] Reed | Triangular mesh methods for the neutron transport equation[END_REF] for solving the wave equation. More precisely, Hou10ni2D relies on Interior Penalty approach [START_REF] Grote | Discontinuous Galerkin finite element method for the wave equation[END_REF][START_REF] Riviere | Discontinuous Galerkin methods for solving elliptic and parabolic equations: theory and implementation[END_REF] to simulate time-dependent elastic wave propagation. Hou10ni simulates seismic wave propagation in heterogeneous media and allows the use of di erent polynomial orders in the mesh elements (p-adaptivity in space) [START_REF] Barucq | E cient DG-like formulation equipped with curved boundary edges for solving elasto-acoustic scattering problems[END_REF]. In addition, Hou10ni2D can run on parallel architectures.

We apply absorbing boundary conditions on the two lateral sides and the bottom of the cross-section, whereas the top of the domain is defined as a free surface. We choose the source time-function as a Ricker wavelet with a dominant frequency f = 15 Hz. As a consequence, the minimum wavelength is Model characteristics (Contact depth, number of triangles, worst element quality, simulation time and number of connected reservoir parts). The quality is the minimal triangle height (2). Simulations are run on a machine with a processor Intel(R) Xeon(R) CPU E5-1620 v2 @ 3.70GHz, 8 Cores.

min Ù v s 3f , (3) Depth (m
where v s is the S-wave velocity. The source is located at depth z = *775 m, two kilometers from the left side. We place a line of 150 receivers (one every 100 m) at the top of the domain (dashed green line in Figure 12). Simulations are run for T = 4 s after the shock. The displacement is computed for all the receivers at each time step.

Comparison of simulation results with a varying gas-water contact depth

Set of models We built 20 models by inserting the gas-water contact at several depths in the initial mesh Figure 13, using our method presented in Section 3. The process takes about 10 seconds per model. Depths vary from z min = *2, 845 m to z max = *2, 645 m (between the two dashed red lines in Figure 12), with a constant step z equal to 10 m.

Each model is characterized (Table 2) by the contact depth, number of triangles, mesh quality, simulation time, and number of connected reservoir parts (Nb CC in Table 2). To evaluate the quality, we chose the minimal of all triangle heights (2), as it is directly linked to the software heuristics for the determination of the time step and the level of the p-refinement. This value has a strong impact on the simulation run time of simulation as can be seen in Figure 14.

Impact of quality on simulation time

The mesh quality depends on the model geometrical configuration. There is no direct link between the worst element quality, the depth of the contact or the number of connected components of the gas-saturated layer (#6 in Figure 12). Poor quality triangles are located around the ends of the gas-water contact line.

They often are small and elongated triangles. The existence of these ill-shaped triangles is partly a consequence of low angles between the contact and stratigraphic boundaries (close-up views in Figure 14b and Figure 14c). However, the position of the new corners is the parameter that most impacts the mesh quality and the simulation time. If the corner is very close to a vertex of an adjacent triangle of the outer region (#5 in Figure 12), the quality of the triangles created is low (Figure 14a andb). The approximation of the scalar field described in Section 3.2 and Figure 10 avoids the creation of very poor quality triangles. Table 3 Impact of elongated triangles deletion: the simulation time is between 12 and 95 times shorter. Simulations are run on a machine with a processor Intel(R) Xeon(R) CPU E5-1620 v2 @ 3.70GHz, 8 Cores.

The scalar field has been locally approximated for two models: z = *2, 775 m and z = *2, 705 m with ✏ = 0.5 m.

These small modifications lead to a large quality improvement and a significant reduction of the computation time (Table 3).

Comparison with the reference model For each modified model M z , a simulation is run using the corresponding velocity model. Then the relative di erence d i,z at each receiver i between the seismogram u i,ref of the reference model M ref (Figure 12) and the seismogram u i,z simulated in each model M z is calculated [START_REF] Geller | A new method for computing highly accurate DSM synthetic seismograms[END_REF]:

d i,z = r r r p î T 0 Ù Ù Ù u ref (t) * u i,z (t) Ù Ù Ù 2 dt î T 0 Ù Ù Ù u ref (t) Ù Ù Ù 2 dt a s s s q 1 2 (4)
The results are di cult to interpret because the simulated wavefield is complex. From Figure 15, we see that the largest waveform di erence ( d i,z > 20 %) is obtained between x = 5, 200 m and x = 13, 400 m. On average, the seismogram error increases when the contact depth error z increases. Horizontally, however, Figure 15 shows that some receivers are not sensitive to contact depth uncertainty, e.g., from x = 0 m to x = 2000 m and from x = 13, 600 m to x = 15, 000 m. This corresponds to receivers for which the possible contacts are not illuminated by the wavefield.

Interestingly, we also observe that the local error decreases with increasing depth di erences, e.g. from x = 6, 500 m and x = 7, 500 m with z > 110 m. This highlights the importance of the starting model to solve inverse problems from sparse seismic measurements to avoid being trapped in local minimum when only a few receivers are used. Due to the complexity of these results, a comparative analysis on a single receiver is delicate, and it is more relevant to take into account the seismograms at a large set of receivers to analyze the results.

To more easily compare simulations we choose a unique indicator, namely the mean z of the errors on the n receivers:

z = n ≥ i=0 d i,z n (5)
We observe in Figure 16 a correlation between the mean error and the depth error z : the larger the distance from the reference depth contact, the higher the mean of L2-norm.

The inverse problem we consider consists in finding a model M z containing a contact at the depth z, such that its seismic response s(z) is "close" to a reference seismic response, obtained from the model M ref defined as s(z) = s(z ref ). The di erence in seismic response is evaluated using the L2-norm (4 and 5). The inversion loop to determine the best model that fits the seismic response of the reference model is, in general, a di cult process. To solve a Bayesian inverse problem, many di erent models have to be generated and seismic simulations are required for each of them.

Here, the inversion only concerns one parameter (the contact depth), so we use a simple grid search (brute-force) approach using the previous results (Figure 15 and Figure 16).

Quantification of uncertainties

To test the impact of the modifications on the simulations, and to study the feasibility of integrating our approach in an inversion loop, we compute the density of probability describing the depth of the contact. Bayes theorem (e.g., [START_REF] Mackay | Information theory, inference and learning algorithms[END_REF] is used to combine prior information on the structural model with the observed data to give the posterior probability density function

⇢(zD) = ⇢(z)L(D, z) î z ⇢(z)L(D, z)dz (6)
where ⇢(z) is the prior probability density, i.e., the information we have on the contact depth before considering the data D. In this study, we choose a non-informative (uniform) prior probability density

⇢(z) = h n l n j 1 z max *z min if z À ⌅ z min ; z max ⇧ 0 otherwise . ( 7
L(D, z) = ⇢(Dz) = e * 1 2 Ò D(z)*D ref Ò 2 2 ˘2⇡ , with = ÒDÒ (8)
This makes it possible to use the values of L2-norm introduced in the previous section to calculate the posterior density of probability. By using (4) to simplify (8) we obtain

Ù Ù Ù D(z) * D ref Ù Ù Ù 2 2 = … i ( d i,z ) 2 = d 2 z , (9)
(10) and then

L(D, z) = ⇢(Dz) = e * 1 2 d 2 z ˘2⇡ (11)
The posterior probability density function is computed for the set of models and we obtain the results presented in Figure 17. We can see that the density of probability is maximum around the reference depth. More generally, such a density of probability quantifies the posterior uncertainties on the depth of the contact.

Discussion and conclusion

In this paper, we have proposed a method to locally modify a meshed geomodel by inserting a line implicitly defined. This avoids to globally rebuild the structural model and create a new mesh, from scratch, which takes time and e ort (e.g., [START_REF] Pellerin | Elements for measuring the complexity of 3d structural models: Connectivity and geometry[END_REF][START_REF] Zehner | Workflows for generating tetrahedral meshes for finite element simulations on complex geological structures[END_REF]. Local editing techniques on geological models allow for small adjustments to the models without having to go through the entire mesh construction process. Although the principles are similar to those highlighted by [START_REF] Suter | A Novel Method for Locally Updating an Earth Model While Geosteering[END_REF], we have proposed a strategy to maintain the mesh conformity, which is a prerequisite for many physical modeling codes. The implicit formulation is interesting for mesh simplifications and generalizations to more complex interfaces as long as they can be represented on the input mesh by a piecewise linear scalar field. This formulation allowed us to demonstrate the practical impact of model updating on wave propagation in a simple example of gas-water contact insertion into a cross-section. The simulation results permit to solve a Bayesian inverse problem, leading to a probabilistic estimation of the contact depth.

Some of the modified models have poor quality triangles, often located near the newly created corners. These triangles strongly impact the simulation time. In order to automate the reversal process, a quality threshold could be defined to automatically adjust the tolerance ✏ used in Section 3.3 and avoid excessively long simulations (Figure 18).

Although this approach locally approximates the target geometry by projecting the level-set on nodes of the input mesh, it allows to automatically and robustly reduce mesh quality issues. As compared to the simplification methodology proposed by [START_REF] Anquez | Automatic correction and simplification of geological maps and cross-sections for numerical simulations[END_REF], this approach is more specific as it considers one interface at a time, but it is significantly easier to implement. So whereas the initial target of this work was not global model simplification as in the case of [START_REF] Anquez | Automatic correction and simplification of geological maps and cross-sections for numerical simulations[END_REF], we believe that a carefully designed incremental insertion of implicit interfaces could provide an alternative way to robustly address this problem by locally merging interfaces which are too close one from another. One may argue that such approximations should be avoided, but we consider that some of them are inevitable to obtain physical solutions in acceptable time. As discussed by [START_REF] Anquez | Automatic correction and simplification of geological maps and cross-sections for numerical simulations[END_REF], due to the limited resolution of the seismic wavefield, some mild interface modifications will only introduce negligible errors in the simulation results.

Coupling the local modifications with the evaluation of their impact by comparing to a reference model (by using wave propagation simulations), leads to the fact that the integration in an inversion code could be done in more complex settings (e.g., detecting fluid substitution during gas storage, CO 2 sequestration, or reducing the horizon position uncertainty [START_REF] Bodin | Seismic tomography with the reversible jump algorithm[END_REF]).

More generally, the local insertion of an implicit interface into an existing mesh opens interesting perspectives to address other types of geological uncertainties, such as the simulations of undetected fractures and faults, or the management of stratigraphic uncertainties. As faults and fractures do not necessarily terminate against pre-existing boundaries, the definition of the region to be modified is essential. This particular region could be defined according to the object that we want to modify or insert. In two dimensions, this could be achieved by defining a circle centered on the fracture to be inserted or more generally by a second level-set [START_REF] Cherpeau | Method for Stochastic Inverse Modeling of Fault Geometry and Connectivity Using Flow Data[END_REF][START_REF] Moës | Non-planar 3D crack growth by the extended finite element and level sets-Part I: Mechanical model[END_REF]. Finally, the zone to be remeshed could be defined by a number of rows of triangles around the fracture, so as to preserve the mesh quality during the insertion.

Although the principles discussed in this article for the local modification of a meshed geomodel in 2D could be extended in 3D, the transition from 2D to 3D geometric operations is seldom trivial. Indeed, the 2D geometric model editing operations are well understood and easily implemented, but the corresponding 3D operations often generate more complex configurations which are di cult to handle. Nonetheless, we are confident that the use of piecewise linear level-sets for local 3D mesh updating are very promising. As a matter of fact, it is already possible to discretize surfaces in 3D multi-material models using a similar workflow relying on the mmg3d code (Figure 19).

Figure 1 :

 1 Figure 1: (1) From several level-sets defining the possible interface positions between water-saturated sandstone and

Figure 2 :

 2 Figure 2: 2D triangulated models are defined by their Corners (black dots), Lines, and Surfaces (a) whose geometry is supported by points, segments and triangles (b)

Figure 3 :

 3 Figure 3: Example of a fault throw modification in which the model is completely erased and meshed again.

Figure 4 :

 4 Figure 4: (Left) 2D level-set function; (right) a corresponding meshed description: the 2D domain is equipped with a

Figure 5 :

 5 Figure 5: Insertion of an implicitly defined salt dome (b) in a geological model (a). The iso-values of the signed distance scalar field represent possible salt dome boundaries; (c) and (d) show two resulting models for scalar values 0 and -70, respectively.

Figure 6 :

 6 Figure 6: The flexibility of the proposed approach demonstrated on the insertion of an arbitrary shape in an edited region composed of surface units A 1 , B 2 , B 3 and A 4 (light blue).

Figure 7 :

 7 Figure 7: Two steps level-set discretization method in Mmg: (a) If the level-set goes exactly through one vertex, the triangle is cut into two triangles; (b) if the level-set cuts the interior of two edges, the triangle is cut into one triangle (red) on one side and two triangles (green and blue) on the other side

Figure 8 :Figure 9 :

 89 Figure 8: Two steps level-set discretization method in Mmg: (a) insertion of the line by creating new nodes and triangles; (b) remeshing of the region

Figure 11 :

 11 Figure 11: Results of MMG remeshing by default (left) and with the option -nofem activated (right). With this option,

Figure 12 :

 12 Figure 12: Reference model: the gas-water depth is z = *2, 785 m. The dashed green line represents a line of 150 receivers (one every 100 m) at the top of the domain. The lithology is presented in Table 1

Figure 14 :

 14 Figure 14: The mesh quality impacts the computation time: the highlighted triangles with the lowest quality (right) for 4 models (left) have elongated form, leading to decreasing time steps and increasing simulation run time (top).

Figure 15 :

 15 Figure 15: L2-norm waveform difference d i,z (4) for each receiver i at position x = i < 100 m as a function of contact depth difference z (results are in percentages). The color represents the L2-norm value at one receiver for one z . The z values are negative when the contact of a model is below the contact of the reference model.

Figure 16 :

 16 Figure 16: Mean of L2-norm z (5) as a function of z

)

  The likelihood function L(D, z) quantifies how equivalent a candidate model M is to the reference model M ref by comparing their seismic simulation results: D(z) at all receivers for the candidate model and D ref = D(z ref ) for the reference model. In this work, we use a classical Gaussian assumption:

Figure 17 :

 17 Figure 17: Posterior probability density function. The density is maximal around the reference depth z = *2, 785 m. The most probable gas-water contact depth can therefore be determined in the proposed settings.

Figure 18 :

 18 Figure 18: A quality threshold can be defined to avoid a high top simulation time. If the quality Q (2) is smaller than 3, the simulation time significantly increases (on a machine with a processor Intel(R) Xeon(R) CPU E5-1620 v2 @ 3.70GHz, 8 Cores).

Figure 19 :

 19 Figure 19: Insertion of an horizontal surface in multi-material model: (a, c) Initial multi material model; (b, d) The surface is inserted in the model only in the outer domain (1) of (a). The elephant is still composed of a unique region. The elephant surface is remeshed to maintain the conformity of the model.

Table

Table 1

 1

	Color map Lithology	P-wave velocity (m/s) S-wave velocity (m/s)
	1	Shale	3000	2500
	2	Sandstone	3400	2260
	3	Shale	3200	2650
	4	Sandstone	3500	2330
	5	Shale	3300	2750
	6	Sandstone (Gas sat.)	3700	2740
	7	Sandstone(W atersat.)	3800	2375
	8	Shale	3900	2600
	9	Limestone	4000	2900
	10	Sandstone	4500	3000

Table 2

 2 ) Nb triangles Quality Time (s) Nb CC

	-2,645	21,997	8.54	17,321	3
	-2,655	21,989	2.00	46,376	3
	-2,665	22,007	6.37	16,973	3
	-2,675	21,997	6.1	18,231	3
	-2,685	21,995	2.35	42,651	3
	-2,695	22,023	4.85	21,300	3
	-2,705	22,003	4.76	21,759	2
	-2,715	21,971	5.1	20,847	2
	-2,725	22,003	7.99	18,469	2
	-2,735	22,011	2.95	33,718	2
	-2,745	21,991	6.01	19,620	1
	-2,755	21,967	25.51	11,348	1
	-2,765	21,997	2.77	35,271	1
	-2,775	21,981	22.76	10,809	1
	-2,785	21,959	7.36	20,298	1
	-2,795	21,969	17.95	11,508	1
	-2,805	21,971	24.84	11,514	1
	-2,815	21,939	11.31	13,674	1
	-2,825	21,955	6.89	17,883	1
	-2,835	21,951	2.57	38,380	1
	-2,845	21,959	22.55	11,015	1

https://github.com/ring-team/LUMOS2D

https://github.com/MmgTools/mmg

https://gitlab.onelab.info/gmsh/gmsh

Acknowledgments

This work was performed in the frame of the RING project (http://ring.georessources.univ-lorraine.fr) at Université de Lorraine. We would like to thank for their support the industrial and academic sponsors of the RING-GOCAD Consortium managed by ASGA. We also thank TOTAL for additional support for the PhD of Capucine Legentil. The software corresponding to this paper is available as a python and C++ code, which significantly relies on Mmg and Gmsh. We also acknowledge Inria and the Mmg Consortium for the hou10ni and Mmg libraries. The authors would also like to thank Hélène Barucq, Julien Diaz for their help.

Code availability section

Name of the code/library : LUMOS2D Contact: capucine.legentil@univ-lorraine.fr License: GNU GPL v.3

Hardware requirements: Linux System Program language: Python Software required: python 3.6 or 3.7, Mmg, Gmsh, shapely, VTK, meshio (to run the jupyter notebook)

Program size: 5Mo

The source codes are available for downloading at the link: https://github.com/ring-team/LUMOS2D


A Filleron

M Simon

S Hantova

G Cambonie

T A Tran

H Marchandin

E Jumas-Bilak

Detection of coagulase-negative staphylococci in pasteurized breast milk

 tuf-PCR-TTGE All positive for CoNS  tuf-PCR-TTGE after DNAse digestion of free DNA  Detection of CoNS DNA from 84/211 samples 39.8% of samples  Number of species by sample: 1 to 7

 Identification  S. epidermidis: major detected species 64%  S. haemolyticus 13%.

 Conventional cultures All negative for

BACKGROUND MAIN OBJECTIVE MATERIALS AND METHODS

DISCUSSION -CONCLUSIONS CoNS:

 most represented group of bacteria in raw breast milk and stools of very preterm neonates [1,2]  leading cause of late-onset sepsis in preterms [3]  Similarity between CoNS isolated from stools and neonate blood cultures suggests a causal relationship [4]  Preterm infants are frequently fed with pasteurized breast milk (PBM) from milk bank  Pasteurization process should completely eliminate remaining CoNS Search for CoNS in pasteurized breast milk by cultivation and cultivation-independent methods to investigate the possible role of PBM as a vehicle for the transmission of CoNS to the neonates fed with PBM.

 603 samples of PBM from milk bank  Conventional and optimized culture of PBM, and rpoB gene-based identification of CoNS [5]  1 mL of sample inoculated onto sheep blood agar (bioMérieux) and incubated at 37°C during 48 h  samples inoculated on trypticase soy agar (TS, Difco) and Schaedler broth (bioMérieux) incubated at 37°C for 10 days, TS agar slant being daily inoculated by turning, and minimal medium (Nutrient agar, Difco) into a tube, incubated for 20 days at 37°C (500 μl of milk in each of the media).

 Cultivation-independent method after DNAse digestion of free DNA: tuf gene based-PCR followed by temporal temperature gradient gel electrophoresis (tuf-PCR-TTGE) [6] Keywords


Serena E Dool
email: serena.dool@gmail.com

Gabriele Uhl
email: gabriele.uhl@uni-greifswald.de

INRAE Cbgp

Knight: First record of a false widow spider Steatoda nobilis preying on a pipistrelle bat in Britain"

A recent publication in Ecosphere claimed to report the first case of a predation event by a theridiid spider on a bat globally and the first case of a member of the genus Steatoda preying on a mammal globally, and the paper concludes with possible implications for public health. Predation is typically understood to mean "capture, kill, and eat". However, none of these aspects had been observed. Moreover, key aspects of bat biology were not taken into account, including that the bat in question was a flightless newborn, the capture of which would require some explanation. We revisit this interesting observation considering both bat and spider ecology.

1 Introduction Dunbar et al. (2022a) report an interesting observation of a young dead bat, as well as subsequently a live adult bat of the same species, found in the web of a false noble widow spider (Steatoda nobilis, Theridiidae) in the UK. The web was situated below the entrance to a bat (Pipistrellus sp.) maternity roost. The authors state that the dead young bat was "evidently wrapped in silk" and "a region of dark purple coloration and [being] slightly shrivelled is indicative of the spider feeding on the bat". The adult bat was found entangled in the web with "no clear sign of interference by the spider". The authors state that at least the first observation indicates predation. However, there is no evidence presented to support this claim since the capturing, killing, and consuming of the bat were not documented. Considering bat biology, it seems somewhat improbable that predation occurred, although this remains a possibility.

Definition of predation

Although some very broad definitions of predation exist, e.g. "the killing by one living organism of another for food", which thus includes herbivory, the more usual sense of the word is the act of an animal killing and eating other animals, e.g. "The natural preying of one animal on others; the behaviour of a predator", where a predator is "An animal that naturally preys on others; an animal that habitually catches and eats prey" (Shorter Oxford English Dictionary, 2007); "the interaction between species' populations in which one organism, the predator, obtains energy (as food) by consuming, usually killing, another organism" [START_REF] Allaby | A Dictionary of Ecology[END_REF]; or "an interaction between two populations of animals in which one (the predator) hunts, captures, and kills the other (the prey) for food" [START_REF] Hine | A dictionary of biology. 8th edition[END_REF]. [START_REF] Nyffeler | Bat predation by spiders[END_REF] comprehensively reviewed known cases of bats preyed upon by spiders, in which a strict distinction was made between intentionally killing and eating (predation) and cases where bats became entangled in spider webs and died of exhaustion, starvation, dehydration, or hypothermia (non-predation deaths). The latter cases might result in scavenging or might not. Here we assess the observations reported in Dunbar et al. (2022a) for the three aspects of predation (capture, killing, eating) with respect to bat and spider biology.

Published by Copernicus Publications on behalf of the European Ecological Federation (EEF).

Bat biology

In the six-page article, neither bat nor spider biology is substantially discussed. However, this is important context for interpreting what was observed and surmising what may have taken place. Temperate species of bat form maternity colonies during the summer months. These groups are typically formed by adult females and young bats of both sexes who are not yet reproductively active. In late June-early July the pregnant females give birth, usually to a single pup [START_REF] Dietz | Bats of Britain, Europe & Northwest Africa[END_REF]. Bat pups are born with closed eyes and without hair. In Pipistrellus sp. (P. pipistrellus/P. pygmaeus) the eyes first open at 4-6 d old and the first downy fur appears at 6-8 d when wing-stretching behaviour also begins [START_REF] Hughes | Ontogeny of 'true' flight and other aspects of growth in the bat Pipistrellus pipistrellus[END_REF][START_REF] Kleiman | Maternal care, growth rate, and development in the noctule (Nyctalus noctula), pipistrelle (Pipistrellus pipistrellus), and serotine (Eptesicus serotinus) bats[END_REF]. Newborns are not fully furred until 2 weeks of age, when they begin their first attempts at flight [START_REF] Hughes | Ontogeny of 'true' flight and other aspects of growth in the bat Pipistrellus pipistrellus[END_REF][START_REF] Kleiman | Maternal care, growth rate, and development in the noctule (Nyctalus noctula), pipistrelle (Pipistrellus pipistrellus), and serotine (Eptesicus serotinus) bats[END_REF]. It takes 1 month before they are close to adult size and weight and can fly in a relatively accomplished and independent manner (i.e. sustained flight, with some manoeuvrability) [START_REF] Racey | Feeding ecology of Pipistrellus pipistrellus (Chiroptera: Vespertilionidae) during pregnancy and lactation, I. Foraging behaviour[END_REF]. However, the "clumsy" and "uncoordinated" nature of such early flights is often noted during at least the first 2 weeks of active flying [START_REF] Hughes | Ontogeny of 'true' flight and other aspects of growth in the bat Pipistrellus pipistrellus[END_REF][START_REF] Kunz | On the Efficiency of the Tuttle Bat Trap[END_REF][START_REF] Racey | Feeding ecology of Pipistrellus pipistrellus (Chiroptera: Vespertilionidae) during pregnancy and lactation, I. Foraging behaviour[END_REF]. Therefore, young bats, including Pipistrellus sp., do not fly independently until at least 1 month of age, by which time they are fully furred. Before they are capable of flying, newborns remain in the roost or attached to their mother. Figure 2 in Dunbar et al. (2022a) clearly shows a naked newborn pup which would be completely incapable of flying. There is no explanation given of how the capture of a flightless newborn could occur. It is likely that no explanation was given as the authors mistakenly believed that the newborn was capable of flying, as stated in a video interview related to the article: "The first one to be caught was actually a juvenile specimen, so fairly small, that was probably trying to learn to fly for the first time" (University of Galway, 2022).

Aside from being born naked, very young bats are incapable of thermoregulation [START_REF] Racey | Life-history and reproductive strategies of bats[END_REF], which is why mothers choose warm sites for maternity colonies, exhibit clustering behaviour, and either return often to feed and warm their young or bring them with them on foraging trips. Mothers carrying flightless young during foraging bouts or when switching roosts is a well-documented behaviour, including in Pipistrellus sp. [START_REF] Bartonička | Roost switching and activity patterns in the soprano pipistrelle, Pipistrellus pygmaeus, during lactation[END_REF][START_REF] Hughes | The flight of pipistrelle bats Pipistrellus pipistrellus during pregnancy and lactation[END_REF]. Pups of this age that are abandoned or lost are very vulnerable and quickly perish. No evidence is presented that the spider killed the newborn bat. Given the brief duration of such events, it may be that they are rarely witnessed (see [START_REF] Nyffeler | Bat predation by spiders[END_REF]. However, considering the extreme vulnerability of a newborn bat outside without its mother, a rapid death would be expected from exposure alone. Dunbar et al. (2022a) state that the slightly shrivelled appearance and dark purple colouration of the pup are indica- tions that the spider fed on the bat. However, the natural appearance of such pups is wrinkled and pink-purple-black in colour (Fig. 1). A pup left exposed for even a short period might very well curl up tightly as a means to prevent loss of heat and dehydration and indeed appear very wrinkled/shrivelled once dead. The photo quality (Fig. 2b in Dunbar et al., 2022a) does not make it easy to assess the state of the bat pup with great clarity, nor can it be verified that the bat pup was wrapped in silk. Overall, there does not appear to be anything atypical about the position, state, or colouration of the dead newborn bat in the photos provided. No evidence is presented that the spider fed on the newborn bat, and the spider was not seen near the bat. The extraintestinal feeding style of spiders requires that the spider punctures the skin of the bat with its chelicerae and releases saliva to liquefy the bat's underlying tissue and then ingest it. This likely requires several hours of handling time as was observed for an orbweaver feeding on a proboscis bat [START_REF] Timm | Orb-Weaving Spider, Argiope Savignyi (Araneidae), predation on the proboscis bat Rhynchonycteris naso (Emballonuridae)[END_REF] the very same feeding event was observed by Gabriele Uhl on 25 July 2005). Puncture marks were not reported in the current account.

Spider biology

As noted by [START_REF] Nyffeler | Bat predation by spiders[END_REF] [START_REF] Dugon | Occurrence, reproductive rate and identification of the non-native noble false widow spider Steatoda nobilis (Thorell, 1875) in Ireland[END_REF]. However, Steatoda silk has been observed to be strong [START_REF] Snazell | The theridiid spider Steatoda nobilis (Thorell, 1875) in Britain[END_REF], and the analysis of the properties of silk of another theridiid spider, Latrodectus hesperus, demonstrated that the silk is at the upper range of extensibility, toughness, and strength (Blackledge et al., 2005;[START_REF] Swanson | Spider dragline silk: correlated and mosaic evolution in highperformance biological materials[END_REF]. The three-dimensional web structure of theridiid spiders consists of a cobweb with sticky gumfooted threads that lead to the ground (used to catch walking prey) and a supporting structure (used to catch flying prey) (Argintean et al., 2006;[START_REF] Benjamin | Steatoda triangulosa and comments on phylogenetic implications (Araneae: Theridiidae)[END_REF]. The gumfoot threads allow these spiders to capture prey much larger than themselves. The gumfoot threads detach from the substrate when an animal walks by and comes into contact with them. When the prey is small, the gumfoot thread pulls it upwards. Bigger prey is not lifted by a single thread. Rather, it has been observed that theridiid spiders run towards large prey and add more threads before biting the prey, which has also been reported for Steatoda nobilis [START_REF] Snazell | The theridiid spider Steatoda nobilis (Thorell, 1875) in Britain[END_REF]. Biting before wrapping was observed in a Steatoda triangulosa dealing with a 6 cm long geckonid lizard [START_REF] Vitkauskaite | Vertebrate prey capture by Latrodectus mactans (Walckenaer, 1805) and Steatoda triangulosa (Walckenaer, 1802) (Araneae, Theridiidae) provide further insights into the immobilization and hoisting mechanisms of large prey[END_REF]. In an experimental study, it was demonstrated that Steatoda paykulliana and S. triangulosa attach pre-tensioned silk threads to the prey to hoist such large prey off the ground and into the web where they are fed on [START_REF] Greco | How spiders hunt heavy prey: the tangle web as a pulley and spider's lifting mechanics observed and quantified in the laboratory[END_REF]. Some Latrodectus spiders have been reported to catch and consume mice [START_REF] Nyffeler | Black widow spiders, Latrodectus spp. (Araneae: Theridiidae), and other spiders feeding on mammals[END_REF]. In conclusion, the specific silk properties together with the specific behavioural strategies combined with the potent venom found in the Latrodectus and Steatoda species studied thus far (Dunbar et al., 2022b;Garb and Hayashi, 2013) render it likely that these spiders do prey on vertebrates, including capture, killing, and consumption at least occasionally. Nevertheless, the flightless newborn bat neither walked under nor flew into the web and killing was not documented, nor were bite marks, envenomation, or eating reported.

(Re-)interpretation of the observations

A dead newborn bat pup was found in a spider web one morning during the first week of July 2021. The following morning, the dead pup had fallen to the ground beneath the web and a live bat was found entangled in the web, was released by the observer, and thence returned to the maternity colony. It is not certain that the second bat in the spider web was an adult. No evidence to support this was presented (e.g. fused phalangeal epiphyses, reproductive status). It is plausible that it was an older juvenile. A plausible scenario which could account for the newborn bat in the web is the following: a mother bat left the roost during the night with her newborn. She became encumbered by the spider web and dropped her pup into it. She may or may not have tried to retrieve her pup, a behaviour which is relatively common in bats [START_REF] Kunz | Parental care and postnatal growth in the Chiroptera[END_REF]. Other possibilities exist, including that the pup died in the roost shortly after its birth and simply tumbled out of the roof space into the web. It might have been an obstacle in the web, and the spider could have discarded it during the following night without feeding on it. It also remains possible that the bat fell in the web (from the roost/from its mother) and was killed by the spider and that the spider fed on it. We do not exclude the possibility that feeding occurred on the newborn bat, possibly as by-catch. However, no compelling evidence was presented to support a predation event. Gut content metabarcoding in combination with detailed behavioural observations would help to clarify if theridiid spiders living in close proximity to a bat roost regularly prey on bats.

Figure 1 .

 1 Figure 1. Adult female common pipistrelle (Pipistrellus pipistrellus) with her newborn pup (Burgdorf, Germany, June 2014). Photo by Bernd Rose, NABU Burgdorf, used with permission.

  Dool and G. Uhl: Bat predation by Steatoda nobilis is not verified 99 whose main prey are flying insects, with vertebrate taxa occurring occasionally as by-catch. The diet of Steatoda nobilis is no exception, with invertebrate prey being the main food type

, the documented cases of spiders preying on bats involve species Web Ecol., 22, 97-100, 2022 https://doi.org/10.5194/we-22-97-2022 E.

Acknowledgements. We are grateful to the editors for providing us with the opportunity to submit our comment to Web Ecology. The journal which published the original article offered to publish our comment but demanded a EUR 1700 article processing charge (APC) for it.

Review statement. This paper was edited by Daniel Montesinos and reviewed by Roland Brandl and one anonymous referee.

Data availability.

No data sets were used in this article.

Author contributions. Both authors wrote the article together.

Competing interests. The contact author has declared that neither of the authors has any competing interests.

Disclaimer. Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.