Solver and run files for the paper "Numerical simulation of magma-rock interaction at Krafla volcano using OpenFOAM software and a simplified thermal model"

Tutorial for OpenFOAM solver multiMeltInterFoam Developed during Aurélie LOUIS-NAPOLEON’s PHD (alouisnapole[email protected]om) Additional input by Muriel Gerbault (muriel.gerb[email protected]) September 29, 2025 Content Introduction 3 1 Presentation of multiMeltInterFoamv2 4 1.1 Solved Equations .................................... 4 1.1.1 Partial melting ................................. 4 1.1.2 Spatial Heterogeneities ............................. 4 1.1.3 Effective Density and viscosity ........................ 5 2 Installation of OpenFOAM and the Solver 6 2.1 Installing OpenFOAM on a personal PC - Ubuntu ................. 6 2.2 Installing OpenFOAM on a personal PC - OpenSuse ................ 6 2.3 Installing OpenFoam on a lab server like IMFT ................... 7 2.4 Installing the solver multiMeltInterFoamv2 ...................... 7 3 Using multiMeltInterFoamv2 8 3.1 General remarks .................................... 8 3.2 Overview of files for the Test Case .......................... 9 4 Pre-Processing 10 4.1 Step 1 : Form the mesh geometry ........................... 10 4.2 Step 2 : Initiate the geotherm ............................. 10 4.3 Step 3 : Set the phases initial geometry (inclusions) ................. 12 4.3.1 Case of spherical inclusions .......................... 12 4.3.2 Case of layers of inclusions ........................... 13 4.4 Set the rheological properties and adimensionalize the problem .......... 17 4.5 Extra step : cooling and freezing of the inclusions .................. 19 4.5.1 Cooling ..................................... 19 4.5.2 Freezing the densities of the inclusions .................... 20 5 Processing 21 5.1 Step 4 : Run the calculation ............................. 21 6 Post-Processing 21 6.1 Step 5 : Visualise the results ............................. 21 6.1.1 Time ....................................... 21 6.1.2 Flow markers .................................. 22 6.1.3 Vertical average of quantities over the entire calculation domain ...... 23 6.1.4 Paraview macros ................................ 23 6.2 Step 6 : Data processing using Python codes .................... 23 6.2.1 Generate good quality videos ......................... 24 6.2.2 Viscosity and temperature profiles, velocity and melting front over time . 24 6.2.3 Exercise : apply these Python codes to the other cases tested in this tutorial 24 6.3 Summary of the main commands ........................... 24 7 To practice 24 7.1 Results .......................................... 25 7.2 Solution ......................................... 26 1 8 Parallel computation 26 8.0.1 The file system/decomposeParDict ...................... 26 8.0.2 Launch a parallel Run on your computer ................... 27 8.0.3 Launch a Run on a Distant server like Olympe on CALMIP ........ 27 9 3D models 29 9.1 Installing the solver on CALMIP ........................... 29 9.2 Run on CALMIP .................................... 29 9.3 Post-processing ..................................... 29 9.3.1 Paraview .................................... 30 9.3.2 Markers ..................................... 30 9.3.3 2D sections fields ............................... 30 10 Solver modifications, version V2 31 10.1 Correction of the viscosity in Pinkerton’s law .................... 31 10.2 Correction of ∆Cp in Wittington’s law ........................ 32 10.3 Modification of the latent heat ............................ 34 10.4 Taking into account the variation of ρin the heat equation ............. 35 10.5 Modifications in the Test Case input files ...................... 36 11 Solver modifications (v4): Darcy-Brinkman-Stokes formulation - in process 36 12 Solver modifications (v3): Krafla study - 2025 36 12.1 Equations and properties ................................ 37 12.2 Study summary and model setup ........................... 37 A Appendices 39 A.1 Setting up the folder 0................................. 39 A.2 Setting up the folder constant ............................ 41 A.3 Setting up the folder system ............................. 44 A.3.1 File controlDict ................................. 44 A.3.2 File fvSchemes ................................. 45 A.3.3 File fvSolution ................................. 45 A.3.4 Additional library for initial and boundary conditions: swak4Foam . . . . 46 A.4 Postprocessing ..................................... 48 A.4.1 Paraview .................................... 48 A.4.2 sampleDict ................................... 48 A.5 Reminder of basic Linux commands .......................... 49 A.5.1 Terminal commands .............................. 49 A.5.2 Aliases ...................................... 49 Figures list 50 2 Introduction This manual is linked to the PHD thesis (https: // www. theses. fr/ 2020INPT0124 ): Modeling gravity instabilities within the partially molten crust using a Volume-Of-Fluid method. The objective of this thesis was to study crustal partial melting in order to explain the formation of migmatite domes. Therefore, we used the OpenFOAM code and created a solver called multiMeltInterFoamv2. This solver combines the VOF method, which makes it possible to finely follow the interfaces of heterogeneities, with rheological laws that allow to model partial melting (presence of a melt rate M, and specific densities and viscosities laws). These heterogeneities correspond to more or less molten rocks which float in a partially molten environment (called BED). The studied system corresponds to figure 1. μNF M=MS L g H μNF M=MS ρ1 ρ3 ρ2 r TC=300°C TH+=1000°C Milieu non fondu Milieu partiellement fondu } MGE μ3 PF μ2 PF μ1 PF Inclusion lourde } } Inclusion légère ρ Figure 1: Model Setup, to simulate a partially melting crustal domain. The purpose of this manual is to explain how to run this solver. The first part presents the equations solved, the second part explains how to compile OpenFOAM and multiMeltInterFoamv2. The 3rd and 4th parts describe step by step how to start the calculation with a little practical exercise. The 5th (processing) and 6th parts explain how to visualize the calculation (post-processing with small exercises). Part 7 indicates how to carry the calculation in parallel both on your own machine and on the CALMIP Olympe server. Section 8 presents the pecularities of the solver for a 3-dimensional problem. Section 9 details the changes made to version V2 of the solver (2023) compared to V1 (2020)[Louis-Napoléon et al., 2020a, Louis-Napoléon et al., 2020b,Louis-Napoléon et al., 2022]. 3 1 Presentation of multiMeltInterFoamv2 1.1 Solved Equations multiMeltInterFoamv2 solves the Navier-Stokes equations : ∂Ci ∂t +U·∇Ci=−∇·(UrCr),(1) ∇·U= 0,(2) ρref ∂U ∂t +ρref U·∇ U=−∇P+ ˜ρg+∇·[˜µ(∇U+ (∇U)T)],(3) ρref Cp∂T ∂t +U·∇T=1 CA[∇·(ρref Cpκ∇T) + Hr],(4) with κ(T < 846K) = 567.3 T−0.062·10−6 κ(T > 846K) = (0.732 −1.35 ×10−4T)·10−6(5) Cp(T < 846K) = 902.7+0.387T−22.6×106T−2 Cp(T > 846K) = 1037.6+0.146T−216.7×10−6T−2(6) Only viscosity ˜µand density ˜ρwill be affected by partial melting and heterogeneities. The heat equation here takes into account internal heating, with Hr(en W.kg−1), the rate of radioactive heat production as well as variations in thermal diffusion (κ(T), Eqs. 5in m2.s−1) heat capacity (Cp(T), Eqs. 6, in m2.s−2.K−1) with temperature [Whittington et al., 2009]. allows the contribution of latent heat to be taken into account Lwhen there is partial melting (Mthe melting rate which itself depends on T, cf. next section) following the formalism of [Ulvrová et al., 2012]. Evaluation of CA(T)is expressed in terms of a Stefan number which is only non-zero in the interval Tsol ≤T≤Tliq, so that CA∼2in that interval and CA= 1 outside. 1.1.1 Partial melting In order to model partial melting, we make a first order approximation by introducing a single parameter, the melt fraction (at constant pressure) 0≤M≤1(at constant pressure). This parameter is defined to vary linearly with temperature: M=T−Tsol Tliq −Tsol ,(7) where, throughout the study, we fix Tsol = 600oCet Tliq = 1000oC. We also introduce a critical value of the melt fraction, noted MS, beyond which the viscosity decreases sharply. 1.1.2 Spatial Heterogeneities The solver can model several fluids with different physical properties (density, viscosity, etc.). In the studied cases, we consider 3 fluids, therefore three volume fractions 0≤Ci≤1(i= 1,2,3) : -Fluid 1 = equivalent bedrock domain (BED) (C1= 1,C2=C3= 0), 4 -Fluid 2 = light inclusions (C1= 0,C2= 1,C3= 0), -Fluid 3 = heavy inclusions (C1=C2= 0,C3= 1). This allows us to distinguish in the system a ’continuous’ phase called BED which may or may not be melted, and two ’discrete’ phases (in the sense that they are initially in the form of inclusions separated from each other). 1.1.3 Effective Density and viscosity Density ˜ρdepends on both temperature T(via the fluid state equation), the phase considered and the melt fraction. Note, however, that ˜ρis independent of the threshold MSunlike viscosity. Density is thus modeled by: ˜ρ=ρref ×[1 −α(T−Tref )],avec ρref = 3 X i=1 Ci×[ρS0 i(1 −M) + ρL0 iM],(8) where ρS0 iis the density of the unmolten rock in phase i,ρL0 ithat of the melting rock from this same phase and Tref a reference temperature for the fluid state equation (here, Tref = 300oC). For the viscosity ˜µ, two laws are introduced depending on the relative value of the melt fraction Mwith respect to the threshold MS. In one case, ˜µdepends on temperature and strain rate ˙ε, in the other case ˜µdepends on the melt fraction and the phase considered. Ultimately, the viscosity in the system is modeled by: ˜µ= µNF = 0.25 ×106(0.75A)−1 n(max( ˙εmin,˙ε)) 1 n−1·exp(Q nRT )si M < MS µPF = 3 X i=1 Ciµ0 i·exp " 2.5 + 1−M M0.48!(1 −M)#si M≥MS (9) where µNF (resp. µP F )refers to viscosity in unmelted rock (resp. partially melted), Ais the material constant in Pa−n.s−1,˙ε=p(tr(2)−tr2()/3) /2is the second invariant of the strain rate tensor = (∇U+(∇U)T)/2en s−1,Qis the activation energy in J.mol−1,nis the exponential factor Ris the universal gaz constant (R= 8.314 J.mol−1.K−1), and µ0 iis a characteristic viscosity for phase iwhen it is partially molten. Note that the equations 4and 9here are corrected with respect to those appearing in the thesis and in the publications [Louis-Napoléon et al., 2020a,Louis-Napoléon et al., 2022]. A 3rd publication (Louis-Napoléon et al., part II, in press at GJInt) includes these corrections. In order to ensure a physical numerical solution and limit the calculation time, we define a numerical Archimedes number ArN≤1, proportional to the real Archimedes number. The equations above are therefore rescaled according to a characteristic time Tand a characteristic distance H: t∗=t T,T=q−1=2Keff ∆ρgH n =2Keff ραg∆T H n , x∗ i=xi H, Ar =2ρq2−1 nH2 Keff .(10) This ratio of numerical vs. real Ar leads to a numerical vs. real viscosity of the system, cf. details in Section 4.4. 5 2 Installation of OpenFOAM and the Solver First, make sure that some packages are preinstalled. As of january 2024, please check https: //www.openfoam.com/documentation/system-requirements, for Ubuntu and Opensuse. Prosaically, specific libraries are recommended like : bear,bison,build-essential, flex,gitcore, libqt5,libncurses,libxt, libgmp,libglu1,openmpi,python,readline,zlib1g,... It may be useful to set the appropriate mpi environment in your home/.bashrc: export LD_PRELOAD ... export LIBRARY_PATH ... (cf . system - requirements webpage above) 2.1 Installing OpenFOAM on a personal PC - Ubuntu Follow the instructions in the following link to install OpenFOAM-v2206 on Ubuntu : https: //develop.openfoam.com/Development/openfoam/-/wikis/precompiled/debian. The "v2***" version branch should be preferred over branch 8,9, or 10. To test the installation under Ubuntu and learn the main steps, you can carry out the test case breaking of a dam, accessible via the online tutorial : https://www.openfoam.com/documentation/ tutorial-guide/4-multiphase-flow/4.1-breaking-of-a-dam. 2.2 Installing OpenFOAM on a personal PC - OpenSuse Since the OpenSuse Leap 15.2 distribution, Openfoam is provided precompiled in the specific repository science. For certain Linux distributions, follow the instructions via this link: https: //develop.openfoam.com/Development/openfoam/-/wikis/precompiled/. For example in October 2022, version OpenFoam2206 is installed in /usr/bin and /usr/ lib. We can use the Yast2 graphical tool and install openfoam2206 (openfoam2206)- common,-devel,-tools,-tutorials,-selector (the latter chooses the version to install itself), or use the shell command zypper: > zypper ar -f <URL > <alias > // Add the science repository > sudo zypper refresh // Refresh > sudo zypper install openfoam // Install package Note: You must be admin for the install. Then you can be user ("login"). Create the working directory OpenFOAM/login-version from home/login if it was not created by default (replace "version" with "2206" or the current version). We will name this path login-path in the following. Each time you use OpenFoam you must source the environment (that of the version you want to use if you have several), either automatically via the .bashrc or via an alias (in file /home/login/.alias). Another alias can allow you to place yourself in the working directory: > source /usr / lib/ openfoam / openfoam2206 / etc/ bashrc // source the environment > openfoam // the prompt ’>’ becomes ’$’ $ of06 // alias going to ’login -path ’ (" cd" command ) You can check the installation with the following commands, go to login-path and type: > foamSystemCheck // should display :.." Can continue to OpenFOAM installation ." > foamInstallationTest // displays paths and ...".. systems ok ." > foamTestTutorial -full incompressible / simpleFoam / pitzDaily // " Passed all 1 test" > cp -R /usr /lib / openfoam / openfoam2206 / tutorials / multiphase / interFoam / laminar / damBreak/damBreak . > cd damBreak > cp -R 0. orig 0 6 > ./ Allrun > paraFoam // should display a rainbow shape pierced with a rectangle . In order to visualize the results (post-processing), OpenFoam has a specific command "paraFoam", who basically sets Paraview in the OpenFoam environment. Previously, you had to install a set of libraries and packages grouped under the name "ThirdParty" yourself. If this directory still exists, it is probably empty, and it is therefore sufficient that paraview be pre-installed (via Yast2 or zypper). If a problem of the following type appears: ParaView_DIR =/ home / login / OpenFOAM / ThirdParty -4. x/ platforms / linux64Gcc / ParaView -5.0.1 sed : can t read / home / login / OpenFOAM / ThirdParty -4. x/ ParaView -5.0.1/ VTK / Utilities / kwsys / kwsysDateStamp . cmake : No such file or directory Do : sed -ibck -e ’s /\[3456\]/[34567]/ ’ ParaView -5.0.1/ VTK / CMake /{ GenerateExportHeader , vtkCompilerExtras }. cmake There may also be version conflicts between Paraview, OpenGL and common Linux kernel tools (for example Paraview v5.10 is available under Leap 15.3, but it requires OpenGL v3.2 which would not be compatible with the Linux kernel of Leap 15.3. Tip: Downgrade the version of Paraview (and at worstOpenFoam). 2.3 Installing OpenFoam on a lab server like IMFT OpenFOAM is installed at IMFT in different versions. To use v****, set : export LMOD_DEFAULT_COLLECTION = All source / PRODCOM / bin / config . sh module load Collection / All module load Collection /. Intel -2019.0 module load Collection /. Intel -2019.0 - MPICH module load PETSc /3.9.4 - intel -2019.0 - mpich -3.2.1 module load MUMPS /5.1.2 - intel -2019.0 - mpich -3.2.1 - parmetis echo $EBROOTPETSC / PRODCOM / Ubuntu18 .04/ Intel -2019.0 - mpich -3.2.1/ PETSc /3.9.4 - intel -2019.0 - mpich -3.2.1 module load openfoam / openfoam % alias of41=" source ␣/ PRODCOM / bin /u18 / openfoam -4.1. sh " Then source the right one ../etc/bashrc in a terminal. To check that it works, you can type : simpleFoam -help. 2.4 Installing the solver multiMeltInterFoamv2 multiMeltInterFoamv2s a solver which is not in the basic version of Open-FOAM, so you have to compile it yourself on your machine. The solver and the associated test case can be downloaded here : https://gitlab.com/AurelieLN/multiMeltInterFoam.git (vsee .zip). It is essential to compile the solver at the right place, and there are 2 options: — First, make sure you have created the folder ’login-path’ somewhere (cf. 2.2, and A.5.1 for a Linux command guide), if it has not already been created : mkdir OpenFOAM && cd OpenFOAM mkdir $WM_PROJECT_USER // if this doesn ’t work , the directory must be named login -version , example : alouis - v2206 cd $WM_PROJECT_USER //(or cd login - path) 7 — Then install the solver: mkdir -p applications / solvers / multiphase cd applications / solvers / multiphase // set multiMeltInterFoam directory here cd multiMeltInterFoamv2 ./ Allwmake — Alternatively, put the folder multiMeltInterFoam in the local OpenFoam install directory and recompile it there. If during this compilation errors appear : 1) phase / phase .C :174:31: error : call of overloaded autoPtr ( NULL ) is ambiguous →Repaired by replacing in phase.C:’autoPtr’ by ’null’. cf. section 8.0.3. 2) / usr /lib /openfoam / openfoam2206 / src/ finiteVolume / lnInclude / MULESTemplates .C :392:16: error : no matching function for call to min (Foam ::tmp < Foam :: Field < double> >, const int &) →Repaired by replacing in multiMeltMixture.C lignes ∼630 (loop for MULES): "1,0" by "oneField(), zeroField()",then "phasei++" by "++phasei". 3 Using multiMeltInterFoamv2 3.1 General remarks OpenFOAM is used in three main steps: •pre-processing (mesh and geometry, launched by commands blockMesh, setFields,..) •the Run (launched with a solver’s name which ends with foam, here multiMeltInterFoamv2) •post-processing (comparison and analysis of the results, here paraFoam). OpenFOAM offers multiple tutorials and test cases on various types of fluid flow problems. Each test case consists of three directories: -a directory 0which contains the initial data of the physical quantities used (temperature, pressure, speed, etc.), -a directory constant which corresponds to the directory of calculation parameters (viscosity, density, etc.) and the mesh (geometry and edge conditions), -a directory systemwhich contains information such as the numerical schemes and solver used, time steps... In order to execute the commands necessary for the calculation to proceed, you must source with an alias command like of41 or run openfoam with openfoam2206 (which does the same), cf. section 2.3,2.2). All commands must be executed IN the case folder. Each OpenFOAM file is written in C++ and has a header which specifies that this file is part of OpenFOAM and which version, in which folder it is located and its name. 1/* -------------------------------*- C++ -*---------------------------------*\ 2| ========= | | 8 42 43 write (m, ∗)’boxToCell ’ 44 write (m, ∗)’{’ 45 write (m,FMT=’(A9 ,F0 .3,A9 ,␣ F0 .3, A4 )’ )’␣box␣(-1␣’ , y , ’␣ -1␣)␣(1␣’ ,y+r , ’␣1);’ ! on met au bon format l e s nombres pour e v i t e r l e s decimales non d e s i r e e s 46 write (m, ∗)’fieldValues␣(␣␣volScalarFieldValue ␣alpha.shamp␣0’ 47 write (m, ∗)’␣ volScalarFieldValue ␣alpha.water␣0’ 48 write (m, ∗)’␣ volScalarFieldValue ␣alpha.oil␣1’ 49 write (m, ∗)’);’ 50 write (m, ∗)’}’ 51 write (m, ∗)’boxToCell ’ 52 write (m, ∗)’{’ 53 write (m,FMT=’(A9 ,F0 .3,A9 ,␣ F0 .3, A4 )’ )’␣box␣(-1␣’ , y+2∗r , ’␣ -1␣)␣(1␣’ ,y+3∗r , ’␣1); ’ 54 write (m, ∗)’fieldValues␣(␣␣volScalarFieldValue ␣alpha.shamp␣1’ 55 write (m, ∗)’␣ volScalarFieldValue ␣alpha.water␣0’ 56 write (m, ∗)’␣ volScalarFieldValue ␣alpha.oil␣0’ 57 write (m, ∗)’);’ 58 write (m, ∗)’}’ 59 60 y = y + 4∗r 61 compt=compt+1 62 end do 63 64 write (m, ∗)’);’ 65 write (m, ∗)’//␣ ********************************************************************//’ 66 END program creafunkySETfIELDS Interfaces with white noise - Here we first define layers with setFields, then we insert white noise on in between each layers. Since we want white noise on the mille-feuilles layering created above, oa use here a fortran program to generate white noise on all interfaces (create_funky_bruitblanc_millefeuille). First we set white noise on the light inclusions interfaces surrounded by BED (upper and lower interfaces), then we suppress the white noise of the BED. We do the same for the heavy phase. 1program creafunkySETfIELDS 2i m p l i c i t none 3i n t e g e r : : m, i , hauteur , k , zi , toto , kprime , compt 4r e a l : : Lx , Lz , Lcouche , s ur fa ce , r , l a r g e u r t r i o , nombretriodsLx , Ly , x , y , deltay , b u l l e l a r g e ur , z 5open (1 , f i l e="funkySetFieldsDict") 6m = 1 7 8write (m, ∗)’/*--------------------*-␣C++␣ -*---------------------*\’ 9write␣(m,*)␣’ | ========= | | ’ 10 write␣(m,*)␣’ | \\ / F i e l d | OpenFOAM: The Open Source CFD Toolbox |’ 11 write␣(m,*)’ | \\ / O p erati on | Version : 4.0 |’ 12 write␣(m,*)’ | \\ / A nd | Web: www.OpenFOAM. org |’ 13 write␣(m,*)’ | \\/ M an ip ul at io n | | ’ 14 write␣(m,*)’\∗−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−∗/’ 15 write␣(m,*)’FoamFile’ 16 write␣(m,*)’{’ 17 ␣␣write␣(m,*) ’ v er sion 2 . 0 ; ’ 15 18 ␣␣␣write␣(m ,*)’ format a s c i i ; ’ 19 ␣␣write␣(m,*) ’ class dictionary ;’ 20 ␣␣write␣(m,*) ’ location "system" ;’ 21 ␣␣␣write␣(m ,*)’ object "funkySetFieldsDict";’ 22 write␣(m,*)’}’ 23 write␣(m,*)’//**************************//’ 24 write (m, ∗)’defaultFieldValues’ 25 write (m, ∗)’(’ 26 write (m, ∗)’␣␣volScalarFieldValue ␣alpha.shamp␣0’ 27 write (m, ∗)’␣ volScalarFieldValue ␣alpha.water␣1’ 28 write (m, ∗)’␣ volScalarFieldValue ␣alpha.oil␣0’ 29 write (m, ∗)’);’ 30 write (m, ∗)’expressions’ 31 write (m, ∗)’(’ 32 33 r = 0.006 34 Lz = 0. 5 35 Lx = 0 .5 36 Ly = 0.35 37 deltay= 0.001 != deltay /2 38 y = r 39 compt= 0 40 do while ( compt < 7) 41 write (m, ∗)’␣ noise ’ 42 write (m, ∗)’␣{␣’ 43 write (m, ∗)’field␣alpha.oil␣;’ 44 write (m, ∗)’expression ␣"␣ rand (0) "; ’ 45 write (m,FMT=’(A22 , F0 .3 ,A16 ,␣ F0 .3, A4)’ )’condition␣␣"( pos().y<’ ,y , ’)␣&&␣(pos() .y>=’ ,y−deltay ,’)";␣’ 46 write (m, ∗)’keepPatches␣1;’ 47 write (m, ∗)’}’ 48 49 write (m, ∗)’␣␣noise␣’ 50 write (m, ∗)’{␣’ 51 write (m, ∗)’field␣alpha.water␣;’ 52 write (m, ∗)’expression ␣"1-␣ rand (0) ";’ 53 write (m,FMT=’(A22 , F0 .3 ,A16 ,␣ F0 .3, A4)’ )’condition␣␣"( pos().y<’ ,y , ’)␣&&␣(pos() .y>=’ ,y−deltay ,’)";␣’ 54 write (m, ∗)’keepPatches␣1;’ 55 write (m, ∗)’}’ 56 57 write (m, ∗)’␣ noise ’ 58 write (m, ∗)’␣{␣’ 59 write (m, ∗)’field␣alpha.oil␣;’ 60 write (m, ∗)’expression ␣"␣ rand (0) "; ’ 61 write (m,FMT=’(A22 , F0 .3 ,A16 ,␣ F0 .3, A4)’ )’condition␣␣"( pos().y<’ , y+r , ’)␣&&␣(pos ().y>=’ ,y+r−deltay ,’)";␣’ 62 write (m, ∗)’keepPatches␣1;’ 63 write (m, ∗)’}’ 64 65 write (m, ∗)’␣␣noise␣’ 66 write (m, ∗)’{␣’ 67 write (m, ∗)’field␣alpha.water␣;’ 68 write (m, ∗)’expression ␣"1-␣ rand (0) ";’ 69 write (m,FMT=’(A22 , F0 .3 ,A16 ,␣ F0 .3, A4)’ )’condition␣␣"( pos().y<’ , y+r , ’)␣&&␣(pos ().y>=’ ,y+r−deltay ,’)";␣’ 70 write (m, ∗)’keepPatches␣1;’ 71 write (m, ∗)’}’ 72 73 write (m, ∗)’␣ noise ’ 16 74 write (m, ∗)’␣{␣’ 75 write (m, ∗)’field␣alpha.shamp␣;’ 76 write (m, ∗)’expression ␣"␣ rand (0) "; ’ 77 write (m,FMT=’(A22 , F0 .3 ,A16 ,␣ F0 .3, A4)’ )’condition␣␣"( pos().y<’ , y+2∗r , ’)␣&&␣( pos ().y>= ’ , y+2∗r−deltay ,’)";␣’ 78 write (m, ∗)’keepPatches␣1;’ 79 write (m, ∗)’}’ 80 81 write (m, ∗)’␣␣noise␣’ 82 write (m, ∗)’{␣’ 83 write (m, ∗)’field␣alpha.water␣;’ 84 write (m, ∗)’expression ␣"1-␣ rand (0) ";’ 85 write (m,FMT=’(A22 , F0 .3 ,A16 ,␣ F0 .3, A4)’ )’condition␣␣"( pos().y<’ ,y+2∗r , ’)␣&&␣( pos ().y>= ’ , y+2∗r−deltay ,’)";␣’ 86 write (m, ∗)’keepPatches␣1;’ 87 write (m, ∗)’}’ 88 89 write (m, ∗)’␣ noise ’ 90 write (m, ∗)’␣{␣’ 91 write (m, ∗)’field␣alpha.shamp␣;’ 92 write (m, ∗)’expression ␣"␣ rand (0) "; ’ 93 write (m,FMT=’(A22 , F0 .3 ,A16 ,␣ F0 .3, A4)’ )’condition␣␣"( pos().y<’ , y+3∗r , ’)␣&&␣( pos ().y>= ’ , y+3∗r−deltay ,’)";␣’ 94 write (m, ∗)’keepPatches␣1;’ 95 write (m, ∗)’}’ 96 97 write (m, ∗)’␣␣noise␣’ 98 write (m, ∗)’{␣’ 99 write (m, ∗)’field␣alpha.water␣;’ 100 write (m, ∗)’expression ␣"1-␣ rand (0) ";’ 101 write (m,FMT=’(A22 , F0 .3 ,A16 ,␣ F0 .3, A4)’ )’condition␣␣"( pos().y<’ , y+3∗r , ’)␣&&␣( pos ().y>= ’ , y+3∗r−deltay ,’)";␣’ 102 write (m, ∗)’keepPatches␣1;’ 103 write (m, ∗)’}’ 104 105 y = y + 4∗r 106 compt=compt+1 107 108 end do 109 110 write (m, ∗)’);’ 111 write (m, ∗)’//*************************************************************//’ 112 END program creafunkySETfIELDS 4.4 Set the rheological properties and adimensionalize the problem We have two viscosity laws and an interpolated density law, cf. 1.1. The fluid parameters associated with these laws are stored in constant>transportProperties. Some parameters are not the same numerically (N) and in reality (R). The format of the physical properties is the same for the 3 fluids (water, oil and shamp), and only water is presented below. At the end of the file are provided the threshold value MS, activation or not of the thermal capacity Cp and thermal diffusion kappa, depending or not on the temperature. There is also the minimum strain rate ˙εmin, imposed to calculate an initial viscosity µNF , and Stefan number provided for the latent heat (Lt/(Cp.∆T)∼1.0), cf. Section 1.1. 17 1phases 2( 3water 4{ transportModel Newtonian ; // peu importe 5nu nu [ 0 2 −1 0 0 0 0 ] 6 . 4 ; //useless 6rhoL rhoL [ 1 −3 0 0 0 0 0 ] 2700; 7rhoS rhoS [ 1 −3 0 0 0 0 0 ] 2800; 8Tref Tref [ 0 0 0 1 0 0 0 ] 573; 9beta beta [ 0 0 0 −1 0 0 0 ] 3 . e−5 ; 10 nu0 nu0 [ 0 2 −1 0 0 0 0 ] 1 . 8 5 e −3;// nu0N=mu0R /(2700*10^15) 11 Arrho0 Arrho0 [ 0 0 0 0 0 0 0 ] 1. e18 ; // ArN= ArR *(10^5) ^n*1. e10 12 nr nr [ 0 0 0 0 0 0 0 ] 2 . 3 ; 13 Er Er [ 1 2 −2 0 −1 0 0 ] 154000.0; 14 Tsol Tsol [ 0 0 0 1 0 0 0 ] 873.0 ; 15 Tliq Tliq [ 0 0 0 1 0 0 0 ] 1273.0 ; 16 A0 A0 [ 1 −1−3 0 0 0 0 ] 1.75 e4 ; // A0N = HrN = 10^10* HrR 17 Cp Cp [ 0 2 −2−1 0 0 0 ] 10 0 0 . 0 ; // used if Cp constant 18 kappa kappa [ 0 2 −1 0 0 0 0 ] 1 . e −6; //used if kappa constant 19 } 20 // others fluides .... 21 valueM 1.0 ; // Useless ! 22 Mseuil 0.3 ; // Mseuil=1 when no melt(law UF), Mseuil=0 when melt(law PF) 23 KTok 1.0 ; // 1.0 if kappa(T), 0 if kappa constant 24 CpTok 1.0 ; // 1.0 if Cp(T), 0 if Cp constant 25 epsmin 1.e−6; // epsminN = 10^10 epsminR , minimal value of strain -rate 26 St 1 . 0 ; // Stefan number for latent heat: Lt/(Cp.DT)~4e5 /(1e3.H) The meaning of these parameters is : Parameters Code PHD numerical value Field Value Unit Rocks density PF rhoL ρL0 i2700 kg.m−3 Rocks density NF rhoS ρS0 i2800 Reference temperature Tref Tref 573 K Thermal expansion beta β3×10−5K−1 Viscosity prefactor µPF /ρ nu0 ν0 i1.85 ×10−31.85 ×1012 m2.s−1 Parameters Arrho0 A1018 3.2×10−4Pa−n.s−1 for viscosity nr n2.3 µNF Er Q154 kJ.mol−1 Solidus temperature Tsol Tsol 873 K Liquidus temperature Tliq Tliq 1273 K Radioactive heat production A0 Hr 1.75 ×1041.75 ×10−6W.m3 Thermal capacity Cp Cp 1000 m2.s−2.K−1 Thermal diffusion kappa κ10−6m2.s−1 Melt Threshold MMseuil MS0.3 Min. Strain-rate epsmin ˙εmin 10−610−16 s−1 Stefan number St St 1 Table 2: Equivalence between physical parameters in the code and in reality. Be aware that OpenFoam defines its numerical convergence parameters according to the FIRST phase listed in file transportproperties. It is thoughtful to define this 1st phase as the leading phase, which controls the model dynamics, hence the following a-dimensionalisation. 18 Note that the code uses the kinematic viscosity ν=µ/ρ (instead of µthe dynamic viscosity). The choice of the models dimensions made for the thesis was to consider that : HN HR= 10−5, and that the numerical Archimedes number ArN=Ar0 rho = 1. From there, we can calculate the different numerical quantities (N) different from those of reality (R) : •AN=Arrho0 = AR×(105)n×1010, •A0N=HrN= 1010 ×HrR, •νN 0=µR 0 ρ×1015 with ρ= 2700 kg.m−3. As another example (Malissa Bakouche, 2023), we consider a physical domain HN HR= 1 of size 50 ×50 cm and a real viscosity µR 0= 1010 Pa.s (Pinkerton law with K=µ0and n= 1). So we first take the definition of the physical characteristic time : t= 1/q = ( 2K ∆ρgH )n=1,hence with Boussinesq: =2K αρg∆T H ,(11) to obtain t∼16 yrs with α= 3.10−5and ∆T= 100,ρ= 2500. Also, q= 1.875 ·10−9(1/s). But in case of no thermal gradient ∆T, hence no Boussinesq approximation, one must use Has the characteristic distance in between ∆ρheterogeneities inside the model domain. We then express the Real Archimedes number, with K∼µ0: ArR=2ρq2−1/nH2 K,hence with Boussinesq : =αρ2g∆T.H3 µ2 0 (12) Since here all the numerical and real values are identical except the viscosities, we have : ArR ArN=(µN 0)2 (µR 0)2⇒µN 0=µR 0·rArR ArN.(13) If we want ArN≤1to ensure physical numerical solutions, we obtain µN 0=µR 0·(AR)1/2=ρ.Hpαg∆T.H, hence µN 0≥ρ.Hpαg∆T.H (14) By transforming back into kinematic viscosity, νN 0=µR 0/ρ ≥H√αg∆T H. The numerical application for this case provides νN 0≥0.05 (for the leading BED phase). 4.5 Extra step : cooling and freezing of the inclusions 4.5.1 Cooling To modify the basal temperature over time, we use a Fortran script which writes a timetemperature list: 1PROGRAM dimT 2i m p l i c i t none 3i n t e g e r : : i , n 4r e a l : : a , T, ti , Trf , Tchaud , tf , b , duree 5 6Open (10 , f i l e =’Temperatures_parentheses.txt’) 7Trf = 873.15 19 8Tchaud = 1000.0+273.15 9t i = 14000 10 duree = 5 11 duree = duree /0.000317 12 t f = t i + duree 13 a = ( Trf−Tchaud) /( tf −t i ) 14 b = Tchaud −a∗t i 15 Write (10 ,∗)’type␣uniformFixedValue ␣;’ 16 Write (10 ,∗)’uniformValue␣table’ 17 Write (10 ,∗)’(’ 18 do while ( t i < t f ) 19 T = a∗t i +b 20 Write (10 ,∗)’(␣’ , t i , ’␣’ ,T , ’)␣’ 21 t i = t i + 1.0 22 end do 23 Write (10 ,∗)’);’ 24 c l o s e (10) 25 END PROGRAM dimT which writes a file Temperature_parenthese.txt 1type uniformFixedValue ; 2uniformValue table 3( 4( 14000.0000 1273.14990 ) 5( 14001.0000 1273.12463 ) 6( 14002.0000 1273.09924 ) 7. . . 8( 29772.0000 873.172058 ) 9) ; We copy it into 0/T and check the bottom boundary condition : 1lowerWall 2{ type uniformFixedValue ; 3uniformValue table 4(... // Temperatures 5) ; 6}; 4.5.2 Freezing the densities of the inclusions To freeze the inclusions densities, we make them independent of the melt; in constant/transportProperties,we change the densities by averaging the previous densities (M= 0.5) : 1o i l 2{ transportModel Newtonian ; 3rhoL rhoL [ 1 −3 0 0 0 0 0 ] 2650; // avant 2500 4rhoS rhoS [ 1 −3 0 0 0 0 0 ] 2650; //avant 2800 5} 6shamp 7{ transportModel Newtonian ; 8rhoL rhoL [ 1 −3 0 0 0 0 0 ] 2850; // avant 2900 9rhoS rhoS [ 1 −3 0 0 0 0 0 ] 2850; // avant 2800 10 } 20 5 Processing 5.1 Step 4 : Run the calculation Type, in the terminal shell: multiMeltInterFoamv2 The file system>controlDict allows you to manage recordings of calculation time steps, restarting, etc. The options of controlDict are detailed in appendix A.3.1. For parallel computation see section 8. 6 Post-Processing 6.1 Step 5 : Visualise the results Once the calculation has been carried out, files appear as the calculated time step folders : -rhoB, for the density which varied with temperature(Eq. 8), -nu, for the viscosity ν, -epsilon, for the strain-rate ˙ε, -M, for the volumeric melt fraction M. These are not initially present in the directory 0. To avoid problems with visualization, I suggest copying these new quantities that appeared in the first time step folder 0. So that these quantities are closer to their value at time 0, we can launch a calculation at time 0.001. If you want to add or delete output files for each time step (in order, for example, to reduce their weight), you must modify the solver multiMeltinterFoamv2/createFields.H. For each quantity to be written/deleted, you must activate/deactivate, then recompile the solver : 1IOobject : :AUTO_WRITE / NO_WRITE The following paragraphs propose useful indications to produce figures: 6.1.1 Time The calculation time is in seconds. For the PHD cases, conversion from numerical time to real time is tR= 1010tNwith tin seconds. Shifting to million years (1 Ma ∼3.151013s) leads to tR(Ma)=0.000317 ×tN(s). Filters>AnnotateTimeFilter , then adjust Format and Scale according to indications given in Fig. 5. One might have to change the syntax Format depending on Paraview version, to t=time:0.1f Ma. In another example (that of Malissa with HR=HN= 0.5m), only viscosity is associated with the control of the numerical vs. real Archimedes number, and with the relations eqs. 13, we express tR tN=qN qR=KR KN=µR 0 µN 0 =1 5×10−9so that tR= 2 ×108tN. To convert to millions years we divide by 3.15 ×1013, so that tR(Ma)=6.345 ×10−6×tN(s). 21 Figure 5: Adapting the numerical time to the real time (here tN= 14900 s⇒tR= 4.7Ma) 6.1.2 Flow markers Markers in the flow illustrate the duration of a convective cycle. Proceed as follows : 1.Start at timestep 0, 2.Sources>Point Source > 0.246 0.09 0 and Point Size at 0, 3.Filters>Temporal>ParticleTracer Input: CasTest.OpenFOAM and Seed Source: PointSource, then select Static Seeds and Static mesh, unselect Compute Vorticity and set Force Reinjection NSteps to 1 (cf. Fig. left), 4.to the right of RenderView1, click on the double vertical rectangle then select SpreadSheetView>ParticleTracer1, 5.Select the cell like shown on the figure below, 6.Plot Selection Over Time, displays new curves, 7.Select the desired curve in the menu list Properties, then store by selecting File > Save data 22 Details for each step: 1.The square area activated at x=0.249, y=0.09, z=0, is displayed on the 2D plot and we can display its properties or simply a color (’Solid Color’). We can modify the units beforehand by using Calculator. 2.In the menu Pipeline, 3 lines of ’ParticleTracer1’ appear, it is of no importance and just select the first line, cf. area highlighted in blue in the Figure. 3.The ’Excel’ lines that appear correspond to distinct time steps. 4.It is important to select the first line, which corresponds to the initial state 5.Select/eliminate the desired curves in Properties, as well as line colors, axis labels, etc., as in a usual PlotOverLine. These outputs remain difficult to manipulate: the curves are displayed according to the time step, and not time. A Python code inclusions.py offers greater freedom of display, cf. section 6.2. This code also allows to follow the trajectory of the marker in space. 6.1.3 Vertical average of quantities over the entire calculation domain Paraview doesn’t seem to be able to give the vertical average of a quantity. To do this, I used a Fortran program average.f95 which calculates the vertical average directly from the output files. Additionally, it calculates the domain average of vertical velocity over time. To use it, you must modify the numbers of cells according to axes xand y. You can also modify the time steps over which to perform the averages. 6.1.4 Paraview macros Once the visualization is carried out on paraview is appropriate , it can be solved in a file .pvsm with File>save state> fichierVisu.pvsm. To re-use it, click on File>load state> fichierVisu.pvsm, choose the tab Choose file and select CasTest.OpenFOAM. 6.2 Step 6 : Data processing using Python codes If Paraview is useful to get a first idea of the evolution of many fields and to visualize the system in 3D, results rendering are often not great. Small python codes (folder Codes_Python) allow to process data quickly and easily and to make figures in .pdf. Therefore copy the scripts .py in the run directory (CasTest here). To apply these python files to your specific case, we recommend you modify only : -the paths of the case files, the names of the cases, -the times measured, the size of the scales, the location of the legends, -the size of the domain (Lx, Ly), the number of cells along the x and y axes (nx, ny), a +1 to add to the size of an array (fluids 2 and 3 after reshape, 22→23, l.137 of Films.py). 23 6.2.1 Generate good quality videos Use Films.py (copied in CasTest), which generates .jpg figures, in a directory CasTest/images that will have been created beforehand. The film can then be creatd by going into images and by typing (here we create two movies for two cases I1M1H2 and I1M1H3 - remove the loop if there is only one case). One might have to add +1 to the arrays fluid2_filt, l.219: mkdir images && python3 Films . py && cd images for (( j=2; j <= 3; j++ )) do ffmpeg -f image2 -i hqI1M1H "$j" -%04d. jpg -r 24 - vcodec mpeg4 -b 15000 k - filter:v " setpts =4.0* PTS " I1M1H "$j". mp4 done or, to build a gif (here for a single case): ffmpeg -f image2 -i hqI1M1HX .%04d. jpg hqI1M1HX . gif 6.2.2 Viscosity and temperature profiles, velocity and melting front over time Use the script Traces_champs_frontfusion_vitessmoy.py to generate figres like in publication Louis-Napoleon et al., 2022. 6.2.3 Exercise : apply these Python codes to the other cases tested in this tutorial . 6.3 Summary of the main commands // go into the case Folder and type: of41 // or other alias to source OpenFOAM ( source etc / bashrc ) // Pre - processing blockMesh // Create the geometry and the mesh setFields // Create the inclusions funkySetFields OR champT .f => 0/ T // Create the geotherm // Processing multiMeltInterFoamv2 // Run the solver // Post - Processing paraFoam // Display the results on paraview 7 To practice Copy the folder Cas Test and modify it to build a new case with the following characteristics : -Mesh : 247 ×175 ×1(x×y×z), -Inclusions radius: 600 m and they touch each other, -µ0 1= 5 ×1015 Pa.s, -µ0 2= 1015 Pa.s, -µ0 3= 25 ×1015 Pa.s. After 5000 seconds, things start to happen. Run until 150000 s, which may take about 2h15. 24 16 } 17 i n t e r p o l a t e false ; 18 } 19 ) ; 20 fields 21 (U T alpha . o i l alpha . water alpha . shamp nucreate ) ; If this is done before the calculation, post-processing is automatic. Be careful, however, it may not work well if the calculation is parallelized and you will therefore have to run the command: solveurName -postProcess. Again, pay attention to the OpenFOAM version used ! The post-processing files will then be stored in a folder /postProcessing which will be used with the Python code Codes_Python/Coupe_Champ_3D.py to produce 1, 3, or 6 .pdf figures. Note: if you wish X axis sections instead of Z axis, i) replace "zCut" by "xCut" in system/singleGraph, and ii) adapt the axes ny,nx to ny,nz in the .py code. 10 Solver modifications, version V2 Several changes have been made to the solver for paper "part II" (Louis-Napoléon et al., 2023, in press) compared to the previous ones. They are listed below. In order to compile this version 2 or any other modified version, you must copy the new directory multiMeltInterFoamv2 (from GitLab) in the OpenFoam build folder (currently OpenFoam-v2106). Then go to this directory, and type ./Allwmake. 10.1 Correction of the viscosity in Pinkerton’s law In /multiMeltInterFoamv2/multiMeltMixturev2/multiMeltMixturev2.C 1Foam : : tmp<Foam : : vo l Sc al ar F ie ld > 2Foam:: multiMeltMixture :: nu() const // viscosity 3{ 4dimensionedScalar R("R" , dimensionSet(1, 2, −2, −1, −1, 0 , 0) , 8 .31 4) ; 5dimensionedScalar T0("T0" , dimensionSet (0 , 0 , 0 , 1 , 0 , 0 , 0) , 27 3.0) ; 6dimensionedScalar Adim("Adim" , dimensionSet(1, −1, −3, 0 , 0 , 0 , 0) ,1.0) ; 7// dimensionedScalar visc ("visc", dimensionSet (0, 2, -1, 0, 0, 0, 0) ,370.); 8// dimensionedScalar ViscMin (" ViscMin ", dimensionSet (0 , 2, -1, 0, 0) ,0.0037); 9PtrDictionary<phase >:: c onst_iter at or i t e r = phases_ . begin ( ) ; 10 11 const v o l S c a l a r F i e l d& Mrho = U_. mesh ( ) . lookupObject<v ol Sc al ar Fi e ld >(" Mrho" ) ; 12 const v o l S c a l a r F i e l d& M = U_. mesh ( ) . lookupObject<v ol Sc al a rF i e ld >("M" ) ; 13 const v o l S c a l a r F i e l d& s tr ai nR at e = U_. mesh ( ) . lookupObject< vo lS ca l ar Fi el d >(" strainRate" ) ; 14 tmp<v ol Sc al a rF ie ld > tnu = i t e r ( ) ∗( i t e r () . nu0 ( ) ∗exp ((2.5+ (pow ( ( 1 . 0 −Mrho) / (max(Mrho , 0 . 1 ) ) , 0 . 48 ) ) ) ∗( 1 . 0 −Mrho) ) ∗M + ( i t e r () ∗0.25 ∗v i s c ∗pow(( i t e r ( ) . Arrho0 ( ) ∗0. 75 ) , −1.0 / i t e r () . nr () ) ∗exp ( i t e r () . Er () / ( i t e r () . nr () ∗R∗(max(T_, T0) ) ) ) ∗pow( strainRate , ( ( 1 . 0 / i t e r ( ) . nr ( ) ) −1 . 0 ) ) ) ∗(1 −M) ) ; 15 16 v o l S c a l a r F i e l d &nu = tnu . r e f ( ) ; 17 18 for (++ i t e r ; i t e r != phases_ . end ( ) ; ++i t e r ) 19 { 31 20 nu += i t e r ( ) ∗( i t e r () . nu0 ( ) ∗exp ((2.5+ (pow ( ( 1 . 0 −Mrho) / ( max(Mrho , 0 . 1 ) ) , 0 . 48 ) ) ) ∗( 1 . 0 −Mrho) ) ∗M + max( ViscMin , i t e r ( ) ∗0.25 ∗v i s c ∗pow(( i t e r ( ) . Arrho0 ( ) ∗0. 7 5) , −1.0 / i t e r () . nr ( ) ) ∗exp ( i t e r ( ) . Er ( ) / ( i t e r ( ) . nr ( ) ∗R∗ (max(T_, T0) ) ) ) ∗pow( s tra in Rat e , ( 1 . 0 / i t e r ( ) . nr ( ) ) − 1 . 0 ) ) ∗(1 −M) ) ; 21 } 22 return tnu ; 23 } 24 [...] 25 tmp<s c a l a r F i e ld > tnu = i t e r () ∗( i t e r () . boundaryField ( ) [ patchi ] ∗i t e r () . nu0 () . value ( ) ∗exp ( ( c oe f+ (pow ( ( 1 . 0 −Mrho . boundaryField ( ) [ pat ch i ] ) / (max(Mrho . boundaryField ( ) [ p at ch i ] , 0 . 1 ) ) , 0 . 48 ) ) ) ∗ ( 1 . 0 −Mrho . boundaryField () [ patchi ] ) ) ∗M. boundaryField () [ patchi ]+ max( ViscMin , i t e r () . boundaryField ( ) [ patchi ] ∗0.25 ∗v is c 2 ∗pow (( i t e r ( ) . Arrho0 () . value ( ) ∗0 .7 5) , −1.0 / i t e r () . nr ( ) . value () ) ∗ exp ( i t e r () . Er () . value () / ( i t e r () . nr ( ) . value ( ) ∗R2 ∗(max(T_. boundaryField () [ patchi ] , T02) ) ) ) ∗pow( strainRate . boundaryField ( ) [ pat ch i ] , ( 1 . 0 / i t e r ( ) . nr ( ) . va lu e ( ) ) −1 . 0 ) ) ∗(1 −M. boundaryField () [ patchi ] ) ) ; 26 27 s c a l a r F i e l d &nu = tnu . r e f () ; 28 29 for (++ i t e r ; i t e r != phases_ . end () ; ++i t e r ) 30 { 31 nu += i t e r ( ) ∗( i t e r () . boundaryField () [ patchi ] ∗i t e r () . nu0 ( ) . value() ∗exp ( ( c oef + (pow ( ( 1 . 0 −Mrho . boundaryField () [ pat ch i ] ) / (max(Mrho . boundaryField ( ) [ p at ch i ] , 0 . 1 ) ) , 0 . 48 ) ) ) ∗( 1 . 0 −Mrho . boundaryField () [ patchi ] ) ) ∗M. boundaryField () [ patchi ] + max( ViscMin , i t e r ( ) . boundaryField () [ patchi ] ∗0.25 ∗vi sc 2 ∗pow( ( i t e r () . Arrho0 () . value ( ) ∗0 .7 5 ) , −1.0 / i t e r () . nr ( ) . value ( ) ) ∗ exp ( i t e r () . Er () . value () / ( i t e r () . nr ( ) . value ( ) ∗R2 ∗(max (T_. boundaryField ( ) [ patchi ] , T02) ) ) ) ∗pow(strainRate . boundaryField ( ) [ p at ch i ] , ( 1 . 0 / i t e r ( ) . nr ( ) . va lu e ( ) ) − 1 . 0 ) ) ∗(1 −M. boundaryField ( ) [ patchi ] ) ) ; 32 } 33 return tnu ; 10.2 Correction of ∆Cp in Wittington’s law In /multiMeltInterFoamv2/multiMeltMixturev2/multiMeltMixturev2.H 1//- Return 2tmp<v ol Sc al a rF ie ld > Cpt () const ;//// 3 4//- Return the 5tmp<s c a l a r F i e ld > Cpt( const l a b e l patch i ) const ;//// In /multiMeltInterFoamv2/multiMeltMixturev2/multiMeltMixturev2.H 1seuilTCp_("seuilTCp" ,dimTemperature , 846.) , 2[...] 3// ----Cpt 4 5Foam : : tmp<Foam : : vo l Sc al ar F ie ld > 6Foam : : multiMeltMixturev2 : : CptInter ( ) const // radioactive term /CP(T) 32 7{ 8dimensionedScalar Adim("Adim" , dimensionSet (0 ,0 , 0 , −1, 0 , 0 , 0) , 1) ; 9tmp<v ol Sc al a rF ie ld > tCptInter = (min( double ( 1 . 0 ) , max( double ( 0 . 0 ) , ( (T_−i t e r () . Tsol () ) / ( i t e r () . Tliq () −i t e r ( ) . Tsol () ) ) ) ) ) ; //INit +CB 10 11 v o l S c a l a r F i e l d &CptInter = tCp tIn ter . r e f ( ) ; 12 13 for ( l a b e l c e l l I =0; c e l l I <mesh_ .C( ) . s i z e () ; c e l l I ++) 14 { 15 if (T_[ c e l l I ] < seuilTCp_ . value ( ) ) 16 { CptInter [ c e l l I ] =0;} 17 else 18 { CptInter [ c e l l I ] =1;} 19 } 20 CptInter . correctBoundaryConditions () ; 21 return tCptInter ; 22 } 23 24 Foam : : tmp<Foam : : vo l Sc al ar F ie ld > 25 Foam : : multiMeltMixturev2 : : Cpt ( ) const // radioactive term /CP(T) 26 { 27 const v o l S c a l a r F i e l d& CptInter= U_. mesh ( ) . lookupObject<v ol S c al ar Fi e ld >(" CptInter" ) ; 28 PtrDictionary<phase >:: c onst_iterat or i t e r = phases_ . begin ( ) ; 29 dimensionedScalar coef6("coef6" , dimensionSet(0 , 2, −2, −1, 0 , 0 , 0) , 199.5∗1 . e3 /221.78) ; 30 dimensionedScalar coef8("coef8" , dimensionSet(0 , 2, −2, −2, 0 , 0 , 0) , 0.0857∗1 . e3 /221.78) ; 31 dimensionedScalar coef9("coef9" , dimensionSet(0 , 2, −2, 1 , 0 , 0 , 0) , 5.0 e6∗1 . e3 /221.78) ; 32 33 dimensionedScalar coef61(" coef61 " , dimensionSet(0, 2, −2, −1, 0 , 0 , 0) , 229.32∗1 . e3 /221.78) ; 34 dimensionedScalar coef81(" coef81 " , dimensionSet(0, 2, −2, −2, 0 , 0 , 0) , 0.0323∗1 . e3 /221.78) ; 35 dimensionedScalar coef91(" coef91 " , dimensionSet(0, 2, −2, 1 , 0 , 0 , 0) , 47.9e−6∗1 . e3 /221.78) ; 36 37 return (coef6 + coef8 ∗T_ −coe f9 / (T_ ∗T_) )∗(1−CptInter ) + ( co ef61 + coef81 ∗T_−coef91 / (T_ ∗T_) ) ∗CptInter ; // tCpt; 38 } In /multiMeltInterFoamv2/TEqn.C 1{ 2fvScalarMatrix TEqn 3( 4Cpt∗(fvm : : ddt ( rhoC ,T) 5+ fvm : : div ( rhoPhi ,T) ) 6−fvm : : l a p l a c i a n ( kappa∗rhoC∗Cpt ,T) 7== Hr 8) ; 9TEqn . re l ax ( ) ; 10 fvOptions . co n st ra i n (TEqn) ; 11 TEqn . so lv e ( ) ; 12 fvOptions . c o r r e c t (T) ; In /multiMeltInterFoamv2/createFields.C : 33 1v o l S c a l a r F i e l d CptInter 2( 3IOobject 4( 5"CptInter" , 6runTime . timeName () , 7mesh , 8IOobject : : READ_IF_PRESENT, 9IOobject : :NO_WRITE 10 11 ) , 12 mixture . CptInter () 13 ) ; 14 CptInter . oldTime () ; 15 volScalarField Cpt 16 ( 17 IOobject 18 ( 19 "Cpt" , 20 runTime . timeName () , 21 mesh , 22 IOobject : : READ_IF_PRESENT, 23 IOobject : :NO_WRITE 24 ) , 25 mixture . Cpt ( ) 26 ) ; 27 Cpt . oldTime ( ) ; 10.3 Modification of the latent heat In /multiMeltInterFoamv2/multiMeltMixturev2/multiMeltMixturev2.H : 1tmp<v ol Sc al a rF ie ld > CA() const ;//// 2 3//- Return the 4tmp<s c a l a r F i e ld > CA( const l a b e l patc hi ) const ;//// In /multiMeltInterFoamv2/multiMeltMixturev2/multiMeltMixturev2.C : 1TM_( "TM" , dimTemperature , 8 7 3 . ) , 2[...] 3Foam : : tmp<Foam : : vo l Sc al ar F ie ld > 4Foam : : multiMeltMixturev2 : :CA() const 5{ 6volScalarField St2( 7IOobject( 8"St2" , 9mesh_ . time () . timeName () , 10 mesh_) , 11 mesh_ , 12 dimensionedScalar( 13 "St2" , 14 dimensionSet (0 , 0 , 0 , 0 , 0) , 15 1 . 0 ) ) ; 16 tmp<v ol Sc al a rF ie ld > tCA = 0+St2 ; //pour initialiser + CB 17 v o l S c a l a r F i e l d &CA = tCA . r e f ( ) ; 34 18 19 for ( l a b e l c e l l I =0; c e l l I <mesh_ .C( ) . s i z e () ; c e l l I ++) 20 { 21 if (T_[ c e l l I ] < TM_. value ( ) ) 22 {CA[ c e l l I ] = 1. ;} 23 else 24 {CA[ c e l l I ] =2;} 25 } 26 CA. correctBoundaryConditions () ; 27 return tCA ; 28 } In /multiMeltInterFoamv2/TEqn.C : 1{ 2fvScalarMatrix TEqn 3( 4CA∗Cpt∗(fvm : : ddt (rhoC ,T) 5+ fvm : : div ( rhoPhi ,T) ) 6−fvm : : l a p l a c i a n ( kappa∗rhoC∗Cpt ,T) 7== Hr 8) ; 9TEqn . re l ax ( ) ; 10 fvOptions . co n st ra i n (TEqn) ; 11 TEqn . so lv e ( ) ; 12 fvOptions . c o r r e c t (T) ; In /multiMeltInterFoamv2/createFields: 1volScalarField CA 2( 3IOobject 4( 5"CA" , 6runTime . timeName () , 7mesh , 8IOobject : : READ_IF_PRESENT, 9IOobject : :AUTO_WRITE 10 ) , 11 mixture .CA() 12 ) ; 13 CA. oldTime () ; In /multiMeltInterFoamv2/multiMeltMixturev2/multiMeltInterFoamv2.C : 1CA = mixture .CA() ; 10.4 Taking into account the variation of ρin the heat equation In /multiMeltInterFoamv2/TEqn.C : 1{ 2fvScalarMatrix TEqn 3( 4CA∗Cpt∗rho∗( fvm : : ddt ( 1 . 0 ,T) 35 5+ fvm : : div ( rhoPhi ,T) ) 6−fvm : : l a p l a c i a n ( kappa∗rho∗Cpt ,T) 7== Hr 8) ; 9TEqn . re l ax ( ) ; 10 fvOptions . co n st ra i n (TEqn) ; 11 TEqn . so lv e ( ) ; 12 fvOptions . c o r r e c t (T) ; 10.5 Modifications in the Test Case input files Boundary conditions in 0/U 1lowerWall 2{ 3type s l i p ; 4} 5atmosphere 6{ 7type fixedValue ; 8value uniform (0 0 0) ; 9} µ0for phase 1 (the two others are divided or multiplied by 10) in transportProperties 1nu0 nu0 [ 0 2 −1 0 0 0 0 ] 0 . 0 3 23 ; 11 Solver modifications (v4): Darcy-Brinkman-Stokes formulation - in process This section briefly describes the evolution of Aurélie’s solver towards solving Darcy-BrinkmanStokes equation, by adding an effective permeability term to the Stokes equation, and considering mineral phase changes and porosity change with two additional conservation equations. In process. 12 Solver modifications (v3): Krafla study - 2025 This section describes the modifications made for the Krafla study presented in [Gerbault et al., 2025] : modifications were made on the temperature-dependent melt relationship and thermo-mechanical properties in order to better account for rhyolite and basalt behavior in the context of Krafla’s magma sill intrusion. In this context, the dimensions of the 2D problem are of a few hundred meters and a few tens to a several hundred years. In practice, one must install the multiMeltInterFoamv3 solver and download the input files available on the GitLab repository and use the input files deposited in directory /Krafla_study, for a case simulating a basalt intrusion with cell sizes 0.45 m and minimum viscosity 2 m/s). 36 12.1 Equations and properties The physical problem is described with the same equations 1to 4and 8, but now, melt fraction Mis parametrized as : M= (1 + ai+biTr+ciT2 r+diT3 r)−1. Tr=T/Trefi, Tref = 1000oC. (15) Parameters in equation (15) were obtained by fitting simulations results of the crystallization of basaltic and rhyolitic magmas with typical Krafla composition [Borisova et al., 2023]. The thermal conductivity κand heat capacity Cpare prescribed as: Cp=Cp0+L∗ dM dT , κ =κ0 1 + κTT,(16) with L∗the latent heat of crystallization, κ0and κTconstants given in [Gerbault et al., 2025]. Dynamic viscosity µ=µmη(φ, )is a product of melt viscosity µmgiven by [Giordano et al., 2008] and a relative viscosity ηdue to the presence of crystals, that depends on the melt fraction M and the strain rate according to [Costa et al., 2009] : log µm=Ai+Ei T−Ci, η(φ, ) = 1+φδ [1−F(ϕ,,γ)]Bφ∗ F= (1 −ξ)·erf h√π 2·(1−ξ)ϕ·(1 + ϕγ)iwith ϕ=φ φ∗, φ = 1 −M. (17) Constants Ai, Ei, Cidiffer between rhyolite and basalt phases (i), parameters δ, φ∗, B, ξ, γ are strain-rate dependent. Further details are provided in Appendix 2 of [Gerbault et al., 2025], in revision. 12.2 Study summary and model setup In this study we reproduce the main heat transfer processes between a hot molten rhyolite (scenario 1) or basalt (scenario 2) intrusion and a colder rhyolite host rock domain at about 2 km depth, cf. Fig.11. We test the influence of mesh size and viscosity ranges on the results. The host rock domain in contact with the intrusion melts over a thickness that increases over time, enough to generate convective heat transfer and extremally sharp temperatures at the melting front, which simulate the sharp temperature jump encountered during the 2009 IDDP-1 drilling project at Krafla (Iceland). More details in [Gerbault et al., 2025]. Animations can be found on the Gitlab repository. The models setup is a 2D square domain of dimension H×Hmade of rhyolite crust (r subscripts parameters) at a uniform temperature of 400°C, in the middle of which a hot basalt intrusion (bsubscript parameters ) of thickness D and Temperature Tois injected. In scenario 1, the rhyolite intrusion has a thickness D= 300 m at temperature To= 980°C, and the domain’s size is H= 750 m. In scenario 2, the intrusion is composed of basalt of thickness D= 100 m at To= 1200°C and the domain’s size is H= 225 m. The top and bottom boundaries have a free-slip condition and temperatures are maintained there, fixed at 400°C. Velocities and temperatures at the lateral borders are set periodic. 37 Figure 11: Convective heat transfer between hot basalt and cold rhyolite at Krafla. Fig. 5 from [Gerbault et al., 2025]. 38 A Appendices A.1 Setting up the folder 0 Each variable file (U,T, ...) begins with a dimension line which informs of its unit : dimensions [kg m s K mol A cd ]; A 0 will correspond to the fact that the dimension is not present in the magnitude, a 1 that it is, a 2 that it is squared... Then, it is necessary to enter in the size file the boundary conditions which correspond to those of the file blockMeshDict : boundaryField { sides { type zeroGradient; } empty { type empty ; } } [of OpenFOAM, 2015a] Code Meaning type fixedValue ; value uniform valeur ; Impose a constant condition of type Dirichlet type uniformFixedValue ; uniformValue table ( (t0 valeur0 ) (t1 valeur1 ) ); Impose a time dependent variable condition type fixedValue ; value nonuniform List <scalar > ( valeur1 valeur2 ); Insert a field of values (fortran program) on chosen cells type zeroGradient ; Neumann condition : ∂grandeur ∂n = 0 type fixedFluxPressure ; value uniform valeur0 ; The pressure gradient is regulated in a way that the flow at the border corresponds to the velocity boundary condition Table 3: Table of the main boundary conditions used 39 Example : Velocity U 1/* ----------------------------*- C++ -*-------------------------------*\ 2========= | 3\\ / F ield | OpenFOAM: The Open Source CFD Toolbox 4\\ / O peration | Website: https :// openfoam .org 5\\ / A nd | Version : 6 6\\/ M anipulation | 7\*--------------------------------------------------------------------*/ 8FoamFile 9{ 10 ve rs io n 2 . 0 ; 11 format a s c i i ; 12 class volVectorField ; 13 location "0" ; 14 object U; // Velocity 15 } 16 //***********************************// 17 dimensions [0 1 −1 0 0 0 0 ] ; 18 i n t e r n a l F i e l d uniform (0 0 0) ; // Vitesse initiale dans le domaine 19 // Initial velocity inside the domain 20 boundaryField 21 { 22 le f t W a l l 23 { 24 type c y c l i c ; // Condition p r i o d i q u e 25 } 26 rightWall 27 { 28 type c y c l i c ; 29 } 30 lowerWall 31 { 32 type fixedValue ; 33 value uniform (0 0 0) ; // Adherence - No slip 34 } 35 atmosphere 36 { 37 type s l i p ; // Glissement - Slip 38 } 39 cote 40 { 41 type c y c l i c ; 42 } 43 cote1 44 { 45 type c y c l i c ; 46 } 47 } 48 // ********************************************************************* // 40 Sinusoidal interfaces - If for example we want a basal "water" layer with sinusoidal disturbance, we can use swak4Foam with the file funkySetFieldsDict, which looks like: 1defaultFieldValues 2( 3vol Sca lar Fie ldV alu e alpha . shamp 1 4vol Sca lar Fie ldV alu e alpha . water 0 5volScalarFieldValue alpha. oil 0 6) ; 7expressions 8( 9noise1 10 { 11 f i e l d alpha . water ; 12 expression "1" ; 13 condition "(pos () .y -0.01* sin (2*3.141* pos ().x /0.1) -0.1) <=0. " 14 keepPatches 1; 15 } 16 noise2 17 { 18 f i e l d alpha . shamp ; 19 ex pr es si on "0" ; 20 condition "(pos () .y -0.01* sin (2*3.141* pos ().x /0.25) -0.1) <=0. " // what is in water , mustbe removed in shamp to avoid phases superposition! 21 keepPatches 1; 22 } 23 . . . // repete noise ... with expression "0" and"1" for the condition >0 24 ) ; Please note, you must define the expression "0" or "1" for each of the 3 phases predefined by default, so in total there are 6 expressions (noise1...noise6, which are ad-hoc terms chosen by the user ; for example "1" for the condition <= 0 then "0" for the condition >0). Check the logical syntax sequence : expressions, noise, field, expression, condition, keepPatches. Warning 2, if we carry out specific tests on this geometry of the phases, remember to reinitiate the files 0/alpha*, copy the reference files (for example "cp CasTest/0.org/alpha* CeRun/0/."), then run the command : funkySetFields -time 0. 47 A.4 Postprocessing A.4.1 Paraview To display OpenFOAM data, type the command paraFoam, which launches paraview and creates a file..Openfoam. This file is deleted when Paraview is stopped. Some useful options: -paraFoam -builtin -paraFoam -touchAll (keeps the display file when paraview is stopped, which allows different cases to be displayed at the same time on the same paraview window) Main Paraview functions that can be used : -slice : to draw a 2D section -isovolume : to see an interval of specific values of a variable -plotOverLine : to plot a variable along a line at a given time. The resolution will allow more points to be recorded -plotSelectionOverTime : to plot the evolution of a variabl at a given point as a function of time (for this, we will activate the box ...) -Contour : to outline specific values (isotherms) If you want to follow the temporal and spatial evolution of a marker: -Point source : mettre les coordonnées du marqueur -ParticleTracer -plotGlobalVariableOverTime You can save the data in csv format (File-Save data). The option ’Precision’ indicates the number of significant digits to be saved. A.4.2 sampleDict You can also calculate quantities after the Run, therefore you must create a file sampleDict in /system and type the command sample. interpolationScheme cellPoint ; // choice of variable interpolation // cellPoint : linear interpolation based on cell size setFormat surfaceFormat ; // choice of 1D writing format : raw - ASCII format in column sets the name of // the directory in which the data will be written ( line_y // / write name of files { type uniform ;// uniform distribution of the 146 points in line axis y; // reference coordinate start ( 0.1 0 0.0085) ;// coordinate of the first point in line end ( 0.1 0.1 0.0085) ; // coordinate of the last point in line nPoints 146; // number of points in the line } ); fields (T); // list of variables to record on the lines //************************************************************************** // [of OpenFOAM, 2015g] 48 A.5 Reminder of basic Linux commands A.5.1 Terminal commands - Distant connexion Effect Command Connecter to a distant server ssh login@nom_serveur_distant Graphical interface utiliser X2go Modify a file ’ emacs ␣ file_name ’ to modify Transfer files scp fichier login@nom_serveur_distant :/ home / login Retrieve files scp login@nom_serveur_distant :/ home / login / fichiers . Compress zip nom . zip nom_a_zipper / Uncompress unzip nom.zip - files, folders, navigation Effect Command Go into a folder cd nom_dossier Rename a file mv nom_fichier nouveau_nom Create a folder mkdir dossier Copy a file cp nom_fichier nom_nouveau_fichier Copy a folder elsewhere into the present folder cp -r chemin / dossier_a_copier . Suppress a file or a folder \rm -r nom_dossier The & option sets the job in background (and frees the prompt shell). The option -r is for recursive action on the files and folders inside a folder. - Run and Jobs To interrupt a job or a run (if the run hasbeen sent on multiple processers, interrupting one processor (PID) will interrupt them all) : kill -9 nom_PID To display the running jobs : top . A.5.2 Aliases In order to simplify the general processing of files, we can create aliases or shortcuts to use in a shell terminal. Therefore these aliases are set either in the .bashrc file or in a specific .alias file. We first open the .bashrc file in an editor: gedit $HOME /. bashrc And we add the aliases inside it : alias para = ’ paraFoam ␣&’ alias block = ’blockMesh ’ alias check = ’checkMesh ’ alias set=’setFields ’ alias dec=’decomposePar’ alias rec=’reconstructPar’ alias doc=’cd␣ Documents ’ alias op=’cd ␣ OpenFOAM ’ 49 List of Figures 1 Model Setup, to simulate a partially melting crustal domain. ............ 3 2 Paraview : Test Case mesh and white background configuration. ......... 10 3 Paraview : Initial temperatuer field in 2D and along a vertical proifle. ...... 12 4 Paraview : Display the inclusions and the isotherms. ................ 13 5 Adapting the numerical time to the real time (here tN= 14900 s⇒tR= 4.7Ma)22 6Inclusions distribution at initial stage. .......................... 25 7Inclusions distribution, viscosity and isotherms after t = 4.8 Ma. ............ 25 8Profiles of temperature, depth and velocity for one marker, with paraview........ 25 9Vertical profiles of the average viscosity, temperature, with Moyenne.f95, for 2 time-steps. 25 10 Mean domain vertical velocity over time. ......................... 25 11 Convective heat transfer between hot basalt and cold rhyolite at Krafla. Fig. 5 from [Gerbault et al., 2025]. .................................. 38 12 Examples of grading .................................. 41 13 Nodes numbering in a 2D mesh. ............................ 41 14 Nodes numbering in a 3D mesh. ............................ 42 15 Connectivity of adjacent faceins of the mesh. .................... 42 50 References [Behrens, 2009]Behrens, T. (2009). Openfoam’s basic solvers for linear systems of equations. http://www.tfd.chalmers.se/~hani/kurser/OS_CFD_2008/TimBehrens/ tibeh-report-fin.pdf. [Borisova et al., 2023]Borisova, A. Y., Melnik, O. E., Gaborit, N., Bindeman, I. N., Traillou, T., Raffarin, M., Stefansson, A., Laurent, O., Leisen, M., et al. (2023). In situ probing of the present day zircon bearing magma chamber at Krafla, northeastern iceland. Frontiers in Earth Science, 11:1307303. [Costa et al., 2009]Costa, A., Caricchi, L., and Bagdassarov, N. (2009). A model for the rheology of particle-bearing suspensions and partially molten rocks. Geochemistry, Geophysics, Geosystems, 10(3). [Gerbault et al., 2025]Gerbault, M., Melnik, O., and Borisova, A. (2025). Numerical simulation of magma-rock interaction at krafla volcano using openfoam software and a simplified thermal model. EGUsphere, 2025:1–30. [Giordano et al., 2008]Giordano, D., Russell, J. K., and Dingwell, D. B. (2008). Viscosity of magmatic liquids: a model. Earth Planet. Sc. Lett., 271(1-4):123–134. [Louis-Napoléon et al., 2022]Louis-Napoléon, A., Bonometti, T., Gerbault, M., Martin, R., and Vanderhaeghe, O. (2022). Models of convection and segregation in heterogeneous partially molten crustal roots with a vof method–i: flow regimes. Geophysical Journal International, 229(3):2047–2080. [Louis-Napoléon et al., 2020a]Louis-Napoléon, A., Gerbault, M., Bonometti, T., Thieulot, C., Martin, R., and Vanderhaeghe, O. (2020a). 3-d numerical modelling of crustal polydiapirs with volume-of-fluid methods. Geophysical Journal International, 222(1):474–506. [Louis-Napoléon et al., 2020b]Louis-Napoléon, A., Gerbault, M., Bonometti, T., Thieulot, C., Martin, R., and Vanderhaeghe, O. (2020b). 3-d numerical modelling of crustal polydiapirs with Volume-Of-Fluid methods. Geophysical Journal International, 222(1):474–506. [of OpenFOAM, 2015a]of OpenFOAM, T. A. (2015a). Boundaries. http://cfd.direct/ openfoam/user-guide/boundaries/. [of OpenFOAM, 2015b]of OpenFOAM, T. A. (2015b). Mesh generation with the blockmesh utility. http://cfd.direct/openfoam/user-guide/blockmesh/#x25-1480005.3.2. [of OpenFOAM, 2015c]of OpenFOAM, T. A. (2015c). Numerical schemes. http://cfd.direct/ openfoam/user-guide/fvschemes/. [of OpenFOAM, 2015d]of OpenFOAM, T. A. (2015d). Running applications in parallel. http: //cfd.direct/openfoam/user-guide/running-applications-parallel/. [of OpenFOAM, 2015e]of OpenFOAM, T. A. (2015e). Solution and algorithm control. http: //cfd.direct/openfoam/user-guide/fvsolution/. [of OpenFOAM, 2015f]of OpenFOAM, T. A. (2015f). Time and data input/output control. http://cfd.direct/openfoam/user-guide/controldict/. [of OpenFOAM, 2015g]of OpenFOAM, T. A. (2015g). Tsampling data. http://cfd.direct/ openfoam/user-guide/sample/. 51 [Ulvrová et al., 2012]Ulvrová, M., Labrosse, S., Coltice, N., Råback, P., and Tackley, P. (2012). Numerical modelling of convection interacting with a melting and solidification front: Application to the thermal evolution of the basal magma ocean. Physics of the Earth and Planetary Interiors, 206:51–66. [Whittington et al., 2009]Whittington, A. G., Hofmeister, A. M., and Nabelek, P. I. (2009). Temperature dependent thermal diffusivity of the earth’s crust and implications for magmatism. Nature, 458(7236):319–321. 52