dynamic-R/RTM
Reaction Transport Modelling: Course Tutorials and Exercises

Package index
Search the dynamic-R/RTM package
Vignettes
Functions
8
Source code
1
Man pages
0
Home
/
GitHub
/
dynamic-R/RTM
/

In dynamic-R/RTM: Reaction Transport Modelling: Course Tutorials and Exercises 

knitr::opts_chunk$set(echo = TRUE)

library(ReacTran)
library(marelac)

If you are a total beginner or need to get started from scratch

You will need to install couple of programs on your computer to be able to do all the learning, modelling and reporting work in the course. The software we will use is free of charge and used by millions of people all around the world. In this section we will explain --- step by step --- how to get started with the installation to make sure that everything works on your computer as required.

We assume here that you use Microsoft Windows. The sequence of steps for Linux or Mac-OS users is the same, although their execution may differ in specific details.

Core R software

R is the programming language that we will use in the course. It can be downloaded from \url{http://www.r-project.org/}. Click on the download R link, choose a mirror, and download the precompiled binary distribution for your system. On this website, you will also find useful documentation.

First, install the base R (e.g., installation file R-4.1.0-win.exe or some later version). Subsequently, install Rtools (e.g., installation file rtools40v2-x86_64.exe). The later step will allow you to compile R-packages on your computer, which will be quite important later on. Ensure that you complete the instructions described in the section Putting Rtools on the PATH after you have installed Rtools.

Rstudio

To execute R code in a nice and comprehensive environment, we will use Rstudio. Go to \url{http://rstudio.org}, choose Download, choose Rstudio Desktop Free, download the installation package for your system (e.g., RStudio-1.4.1717.exe or some later version), and install it.

Open Rstudio and test whether you can compile an R-package from source, as described in the instructions for installing Rtools (see above).

R packages

Many of the useful functions are not part of the R base program, but are available as R packages. They are not necessarily made by the R core-development team, but are developed and shared with the R community by experts in other fields, e.g., in the field of modelling.

To install an R-package in Rstudio, click on the Packages tab (lower-right window), click on install, type the package name, make sure you have checked install dependencies, and click install. When asked about compilation for packages that need to be installed from sources, click yes. Be patient, the package compilation step may take a while.

In this course, we will use several packages. We recommend that you install them in this order:

devtools::install_github("dynamic-R/RTM", depend=TRUE)

More information about the RTM package is provided in Section 2 of this document.

LaTeX distribution

We will use the typesetting program \LaTeX\ to produce beautifully typeset documents. This program is extremely powerful and versatile, and it is widely used especially by people who want to write equations in their documents. In this course we will use it to convert documents written in markdown into PDF. To install \LaTeX, go to \url{https://miktex.org/}, choose Download, download the installation file (e.g., basic-miktex-21.6-x64.exe), and install it.\footnote{The MacTex distribution (\url{https://tug.org/mactex/}) is recommended for Mac-OS users, the TexLive distribution is recommended for Linux users. Follow the installation instructions on the respective websites.} During the installation process, select yes when asked Install missing packages on-the-fly. Also, run the MikTeX console from the star menu just after the installation has finished, select check updates, and update any packages suggested by this check. At the end, restart Rstudio if it was running before \LaTeX\ installation.

Check whether everything is working

To check whether everything is working at this point, open Rstudio, and create a new R Markdown file. This is done from the menu by selecting File $\rightarrow$ New file $\rightarrow$ R Markdown $\rightarrow$ Document. Type the title, your name in the author box, check the PDF checkbox, and click Ok. If everything goes as expected, your Rstudio session should look similar as shown in Figure 1.

Screenshot of Rstudio with the default Rmd file.{width=98%}

The top-left panel contains the editor, where you will write your text and R-code. In the top-right panel, you can view several things, including the Environment (which shows values of your R-variables) and Tutorials. The bottom-left panel contains the R-console, where you can type R-commands, and the R Markdown console, where you can see the progress of your Rmd file compilation (see next). Finally, the bottom-right panel allows you to manage your R-packages, view help files, or see plots.

Once you create a new Rmd file, it is best to first save it under a name of your choice (e.g., test1.Rmd). Then, after you are satisfied with its content, click on the Knit button in the top part of the editor panel (useful short-cut: press Ctrl+Shift+K). This will start a compilation process, which will convert the Rmd file into a nicely typeset document in one of three possible formats: PDF, DOC, or HTML (Figure 2). Read the following section if you experience problems during compilation of your Rmd file into PDF.

Nicely formatted output produced from the Rmd file. Left: PDF, viewed in SumatraPDF. Middle: DOC, viewed in Microsoft Word. Right: HTML, viewed in Firefox.

Possible issues

Note that if you are converting the Rmd file (e.g., \texttt{test1.Rmd}) into PDF for the very first time, the compilation may take a while or it may not work at all. This is because some packages required by the knitting process are missing and need to be installed on your system. Try the following steps to debug this issue.

If you have selected \emph{yes} when asked \emph{Install missing packages on-the-fly} during \emph{MikTeX} installation, the compilation into PDF should finish at some point; you only need to be patient.

If, however, you did not select \emph{yes}, you need to stop the knitting process by clicking the red stop button and fix things manually. To do this, open your file browser and navigate to the folder where you saved the \texttt{test1.Rmd} file. If the knitting of the file \texttt{test1.Rmd} did not succeed, this folder should contain a file called \texttt{test1.tex}. Open it in the native \TeX\ editor that comes with the \emph{MikTeX} distribution (most likely \texttt{TeXworks}) and compile it with pdf\LaTeX\ from there.\footnote{To do this, select pdf\LaTeX\ in the drop-down menu at the top and click on the green arrow.} In this process, you should be prompted to install missing \LaTeX\ packages. Be patient and wait until the missing packages are installed and the compilation process is finished.\footnote{At the time of writing this Reader, it took about 2 min until the first-time knitting and \LaTeX\ compilation were complete.} It may happen that errors due to missing files such as geometry.sty or fancyvrb.sty occur. Under such circumstances, run the MikTeX console and install these packages from there.

After completing the debugging steps above, everything should work well and you should find a file test1.pdf in the same folder as the file test1.Rmd. You can view it in any PDF reader. Note, however, that \emph{Acrobat Reader}, which is a commonly used PDF viewer and may be set as the default PDF viewer on your system, may not be the best choice. This is because the PDF file is \emph{locked} for writing when it is open by \emph{Acrobat Reader}. This means that if you want to re-knit your Rmd file into a PDF file, you may end up with a writing error rather than an updated PDF file. To avoid this error, you will need to close the PDF file in \emph{Acrobat Reader} before clicking the \emph{knit} button in \emph{Rstudio}. We recommend a better solution, though: use a PDF viewer that does \emph{not} lock the viewed PDF file. One such program is called \emph{SumatraPDF}, which can be downloaded from \url{http://www.sumatrapdfreader.org}. Others include \emph{Chrome} or \emph{Edge}.

Learning reaction-transport modelling using the package RTM

Now that you have installed the required software, you can proceed by studying other material prepared for the reaction-transport modelling course. This material is put together in the R-package called RTM. You can install this package in Rstudio by entering in the R-console:

devtools::install_github("dynamic-R/RTM", depend=TRUE)

The RTM package contains tutorials, exercises and extra readers. This division reflects the main course activities, as described in this section.

RTM tutorials

The tutorials contain short videos about a specific modelling topic. The topics range from basic principles (e.g., conceptual diagrams, mass balance equations) to more advanced (e.g., rate laws) and very advanced topics (e.g., pH modelling, reaction-transport modelling in porous media). The tutorials also contain simple tasks and exercises, which allow you to test your knowledge gained from the videos and prepare you for the exam.

We recommend that you study the tutorials in advance, i.e., before the class. You try to learn as much as possible by watching the short videos (once, but preferably multiple times), consult the textbook, and solve the tasks and exercises within the individual tutorials. In the process, you note questions about topics that you have not fully understood. You will be able to consult these questions with your peers and with the lecturers during the class.

If you installed the RTM package, you will find the tutorials in the Tutorial tab of the Rstudio program (Figure 3). Run the tutorial by clicking the Start Tutorial button.

Left: Rstudio with the Tutorial tab selected. Right: RTM tutorial *conceptual* run as a web-based app, created using the *learnr* package.

Alternatively, enter the following command in the R-console to see the complete list of available tutorials (output shown in table below):

RTMtutorial("?")

#RTM:::RTMtutorial("?")
knitr::kable(RTM:::RTMtutorial("?"))

Then, choose a tutorial you want to view (e.g., conceptual), and start it by entering the following command in the R-console:

RTMtutorial("conceptual")

Note that the tutorials are sorted from basics to intermediate and advanced topics. By studying the tutorials in the same order, you will be able to build up your knowledge of modelling in R in a gradual way. We will adhere to the same order during the course.

RTM exercises

Each tutorial comes with a set of exercises, which contain more complicated problems that you should be able to solve once you have learned the basics from the tutorial. These exercises will be solved during the class. This means that you will be able to consult your questions, ideas, and solutions with the lecturers in a direct and interactive way.

The exercises recommended for a given tutorial are listed at the end of the tutorial (section More exercises). For example, the tutorial on Reversible reactions and equilibrium chemistry recommends three additional exercises (equilibriumNH3, equilibriumHCO3, equilibriumOMD) (Figure 4).

Screenshot of an RTM tutorial *Reversible reactions and equilibrium chemistry*, showing the suggestions for extra exercises to solidify the knowledge gained in the tutorial.

Enter the following command in the R-console to see the complete list of available exercises (output shown in table below):

RTMexercise("?")

#RTM:::RTMexercise("?")
knitr::kable(RTM:::RTMexercise("?"))

Then, choose an exercise you want to solve (e.g., diagenesis) and generate the corresponding PDF file by entering the following command in the R-console:

RTMexercise("diagenesis", output="PDF")

If you want to export the exercise in a HTML or DOC format, change the value of the parameter output to HTML or DOC, respectively.

RTM readers

Finally, the RTM package contains extra readers. These readers explain topics that are not part of the course but are nevertheless useful in the context of reaction-transport modelling. They explain, e.g., how to visualise multi-dimensional data generated by the models, implement events and forcing functions based on data, fit data with a model, create interactive model applications, analyse the response times of perturbed systems, or model pH gradients.

Enter the following command in the R-console to see the complete list of available readers (output shown in table below):

RTMreader("?")

knitr::kable(RTM:::RTMreader("?"))

To generate a PDF file of the reader (e.g., visualisation), enter the following command in the R-console:

RTMreader("visualisation", output="PDF")

Scientific programming practice

Scientific research often involves computation, e.g., to convert numbers to different units, to calculate derived quantities, to perform statistical analyses on data, or to create and solve models. More and more of these computations are done by writing computer codes, or scripts. These scripts consist of a sequence of commands that tell the computer what it should do.

It is a good programming practice to write computer codes that are easy to understand. This not only facilitates exchange of these codes with other people (e.g., colleagues, your supervisor or lecturer), but it will also give you a head-start if you want to use these codes in the future. In addition, a well-documented code is an essential step to reproducible and reusable science.

The first coding attempt in R

Consider the following R-code:

head(airquality,n=2)
par(mfrow=c(1,2))
plot(airquality$Solar.R,type="l",xlab="time",ylab="Solar.R")
plot(airquality$Wind,type="l",xlab="time",ylab="Wind")

While we can more or less guess what this code does, it is not clear where the data come from and what the purpose of this code is. Also, the statements themselves are not extraordinarily legible.

Adding structure and comments in the code

Legibility of the code can easily be increased by adding a visual structure to the code: alignment of comparable sections, use of spaces, etc.

One way to document the above code is to add comments. In R, comments are preceded by a hash-sign (#).

The following code is self explanatory. If we execute this code, then the first two lines of the dataset will be printed to the console, and the graphs will be created and included in a figure within the document (Figure 5). However, the output will be separate from the code that generated it.

# The R dataset "airquality" contains daily air quality measurements in New York,
# from May to September 1973. 
# In the code below, we first look at the first two rows of this dataset
# and then plot the solar radiation and wind data, in two figures next to one another.
# Source code written by Karline Soetaert.

head(airquality, n = 2)  # show first two lines of the dataset

par(mfrow = c(1, 2))     # figures aligned in one row, two columns
plot(airquality$Solar.R, type = "l", xlab = "time", ylab = "Solar.R")
plot(airquality$Wind,    type = "l", xlab = "time", ylab = "Wind"   )

R Markdown --- creating fully integrated documents

Even better is to have the text, the code and the output in one document. This is what R Markdown does. When the final document is "created", the R-code is executed and the results are merged in the document. Figure 6 illustrates how this works in Rstudio.

Screenshot of *Rstudio* showing an R Markdown code (left) and the resulting document produced via *knit* (right).

There are several parts in the R Markdown document:

Finally, the Knit button (encircled in blue in Figure 6) will compile the R Markdown document and generate the document on the right. This final document contains both the R-code and the results. It fully documents the analysis.

Getting started with R Markdown

The easiest way to get started is to create a first R Markdown document using Rstudio. From the menu choose File $\rightarrow$ new File $\rightarrow$ R markdown. This will open an R Markdown file that already contains some R code and text. Do this, and look at the contents of this document.

Two useful information sources are accessible from the Rstudio menu: help/cheatsheets/R markdown Cheat Sheet and help/cheatsheets/R markdown Reference Guide.

Short help on the text format that can be used in Markdown documents can be found in help/Markdown quick reference.

Exercises --- Creating and exporting an Rmd document

Open the R Markdown file called Rmarkdown_small.Rmd in Rstudio. The file is available via Blackboard and was used to generate Figure 5.

A quick overview of R and Rstudio

An R-code is easy to read once you realise that:

Console versus R-scripts versus R Markdown

There are three ways in which to work with R.

Console

We can type commands into the R console window at the command prompt (>) and use R as a powerful scientific calculator. For instance, enter in the console window:

 pi*0.795^2 ; 25*6/sqrt(67) ; log(25) ; log10(25)

Here sqrt, log and log10 are built-in functions in R, pi is a built-in constant, and the semi-colon (;) is used to separate R-commands entered on one line.

In the console window, the "UP" and "DOWN" arrow keys can be used to navigate through previously typed commands (command history). Note that typing commands in the console is not the most efficient way of using R, especially if you want to debug and save your work for later use (see next).

Throughout these notes, the following convention is used:

```{R, eval = FALSE }

3/2

  denotes an R-code (``>`` is the prompt), and 

```{R, eval = FALSE }
[1] 1

is an R output, as written in the console window.

R-scripts

A better way of using R is by creating R-scripts in Rstudio's editor and save them in a file (e.g., "MyModel1.R") for later re-use.

R-scripts are sequences of R-commands and expressions. These scripts should be submitted to R before they are executed. This can be done in several ways:

```{R, eval = FALSE } source ("MyModel1.R") 

*   by opening the file in a text editor, copying the R-script to the clipboard (ctrl+C) and pasting it (ctrl+V) into the R-console;
* if you use *Rstudio*, which we recommend, you can either execute the current line (Ctrl+Enter), a section, or the entire file (Ctrl+S).

### R Markdown

The most integrated way to work with R is by using *R Markdown*. This is what we will do in the modelling course. This means that you will make, for each modelling project or exercise, a document that contains both the R-code *and* the text that provides information on the 'how' and 'why'. 

Typically you will use the Rmarkdown document to include everything that you want to show to us, or to keep for later use. You will still use the console to do some quick calculations, or to produce a quick graph, to see whether you are on the right track. 

## Getting help, examples, demonstrations

R has an extensive help facility. Apart from the Help window launched from the Help menu, it is also available from the command line prompt.

For instance, typing

```{R, eval = FALSE }
 ?log
 ?sin
 ?sqrt
 ?round
 ?Special

will explain about logarithms and exponential functions, trigonometric functions, and other functions.

```{R, eval = FALSE } ?Arithmetic

lists the arithmetic operations in R.

```{R, eval = FALSE }
 help.search("steady")

will list occurrences of the word "steady" in R-commands.

Sometimes the best help is provided by the very active mailing list. If you have a specific problem, just type R: your problem on your search engine. Chances are that someone encountered the problem and it was already solved.

Most of the help files also include examples. You can run all of them by using R-statement example.

For instance, typing into the console window:

```{R, eval = FALSE } example(pairs)

will run all the examples from the pairs help file.
(Try this, ``pairs`` is a very powerful way of visualizing pair-wise relationships.)

Alternatively, you may select one example from the help file, copy it to the clipboard (ctrl-C for  windows users) and then paste it (ctrl-V) in the console window. 
In addition, the R main software and many R-packages come with demonstration material. Typing 

```{R, eval = FALSE }
 demo()

will give a list of available demonstrations in the main software.

```{R, eval = FALSE } demo(transport1D)

will demonstrate some modelling output generated with package  *ReacTran*.

## Small things to remember

### Slash vs. backslash

Pathnames in R are written with forward slashes (``/``). Note that this contrasts to the convention used in Windows, which uses a backslash (``\``) to separate folders in the paths.

To set a working directory in R (e.g., to ``C:\Dirk\MyRcodes\``), use this command:

```{R, eval = FALSE }
setwd("C:/Dirk/MyRcodes/")

Within Rstudio, you can do this via menu Session/Set Working Directory.

Finished vs. unfinished R-commands

If an R-command on one line is syntactically correct, R will execute it, even if the intention was that the command proceeds on the next line.

For instance, if we write

```{R, eval = FALSE } A <- 3 + cos(pi) - sqrt(5)

then A will be assigned the value of $3+cos(\pi)$ and R will display the value of $-\sqrt 5$ ($-2.236068$) in the console.

In contrast, in the following lines,

```{R, eval = FALSE }
A <- 3 + cos(pi) -
     sqrt(5)
A

R will assign to A the value of $3+cos(\pi)-\sqrt 5$, as the command on the first line was not syntactically finished. In this case, R has (correctly) assumed that the command continued on the next line.

Be careful if you want to split a complex statement over several lines! These errors are very difficult to trace, so it is best to avoid this practice.

Exercises --- Using R as a calculator

It is very convenient to use R as a powerful calculator. This can best be done from within the R-console.

Use the console to calculate the value of:

Hint: you may need to look up some help for some of these functions. Typing ```{R, eval = FALSE } ?"+" ?log

will display help for the common arithmetic operators and about the built-in R-function ``log``.


# R-variables

R calculates as easily with *vectors*, *matrices* and *arrays* as with *single numbers*.

R also includes more complex structures such as *data frames* and *lists*, which allow to combine several types of data.

Learning how to create these variables, how to address them and modify them is essential if you want the make good use of the R software.

## Numbers, vectors, matrices and arrays

### Value assignment

When variables are used, they need to be initialised with numbers. Here is an example.

```r
A <- 1
B <- 2
A + B

R can take as arguments for its functions single numbers, vectors, matrices, or arrays.

V <- exp(-0.1)

calculates the exponential of $-0.1$ ($e^{-1}$). The operator <- assigns the result of this calculation to variable V.

V can then be used in subsequent calculations:

log(V)

Note that the assignment of a value to V does not display it.

To display V we simply write in the console:

V

Alternatively, we may assign the result of calculations to a variable AND view the results, by embracing the statement between parentheses:

( X <- sin(3/2*pi) )

Apart from integers, real and complex numbers, R also recognizes infinity (Inf) and Not a Number (NaN). Try:

```{R, eval = FALSE } 1/0 # yields Inf 0/0 # yields NaN 2.3e-8 * 1000 # it is recommended to use 2.3e-8 rather than 2.3*10^-8

Note that ``e-8`` is a preferred syntax for $10^{-8}$, as the latter one is prone to confusion or unexpected errors.

### Vectors

Vectors can be created in many ways; most often we will use:

*   The function ``c()`` combines numbers into a vector.\footnote{This is perhaps THE most important function in R.}
*   The operator ``:`` creates a sequence of values, each by 1 larger (or smaller) than the previous one.
* A more general sequence can be generated by R-function ``seq``.

For instance, the commands

```{R, eval = FALSE }
c(0, pi/2, pi, 3*pi/2, 2*pi)
seq(from = 0, to = 2*pi, by = pi/2)
seq(0, 2*pi, pi/2)
seq(to = 2*pi, from = 0, by = pi/2)

will all create a vector, consisting of: 0, $\pi/2, \pi, \ldots, 2\pi$.

Note that R-function seq takes as input (amongst others) parameters from, to and by (second line). If the order of these arguments is kept, they do not need to be specified by name (third line). But we recommend that you always use the names of the input argument if you want to specify their value. In this case, you can change the order of the input arguments as you like (fourth line).

The next command calculates the sine of this vector and outputs the result:

sin( seq(0, 2*pi, pi/2) )

The next statements create a sequence of integers between 1 and 10 and take the square root of all of them, displaying the result to the screen. The operator <- assigns the sequence to V.

V <- 1:10
sqrt(V)

A peculiar feature of R is that the elements of a vector can also be given names:

( Ocean <- c(total.mass = 1.35e25, volume = 1.34e18, mean.depth = 3690) )
names(Ocean)

This will be a very useful feature in our models.

Matrices

Matrices can also be created in several ways; most often we will use the R-function matrix:

The statement:

A <- matrix(nrow = 3, data = c(1, 2, 3, 4, 6, 8, 10, 12, 14))

creates a matrix A, with three rows, and, as there are nine elements, three columns. Note that the data are input as a vector (using the c() function), and these values are sorted into the matrix column after column (by default). The following commands display the matrix and the square root of its elements.

A
sqrt(A)

Selecting and extracting elements

To select subsets of vectors or matrices, we can either

Simple indexing

The elements of vectors, matrices and arrays are indexed using the [] operator:

To show only the volume of our vector Ocean:

Ocean["volume"]

The following statement takes the elements on the $1^{st}$ and $3^{rd}$ row and on the first two columns of matrix A.

A[c(1, 3), 1:2]

If an index is omitted, then all the rows ($1^{st}$ index omitted) or columns ($2^{nd}$ index omitted) are selected. In the following, the elements on the first three rows of A are multiplied with $2$ and reassigned to the same elements of the matrix.

A[1:3, ] <- A[1:3, ] * 2
A

Similar selection methods apply to vectors:

V[1 : 10]
V[seq(from = 1, to = 5, by = 2)]

Logical expressions

Logical expressions are often used to select elements from vectors and matrices that obey certain criteria. R distinguishes logical variables TRUE and FALSE, represented by the integers $1$ and $0$.\footnote{Note that variables \textit{T} and \textit{F} are reserved in \textit{R} to represent the logical \textit{TRUE} and \textit{FALSE}, respectively. Therefore, if you want to avoid unexpected consequences, you should \textit{not} use \textit{T} and \textit{F} as variables to which you assign values. For example, it might be tempting to use \textit{T} to denote temperature, or \textit{F} to denote flux. It is recommended that you do \textit{not} do this, but use instead more intuitive names such as \textit{Temp} and \textit{Flux}.}

The following will return TRUE for values of sequence V that are greater than 1:

(V <- seq(from = -2, to = 2, by = 0.5))
V > 1

while

V [V > 1]

will select the values from V that are greater than 1, and 

V [V > 1] <- 0

will set to zero all elements in V that are greater than 1.

Removing elements

When the index is preceded by a -, the element is removed. For example,

 A[ ,-1]

will show the contents of matrix A, except the first column, while the command

V[- V < 0]

will only show the positive or zero elements of V (i.e., elements that are not negative).

More complex data structures

Frequently used data structures that are more complex than vectors or matrices are data.frames and lists.

Data.frames

A data.frame superficially looks like a matrix, but its columns may contain different types of elements, e.g., one column may contain strings, another integers, etc.

For instance, the data set iris is of class data.frame:

head(iris)

The data set contains strings and numbers. Each column has a name, and can be accessed by its name using the dollar-sign ($):

mean(iris$Petal.Width)

Lists

A list is a combination of several objects; each object can be of different length. For instance:

LL <- list(Vector = V, Matrix = A)

will combine the previously defined vector V and matrix A in a list called LL. 

names(LL)
LL

Selecting elements from data frames and lists

Lists and data frames can be accessed by their names (using $), or by the square bracket ([ ]) or double square bracket ([[ ]]) operators.

Note: The object resulting from a selection using single brackets [ ], will be a data.frame respectively a list itself; with double brackets [[ ]], one obtains a vector (data.frames) or a variable data-type (lists).

For instance:

LL$Vector * 2

will multiply all values of element Vector with 2.

mean(LL[[1]])

will calculate the mean of element Vector.

Exercises --- Vectors and sequences

Mean of a vector

Sediment depth profiles

       plot(depth, porosity)

Estuarine morphology

Plotting observed data

     (210, 250, 260, 289, 280, 260, 270, 260).

R-functions

One of the strengths of R is that one can make user-defined functions that add to R's built-in functions.

Function definition

Typically, complex functions are written in R-script files or in an R-markdown document, as you will want to use the function several times. For instance,

options(prompt = " ")

Circlesurface <- function (radius) 
    return(pi*radius^2)

defines a function (called Circlesurface), which takes as an input argument a variable called radius and returns the value $\pi*radius^2$, which is the surface of a circle.

After submitting this function to R, we can use it to calculate the surfaces of circles with a given radius:

Circlesurface(10)
Circlesurface(1:10)

The latter statement will calculate the surface of circles with radii 1, 2, ..., 10.

More complicated functions may return more than one element:

Sphere <- function(radius) {
 V  <- 4/3 *pi*radius^3
 S <- 4 *pi*radius^2
 return( list(volume = V, surface = S) )
}

Here, we recognize

The Earth has an approximate radius of $6371~km$, so its volume ($km^3$) and surface ($km^2$) are:

Sphere(6371)

The next statement will only display the volume of spheres with radius 1, 2, $\ldots$ 5:

Sphere(1:5)$volume

Sometimes it is convenient to provide default values for the input parameters.

For instance, the next function estimates the density of "standard mean ocean water" (in $kg~m^{-3}$) as a function of temperature in $^{\circ}$C, TC, and for salinity=0 and pressure = $1~atm$ (Millero, 1981).

The input parameter TC is by default equal to 20 $^{\circ}$C:

Rho_W <- function(TC = 20) {
  rho <- 999.842594 + 0.06793952 * TC - 0.00909529 * TC^2 +
         0.0001001685 * TC^3 - 1.120083e-06 * TC^4 + 6.536332e-09 * TC^5
  return(rho)
}

Note that, within the function body, we ended the first line with a + in order to make clear that the statement is not finished and continues on the next line. It would have been wrong to put the + on the next line (and very difficult to trace this error).

Calling the function without specifying the temperature (input argument TC) uses the default value ($TC=20$):

Rho_W()
Rho_W(20)
Rho_W(0)
Rho_W(TC=0)

Functions in R-packages

For the modelling class you will use many functions from two packages: marelac and ReacTran.

The marelac package

The R package marelac contains many functions useful for the aquatic sciences.

Try:

```{R, eval = FALSE } ?marelac

to see what it contains.

For instance, its function ``sw_dens`` estimates the density of seawater as a function of salinity (``S``), temperature(``t``) and pressure (``p``), and using three different functions (``method``).

To see how it is used, type:

```{R, eval = FALSE }
?sw_dens

To estimate the density for a salinity ranging from 0 to 10, and a temperature of $15~^\circ$C, write:

library(marelac)
sw_dens(S = 0:10, t = 15)

The ReacTran package

The ReacTran package has been written specifically for solving reaction-transport models. It contains many functions to make this type of modelling simple.

For reaction-transport modelling in sediments, for instance, one not only needs a value for the porosity in the centre of the boxes, but also at the box interface. As this requires quite a complicated book-keeping, the R-package ReacTran contains two functions that do that for you.

For instance, to subdivide a $5~cm$ thick sediment column into 10 layers, we write:

sed <- setup.grid.1D(L = 5, N = 10)

This function returns a list that contains many elements that are needed in reaction-transport models:

sed

It is now simple to define the porosity on this grid. First, we define a function that estimates, for a certain sediment depth, $x$, the corresponding porosity. Then, we use the ReacTran function setup.prop.1D to calculate the porosity on this grid:

porfunc <- function (x) return(0.7 + 0.3*exp(-x))
porosity <- setup.prop.1D(func = porfunc, grid = sed)

Porosity is now defined, both in the middle of slices and at the slice interfaces: 

porosity

Exercises --- R-functions

R-function to estimate saturated oxygen concentrations

The saturated oxygen concentration in water ($\mu mol~kg^{-1}$), also called oxygen solubility, can be calculated based on an empirical formula $SatOx = e^A$, where

$$A = -173.9894 + 25559.07/T + 146.4813 \times log_e(T/100) -22.204 \times T/100 + $$ $$ S \times (-0.037362+0.016504 \times T/100-0.0020564 \times T/100 \times T/100)$$

T is temperature in Kelvin (Tkelvin = Tcelsius + 273.15), and S is salinity (reported as unitless, but meaning a value in g/kg).

Molecular diffusion coefficient

Package marelac contains a function to calculate molecular diffusion coefficients.

Note: the diffcoeff function from the marelac package returns a data.frame or a list. For plotting, it is easiest to subset this so as to have a vector. Thus, 

diffcoeff()$O2

will provide the output in a format that is easy to work with.

R-function sphere

Organisms can have many shapes, from spherical to rod-like to amorphous.

In order to create reaction-transport models that describe for instance the oxygen concentration in the body of these organisms, we need to estimate the surface area at certain distances from the centre of their body.

Assume a spherical ciliate with a diameter of 100 $\mu m$. For modelling purposes, assume that this ciliate consists of concentric spheres, each 1 $\mu m$ thick.

Porosity profile and estuarine morphology as a function

Estuarine morphology using ReacTran

Solving ordinary differential equations in R

Throughout the modelling course, you will specify many models in the form of differential equations. Here, it is shown how differential equations are implemented and solved in R.

Consider the following set of two ordinary differential equations: $$ \frac{{dA}}{{dt}} = r \cdot (x - A) - k \cdot A \cdot B $$ $$\frac{{dB}}{{dt}} = r \cdot (y - B) + k \cdot A \cdot B $$ Here, $A$ and $B$ are called the state variables, $\frac{dA}{dt}$ and $\frac{dB}{dt}$ are the time derivatives (also called the rates of change), while $r$, $x$, $y$ and $k$ are constant parameters.

Specifying the differential equation model

The first step to solving this set of differential equations in R is to define the model function, which specifies the right-hand side of the differential equations.\footnote{Note that in R we use "*" for multiplication, and not $\cdot$ or $\times$ as in the mathematical formalism.}

This function has three different arguments as input: the actual time (t), the values of the state variables (state) and the values of the parameters (parameters). 

model <- function(t, state, parms) {

with( as.list(c(state, parms)), {

 dA.dt <- r*(x-A) - k*A*B
 dB.dt <- r*(y-B) + k*A*B
 return (list(c(dA.dt, dB.dt), sum = A+B))

 })  # end of the *with* function
}

This function simply calculates the time derivatives of the state variables (dA.dt and dB.dt) and an output variable called "sum". The derivatives are combined in a vector, and both this vector and the output variable are returned as a list.

The R-statement with( as.list(c(state, parms)), { ensures that the state variables and parameters can be addressed by their names. This statement embraces all other statements in the function --- it ends at the line that says }).

Solving the differential equation model

Before we can actually solve this model, we need to:

pars      <- c(x = 1, y = 0.1, k = 0.05, r = 0.05)
state.ini <- c(A = 1, B = 1)  # the same order as dA.dt and dB.dt in the model function!!
time.seq  <- seq(from = 0, to = 300, by = 1)

The model can now be solved. To do so, we use the integration routine ode, which can be found in the R-package deSolve (Soetaert et al., 2010). This package is loaded first (using function require).

require(deSolve)

The routine ode will calculate a value for the state variables A and B at each time value specified in the vector time.seq. It does so by numerical integration. The name (ode) hints at the type of differential equations that this function solves, which are Ordinary Differential Equations.

The actual numerical solution of our ODE model is obtained by the following single R-statement:

out <- ode(y = state.ini, times = time.seq, func = model, parms = pars)

The output is stored in a matrix, called out. All we need to do now is to plot this model output. Before we do so, we can have a look at the top few rows of the output matrix:

head(out,4)

The data are arranged in three columns: first the time, then values of the state variables A and B at each time point, followed by the output variable called sum.

We can plot the output in several ways. The easiest is to plot the entire object at once, either plotting each variable in a separate graph,

plot(out, xlab = "time", ylab = "concentration", lwd = 2, type = "l", mfrow=c(1,3))

or all variables in one graph.

matplot.0D(out, lty = 1, lwd = 2, main=NA)

We can run the model with different values of the parameter k, store the output in a matrix out2, and plot the first and second run at the same time. In the code below, we first take a copy of the parameter vector (parms2), and then change the parameter named ("k"); we then solve the model, passing the updated parameter vector (parms = parms2). We can plot the outcome of the two runs at once. Note that by the argument mfrow = c(1,3), we force the output to be in one row and 3 columns (mfrow stands for multiple figures in a row).

pars2      <- pars
pars2["k"] <- 0.5
out2 <- ode(y = state.ini, times = time.seq, func = model, parms = pars2)
plot(out, out2, lty = 1, lwd = 2, mfrow = c(1, 3))

For completeness, the entire R-code that solves and plots the above set of differential equations is listed here.

model <- function(t, state, parms) {
  with (as.list(c(state, parms)), {

    dA.dt <- r*(x-A) - k*A*B
    dB.dt <- r*(y-B) + k*A*B
    return (list(c(dA.dt, dB.dt), sum = A+B))

 })
}

pars      <- c(x = 1, y = 0.1, k = 0.05, r = 0.05)
state.ini <- c(A = 1, B = 1)
time.seq  <- seq(from = 0, to = 300, by = 1)

require(deSolve)
out <- ode(y = state.ini, times = time.seq, func = model, parms = pars)
plot(out, xlab = "time", ylab = "concentration", lwd = 2, type = "l")

Steady-state conditions of differential equations

Sometimes we are interested in the conditions of a differential equation where the state variables do not change anymore. One way to obtain this is to use the function steady from the R-package rootSolve (Soetaert, 2009):

require(rootSolve)
STD  <- steady(y = state.ini, func = model, parms = pars)
STD$y; STD$sum

Exercises --- Solving ordinary differential equations in R

The Lotka-Volterra model

The Lotka-Volterra model is a famous model that describes predator-prey interactions or competitive interactions between two species. A. J. Lotka and V. Volterra formulated the model almost simultaneously in the 1920's.

The Lorenz Butterfly

The Lorenz equations represents the first set of differential equations in which chaotic behaviour was discovered. These three differential equations represent an idealized model for the circulation of air within the atmosphere of the Earth.

$$ \frac{{dX}}{{dt}} = - \frac{8}{3} \cdot X + Y \cdot Z$$ $$\frac{{dY}}{{dt}} = - 10 \cdot (Y - Z) $$ $$\frac{{dZ}}{{dt}} = - X \cdot Y + 28\cdot Y - Z $$

References

R Core Team (2018). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. \url{https://www.R-project.org/}.

Soetaert, Karline and Meysman, Filip, 2012. Reactive transport in aquatic ecosystems: Rapid model prototyping in the open source software R Environmental Modelling & Software, 32, 49--60.

Soetaert K. (2009). rootSolve: Nonlinear root finding, equilibrium and steady-state analysis of ordinary differential equations. R-package version 1.6

Karline Soetaert, Thomas Petzoldt, R. Woodrow Setzer (2010). Solving Differential Equations in R: Package deSolve. Journal of Statistical Software, 33(9), 1--25. \url{http://www.jstatsoft.org/v33/i09/} DOI 10.18637/jss.v033.i09

Karline Soetaert and Thomas Petzoldt (2018). marelac: Tools for Aquatic Sciences. R package version 2.1.7.

Andreas F. Hofmann, Karline Soetaert, Jack J. Middelburg and Filip J. R. Meysman, 2010. AquaEnv - an Aquatic acid-base modelling Environment in R. Aquatic Geochemistry. DOI: 10.1007/s10498-009-9084-1.



dynamic-R/RTM documentation built on Feb. 28, 2025, 1:23 p.m.

R Package Documentation

rdrr.io home R language documentation Run R code online

Browse R Packages

CRAN packages Bioconductor packages R-Forge packages GitHub packages

We want your feedback!

Note that we can't provide technical support on individual packages. You should contact the package authors for that.
GitHub issue tracker
ian@mutexlabs.com
Personal blog

 
Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.

Close