Temporal Evolution of Cooling by Natural Convection in an Enclosed Magma Chamber

Abstract

This research numerically studies the transient cooling of partially liquid magma by natural convection in an enclosed magma chamber. The mathematical model is based on the conservation laws for momentum, energy and mass for a non-Newtonian and incompressible fluid that may be modeled by the power law and the Oberbeck–Boussinesq equations (for basaltic magma) and solved with the finite volume method (FVM). The results of the programmed algorithm are compared with those in the literature for a non-Newtonian fluid with high apparent viscosity (10–200 Pa s) and Prandtl (Pr = 4 × 10⁴) and Rayleigh (Ra = 1 × 10⁶) numbers yielding a low relative error of 0.11. The times for cooling the center of the chamber from 1498 to 1448 K are 40 ky (kilo years), 37 and 28 ky for rectangular, hybrid and quasi-elliptical shapes, respectively. Results show that for the cases studied, natural convection moved the magma but had no influence on the isotherms; therefore the main mechanism of cooling is conduction. When a basaltic magma intrudes a chamber with rhyolitic magma in our model, natural convection is not sufficient to effectively mix the two magmas to produce an intermediate SiO₂ composition.

1. Introduction

To understand the processes that occur in magma chambers it is necessary to use multidisciplinary knowledge from geology, fluid dynamics and engineering [1]. Rheology governs processes such as the rise of magma, the mixing of various batches, the settling of crystals and the transport of xenoliths [2]. Magmas have a range of temperatures and chemical compositions that affect fluid convection and elemental diffusion [3]. Basaltic magma is usually treated as a quasi-Newtonian fluid, but several laboratory experiments have shown that, even for low viscosity, quasi-Newtonian behavior cannot be applied [4]. Instead, basaltic magma behaves like a Bingham fluid [5,6] or a power law fluid [4]. A Bingham fluid behaves like a rigid body when subjected to low stresses but flows like a viscous fluid under higher stresses. In this paper, we include a non-Newtonian model for the magma.

Modeling convective cooling in a magma chamber must take into account an estimate of its size and shape, as well as other geometric, thermodynamic, thermophysical and compositional factors that influence the history of pluton cooling [7,8]. Rectangular shapes approximate a sill-like intrusion or a laccolith. The oblate ellipsoid is a common form for many large chambers [9]. In this paper we study both simple rectangular and elliptical shapes.

Convection in a magma chamber [10,11] due to heating from below and cooling from above is a special case of natural convection called Rayleigh–Bernard (RB) convection [12]. For Rayleigh values higher than 1703, the conditions are thermally unstable and advection occurs. For Rayleigh numbers above 5 × 10⁴, fluid motion evolves through many different patterns before becoming turbulent [13]. Correctly calculating a Rayleigh number based on the temperature (Ra_T) due to thermal coupling between the magma and host is complicated. A more realistic method uses the flux-based “magmatic Rayleigh number” (Ra_F) [14] for natural magma convection in a chamber. In this work the Ra_F is used.

A typical Rayleigh number for real chambers may be in the range of 1 × 10⁹–1 × 10¹⁷, or 1 × 10⁶–1 × 10⁹ [15] for smaller chambers. Magmas have high Prandtl numbers of 1 × 10³–1 × 10⁸ [16], implying that the momentum diffusion is much greater than thermal diffusion. Magmas have an even higher Lewis number of 1 × 10⁴–1 × 10¹³ [17], which means that heat flow is much greater than chemical molecular diffusion. In this work, these observations about the dimensionless numbers will be very important in analyzing the simulation results.

This work studies the cooling of basaltic magma by RB natural convection inside a closed chamber located in the Peruvian Coastal Batholith (PCB). The novelty of this study has two main aspects: (1) the development of a mathematical model of RB convection that includes experimental data for rheology of basaltic magma as it just begins to crystallize and (2) computational simulations of magma chambers that take into account geological aspects of the Humay unit of the PCB. The study describes the magma flow pattern and temperature distribution inside the chamber when the crystal volume in the magma is low. In this work we use a power law non-Newtonian model. The constants for our mathematical model of non-Newtonian viscosity are obtained from experimental data for magma cooling over the range 1498–1448 K [4]. Calculations at ambient pressure indicate that the liquidus and solidus temperature of basaltic magma are approximately 1503 and 1223 K respectively [18]; therefore, if the magma is in equilibrium, it has already begun to solidify at 1498 K. The temperature range of 1498–1448 K is in a higher range where the magma is mostly liquid and has just begun to crystalize. The mathematical model includes the interaction of the magma with the host rock where only conduction is considered. The simulations were run until the temperature at a central point located within the chamber dropped to 1448 K.

The geological aspects of the simulation are considered in detail in Section 2.1. In Section 3.2.1 and Section 3.2.2, the cooling time of magma as a function of the shape of the chamber is studied. In Section 3.2.3, a hypothetical case where a basaltic magma intruded into a rhyolitic magma is also simulated.

Our mathematical model is focused on solving the movement of the magma and does not include the phase change enthalpy in the energy equation. The phase change enthalpy normally is included in mathematical models as an additional term in the energy equation, which is a function of phase change [15,19,20]. Cooling times that do not include the phase change enthalpy are always shorter so the cases solved in this work yield the minimum cooling time for liquid in a magma chamber. These considerations mean that there are two assumptions made on which our analysis are based: (1) the streamlines and the velocities with which the magma moves are always the highest that can be found for each chamber studied and therefore, (2) since the phase change function and the enthalpy of solidification are not considered, the cooling times for each chamber studied are the shortest. This study does not consider water or vapor advection in the host rock which could decrease the cooling times of the magma.

2. Materials and Methods

2.1. The Physical and Geological Situation

Three chamber shapes are simulated: rectangular (case 1), hybrid (case 2) and quasi-elliptical (case 3) [21]. In all cases, the chamber has the same volume of 18 km³, but varying widths. Figure 1 displays the magma chamber for the most complex case as suggested by previous work in the PCB (Humay unit) [22,23]. Geological data indicate that PCB plutons are commonly 5 km thick or less [24,25]. Analysis of saline inclusions in the PCB Linga complex suggests magma crystallization pressures between 800 and 900 bar, indicating that the depth to the top of the magma chamber was between 2.9 and 3.2 km from the Earth’s surface when crystallizing [22]. In this work we locate the top of the chamber at 3 km from the surface and use a maximum thickness for the chamber of 3 km (L_c = 3000 m).

The basaltic magma inside the chamber is initially stationary with initial temperature of $T_{f,0}=1498 \mathrm{K}$.

This paper uses an initial temperature gradient in the host rock of 25 K/km [26] which corresponds to T_s,0 = 298 K at the upper boundary (surface) and T_s,0,b = 448 K at the lower boundary. The lower edge of the computational domain (chamber and host rock) is considered adiabatic. This could occur if the cooling below the bottom boundary of the model together with the cooling of the modeled magma chamber would lead to a small temperature gradient near the bottom boundary.

The host rock for the PCB is primarily basaltic andesite with SiO₂ silica concentration of 50–55 wt%; thus, C_s = 0.5 is used in all simulations. It is common for basaltic magma to have an SiO₂ concentration of 45–50 wt% and rhyolitic magma of 65–75 wt% SiO₂ [18]. Section 3.2.3 presents the simulation results of a chamber with two magma compositions. Initially SiO₂ in the lower zone of the chamber is $C_{f,0,b}=0.47$ and in the upper zone $C_{f,0,a}=0.7$, corresponding to basaltic magma intruding rhyolitic magma. This means that at the beginning of the intrusion simulation the chamber contains mostly rhyolitic magma, but with basaltic magma at L_intr = 1000 m at the bottom, with all at the same temperature.

2.2. Mathematical Model

A mathematical model for transient two-dimensional natural convection for non-Newtonian fluid flow is now presented. The relation between density and pressure is simplified by the so-called Oberbeck–Boussinesq approximation [27] which means that density can only change with temperature and hence flow is incompressible, assuming no volatiles are present. The mathematical model assumes that the crystals in the magma follow the flow path of the fluid and there is no settling. The magma stops moving when there is no more convection, which occurs when it reaches around 40% crystals (mush state). The system of equations presented corresponds to the one used generally for laminar flow and is suggested for RB convection for Ra < 1 × 10¹⁴ by several authors [15,27]. The partial differential equations (PDEs) for mass conservation and momentum for a melt are:

$$\frac{\partial u_f}{\partial x}+\frac{\partial v_f}{\partial y}=0$$

(1)

$${\rho }_f\left[\frac{\partial u_f}{\partial t}+u_f\frac{\partial u_f}{\partial x}+v_f\frac{\partial u_f}{\partial y}\right]=-\frac{\partial p_f}{\partial x}+\left[\frac{\partial {\sigma }_{xx}}{\partial x}+\frac{\partial {\sigma }_{yx}}{\partial y}\right]$$

(2)

$${\rho }_f\left[\frac{\partial v_f}{\partial t}+u_f\frac{\partial v_f}{\partial x}+v_f\frac{\partial v_f}{\partial y}\right]=-\frac{\partial p_f}{\partial y}+\left[\frac{\partial {\sigma }_{xy}}{\partial x}+\frac{\partial {\sigma }_{yy}}{\partial y}\right]+{\rho }_{f}g{\beta }_f\Delta T_{ref}$$

(3)

where ${\rho }_f$ is the density of the reference magma, $g$ is gravity and ${\beta }_f$ is the volumetric expansion coefficient. The $\Delta T_{ref}$ is the reference temperature variation which is related to the average magma temperature in the y direction. This difference is calculated numerically with the equation, $\Delta T_{ref}={\displaystyle \sum _{j=1}^{ny}{(\Delta T/\Delta y)}_j}$, where ny is the total node number for the magma in the y direction. In these Equations (1)–(3) and the following equations, the subscript f is used for the fluid. The energy equation is given by:

$${\rho }_{f}C{p}_f\left[\frac{\partial T_f}{\partial t}+u_f\frac{\partial T_f}{\partial x}+v_f\frac{\partial T_f}{\partial y}\right]=k_f\left[\frac{{\partial }^2{T}_f}{\partial x^2}+\frac{{\partial }^2{T}_f}{\partial y^2}\right]$$

(4)

The equation for mass concentration of SiO₂ in magma is:

$$\frac{\partial C_f}{\partial t}+u_f\frac{\partial C_f}{\partial x}+v_f\frac{\partial C_f}{\partial y}=D\left[\frac{{\partial }^2{C}_f}{\partial x^2}+\frac{{\partial }^2{C}_f}{\partial y^2}\right]$$

(5)

In simulations, these equations are solved by coupling them together. This means that in each time step the fluid velocity is used to calculate the temperature, which in turn is used to calculate the non-Newtonian viscosity and the velocity in the next iteration until the convergence criterion imposed for the time step is fulfilled. Section 2.3 gives more details of the procedure to solve the governing equations for the magma.

The thermal diffusion equation is used to calculate the conduction in the host rock:

$${\rho }_{s}C{p}_s\frac{\partial T_s}{\partial t}=k_s\left[\frac{{\partial }^2{T}_s}{\partial x^2}+\frac{{\partial }^2{T}_s}{\partial y^2}\right]$$

(6)

The subscript s is used for the solid and indicates properties of the host rock. The silica content SiO₂ of the host rock is considered to be constant at C_s = 0.5.

According to the power law model, the equations for the shear stress and shear rate for a non-Newtonian liquid are [28]:

$${\sigma }_{xx}=\eta \left[2\frac{\partial u}{\partial x}\right] ,{\sigma }_{yy}=\eta \left[2\frac{\partial v}{\partial y}\right], {\sigma }_{xy}={\sigma }_{yx}=\eta \left[\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right]$$

(7)

$$\dot{\gamma }={\left\{2\left[{\left(\frac{\partial u}{\partial x}\right)}^2+{\left(\frac{\partial v}{\partial y}\right)}^2\right]+{\left[\frac{\partial v}{\partial x_{}}+\frac{\partial u}{\partial y_{}}\right]}^2\right\}}^{0.5}$$

(8)

Apparent viscosity $\eta$ was defined by the Oswald de Waele power law in terms of the power law index n [29]:

$$\eta =K\times {\dot{\gamma }}^{n-1}Exp\left[\frac{n\times \Delta E(T_f)}{R\times T_f}\right]$$

(9)

where K is a constant called the consistency index and $\Delta E(T_f)$ is the activation energy defined by a temperature function calculated as follows:

$$\Delta E(T_f)=E_0+a\times T_f$$

(10)

The constants E₀ and a were calculated from experimental data for basaltic magma at four temperatures: 1498, 1473, 1458 and 1448 K [4]. The resulting values used in the non-Newtonian model of Equations (9) and (10) for cooling of basaltic magma are: E₀ = 1.5 × 10⁶ J/mol, a = −994.4827 J/mol K, K = 15 Pa s and n = 0.58. These parameters used in the range 1498–1448 K are constants for the non-Newtonian power law model and do not depend on either time or temperature. This means that the effects of crystals in the movement of the basaltic magma are included.

The velocity scale in the fluid, which coincides with the maximum velocity of the magma, is defined as:

$$U_0=\sqrt{g\times \beta \times L_{cl}\times \Delta T_{cl}}$$

(11)

where $\Delta T_{cl}$ is the maximum temperature variation which is related to thermal boundary layer L_cl (see Equation (22)), g is gravity and $\beta$ (1/K) is the volumetric expansion coefficient. Independent and dependent dimensionless variables are defined as:

$$X=\frac{x}{L_c};Y=\frac{y}{L_c};\tau =\frac{t\times U_0}{L_c};U=\frac{u_f}{U_0};V=\frac{v_f}{U_0};P=\frac{p_f\times L_c}{{\eta }_{ref}\times U_0}$$

(12)

$$\begin{array}{l}\theta =\frac{T-T_{f,0}}{T_{s,0}-T_{f,0}}; {\theta }_m=\frac{\Delta T_{ref}}{\Delta T_{cl}}; {\alpha }^{*}=\frac{{\alpha }_f}{{\alpha }_s}; {\eta }^{*}=\frac{\eta }{{\eta }_{ref}}; {\dot{\gamma }}^{*}=\frac{L_c}{U_0}\dot{\gamma }; k^{*}=\frac{k_f}{k_s};\\ C^{*}=\frac{C_f-C_{f,0}}{C_{s,0}-C_{f,0}}; L^{*}=\sqrt{\frac{L_c}{L_{cl}}}; \Delta T_{ref}=(T_f-T_{ref})\end{array}$$

(13)

The dimensionless Prandtl, Lewis and Rayleigh numbers (magmatic) are calculated as follows:

$$\mathrm{Pr}=\frac{{\eta }_{ref}C{p}_f}{k_f};Le=\frac{{\alpha }_f}{D};Ra=\frac{{\rho }_{f}g{\beta }_f{L}_c^3\Delta T_{cl}}{{\eta }_{ref}{\alpha }_f}$$

(14)

The dimensionless PDE system for the conservation of mass, momentum, energy and SiO₂ are:

$$\frac{\partial U_f}{\partial X}+\frac{\partial V_f}{\partial Y}=0$$

(15)

$$\sqrt{\frac{Ra}{\mathrm{Pr}}}\left[\frac{\partial U_f}{\partial \tau }+U_f\frac{\partial U_f}{\partial X}+V_f\frac{\partial U_f}{\partial Y}\right]=-L^{*}\frac{\partial P_f}{\partial X}+{\eta }^{*}{L}^{*}\left[\frac{{\partial }^2{U}_f}{\partial X^2}+\frac{{\partial }^2{U}_f}{\partial Y^2}\right]$$

(16)

$$\sqrt{\frac{Ra}{\mathrm{Pr}}}\left[\frac{\partial V_f}{\partial \tau }+U_f\frac{\partial V_f}{\partial X}+V_f\frac{\partial V_f}{\partial Y}\right]=-L^{*}\frac{\partial P_f}{\partial Y}+{\eta }^{*}{L}^{*}\left[\frac{{\partial }^2{V}_f}{\partial X^2}+\frac{{\partial }^2{V}_f}{\partial Y^2}\right]+\sqrt{\frac{Ra}{\mathrm{Pr}}}\times L^{*2}\times {\theta }_m$$

(17)

$$\sqrt{Ra\times \mathrm{Pr}}\left[\frac{\partial {\theta }_f}{\partial \tau }+U_f\frac{\partial {\theta }_f}{\partial X}+V_f\frac{\partial {\theta }_f}{\partial Y}\right]=L^{*}\left[\frac{{\partial }^2{\theta }_f}{\partial X^2}+\frac{{\partial }^2{\theta }_f}{\partial Y^2}\right]$$

(18)

$$Le\sqrt{Ra\times \mathrm{Pr}}\left[\frac{\partial C^{*}{}_f}{\partial \tau }+U_f\frac{\partial C^{*}{}_f}{\partial X}+V_f\frac{\partial C^{*}{}_f}{\partial Y}\right]=L^{*}\left[\frac{{\partial }^2{C}^{*}{}_f}{\partial X^2}+\frac{{\partial }^2{C}^{*}{}_f}{\partial Y^2}\right]$$

(19)

The dimensionless transient heat diffusion equation for the container (host rock) is:

$${\alpha }^{*}\sqrt{Ra\times \mathrm{Pr}}\frac{\partial {\theta }_s}{\partial \tau }=L^{*}\left[\frac{{\partial }^2{\theta }_s}{\partial X^2}+\frac{{\partial }^2{\theta }_s}{\partial Y^2}\right]$$

(20)

In order to find the thickness of the thermal boundary layer L_cl and the magmatic Rayleigh number Ra, the method proposed by [14] for cooling lava lakes was used. This method assumes that the host regime determines the total flux of the host rock/chamber system. The flux-based Rayleigh number Ra_F is related to the temperature-based version Ra_T by the Nusselt number Nu, as follows:

$$R{a}_F=Nu\times R{a}_T; R{a}_F=\frac{\beta g{L}_c{}^{4}F}{{\alpha }_f{\upsilon }_f{k}_s}; F=\frac{Nu{k}_s\Delta T_c}{L_c}; {\upsilon }_f=\frac{{\eta }_{ref}}{{\rho }_f}$$

(21)

The calculated heat flux F is related to the total flux by conduction between the upper and lower surfaces of the magma chamber. The Nusselt number Nu can be used to relate the total heat flux by conduction and convection in the host rock/chamber system. This dimensionless number is the ratio of the total heat flux F across a convective system characterized by a particular temperature difference ($\Delta T_c$) to the heat flux that would be realized by conduction assuming the same temperature difference. The Nusselt number and the flux-based Rayleigh number Ra_F can be related by an empirical correlation for the RB case ($Nu=0.374R{a}_F^{1/4}$) [14]. The earlier equations allow one to obtain the total thermal boundary layer thickness L_cl, their temperature difference $\Delta T_{cl}$ and the magmatic Rayleigh number Ra.

$$L_{cl}=\frac{L_c}{Nu}; \Delta T_{cl}= \frac{F\times L_{cl}}{k_s}; Ra=\frac{{\rho }_{f}g{\beta }_f\Delta T_{cl}}{{\eta }_{ref}}\frac{L_c^3}{{\alpha }_f}$$

(22)

The properties and parameters used in the simulations are listed in

Table 1. Magma and host rock properties and parameters used in the mathematical model.

Transport Properties	Magma [18]	Host Rock [30]	Param.		Param.
$\rho ({\mathrm{kg}/\mathrm{m}}^3)$	2600	2670	$R{a}_F$	1.62 × 1017	Pr	2.42 × 105
$C_p (\mathrm{J}/\mathrm{kg} \mathrm{K})$	1450	1000	$\Delta T_c (\mathrm{K})$	75	Le	5.3 × 104
$k (\mathrm{W}/\mathrm{m} \mathrm{K})$	0.6	2.65	$Nu$	7426	Tf,0 (K)	1498
$\alpha ({\mathrm{m}}^2/\mathrm{s})$	1.6 × 10−7	9.93 × 10−7	$L_{cl} (\mathrm{m})$	0.4	Tref (K)	1473
$\beta (1/\mathrm{K})$ *	5 × 10−5	-	$\Delta T_{cl} (\mathrm{K})$	1.01 × 10−2	Cs,0	0.5
${\eta }_{ref} (\mathrm{Pa} \mathrm{s})$	100	-	$Ra$	2.19 × 1013	Cf,0,a	0.47
$D ({\mathrm{m}}^2/\mathrm{s})$	3 × 10−12	-	$U_0 (\mathrm{m}/\mathrm{s})$	1.41 × 10−3	Cf,0,b	0.7

* Obtained from [31].

. It is interesting to note the close match between the values for maximum velocity calculated in this work ($U_0$ = 1.4 × 10⁻³ m/s) and those calculated by the correlation proposed by [17] for magma chambers heated from below ($U_0$ = 0.9 × 10⁻³ m/s).

2.3. Numerical Procedure

The coupled PDEs system is solved numerically with the classical finite volume method proposed originally by [32] and used widely [33,34,35]. The algorithm programmed in FORTRAN 77 computer code solved each governing equation which was treated in the generalized form (23). This equation contains non-steady state, convection, diffusion and linearized source terms:

$$\frac{\partial (\rho \phi )}{\partial t}+div\left(\rho \overrightarrow{u}\phi \right)=div\left(\Gamma \times grad\phi \right)+Sc+Sp\times \phi$$

(23)

The integration over time is performed with an explicit Euler scheme.

$$\frac{\partial \phi }{\partial \tau }=\frac{{\phi }_{}^{\tau +\Delta \tau }-{\phi }_{}^{\tau }}{\Delta \tau }$$

(24)

At each time step $\Delta \tau$, the system of discretized nodal equations for each main variable (φ = U, V, $\theta$, C*) is solved iteratively (internal iterations) with a combination of the alternating tridiagonal matrix algorithm (TDMA) and the Gauss–Seidel method. The coupling of these main variables (external iterations) uses the SIMPLE (Semi-Implicit Linked Equations) method. In order to improve the convergence, this procedure uses under-relaxation coefficients that convert the implicit in a semi-implicit method. For the cases studied in this paper, relaxation coefficients of 0.25 are adopted for the velocity components U and V, while 0.8 is adopted for pressure. For temperature and concentration, it is not necessary to relax variables. For more details about the finite volume method and the SIMPLE algorithm, consult [36].

Convergence criterion for the residuals of the mass conservation equation of 1 × 10⁻⁴ was used in all simulations. The convergence criteria for the φ variables for the internal iteration k, in the control volume (i,j), in eachdimensionless time step ($\Delta \tau$ = 10), are:

$$\left|{\phi }_{i,j}^k-{\phi }_{i,j}^{k-1}\right|\le {10}^{-4}$$

(25)

The calculations in the first (external) iteration were carried out using Newtonian behavior in order to find a velocity field that would allow the calculation of the apparent viscosity, find its average value and consequently adjust the pressure gradient, as proposed by [32].

$${\frac{\partial p}{\partial y}|}_{\tau }={\frac{\partial p}{\partial y}|}_{\tau -1}\frac{\overline{\eta }}{{\eta }_{ref}}$$

(26)

The programmed algorithm required that when the control volume had the thermal properties of a solid, the viscosity was infinite (1 × 10³⁰ Pa s) and the velocities were zero (m/s). No convection is assumed in the volume occupied by the host rock nor in the solidified portion of melt below the solidus temperature (1223 K), and diffusion becomes the only transport mechanism.

In all three cases studied, the simulation ends when the control volume in position x = 12,000 m and y = −4500 m (Pm) meets the condition, $T_f\le 1448 \mathrm{K}$. In order to compare reference and simulations results the following equation for the relative error (RE) was used [37]:

$$RE=\frac{\left|{\displaystyle \underset{y_p=0}{\overset{y_p=\max }{\int }}f{(x_p,y_p)}_{\tau ,FVM}d{y}_p}-{\displaystyle \underset{y_p=0}{\overset{y_p=\max }{\int }}f{(x_p,y_p)}_{\tau ,ref}d{y}_p}\right|}{\left|{\displaystyle \underset{y_p=0}{\overset{y_p=\max }{\int }}f{(x_p,y_p)}_{\tau ,ref}d{y}_p}\right|}$$

(27)

In Equation (27), x_p and y_p are the abscissa and ordinate of the graphs used to obtain the RE.

The calculation of the relative error is carried out by numerical integration (Equation (27)) on the functions obtained by numerical interpolation with the discrete values from the finite volumes algorithm. This is a more precise and rigorous method for calculating the REs between results when compared to the one that is generally used and that will consider only a discrete number of comparison points. In this work, REs were up to 30% higher when the results of Equation (27) were compared with those of the classic RE equation that uses discrete points without integration.

2.4. Mesh Study

In order to optimize the calculation times with a mesh that had no effect on the results, a mesh study with four meshes was carried out. The meshes used and the results are shown in

Table 2. Different meshes studied and the relative errors (RE) between V velocity profiles in position V (2.8–5.2; 0.5), compared with the more refined mesh 1.

Mesh	Nodes X × Y	Total Nodes	N° Nodes in Chamber Zone	N° Nodes in Host Rock	RE	t (s)
1	300 × 90	27,000	2800	24,200	0	13,180
2	280 × 70	19,600	2450	17,150	0.0078	9836
3	196 × 58	11,368	2450	8918	0.1283	4393
4	156 × 60	9360	1800	7560	1.1756	3440
5	130 × 50	6500	1250	5250	2.0933	2960

.

All meshes are uniform except mesh 3, which was refined in the chamber area. The physical situation corresponds to case 1 with a rectangular chamber. Some boundary conditions to reach the steady state were simplified by using a constant dimensionless temperature of $\theta =-1$ in the host rock, an initial temperature of $\theta =0$ in the magma chamber and a constant time step of $\Delta \tau$ = 1 × 10⁻². The mathematical model assumed $L_{cl}=L_c$ and $\Delta T_{cl}= 1$. Since the driving force for natural convection is gravity acting on the Y coordinate, it is assumed that the velocity gradients would be greater in the V component of the velocity and hence the mesh convergence is studied analyzing the V profile once all variables ($\phi =U,V,T$) have converged. This is verified after $\tau =7\times {10}^4$. To calculate the RE, the dimensionless profiles in position V (2.8–5.2; 0.5) obtained with each mesh were compared with those of the denser mesh (mesh 1). The results shown in Table 2 indicate clear convergence when the mesh is denser. Mesh 5 with 6500 nodes distributed over the entire computational domain has a RE of 2.0933. By increasing the nodes to 9360 (mesh 4) the RE decreases to 1.1756, which is approximately 10 times smaller than results found for mesh 5. The greatest gradients occur in the chamber where there is natural convection of the magma, while in the host rock the diffusion produces smaller gradients. This leads us to formulate a mesh with a higher number of nodes in the chamber zone and a lower number of nodes in the host rock. Mesh 3 keeps the same quantity of nodes from mesh 2 in the chamber area, but decreases its number in the host rock zone. The RE is only 0.1283 even when it has about one third of the total nodes of mesh 1. The same happens with the computational time which decreases from 13,180 to 4393 s. Hence, a non-uniform mesh (mesh 3) of 196 × 58 was used in the simulations.

3. Results

3.1. Validation

In order to validate our computer code, we here reproduced the results obtained from [38]. Transient natural convection cooling (Ra_T = 1 × 10⁶) of a highly non-Newtonian fluid with a high Prandtl number (4 × 10⁴) inside a square cavity surrounded by a solid was studied numerically. The calculated domain was 0.2 m in x and y directions that included a Plexiglas contour of 0.01 m thickness. Initially the fluid was static (U = V = 0) and both fluid temperatures and Plexiglas were constant at T₀ = 45 °C ($\theta =0$). In order to compare published results with those obtained using our algorithm, two cases were used. In the first case a constant temperature of T_w = 25 °C ($\theta =-1$) outside the container was used. In the second case the same temperature T_w = 25 °C ($\theta =-1$) on three boundaries was imposed, but an adiabatic condition at the bottom boundary was used. The dimensionless mathematical model introduced in Section 2.2 was used to simulate the mathematical model (using $L_{cl}=L_c$ and $\Delta T_{cl}=1$) and FVM proposed in this work. In both cases, the calculated domain was discretized with five uniform meshes of 40 × 40, 180 × 180, 120 × 120, 60 × 60 and 40 × 40 nodes and three-time steps $\Delta \tau$, of 1 × 10⁻⁴, 1 × 10⁻³ and 1 × 10⁻² (in dimensionless time). Simulations were carried out until reaching 1785s ($\tau$ = 1165). To calculate the apparent viscosity, the following equation was used: 10⁻⁴

$$\begin{array}{l}\eta ={\eta }_{ref}\exp \left(\frac{1}{A{T}^2+BT+C}\right); {\eta }_{ref} =4.485E-8 \mathrm{P}\mathrm{a} \mathrm{s}; \\ A=-7.5907\times {10}^{-7}; B=3.8968\times {10}^{-4}; C=4.0130\times {10}^{-2}\end{array}$$

(28)

where the temperature must be entered in °C. The apparent viscosity varied from approximately 10 to 200 Pa s for a cooling range of 45–25 °C. These viscosity values are close to those of basaltic magma at temperatures between 1498 and 1473 K [4] and a Prandtl number of 4 × 10⁴ is close to the basaltic magma value [39]. It is important that the Prandtl number is close to that of basaltic magma since the theory of similarity between dimensionless numbers allows prediction of the same behavior of both fluids in terms of the amount of movement and heat transfer. This is very important for the validation of magmatic flows since the high temperatures at which they occur make their study very difficult in the laboratory. The two previous characteristics of the reference problem (viscosity and Pr number) allow us to confirm that the dimensionless solution for speeds and temperatures is comparable to those that occur in basaltic magmas. In Figure 2 velocity vectors and isotherms from [38] are compared to those obtained in the present work after 1785 s using a mesh with 180 × 180 nodes.

Both the streamlines and the temperature values behaved similar to those in [38] in the two cases compared. The effects of the buoyancy term on current lines and temperatures can be clearly seen. Figure 3 displays the dimensionless average temperatures of the fluid as a function of time for five meshes for the first case study in this section. The results using less refined meshes show RE above 0.16. There is a tendency toward a single result that coincides with that obtained by [38] when the mesh is refined.

The calculated average values with 200 × 200 and 180 × 180 nodes follow the trend obtained from the reference, but there are small deviations which increase slightly at the end of the calculation period. The RE of 0.057 which is the same for both meshes indicates a very good approximation and that the algorithm has converged to its final solution.

The code used to obtain results for the first case in this section can be obtained from the Supplementary Materials (https://bclausen.net/Zambra/Program_Validation.txt. Last access: 20 December 2021).

3.2. Results and Discussion

Three chamber shapes with the same 18 km³ volume are simulated (where the third dimension z is considered to be of unitary magnitude). The cases studied are: rectangular (case 1), hybrid (case 2) and quasi-elliptical (case 3). The last two could have been formed by geological phenomena that have been detailed by [21]. In all cases the simulation was carried out until the temperature at point Pm (12,000 m; −4500 m) was lower than 1448 K. This point is in the center of the rectangular chamber in case 1.

Code used to obtain results for case 3 can be obtained from the Supplementary Materials (https://bclausen.net/Zambra/Program_Chamber_case3.txt. Last access: 20 December 2021).

3.2.1. Temperature Results

Figure 4 shows the temperature distribution in the host rock and within the chamber at two different times (5 and 22 ky) for three cases. Temperatures decrease from the center to the edges of the chamber except at the bottom edge. The shape of the chamber directly affects the rate of heat loss. In case 1, the heat is retained better and concentrated in the lower central zone. In contrast, case 3 (with a quasi-elliptical shape) enhances heat elimination. This is due to the extreme lateral zones of the chamber where temperatures decrease rapidly (see case 3, t = 22 ky).

Up to 5 ky of cooling, the isotherms in the host rock follow the contour of the chamber. After 22 ky of cooling, the isotherms lose this characteristic and tend to show symmetry. After 22 ky the high temperature isotherms (above 1016 K) are greater in cases 1 and 2, which produce higher average temperatures inside the chamber as presented in Figure 5.

It can be seen that the values of the isotherms decrease regularly, in layers from the edges with lower temperature and towards the center. This occurs on all edges except the bottom edge, where it is adiabatic. This condition is typical of conduction-driven processes. Some examples can be seen in the isotherms of the wall in the cooling of Golden Syrup [38], in a mold in the case of solidification of an aluminum alloy with an adiabatic edge [40] or in the fluid within a cavity when the Rayleigh number is low [41]. The magma is gravitationally stable and the main mechanism of thermal energy transport is conduction. This behavior can be observed throughout the entire simulation time. One of the reasons for this behavior is the high magma apparent viscosity of 9 to 1.8 × 10³ Pa s according to Equations (9) and (10) for the range between 1498 to 1448 K, which produce very low velocities (see Figure 6 and Figure 7) by natural convection in the chamber. Laboratory studies have shown that natural convection may not have visible effects on isotherms when a cavity is heated from the bottom and cooled from the top [41]. The results of the simulations indicate that the velocities produced by natural convection in the magma are not high enough to alter the isotherms. Conduction between the magma and host rock produces a rapid decrease in the temperature at the edges of the chamber. In general, cooling rates are minimized when the temperature difference between magma and host is small [42].

Figure 5 shows (a) the average temperature variation with time inside the chamber until the temperature drops below 1448 K at point Pm and (b) the temperature profile at x = 12,000 m for 0.1 ky and 22 ky. Figure 5a shows that average temperatures for case 1 decrease the slowest and case 3 cools down the fastest. The average temperatures after 40, 37 and 28 ky were 1215, 1156 and 1080 K for cases 1, 2 and 3, respectively. Cooling times are related to heat flow at the edges in contact with the host rock, which depends on the shape of the chamber. These contact areas are 18, 22.3 and 28.2 km² for cases 1, 2 and 3, respectively. This in itself could explain why the cooling times in Figure 5a are different, but there is another reason. If the same boundary conditions were applied, similar heat losses should occur. However, the chambers have different shapes and this causes the edges of the chamber to come in contact with different temperatures of the host rock where the temperature increases by 25 K/km of depth. Fourier’s law indicates that the greatest heat flux at the host rock/chamber interface is found where there is the highest temperature gradient. In our case, this occurs when the host rock has lower temperatures (where the magma has a uniform initial temperature of T_f,0) at a shallower depth. Case 1 has less area exposed to high gradients, while cases 2 and 3 have larger exposed areas where the temperature gradient is high. These two reasons explain the large differences in temperature between the shapes studied and highlight the effect of chamber shape on the magma’s cooling times. Small batches of magma moving from hotter to colder regions in the upper crust can achieve thermal equilibrium in 1 ky, while larger igneous systems with spatial scales of 10–100 km may have lifetimes of 100 to 10,000 ky [43].

Recent field and laboratory studies for the mafic Humay pluton of the PCB suggest that cooling times for the crystallization of zircons (~1200–1100 K) may be on the order of a million years or more [44]. For comparison to simulations in this section, cooling times at the center of the pluton during the fully liquid stage (1498–1448 K) are 20–40 ky. When the simulations were run until the temperature dropped to 1173 K, the cooling times increased to 419–460 ky depending on the shape of the chamber. The mathematical model for simulations in this study gives a relatively short cooling time, since it does not consider the enthalpy of phase change. A model closer to reality, but not carried out in this work, would incorporate the phase change, thus considerably increasing the cooling time closer to that presented by [44]. Additionally, the Humay pluton could be up to 5 km thick [23], instead of 3 km simulated here, which would also increase the estimated cooling time.

In general, the temperature profiles in Figure 5b show similar variations for the three cases and the two times, except for case 3 (near −6000 m) where the profile at x = 12,000 m is slightly displaced from the area where the chamber has its greatest depth (see Figure 4, case 3). After 22 ky the temperatures inside the chamber first drop in the extreme upper zone to values around 1000 K. The temperatures of the host rock have increased significantly to values around 375 K at a depth of one km from the surface. This is a little more than 50 K higher than the original temperature due to the geothermal depth gradient. At 300 m depth, the temperatures increase by about 10 K.

These simulations can be used to estimate the presence, depth and availability of geothermal reservoirs, which require geothermal fluids at temperatures higher than 433 K to be exploited for high-temperature heat applications at the surface [45]. More accurate simulations of the temperature in the host rock above the magma chamber would require consideration of 3D chamber geometry, latent heat and groundwater/volatile circulation, unless the host rock has very little porosity and/or water content.

3.2.2. Velocity Results

Figure 6 shows the streamlines and velocities in the vertical v-direction within the chamber at two different times for the cases studied. Streamlines in Figure 6 are in the temperature zones above 1223 K (solidus temperature). In all three cases, natural convection produces movement of the magma. For the two times presented, the v velocities do not exceed 1.4 × 10⁻³ m/s or 31.5 km/y. Maximum velocities are observed near the lateral wall. In general, the lower local resultant velocities are in the range of 0.32–0.032 km/y. This means that in 100 years an element of the magma travels at least one complete revolution in the chamber. Case 1 is the one that retains the most heat and therefore has the most movement. Two big vortices are produced in the fluid zone in all cases. In cases 2 and 3, streamlines show the effects of the chamber shape. This is best seen at 5 ky when the fluid area is still large.

Figure 7a,b show the velocity profiles at x = 12,000 m. The absolute value of the velocities u and v are low at 9 × 10⁻⁴ m/s or less. Velocities within the chamber in each case are distributed in a different form. Cases 2 and 3 show positive and negative values of u and v. The velocities are displaced with time due to the solidification front for case 3. The velocity profile in case 3 has no contact with the lower part of the chamber; this produces velocities of zero near the bottom of the chamber. For this profile, case 1 displays the lowest velocities of the three studied cases. The time evolution of the average resulting velocity (V_res) in the chamber is shown in Figure 7c. The V_res increases rapidly during the first years of cooling. Only in case 1 is the V_res value (1.1 × 10⁻³ m/s) near the maximum calculated by Equation (11) (U₀ = 1.4 × 10⁻³ m/s). The maximum V_res are 5.5 × 10⁻⁴ and 1.6 × 10⁻⁴ m/s for cases 2 and 3, respectively. The V_res decreases with: case 1 > case 2 > case 3, which is consistent with the temperature decrease for the studied chambers.

The movement produced by natural convection could transport crystals to areas with different temperatures, where the addition of minerals other than those initially present could occur and thus generate the chemical differentiation that can be observed in plutonic rocks [17]. In the Humay unit of the PCB, crustal assimilation and magmatic recycling provide evidence for magma movement [44] which could result from the natural convection simulated in this paper.

3.2.3. Concentration Results

As illustrated in Figure 8, a simulation was run to determine the evolution of SiO₂ concentration when basaltic magma is intruded beneath a rhyolitic magma chamber. Both magma bodies initially had the same temperature. The basaltic magma initially had 47 wt% SiO₂, the rhyolitic magma 70 wt% and the host rock 50 wt%. Using these parameters and a typical SiO₂ diffusion coefficient at atmospheric pressure of D = 3 × 10⁻¹² m²/s [18], no variation was observed in the concentration of SiO₂ for the three cases (see Figure 8a for case 3). Additional tests were made with similar results using a self-diffusivity value of D = 6.7 × 10⁻¹¹ m²/s obtained at 1 GPa of pressure and 1873 K for SiO₂ in anhydrous basaltic liquid [46]. These diffusivity values are too low to observe the SiO₂ transport by molecular diffusion. Advection is not capable of mixing or mingling the two contrasting magma types. It is important to remember that for the magma chambers studied here the Rayleigh, Prandtl and Lewis numbers indicate that the diffusion of momentum is much greater than that of heat, which is much greater than the molecular diffusion of SiO₂. Therefore, it was expected that natural convection would not produce changes in the iso-concentration of SiO₂.

Our simulation results agree with experimental results that can be found elsewhere [47,48]. These authors concluded that if a new fluid is intruded into a chamber that originally had a fluid with different and homogenous concentration, the new fluid forms a separate layer at the base of the chamber. This is true even if vigorous convection in both fluids is occurring. At the same time there is a negligible compositional flux across the interface between the fluids.

A sensibility analysis was made by carrying out several simulations using artificial diffusion coefficients of 1 × 10⁻¹¹ to 1 × 10⁻⁶ m²/s. Results show that the minimum diffusivity necessary to observe diffusion of silica in the chamber is around 1 × 10⁻⁸ m²/s (Le = 16). The result obtained using a hypothetical value for the diffusivity coefficient of D = 1 × 10⁻⁸ m²/s is presented at 22 ky in Figure 8b. The advection–diffusion transport mechanisms change the concentration of SiO₂ while the magma is in the liquid phase. This occurs mainly in the first 10 ky of cooling. Figure 8b shows the variation in concentration level near the boundary between basalt and rhyolite. The SiO₂ content increases in the zone where initial basaltic magma intrudes. The rhyolite composition decreases from 0.70 to 0.60 and the basalt composition increases from 0.47 to 0.50 near the boundary interface between rhyolite and the basalt intrusion.

Molecular double diffusion or the existence of magnetic compounds that can be found in multicomponent fluids could increase convection and produce molecular stratification of minerals within the chamber [49,50]. These phenomena cannot be reproduced with our mathematical model.

4. Conclusions

Based on the cases proposed in this work for a magma chamber that formed a pluton in the Humay unit of the PCB, the simulations show that the main mechanism of heat transfer for magma cooling is conduction. The magma velocities produced by natural convection within the chamber do not transfer heat fast enough to observe an influence on the isotherms. That does not imply that natural convection does not participate in magmatic differentiation, since the movement of the fluid could transport crystals to areas with different temperatures in which other minerals can be added and thus form magma and rocks with different chemical composition from the original. Molecular diffusivity of SiO₂ under the conditions examined in this study produced little chemical differentiation of magma.

Shapes of the magma chamber and the interaction with the host rock determine the cooling rates of the magma chambers. Shapes that have a larger heat transfer area with the host rock will tend to cool faster. Additionally, wider shapes in the upper zone of the chamber will tend to cool faster than taller shapes due to the area exposed to lower temperatures of the host rock.

Based on the results reported in this work, it would seem adequate to affirm that a simpler mathematical model of heat conduction could be used to estimate the cooling times of a magma chamber with these specific conditions. However, a conduction model would not provide any information regarding the movement of the magma in the chamber, which can be observed in field and laboratory studies showing that there was convection. A coupled model of momentum, heat and mass such as the one proposed in this work is more suitable to study the cooling of chambers where magma flows and to explain the origin of what is observed in the field in plutons of the PCB.