Coupling MEDUSA Sediment Model to iLOVECLIM

How to install and run iLOVECLIM including MEDUSA

!! This section is an update from April 2020. All the other sections below are anterior and have not been verified. !!

The version of iLOVECLIM including MEDUSA can be installed only with the new installation procedure. It is preferable that you extraxt iloveclim to a new directory, here called my_iloveclim_medusa.

svn co svn+ssh:// my_iloveclim_medusa/.

First, copy the MEDUSA directory with the corresponding libraries on obelix to the directory where you want MEDUSA to be (preferentially OUTSIDE the iloveclim main directory):

cp  -vpfR /home/acclimate/roche/medusa-svn /my_path/my_MEDUSA_directory 

Second, modify the path of the library in my_iloveclim_medusa/config/medusa.path to your /my_path/my_MEDUSA_directory

Then you can run a full install of iloveclim using the external libraries intel-fcomp/2013 netcdf/4

module purge
module load intel-fcomp/2013 netcdf/4

./ -E

If this works, you should theoretically be able to lauch a simulation with MEDUSA using the following flags in choixcomposantes.h, similarly to the Carbon Cycle CarbonCyclePage :

#define IMSK 1
#define COMATM 1
#define ROUTEAU 1
#define EVAPTRS 1
#define EVAPSI 1
#define CLAQUIN 0
#define F_PALAEO 0

All the flags related to GRISLI should be set to 0 for the moment.

#define CYCC 2
#define OCYCC 1
#define OLDC14 0
#define KC14 1
#define KC14P 0
#define OXNITREUX 0
#define WINDINCC 0
#define O2ATM 0
#define N2OATM 0
#define MEDUSA 1
#define PATH 1

Importantly, for the moment only the flag BATHY 0 is supported. Later development should include the possibility to use GRISLI and thus BATHY 1 or BATHY 2.

Then you can run your simulation with the options you want as usual, using the new command line ./bin/run-iloveclim.

!! This section is an update from April 2020. All the following sections below are anterior and have not been verified. !!


The new iLOVECLIM v1.1β coupled with MEDUSA sediment model is an extended version of the iLOVECLIM carbon cycle version presented and described in Bouttes et al. (2015) and described here < >. As established in this publication, the long-term objective is to study past and future carbon cycle changes over timescales of a few thousand to hundred of thousand years, typical of glacial–interglacial changes and, obviously, sedimentary processes are relevant to such timescales.

MEDUSA sediment model has been coupled to iLOVECLIM to allow the direct interaction between the ocean bottom layer and a sediment reactive layer, as well as to create new sediment layers or erode existing ones.

In our current version, fluxes of particulated organic matter (POC) and calcite are calculated and averaged every 10 years of simulation and used to contribute in the creation of sediment layers in MEDUSA every time step. At the same time, fluxes of dissolved DIC, DOCS, oxygen, alkalinity, nitrate and phosphate are exchanged between the bottom ocean layer and the surface reactive sediment layer every time step of the iLOVECLIM carbon cycle version. For the dissolved components, net fluxes are calculated and considered positive in the direction of the sediments, which means that, if the flux is positive, the flux must be substracted from the ocean and added to the sediment.

Due to the accumulation of particulated material on the bottom of the ocean, new sediment layers are generated. Therefore, the deeper sediment layers are not accessible any more by the ocean bottom for exchange. This "virtually lost" layers are used as LOOP BACK FLUXES which emulate the riverine and eolian compensation in the global ocean carbon cycle.

Note: iLOVECLIM - OCYCC uses DIC variable but exchange between MEDUSA and iLOVECLIM is done in terms of HCO3 , CO3 and CO2(dis)

MEDUSA Description

MEDUSA (Model of Early Diagenesis in the Upper Sediment) is a transient one-dimensional advection-diffusion-reaction model that describes the coupled early diagenesis processes of carbonates and organic matter in the surface sediment. MEDUSA sediment model has been developed by Guy Munhoven from the University of Liège (Belgium) < >

In summary,MEDUSA is a transient early diagenesis model originally fully coupled to a multibox model (MBM) of the carbon cycle, as described in Munhoven (2007). This model combines advection, diffusion and reaction processes to describe the early diagenesis processes of organic and inorganic (carbonates) carbon components. Ocean water columns interact directly with a reactive layer and accumulation generates an historical zone. Model includes loopback fluxes to keep a conservative mass balance and emulate riverine and eolian compensation.

Running a simulation with MEDUSA

Running a simulation of the full iLOVECLIM carbon cycle (OCYCC) with MEDUSA requires to configure the original iLOVECLIM carbon cycle model as well as several input files (namelists) for MEDUSA model. The iLOVECLIM version coupled with MEDUSA works in the same way that the original iLOVECLIM for Carbon Cycle except for the new input and configuration files for running MEDUSA.

  • Modify the file: lbm/sources/choixcomposantes.h to have the same options than for the Carbon Cycle procedures and also MEDUSA 1, please check < > for more details. Do not forget to establish your initial conditions for restarting using, for example, -s 50800 and -S data/default/ic050800 (to start from the initial conditions)
  • Prepare your restarting files for MEDUSA using reaclay, flx and sedcore files from a previous simulation or, in case you don't want to use any initial condition for MEDUSA, set the corresponding options of medusa_filelist.nml file to "/dev/null"
  • Configure your MEDUSA input files and namelists in medusa-svn/inputdata-cpl: medusa.rrp, medusa.tsi, medusa_filelist.nml, medusa_seafloor_init.nml
    • It is specially important to pay attention to medusa_filelist.nml, cfn_ncin_sedcore file is the main file for doing a restart in the MEDUSA simulation. In MEDUSA, it is possible to read this input file and create a different output file cfn_ncout_sedcore file. However, there is a small bug and this option does not work in our current version. Therefore, the procedure suggested is to make a copy of the original sedcore input nectdf file, re-name it and use this file as the same file for input and output (that works perfectly!)

Here you can find an example of medusa_filelist.nml

cfn_nmlin_init     = "medusa_seafloor_init.nml"
cfn_ncin_init      = "/MEDUSA_RESTART_DIR/"
cfn_ncin_flx       = "/MEDUSA_RESTART_DIR/"
cfn_ncin_sedcore   = "/MEDUSA_RESTART_DIR/"

cfn_ncout_reaclay  = "outputdata/sediments/"
cfn_ncout_reaction = "outputdata/sediments/"
cfn_ncout_bc       = "outputdata/sediments/"
cfn_ncout_flx      = "outputdata/sediments/"
cfn_ncout_sedcore  = "/MEDUSA_RESTART_DIR/"
ctitle_ncfiles     = "Short description"
cfn_ncout_aux      = "outputdata/sediments/"

The other important configuration file is medusa.tsi

! initial time of the simulation
ayears_medusaini = 0.0D+00
! Medusa time step length in days
ndays_medusastep = 3600
! Frequency of NetCDF output (every nn Medusa time steps)
nmedusasteps_ncoutfrequency = 10

which allows to set the number of days of the MEDUSA time step and the frequency of the NETCDF Output

Finally the current default files form medusa.rrp and medusa_seafloor_init.nml are:

k_1   = 365.25D+00
n_1   = 4.5D+00

k_2   = 0.06D+00
khs_2 = 0.020D+00

! Template for medusa_seafloor_init.nml, automatically generated by
! CREATE_MOD_SEAFLOOR_INIT from the MEDUSA configuration
! utility MedusaCoCoGen.
! Solids'  (prefix "solid_") are expected to be given in %dry-weight
! Porewater solutes' (prefix "porew_") are expected to be given in mol/m3
solid_pcent_clay =  100.00D+00
solid_pcent_calc =    0.00D+00
porew_molm3_co3 =    81.76D-03
porew_molm3_hco3 = 2307.01D-03
porew_molm3_co2 =    31.23D-03
solid_pcent_om =      0.00D+00
porew_molm3_o2 =    100.00D-03

Coupling description

Main aspects of general MEDUSA Coupling are described in detail in MEDUSA-- Guide to Coupling by Guy Munhoven. In this section, we will describe the main aspects of the coupling of MEDUSA specifically with iLOVECLIM. Most of the Fortran code for the coupling is contained in three module files: medusa_wrapper.f90, flux_from_sediments_mod.f90 and flux_to_sediments_mod.f90.

The last two modules contain the fuctions and subroutines for averaging fluxes of particules from iLOVECLIM to MEDUSA, the lines of code for the exchange of fluxes between OCYCC (ocean carbon cycle model) and MEDUSA (sediments) and the lines that modify the corresponding state variables in iLOVECLIM (e.g., DIC, PO4, NO3, O2) due to the effect of sediments.

We will describe more in detail fluxes, loopback fluxes and variables in sediments in the current version, as well as, the files and modules in the code.

<Link to the two documents describing the coupling and compilation of MEDUSA>

Units in MEDUSA

Particle fluxes (organic matter and calcite) need to be expressed in MEDUSA as [kg/m2/yr]. The conversion is performed in mod_iloveclim_o2s.F assuming that the units of the particle fluxes coming from OCYCC are [molC/m2/yr]. In flux_to_sediments_mod.f90, particle fluxes are passed to MEDUSA in the variables TPP_mafond and caco3_mafond, which are obtained averaging TPP_ma and caco3_ma, expressed in [molC/m2/day] and multiplying by 360 days to get the values per year.

TPP_ma is calculated in line 105 of maphot.f taking TPP_m in Tmol and dividing this value by SQRO2 (area in m2) and then dividing the result by a factor of 1.0d-12 which precisely would convert TmolC/m2/year into molC/m2/year, as requested by MEDUSA code.

Same procedure is repeated in line 262 in maphot.f with caco3_ma.

More details about units and scale factors in the Carbon Cycle Module can be found in:

< >

Future work

Including specific aspects of isotopes in MEDUSA.

Last modified 18 months ago Last modified on 04/28/20 18:05:43