Using the transport model classes

Class hierarchy

BaseTransportModel1D is the abstract base class for all other classes. It includes functions i.e. to retrieve results, set up the model, etc., as far as they are not solver-specific.

Derived model classes are:

• EulerBackwardsModel1D, which uses the implicit Euler (first order) time evolution scheme and

• CrankNicolsonTransportModel1D, which builds on some routines of the former and implements the Crank-Nicolson scheme.

The classes in ensemble.py are wrappers to run multiple submodels of a given transport problem in parallel. A submodel consists of the transport model for a complete set of radionuclides, that is, a set of radionuclides in which no radionuclide decays into a radionuclide that is not included in the submodel

General remarks

The life cycle of transport model objects consists of instantiation/parametrization, running the model, and analyzing the results. Parametrization always happens at instantiation; i.e. parameters should not be changed after instantiation. Similarly, running the same model object twice will usually not work.

Instantiation/Parametrization

All model classes share, as far as possible, the same interface. For example, the EulerBackwardsModel1D has the following arguments:

def __init__(self,
             species,
             parameters_rock,
             parameters_transport_material=None,
             geometry=None,
             bcs=None,
             time_interval=1e6 * u.yr,
             options=None,
             parameters_source=None
             ):

species is just a list of the radionuclides considered in the model. The physical parameters of the model domain are provided via two arguments, parameters_rock and parameters_transport_material. The former is structured as a list of geological units, where each list entry is a dictionary. An example is shown below

import astropy.units as u
parameters_rock = [{'name': "unit1", 
                    'thickness': 1000 * u.m,
                    'porosity': 0.1,  
                    'head_gradient': 0.0, 
                    'hydraulic_conductivity': 2e-13 * u.m / u.s, 
                    'bulk_density': 2650 * u.kg / u.m**3}]

This defines a setup with only one layer, named unit1. Note the usage of the astropy.units package. We make use of it in many places. parameters_transport_material is another list of dictionaries, with each entry defining the transport parameters for one geological unit. An example is again shown below.

parameters_transport_material = [
                                  { 
                                   'I-129': 
                                     {'effective_diffusion_coefficient': 2.4e-12 * u.m**2 / u.s, 
                                      'sorption_coefficient': 0.0 * u.m**3 / u.kg},
                                    'Ca-41': 
                                     {'effective_diffusion_coefficient': 1e-11 * u.m**2 / u.s, 
                                      'sorption_coefficient': 6.3e-4 * u.m**3 / u.kg}
                                   }
                                 ]

In this case, sorption and diffusion coefficients are given for unit1, for the radionuclides I-129 and Ca-41. Both parameters_rock and parameters_transport_material need to have the same length.

geometry specifies the grid spacing and/or the grid. Internally, this is treated as an GeometryArgs object, but it can also be provided as a dict object, e.g.

geometry = {'dx':1*u.m}

This will generate a grid with a grid spacing of 1 m.

bcs will be handled as a BCArgs object, but can be provided a dictionary as well. If None is provided, Dirichlet boundary conditions at 0 mol/m^3 will be assigned per default.

time_interval defines the total time period for the simulation. options will be converted to an OptionsArgs object, but can again be provided as a dictionary. It can be used to set various options, like flags to include/not include specific processes, etc. parameters_source will evaluate to a SourceArgs object and can be provided as a dictionary. It controls how source terms are included in the calculation.

Running the model

A parametrized model can be run like this:

model.run(ics, t_eval)

where we have omitted some optional arguments. ics is a data structure that provides the initial conditions, usually a list or a dictionary. If it is a list, it has to provide the initial concentration vector for each species. If it is a dict, it has to provide the initial concentration vector for each species as key:value pair. All concentration vectors must have astropy units. Units can be kg/m^3 or mol/m^3.

t_eval is a vector of the times at which outputs of the simulation should be stored, and needs to have units.

Initial conditions need to specify the total mass of radionuclides per rock volume, including the sorbed fraction of radionuclides.

Analyzing the results

The calculated concentrations are always the total concentration of radionuclides (including any sorbed species), either in mass or amount of substance, per rock volume or per volume with the functions below. Internally, the concentrations are stored as total concentrations in mol/m^3 per fluid volume.

Access to the model results can be obtained by the functions get_curves_t, get_curves_x, and some other derived functions. The former allows accessing the spatial distribution of concentrations at a fixed time, whereas the latter returns the temporal evolution of the system at a fixed location. Some more complex routines are already implemented, e.g. plot_overview_c_x in ensemble.VSGEnsembleModel1D.

For example, the following code retrieves the spatial distribution of I-129 (as c), and the location of grid nodes (as x), at the last time stored:

x, c = model.get_curves_x("I-129", itime=-1, quantity="molar density", units=u.mol/u.l, per_rock_volume=True)