Thermo-mechanical modeling of pancakelike domes on Venus

Calusi, Benedetta; Farina, Angiolo; Fusi, Lorenzo; Rosso, Fabio

doi:10.1063/5.0209674

In this paper, we present a mathematical model aimed at describing both the effusive and relaxing phase of pancakelike lava domes on the Venus surface. Our model moves from the recent paper by Quick et al. [“New approaches to inferences for steep-sided domes on Venus,” J. Volcanol. Geotherm. Res. 319, 93–105 (2016)] but generalizes it under several respects. Indeed, we consider a temperature field, playing a fundamental role in the flow evolution, whose dynamics is governed by the heat equation. In particular, we suggest that the main mechanism that drives cooling is radiation at the dome surface. We obtain a generalized form of the equation describing the dome shape, where the dependence of viscosity on temperature is taken into account. Still following Quick et al. [“New approaches to inferences for steep-sided domes on Venus,” J. Volcanol. Geothermal Res. 319, 93–105 (2016)], we distinguish an isothermal relaxing phase preceded by a non-isothermal (cooling) effusive phase, but the fluid mechanical model, developed in an axisymmetric thin-layer approximation, takes into account both shear thinning and thermal effects. In both cases (relaxing and effusive phase), we show the existence of self-similar solutions. In particular, this allows to obtain a likely scenario of the volumetric flow rate which originated this kind of domes. Indeed, the model predicts a time varying discharge, which is maximum at the beginning of the formation process and decreases until vanishing when the effusive phase is over. The model, in addition to fitting well the dome shape, suggests a possible forming scenario, which may help the largely debated questions about the emplacement and lava composition of these domes.

I. INTRODUCTION: GENERAL CONSIDERATIONS ABOUT TERRESTRIAL AND EXTRA-TERRESTRIAL LAVA DOMES

The mathematical modeling of volcanic domes is an extremely challenging subject due to the great variety of morphological shapes observed at the end of the effusion phase. Indeed, the erupted lava is a very complex fluid, being its viscosity and density strictly depending on the chemical composition, water content, and temperature, which in turn affects the degree of crystallization.^2–5 Dome shapes may vary enormously, although the common origin (eruption from a vent) leads, usually, to an axisymmetric structure, whose aspect ratio can range from 0.1 to 0.05 or even less. The most massive shield volcanoes on Earth are of this type, with heights of order thousands of meters and diameters of order tens of kilometers.

The leading physical mechanisms that guide dome formation are gravity and temperature, and most of the mathematical research devoted to presenting a model for dome evolution is based on the equations of fluid mechanics in the framework of the lubrication approximation. A pioneer of this line of investigation applied to dome growth was Huppert in a series of papers^6–8 although other authors before him dealt with the dynamics of gravity currents.⁹ Many other investigators have approached this topic since then.^10–16 These contributions consider both isothermal and non-isothermal conditions and Newtonian as well as non-Newtonian rheology.

Some domes on Venus (known in the literature with the name of pancake domes) received special attention in Refs. 1, 15, and 17–27 because of certain special characteristics of these volcanic structures. The aim of the present paper is to make a contribution, namely, in this latter case.

Before going into details, we remind that volcanic domes are observable in all planets and satellites of the solar system. While on Earth most physical parameters are more or less easily accessible, the situation is totally different for non-terrestrial domes. In the latter case, the best we can know is the present shape of the dome, and in some of these cases, even the shape is only partially known with sufficient precision. Domes on Venus are just like that; indeed, its very dense atmosphere precludes direct visual observations and all the data available come from radar measurements by dedicated spacecraft. In particular, for any Venusian dome, we have no indication on its age, that is on the duration of the growing or effusion phase [corresponding to an outflow at the vent $Q (t) \neq 0$ ] and that of the relaxing phase [i.e., Q(t) = 0], when the cooling process competes with gravity and ambient pressure to produce the observed shape. The characteristic times of the two phases are, on the Earth, very different. However, as we shall see in the sequel, the situation does not seem to be the same on Venus.

The case of Venus should be viewed as a case study. Indeed, as emphasized in Ref. 28, the formation of lava domes on Earth requires viscous magma, and, as a result, dome-building volcanoes are primarily found at plate boundaries.²⁹ Although the formative processes are sometimes unclear.^24,27,30 volcanic domes observed on other planetary bodies suggest a similar evolution, including the existence of plate tectonics.³¹ In particular, volcanism also suggests the process of volatile degassing, a phenomenon that has been linked to the development of the nitrogen-enriched atmosphere required to create a habitable planet.³² Some authors^18,19 suggested, in addition, a possible analogy between Hawaiian submarine domes and Venus pancake domes. Indeed, the control of volatile exsolution by pressure on Venus and Earth seafloor can cause lavas to have similar viscosities and densities, although the latter will be counteracted by high buoyancy underwater. In any event, analogous effects of the Venusian and seafloor environments alone are probably not sufficient to produce similar volcanoes. Other points, still under debate, concern both which kind of lava, basaltic^17,18 or rhyolitic^27,30 allows the observed shapes and whether the observed morphology is closer to a single continuous emplacement of magma or to multiple events.

All the indications of geologists and planetary experts mentioned up until now, appear, by far, neither convincing nor conclusive. Indeed, it is quite evident that, lacking precise chemical-geophysical data, the best one could do is to construct a “morphology-oriented” mathematical model, sufficiently flexible in the key parameters to be able to justify the observed shapes, which is exactly the aim of the present paper.

In the early 1990s, the Magellan spacecraft (also known as the Venus Radar Mapper) used SAR (Synthetic Aperture Radar) imaging to produce detailed topography data for most of the planet's surface. It turned out that about 90% of Venus's surface was covered by lava flows and broad shield volcanoes. On Venus, large shield volcanoes are an impressive 700 km wide at the base but are only about 5.5 km in height (aspect ratio about 0.007), thus being several times as wide as those on Earth but with a much gentler slope (for example, the Vesuvio aspect ratio is about 0.025). A possible reason for this could be that lava's mobility might be enhanced by the planet's average surface temperature (⁠ $\approx 700$ K). Indeed, the high temperature contributes to maintaining more of the volcanic material in its fluid state, preventing solidification and favoring the length of the relaxing time after the end of the volume-growing phase.

The Magellan exploration mission also allowed to identify the previously mentioned special class of volcanoes called pancake domes since then. They appear to be steep-sided and flat-topped, both features of enigmatic origin. Pavri et al.²⁶ reported 145 domes of this type on the planet, with diameters ranging from 7 to 94 km and heights over 4 km for up to 60 of them. These numbers have been updated by Stofan et al.¹⁷ These domes are thought to form by the extrusion of viscous lava. However, the composition of these domes is still puzzling. Indeed, although their smooth surfaces are consistent with basaltic compositions,^17,19,20 their steep sides suggest that they are composed of viscous lava, consistent with high-silica, rhyolitic compositions.^26,27,33 In any event, the concentric and radial fracture pattern shown by many domes is consistent with the stretching of dome surfaces during their formation as lava wells up from interior vents and spreads laterally on the surface.^17,33

The first data provided by the Magellan mission were SAR data (see Fig. 1). Only some years later, Gleason et al.²⁵ applied stereo radargrammetric techniques to Magellan SAR data (⁠ $\approx 100$ m resolution) to generate DEMs (Digital Elevation Models), $\approx 1$ km horizontal resolution, for 23 steep-sided domes. The heights of the domes range from 390 to 1590 m, aspect ratios (height/diameter) range from 0.01 to 0.07, and volumes range from 73 to 1790 ${km}^{3}$ ⁠. Figure 2 shows one of these pancakelike domes in the Rusalka Planitia on Venus. Figure 3 reproduces the SAR altimetry profiles obtained through two different orbits of the Magellan mission for the dome shown in Fig. 2,²⁷ while Fig. 4 shows two different DEM profiles of the same dome.

FIG. 1.

View large Download slide

High-resolution Magellan SAR image of a group of three “pancakelike” domes in Carmenta Farra Planitia on Venus (NASA, public domain).

FIG. 2.

Left panel: Magellan SAR image of a volcanic “pancakelike” dome, located 3 o S, 151 o E on Venus. Right panel: DEM model. This dome is ≈ 30 − 35 km in diameter and about 1370 − 1560 m tall. The two transects refer to the N-S and SW-NE section from which the two profiles shown in Fig. 4 are derived.25 Reproduced with permission from Gleason et al., “Analysis of Venusian steep-sided domes utilizing stereo-derived topography,” J. Geophys. Res.: Planets 115, E06004 (2010). Copyright 2010 John Wiley and Sons, License No. 5770120388193.

View large Download slide

Left panel: Magellan SAR image of a volcanic “pancakelike” dome, located $3^{o}$ S, $151^{o}$ E on Venus. Right panel: DEM model. This dome is $\approx 30 - 35 km$ in diameter and about $1370 - 1560 m$ tall. The two transects refer to the N-S and SW-NE section from which the two profiles shown in Fig. 4 are derived.²⁵ Reproduced with permission from Gleason et al., “Analysis of Venusian steep-sided domes utilizing stereo-derived topography,” J. Geophys. Res.: Planets 115, E06004 (2010). Copyright 2010 John Wiley and Sons, License No. 5770120388193.

FIG. 3.

View large Download slide

Altimetry profile from Magellan SAR data of the dome in Fig. 2 corresponding two different orbits. The plot-markers indicate altimetry data from orbit Nos. 1277 and 3067.²⁷ These data have an accuracy of 100–150 m as described in Ref. 34.

FIG. 4.

Digital elevation profiles along the N-S and SW-NE transection of the dome shown in Fig. 2.25 Reproduced with permission from Gleason et al., “Analysis of Venusian steep-sided domes utilizing stereo-derived topography,” J. Geophys. Res.: Planets 115, E06004 (2010). Copyright 2010 John Wiley and Sons, License No. 5770120388193.

View large Download slide

Digital elevation profiles along the N-S and SW-NE transection of the dome shown in Fig. 2.²⁵ Reproduced with permission from Gleason et al., “Analysis of Venusian steep-sided domes utilizing stereo-derived topography,” J. Geophys. Res.: Planets 115, E06004 (2010). Copyright 2010 John Wiley and Sons, License No. 5770120388193.

The aim of the present paper is to generalize the model proposed by Quick et al.,¹ but allowing a more consistent way to take into account the influence of temperature on the dome evolution. First, lava is a strongly non-Newtonian fluid with highly temperature dependent material properties.^2–5,35 The fluid rheology is complicated yet further by the phase change of solidification. For lava, the dome is unlikely to be vertically isothermal because the diffusivity is small. Instead, the cooling of the surface creates a thermal boundary layer that slowly migrates into the fluid interior. However, in the special limit of the thin layer approximation (the domes we consider have a low aspect ratio), the evolution is rather slow, and thermal diffusion acts sufficiently quickly to make the fluid isothermal in vertical cross sections. For this reason, we consider an idealized problem: the emplacement of a cooling non-Newtonian fluid on a plane, modeling its “solidification” as the increasing of viscosity as temperature decreases. We do this both in the interest of simplicity, and because we would first like to understand the dynamics of purely fluid Venusian domes. We allow a more general rheological model with respect to that considered by Quick et al.¹ Indeed, we follow the approach by Rajagopal et al.,³⁶ adopting a shear-thinning power-law rheology within the lubrication approximation framework. Nevertheless, our approach differs from that Rajagopal et al.,³⁶ since they assume viscosity to be pressure-dependent, while, in the present context, a temperature dependence is more appropriate. We point out that having considered a non-Newtonian constitutive model represents a step forward with respect to the theory illustrated in Ref. 1. Indeed, the non-Newtonian rheology, introducing a further parameter besides viscosity, allows to fit different profiles of the same dome (e.g., SW-NE and N-S) with the same set of parameters. On the contrary, this result is not always possible if a purely Newtonian model is considered (as it will be emphasized in Fig. 6).

As in Ref. 1, we look for similarity solutions of the nonlinear partial differential equation for the dome shape, in both isothermal and non-isothermal conditions. In addition to allowing a more general rheology, the main difference with respect to Ref. 1 is that we couple the momentum equations with the heat equation (although in a simplified form). More precisely, we assign the dependence of viscosity on temperature via an Arrhenius-type law, but, unlike in Ref. 1, the change of temperature with time is governed by the heat equation. Indeed, through a careful analysis of the timescales involved in the temperature evolution equation, we suggest that the main mechanism that drives cooling is radiation at the dome surface. From the point of view of including the thermal effects during the dome evolution, we mention the work by Balmforth and Craster.¹¹ However, their paper does not consider the possibility of the existence of similarity solutions for the differential equation which governs the evolution of the dome surface. Furthermore, their analysis is more addressed to compare the dome evolution as predicted by their model to laboratory experiments with experimental slurry domes, cooled by conduction and convection in overlying water. Our goal is rather different: just as Quick et al.,¹ the only comparison that can be made is with the present shape of the dome as it appears from radar images. Moreover, as in Ref. 1, we are more interested in linking the observed shape to the main phases (effusive and relaxing) of the dome evolution, trying to give some estimates of the timescales involved. Furthermore, our model suggests a possible flow rate at the vent, which, differently from Ref. 1, decreases monotonically with time. Indeed, by using the mass conservation law, we identify a specific decaying law for Q(t), characterized by its vanishing at the end of the effusion phase.

The plan of the paper is as follows: in Sec. II, we present the general rheological model of a power-law fluid flowing on a plane in Cartesian coordinates. Then, we introduce some scale factors and rewrite the main equations in dimensionless form. Following the procedure illustrated in Ref. 36 we obtain a generalized form of the nonlinear partial differential equation for the evolution of the dome surface boundary [see Eq. (19)], where the dependence of viscosity on temperature is taken into account. In Sec. III, we discuss the relaxing phase in the framework of lubrication theory and in the axisymmetric geometry. First, we assume that relaxing goes isothermally. This leads to a simpler form of the equation describing the evolution of the dome top surface h(r, t), which coherently reduces to its classical form in the case of a Newtonian rheology.⁶ Subsequently, we adhere to the same guidelines as Quick et al.¹ to distinguish the forming phase from the relaxing one. This involves to set the initial instant of the relaxing phase to be a new parameter τ, which can be identified by comparing the present dome shape with the one obtained through simulations. Our solution is fully explicit and provides h(r, t), the terminal front R(t) of the dome and the relation between the dimensionless volume of the dome with respect to all the relevant parameters. We remark that, in the case of a Newtonian rheology, our analysis reproduces exactly the same solutions obtained in Ref. 1. Section IV is devoted to comparing the evolution of the dome as provided by the model, with the radar profile of the dome shown in Fig. 2. The aim is to identify the main dimensionless parameters of the model, as the duration of both the effusive and relaxation phases. We rescale the parameters to their dimensional values by using the literature data, and we find a satisfactory agreement with the values reported in Refs. 25 and 27. One important result is that the procedure provides an estimate of the duration of the effusive phase with respect to the relaxing phase. To get more insights on these time intervals, a further discussion on the choice of the right timescale is needed. This topic is illustrated in Sec. V. In particular, we estimate the values of the Reynolds and Froude numbers to check that all simplifying assumptions made to develop our model are effectively largely satisfied. Section VI is devoted to couple the flow with the thermal field. By a careful identification of the characteristic times involved by diffusion, convection, and radiation, we obtain a simplified version of the heat equation, which should capture the most important effects of the thermal field during the effusive phase. In Sec. VII, we consider the case of spatially uniform temperature $θ (t)$ ⁠. In this case, although explicit solutions for h(r, t) are no longer available, we have been able to identify a similarity solution. Indeed, the equation for h(r, t) can be transformed into a suitable equation for $h (r, φ)$ where $φ = φ (θ)$ is an auxiliary function. This latter has to solve a particular differential equation, while $θ = θ (t)$ is obtained from the heat equation when radiation is the only cooling mechanism. In Sec. VIII, we solve the equations for the coupled problem, integrating the evolution equation backward in time, using as initial data those obtained from the analysis of the purely relaxing phase. Section IX is devoted to some conclusions.

II. FORMULATION OF THE GENERAL MODEL

Let us consider an incompressible continuum flowing on a horizontal plane where

O, x, y, z

is the Cartesian frame of reference, and z is the coordinate orthogonal to the plane. The velocity field is

v = u (x, y, z, t) e_{x} + v (x, y, z, t) e_{y} + w (x, y, z, t) e_{z},

and the domain occupied by the flowing fluid is given by

Ω = {(x, y, z) \in ℝ^{3} | (x, t) \in D_{t}, 0 \leq z \leq h (x, y, t)},

(1)

where

D_{t} \subset ℝ^{2}

is portion (evolving in time) of the (x, y) plane “wet” by the material and

z = h (x, y, t)

is the top surface. In particular,

{h (x, y, t) |}_{(x, y) \in \partial D_{t}} = 0

⁠.

The continuity equation and the momentum equation are

{\begin{matrix} \nabla \cdot v = 0, \\ ρ \frac{D v}{D t} = - ρ g e_{z} - \nabla p + \nabla \cdot S, \end{matrix}

(2)

where

\frac{D}{D t}

denotes the material derivative, ρ is the fluid density, g is the gravity acceleration,

p

denotes the pressure, and

S

represents the deviatoric part of the Cauchy stress tensor, whose constitutive equation is

S = S (θ, | | D | |) D,

(3)

where

S (θ, | | D | |)

is the apparent viscosity of the fluid, θ its absolute temperature and

D = [\nabla v + \nabla v^{T}] / 2

the strain rate tensor. In particular, we assume

S (θ, | | D | |) = 2 μ_{c} η (θ - θ_{amb}) {‖ D ‖}^{(1 - λ) / λ},

(4)

with

θ_{amb}

external temperature,

η (θ - θ_{amb}) > 0

a sufficiently regular, dimensionless and monotone decreasing function of θ, μ_c a physical parameter with dimension

Pa \cdot s^{1 / λ}

⁠, and where λ is a strictly positive constant. The constitutive model allows for these behaviors: shear-thinning if

λ > 1

⁠, shear-thickening if

0 < λ < 1

⁠, Newtonian if λ = 1.

Since the outflow at the vent vanishes (we are indeed modeling the relaxing phase), we consider these boundary conditions for v and

S

⁠:

{\begin{matrix} v = 0 & at z = 0, \\ S n = 0 & at z = h (x, y, t), \end{matrix}

(5)

where n is the outer normal to the top surface.

Concerning the top surface, we follow Ref. 36 recalling at first that

z = h (x, y, t)

is a material surface, so that

{w |}_{z = h} = \frac{\partial h}{\partial t} + {u |}_{z = h} \frac{\partial h}{\partial x} + {v |}_{z = h} \frac{\partial h}{\partial y} .

(6)

We then integrate ₁ over the domain depth and exploit and so as to get

\frac{\partial h}{\partial t} + \frac{\partial (h \bar{u})}{\partial x} + \frac{\partial (h \bar{v})}{\partial y} = 0,

(7)

where

h \bar{u} = \int_{0}^{h} u d z, h \bar{v} = \int_{0}^{h} v d z

⁠.

We introduce the following dimensionless quantities:

\begin{matrix} u^{*} = \frac{u}{U_{ref}}, v^{*} = \frac{v}{U_{ref}}, w^{*} = \frac{w}{ε U_{ref}}, \\ t^{*} = \frac{t}{t_{ref}}, x^{*} = \frac{x}{L_{ref}}, y^{*} = \frac{y}{L_{ref}}, z^{*} = \frac{z}{H_{ref}}, h^{*} = \frac{h}{H_{ref}}, \\ θ^{*} = \frac{θ - θ_{amb}}{θ_{ref} - θ_{amb}}, p^{*} = \frac{p - p_{0}}{ρ g H_{ref}}, D^{*} = \frac{L}{U} D, \end{matrix}

(8)

where

U_{ref}

is the characteristic planar velocity,

ε = \frac{H_{ref}}{L_{ref}} ≪ 1,

is the domain aspect ratio (⁠

H_{ref}

and

L_{ref}

are the characteristic thickness and length of the domain Ω),

t_{ref}

is the characteristic timescale (which will be investigated in Sec. V), p₀ and

θ_{amb}

are the external pressure and temperature, respectively, and

θ_{ref}

is a reference temperature. Therefore, after some algebra, system reduces to (in the sequel, we omit the “*” to keep notations as light as possible)

{\begin{array}{l} \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} + \frac{\partial w}{\partial_{z}} = 0, \\ ε Re (\frac{t_{tr}}{t_{ref}} \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y} + w \frac{\partial u}{\partial z}) \\ = - ε \frac{Re}{F r^{2}} \frac{\partial p}{\partial x} + ε \frac{\partial S_{x x}}{\partial x} + ε \frac{\partial S_{x y}}{\partial y} + \frac{\partial S_{x z}}{\partial z}, \\ ε Re (\frac{t_{tr}}{t_{ref}} \frac{\partial v}{\partial t} + u \frac{\partial v}{\partial x} + v \frac{\partial v}{\partial y} + w \frac{\partial v}{\partial z}) \\ = - ε \frac{Re}{F r^{2}} \frac{\partial p}{\partial y} + ε \frac{\partial S_{y x}}{\partial x} + ε \frac{\partial S_{y y}}{\partial y} + \frac{\partial S_{y z}}{\partial z}, \\ ε^{2} Re (\frac{t_{tr}}{t_{ref}} \frac{\partial w}{\partial t} + u \frac{\partial w}{\partial x} + v \frac{\partial w}{\partial y} + w \frac{\partial w}{\partial z}) \\ = - \frac{Re}{F r^{2}} (\frac{\partial p}{\partial z} + 1) + ε \frac{\partial S_{z x}}{\partial x} + ε \frac{\partial S_{z y}}{\partial y} + \frac{\partial S_{z z}}{\partial z}, \end{array}

(9)

where

t_{tr} = \frac{L_{ref}}{U_{ref}}

is the so-called transit time and

Re = \frac{ρ H_{ref}^{1 / λ} U_{ref}^{(2 λ - 1) / λ}}{μ_{c}}, F r^{2} = \frac{U_{ref}^{2}}{g L_{ref}}

(10)

are the Reynolds and Froude number, respectively.

Then, according to Ref. 36, we stipulate that:

The pressure gradient essentially balances the viscous stresses, so that

$U_{ref} = {(\frac{ρ g H_{ref}^{(2 λ + 1) / λ}}{μ_{c} L_{ref}})}^{λ},$

(11)

yielding

$1 = ε \frac{Re}{F r^{2}} .$

(12)
The inertia terms can be neglected,

$ε Re ≪ 1, ε Re \frac{t_{tr}}{t_{ref}} ≪ 1.$

(13)

Thus, by omitting the

O (ε^{n})

⁠,

n \geq 1

⁠, terms and by using , , and , system reduces to

{\begin{array}{l} \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} + \frac{\partial w}{\partial z} = 0, \\ \frac{\partial}{\partial z} [η (θ) {({(\frac{\partial u}{\partial z})}^{2} + {(\frac{\partial v}{\partial z})}^{2})}^{(1 - λ) / 2 λ} \frac{\partial u}{\partial z}] = 2^{(1 - λ) / 2 λ} \frac{\partial p}{\partial x}, \\ \frac{\partial}{\partial z} [η (θ) {({(\frac{\partial u}{\partial z})}^{2} + {(\frac{\partial v}{\partial z})}^{2})}^{(1 - λ) / 2 λ} \frac{\partial v}{\partial z}] = 2^{(1 - λ) / 2 λ} \frac{\partial p}{\partial y}, \\ \frac{\partial p}{\partial z} + 1 = 0. \end{array}

(14)

Now, we write, in dimensionless form, the kinematic condition for the top surface,

\frac{t_{tr}}{t_{ref}} \frac{\partial h}{\partial t} + \frac{\partial (h \bar{u})}{\partial x} + \frac{\partial (h \bar{v})}{\partial y} = 0,

(15)

and the boundary conditions where all

O (ε)

terms have been neglected,

{\begin{array}{l} {u |}_{h} = {v |}_{h} = {w |}_{h} = 0, \\ {p |}_{h} = 0, \\ {\frac{\partial u}{\partial z} |}_{h} = {\frac{\partial v}{\partial z} |}_{h} = 0, \end{array}

(16)

meaning that the external pressure has been rescaled to zero and no tangential stress acts on the top surface.

Generalizing the procedure illustrated in Ref. 36, we obtain the solution of the boundary-value problems and ,

{\begin{array}{l} u = - 2^{(1 - λ) / 2} \frac{\partial h}{\partial x} {‖ \nabla_{x, y} h ‖}^{λ - 1} \int_{0}^{z} {(\frac{h - s}{η (θ - θ_{amb})})}^{λ} d s, \\ v = - 2^{(1 - λ) / 2} \frac{\partial h}{\partial y} {‖ \nabla_{x, y} h ‖}^{λ - 1} \int_{0}^{z} {(\frac{h - s}{η (θ - θ_{amb})})}^{λ} d s, \\ w = 2^{(1 - λ) / 2} \nabla_{x, y} \cdot [{‖ \nabla_{x, y} h ‖}^{λ - 1} \nabla_{x, y} h \int_{0}^{z} (\int_{0}^{s} {(\frac{h - σ}{η (θ - θ_{amb})})}^{λ} d σ) d s], \\ p = h - z, \end{array}

(17)

where

\nabla_{x, y} = \frac{\partial}{\partial x} e_{x} + \frac{\partial}{\partial y} e_{y}

is the gradient in the x, y plane.

Introducing the functional dependence on h and θ,

F (h, θ) = \int_{0}^{h} \int_{0}^{z} {(\frac{h - s}{η (θ - θ_{amb})})}^{λ} d sdz = \int_{0}^{h} (h - s) {(\frac{h - s}{η (θ - θ_{amb})})}^{λ} d s,

(18)

we have

\begin{matrix} h \bar{u} = \int_{0}^{h} u d z = - 2^{(1 - λ) / 2} \frac{\partial h}{\partial x} {‖ \nabla_{x, y} h ‖}^{λ - 1} F (h, θ), \\ h \bar{v} = \int_{0}^{h} v d z = - 2^{(1 - λ) / 2} \frac{\partial h}{\partial y} {‖ \nabla_{x, y} h ‖}^{λ - 1} F (h, θ), \end{matrix}

and so Eq. acquires this form:

\frac{t_{tr}}{2^{(1 - λ) / 2} t_{ref}} \frac{\partial h}{\partial t} = \nabla_{x y} \cdot [F (h, θ) {‖ \nabla_{x, y} h ‖}^{λ - 1} \nabla_{x, y} h] .

(19)

To close the problem, we have to write also an evolution equation for the temperature field θ. The model obviously becomes quite complicated due to the coupling between the equation for h and the one for temperature. In the sequel, we first consider some approximations regarding both geometry and temperature.

III. MODELING THE RELAXING PHASE: ISOTHERMAL AND RADIAL APPROXIMATION

We assume that:

The phenomenon is essentially isothermal, and we take η = 1.
$\frac{t_{tr}}{t_{ref}} = O (1)$ ⁠. Actually, for the sake of simplicity, we take $\frac{t_{tr}}{t_{ref}} = 1$ ⁠. The validity of this hypothesis will be verified in Sec. V.
Cylindrical symmetry applies, which, considering cylindrical coordinates $r, φ, z$ ⁠, is equivalent to require: $h = h (r, t), v = v_{r} (r, z, t) e_{r} + w (r, z, t) e_{z}$ ⁠, and

$Ω = {(r, z, φ) \in ℝ^{3} | 0 \leq r \leq R (t), 0 \leq z \leq h (r, t), 0 \leq φ < 2 π},$

where $r = R (t)$ is the advancing front, namely $h (R (t), t) = 0$ ⁠.

If η = 1,

F

given by becomes

F (h) = \frac{h^{λ + 2}}{λ + 2},

and acquires this form in cylindrical coordinates:

\frac{λ + 2}{2^{(1 - λ) / 2}} \frac{\partial h}{\partial t} - \frac{1}{r} \frac{\partial}{\partial r} [r h^{λ + 2} {| \frac{\partial h}{\partial r} |}^{λ - 1} \frac{\partial h}{\partial r}] = 0.

(20)

We remark that, if λ = 1 (i.e., Newtonian fluid), Eq. reduces to the classical model,^1,6 in cylindrical coordinates.

Now, following a procedure similar to the one illustrated in Ref. 1, we look for a self-similar solution to taking

h (r, t) = β t^{- α} P (ξ) with ξ = γ r^{2} t^{- α},

(21)

with α, β, and γ suitable constant. In particular,

h (R (t), t) = 0

implies

{ξ |}_{R (t)} = γ {(R (t))}^{2} t^{- α} = const .

⁠, which provides the temporal behavior of

R (t)

⁠, i.e.,

R (t) = (constant) \cdot t^{α / 2} .

(22)

After plugging into , lengthy but straightforward algebraic manipulations allow us to recast in the following form (here,

P'

represents the ordinary derivative of P with respect to ξ):

\frac{d}{d ξ} {\frac{λ + 2}{2^{(1 + λ) / 2} (3 + 5 λ)} ξ P + β^{2 λ + 1} γ^{(1 + λ) / 2} ξ^{(1 + λ) / 2} P^{λ + 2} {| P^{'} |}^{λ - 1} P^{'}} = 0,

(23)

provided

α = \frac{2}{3 + 5 λ} .

(24)

Integrating , we obtain

\frac{λ + 2}{2^{(1 + λ) / 2} (3 + 5 λ)} ξ P + β^{2 λ + 1} γ^{(1 + λ) / 2} ξ^{(1 + λ) / 2} P^{λ + 2} {| P^{'} |}^{λ - 1} P^{'} = 0,

(25)

since the integration constant vanishes once the solution is evaluated at

ξ = 0

⁠. We now require that h is monotonically decreasing with respect to

r

⁠, i.e.,

P^{'} < 0

⁠. We thus rewrite as

- \frac{2 λ + 1}{λ} S^{1 / λ} ξ^{(1 - λ) / (2 λ)} = {(P^{(2 λ + 1) / λ})}^{'}

(26)

with

S (γ, β, λ) = \frac{2^{- (1 + λ) / 2} (λ + 2)}{(3 + 5 λ) β^{2 λ + 1} γ^{(1 + λ) / 2}} .

Integrating again ,

P (ξ) = {[C - \frac{2 (2 λ + 1)}{1 + λ} S^{1 / λ} ξ^{(1 + λ) / (2 λ)}]}^{λ / (2 λ + 1)}

with C constant, which gives

h (r, t, λ) = β t^{- 2 / (3 + 5 λ)} {[C - \frac{2 (2 λ + 1)}{1 + λ} S^{1 / λ} γ^{(1 + λ) / (2 λ)} {(\frac{r}{t^{1 / (3 + 5 λ)}})}^{(1 + λ) / λ}]}^{λ / (2 λ + 1)} .

(27)

It is worth highlighting that the latter reduces, when λ = 1, to formula (A13) of Ref. 1 provided we assume

ν_{0} / g = 1

in (A13) itself. We then impose

h (R (t), t, λ) = 0

and obtain with α given by , i.e.,

R (t, λ) = \frac{C (1 + λ)}{2 (2 λ + 1) S^{1 / λ} γ^{(1 + λ) / (2 λ)}} t^{1 / (3 + 5 λ)}

(28)

and

h (r, t, λ) = H (λ) t^{- 1 / (2 λ + 1)} R {(t, λ)}^{(1 + λ) / (2 λ + 1)} {[1 - {(\frac{r}{R (t, λ)})}^{(1 + λ) / λ}]}^{λ / (2 λ + 1)}

(29)

with

H (λ) = 2^{(λ - 1) / (4 λ + 2)} {[\frac{2 λ + 1}{1 + λ} {(\frac{λ + 2}{3 + 5 λ})}^{1 / λ}]}^{λ / (2 λ + 1)} .

(30)

We point out that

R (t, λ)

is defined up to the constant

C

and so is h. For a physical interpretation of this constant, we proceed recalling that provides the height of the dome for each

t \geq 0

⁠. However, we are considering the relaxing phase, that is, the phase in which the lava flow from the vent is over (and the dome volume remains unchanged). Therefore, if the growing or effusion phase ended at the (dimensionless) time τ, model applies for

t > τ

⁠. So, following Ref. 1, we set

t = (t - τ) + τ = \tilde{t} + τ

⁠, with

\tilde{t} \geq 0

⁠, and obtain

h (r, \tilde{t}, τ, λ) = H (λ) {(\tilde{t} + τ)}^{- 1 / (2 λ + 1)} R {(\tilde{t} + τ, λ)}^{(1 + λ) / (2 λ + 1)} {[1 - {(\frac{r}{R (\tilde{t} + τ, λ)})}^{(1 + λ) / λ}]}^{λ / (2 λ + 1)} .

(31)

It is worth noticing that h does not depend on proportionality constant

β

introduced through the similarity solution . In particular, if we introduce

Θ (λ) = \int_{0}^{1} {(1 - ξ^{(1 + λ) / λ})}^{λ / (2 λ + 1)} ξ d ξ,

which reduces to 3/8 when λ = 1, the volume V of the dome at time

\tilde{t}

is

V (\tilde{t}, τ, λ) = 2 π \int_{0}^{R (\tilde{t} + τ, λ)} h (r, \tilde{t}, τ, λ) r d r = 2 π H (λ) Θ (λ) {[\frac{R {(\tilde{t} + τ, λ)}^{(3 + 5 λ)}}{\tilde{t} + τ}]}^{1 / (2 λ + 1)} .

(32)

On the other hand, the aim of the model is to justify someway the estimated profiles of the dome after the effusion phase has ended and the volume has reached its maximum value. In other words, in the present context, the initial time corresponds to a given instant

τ > 0

⁠, which is related to the position of the front R_o and the volume V_o at the end of the effusion phase through the following relation:

V_{o} = 2 π H (λ) Θ (λ) {[\frac{R_{o}^{3 + 5 λ}}{τ}]}^{1 / (2 λ + 1)},

(33)

which provides τ once R_o, λ, and V_o are prescribed, namely,

τ = {(\frac{2 π H (λ) Θ (λ)}{V_{o}})}^{2 λ + 1} R_{o}^{3 + 5 λ} .

(34)

Function

Θ

is not explicit in terms of λ, unless λ = 1. However, in the interval (0,1] it turns out

Θ (λ) = 3 / 8 + 0.68 (1 - λ) + 0.04 {(1 - λ)}^{2} + O [{(1 - λ)}^{3}]

⁠. For λ = 1, implies

τ = {(\frac{3}{4})}^{5} {(\frac{π}{V_{o}})}^{3} R_{o}^{8}

⁠, exactly the same as formula (A29) in Ref. 1. Furthermore, since

V (\tilde{t}, τ, λ) = V_{o}

for

\tilde{t} \geq 0

⁠, formula gives the advancing front, namely,

R (\tilde{t} + τ, λ) = \underset{R_{o} τ^{- 1 / (3 + 5 λ)}}{\underset{}{\underset{︸}{{(\frac{V_{o}}{2 π H (λ) Θ (λ)})}^{(1 + 2 λ) / (3 + 5 λ)}}}} {(\tilde{t} + τ)}^{1 / (3 + 5 λ)} = R_{o} {(1 + \frac{\tilde{t}}{τ})}^{1 / (3 + 5 λ)}, \tilde{t} \geq 0,

(35)

which reduces to formula (A26) in Ref. 1 when λ = 1. We finally remark that allows linking the constant C, appearing in , the dome volume at the end of the growth phase, i.e.,

V_{o},

or, exploiting , to

R_{o}

and τ (Table I).

TABLE I.

Summary of the best fitting values of V_o, ${\tilde{t}}_{end}$ ⁠, λ, and R_o. The value of τ is computed through Eq. (34). We notice that values remain practically unchanged for the four different profiles considered in the present paper.

$Dome profile$	V_o	${\tilde{t}}_{end}$	λ	R_o	Τ
DEM (SW-NE)	2.1	0.43	0.8	0.66	0.04
DEM (N-S)	2.0	0.45	0.8	0.65	0.04
SAR (orbit No. 1277)	2.2	0.45	0.8	0.67	0.04
SAR (orbit No. 3067)	2.3	0.43	0.8	0.67	0.036

$Dome profile$	V_o	${\tilde{t}}_{end}$	λ	R_o	Τ
DEM (SW-NE)	2.1	0.43	0.8	0.66	0.04
DEM (N-S)	2.0	0.45	0.8	0.65	0.04
SAR (orbit No. 1277)	2.2	0.45	0.8	0.67	0.04
SAR (orbit No. 3067)	2.3	0.43	0.8	0.67	0.036

Finally, as shown in Fig. 5, because of (33), τ increases with increasing R_o and V_o, independently of λ.

FIG. 5.

View large Download slide

Contour lines of $τ (R_{o}, V_{o}, λ)$ obtained from (33) for $λ = 0.6 - 0.9$ ⁠.

IV. SOME SIMULATIONS USING BOTH SAR AND DEM ALTIMETRY DATA

In this section, we use data reported in Refs. 25 and 27. Figure 3 shows that the diameter of the dome is $\approx 30 km$ ⁠, the height being $\approx 1.5 km$ so that the aspect ratio $ε \approx 0.05$ ⁠. We make dimensionless the original geometrical data using as $L_{ref}$ and $H_{ref}$ the dome maximum radius and height (whose ranges have been just mentioned), respectively. The specific values of $L_{ref}$ and $H_{ref}$ depend on the profile considered. Those used in our simulations are listed in Table II. We define ${\tilde{t}}_{end}$ as the emplacement time, i.e., the duration of the relaxation phase.

TABLE II.

Dimensional volume $V_{model} = V_{o} H_{ref} L_{ref}^{2}$ of the dome in Fig. 2 according to four different profiles, two based on DEM data, and two on SAR data, respectively. The values provided by the mathematical model are compared with $V_{estim .}$ as reported in the literature. Our result appears to agree very well with those suggested by Gleason et al.²⁵ On the contrary, for the case of SAR data, the values proposed by McKenzie et al.²⁷ agree neither with ours nor with Gleason's et al. ones. It should be noticed, however, that due to the automated methodology used to process the Magellan altimetry data, there are significant differences between the altimetry derived measurements of Pavri et al.^26,27 and the stereo-derived DEMs by Gleason et al. In the last column, we reported the volume estimated by McKenzie et al. for orbit Nos. 1277 and 3067.

$Dome profile$	$H_{ref} (km)$	$L_{ref} (km)$	$V_{model} {(km}^{3})$	$V_{estim .} ({km}^{3})$	Reference
DEM (SW-NE)	1.48	15.9	748	795	25
DEM (N-S)	1.43	15.4	712	706	25
SAR (orbit No. 1277)	1.45	14.9	691	1040	27
SAR (orbit No. 3067)	1.41	13.9	618	820	27

As a first step, we express $h (r, {\tilde{t}}_{end})$ ⁠, given by (31), in terms of the dimensionless parameters V_o, R_o, $λ$ ⁠, and ${\tilde{t}}_{end}$ exploiting Eqs. (34) and (35). Next, we vary the dimensionless parameters V_o, R_o, $λ$ ⁠, and ${\tilde{t}}_{end}$ in order to fit the following profiles: orbit No. 1277, orbit No. 3067, considered in Ref. 27 and the N-S, SW-NE profiles considered in Ref. 25. To further clarify the procedure we followed, some remarks are appropriate:

The fitting parameters V_o, R_o, λ and ${\tilde{t}}_{end}$ are selected within appropriate ranges so that the dimensional values are compatible with the orders of magnitude of the experimentally observed scales.
Once V_o, R_o, λ, and ${\tilde{t}}_{end}$ have been selected, the dimensionless parameter τ is computed through the formula (34). Its value provides an estimate of the duration of the effusive phase with variable volume with respect to the relaxation phase at constant volume.
Finally, ${\tilde{t}}_{end}$ gives an estimate of the duration (dimensionless) of the relaxation phase at constant volume.
Dimensional values for both the effusive and relaxation durations can be estimated once a correct timescale $t_{ref}$ has been identified. This issue will be discussed in Sec. V.

In Table II, we list the corresponding dimensional volumes for each profile analyzed. In this connection, we remark that McKenzie et al. use $(3 π / 4) H_{ref} L_{ref}^{2} \approx 2.36 H_{ref} L_{ref}^{2}$ as a scale factor for the volume, that is by assuming a dome dimensionless profile proportional to ${(1 - r)}^{1 / 2}$ [see Eq. (1) in Ref. 27]. Figures 6 and 7 show the best fitting of the dome profiles based on both the DEM and SAR data, respectively. In particular, in Fig. 6 we also report the best fitting profiles (in red), if we take a Newtonian rheological model (λ = 1). We notice that, in the latter case we obtain a significant change in the value of τ, difficult to be explained physically (the dome is the same). On the contrary, considering a non-Newtonian model, λ provides an extra degree of freedom to fit the dome profile, with practically identical values of V_o, ${\tilde{t}}_{end}$ ⁠, R_o, and τ [calculated through (33)].

FIG. 6.

Best fitting of the dome profiles based on the DEM data (N-E on the left, SW-NE on the right). The estimated values of V o , t ̃ end , R o obtained in the non-Newtonian case ( λ = 0.8) are practically identical [ τ is calculated through (33)]. On the contrary, neglecting shear thinning (λ = 1), the values of τ that realize the best fitting are significantly different.

View large Download slide

Best fitting of the dome profiles based on the DEM data (N-E on the left, SW-NE on the right). The estimated values of $V_{o}, {\tilde{t}}_{end}, R_{o}$ obtained in the non-Newtonian case (⁠ $λ = 0.8$ ⁠) are practically identical [ $τ$ is calculated through (33)]. On the contrary, neglecting shear thinning (λ = 1), the values of τ that realize the best fitting are significantly different.

FIG. 7.

Best fitting of the dome profiles based on the SAR data (orbit No. 1277 on the left, orbit No. 3067 on the right). The estimated values of V o , t ̃ end , λ , R o are very similar [ τ is calculated through (33)].

View large Download slide

Best fitting of the dome profiles based on the SAR data (orbit No. 1277 on the left, orbit No. 3067 on the right). The estimated values of $V_{o}, {\tilde{t}}_{end}, λ, R_{o}$ are very similar [ $τ$ is calculated through (33)].

Figure 8 shows the temporal evolution of the profile of the dome starting from τ, i.e., from $\tilde{t} = 0$ ⁠. As can be seen in Table II, the (dimensionless) duration of the relaxing phase is ${\tilde{t}}_{end} ≃ 0.4$ ⁠, and in this interval, the (dimensional) radius of the dome has grown from approximately $0.67 L_{ref}$ to $L_{ref}$ ⁠.

FIG. 8.

Dome relaxation at constant volume starts at t ̃ = 0. The mathematical model predicts that the dome attained the present shape at t ̃ = t ̃ end ≃ 0.44.

View large Download slide

Dome relaxation at constant volume starts at $\tilde{t} = 0$ ⁠. The mathematical model predicts that the dome attained the present shape at $\tilde{t} = {\tilde{t}}_{end} ≃ 0.44$ ⁠.

V. DEALING WITH THE TIME SCALE

So far, we have not addressed the estimation of the timescale

t_{ref}

⁠, just assuming that, in the isotherm approximation,

t_{ref} / t_{tr} = O (1)

⁠. Which could be a reasonable timescale to evaluate the dimensional values of

{\tilde{t}}_{end}

and τ? This problem is strictly related to the cooling process and its influence on the viscosity, which in turn depends on the real chemical composition of the magma. As pointed out by Huppert,⁶ the observed shape of the dome can be used for an estimate of the apparent viscosity of the magma. According to McKenzie et al.,²⁷

t_{ref} \approx \frac{H_{ref}^{2}}{π^{2} κ},

(36)

where

κ = O (10^{- 6} - 10^{- 7}) m^{2} s^{- 1}

is the thermal diffusivity.

Remark 1. The diffusion coefficient κ is slightly decreasing in the range 50–1000 $° C$ ⁠. For basaltic lava,³⁷ $κ$ remains in the range 7–4 × $10^{- 6} m^{2} s^{- 1}$ ⁠, while for rhyolitic melts³⁸ κ remains in the range 6–4 × $10^{- 7} m^{2} s^{- 1}$ ⁠. We notice that the uncertainty of one order of magnitude for κ implies that $t_{ref}$ can range, in the present case, from hundreds to thousands of years.

Therefore, the dimensional duration of the relaxing phase is

t_{ref} {\tilde{t}}_{end}

and the instant at which the volume of the dome has reached its final (i.e., current) value V_o (the effusion phase is over) is

t_{ref} τ

⁠. It should be said, however, that the scale refers to magma cooling, which is an extremely complex phenomenon^39,40 and to identify a correct timescale would require much more knowledge about the dome chemical composition rather than just its shape. Nevertheless, we can use the above value of

t_{ref}

and the dimensionless values provided by the mathematical model and listed in Table II to obtain an estimate of μ_c. Indeed, writing in a dimensional form and using the selected scales in Table II, we obtain

H_{ref} L_{ref}^{2} V_{o} = 2 π {(\frac{μ_{c}}{ρ g})}^{λ / (2 λ + 1)} H (λ) Θ (λ) {[\frac{{(L_{ref} R_{o})}^{3 + 5 λ}}{τ t_{ref}}]}^{1 / (2 λ + 1)},

(37)

from which we compute μ_c. We recall¹ that, on Venus,

g = 8.87 {ms}^{- 2}

and that

ρ = 2535 {kgm}^{- 3}

is the density difference between the flowing material (lava has a density of

2600 {kgm}^{- 3}

⁠) and the ambient fluid density (Venus atmosphere has a density of

65 {kgm}^{- 3}

⁠).

Our simulations indicate that the fitting of all profiles available occur, on average, for

λ ≃ 0.8, V_{o} ≃ 2.2, {\tilde{t}}_{end} ≃ 0.44, R_{o} ≃ 0.66

⁠. We first notice that, being the transit timescale estimated by

t_{tr} = t_{ref} {\tilde{t}}_{end},

simulations indicate that the assumption

t_{tr} / t_{ref} = O (1)

is essentially correct. In all simulations, the corresponding value of τ provided by is

τ \sim 0.04.

Letting now, on average,

L_{ref} = 1.5 \times 10^{4} m, H_{ref} = 1.45 \times 10^{3} m

⁠, we first notice that during the relaxation phase at constant volume, the dome front advances

4.9 \times 10^{3} m

⁠. If we assume, on average (see Remark 1),

κ ≃ 5 \times 10^{- 6} m^{2} s^{- 1}

⁠, then

t_{ref} ≃ 4.26 \times 10^{10} s

(that is 1350 yr), and, as a consequence, the emplacement time is a little less than a half of that, i.e., about 608 yr. In particular, the total duration of the process is

{\tilde{t}}_{end} + τ ≃ 0.54

⁠. In Table III, we have summarized the values obtained.

TABLE III.

Summary of the dimensionless/dimensional timescales involved in the dome evolution.

Duration	dimensionless	dimensional (⁠ $t_{ref} ≃ 1350$ yr)
Effusion	$τ ≃ 0.04$	$≃ 54$ (yr)
Relaxing	${\tilde{t}}_{end} ≃ 0.45$	$≃ 608$ (yr)
Dome age	$τ + {\tilde{t}}_{end} ≃ 0.49$	$≃ 662$ (yr)

Duration	dimensionless	dimensional (⁠ $t_{ref} ≃ 1350$ yr)
Effusion	$τ ≃ 0.04$	$≃ 54$ (yr)
Relaxing	${\tilde{t}}_{end} ≃ 0.45$	$≃ 608$ (yr)
Dome age	$τ + {\tilde{t}}_{end} ≃ 0.49$	$≃ 662$ (yr)

By definition, the velocity of the front during the relaxation phase is

U_{ref} = \frac{(1 - R_{o}) L_{ref}}{t_{ref} {\tilde{t}}_{end}},

that is

2.58 \times 10^{- 7} {ms}^{- 1}

(that is 8.1 m/yr). Recalling the definition of the scale

U_{ref}

⁠, we can solve for μ_c to get

μ_{c} ≃ 5.4 \times 10^{18} (Pa s^{5 / 4}) .

(38)

One can easily verify that and provide estimates of the same order of magnitude. Indeed, from we get

μ_{c} ≃ 3.63 \times 10^{18} (Pa s^{5 / 4})

⁠. Although quite rough, such an approach is useful as a guide to the likely order of magnitude of the (apparent) viscosity μ_c at the end of the relaxation phase, that is the present viscosity of the material. To get reliable estimates of the behavior of viscosity with time from the onset of the effusive phase requires many more information.^5,41,42 For example, we should know more about the chemical composition of magma, the initial water content and how the crystallization percentage change as temperature decreases during the cooling process. This kind of information is generally stated by adopting a dependency of viscosity on temperature of exponential type (such as the Arrhenius law).

We now use (38) to verify that conditions (12) and (13) are largely satisfied (see Table IV).

TABLE IV.

Estimated values of $U_{ref}$ ⁠, $Re$ ⁠, $Fr$ ⁠, and $ε Re / F r^{2}$ during the relaxation phase, based on simulations and (38).

$U_{ref} ({ms}^{- 1})$	$Re$	$F r^{2}$	$ε = H_{ref} / 2 L_{ref}$	$ε Re / F r^{2}$
$≃ 2.3 \times 10^{- 7}$	$≃ 3.97 \times 10^{- 17}$	$≃ 4.26 \times 10^{- 19}$	$≃ 0.05$	$≃ 4.82$

In Sec. VI, we investigate the effusion phase, relaxing the isothermal assumption. Indeed, we couple the dynamics of h with the cooling process in order to get some insights on the role of temperature on the dome formation process.

VI. MODELING OF THE EFFUSION PHASE: COUPLING THE FLOW WITH THE THERMAL FIELD

We now include in the model the energy balance equation which, still assuming the radial approximation, is written as

ρ c_{p} (\frac{\partial θ}{\partial t} + {\bar{v}}_{r} \frac{\partial θ}{\partial r}) = K \frac{1}{r} \frac{\partial}{\partial r} (r \frac{\partial θ}{\partial r}) - E σ \frac{θ^{4} - θ_{amb}^{4}}{h},

(39)

where c_p is the specific heat capacity at constant pressure, E is the emissivity,

σ = 5.67 \times 10^{- 8} {Wm}^{- 2} K^{- 4}

the Stefan–Boltzmann constant, K the thermal conductivity,

θ_{amb}

the external ambient temperature (Venus average ambient temperature on the dome surface), and

{\bar{v}}_{r} (r, t) = \frac{1}{h (r, t)} \int_{0}^{h (r, t)} v_{r} (r, z, t) d z

(40)

is the mean radial velocity. To justify the division by h in the radiation term at the rhs of , we recall this term must represent the heat loss per unit volume, while

E σ (θ^{4} - θ_{amb}^{4})

gives, by definition, the heat loss per unit area. Introducing the dimensionless temperatures

θ^{*} = \frac{θ}{θ_{ref}}, θ_{amb}^{*} = \frac{θ_{amb}}{θ_{ref}},

where

θ_{ref}

is the temperature of the lava at the effusive vent, and using , Eq. can be rewritten in dimensionless form (remind that we have omitted the “*” to keep notation as light as possible) as

\frac{1}{t_{ref}} \frac{\partial θ}{\partial t} + \frac{1}{t_{tr}} {\bar{v}}_{r} \frac{\partial θ}{\partial r} = \frac{1}{t_{diff}} \frac{1}{r} \frac{\partial}{\partial r} (r \frac{\partial θ}{\partial r}) - \frac{1}{t_{rad}} \frac{θ^{4} - θ_{amb}^{4}}{h},

(41)

where

t_{diff} = \frac{L_{ref}^{2} ρ c_{p}}{K},

(42)

t_{rad} = \frac{ρ c_{p} H_{ref}}{E σ θ_{amb}^{3}}

(43)

are the time scales featuring the diffusive and radiation processes, respectively. The dimensional values of

θ_{ref}

and

θ_{amb}

are reported in Table V which completes the list of physical parameters relevant to the cooling dynamics. In particular, referring to Table V,

θ_{amb} ≃ 0.54

⁠, so that

θ_{amb}^{4} = O (ε) .

TABLE V.

Order of magnitude (or estimated values) of some relevant parameters for Venusian domes^{20,38–40,43,44} and the corresponding characteristic radiation and diffusive timescale.

c_p (⁠ ${Jkg}^{- 1} K^{- 1}$ ⁠)	K (⁠ ${Wm}^{- 1} K^{- 1}$ ⁠)	E	$θ_{ref}$ (K)	$θ_{amb}$ (K)	$t_{rad}$ (s)	$t_{diff}$ (s)
$1.2 \times 10^{3}$	15	0.8	$1.3 \times 10^{3}$	$7 \times 10^{2}$	$5.6 \times 10^{7}$	$4.2 \times 10^{10}$

c_p (⁠ ${Jkg}^{- 1} K^{- 1}$ ⁠)	K (⁠ ${Wm}^{- 1} K^{- 1}$ ⁠)	E	$θ_{ref}$ (K)	$θ_{amb}$ (K)	$t_{rad}$ (s)	$t_{diff}$ (s)
$1.2 \times 10^{3}$	15	0.8	$1.3 \times 10^{3}$	$7 \times 10^{2}$	$5.6 \times 10^{7}$	$4.2 \times 10^{10}$

We now focus on the effusion phase, whose duration,

τ t_{ref}

⁠, has been estimated in Sec. V. In particular,

τ t_{ref} \sim 10^{- 2} t_{ref} \sim 10^{8} s

⁠, with

t_{ref}

given by . Since the characteristic velocity

U_{ref}

during the dome formation can be estimated as

R_{o} L_{ref} / τ t_{ref}

⁠, the corresponding transit time is

t_{tr} = τ t_{ref}

⁠. We thus select

τ t_{ref}

as new timescale and notice that

\frac{τ t_{ref}}{t_{rad}} = O (1)

⁠, while

\frac{τ t_{ref}}{t_{diff}} = O (ε)

⁠. Hence, neglecting

O (ε^{n}), n \geq 1

⁠, Eq. reduces to

\frac{\partial θ}{\partial t} + {\bar{v}}_{r} \frac{\partial θ}{\partial r} = - r \frac{θ^{4}}{h},

(44)

with

r = \frac{τ t_{ref}}{t_{rad}} = O (1) .

Equation is physical consistent only if θ varies in the interval

(θ_{amb}, 1)

⁠, where 1 is the dimensionless temperature of the fluid at the vent. Concerning conditions , these are still largely fulfilled.

Equation must be coupled with the radial version of , i.e.,

\frac{1}{2^{(1 - λ) / 2}} \frac{\partial h}{\partial t} - \frac{1}{r} \frac{\partial}{\partial r} [r F (h, θ) {| \frac{\partial h}{\partial r} |}^{λ - 1} \frac{\partial h}{\partial r}] = 0,

(45)

supplemented with the condition

V (t) = 2 π \int_{0}^{R (t)} h (r, t) r d r,

(46)

where

V (t)

is the dome volume at time t and where we recall that

F (h, θ) = \int_{0}^{h} (h - s) {(\frac{h - s}{η (θ - θ_{amb})})}^{λ} d s = \frac{h^{λ + 2}}{η^{λ} (θ - θ_{amb}) (λ + 2)} .

(47)

To close the problem, we must assign a physically significant function

η (θ - θ_{amb})

⁠, appearing in the constitutive equation . We could assume, for example, the following Arrhenius-type law, that is

η (θ - θ_{amb}) = exp [b {(θ - θ_{amb})}^{δ}],

(48)

where b and δ are the dimensionless parameters that depend on the chemical composition and other characteristics of the fluid, and

θ \in (θ_{amb}, 1)

⁠, with

0 < θ_{amb} < 1

⁠. Viscosity can vary of orders of magnitude, at a given temperature, changing the nature of magma (rhyolitic, basaltic, etc.) and the percentage of water content.^45–47 Our model suggests that

μ_{c} = O (10^{18})

is the present (dimensional) viscosity of the fluid and many researchers argued about the viscosity at the eruption time. Figure 9 shows

μ_{c} η (θ - θ_{amb})

in two cases,

b = - 10, - 15

⁠, with

δ = 1 / 3

⁠. We remark that, as shown in Ref. 47, the precise values of viscosity vary in a wideband depending on the magma water percentage. Therefore, deductions based on , should be viewed more from a qualitative viewpoint than from a quantitative one.

FIG. 9.

View large Download slide

Possible behavior of viscosity with (dimensionless) temperature based on (48) with $δ = 1 / 3$ ⁠.

VII. UNIFORM TEMPERATURE APPROXIMATION

We now assume that:

The temperature is spatially uniform, i.e., $θ = θ (t)$ ⁠, and, inspired by (44), we stipulate that its evolution equation is

$\dot{θ} = - A θ^{4},$

(49)

where $A$ is a suitable positive parameter which we will specify shortly.
At t = 1, namely, at the end of the effusion phase, the fluid has thermalized with ambient, while at the beginning of the process, i.e., at t = 0, the temperature is maximum, i.e.,

$\begin{array}{l} θ (1) = θ_{amb}, θ (0) = 1. \end{array}$

(50)

We integrate backward in time (i.e., for

t \leq 1

⁠) Eq. with the “initial” condition

_{1}

⁠, getting

θ (t) = θ_{amb} {(1 - 3 A θ_{amb}^{3} (1 - t))}^{- 1 / 3}

⁠. We then impose ₂ obtaining

θ (t) = \frac{θ_{amb}}{\sqrt[3]{1 - (1 - θ_{amb}^{3}) (1 - t)}},

so that

A = \frac{1 - θ_{amb}^{3}}{θ_{amb}^{3}} .

(51)

Using , Eq. becomes

η^{λ} (θ - θ_{amb}) \frac{\partial h}{\partial t} - \frac{2^{(1 - λ) / 2}}{(λ + 2) r} \frac{\partial}{\partial r} [r h^{λ + 2} {| \frac{\partial h}{\partial r} |}^{λ - 1} \frac{\partial h}{\partial r}] = 0.

(52)

Next, we remark that θ and t are in a one-to-one correspondence, and we can write

t = t (θ)

⁠, with

t (θ) = 1 - \frac{1}{1 - θ_{amb}^{3}} [1 - {(\frac{θ_{amb}}{θ})}^{3}],

(53)

and look for

\tilde{h} (r, θ) = h (r, t (θ))

⁠, so that

\frac{\partial h}{\partial t} = \frac{\partial \tilde{h}}{\partial θ} \dot{θ} = - A θ^{4} \frac{\partial \tilde{h}}{\partial θ} .

Equation can be rewritten as

- A θ^{4} η^{λ} (θ - θ_{amb}) \frac{\partial \tilde{h}}{\partial θ} - \frac{2^{(1 - λ) / 2}}{(λ + 2) r} \frac{\partial}{\partial r} [r {\tilde{h}}^{λ + 2} {| \frac{\partial \tilde{h}}{\partial r} |}^{λ - 1} \frac{\partial \tilde{h}}{\partial r}] = 0.

We now introduce the auxiliary variable

φ = φ (θ) \Leftrightarrow θ = θ (φ)

defined through the differential equation

- A θ^{4} η^{λ} (θ - θ_{amb}) \frac{\partial φ}{\partial θ} = 1.

(54)

Recalling , gives

φ (θ) = φ_{o} - \frac{3}{A} Λ (θ; λ, b, δ),

(55)

where

Λ (θ; λ, b, δ) = \int_{θ_{amb}}^{θ} s^{- 4} exp [- λ b {(s - θ_{amb})}^{δ}] d s,

and the initial value

φ_{o} = φ (θ_{amb})

has to be conveniently chosen. We require

φ

to be a non-negative function of

θ \in (θ_{amb}, 1)

⁠. Thus, we select (see Fig. 10)

φ_{o} = \frac{3}{A} max_{θ, λ, b, δ} Λ (θ; λ, b, δ) = \frac{3}{A} max_{λ, b, δ} Λ (1; λ, b, δ) .

(56)

FIG. 10.

Behavior of φ o given (56), for λ = 0.8 , δ = 1 / 3, and b ∈ ( − 15 , − 10 ).

View large Download slide

Behavior of $φ_{o}$ given (56), for $λ = 0.8, δ = 1 / 3$ ⁠, and $b \in (- 15, - 10)$ ⁠.

Figures 11 and 12 show Λ and $φ$ ⁠, respectively, as functions of $θ \in (θ_{amb}, 1)$ ⁠, for prescribed values of b, δ, and λ. Figure 13 shows φ, directly, as a function of t.

Clearly, by inverting both $t = t (θ)$ in (53) and $φ = φ (θ)$ in (55) we can also evaluate t as a function of $φ$ and vice versa.

FIG. 11.

Behavior of φ ( θ ) for λ = 0.8 , δ = 1 / 3, and b = − 15 , − 10, with θ ∈ ( θ amb , 1 ) ≃ ( 0.54 , 1 ).

View large Download slide

Behavior of $φ (θ)$ for $λ = 0.8, δ = 1 / 3$ ⁠, and $b = - 15, - 10$ ⁠, with $θ \in (θ_{amb}, 1) ≃ (0.54, 1)$ ⁠.

FIG. 12.

Representation of Λ ( θ ) for λ = 0.8 , δ = 1 / 3, and b = − 15 , − 10 for θ ∈ ( θ amb , 1 ) ≃ ( 0.54 , 1 ).

View large Download slide

Representation of $Λ (θ)$ for $λ = 0.8, δ = 1 / 3$ ⁠, and $b = - 15, - 10$ for $θ \in (θ_{amb}, 1) ≃ (0.54, 1)$ ⁠.

FIG. 13.

View large Download slide

Behavior of $φ (θ (t))$ for $λ = 0.8, δ = 1 / 3$ ⁠, and $b = - 15, - 10$ ⁠.

Then, we look for

\hat{h} (r, φ) = \tilde{h} (r, θ (φ)) = h (r, t (θ (φ)))

⁠, whose equation is

\frac{\partial \hat{h}}{\partial φ} - \frac{2^{(1 - λ) / 2}}{(λ + 2) r} \frac{\partial}{\partial r} [r {\hat{h}}^{λ + 2} {| \frac{\partial \hat{h}}{\partial r} |}^{λ - 1} \frac{\partial \hat{h}}{\partial r}] = 0.

(57)

To solve Eq. , differently from where

h (r, t) = β t^{- α} P (ξ)

⁠, we now seek the solution in the form

\hat{h} (r, φ) = φ^{n - α} P (ξ) with ξ = γ r^{2} φ^{- α},

(58)

with

γ, n, α

positive constants to be conveniently chosen (as it will be specified in the sequel). We remark that

φ \geq 0

⁠, because of . In particular, since by definition

h (R (t), t) = 0

for every t, we impose

P (1) = 0

⁠, yielding

R (φ (θ (t))) = \sqrt{\frac{φ^{α} (θ (t))}{γ}} .

(59)

On the other hand, exploiting , the dimensionless volume of the dome writes

V (t) = 2 π \int_{0}^{R (t)} h (r, t) r d r = \frac{π φ {(θ)}^{n}}{γ} \int_{0}^{γ R {(φ)}^{2} φ {(θ)}^{- α}} P (ξ) d ξ .

(60)

In particular, at the end of the effusion phase, i.e., at t = 1, we have

θ = θ_{amb}, φ_{o} = φ (θ_{amb})

⁠, and so,

V (θ_{amb}) = \frac{π φ_{o}^{n}}{γ} \int_{0}^{γ R {(θ_{amb})}^{2} φ_{o}^{- α}} P (ξ) d ξ,

(61)

where we recall that, according to Table I,

R (θ_{amb}) = R_{o} ≃ 0.65

⁠, and

V (θ_{amb}) = V_{o} ≃ 2.2

⁠. We can simplify choosing

γ = R_{o}^{- 2} φ_{o}^{α},

(62)

so that

V (θ) = 2 π \int_{0}^{R (θ)} h (r, t) r d r = π φ^{n - α} (θ) R_{o}^{2} \int_{0}^{1} P (ξ) d ξ,

(63)

yielding

\int_{0}^{1} P (ξ) d ξ = \frac{φ_{o}^{α - n}}{π R_{o}^{2}} V_{o},

(64)

which plays the role of a side condition for the unknown function

P (ξ)

⁠. In addition, rewrites

R (φ (θ (t))) = R_{o} {[\frac{φ (θ (t))}{φ_{o}}]}^{α / 2} .

(65)

Plugging into , we obtain the following differential equation in ξ:

- \frac{n (λ + 2)}{2^{\frac{3 + λ}{2}}} P + \frac{d}{d ξ} [\frac{λ + 2}{2^{\frac{λ + 1}{2}}} α ξ P + γ^{1 + 2 λ} ξ^{\frac{λ + 1}{2}} P^{λ + 2} {| P^{'} |}^{λ - 1} P^{'}] = 0,

(66)

provided

α = \frac{2 (1 + n + n λ)}{3 + 5 λ} .

Let us define

A_{1} = - \frac{n (λ + 2)}{2^{(3 + λ) / 2}}, A_{2} = \frac{λ + 2}{2^{(λ + 1) / 2}} α, A_{3} = γ^{1 + 2 λ},

with γ given by . As we already noticed, h is monotonically decreasing with respect to r, i.e.,

P' < 0

⁠. Therefore, we rewrite as follows:

\frac{A_{1}}{A_{3}} P + \frac{d}{d ξ} [\frac{A_{2}}{A_{3}} ξ P - ξ^{\frac{λ + 1}{2}} P^{λ + 2} {| P^{'} |}^{λ}] = 0.

(67)

The boundary conditions to be coupled to are

{\begin{array}{l} P (1) = 0, \\ \int_{0}^{1} P (ξ) d ξ = \frac{φ_{o}^{α - n}}{π R_{o}^{2}} V_{o} . \end{array}

(68)

VIII. NUMERICAL SIMULATIONS

Clearly, to integrate the differential equation for P is the key point. We verified, “a posteriori,” that changing either n, δ, or λ has a very modest effect on the shape of P. On the other hand, changing b modifies significantly the value of the integral of P, and, consequently, the parameter b is crucial in order to satisfy the non-local condition

_{2}

⁠. Indeed, to solve and , we applied a shooting method allowing large values of

| P' (1) |

(which is theoretically infinite) and changing

P' (1)

until

_{2}

is met. To limit the variability of the parameters involved, we fixed

λ = 0.8

⁠, n = 0.3, and

δ = 1 / 3

⁠, allowing b to vary of one order of magnitude at least. Then, we identified, for each b, the value of

P' (1)

such that

Δ = | \frac{φ_{o}^{α - n}}{π R_{o}^{2}} V_{o} - \int_{0}^{1} P (ξ) d ξ | \leq 10^{- 4} .

(69)

Table VI shows, for each choice of b, the value of

| P' (1) |

needed to fulfill condition

_{2}

within the approximation given by .

TABLE VI.

Maximum value of $| P' (1) |$ ⁠, for any given b, needed to fulfill condition (68) $_{2}$ within the approximation given by (69).

b	−1	−5	−10	−15	−25
$P' (1)$	$3.08 \times 10^{5}$	$6.72 \times 10^{5}$	$2.99 \times 10^{6}$	$1.56 \times 10^{7}$	$5.27 \times 10^{7}$

b	−1	−5	−10	−15	−25
$P' (1)$	$3.08 \times 10^{5}$	$6.72 \times 10^{5}$	$2.99 \times 10^{6}$	$1.56 \times 10^{7}$	$5.27 \times 10^{7}$

Table VI suggests that the larger is $| b |$ the larger is the value of $| P' (1) |$ ⁠, thus approximating the theoretical condition $| P' (1) | = + \infty$ ⁠.

It is interesting to notice that large values of b suggest the idea of a very viscous fluid, typical of a rhyolitic magma. For a rhyolitic magma, we can choose $κ = 5 \times 10^{- 7} m^{2} s^{- 1}$ (see Remark 1). Then, the corresponding emplacement time are quite similar to those estimated by Loewen et al.⁴⁸ to the case of Summit Lake at Yellowstone (USA). Indeed, this terrestrial dome is a prime example of an axisymmetric lava flow propagated over a horizontal surface. This circumstance suggests to the authors to apply the similarity solution for the case of a Newtonian fluid flow of Huppert⁶ to estimate discharge and the duration of its emplacement (see Table I in Ref. 48). However, the same authors point out that, as a lava flow cools, the yield strength of a thickening cooled crust will exert primary control on lava flow behavior, requiring application of a viscoelastic model^11,49 but also resulting in more complex compound flow. Of course, although fascinating and suggestive, this conclusion cannot be considered in any manner a conclusive indication about the largely debated question of the physical mechanism that originated the pancake domes on Venus.

In Fig. 14, we display $P (ξ)$ for various values of b and n. Note that the solution is singular at the origin, reflecting the fact that fluid is being introduced at a non-zero rate at r = 0.

FIG. 14.

View large Download slide

Left panel: the solution $P (ξ)$ of the differential equation (67) for $b = - 25, - 15, - 10, - 5$ ⁠, and $δ = 0.33, λ = 0.8$ ⁠, n = 0.3. Right panel: the solution $P (ξ)$ of the differential equation (67) for $n = 0.1, 0.3, 0.5, 0.9$ ⁠, and $δ = 0.33, λ = 0.8$ ⁠, b = −10. For lower or greater values of b, the four plots remain practically coincident.

Figure 15 shows R as a function of both temperature and time, during the effusion phase. We notice that the greater is $| b |$ (that is the smaller is the viscosity of the fluid), the faster is the dome radius to approach its final value at the end of the effusion phase. For instance, with b = −15, $R (t) = 0.95 R_{o}$ for t = 0.23.

FIG. 15.

View large Download slide

Evolution of R vs time for $b = - 1, - 5, - 10, - 15$ ⁠, n = 0.3, $δ = 0.33$ ⁠, and $λ = 0.8$ ⁠. At t = 1 the effusion phase is over, and the relaxation phase starts when the front edge is in R_o.

Figure 16 shows h as a function of temperature over the interval $(0.2, R_{o}) = (0.2, 0.65)$ ⁠, during the effusion phase, for $b = - 10, - 25$ ⁠. A representation closer to r = 0 is meaningless because of the singularity. Indeed, at the vent, the asymptotic analysis based on the thin layer approximation formally breaks down during the effusive phase.⁶

FIG. 16.

View large Download slide

Left panel: evolution of $h (r, θ)$ over the interval $(0.2, R_{o}) = (0.2, 0.65)$ as temperature decreases from 1 to $θ_{amb}$ for n = 0.3, b = −10, $δ = 0.33$ and $λ = 0.8$ ⁠. Right panel: Evolution of $h (r, θ)$ over the interval $(0.2, R_{o}) = (0.2, 0.65)$ as time increases from 0 to 1, for n = 0.3, b = −25, $δ = 0.33$ and $λ = 0.8$ ⁠.

The volume flux (discharge) at the vent during the effusion phase can be evaluated simply by differentiating V(t). It turns out

Q (t; n, λ, b, δ) = V_{o} φ_{o}^{- n} \frac{d}{d t} [φ^{n} (θ (t), λ, b, δ)] .

(70)

It must be emphasized that Q(t) is uniquely determined by $φ$ ⁠, that is, by the function which links viscosity to temperature during the evolution of the similarity solution. This suggests that the parameter b must have a non-negligible effect on Q(t). Figure 17 shows that Q(t) decreases monotonically, practically vanishing at the end of the effusion phase, consistently with the physical expectation. Indeed, for b decreasing from −15 to −25, Q(1) reduces from $O (10^{- 2})$ to $O (10^{- 6})$ ⁠.

FIG. 17.

View large Download slide

Evolution of Q over the time interval (0, 1) for n = 0.3, $b = - 10, - 15, - 25, δ = 0.33$ ⁠, and $λ = 0.8$ ⁠.

We find the dimensional discharge to be

Q (t) = \frac{H_{ref} L_{ref}^{2}}{τ t_{ref}} \int_{0}^{1} Q (t) d t = \frac{H_{ref} L_{ref}^{2}}{τ t_{ref}} V_{o} ≃ 6 \times 10^{3} m^{3} s^{- 1} .

It is interesting to note that the above value is within the range estimated by Loewen et al. in Ref. 48 analyzing the case of a point source eruption of rhyolitic lava which formed the Summit Lake, just mentioned before.

IX. CONCLUSIONS

Radar backscatter data from many steep-sided domes on Venus indicate that they have smooth, relatively unbroken surfaces. Their small aspect ratio allows applying a mathematical approach based on the thin-layer approximation. Our research does not address questions about the physical mechanism that leads to the observed shapes, such as those considered by Stofan et al.,¹⁷ but criticized by others.¹ Similarly, we do not discuss the question of which type of lava might have formed these domes.²⁷

Our primary goal is to generalize the approach of Quick et al.,¹ establishing a mathematical model for a non-Newtonian viscous fluid cooling and spreading over a rigid surface, capable of capturing the patterns observed in the SAR images and altimeter profiles. Our model, indeed, generalizes the one by Quick et al.,¹ under several respects, listed in the sequel:

Lava is considered a non-Newtonian fluid of power-law type. In particular, we adapt the approach of Rajagopal et al.,³⁶ allowing the viscosity to depend on temperature (instead of pressure). The non-Newtonian constitutive model represents a step forward with respect to the theory illustrated in Ref. 1. Indeed, the introduction of a further parameter besides viscosity allows to fit two different profiles of the same dome (e.g., SW-NE and N-S) with the same set of parameters. Figure 6 shows that, with λ = 1, the two profiles of the dome are fitted with two significantly different values of τ (duration of the relaxing phase), a circumstance difficult to be justified from the physical point of view.
As in Quick et al.,¹ we distinguish the two main phases of the dome evolution, the effusive and relaxing ones. However, unlike Ref. 1, the influence of the temperature is coherently included in the model by using the heat equation in the purely radiative regime, since a careful dimensional analysis shows that diffusion and convective effects can be neglected.
Our simulations and subsequent comparison with the available data, suggest a precise variable flow rate at the vent Q(t), based on the mass conservation law. In our model, Q(t) turns out monotonically decreasing with time. From this point of view, we are far from the approach in Ref. 1: there, the authors allow a constant flow rate obtained by applying the model of Babu and van Genuchten^50,51 which we believe not fully consistent in the present context. Indeed, the approach by Babu and van Genuchten [see Eq. (a) in Ref. 51 and compare it with Eq. (52) in the present paper] concerns an isothermal Newtonian flow taking place in radial direction into an initially dry homogeneous isotropic soil under the assumption of a constant injection rate at the origin.
The values of the dimensional flow rate provided by our model are compatible with those suggested by Ref. 48 for Summit Lake flow, a rhyolitic eruption in Yellowstone, applying the classical Newtonian model by Huppert.^6,7 Although not conclusive, this circumstance suggests a possible answer to the intriguing and largely debated question of the physical process that originated the pancake domes on Venus.
Finally, we remark that our model is stated in a fully dimensionless form, but when we introduce a suitable length scale, the calculated physical quantities appear to match significantly those present in the current literature. Nevertheless, a lot of caution should be taken with these values, since, we remind, nothing is known about the chemical composition of the lava which formed the dome considered.

ACKNOWLEDGMENTS

This work has been done under the framework of GNMF (Indam), PRIN 2022 Project “Mathematical modelling of heterogeneous systems,” and National Recovery and Resilience Plan, Mission 4 Component 2—Investment 1.4—National Center for HPC, Big Data and Quantum Computing—funded by the European Union—Next Generation EU—CUP B83C22002830001.

AUTHOR DECLARATIONS

Conflict of Interest

The authors have no conflicts to disclose.

Author Contributions

Benedetta Calusi: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Validation (equal); Writing – review & editing (equal). Angiolo Farina: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Validation (equal); Writing – review & editing (equal). Lorenzo Fusi: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Validation (equal); Writing – review & editing (equal). Fabio Rosso: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Software (lead); Validation (equal); Writing – original draft (lead); Writing – review & editing (equal).

DATA AVAILABILITY

Data sharing is not applicable to this article as no new data were created or analyzed in this study.

REFERENCES

1.

L. C.

Quick

,

L. S.

Glaze

,

S. M.

Baloga

, and

E. R.

Stofan

, “

New approaches to inferences for steep-sided domes on Venus

,”

J. Volcanol. Geotherm. Res.

319

,

93

–

105

(

2016

).

https://doi.org/10.1016/j.jvolgeores.2016.02.028

Google Scholar

Crossref

2.

N. S.

Bagdassarov

, “

Non-Newtonian and viscoelastic properties of lava flows

,”

AGU Fall Meeting Abstracts

(

NASA

,

2004

), pp.

V51D

–

V504

.

Google Scholar

3.

I.

Sonder

, “

Non-Newtonian properties of magmatic melts

,” Ph.D. thesis (

Universität Würzburg

,

2010

).

Google Scholar

4.

M.

Dragoni

, “

Modelling the rheology and cooling of lava flows

,” in

Active Lavas

(

Routledge

,

2022

), pp.

235

–

261

.

Google Scholar

Crossref

5.

I.

Sonder

,

B.

Zimanowski

, and

R.

Büttner

, “

Non-Newtonian viscosity of basaltic magma

,”

Geophys. Res. Lett.

33

(

2

),

L02303

, https://doi.org/10.1029/2005GL024240 (

2006

).

Google Scholar

Crossref

6.

H. E.

Huppert

, “

The propagation of two-dimensional and axisymmetric viscous gravity currents over a rigid horizontal surface

,”

J. Fluid Mech.

121

,

43

–

58

(

1982

).

https://doi.org/10.1017/S0022112082001797

Google Scholar

Crossref

7.

H. E.

Huppert

,

J. B.

Shepherd

,

H.

Sigurdsson

, and

R. S. J.

Sparks

, “

On lava dome growth with application to the 1979 lava extrusion of the Soufière of St.Vincent

,”

J. Volcanol. Geotherm. Res.

14

,

199

–

222

(

1982

).

https://doi.org/10.1016/0377-0273(82)90062-2

Google Scholar

Crossref

8.

H. E.

Huppert

, “

Gravity currents: A personal perspective

,”

J. Fluid Mech.

554

,

299

–

322

(

2006

).

https://doi.org/10.1017/S002211200600930X

Google Scholar

Crossref

9.

T. B.

Benjamin

, “

Gravity currents and related phenomena

,”

J. Fluid Mech.

31

,

209

–

248

(

1968

).

https://doi.org/10.1017/S0022112068000133

Google Scholar

Crossref

10.

G.

Algwauish

and

S.

Naire

, “

The thermo-viscous fingering instability of a cooling spreading liquid dome

,”

Phys. Fluids

35

(

11

),

112109

(

2023

).

https://doi.org/10.1063/5.0173902

Google Scholar

Crossref

11.

N. J.

Balmforth

and

R. V.

Craster

, “

Dynamics of cooling domes of viscoplastic fluid

,”

J. Fluid Mech.

422

,

225

–

248

(

2000

).

https://doi.org/10.1017/S002211200000166X

Google Scholar

Crossref

12.

N. J.

Balmforth

,

A. S.

Burbidge

,

R. V.

Craster

,

J.

Salzig

, and

A.

Shen

, “

Visco-plastic models of isothermal lava domes

,”

J. Fluid Mech.

403

,

37

–

65

(

2000

).

https://doi.org/10.1017/S0022112099006916

Google Scholar

Crossref

13.

N. J.

Balmforth

,

R. V.

Craster

, and

R.

Sassi

, “

Shallow viscoplastic flow on an inclined plane

,”

J. Fluid Mech.

470

,

1

–

29

(

2002

).

https://doi.org/10.1017/S0022112002001660

Google Scholar

Crossref

14.

R. W.

Griffiths

, “

The dynamics of lava flows

,”

Annu. Rev. Fluid Mech.

32

(

1

),

477

–

518

(

2000

).

https://doi.org/10.1146/annurev.fluid.32.1.477

Google Scholar

Crossref

15.

D.

Bercovici

, “

A theoretical model of cooling viscous gravity currents with temperature-dependent viscosity

,”

Geophys. Res. Lett.

21

(

12

),

1177

–

1180

, https://doi.org/10.1029/94GL01124 (

1994

).

Google Scholar

Crossref

16.

M.

Dragoni

,

I.

Borsari

, and

A.

Tallarico

, “

A model for the shape of lava flow fronts

,”

J. Geophys. Res.

110

(

B9

),

B09203

, https://doi.org/10.1029/2004JB003523 (

2005

).

Google Scholar

17.

E. R.

Stofan

,

W. S.

Anderson

,

D. A.

Crown

, and

J. J.

Plaut

, “

Emplacement and composition of steep-sided domes on Venus

,”

J. Geophys. Res.

105

(

E11

),

26757

–

26771

, https://doi.org/10.1029/1999JE001206 (

2000

).

Google Scholar

Crossref

18.

S. E. H.

Sakimoto

and

M. T.

Zuber

, “

The spreading of variable-viscosity axisymmetric radial gravity currents: Applications to the emplacement of Venusian ‘pancake’ domes

,”

J. Fluid Mech.

301

,

65

–

77

(

1995

).

https://doi.org/10.1017/S0022112095003806

Google Scholar

Crossref

19.

N. T.

Bridges

, “

Submarine analogs to Venusian pancake domes

,”

Geophys. Res. Lett.

22

(

20

),

2781

–

2784

, https://doi.org/10.1029/95GL02662 (

1995

).

Google Scholar

Crossref

20.

N. T.

Bridges

, “

Ambient effect on basalt and rhyolite lavas under Venusian, subaerial and subacqueous conditions

,”

J. Geophys. Res.

102

(

E4

),

9243

–

9255

, https://doi.org/10.1029/97JE00390 (

1997

).

Google Scholar

Crossref

21.

M. E.

Borrelli

,

J. G.

O'Rourke

,

S. E.

Smrekar

, and

C. M.

Ostberg

, “

A global survey of lithospheric flexure at steep-sided domical volcanoes on Venus reveals intermediate elastic thicknesses

,”

J. Geophys. Res.: Planets

126

(

7

),

e2020JE006756

, https://doi.org/10.1029/2020JE006756 (

2021

).

Google Scholar

Crossref

22.

L. C.

Quick

,

L. S.

Glaze

, and

S. M.

Baloga

, “

Venusian steep-sided domes: Essential exploration targets for constraining the range of volcanic emplacement conditions

,” in

Workshop on Venus Exploration Targets

(

The Venus Exploration Analysis Group

,

2014

), Vol.

1781

, p.

6014

.

Google Scholar

23.

M. A.

Ivanov

, “

Stratigraphic and geographic distribution of steep-sided domes on Venus: Preliminary results from regional

,”

J. Geophys. Res.

104

,

18907

–

18924

, https://doi.org/10.1029/1999JE001039 (

1999

).

Google Scholar

Crossref

24.

D. K.

Smith

, “

Comparison of the shapes and sizes of seafloor volcanoes on earth and ‘pancake’ domes on Venus

,”

J. Volcanol. Geotherm. Res.

73

(

1

),

47

–

64

(

1996

).

https://doi.org/10.1016/0377-0273(96)00007-8

Google Scholar

Crossref

25.

A. L.

Gleason

,

R. R.

Herrick

, and

J. M.

Byrnes

, “

Analysis of Venusian steep-sided domes utilizing stereo-derived topography

,”

J. Geophys. Res.

115

,

E06004

, http://dx.doi.org/10.1029/2009JE003431 (

2010

).

Google Scholar

26.

B.

Pavri

,

J. W.

Head

,

K. B.

Klose

, and

L.

Wilson

, “

Steep-sided domes on Venus: Characteristics, geologic setting, and eruption conditions from Magellan data

,”

J. Geophys. Res.

97

(

E8

),

13445

–

13478

, https://doi.org/10.1029/92JE01162 (

1992

).

Google Scholar

Crossref

27.

D.

McKenzie

,

P. G.

Ford

,

F.

Liu

, and

G. H.

Pettengill

, “

Pancakelike domes on Venus

,”

J. Geophys. Res.

97

(

E10

),

15967

–

15976

, https://doi.org/10.1029/92JE01349 (

1992

).

Google Scholar

Crossref

28.

C. E.

Harnett

,

M. J.

Heap

, and

M. E.

Thomas

, “

A toolbox for identifying the expression of dome-forming volcanism on exoplanets

,”

Planet. Space Sci.

180

,

104762

(

2020

).

https://doi.org/10.1016/j.pss.2019.104762

Google Scholar

Crossref

29.

E.

Cottrell

, “

Global distribution of active volcanoes

,” in

Volcanic Hazards, Risks and Disasters

(

Elsevier

,

2015

), pp.

1

–

16

.

Google Scholar

Crossref

30.

J. H.

Fink

,

N. T.

Bridges

, and

R. E.

Grimm

, “

Shapes of venusian ‘pancake’ domes imply episodic emplacement and silicic composition

,”

Geophys. Res. Lett.

20

(

4

),

261

–

264

, https://doi.org/10.1029/92GL03010 (

1993

).

Google Scholar

Crossref

31.

D.

Bercovici

and

Y.

Ricard

, “

Plate tectonics, damage and inheritance

,”

Nature

508

(

7497

),

513

–

516

(

2014

).

https://doi.org/10.1038/nature13072

Google Scholar

Crossref

PubMed

32.

S.

Mikhail

and

D. A.

Sverjensky

, “

Nitrogen speciation in upper mantle fluids and the origin of earth's nitrogen-rich atmosphere

,”

Nat. Geosci.

7

(

11

),

816

–

819

(

2014

).

https://doi.org/10.1038/ngeo2271

Google Scholar

Crossref

33.

J. W.

Head

,

D. B.

Campbell

,

C.

Elachi

,

J. E.

Guest

,

D. P.

McKenzie

,

R. S.

Saunders

,

G. G.

Schaber

, and

G.

Schubert

, “

Venus volcanism: Initial analysis from magellan data

,”

Science

252

(

5003

),

276

–

288

(

1991

).

https://doi.org/10.1126/science.252.5003.276

Google Scholar

Crossref

PubMed

34.

S.

Hensley

and

S.

Shaffer

, “

Automatic dem generation using Magellan stereo data

,” in

Proceedings of IGARSS '94 – 1994 IEEE International Geoscience and Remote Sensing Symposium

(

IEEE

,

1994

), Vol.

3

, pp.

1470

–

1472

.

Google Scholar

Crossref

35.

A. R.

McBirney

and

T.

Murase

, “

Rheological properties of magmas

,”

Annu. Rev. Earth Planet. Sci.

12

(

1

),

337

–

357

(

1984

).

https://doi.org/10.1146/annurev.ea.12.050184.002005

Google Scholar

Crossref

36.

K. R.

Rajagopal

,

G.

Saccomandi

, and

L.

Vergori

, “

Flow of fluids with pressure- and shear-dependent viscosity down an inclined plane

,”

J. Fluid Mech.

706

,

173

–

189

(

2012

).

https://doi.org/10.1017/jfm.2012.244

Google Scholar

Crossref

37.

P.

Hartlieb

,

M.

Toifl

,

F.

Kuchar

,

R.

Meisels

, and

T.

Antretter

, “

Thermo-physical properties of selected hard rocks and their relation to microwave-assisted comminution

,”

Miner. Eng.

91

,

34

–

41

(

2016

).

https://doi.org/10.1016/j.mineng.2015.11.008

Google Scholar

Crossref

38.

W. L.

Romine

,

A. G.

Whittington

,

P. I.

Nabelek

, and

A. M.

Hofmeister

, “

Thermal diffusivity of rhyolitic glasses and melts: Effects of temperature, crystals and dissolved water

,”

Bull. Volcanol.

74

,

2273

–

2287

(

2012

).

https://doi.org/10.1007/s00445-012-0661-6

Google Scholar

Crossref

39.

J. V.

Smith

, “

Shear thickening dilatancy in crystal-rich flows

,”

J. Volcanol. Geotherm. Res.

79

(

1–2

),

1

–

8

(

1997

).

https://doi.org/10.1016/S0377-0273(97)00020-6

Google Scholar

Crossref

40.

S. L.

Webb

and

D. B.

Dingwell

, “

Non-Newtonian rheology of igneous melts at high stresses and strain rates: Experimental results for rhyolite, andesite, basalt, and nephelinite

,”

J. Geophys. Res.

95

(

B10

),

15695

–

15701

, https://doi.org/10.1029/JB095iB10p15695 (

1990

).

Google Scholar

Crossref

41.

L.

Liu

and

R. P.

Lowell

, “

Modeling heat transfer from a convecting, crystallizing, replenished silicic magma chamber at an oceanic spreading center

,”

Geochem. Geophys. Geosyst.

12

(

9

),

Q09010

(

2011

).

https://doi.org/10.1029/2011GC003612

Google Scholar

Crossref

42.

M. R.

James

,

N.

Bagdassarov

,

K.

Müller

, and

H.

Pinkerton

, “

Viscoelastic behaviour of basaltic lavas

,”

J. Volcanol. Geotherm. Res.

132

(

2–3

),

99

–

113

(

2004

).

https://doi.org/10.1016/S0377-0273(03)00340-8

Google Scholar

Crossref

43.

C. E.

Lesher

and

F. J.

Spera

, “

Thermodynamic and transport properties of silicate melts and magma

,” in

The Encyclopedia of Volcanoes

(

Elsevier

,

2015

), pp.

113

–

141

.

Google Scholar

Crossref

44.

C. H.

Donaldson

, “

Calculated diffusion coefficients and the growth rate of olivine in a basalt magma

,”

Lithos

8

(

2

),

163

–

174

(

1975

).

https://doi.org/10.1016/0024-4937(75)90023-7

Google Scholar

Crossref

45.

S.

Diniega

,

S. E.

Smrekar

,

S.

Anderson

, and

E. R.

Stofan

, “

The influence of temperature-dependent viscosity on lava flow dynamics

,”

J. Geophys. Res.: Earth Surf.

118

(

3

),

1516

–

1532

(

2013

).

https://doi.org/10.1002/jgrf.20111

Google Scholar

Crossref

46.

H.

Blatt

,

R.

Tracy

, and

B.

Owens

,

Petrology: Igneous, Sedimentary and Metamorphic

, 3rd ed. (

Freeman

,

New York

,

2006

).

Google Scholar

47.

B.

Clarke

,

E. S.

Calder

,

F.

Dessalegn

,

K.

Fontijn

,

J. A.

Cortés

,

M.

Naylor

,

I.

Butler

,

W.

Hutchison

, and

G.

Yirgu

, “

Fluidal pyroclasts reveal the intensity of peralkaline rhyolite pumice cone eruptions

,”

Nat. Commun.

10

(

1

),

2010

(

2019

).

https://doi.org/10.1038/s41467-019-09947-8

Google Scholar

Crossref

PubMed

48.

M. W.

Loewen

,

IN.

Bindeman

, and

O. E.

Melnik

, “

Eruption mechanisms and short duration of large rhyolitic lava flows of Yellowstone

,”

Earth Planet. Sci. Lett.

458

,

80

–

91

(

2017

).

https://doi.org/10.1016/j.epsl.2016.10.034

Google Scholar

Crossref

49.

S.

Blake

and

B. C.

Bruno

, “

Modelling the emplacement of compound lava flows

,”

Earth Planet. Sci. Lett.

184

(

1

),

181

–

197

(

2000

).

https://doi.org/10.1016/S0012-821X(00)00278-8

Google Scholar

Crossref

50.

D. K.

Babu

and

M. T.

Van Genuchten

, “

A similarity solution to a nonlinear diffusion equation of the singular type: A uniformly valid solution by perturbations

,”

Quart. Appl. Math.

37

(

1

),

11

–

21

(

1979

).

https://doi.org/10.1090/qam/530666

Google Scholar

Crossref

51.

D. K.

Babu

and

M. T.

Van Genuchten

, “

A perturbation solution of the nonlinear Boussinesq equation: The case of constant injection into a radial aquifer

,”

J. Hydrol.

48

(

3–4

),

269

–

280

(

1980

).

https://doi.org/10.1016/0022-1694(80)90119-5

Google Scholar

Crossref

Published open access through an agreement with Università degli Studi di Firenze

Thermo-mechanical modeling of pancakelike domes on Venus

I. INTRODUCTION: GENERAL CONSIDERATIONS ABOUT TERRESTRIAL AND EXTRA-TERRESTRIAL LAVA DOMES

II. FORMULATION OF THE GENERAL MODEL

III. MODELING THE RELAXING PHASE: ISOTHERMAL AND RADIAL APPROXIMATION

IV. SOME SIMULATIONS USING BOTH SAR AND DEM ALTIMETRY DATA

V. DEALING WITH THE TIME SCALE

VI. MODELING OF THE EFFUSION PHASE: COUPLING THE FLOW WITH THE THERMAL FIELD

VII. UNIFORM TEMPERATURE APPROXIMATION

VIII. NUMERICAL SIMULATIONS

IX. CONCLUSIONS

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

Thermo-mechanical modeling of pancakelike domes on Venus Open Access

I. INTRODUCTION: GENERAL CONSIDERATIONS ABOUT TERRESTRIAL AND EXTRA-TERRESTRIAL LAVA DOMES

II. FORMULATION OF THE GENERAL MODEL

III. MODELING THE RELAXING PHASE: ISOTHERMAL AND RADIAL APPROXIMATION

IV. SOME SIMULATIONS USING BOTH SAR AND DEM ALTIMETRY DATA

V. DEALING WITH THE TIME SCALE

VI. MODELING OF THE EFFUSION PHASE: COUPLING THE FLOW WITH THE THERMAL FIELD

VII. UNIFORM TEMPERATURE APPROXIMATION

VIII. NUMERICAL SIMULATIONS

IX. CONCLUSIONS

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Related Content

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only

Thermo-mechanical modeling of pancakelike domes on Venus