Elodie Lesage

Frédéric Schmidt

François Andrieu

Hélène Massol

Constraints on effusive cryovolcanic eruptions on Europa using topography obtained from Galileo images

Keywords: Europa, Icy moon, cryovolcanism, DEM, shape-from-shading, uncertainties

Images of Europa's surface taken by the Galileo Solid State Imager (SSI) show smooth features measuring a few kilometers, potentially resulting from eruptions of low-viscosity material such as liquid cryomagma.

We estimated the volume of four of these smooth features by producing Digital Elevation Models (DEMs) of four Galileo/SSI images. We used the shape-from-shading technique with special care to estimate the uncertainties on the produced DEMs and estimated feature volumes to be between (5.7±0.9)×10 7 m 3 and (2.7±0.4)×10 8 m 3 . We discussed the implications for putative sub-surface liquid reservoir dimensions in the case of eruptions induced from freezing reservoirs. Our previous cryovolcanic eruption model was improved by considering a cycle of cryomagma freezing and effusion and by estimating the vaporized cryolava fraction once cryolava spreads onto Europa's surface. Our results show that the cryomagma reservoirs would have to be relatively large to generate these smooth features (1 to 100 km 3 if the flow features result from a single eruption, and 0.4 to 60 km 3 for the full lifetime of a reservoir generating cyclic eruptions). The two future missions JUICE (ESA) and Europa Clipper (NASA) should reach Europa during the late 2020s. They shall give more information on those putative cryovolcanic regions which appear as interesting targets that could provide a better understanding of the material exchanges between the surface, sub-surface and ocean of Europa.

Introduction

The Jovian moon Europa is believed to hide a global liquid water ocean under its ice crust [START_REF] Khurana | Induced magnetic fields as evidence for subsurface oceans in Europa and Callisto[END_REF][START_REF] Pappalardo | Does Europa have a subsurface ocean? Evaluation of the geological evidence[END_REF]. This ocean is predicted to be in contact with a silicate mantle, which could allow the chemical exchanges needed to create a rich habitable environment [START_REF] Kargel | Brine volcanism and the interior structures of asteroids and icy satellites[END_REF][START_REF] Kargel | Europa's crust and ocean: Origin, composition, and the prospects for life[END_REF]. The habitability of Europa's ocean depends on chemical conditions and equilibrium in it [START_REF] Vance | Geophysical investigations of habitability in ice-covered ocean worlds[END_REF], but for now, it remains impossible to directly sample. Also, water coming directly from the ocean seems unlikely to erupt at its surface because of the very high pressure required for it to ascend through the whole ice crust [START_REF] Gaidos | Tectonics and water on Europa[END_REF][START_REF] Manga | Pressurized oceans and the eruption of liquid water on Europa and Enceladus[END_REF][START_REF] Rudolph | Fracture penetration in planetary ice shells[END_REF]. The strong tides generated by Jupiter in Europa's ice crust could be at the origin of ice local melting [START_REF] Tobie | Tidal dissipation within large icy satellites: Applications to Europa and Titan[END_REF][START_REF] Mitri | A model for the temperature-dependence of tidal dissipation in convective plumes on icy satellites: Implications for Europa and Enceladus[END_REF][START_REF] Vilella | Tidally heated convection and the occurrence of melting in icy satellites: Application to Europa[END_REF]. Such melted reservoirs appear to be good candidates to host life forms if they were able to remain dormant in the ice between melting episodes [START_REF] Gaidos | Tectonics and water on Europa[END_REF].

Identifying the geological features emplaced during eruptions of liquid water could give information on the location of the terrains that are most likely to show biosignatures, which can be useful to the two upcoming missions JUICE (ESA) and Europa Clipper (NASA). In [START_REF] Lesage | Cryomagma ascent on Europa[END_REF], we demonstrated the possibility of erupting liquid water from sub-surface freezing reservoirs. Here, we propose to test an improved version of our previous model against Galileo data. First, we select images showing geological features that could result from the eruption of liquid water/brines at the surface, i.e. the smooth plains. We then generate the Digital Elevation Model (DEM) of the chosen features and measure their volumes. We then use these results to constrain the volume and depth of the source reservoirs.

The highest resolution images of Europa were acquired during the Galileo mission with the Solid State Imager (SSI) [START_REF] Belton | The Galileo Solid-State Imaging experiment[END_REF]. These images have shown a geologically active surface, characterized by a wide variety of features [START_REF] Greeley | Europa: Initial Galileo geological observations[END_REF][START_REF] Greeley | Geologic mapping of Europa[END_REF]. The low crater density found on Europa demonstrates the vigorous resurfacing processes taking place on this moon, making its surface one of the youngest in the solar system, with an age under 90 Myrs [START_REF] Zahnle | Cratering rates in the outer Solar System[END_REF]. Plate tectonics-like processes were thought to frequently recycle the icy surface [START_REF] Greeley | Europa: Initial Galileo geological observations[END_REF][START_REF] Kattenhorn | Evidence for subduction in the ice shell of Europa[END_REF] but numerical modeling of the ice shell shows that a global Earth-like plate tectonics is unlikely on Europa [START_REF] Howell | Can Earth-like plate tectonics occur in ocean world ice shells?[END_REF]. Cryovolcanic activity may therefore contribute covering some older terrains with fresh material.

A wide range of local-scale geological features is observed at Europa's surface, such as chaos, lenticulae, domes, pits, and ridges [START_REF] Greeley | Europa: Initial Galileo geological observations[END_REF][START_REF] Greeley | Geologic mapping of Europa[END_REF]. Several formation mechanisms have been proposed to explain the formation of these features and invoke, in most cases, a diapiric ascent of warm ice [START_REF] Head | Europa: Morphological characteristics of ridges and triple bands from Galileo data (E4 and E6) and assessment of a linear diapirism model[END_REF][START_REF] Sotin | Europa: Tidal heating of upwelling thermal plumes and the origin of lenticulae and chaos melting[END_REF][START_REF] Fagents | Considerations for effusive cryovolcanism on Europa: The post-Galileo perspective[END_REF][START_REF] Schenk | Topographic variations in chaos on Europa: Implications for diapiric formation[END_REF][START_REF] Quick | Heat transfer of ascending cryomagma on Europa[END_REF], or a direct link with the ocean [START_REF] Greenberg | Chaos on Europa[END_REF][START_REF] Greenberg | Europa's dynamic icy crust[END_REF]. More recent studies show that the formation of some of these features could be related to the presence of sub-surface liquid reservoirs, as for lenticulae (Michaut and Manga, 2014;[START_REF] Manga | Formation of lenticulae on Europa by saucer-shaped sills[END_REF], chaos [START_REF] Schmidt | Active formation of 'chaos terrain' over shallow subsurface water on Europa[END_REF] and double ridges [START_REF] Johnston | Formation of ridges on Europa above crystallizing water bodies inside the ice shell[END_REF]. Numerical models showed the possibility of generating locally melted zones within Europa's ice crust due to the combination of convection and tidal heating [START_REF] Sotin | Europa: Tidal heating of upwelling thermal plumes and the origin of lenticulae and chaos melting[END_REF][START_REF] Mitri | A model for the temperature-dependence of tidal dissipation in convective plumes on icy satellites: Implications for Europa and Enceladus[END_REF][START_REF] Han | Coupled convection and tidal dissipation in Europa's ice shell[END_REF][START_REF] Vilella | Tidally heated convection and the occurrence of melting in icy satellites: Application to Europa[END_REF], which could explain the formation of sub-surface molten reservoirs. [START_REF] Fagents | Considerations for effusive cryovolcanism on Europa: The post-Galileo perspective[END_REF], [START_REF] Quick | Cryovolcanic emplacement of domes on Europa[END_REF], and [START_REF] Nunez | Developing a database for candidate cryovolcanic domes on Europa[END_REF] suggest that a subset of domal lenticulae have been put in place by the eruption of cryomagma reservoirs. Thank to DEMs and numerical modeling, [START_REF] Quick | Cryovolcanic emplacement of domes on Europa[END_REF] demonstrated that domes could by put in place the effusion of viscous cryomagma. Nevertheless, thiner and smoother features observed on Galileo images have received less attention.

In their exhaustive classification of Europa's geological features, [START_REF] Greeley | Geologic mapping of Europa[END_REF] introduced the so-called "smooth plains" units, which are defined as smooth surfaces, with no visible texture, that embay or overprint preexisting terrains. [START_REF] Greeley | Geologic mapping of Europa[END_REF] proposed two models of formation for smooth plains, which are (1) the cryovolcanic emplacement of low-viscosity material (such as liquid water-based mixture) and (2) the melting of the surface due to a local heat source. Nevertheless, as we will show it in ths study, the smooth plains morphology is not consistent with surface melting only. Moreover, because of the shape of the smooth plains, which are very thin and topographically constrained, these features have been widely interpreted as the result of liquid flows on the surface [START_REF] Pappalardo | Does Europa have a subsurface ocean? Evaluation of the geological evidence[END_REF][START_REF] Fagents | Considerations for effusive cryovolcanism on Europa: The post-Galileo perspective[END_REF][START_REF] Miyamoto | Putative ice flows on Europa: Geometric patterns and relation to topography collectively constrain material properties and effusion rates[END_REF]. [START_REF] Miyamoto | Putative ice flows on Europa: Geometric patterns and relation to topography collectively constrain material properties and effusion rates[END_REF] modeled liquid flow under Europa's surface conditions and have shown that the effusion of a low viscosity material such as water or brines may create flow-like features before freezing that are consistent with the morphologies of some smooth plains. In this study, we focus on small-scale smooth plains a few kilometers wide that we call "smooth features" hereafter.

In [START_REF] Lesage | Cryomagma ascent on Europa[END_REF], we modeled the eruption of liquid cryolava from a freezing sub-surface reservoir.

In the case of a single eruption, we showed that the erupted volume of cryolava mainly depends on the reservoir volume and depth. Assuming that smooth features may be formed by the effusion of liquid cryolava to the surface, measuring the volume of cryolava flows will provide constraints on the volume and depth of the cryomagma reservoirs using the framework and results of [START_REF] Lesage | Cryomagma ascent on Europa[END_REF]. To estimate the volume of putative flow features, we generate DEMs of Europa's surface using the AMES Stereo Pipeline (ASP) [START_REF] Beyer | The Ames Stereo Pipeline: NASA's open source software for deriving and processing terrain data[END_REF]. This open-source tool has been used previously by several authors to generate DEMs of Europa's surface: for instance, [START_REF] Schenk | Topographic variations in chaos on Europa: Implications for diapiric formation[END_REF] and [START_REF] Schmidt | Active formation of 'chaos terrain' over shallow subsurface water on Europa[END_REF] have generated DEMs of chaotic terrains; [START_REF] Dameron | Europan double ridge morphology as a test for hypothesized models of formation[END_REF] proposed a statistical study of double ridge morphology; a large database that includes DEMs of putative cryolava domes is currently being built by [START_REF] Nunez | Developing a database for candidate cryovolcanic domes on Europa[END_REF]. Here we develop our own cryovolcanic features database based on two morphologic criteria (see sec. 2.1): (i) a thin feature confined by the surrounding ridges and (ii) a smooth appearance with no blocks showing a ridged texture.

We selected four images with smooth features that are compatible with liquid flows, and we calculate the volume of these features. We also calculate the uncertainty on the DEMs based on a sensitivity study (see supplementary materials for details).

Methods

Selection of cryovolcanic features : criteria

To obtain the volume of cryovolcanic eruptions putative products, we generate DEMs of features that might have possibly been put in place during the effusion of liquid water at the surface. These features, called hereafter "smooth features", are a few kilometers wide and have common characteristics with smooth plains previously defined by [START_REF] Greeley | Geologic mapping of Europa[END_REF]. To select the most relevant features in the Galileo dataset and avoid confusion with other features, such as chaos, we choose two major criteria to define the putative cryovolcanic features:

1. Thin features occupying topographic lows. Because we are interested in features that may have been emplaced during the effusion of liquid onto the surface, and as modeled by [START_REF] Miyamoto | Putative ice flows on Europa: Geometric patterns and relation to topography collectively constrain material properties and effusion rates[END_REF] ("A-type" features in their study), the resulting features are expected to be thin, with no particular domal shape. These features occupy the topographic lows, and are constrained by pre-existing relief.

On Europa, these features are typically delimited by ridges.

2.

A smooth appearance with no visible blocks of older surface. Some features at Europa's surface present a smooth matrix, but also contain blocks presenting a texture similar to the pre-existing, ridged plains. These features, named chaos, are identified as the result of local melting and disruption of the surface, followed by its freezing [START_REF] Greenberg | Europa's dynamic icy crust[END_REF][START_REF] Figueredo | Geology and origin of Europa's "Mitten" feature (Murias chaos)[END_REF][START_REF] Schmidt | Active formation of 'chaos terrain' over shallow subsurface water on Europa[END_REF].

Here we study smooth features assumed to have been emplaced by effusion of liquid cryolava onto the surface, which is a completely different process. Nevertheless, we do not rule out the possibility of forming new blocks of ice, either during the flowing of water at the surface because of the very low temperature at Europa's surface, or due to the transport of small blocks of ice due to stopping as fluid ascends toward the surface. Hence, we do not exclude features containing relatively small blocks (covering a surface fraction less than ∼ 10% of the feature), as long as they do not have a ridged surface.

We looked at all SSI images with a resolution higher than ∼100 m/px and found 4 regions fulfilling these criteria (see Fig. 1). A fifth interesting smooth plain was identified (image 8613r) but unfortunately, we did not manage to produce its DEM because of a large projected shadow on the image.

DEM generation

A few steps are required to obtain a ready-to-use DEMs . After calibration, noise-filtering, and mapprojection of the images using ISIS 3 (https://isis.astrogeology.usgs.gov, see details in supplementary materials, section 1 and in the "ISIS" part of flowchart in Fig. 2), we use the NASA AMES StereoPipeline to generate the DEMs [START_REF] Beyer | The Ames Stereo Pipeline: NASA's open source software for deriving and processing terrain data[END_REF], "AMES StereoPipeline" part of the flowchart in Fig. 2). These DEMs are then used to calculate the smooth feature volumes ("QGis" part of the flowchart in Fig. 2).

Two main tools can be used with ASP in order to produce DEMs: the Stereo tool, based on the stereoscopic analysis, and the Shape from Shading (SfS) tool, based on the photoclinometric principle. The Stereo tool produces robust DEMs because this method is based on the correspondences between pixels of two (or more) images [START_REF] Beyer | The Ames Stereo Pipeline: NASA's open source software for deriving and processing terrain data[END_REF]. Stereoscopy requires at least two images of the same terrain, each taken from a different point of view (at least a few degrees of difference) and illuminated from a similar direction. These two criteria are very limiting for the use of the Galileo SSI data: due to the limited number of high-resolution images, image pairs that satisfy these two criteria are extremely rare. For this reason, we cannot use the Stereo tool to produce DEMs of the smooth plains, and we focus on the SfS. The ASP Stereo tool is expected to give more robust results, but SfS produces 3 to 5 times better resolved DEMs [START_REF] Nimmo | Stereo and photoclinometric comparisons and topographic roughness of Europa[END_REF], which allows a more precise study of small scale elevation variations. Also, even though SfS is only able to give relative heights and cannot be used to infer an absolute elevation, this is not limiting in our case as we only want to know the height of cryovolcanic features relative to the surrounding terrain.

To produce the DEMs presented here, we use the SfS tool based on the photoclinometric principle [START_REF] Alexandrov | Multiview Shape-From-Shading for planetary images[END_REF]. SfS uses light intensity of each pixel of an image to infer the surface topography.

In fact, the mean brightness of each surface facet (i.e. the footprint of each pixel on the real surface) is a function of the angle between the sunlight incidence direction and the direction normal to the facet. Based on a known solar elevation and azimuth, SfS computes the slope of each image pixel from its brightness.

Then, it integrates pixel slopes to give the terrain shape. Numerically, this process is done by minimizing the following cost function [START_REF] Alexandrov | Multiview Shape-From-Shading for planetary images[END_REF]:

[I (h (x, y)) -T.A.R (h (x, y))] 2 + µ 2 ∇ 2 h (x, y) 2 +λ 2 [h (x, y) -h 0 (x, y)] 2 dx dy (1)
where h(x, y) is the function describing theterrain topography, I (h(x, y)) is the image viewed by the camera reconstructed using the terrain h(x, y), T is the image exposure, A is the terrain albedo (considered to be constant on the whole image) and R (h (x, y)) is the reflectance. h 0 (x, y) is an "initial guess" terrain (also called "apriori" in the SfS tool) from which the minimization algorithm starts to calculate the real shape of the terrain. This initial guess is provided by the user in the form of a geolocalized DEM. The best initial guess would be a DEM generated with the stereoscopic technique, which could be refined by the SfS tool.

Here we do not possess such a DEM, so we provide a flat DEM to SfS to start the minimization. Finally, µ and λ are two positive coefficients chosen by the user and controlling respectively the smoothing of the DEM and the weight of the initial guess terrain.

The first term of the cost function constrains the brightness and ensures that the simulated light intensity fits with the image recorded by the camera. This term depends on the reflectance model used, which can be chosen by the user. The Lambertian, LunarLambert and Hapke models are available to model icy surfaces.

In accordance with a recent photometric study [START_REF] Belgacem | Regional study of Europa's photometry[END_REF], we use the Hapke model with the following parameters: ω=0.9, b=0.35, c=0.65, B 0 =0.5, h=0.6. The influence of the photometry on the DEMs is showed in supplementary materials (section 2). The relative uncertainty due to photometry has been estimated at ±10% of the measured volume (see supplementary materials).

The second term of the cost function controls the smoothness of the output DEM by minimizing the second-order derivative of the slope on each point. The smoothness coefficient µ is chosen by the user.

Theoretically, the higher µ, the smoother the DEM, making the small-scale details less visible and flattening higher relief features. Nevertheless, we noticed that very small values of µ also produce flattened DEMs. To avoid an extremely flattened DEM, we typically choose 1<µ<10. We tested several values of this coefficient for each image and concluded that there is no ideal value of µ that can be used for all the images. In fact, the smoothing effect controlled by µ depends on the terrain roughness and therefore differs for each studied image [START_REF] Alexandrov | Multiview Shape-From-Shading for planetary images[END_REF]. For each image, we need to test several µ values to keep the most appropriate one. The specific effects of µ variation on the produced DEM are shown in the supplementary materials (section 2). For small features, we found that the DEM is always consistent with the feature heights estimated from shadows (see supplementary materials section 2), in a large domain of the smoothness parameter µ. The relative uncertainty due to the smoothness coefficient has been estimated at ±5% of the measured volume.

Finally, the third term of the cost function describes the difference between the calculated DEM and the initial guess h 0 (x, y) given to SfS. Here, we do not have an initial guess of the terrain elevation and SfS is the only tool used to generate the DEMs. Hence, we use a flat terrain at 0 elevation as the initial guess, and we set the parameter λ=0. By doing so, SfS should not depend on this flat initialization DEM during the minimization iterations.

Finally, the total relative uncertainty of the volumes measured on the DEMs is ±15% (see supplementary materials).

Volume estimation

In order to estimate the cryolava volume erupted during smooth plain emplacement, the first step is to calculate the total smooth feature volume from the DEM. This process is detailed in the next section. With the example of image 5452r, we show in section 2.3.2 that it is necessary to invoke subsidence to explain the morphology of the smooth features. To obtain a more realistic cryolava volume, we also take into account the estimated underlying pre-existing topography (ridges in our case, coupled with terrain subsidence). The hypothesis and methods used to estimate the pre-existing topography are detailed in section 2.3.3.

GIS processing: simple approach

We first calculate the smooth feature volume using QGis (QGIS Development Team, 2019) based on a simple approach. The idea is to subtract the pre-existing terrain beneath the cryolava flow from the feature itself in order to obtain the volume of cryolava. A few steps are necessary to obtain this result and are summarized in the "QGis" part of the flowchart given in Fig. 2. First, the smooth feature needs to be delimited. We draw the shape of the feature and exclude ridges that might intersect it. Then, the preexisting terrain beneath the cryovolcanic feature is inferred from topography surrounding the feature. In order to estimate it, we choose several reference points in the valleys around the smooth feature, in which the cryolava has flowed, and we use the interpolation tool in QGis to create a layer that approximates the terrain elevation under the feature. Once the pre-existing terrain elevation is subtracted from the DEM, we sum the height of all the pixels composing the feature and multiply it by the area of a pixel, so we obtain the volume of the cryolava flow. This volume is called "measured volume" (noted V measured ) in the following sections.

Subsidence and thermal erosion: example of image 5452r

Image 5452r was taken by the Galileo spacecraft during its fourth orbit with a resolution of approximately 30 m/px. It shows a very smooth circular feature. This feature was first presented by [START_REF] Head | Cryovolcanism on Europa-Evidence for the emplacement of flows and related deposits in the E 4 region(5 N, 305 W) and interpreted eruption conditions[END_REF] and described by [START_REF] Pappalardo | Does Europa have a subsurface ocean? Evaluation of the geological evidence[END_REF] as a "smooth deposit, probably emplaced as a cryovolcanic eruption of low-viscosity material, perhaps liquid water". This feature interpretation is now accepted in literature because of its morphology [START_REF] Fagents | Considerations for effusive cryovolcanism on Europa: The post-Galileo perspective[END_REF][START_REF] Miyamoto | Putative ice flows on Europa: Geometric patterns and relation to topography collectively constrain material properties and effusion rates[END_REF].

Some impact craters are visible on the smooth feature, especially a large one, almost centered. We assume that all craters have been formed after the emplacement of the smooth feature and have no relation with its formation. The first reason is because other craters are also present around the smooth feature and do not seem to interact with the older terrain (no melting, no particular ejecta, etc.), and the second reason is because if the smooth feature has been formed by melting after the impact, the impact crater itself would not be visible.

Image 5452r and its DEM generated with SfS are shown respectively in Fig. 3a and3b. We determine the smooth feature edges (in blue on Fig. 3b) using the DEM. In fact, in the particular case of image 5452r, the sunlight comes from the east side of the image, which is the direction perpendicular to the surrounding ridges and thus makes the east-west oriented lobes hardly visible on the raw image. Nevertheless, the edges of the feature are well visible on the DEM, more specifically in the valleys between the ridges. We also choose five points located in valleys adjacent to the feature (see the black crosses on Fig. 3b). These points are used to estimate the pre-existing terrain elevation beneath the feature (see dashed lines in Fig. 4). By interpolation, we obtain a surface that is subtracted from the smooth feature. By doing this, we can calculate the feature volume itself. DEMs of the three other images used in this study are given in supplementary materials, section 3.

We propose some geomorphological interpretations of this smooth feature. First of all, the smooth feature, mostly in white on the DEM, has an elevation of around 0 m, which is higher than the surrounding valley bottoms, that have a negative elevation. At the locations indicated by the blue arrows in Fig. 3, one can see the feature edges, filling the bottom of the valleys. This is in agreement with the hypothesis that the smooth feature is possibly made of low-viscosity material added on the preexisting terrain, flowing between the ridges.

Moreover, on the north, south, and west of the smooth feature, we can see that some ridges are not covered by the smooth material, which indicates that the putative flow had a relatively low viscosity [START_REF] Miyamoto | Putative ice flows on Europa: Geometric patterns and relation to topography collectively constrain material properties and effusion rates[END_REF].

Finally, at the central region of the smooth feature, ridges are not visible. This was unexpected because a flow that covers a pre-existing terrain should theoretically have an elevation higher than the pre-existing terrain. But here, the surrounding ridges have an elevation higher than the smooth feature. Two scenarii can be put forward to explain this effect and are illustrated in Fig. 4 (see the two topographic profiles numbered 1 and 2 in Fig. 3b, 3c and 3d). The emplaced material could generate a local subsidence of the ice crust (model A) or it could melt/erode the surface as it flows (model B). Both cases, or a mix of them (model C), might explain the lower topography of the pre-existing terrain especially near the center.

Model A detailed in Fig. 4a could result from the local subsidence of the ice crust due to the presence of a liquid reservoir at depth centered on the feature. The required condition to create a few kilometers of large depletion is a thin elastic layer of less than a few hundred meters beneath the reservoir [START_REF] Manga | Formation of lenticulae on Europa by saucer-shaped sills[END_REF]. With this formation model, the real pre-existing terrain beneath the smooth feature has an elevation lower than the pre-existing terrain estimated with the DEM, which is the mean level of the surrounding valleys. This has to be taken into account in the volume measurement (see section 2.3.3).

In the case of model B, as shown in Fig. 4b, the ridges are subject to thermal erosion only. This could happen if a warm liquid flows onto the surface (see thermal erosion experiment from [START_REF] Kerr | Thermal erosion by laminar lava flows[END_REF]. In this case, the flowing liquid would be composed of a mixture of cryolava and molten terrain. After freezing, the molten terrain would return to a frozen state with a density similar to its original state, so finally, the net volume change due to thermal erosion would be null. In the case of image 5452r, the smooth feature is around 0 m elevation, which is also the mean elevation of the surrounding terrain based on the initial guess terrain chosen for the DEM generation. This would mean that in the case of terrain melting and freezing, no material coming from the interior was added onto the surface and thus only heat transfer is responsible for the feature emplacement. This is not in agreement with our volume measurement results (see section 4.1), so terrain melting only cannot explain the smooth feature morphology, and we do not consider model B hereafter.

Model C is a combination of thermal erosion and local subsidence, as illustrated in the sketch from Fig. 4c.

As for model A, the real pre-existing terrain is lower than the estimated one, implying a difference between the flow volume measured from the DEM and the real cryolava volume coming from the reservoir. This has to be taken into account and is detailed hereafter.

From measured volumes to cryolava volumes

On Europa, most of the surface is ridged (see images in Fig. 1 and [START_REF] Greeley | Geologic mapping of Europa[END_REF]. Thus, the simple method proposed in section 2.3.1 may overestimate the volume of actual erupted material as it also takes ridges volume into account. As we are interested in the actual volume of cryolava erupted onto the surface, we propose to use a volume factor α V which expresses the actual cryolava volume V cryolava with respect to the apparent volume V measured measured with the DEMs:

V cryolava = α V V measured (2)
The calculation of α V is illustrated in Fig. 5 and described hereafter. We use a topographic profile AB extracted from a nearby ridged plain from image 5452r to estimate α V (Fig. 5) The situation presented in Fig. 5 b and c on the left might not be realistic because, as discussed above, a cryolava flow on the surface cannot produce the observed features without local subsidence (see Fig. 4).

In
Such subsidence may significantly impact α V factor, as estimated in Fig. 5b and c on the right. We simulate the subsidence of cross-section AB to calculate the associated α V . We study two extrema to estimate the 3). The grey zone represents the pre-existing ridged terrain whereas the light blue zone stands for the smooth feature. The dark blue points delimit the edge of the smooth feature corresponding to the dark blue line in Fig. 3. The black crosses show the bottom of the nearest ridges which are used as reference elevations to interpolate the pre-existing terrain mean elevation. This interpolated pre-existing terrain elevation, represented by the black dotted line, is used to calculate the feature volume. The actual pre-existing terrain of the feature may be different from the interpolated one and is plotted with a dotted grey line. Three mechanisms are sketched: (a) surface subsidence under the feature (model A), (b) melting or thermal erosion of the pre-existing terrain (model B) and (c) a mix of these two processes (model C). A difference between the interpolated pre-existing terrain and the real one exists in the case of Model A and C. range of possible cases: (i) shallow subsidence of only the necessary height for the ridges to be embayed in the smooth feature (∼5 m, Fig. 5b), (ii) deeper subsidence of 40 m, which is the maximum subsidence modeled by [START_REF] Manga | Formation of lenticulae on Europa by saucer-shaped sills[END_REF] for subsurface reservoirs less than 10 km wide (Fig. 5c). Hence, we calculate the cross-sectional area A cryolava , which is filled with cryolava after the eruption and subduction of the terrain, represented in dotted in Fig. 5 b and c on the right. We finally calculate α V = A cryolava /A measured using these two areas. This 2D approximation is justified by the homogeneity of the ridged terrains. Based on cross-sectional areas ratio A cryolava /A measured measured on Fig. 5b andc, we obtain α V 0.8 for shallow subsidence and α V 5 for maximum subsidence. These two factors will be taken into account further.

One should also note that the putative melting or thermal erosion of the terrain during the eruption should not modify the α V value. Indeed, this would only transfer some eroded material from the underlying terrain to the smooth feature, with a null net weight balance (as explained in section 2.3.2). For this reason, the results presented within this study are not affected by putative melting or thermal erosion of the underlying terrain.

Eruption model

Pressurization by cryomagma freezing

The model described in this section was presented previously by [START_REF] Lesage | Cryomagma ascent on Europa[END_REF]. In this study, we tested the feasibility of a cryomagmatic eruption model proposed by [START_REF] Fagents | Considerations for effusive cryovolcanism on Europa: The post-Galileo perspective[END_REF] in which a cryomagma reservoir is pressurized by its freezing. We modeled a cryomagma reservoir as a spherical cavity in Europa's ice crust filled with liquid. Because of the temperature gradient between the cryomagma and the surrounding ice, cryomagma freezes from the reservoir wall toward its center, and we model the solidification front position as a function of time by solving the Stefan problem. The density contrast between liquid and solid cryomagma generates overpressure in the reservoir, thus tangential stress on the wall. When the overpressure in the reservoir is high enough, the wall breaks and a fracture can propagate toward the surface [START_REF] Lister | Fluid-mechanical models of crack propagation and their application to magma transport in dykes[END_REF][START_REF] Rubin | Tensile fracture of rock at high confining pressure: Implications for dike propagation[END_REF]. Then, cryolava can flow onto the surface until the overpressure in the reservoir has been released.

Based on this model, we can calculate the cryolava volume emitted at the surface during an eruption and the total duration of this event. These results are obtained as a function of the cryomagmatic reservoir parameters (such as its volume V and depth H) and its environment (temperature gradient in the ice crust, ice and cryomagma compositions). Results and details of this model are available in [START_REF] Lesage | Cryomagma ascent on Europa[END_REF]. andc) We measure the cross-sectional area of the smooth feature obtained with the simple approach A measured (in grey), i.e. the difference between the mean elevation of the smooth feature's top and the elevation of the interpolated pre-existing terrain multiplied by length AB. We also measure a more realistic cross-sectional area which should be filled with cryolava A cryolava (dotted). This cross-sectional area is defined as the difference between the mean elevation of the smooth feature's top and the elevation of a more realistic pre-existing terrain, i.e. a ridged terrain with subsidence of 5 to 40 m, multiplied by the length AB. We then calculate α

V = A cryolava /A measured .
We investigate the impact of putative terrain subsidence on α

V . We simulate a subsidence of the terrain in 2 cases: (b) a shallow subsidence, where the ridges tops are showing on the surface (-5 m, profile A'B'), and (c) a maximum subsidence value after Manga and Michaut (2017) (-40 m, profile A"B"). The reservoir is filled with pure or briny liquid water at isostatic pressure P 0 . (b) An initial liquid volume V i freezes up to a volume V f of ice, inducing an overpressure ∆P in the reservoir. (c) When the pressure reaches a critical value ∆Pc, the wall fractures and the pressurized liquid rises to the surface through an H long fracture. (d) Once a certain amount of liquid has erupted at the surface and the overpressure in the reservoir is released, the eruption ends. The liquid in the fracture freezes, which seals the reservoir. Finally, the reservoir returns to initial condition similar to situation (a) but with a smaller radius R i+n . Freezing of the reservoir continues such that these 4 steps repeat and form an eruptive cycle, leading to several eruptions during the reservoir's lifetime.

a) b) c) V P 0 P 0 P 0 V V P 0 liquid water surface V i V f ice P 0 +∆P ∆P<∆P c H R i d) V P 0 R i+n P 0 R i = R i+n
The model is derived for two cryomagma compositions: 1) pure liquid water and 2) a briny mixture of 81 wt% H 2 O + 16 wt% MgSO 4 + 3 wt% Na 2 SO 4 that is predicted to be close to the Europa's ocean and ice composition [START_REF] Kargel | Brine volcanism and the interior structures of asteroids and icy satellites[END_REF]. This briny mixture is assumed to be an eutectic composition in the model: when it freezes, the ice has the same composition and salt content. The recent detection of chlorides such as NaCl [START_REF] Trumbo | Sodium chloride on the surface of Europa[END_REF] or Mg-bearing chlorinated species [START_REF] Ligier | VLT/SINFONI observations of Europa: new insights into the surface composition[END_REF]) on Europa's surface indicates that the actual cryomagma composition may be a bit different from the one considered here. Nevertheless, the model results are functions of the density contrast between the liquid cryomagma and the ice, which is very similar for sulfates and chlorides [START_REF] Kargel | Brine volcanism and the interior structures of asteroids and icy satellites[END_REF][START_REF] Hogenboom | Magnesium sulfate-water to 400 MPa using a novel piezometer: Densities, phase equilibria, and planetological implications[END_REF][START_REF] Mccarthy | Solidification and microstructures of binary ice-I hydrate eutectic aggregates[END_REF][START_REF] Quick | Heat transfer of ascending cryomagma on Europa[END_REF][START_REF] Lesage | Cryomagma ascent on Europa[END_REF]. So finally, in our model, the salt content is more important than cryomagma exact composition. Here we take a mixture which contains 19% of sulfates.

Cyclic eruptions

In [START_REF] Lesage | Cryomagma ascent on Europa[END_REF], a single cryovolcanic eruption was modeled. Nevertheless, a reservoir might trigger several eruptions during its lifetime, as shown in Fig. 6. At the end of an eruption, the liquid in the fracture may freeze or tectonic stress may close the reservoir. Nevertheless, the solidification continues and the freezing front progresses toward the reservoir's center, which pressurizes the cryomagma once more and can lead to a second eruption. A reservoir might hence be able to erupt several times. Fagents (2003) predicted that cyclic eruptive events could result in morphologies where multiple flow lobes are present. Cyclic eruptions might explain the lobate morphology that can be seen at the center of the smooth feature from image 9352r (see Fig. 1d, blue arrow). As this lobate feature is located on top of the smooth plain, it might have been emplaced after the smooth plain. We do not clearly see several lobate forms on the other features, which may be explained by the limited resolution of the Galileo images. As an example, for image 5452r, the resolution is approximately 25 m/px. The solar incidence angle is ∼75°, so only features higher than ∼11 m can be seen on the image with their projected shadow, assuming that the shadow is projected on a flat terrain. Taking into account the complex processes expected to affect liquid water flowing onto Europa's surface, such as the possibility of endogenous cryolava flow and the competition between freezing and vaporizing (see [START_REF] Allison | Ice-covered water volcanism on Ganymede[END_REF][START_REF] Fagents | Considerations for effusive cryovolcanism on Europa: The post-Galileo perspective[END_REF][START_REF] Quick | Cryovolcanic emplacement of domes on Europa[END_REF] and section 3.3 for details), it is hard to determine whether the smooth features could result preferentially from a single or multiple eruptions, so here we consider these two possibilities.

To take into account the possibility of cyclic eruptions, we compare the smooth feature volumes with the total erupted volume during the lifetime of a cryomagma reservoir producing cyclic eruptions. To obtain this volume, we iterate the cryomagmatic eruption model discribed in section 3.1 as illustrated in 6. Only the reservoir fraction that remains liquid at the end of each freezing is taken into account to calculate the erupted cryolava volume. Fig. 7a shows the volume erupted at the surface during the activity lifetime of a reservoir of 10 9 m 3 located 4 km beneath the surface and filled with briny cryomagma. Each step in Fig. 7a represents an eruptive cycle, i.e. cryomagma freezing, pressurization, and eruption. Fig. 7b shows the volume of cryolava erupted at the surface during each eruption. It follows a decreasing logarithmic trend (note that y-axis is logarithmic). We can predict the volume V #i erupted during each eruption as following. For the first eruption, the erupted volume is equal to the increase in the frozen cryomagma volume. The freezing cryomagma volume is equal to nV where n is the cryomagma fraction necessary to freeze in order to trigger an eruption [START_REF] Fagents | Considerations for effusive cryovolcanism on Europa: The post-Galileo perspective[END_REF][START_REF] Lesage | Cryomagma ascent on Europa[END_REF] and V is the reservoir initial volume. The cryomagma volume that does not freeze is defined as

V 0i = V (1 -n).
The newly formed ice has a volume n ρ l ρs V where ρ l is the cryomagma density and ρ s is the ice density, so the liquid part of the reservoir is compressed to a volume Lesage et al., 2020 for details). Finally the volume change of the liquid part in the reservoir after freezing is [START_REF] Fagents | Considerations for effusive cryovolcanism on Europa: The post-Galileo perspective[END_REF][START_REF] Lesage | Cryomagma ascent on Europa[END_REF]:

V 0f = V (1 -n ρ l ρs ) (see
V #0 = V 0i -V 0f = nV ρ l ρ s -1 (3)
Then, for the second eruption, we calculate V 1i and V 1f from the remaining cryomagma after the first eruption, which has a volume

V 1i = V 0f (1 -n).
The erupted volume V #1 is:

V #1 = V 1i -V 1f = V 0f (1 -n) -V 0f 1 -n ρ l ρ s = V #0 1 -n ρ l ρ s (4)
Eq. ( 4) can be generalized for all the following eruptions, so finally the erupted volume at eruption #i can be written:

V #i = V #0 1 -n ρ l ρ s #i (5)
or: Fig. 7c demonstrates that the time between each eruption initially increases before decreasing. In fact, this freezing time is a competition between two phenomena: firstly, the thermal transfer between the warm reservoir and the cold surrounding ice slows down over time as the reservoir gradually cools, so the thermal wave propagates slower; secondly, the liquid volume decreases, which reduces the amount of cryomagma that needs to freeze to trigger an eruption. For the example given in Fig 7c, the first phenomenon dominates from eruptions #1 to #12, then the second one becomes preponderant. To explain this trend, one must compare two caracteristic velocities: the first one is the solidification front progression, and the second one is the heat transfer in the ice. These two velocities are obtained by solving the Stefan problem at the reservoir wall and thus are different for each eruption.

V #i = nV ρ l ρ s -1 1 -n ρ l ρ s #i (
Fig. 8 shows the total activity lifetime of a reservoir, considering that one cryomagmatic reservoir may erupt several times. The reservoir activity lifetime is calculated for 5 different reservoir volumes ranging from 10 6 to 10 11 m 3 and for 6 different depths ranging from 1 to 10 km under the surface. Lifetime increases with the reservoir volume, for both pure and briny cryomagma, ranging from 0.4 years to 10 5 years. In addition, reservoir lifetime is approximately 10 times larger for reservoir at 10 km depth in comparison to 1 km depth.

Reservoir lifetime depends on reservoir depth because of the temperature gradient in the ice crust: reservoirs near the surface are located in a colder environment, which makes them freeze faster.

In the model presented above, we considered that a cryomagma reservoir may remain active as long as some liquid remains in its interior. Nevertheless, as cryomagma freezes, the crystal concentration in the liquid might increase, which increases the effective cryomagma viscosity [START_REF] Roscoe | The viscosity of suspensions of rigid spheres[END_REF]. The Einstein-Roscoe law [START_REF] Roscoe | The viscosity of suspensions of rigid spheres[END_REF] gives the effective viscosity µ ef f of a fluid and crystals mixture as a function of the fluid viscosity µ and the crystal concentration φ assuming spherical crystals: [START_REF] Marsh | On the crystallinity, probability of occurrence, and rheology of lava and magma[END_REF] studied the effect of crystals in terrestrial magmas, and based on the Einstein-Roscoe law, calculated that eruption of lava should stop once 50 to 60% of crystallinity is reached due to a large increase in viscosity. [START_REF] Quick | Cryovolcanic emplacement of domes on Europa[END_REF] proposed to apply this criterion to cryomagma on Europa and consider that cryomagma eruption is not possible for crystal content above ∼55%. Nevertheless, cryomagma on Europa is expected to be a water-based mixture [START_REF] Kargel | Brine volcanism and the interior structures of asteroids and icy satellites[END_REF], with a viscosity lower than terrestrial magmas of several orders of magnitude. To calculate the maximum crystallinity allowing cryomagma to erupt, we use 1987 ;[START_REF] Fagents | Considerations for effusive cryovolcanism on Europa: The post-Galileo perspective[END_REF][START_REF] Quick | Cryovolcanic emplacement of domes on Europa[END_REF]. [START_REF] Allison | Ice-covered water volcanism on Ganymede[END_REF] studied the flow of liquid water on Ganymede's surface and found that a ∼ 0.5 m thick ice crust is enough to prevent the flow from boiling, which allows the underlying liquid to flow onto the surface. [START_REF] Quick | Cryovolcanic emplacement of domes on Europa[END_REF] calculated crust thickness as a function of the time for erupted brines on Europa and found that a 0.5 m thick crust would form in ∼7.5 days. Before the formation of a solid crust on top of the flow, the liquid would boil violently. Moreover, once the crust is formed, it would stabilize the interior portion of the flow [START_REF] Allison | Ice-covered water volcanism on Ganymede[END_REF][START_REF] Quick | Cryovolcanic emplacement of domes on Europa[END_REF], but the expanding edges are supposed to look like a mix of boiling water and ice blocks being pushed by the liquid [START_REF] Allison | Ice-covered water volcanism on Ganymede[END_REF]. It is thus necessary to estimate the vaporized fraction of cryolava during the eruption.

µ ef f = µ (1 -1.35φ) -2.5 (7)
The water fraction being vaporized at the surface during the cryolava flow is a key parameter to link the volume of the smooth features and the initial erupted volume of liquid water, hence, it must be evaluated.

This quantity was roughly estimated by [START_REF] Porco | Cassini observes the active south pole of enceladus[END_REF], who calculated the vaporized fraction of water during the opening of cracks on Enceladus. One should note that this study does not take into account the formation of an ice crust as discussed above, so the amount of material vaporized used in this study can be considered as an upper boundary. [START_REF] Porco | Cassini observes the active south pole of enceladus[END_REF] assumes that the latent heat of fusion L f generated by the freezing part of liquid water is used as latent heat of vaporization L v by the vaporized part of the fluid, so they calculate that a fraction x = L f /(L f + L v ) = 0.13 of liquid is vaporized. This simple calculation from [START_REF] Porco | Cassini observes the active south pole of enceladus[END_REF] gives an idea of the quantity of liquid turned into vapor when it reaches the icy moon's surface, but to have a better knowledge of the transformation occurring in the two-phase region it is necessary to use the phase diagram of water. To assess the water behavior in extreme environments such as Europa's surface, it is relevant to use a temperature-entropy (T-s) phase diagram [START_REF] Lu | Thermodynamics and mass transport in multicomponent, multiphase H2O systems of planetary interest[END_REF] as the flow might be considered as isentropic (adiabatic and reversible) as discussed by [START_REF] Kieffer | Isentropic decompression of fluids from crustal and mantle pressures[END_REF]. An isentropic process is a vertical line on a T-s phase diagram, so it is easy to deduce flow properties based on this representation, and one can directly read the mass ratio of gas/solid on such a diagram. Fig. 10 shows the T-s diagram of water adapted from [START_REF] Lu | Thermodynamics and mass transport in multicomponent, multiphase H2O systems of planetary interest[END_REF]. On this diagram, the blue arrow is the isentropic depressurizing process from liquid water at the triple point (T=273 K, P=612 Pa) to the exit state, i.e. the conditions at Europa's surface (low pressure and low temperature). Figure 10: Temperature-entropy diagram of pure water. S stands for "solid", L for "liquid" and G for "gas". x is the vapor ratio in the solid + gas mixture. The blue arrow shows the path followed by liquid water at 273 K coming from a sub-surface reservoir (T=273 K, P=612 Pa the pressure at the triple point) and flowing on the surface of an icy moon (low pressure and temperature). This process takes place between the x=0.1 and x=0.15 vapor ratio lines. This diagram is adapted from [START_REF] Lu | Thermodynamics and mass transport in multicomponent, multiphase H2O systems of planetary interest[END_REF].

These results show that, without taking into account the formation of an ice crust on top of the flow, the vaporized fraction of water should range between 10 and 15% of the total erupted cryolava. The ice crust formed on top the of the flow plays a role in protecting the well-developed flow from boiling, however, the flow edges are still subject to boiling and freezing counteractions. We thus consider a vaporized fraction of 13±3% after [START_REF] Porco | Cassini observes the active south pole of enceladus[END_REF] and [START_REF] Lu | Thermodynamics and mass transport in multicomponent, multiphase H2O systems of planetary interest[END_REF] to calculate the liquid volume that may be at the origin of the smooth features, keeping in mind that this calculated volume is a lower limit of the erupted volume. Also, the addition of salts or impurities in the cryomagma could slightly modify this result as it lowers the vapor pressure on Europa from ∼600 to 300-500 Pa depending on the salt content [START_REF] Quick | Cryovolcanic emplacement of domes on Europa[END_REF].

Another parameter that could affect the density contrast between the liquid and solid cryomagma phases is the formation of hydrates during the solidification process. In our previous work [START_REF] Lesage | Cryomagma ascent on Europa[END_REF], we considered a briny cryomagma composed of the following mixture: 81 wt% H 2 O + 16 wt% MgSO 4 + 3 wt% Na 2 SO 4 , which is the composition of Europa's ocean and ice predicted by [START_REF] Kargel | Brine volcanism and the interior structures of asteroids and icy satellites[END_REF][START_REF] Mccarthy | Solidification and microstructures of binary ice-I hydrate eutectic aggregates[END_REF] show that hydrates of MgSO 4 and Na 2 SO 4 form for concentrations above respectively 17.3 and Nevertheless, hydrates formation must be taken into in the case of a higher salt concentration.

Results

Measured and erupted volumes

Table 1 summarizes smooth feature volumes measured on the four images shown in Fig. 1. Measured volumes V measured are extracted directly from the DEM (see the method in section 2.3.1 and flowchart in Fig. 2). As detailed previously, the measured volume is not equal to the volume of cryomagma erupted at the surface (V erupted ) during the eruption that created the feature. To obtain the erupted volume, we multiply the measured volume by the α V factor to take into account the ridges on the covered surface and putative subsidence (see section 2.3.3 and Fig. 5). We calculate the results for two extreme values of α V : α V = 0.8, which is the case of shallow subsidence, and α V = 5, which describes the maximum subsidence possibly induced by a liquid subsurface reservoir of 5 km radius according to [START_REF] Manga | Formation of lenticulae on Europa by saucer-shaped sills[END_REF].

We also multiply the volume V measured by a factor 1 /1-x (where x is the cryomagma vaporized fraction) to take into account the vaporization of the erupted liquid. Finally, we multiply V measured by a factor ρs ρ l ( 1130 1180 for briny cryomagma or 920 1000 for pure water) to estimate the liquid volume before expansion due to phase change. This process is summarized in the following equation:

V erupted = α V V measured 1 (1 -x) ρ s ρ l (8)
where V erupted is the fluid volume erupted from the reservoir, V measured is the smooth feature volume measured from the DEM using the simple approach, α V is a volume factor to take into account the underlying terrain (see section 2.3.3), x = 0.13 is the fluid vaporized fraction, ρ l is the density of the liquid cryomagma and ρ s is the density of the corresponding ice.

We propagate the uncertainties from the DEM in order to take into account the two main uncertainty sources, i.e. the smoothness coefficient and the reflectance model (and the associated albedo). Uncertainty calculation using mean deviation are detailed in supplementary materials, section 2. We found an uncertainty of ±15% on the volumes measured from the DEM. Moreover, we added a ±3% uncertainty on the calculation of the erupted volume due to the uncertainty on the vaporized fraction (see section 3.3). The DEM of each smooth feature used to calculate the volumes are showed in supplementary materials, section 3. Finally, we obtain the results given in Table 1.

Image

Measured volume (m 3 ) Erupted volume for pure water (m 3 ) (see Eq. 8) α V = 0.8 α V = 5 5452r

(5.7±0.9)×10 7 (4.9±0.9)×10 7 (3.0±0.5)×10 8 0713r (6.6±1.0)×10 (1.5±0.3)×10 9

Table 1: Summary of the measured volumes V measured using the DEMs, and the corresponding erupted volumes V erupted calculated with Eq. ( 8) for 5 m subsidence (α V = 0.8) and 40 m subsidence (α V = 5).

Volume of cryomagmatic reservoirs

We previously obtained cryolava volume necessary to explain the emplacement of the four smooth features from Fig. 1. We can now deduce the reservoir volume required to generate this cryolava amount from our eruption model.

In [START_REF] Lesage | Cryomagma ascent on Europa[END_REF], we obtained the volume of cryolava erupted at the end of a single eruptive event as a function of reservoir volume and depth, for two different cryomagma compositions: pure water and a briny mixture of water and salts: 81 wt% H 2 O + 16 wt% MgSO 4 + 3 wt% Na 2 SO 4 [START_REF] Kargel | Brine volcanism and the interior structures of asteroids and icy satellites[END_REF]. These volumes are shown in Fig. 11 for a reservoir depth ranging from 1 to 10 km and a reservoir volume ranging from 10 8 to 10 12 m 3 , which corresponds to reservoir radius between ∼0.3 and 6.2 km. In Fig. 11, we compare these results with the erupted volumes measured from the four DEMs (uncertainties on the erupted volumes are taken into account). One can see in Fig. 11 that the eruption of a 3×10 9 to 10 11 m 3 reservoir (0.9 to 2.9 km radius) is necessary to explain the formation of smooth features from a single cryovolcanic event for shallow terrain subsidence (α V = 0.8). For deeper subsidence (α V = 5), a 2×10 10 to 10 12 m 3 reservoir (1.7 to 6.2 km radius) is required. The range of erupted volumes is nearly identical for these two compositions as the ρice ρ liq factor does not differ significantly between pure water and a water-based mixture. The erupted volume necessary to produce the observed features depends mostly on the α V ratio. One should note that larger reservoirs are expected to create deeper subsidence of the surface [START_REF] Manga | Formation of lenticulae on Europa by saucer-shaped sills[END_REF], so it is worth it to consider both solutions.

In addition, we also consider in this study the case of several eruptive cycles instead of one single eruption V e (m³)

Figure 11: Volume erupted at the surface during one cryovolcanic eruption after the model from [START_REF] Lesage | Cryomagma ascent on Europa[END_REF]. The smooth feature volume measured on the four analyzed images (see Fig. 1) corresponds to reservoir volumes ranging in between the two dashed lines for α V = 0.8 and between the two solid lines for α V = 5 (see Table 1).

Pure water

Briny cryomagma α V = 0.8

Min. reservoir size (image 5452r) 4×10 8 (R = 0.5 km) 10 9 (R = 0.6 km) Max. reservoir size (image 9352r) 4 ×10 9 (R = 1 km) 10 10 (R = 1.3 km)

α V = 5
Min. reservoir size (image 5452r) 2×10 9 (R = 0.8 km) 6×10 9 (R = 1.2 km) Max. reservoir size (image 9352r) 2×10 10 (R = 1.6 km) 6 ×10 10 (R = 2.4 km)

Table 2: Size of the reservoir required to erupt cryolava amount necessary to generate smooth features observed depending on cryomagma composition and α V ratio. We give these results for the smallest and largest features, respectively from image 5452r and 9352r. We use the following composition for the briny cryomagma: 81 wt% H 2 O + 16 wt% MgSO 4 + 3 wt% Na 2 SO 4 [START_REF] Kargel | Brine volcanism and the interior structures of asteroids and icy satellites[END_REF].

(see section 3.2). In this case, the same reservoir will produce several eruptions during its whole lifetime.

According to Fig. 8, this kind of reservoir may erupt 5×10 6 to 7×10 9 m 3 of cryolava if we consider that the reservoir is active until 100% of the cryomagma is frozen. One should note that this amount may be overestimated as eruption may stop when 55% of the reservoir is frozen [START_REF] Quick | Cryovolcanic emplacement of domes on Europa[END_REF]. This could divide cryolava erupted volume by a fator ∼2, but is not taken into account in the following results as cryolava without crystal is in better agreement with the smooth feature morphologies. We summarize in Table 2 the reservoir volumes required to obtain the minimum and maximum erupted volumes given in Table 1 (respectively for image 5452r and 9352r), taking uncertainties into account.

Discussion and conclusions

We identified four images from Galileo SSI data presenting smooth features that may have formed during one or several eruptions of cryolava at Europa's surface. We produced DEMs of those smooth features using the Shape from Shading tool of the AMES Stereo Pipeline, and we measured their volumes and associated uncertainties. The major uncertainties on the DEM generation and volume measurement come from the user parameters of SfS, but their effect is relatively low for small scale features. We estimated the uncertainty on the volumes at around ±15% (see supplementary materials, section 2). Volume measurements of the four selected smooth features gave results ranging from 5×10 7 to 5×10 8 m 3 . Shape from Shading tool allowed us to generate DEMs from single images, with high precision at low scales [START_REF] Nimmo | Stereo and photoclinometric comparisons and topographic roughness of Europa[END_REF]. Nevertheless, SfS presents its own limitations. First of all, it does not manage albedo or photometric heterogeneities in a single image. In fact, the reflectance model chosen by the user is applied to the whole image, but the surface properties can differ at regional or local scales, as shown by [START_REF] Jiang | Fusion of photogrammetric and photoclinometric information for high-resolution DEMs from Mars in-orbit imagery[END_REF] and [START_REF] Belgacem | Regional study of Europa's photometry[END_REF]. We showed that the reflectance model and the associated chosen coefficients (such as albedo) are not expected to create uncertainties on volume higher than ±10% for small scale features. Nevertheless, as discussed by [START_REF] Jiang | Fusion of photogrammetric and photoclinometric information for high-resolution DEMs from Mars in-orbit imagery[END_REF] for Martian images, too many heterogeneous surface properties in one image can lead to non-convergence of the algorithm. From the images we selected, four were converging (images presented in Fig. 1), but it was not the case of image 8613r, which contains a smooth material darker than the surrounding terrain, maybe due to an albedo heterogeneity. We could not provide a DEM for this image. It could be interesting to investigate the possibility of considering several zones with different surface properties in an image: this could solve the problem of DEM generation of image 8613r and give information on this smooth feature.

Geomorphological interpretations of DEMs are consistent with smooth feature formation by the flow of a fluid on the surface: smooth features have a very thin appearance, are constrained by the surrounding ridges and occupy topographic lows. Moreover, DEMs show that some double ridges have an elevation that decreases in the direction of smooth feature centers. We suggest that this could result from terrain subsidence beneath the smooth features. Thermal erosion could also participate to this particular morphology, but it is not required to explain the formation of smooth features.

To link the volume of the smooth features measured on the DEM (using a simple approach) with the actual volume of cryolava erupted during their emplacement, we propose to take into account a range of possible subsidence depths, as modeled by [START_REF] Manga | Formation of lenticulae on Europa by saucer-shaped sills[END_REF]. We studied the two extreme cases of shallow subsidence of 5 m and maximum subsidence of 40 m. In the first case, approximately 4×10 7 to 4.4×10 8 m 3 of cryolava are required to erupt in order to emplace the observed smooth features. In the second case with deeper subsidence, approximately 2.5×10 8 to 2.7×10 9 m 3 of cryolava are required. These volumes depend slightly on the cryomagma composition, but this does not change significantly the results. These extreme values may be better constrained knowing subsidence height under the smooth features, but unfortunately, this is not possible to infer with current data. Larger reservoirs are expected to create deeper surface subsidence [START_REF] Manga | Formation of lenticulae on Europa by saucer-shaped sills[END_REF], so it is worth considering both of these solutions.

Using our previous model [START_REF] Lesage | Cryomagma ascent on Europa[END_REF], we can predict the cryolava volume erupted at the surface at the end of a cryovolcanic eruption as a function of reservoir volume and depth. These volumes are compared to the volumes measured on the four Galileo images to constrain cryomagmatic reservoir volumes. We found that a 3×10 9 to 10 11 m 3 (0.9 to 2.9 km radius) cryomagma reservoir is required to explain the emplacement plains are the result of one or several eruptions. Higher resolution images expected from future missions will provide better resolved DEMs of the surface and more information on smooth plain morphology. This would help to better constrain putative cryomagma reservoir dimensions.

Figure 1 :

 1 Figure 1: Images selected to generate DEMs. (a) Image 5452r (-1°, 340°), resolution: 27 m/px; (b) Image 0713r (-79°, 124°), resolution: 57 m/px; (c) Image 0739r (-81°, 132°), resolution: 57 m/px; (d) Image 9352r (-28°, 218°), resolution: 60 m/px. All these images are from the Galileo SSI data. Scales are indicative as these images are not projected. Orange arrows show the sunlight direction. The blue arrow points out a lobate feature located on top of the smooth plain of image 9352r.

  Fig. 5 b and c on the left, we show a reference situation where cryolava flows between the ridges without any subsidence. We calculate the cross-sectional area of the smooth feature measured with the simple approach A measured , which is the mean height of the smooth feature top (elevation ∼0 m) above the mean elevation of valley bottoms (∼-7 m) multiplied by the feature's width. This cross-sectional area is filled in grey in Fig. 5 b and c on the right.

Figure 3 :Figure 4 :

 34 Figure 3: (a) Image 5452r from Galileo SSI. Sun light comes from the East (right on the image).(b) The DEM of image 5452r produced with SfS. The blue line indicates the limit of the flow-like feature. Blue arrows show the lobate zones interpreted as filling the valleys between the ridges. The two profiles numbered (c) 1 and (d) 2 are respectively perpendicular and parallel to the surrounding ridges. Schematic representations of feature sections along these two profiles are shown in Fig. 4. Black crosses indicate the chosen points used to calculate the reference level of the pre-existing terrain beneath the smooth feature (see Fig. 4).

Figure 5 :

 5 Figure 5: Illustration of the calculation of the α V = V cryolava

Figure 6 :

 6 Figure6: Schematic representation of a cryomagma reservoir of volume V and radius R = R i , located at depth H below the surface, with liquid cryomagma in white and frozen cryomagma in dashed grey. (a) The reservoir is filled with pure or briny liquid water at isostatic pressure P 0 . (b) An initial liquid volume V i freezes up to a volume V f of ice, inducing an overpressure ∆P in the reservoir. (c) When the pressure reaches a critical value ∆Pc, the wall fractures and the pressurized liquid rises to the surface through an H long fracture. (d) Once a certain amount of liquid has erupted at the surface and the overpressure in the reservoir is released, the eruption ends. The liquid in the fracture freezes, which seals the reservoir. Finally, the reservoir returns to initial condition similar to situation (a) but with a smaller radius R i+n . Freezing of the reservoir continues such that these 4 steps repeat and form an eruptive cycle, leading to several eruptions during the reservoir's lifetime.

Figure 7 :

 7 Figure7: (a) Cumulative volume erupted at the surface during the active lifetime of a reservoir. Each step corresponds to a complete eruptive cycle: the x-axis represents the freezing time (the eruption time is negligible compared to the freezing time), and the y-axis stands for the erupted volume. This result is obtained with a 10 9 m 3 reservoir located 4 km beneath the surface and filled with briny cryomagma. (b) Volume of cryolava erupted during each eruption as a function of the eruption number. Note that y-axis is logarithmic here. (c) Freezing time before each eruption as a function of the eruption number (this time is non-cumulative).

Figure 8 :

 8 Figure 8: Total activity time of a cryomagmatic reservoir producing cyclic eruptions as a function of the reservoir volume and depth. Two cryomagma compositions are tested: (a) pure water and (b) a briny mixture of 81 wt% H 2 O + 16 wt% MgSO 4 + 3 wt% Na 2 SO 4 . The reservoir is active until all the cryomagma is frozen.

Figure 9 :

 9 Figure 9: Total cryolava volume erupted from a reservoir producing cyclic eruptions as a function of reservoir volume. Two cryomagma compositions are used: (a) pure water and (b) a briny mixture of 81 wt% H 2 O + 16 wt% MgSO 4 + 3 wt% Na 2 SO 4 .The reservoir is active until the volume of ice in it reaches the cryomagma crystallinity φ times the reservoir volume V (or simply the reservoir volume V in the case of φ =0).

  This process takes place between the x=0.1 and x=0.15 vapor fraction lines on Fig.10. It indicates that at the end of the process, only 10 to 15% of the liquid water erupted at the surface is vaporized; the major fraction of water freezes to solid.

  Flowchart of the methods used in this study to generate and interpret DEMs. See supplementary materials for details on the pre-processing steps using ISIS. DEM generation using the AMES StereoPipeline is described in section 2.2. Post-processing with QGis and volume calculation are detailed in section 2.3.1.

	INPUT: Raw Galileo SSI image + label text le	Calibrate, correct and project the image	Generate a DEM with SfS	Draw limits of the smooth feature Interpolate terrain elevation the pre-existing	Volume x area calculation: smooth feature -pre-existing terrain	OUTPUT: Smooth feature volume
		ISIS	AMES StereoPipeline		QGis
	Figure 2:

wt%. Thus, in the cases considered in this work, hydrates should not form in the freezing cryomagma.

Acknowledgments

We acknowledge support from the "Institut National des Sciences de l'Univers" (INSU), the "Centre National de la Recherche Scientifique" (CNRS) and "Centre National d'Etudes Spatiales" (CNES) through the "Programme National de Planétologie". We thank Benoit Jabaud for his analysis of the Galileo/SSI images. We also thank Baptiste Journaux for the interesting discussion we had. We gratefully acknowledge the two anonymous reviewers for their comments and suggestions that greatly improved this manuscript.

Eq. ( 7) and effective viscosity measurements of cryovolcanically emplaced domes done by [START_REF] Quick | Cryovolcanic emplacement of domes on Europa[END_REF].

They deduced from those measurements that cryolava should have an effective viscosity µ ef f ranging from 1 to 10 4 Pa.s during the eruption to create dome features. Based on this result, and using the liquid water viscosity µ 10 -3 Pa.s as the cryolava liquid phase viscosity, we can deduce an upper limit of crystallinity of φ 0.74 for erupted cryomagma, assuming that the crystals were the only effect acting to increase the cryolava viscosity at the surface.

As eruptions are expected to stop when cryomagma reaches a maximum crystallinity, this threshold can be considered as a stop condition for our cyclic eruption model. Fig. 9 shows the total volume of cryolava erupted during the reservoir lifetime for three different maximum cryolava crystallinity: φ=0, 0.55 and 0.74.

φ=0 corresponds to cryolava that does not contain crystals. This could occur if crystals are totally separated from cryomagma, if they remain on the reservoir wall for example. In this case, crystals do not ascend with the fluid and cryolava keeps a very low viscosity, close from the pure water viscosity. Because of the very smooth and thin appearance of the features studied here, φ could potentially be close to 0. To model this extreme case, we continue the reservoir freezing until all the cryomagma is turned to ice. φ=0.55 is the maximal cryolava crystallinity suggested by [START_REF] Quick | Cryovolcanic emplacement of domes on Europa[END_REF] according to the study of [START_REF] Marsh | On the crystallinity, probability of occurrence, and rheology of lava and magma[END_REF]. To model this maximum crystal concentration in cryolava, we stop the reservoir activity when 55% of reservoir volume is solid. Finally, φ=0.74 corresponds to the crystallinity calculated from Eq. ( 7) and the maximum cryolava viscosity obtained by [START_REF] Quick | Cryovolcanic emplacement of domes on Europa[END_REF]. We model this threshold by stopping the reservoir activity when 74% of reservoir volume is solid. For reservoir volumes ranging from 10 7 to 10 11 , we finally obtain erupted volumes ranging from 5×10 5 to 7×10 9 m 3 for pure water (Fig. ( 9)a), and from 2×10 5 to 3×10 9 m 3 for the briny cryomagma. To obtain the following results, we use φ=0 as it seems to be in better agreement with the very thin topography of smooth plains studied here. Nevertheless, one should keep in mind that the erupted volumes may be divided by at most a factor two if a more viscous cryomagma is involved.

Vaporized fraction of water

The putative flow of water-based liquid on Europa takes place in a low pressure and temperature environment. The pressure at Europa's surface is near 10 -6 Pa [START_REF] Hall | Detection of an oxygen atmosphere on Jupiter's moon Europa[END_REF] and the mean temperature is around 110 K [START_REF] Spencer | Temperatures on Europa from Galileo photopolarimeter-radiometer: Nighttime thermal anomalies[END_REF], so that the liquid water erupting at the surface is subjected to the counteractions of freezing and boiling. Because of the 160 K difference between the liquid and the environment, several authors previously proposed that an ice crust forms rapidly on top of the flow (Allison and Clifford, of these smooth features from a single eruption in the case of shallow subsidence (Fig. 11, dashed lines). For deep subsidence of 40 m, the required reservoir volumes are 10 times higher (reservoirs up to 6 km in radius, see Fig. 11, solid lines). In the case of cyclic eruptions from the same reservoir, for a pure water reservoir and at the end of its activity lifetime, a 4×10 8 to 4×10 9 m 3 (0.5 to 1 km radius) liquid reservoir is needed in case of shallow subsidence, and a 2×10 9 to 2×10 10 m 3 (0.8 to 1.6 km radius) one in case of a 40 m subsidence.

The cryomagma composition slightly changes these results: two to three times greater reservoir volume is required for a briny cryomagma compared to pure water (up to 2.4 km radius). The total lifetime of such reservoirs ranges from 5×10 3 to 10 5 years if we consider that eruptions stop when 100% of the cryomagma freezes. Nevertheless, as suggested by [START_REF] Quick | Cryovolcanic emplacement of domes on Europa[END_REF], eruptions could become impossible when 55% of crystallinity is reached because cryomagma viscosity at that point becomes too high. In this case, a two times bigger reservoir is needed to explain the smooth features observed.

The recent detection of chlorides such as NaCl [START_REF] Trumbo | Sodium chloride on the surface of Europa[END_REF] or Mg-bearing chlorinated species [START_REF] Ligier | VLT/SINFONI observations of Europa: new insights into the surface composition[END_REF]) on Europa's surface may indicate that the cryomagma composition could be different from the one used here (81 wt% H 2 O + 16 wt% MgSO 4 + 3 wt% Na 2 SO 4 , see [START_REF] Kargel | Brine volcanism and the interior structures of asteroids and icy satellites[END_REF]. Two main quantities might be impacted by the cryomagma composition: (i) the freezing time-scale of the reservoir and (ii) the cryolava erupted volume. The freezing time-scale is a function of the freezing temperature of the solution, which is slightly lower for chloride aqueous solutions than for sulfate ones [START_REF] Quick | Heat transfer of ascending cryomagma on Europa[END_REF]. Thus, the freezing time scale must be slightly larger for chloride solutions. On the other hand, the erupted volume of cryomagma depends on the density contrast between the cryomagma solution and the corresponding ice.

Here we tested a sulfate-based cryomagma and the extreme case of pure water. The erupted volumes obtained with these two compositions are similar for a single eruption and differ by at most a factor two for several eruptions during the whole reservoir lifetime. Hence, the cryomagma composition is not expected to modify the results order of magnitude, at least for a reasonable salt content. Moreover, Eq. ( 6) provides a very simple way to predict the cryomagma volume erupted from a subsurface reservoir depending on cryomagma and ice densities. It is then possible to adapt this model to any cryomagma composition.

We demonstrated in this study that cryomagmatic reservoirs of ∼10 7 to 10 12 m 3 located a few kilometers under Europa's surface may possibly be at the origin of smooth features seen on the Galileo images 5452r, 0713r, 0739r and 9352r. This information could help the two upcoming missions JUICE (ESA) and Europa Clipper (NASA) to target interesting locations to search for biosignatures. In order to better constrain the characteristic size and depth of the source reservoirs, it would be necessary to determine whether the smooth