Tutorial

Overview

Sutra is a python package to calculate the advective removal of microbial organisms (also called ‘pathogens’) from source to end_point.

Main features:
  • Includes database of removal parameters for microbial organisms.

  • Calculate the removal and concentration of the microbial organism over distance, and with time

The code is based on the equations given in chapter 6.7 of BTO 2012.015 “Gezondheidsrisico’s van fecale verontreiniging” [in Dutch], pp 71-74. The equations used are described below.

Background information

During transport in the subsurface, microbial organism removal takes place by both attachment to the soil matrix and by inactivation. The virus concentration ‘C’ [m-3] through steady-state transport of microbial organisms along pathlines in the saturated groundwater can be approximated by:

Equation1: dC/dx + ((k_att + μ_1) / v) * C = 0

Where k_att + mu1 –> equal to removal rate ‘lambda’ ‘k_att’: attachment coefficient [day-1] ‘mu1’: inactivation coefficient [day-1] x: the distance traveled [m] v: the porewater velocity [m day-1] or ‘darcyflux divided by the effective porosity’

Assuming that the background concentration of the relevant microbial organism equals 0, the relative removal ‘C_x/C0’ can be calculated as follows.

Equation2: log(C_x/C_0) = -(((k_att + μ_1))/ln⁡(10)) * (x/v)

The attachment coefficient ‘k_att’ depends on the effective porosity ‘epsilon’, the grain diameter of the sediment ‘d_c’, ‘sticky coefficient’ alpha [day-1], the porosity dependent Happel’s parameter ‘A_s’, diffusion constant ‘D_BM’ [m2 day-1], and the porewater velocity [m day-1].

Equation3: k_att = 3/2 *(((1-ε))/d_c) * α*4*A_s^(1⁄3)*(D_BM/(d_c*ε*v))^(2⁄3) * v

The sticky coefficient alpha is determined by coefficient ‘alpha0’, which depends on both the soil type and the type of organism. Alpha0 is being determined for a reference pH [pH0], e.g. pH=7.5. Alpha relates to alpha0 as follows [corrected for different pH].

Equation4: α = α_0 * 0.9^(((pH - pH_0) / 0.1))

The other parameters are calculated as follows:

Happel’s porosity dependent parameter

Equation5: A_s = 2 * ((1 - γ^5))/((2 - 3γ + 3γ^5 - 2γ^6))

Where:

Equation6: γ = (1 - ε)^(1⁄3)

Boltzmann diffusion coefficient:

Equation7: D_BM = (K_B * (T + 273))/(3π * d_p * μ) * 86400

with Boltzmann constant K_B [1,38 × 10-23 J K-1], organism diameter d_p [m], water temperature T [degr C], and conversion factor 86,400 [s day-1].

The dynamic viscosity ‘mu’ [kg m-1 s-1] depends on the groundwater density ‘rho’. The water density is assumed to be 999.7 [kg m-3], representative for fresh groundwater in the Netherlands under a reference temperature of 12 degrees centigrade.

Equation8: μ = (ρ * 497*10^(-6))/(T + 42.5)^(3⁄2)

Steps

Operating the microbial organism removal involves 2 steps:

#. Run/load the removal_functions.MicrobialRemoval class to retrieve the default microbial (removal) parameters, if present in the database. Otherwise, an empty dataframe is returned. #. Run removal_functions.calc_advective_microbial_removal to calculate the final concentration after a distance and time traveled.

Now, let’s try some examples. First we import the necessary python packages

In [1]: import numpy as np

In [2]: import pandas as pd

In [3]: import os

In [4]: import sys

In [5]: from pathlib import Path

In [6]: import sutra.removal_functions as rf

Scenario A: Calculate removal of a microbial organism using default database parameters.

## Default removal parameters ##
In [7]: organism_name = "carotovorum"

# Redox condition: 3 options ['deeply_anoxic','anoxic','suboxic']
In [8]: redox_cond = 'anoxic'

# organism diameter [m]
In [9]: organism_diam = 1.803e-6

# Starting concentration
In [10]: conc_start = 1.

# Ambient groundwater concentration
In [11]: conc_gw = 0.

# effective porosity
In [12]: por_eff = 0.33

# Sediment grainsize
In [13]: grainsize = 0.00025

# pH of the groundwater
In [14]: pH_water = 7.5

# Water temperature
In [15]: temp_water = 10.

# Water density [kg m-3]
In [16]: rho_water = 999.703

# Distance traveled along pathline [m]
In [17]: distance_traveled = 100.

# Time traveled [days]
In [18]: traveltime = 1.

# Porewater velocity [m day-1]
In [19]: porewater_velocity = distance_traveled / traveltime

First initialize a class for calculating the removal of an organism.

In [20]: mbo_removal_scenA = rf.MicrobialRemoval(organism = organism_name)

In [21]: removal_parameters = mbo_removal_scenA.removal_parameters

# Return the (default) removal parameter values
In [22]: print(removal_parameters)
{'organism_name': 'carotovorum', 'alpha0': {'suboxic': 0.3, 'anoxic': 0.577, 'deeply_anoxic': 0.577}, 'pH0': {'suboxic': 7.5, 'anoxic': 7.5, 'deeply_anoxic': 7.5}, 'organism_diam': 1.803e-06, 'mu1': {'suboxic': 1.2664, 'anoxic': 0.1279, 'deeply_anoxic': 0.1279}}

Calculate final concentration after advective microbial removal

# Calculate final concentration and print it
In [23]: C_final_default = mbo_removal_scenA.calc_advective_microbial_removal(grainsize = grainsize,
   ....:                                     temp_water = temp_water, rho_water = rho_water,
   ....:                                     pH = pH_water, por_eff = por_eff,
   ....:                                     conc_start = conc_start, conc_gw = conc_gw,
   ....:                                     redox = redox_cond,
   ....:                                     distance_traveled = distance_traveled,
   ....:                                     traveltime = traveltime)
   ....: 

In [24]: print(C_final_default)
8.2065781569924e-12

# Print lambda (default): removal rate [day-1]
In [25]: print(mbo_removal_scenA.lamda)
25.526085068992856

Scenario B1: Manual input of removal parameters, not included in default database

## Removal parmeters ##
# Organism name
In [26]: organism_name = "MS2"

# Redox condition: 3 options ['deeply_anoxic','anoxic','suboxic']
In [27]: redox_cond = 'anoxic'

# alpha0: 'sticky coefficient' [-]
In [28]: alpha0 = 0.001

# Reference pH for calculating 'alpha' [-]
# --> if pH == pH0, then collision efficiency alpha equals the value of alpha0
In [29]: pH0 = 7.5

# time dependent inactivation coefficient 'mu1' [day-1]
In [30]: mu1 = 0.149

# organism diameter [m]
In [31]: organism_diam = 2.33e-8

Add remaining ‘ambient’ input parameters

# effective porosity
In [32]: por_eff = 0.33

# Sediment grainsize
In [33]: grainsize = 0.00025

# pH of the groundwater
In [34]: pH_water = 7.5

# Water temperature
In [35]: temp_water = 10.

# Water density [kg m^-3]
In [36]: rho_water = 999.703

# Distance traveled along pathline [m]
In [37]: distance_traveled = 100.

# Time traveled [days]
In [38]: traveltime = 1.

# Porewater velocity [m day-1]
In [39]: porewater_velocity = distance_traveled / traveltime

# Starting concentration
In [40]: conc_start = 1.

# Ambient groundwater concentration
In [41]: conc_gw = 0.

Initialize a class for calculating the removal of an organism

In [42]: mbo_removal_B1 = rf.MicrobialRemoval(organism = organism_name)

Calculate (relative) concentration following advective microbial removal

In [43]: C_final_B1 = mbo_removal_B1.calc_advective_microbial_removal(grainsize = grainsize,
   ....:                                     temp_water = temp_water, rho_water = rho_water,
   ....:                                     pH = pH_water, por_eff = por_eff,
   ....:                                     conc_start = conc_start, conc_gw = conc_gw,
   ....:                                     redox = redox_cond,
   ....:                                     distance_traveled = distance_traveled,
   ....:                                     traveltime = traveltime,
   ....:                                     organism_diam = organism_diam,
   ....:                                     mu1 = mu1,
   ....:                                     alpha0 = alpha0,
   ....:                                     pH0 = pH0 )
   ....: 

# print final concentration
In [44]: print(C_final_B1)
0.38739172625173746

Print the attachment coefficient ‘k_att’ and removal rate ‘lambda’

# k_att, calculated
In [45]: print(mbo_removal_B1.k_att)
0.7993188853572424

# lambda, calculated
In [46]: print(mbo_removal_B1.lamda)
0.9483188853572424

Scenario B2: An alternative way to enter removal parameters and calculate the final concentration Should compare to previous input, be aware to enter the correct redox related values for ‘anoxic’ situation

In [47]: mbo_removal_B2 = rf.MicrobialRemoval(organism = organism_name,
   ....:             alpha0_suboxic=None,
   ....:             alpha0_anoxic=0.001,
   ....:             alpha0_deeply_anoxic=None,
   ....:             pH0_suboxic=None,
   ....:             pH0_anoxic=pH0,
   ....:             pH0_deeply_anoxic=None,
   ....:             mu1_suboxic=None,
   ....:             mu1_anoxic=mu1,
   ....:             mu1_deeply_anoxic=None,
   ....:             organism_diam=organism_diam,
   ....:             )
   ....: 

Read removal parameters from ‘removal_parameters’. Check these values as follows

In [48]: removal_parameters = mbo_removal_B2.removal_parameters

In [49]: print(removal_parameters)
{'organism_name': 'MS2', 'alpha0': {'suboxic': None, 'anoxic': 0.001, 'deeply_anoxic': None}, 'pH0': {'suboxic': None, 'anoxic': 7.5, 'deeply_anoxic': None}, 'organism_diam': 2.33e-08, 'mu1': {'suboxic': None, 'anoxic': 0.149, 'deeply_anoxic': None}}

Calculate the final concentration, removal parameters for redox condition ‘anoxic’ (given by ‘redox_cond’)

# Only include 'ambient'/'physical' parameters (removal parameters loaded preiously)
In [50]: C_final_B2 = mbo_removal_B2.calc_advective_microbial_removal(grainsize = grainsize,
   ....:                                     temp_water = temp_water, rho_water = rho_water,
   ....:                                     pH = pH_water, por_eff = por_eff,
   ....:                                     conc_start = conc_start, conc_gw = conc_gw,
   ....:                                     redox = redox_cond,
   ....:                                     distance_traveled = distance_traveled,
   ....:                                     traveltime = traveltime)
   ....: 

# k_att, calculated
In [51]: k_att = mbo_removal_B2.k_att

In [52]: print(k_att)
0.7993188853572424

# lambda, calculated
In [53]: lamda = mbo_removal_B2.lamda

In [54]: print(lamda)
0.9483188853572424

# Final concentration
In [55]: print(C_final_B2)
0.38739172625173746

Notice that the output concentrations C_final_B1 & C_final_B2 are equal

In [56]: print(C_final_B1 == C_final_B2)
True