Stratification stability of tropical lakes and their responses to
climatic changes: Lake Towuti (Indonesia)
A THESIS
SUBMITTED TO THE FACULTY OF THE GRADUATE SCHOOL
OF THE UNIVERSITY OF MINNESOTA
BY
Tongyao Pu
IN PARTIAL FULFILLMENT OF THE REQUIREMENTS
FOR THE DEGREE OF
MASTER OF SCIENCE
Sergei Katsev
Jan, 2023
© Tongyao Pu 2023
ALL RIGHTS RESERVED
Acknowledgements
There are many people that have earned my gratitude for their contribution to my time in
graduate school.
Dr. Sergei Katsev, being my academic advisor, has provided me with an enormous amount of
help during my master’s. He guided me to see the big research picture; asked me important
questions; answered my questions patiently; and helped review and revise my thesis multiple
times. Sergei is also very supportive for my career development. He advised me in attending
conferences and helped me advertise myself to recruiters.
My committee members, Dr. Sam Kelly, Dr. Chan Lan Chun, and Dr. Byron Steinman came
up with very interesting questions. Their suggestions are inspiring but practical, which greatly
improved the current thesis content compared to the previous version.
Dr. James Russell, Dr. Doug Haffner, and Dr. Sean Crowe have provided me with long term,
high quality data they have collected from Lake Towuti. Without their help, the simulations of
Lake Towuti will be more foggy.
My master’s experience cannot be separated from the amazing courses I had. Physical
Limnology taught by Dr. Jay Austin, Fluid Dynamics taught by Dr. Sam Kelly, Data Analysis
taught by Dr. Rik Gran, and Environmental Modeling taught by Dr. Nate Johnson all helped
me build up my knowledge and skills to support my research work.
i
I also met countless wonderful people who made Duluth feel like a second home to me,
especially my friends from the WRS and the physics department, and my lovely students. I
am also grateful to have the chance to cat-sit Chandler, Sherbert, and Neko. These wonderful
companions greatly reduced my sense of loneliness during long and dreadful holiday seasons.
Unfortunately, Sherbert passed away shortly after my defense. This brought me great sadness.
I will forever remember him.
Last but not least, my parents, Jianying Xia and Zhifeng Pu, are very supportive of my
education, both financially and mentally. My grandma, Jinhua Li, who raised me growing
up, is also very very supportive of my decisions. Unfortunately, we have not been able to get
together for 3.5 years due to COVID. I hope we will see each other again soon!
ii
Abstract
Tropical lakes have different physical dynamics than temperate lakes. They are less likely to
experience thermal convection and can be permanently stratified by temperature gradients alone
(thermogenic meromixis). Climate research suggests that tropical lakes are most susceptible
to increasing surface temperature and stratification with future prediction of increasing air
temperature. Tropical meromictic lakes are also important sites for paleoclimate research, as
the anoxic and stagnant bottom layer preserves laminated sediments. Ancient Lake Towuti,
located in Sulawesi island in Indonesia, is believed to be meromictic sustained without profound
salinity gradient. This project aims to investigate the strength and dynamics of thermogenic
meromixis in tropical lakes. To address this problem, it simulates the physical dynamics of the
lake via the use of a 1-D k-ε hydrodynamic model Simstrat. The simulations reveal that Lake
Towuti may be oligomictic rather than meromictic, as was previously assumed. More generally,
tropical lakes much deeper than 80 m are likely to be oligo- or meromictic, whereas shallower
lakes are likely to be warm monomictic. Air temperature stands out to be the most influential
meteorological variable affecting stratification. Oligomixis may be regulated by a negative
feedback that plays out over multiple years: A full water column mixing event will decrease
the bottom water temperature of a lake and eventually increase surface-bottom temperature
difference and average stability. This leads to increasing difficulties for whole-water column
mixing. The temporal variation in the Lake Towuti stratification is only weakly correlated with
the ENSO index. During the Last Glacial Maximum (LGM), Lake Towuti likely experienced
colder and drier climate, which lowered lake level, and induced more frequent whole-water
column mixing. A 1-D biogeochemical reaction-transport model, which was provided with the
output of hydrodynamic simulations for the intensity of turbulent mixing, further confirmed a
possibility of the whole water column mixing in Lake Towuti as recently as 2008. A weak
iii
mixing event lasting a few weeks during the dry season would be consistent with the observed
chemical profiles of O2, H2S, and Fe2+. Coexisting dissolved iron and hydrogen sulfide at low
concentrations (few µmolL−1) likely indicate observations reflecting a transient state. Findings
suggest that, under the future climate characterized by warmer air temperatures, Lake Towuti is
likely to be more stable. Stronger stability and isolation of monimolimnion water would lead
to further accumulation of dissolved iron. Similar tropical lakes whose stability are controlled
primarily by water temperature are expected to have increasing stability in the future as well.
iv
Contents
Acknowledgements i
Abstract iii
List of Tables viii
List of Figures ix
1 Introduction 1
1.1 Physical Stratification of Lakes . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Chemical Stratification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.3 Hydrodynamic Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.4 Reaction-Transport Biogeochemical Modeling . . . . . . . . . . . . . . . . . . 9
1.5 Tropical Ferruginous Lake Towuti . . . . . . . . . . . . . . . . . . . . . . . . 9
1.6 Research Questions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2 Physical Model 15
2.1 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.1.1 Modeling Framework . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.1.2 Hydrodynamic Model . . . . . . . . . . . . . . . . . . . . . . . . . . 17
v
2.1.3 Stratification Strength Indicators . . . . . . . . . . . . . . . . . . . . . 26
2.1.4 Sensitivity Test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
2.1.5 Correlation with the ENSO Climate Index . . . . . . . . . . . . . . . . 29
2.1.6 Lake Towuti in the Last Glacial Maximum . . . . . . . . . . . . . . . 29
2.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
2.2.1 Stratification of Lake Towuti . . . . . . . . . . . . . . . . . . . . . . . 30
2.2.2 Sensitivity Test Results . . . . . . . . . . . . . . . . . . . . . . . . . . 39
2.2.3 Lake Towuti during the Last Glacial Maximum . . . . . . . . . . . . . 42
2.3 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
2.3.1 Multi-year Stability and Potential Oligomixis of Lake Towuti . . . . . 44
2.3.2 Mechanisms of Stratification and Mixing . . . . . . . . . . . . . . . . 45
2.3.3 Lake Towuti During the Last Glacial Maximum . . . . . . . . . . . . . 52
2.3.4 The Impact of Global Climate Changes . . . . . . . . . . . . . . . . . 52
2.3.5 Correlation with ENSO . . . . . . . . . . . . . . . . . . . . . . . . . . 53
2.3.6 What Makes a Tropical Lake Meromictic . . . . . . . . . . . . . . . . 54
2.3.7 Performance of the Model, Caveats of Interpretations . . . . . . . . . . 56
2.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
3 Biogeochemical Model 58
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
3.2 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
3.2.1 Chemistry Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
3.2.2 Using the Kz from the Hydrodynamic Model . . . . . . . . . . . . . . 60
3.2.3 Simulation Procedure . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
3.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
3.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
vi
3.4.1 Gentle Mixing Reproduces Observed Reduced Water Column . . . . . 69
3.4.2 Limitations of the Chemical Model . . . . . . . . . . . . . . . . . . . 70
3.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
4 Conclusion 72
5 Bibliography 75
Appendix A. Glossary and Acronyms 84
A.1 Acronyms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
Appendix B. Hydrodynamic Modeling Parameters 86
Appendix C. Chemsitry Model 89
C.1 Model Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
C.2 Biogeochemical Model Reactions . . . . . . . . . . . . . . . . . . . . . . . . 94
C.3 Biogeochemical Model Reaction Rates . . . . . . . . . . . . . . . . . . . . . . 95
C.4 Biogeochemical Model Parameters . . . . . . . . . . . . . . . . . . . . . . . . 96
Appendix D. Extra Figures 98
vii
List of Tables
1.1 Selected large tropical lakes (between 23°26’ S and 23°26’ N) and their
stratification status. Continental location is provided : AF - Africa, AS -
Asia, SA - South America. Not all the currently known meromictic lakes are
included in this table. Some small meromictic lakes can be found in Alcocer
(2017) (Mexican lakes) and De Crop and Verschuren (2021) (Uganda lakes) . . 4
2.1 Variations in meteorological variables used in the Monte Carlo analysis. . . . . 28
2.2 Simulation runs for investigating the effect of LGM conditions. . . . . . . . . . 29
2.3 Correlation between ERA5 climate variables and the SSI of Simstrat simulation
run with them as inputs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
2.4 Concluded mixing criteria for tropical lakes without salinity gradient . . . . . . 57
3.1 Parameters defining analytical Kz profiles . . . . . . . . . . . . . . . . . . . . 62
A.1 Acronyms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
B.1 Physical parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
C.1 Chemical model boundary conditions . . . . . . . . . . . . . . . . . . . . . . 91
C.2 Chemical model reactions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
C.3 Chemical reaction rates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
C.4 Chemical model parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
viii
List of Figures
1.1 The bathymetry of Lake Towuti (modified from Russell et al. (2016)). . . . . . 10
1.2 Lake Towuti catchment (including the upstream Lake Matano). . . . . . . . . . 11
1.3 Towuti lake level via satellite altimetry Sulistioadi et al. (2015) . . . . . . . . . 12
2.1 The modeling process. Rounded rectangles represent process / models;
Document symbols are used to indicate inputs and outputs . . . . . . . . . . . 17
2.2 ERA5 reanalysis climate variables from 1989 to 2016, marked with ENSO
events indicated by the MEI.v2 index: El Nin˜o (MEI > 0.5) are marked with
red bands; La Nin˜a (MEI < -0.5) are marked with blue bands. . . . . . . . . . 24
2.3 Comparison of the annual climate data around Lake Towuti from different
sources. ERA5 reanalysis data are presented as mean monthly averages derived
for 2000-2016. Observational data are from weather stations in Kendari and
Palu city (https://www.weather-atlas.com/, year range unknown). . . . . . . . . 25
2.4 Temperature distributions and evolution in Lake Towuti simulated using the
unmodified ERA5 reanalysis meteorological data as model input. . . . . . . . . 30
2.5 Comparison of simulations with the temperature data from a moored instrument
array (Dr. James Russell, pers. comm.; observational data after June 2014 have
less confidence due to instrumentation problems). . . . . . . . . . . . . . . . . 31
ix
2.6 Residuals from the model (model minus mooring temperature). Warm color
indicates warmer simulated temperature, colder color indicates colder simulated
temperature compared to the observation . . . . . . . . . . . . . . . . . . . . . 32
2.7 Comparison between observational CTD temperature profiles and modeled
profiles at the corresponding time. Vertical dashed lines represent 28.5 °C. It is
plotted to provide an easy comparison between different CTD profiles. . . . . . 33
2.8 The simulated time series for the overall stratification indicators: temperature
difference indicator (TDI), maximum buoyancy frequency (MBF), heat content
indicator (HCI), and Schmidt Stability index (SSI). . . . . . . . . . . . . . . . 34
2.9 Heat fluxes across the lake surface calculated by Simstrat in 2012. From top to
bottom: shortwave radiation, longwave radiation reflected from sky into the lake
(longwave in), longwave radiation from lake into the atmosphere (longwave
out), sensible heat, latent heat, and the net heat flux. Red shaded regions mark
a typical Sulawesi wet season, and blue shades indicate a typical dry season. . . 35
2.10 Scatter plots of SSI against model input climate variables (from ERA5
reanalysis). Because the (monthly) data points correspond to monthly averages,
a brief mixing of the water column corresponds to low values of SSI but not
necessarily zero. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
2.11 The MEI.v2 index of ENSO (> 0.5 - El Nin˜o, <−0.5 - La Nin˜a) and the moving
(12 months) mean of the Schmidt Stability (SSI) from the simulation. . . . . . 38
x
2.12 Normalized cross correlation between the MEI.v2 index and unfiltered Schmidt
Stability index (SSI) (all detrended). This figure shows that the high MEI.v2
indexes, i.e. El Nin˜o events, are correlated with weaker stability simulated
in the model, and vice versa for La Nin˜a events. The maximum positive
correlation occurs with a -4 month lag with the correlation coefficient of 0.197.
The maximum negative correlation occurs with a 14 months lag with the
correlation coefficient of 0.134. These are not strong correlations. If ENSO
climate variation is a precursor for the lake status change, the effect occurs
approximately 4 months later. . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
2.13 Correlation between the MBF indicator for the maximum strength of the
thermocline (dots), the WTD indicator for the minimum temperature difference
(circles), and the mean values of climate variables, from 100 runs of Monte
Carlo simulations (2001 - 2008 period). The wind speed and air temperature
have the most significant impact on the stability of stratification. . . . . . . . . 40
2.14 Correlation between MBF (dots), WTD (circles), and the mean values of
climate variables, from 100 runs of Monte Carlo simulations (2009 - 2016
period). Notice that a lot of simulations’ WTD between 2009 - 2016 period is
non-zero. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
2.15 WTD and MBF indicator comparison between pre 2008 and post 2009 from
simulation. With criteria of WTD < 0.4 °C indicating mixing, there are 68
simulations simulated mixing before 2008, but only 3 of them also simulated
mixing between 2009 - 2016 . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
xi
2.16 Simstrat simulation exploring conditions during the Last Glacial Maximum
(LGM) for different water levels (100 m, 150 m, and 190 m) and for different
climate conditions (left column : current ERA5 reanalysis climate conditions,
right column : LGM climate conditions described in subsection 2.1.6). The
simulated processes are similar to the nominal hydrodynamic simulation for the
modern lake, but the outputs are cropped so that it represents a relatively steady
stage of establishing LGM hydrodynamics. . . . . . . . . . . . . . . . . . . . 43
2.17 Schmidt Stability and Heat Content (relative to the water temperature at the
calculated depth). The magnitude of SS indicates the mechanical energy
required to homogenize the water column to the desired depth while the
magnitude of HC indicates the negative thermal energy required to cool down
all the water above that depth. . . . . . . . . . . . . . . . . . . . . . . . . . . 47
2.18 Conceptual diagram of the mixing-stratification dynamics in Lake Towuti.
Starting from a weak thermal gradient, when there is a strong negative surface
heat flux, surface water cools down to a point where the whole water column
mixes. If the negative heat flux is strong, the whole lake cools down by
a significant amount. After this event, the lake gradually reestablishes its
surface temperature, while the bottom water stays cold. Lake stability after
restratification gets stronger, and thus the lake is less vulnerable to mixing.
With enough time, the gradual diffusional process will blend out this strong
stratification, and make it weak again, to its original state. . . . . . . . . . . . . 51
3.1 Turbulent diffusivity of temperature, equivalent to Kz, simulated by the
hydrodynamic model discussed in Chapter 2 . . . . . . . . . . . . . . . . . . 61
3.2 Typical stratification (left) and mixing (right) Kz profiles extracted from the
hydrodynamic model (square and circular markers) and their corresponding
analytical approximation (solid line) used in the biogeochemical model. . . . . 63
xii
3.3 Development of chemical stratification starting from fully oxygenated initial
condition. A typical ”stratification” Kz was applied over the course of the entire
simulation. Red circles are observational profiles (Dr. Sean Crowe, pers. comm.) 66
3.4 Chemical profiles during a mixing event. Starting from the initial conditions
generated from the 40 years of stratification (Fig. 3.3), a typical mixing Kz was
applied over the course of the simulation (2 months). . . . . . . . . . . . . . . 67
3.5 Chemical profiles during the restratification process. Starting from conditions
after mixing (Fig. 3.4), a typical stratify Kz was applied over the course of the
simulation (seven years) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
C.1 Illustration of the important boundary conditions and reactions of the
biogeochemistry model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
D.1 Root mean squared error (RMSE) between model and mooring (2012 - 2015)
temperature data as a function of time. The dashed line indicates an average
RMSE of the dataset. The solid vertical line represents June 2014. After this
time, the confidence in mooring data decreased. This figure illustrates that,
RMSE significantly increased during 2012-2013 wet season. However, this
increasing of error did not accumulate in the proceeding simulation. . . . . . . 98
D.2 RMSE between model and mooring (2012 - 2015) temperature data as a
function of depth. Dashed line indicates the overall RMSE of the dataset.
RMSE is larger than average at the surface, and is minimal at around 100 m.
This indicates that Simstrat simulation errors are mostly observed at the surface
rather than deeper part of the lake and that the modeling of surface heat flux
over Lake Towuti results in discrepancy between simulation and observation.
Bottom temperature simulated by simstrat is more satisfactory. . . . . . . . . . 99
xiii
D.3 RMSE sensitivity with meteorological variables. This figure is generated with
the same data in 2.13, except the y-axis is now plotting the RMSE between the
model and mooring data. This figure illustrates that the current configuration of
meteorological variables is a satisfactory representation that allows simulation
agreement with observation. No better combinations of meteorology can
significantly improve the performance of the model. . . . . . . . . . . . . . . . 100
D.4 Cross correlation between meteorological variables and MEI.v2 index.
negative lag indicates MEI.v2 is prior to the local meteorological variable.
The correlations are mostly weak (with Pearson’s r smaller than 0.5 for all
the climate variables). MEI.v2 index represents the air temperature over Lake
Towuti several months later, but has delays in representing the wind speed,
solar radiation, vapor pressure, and cloud cover over Lake Towuti. . . . . . . . 101
xiv
Chapter 1
Introduction
1.1 Physical Stratification of Lakes
Physical stratification is one of the important characteristics of a lake. The density of water is
affected by temperature and the concentrations of solutes, the distributions of which within the
lake are influenced by mixing rates. A strong density gradient within the water column results
in stratification, creating horizontal layers with different physical and chemical properties.
Density gradients restrict vertical transport of chemicals and organisms in the water column,
modify chemical conditions, affect reaction rates, and influence the way minerals and
biomarkers become preserved in the sediment. Changes in stratification patterns strongly affect
the functioning of lake ecosystems by controlling the transport rates of nutrients and oxygen,
and the distributions of plankton (Woolway et al., 2021). Change of stratification can even lead
to catastrophic events. For example, shoaling of deep water deficient in dissolved oxygen can
lead to fish kills (Fukushima et al., 2017). In rare cases, release of accumulated carbon dioxide
as a result of destabilization of stratification can asphyxiate humans and livestock, such as in
the 1984 Lake Monoun (Sigurdsson et al., 1987) and 1986 Lake Nyos (Kling et al., 1989)
1
2disaster.
Most water bodies stratify vertically for some periods of time, followed by periods of mixing.
Based on the frequency of mixing and the type of stratification, lakes are classified into amictic,
cold monomictic, dimictic, warm monomictic, oligomictic, and polymictic (Hutchinson
and Lo¨ffler, 1956). Amictic lakes are later believed to be a category of meromictic lakes
(Walker and Likens, 1975). Permanent stratification (meromixis) in lakes is a special case of
stratification (’mero-’ stands for ’partial’). In meromictic lakes, the upper part (mixolimnion)
of a lake mixes and stratifies based on the local climate conditions, while the bottom part
(monimolimnion) stays undisturbed. Prolonged stability limits the fluxes of energy and matter
into the monimolimnion (Boehrer et al., 2017). As a result, chemical composition becomes
very different between the surface and the bottom parts of the lake.
Walker and Likens (1975) classified the mechanisms by which meromixis is created into
ectogenesis and endogenesis. Ectogenesis is established by an external inflow of freshwater
at the surface or an inflow of denser water at depth. Endogenesis is generated by a lake’s
internal processes of accumulating salinity. In most meromictic lakes, whether ectogenic or
endogenic, stratification is sustained by a salinity gradient that prevents water mixing between
the mixolimnion and monimolimnion. However, some tropical lakes can be meromictic
without a strong salinity gradient. Katsev et al. (2010) termed the phenomenon ’thermogenic
meromixis’. Deep tropical lakes are known to be stabilized by a thermal gradient alone.
Examples include Lake Matano. Unlike their temperate counterparts, they do not experience
thermal convection that becomes more likely as the surface water approaches the temperature
of maximum density of 4 °C (Katsev et al., 2017).
Table 1.1 lists tropical meromictic lakes, with their surface area A0, maximum depth zm, and
relative depth zr (zr =
50∗ zm ∗
√
pi√
A0
, Walker and Likens (1975)). The relative depth metric was
developed for mid-latitude lakes as a simple predictor for a lake’s mixing status. Lakes with
3zr > 4% were suggested to be likely seasonally or permanently stratified. This categorization,
however, does not work for tropical lakes. Only three out of nine tropical meromictic lakes
mentioned in Table 1.1 have zr greater than 4 %. This table also lists some mechanisms for
establishing permanent stratification, based on available information. It could be seen that
tropical meromictic lakes are dynamic with drastically different shape, size, and properties.
Many large deep tropical lakes are located in tectonic depressions that existed for hundreds
of thousands of years. Such lakes accumulate thick piles of sediments that contain valuable
records of past climates. Meromixis and anoxia in such lakes preserve laminated sediments
(Anderson, 1985), which preserve detailed chronology. Thus, such lakes often become sites of
paleoclimate research. For example, the International Continental Scientific Drilling Program
(ICDP, https://www.icdp-online.org/projects/) has funded projects on the East African lakes
Malawi and Tanganyika, and the Indonesian lake, Towuti.
The stratification status of a lake is often linked to its climate, which determines the heat
and momentum fluxes received by the lake’s surface. The surface heat fluxes are commonly
subdivided to shortwave radiation, longwave radiation, sensible heat transfer (thermal
conduction), and latent heat transfer (evaporation). These heat flux components are linked
to climate variables, such as air temperature, wind speed, solar radiation, vapor pressure,
and cloud cover (Imboden and Wu¨est, 1995). Increases in air temperature, vapor pressure,
and solar radiation would lead to heat fluxes into the lake, resulting in warmer lake surface
and potentially stronger stratification. Increases in wind speed lead to heat loss via negative
sensible and latent heat fluxes. The wind also exerts mechanical stress at the surface of the
lake, transferring momentum, which increases turbulent kinetic energy and leads to mixing
(Imboden and Wu¨est, 1995). Thus, increasing wind speed may destabilize stratification. Cloud
cover has multiple effects that partially cancel each other. Although clouds block sunlight and
inhibit shortwave radiation, they also reflect more of the longwave radiation back to the Earth
4Ta
bl
e
1.
1:
Se
le
ct
ed
la
rg
e
tr
op
ic
al
la
ke
s
(b
et
w
ee
n
23
°2
6’
S
an
d
23
°2
6’
N
)
an
d
th
ei
r
st
ra
tifi
ca
tio
n
st
at
us
.
C
on
tin
en
ta
ll
oc
at
io
n
is
pr
ov
id
ed
:
A
F
-
A
fr
ic
a,
A
S
-
A
si
a,
SA
-
So
ut
h
A
m
er
ic
a.
N
ot
al
lt
he
cu
rr
en
tly
kn
ow
n
m
er
om
ic
tic
la
ke
s
ar
e
in
cl
ud
ed
in
th
is
ta
bl
e.
So
m
e
sm
al
lm
er
om
ic
tic
la
ke
s
ca
n
be
fo
un
d
in
A
lc
oc
er
(2
01
7)
(M
ex
ic
an
la
ke
s)
an
d
D
e
C
ro
p
an
d
V
er
sc
hu
re
n
(2
02
1)
(U
ga
nd
a
la
ke
s)
L
ak
e
N
am
e
L
oc
at
io
n
A
0
[k
m
2 ]
z m
[m
]
z r
[%
]
m
er
om
ix
is
m
ec
ha
ni
sm
re
fe
re
nc
e
1
L
ak
e
N
yo
s
A
F
1.
58
20
8
14
.6
6
di
ss
ol
ve
d
ga
s
K
lin
g
et
al
.(
19
89
)
2
L
ak
e
M
on
ou
n
A
F
0.
31
99
15
.7
6
di
ss
ol
ve
d
ga
s
Si
gu
rd
ss
on
et
al
.(
19
87
)
3
L
ak
e
K
iv
u
A
F
27
00
48
0
0.
82
di
ss
ol
ve
d
ga
s
So
m
m
er
et
al
.(
20
19
)
4
L
ak
e
So
na
ch
i
A
F
0.
18
7
1.
46
sa
lin
ity
gr
ad
ie
nt
M
ac
In
ty
re
an
d
M
el
ac
k
(1
98
2)
5
L
ak
e
Ta
ng
an
yi
ka
A
F
32
90
0
14
70
0.
72
th
er
m
og
en
ic
C
hi
ta
m
w
eb
w
a
(1
99
9)
6
L
ak
e
M
al
aw
i
A
F
29
60
0
70
6
0.
36
te
m
p
an
d
sa
lg
ra
di
en
t
W
ue
st
et
al
.(
19
96
)
7
L
ak
e
M
at
an
o
A
S
16
4
59
0
4.
08
th
er
m
og
en
ic
K
at
se
v
et
al
.(
20
10
)
8
L
ak
e
To
ba
A
S
11
30
50
5
0.
01
3
th
er
m
og
en
ic
Fu
ku
sh
im
a
et
al
.(
20
22
)
9
L
ak
e
A
tit
la
n
SA
13
7
34
1
2.
58
R
ej
m
a´n
ko
va´
et
al
.(
20
11
)
10
L
ak
e
L
an
ao
A
S
35
7
12
2
0.
00
06
w
ar
m
m
on
om
ic
tic
Fr
ey
(1
96
9)
5surface.
The stratification status of a lake can change over time. The change can happen episodically
or regularly. An increasing amount of research has recently focused on the effects of climate
change on lake stratification. Overall, studies suggest warming air temperatures and reduction
in wind speed increase lake stability in lakes around the world (Verburg and Hecky (2009),
Lake Tanganyika, Africa, meromictic; Weinberger and Vetter (2012), Ammersee, Germany,
dimictic; Sahoo et al. (2016), Lake Tahoe, US, warm-monomictic; Woolway et al. (2017),
Vo˜rtsja¨rv, Estonia, polymictic; Oleksy and Richardson (2021), Mohonk Lake, US, dimictic;
Stetler et al. (2021), Crater Lake, US, meromictic; Fukushima et al. (2017), Indonesian
deep tropical lakes; Woolway and Merchant (2019), lakes globally, including tropical lakes).
The general patterns are that lake temperature will be increasing in the future, leading to
strengthening of stratification; Stratification will be shifting from dimictic to monomictic to
meromictic (Woolway and Merchant, 2019).
Although climate warming generally increases lake stratification, individual lakes may respond
differently. Hambright et al. (1994) (Lake Kinneret) and Kraemer et al. (2015) (lakes around the
world) summarized the importance of lake morphology and water temperature in defining the
lake’s response to climate change. Deeper lakes become more strongly stratified, with deeper
thermoclines. The lake surface area serves a less important role. Warmer lakes will experience
greater change in enhancing stratification and decreasing depth of thermocline. Richardson
et al. (2017) also added the effect of water clarity in lakes. More transparent lakes tend to
have greater increases in surface water temperature and stability in response to warming air
temperature. Kraemer et al. (2015) points out that tropical lakes tend to be more influenced by
climate change and experience stronger stratification compared to mid-latitude lakes.
61.2 Chemical Stratification
Prolonged physical stratification limits the exchange of substances between surface and bottom
water. The transport of oxygen, in particular, becomes restricted. Over time, oxygen in the deep
waters becomes depleted by the respiration reactions that oxidize organic matter in the water
column below the density gradient.
Once oxygen is depleted, other chemical species act as electron acceptors in a sequence of
redox reactions that fuel anaerobic microbial metabolisms. The sequence of organic matter
mineralization reactions includes aerobic respiration, nitrate reduction, manganese reduction,
iron reduction, sulfate reduction, and methanogenesis. While the sequence is generally dictated
by the thermodynamic energy gain of the reactions (Froelich et al., 1979, Berner, 1981),
depending on the reaction kinetics and activities of chemicals and the order of the reactions, the
sequence may change. In particular, the overlap between iron reduction and sulfate reduction
is commonly observed (Canfield and Thamdrup, 2009).
In anoxic environments, accumulation of dissolved ferrous iron (Fe2+, a product of iron
reduction) and dissolved hydrogen sulfide (H2S, a product of sulfate reduction) are of great
importance to lake geochemistry. The prevalence of either of these two species define
ferruginous (Fe2+-rich) or euxinic (H2S-rich) conditions. Aqueous Fe2+ and H2S co-precipitate
as poorly crystalline ferrous sulfide (FeS), which matures into more stable pyrite (FeS2). Thus,
only one of the aqueous species can be present at substantial (micromolar) concentrations,
as the less abundant species would precipitate out of the solution. The classification of
geochemical conditions as ferruginous or euxinic varies under different research contexts. To
interpret the past water column environment based on sediment record, geologists would look
at the ratio of pyrite (FeS2) to highly reactive iron (FeHR). Ratios above about 0.7 are taken
as evidence of euxinic conditions, whereas those below this threshold indicate ferruginous
conditions (Poulton and Canfield, 2011).
7In this study, ferruginous conditions are defined as the concentration of Fe2+ being significantly
higher than the concentration of H2S.
Understanding what controls the relative abundances of ferrous iron and hydrogen sulfide is
important in both ecological and geological contexts. In modern settings, sulfide can limit
the amounts of trace metals being available for organisms. A large quantity of sulfide will
generate odor and be toxic to micro and macro organisms (Shakeri Yekta et al., 2014). In
past geological epochs, before oxygen accumulated in the atmosphere, the Earth’s oceans were
stratified and anoxic. Prevalence of ferruginous or euxinic conditions thus defined the chemistry
in the Archean and Proterozoic oceans (Poulton et al., 2004, Guilbaud et al., 2015, Swanner
et al., 2020).
1.3 Hydrodynamic Models
Hydrodynamic modeling simulates the movement of water based on mass, momentum, and
energy conservation laws. By including mechanical forcing, such as wind drag, bottom current
drag, seiche activities, and by including heat fluxes, one can simulate the dynamics of water
mass (density) and temperature.
Hydrodynamic modeling (or similarly, reactive-transport modeling introduced later) can be
performed at different levels of complexity. Box models simplify the system into two or
more homogeneous blocks. One-dimensional models include continuous changes in water
properties along the vertical dimension, but ignore horizontal variability. Two dimensional
or three dimensional models allow more realistic simulations of processes, including seiche
activity, wave propagation, horizontal transport of water, coastal processes etc. Adding more
dimensions escalates the complexity of models and can exhaust computational power quickly.
1-D models are often used for studies of vertical mixing in lakes where horizontal variability is
limited, such as in small and medium-sized lakes.
8Multiple 1-D hydrodynamic lake models have been developed. Some of the popular ones
include Hostetler (Hostetler, 1987), Simstrat (Goudsmit et al., 2002), LAKEoneD (Joehnk and
Umlauf, 2001), LAKE (Stepanenko and Lykossov, 2005), FLake (Mironov, 2008), MINLAKE
(Fang and Stefan, 1996), and CLM4-LISSS (Subin et al., 2012). Lake Model Intercomparison
Project (LakeMIP) (https://www.unige.ch/climate/lakemip/) compares these models in terms
of their simulation mechanisms, accuracy, and computational efficiency (Perroud et al., 2009,
Stepanenko et al., 2010, Thiery et al., 2014). The core difference in their simulation mechanisms
lies in how the turbulence is estimated. The model being used in this project, Simstrat, applies
the k−ε turbulent closure scheme. The details of this scheme are explained in subsection 2.1.2.
Simstrat also parameterizes seiche activities in one dimension, potentially improving simulation
accuracy (Perroud et al., 2009).
Most of these modeling tools have been used for simulating mid-latitude lakes, such as in
Switzerland and the US. Studies of tropical lakes have been rare. Exceptions are African Great
Lakes: Lake Kivu (Kranenburg et al., 2020, Thiery et al., 2014, Ba¨renbold et al., 2022) and
Lake Malawi (Likongwe et al., 2015). The difficulty in simulating tropical lakes often lies in
the scarcity of direct meteorological measurements and lack of long-term continuous in-lake
monitoring data for calibration and validation.
A typical hydrodynamic model application involves tethering the simulations to observational
time series, and often iteratively correcting the model parameters to obtain a better fit to the
highly available data. In the absence of such data, however, using a model in an exploratory
fashion should be a useful way of obtaining at least qualitative understanding about the possible
behaviors of the system. Performing a sensitivity analysis by observing the change in simulated
solutions in response to parameter variations can be further used to understand the sensitivity of
the real system to changes in the lake environment.
91.4 Reaction-Transport Biogeochemical Modeling
By simulating the effects of known biogeochemical reactions, modeling provides an
inexpensive way to understand chemical dynamics in water bodies. To ensure the chemistry
is realistic, the model needs to account properly for the physical transport of chemicals
(Engesgaard and Christensen, 1988, Ajayi and Gupta, 2019). Reaction-transport chemistry
models include such physical dynamics. For meromictic lakes, the main mechanism for
chemical transport is turbulent diffusion, commonly characterized by a vertical eddy diffusivity
parameter (Kz). However, measuring Kz (Jassby and Powell, 1975) is difficult in meromictic
lakes, as it requires data such as continuous measurement of a tracer, or continuous and fine
resolution of temperature profiles within the stagnant layer (Katsev et al., 2010).
Hydrodynamic models make up for the deficiency. The k−ε turbulent closure scheme provides
the simulated dynamics of Kz, which can then be applied to the chemistry model. In return,
the outputs from the chemistry model, when compared to measured observational values, can
inform the validity of the diffusivity values simulated by the hydrodynamic model.
1.5 Tropical Ferruginous Lake Towuti
This project uses Lake Towuti as an example, to investigate the influence of climate conditions
on lake hydrodynamics, and their subsequent effects on redox chemistry dynamics.
Lake Towuti (2.75 °S, 121.5 °E) is the largest lake by surface area in the Malili Lake system on
the Sulawesi island of Indonesia (Vaillant et al., 2011). It is more than 1.2 Ma old (Damanik
et al., 2022). With an area of 560 km2 and a maximum depth of 200 m (Fig 1.1), it lies in
ophiolitic rocks and is surrounded by lateritic soils (Vuillemin et al., 2016, Hasberg et al., 2019).
The climate over Lake Towuti has an annual cycle of distinguishable dry and wet seasons. For
Sulawesi Island, The dry season typically lasts from May to October and the wet season lasts
10
Figure 1.1: The bathymetry of Lake Towuti (modified from Russell et al. (2016)).
from November to next April. The dry season, as its name suggests, has lower precipitation,
lower humidity, and lower cloud cover. It also typically comes with higher wind speed and
lower air temperature. The wet season is hotter and more humid. The air temperature varies
between 20 to 30 °C, with a mean of around 25 °C (Fig. 2.2). Lake Towuti is oligotrophic due
to nutrient limitation, although the photosynthetically active radiation is more abundant than
mid-latitude lakes. Alin and Johnson (2007) predicts that lakes in tropical regions have primary
productivity of around 1000 gC/m2/y. However, the primary productivity on Lake Towuti is
only around 100 gC/m2/y (Dr. Sean Crowe, pers. comm.).
Lake Towuti catchment area is around 4500 km2 (roughly measured with google map, Fig.
1.2). It is not much larger compared to the lake’s surface area (561 km2). Thus, the effect of
catchment runoff is less significant in affecting the hydrodynamic process in Lake Towuti. The
variation of Towuti lake level is also small. Annual fluctuation is within one meter (Fig. 1.3).
Lake Towuti is stratified primarily by a temperature gradient of a few degrees Celsius (varying
seasonally between 0.5 and 3.5 °C between surface and bottom, with only a weak salinity
11
Figure 1.2: Lake Towuti catchment (including the upstream Lake Matano).
gradient. However, since there is only weak seasonality in the equatorial region, the temperature
stratification may allow the lake to remain meromictic for decades, if not centuries or millennia,
as is the case for the neighboring 600 m deep Lake Matano (Katsev et al., 2010). Lake Towuti
has a deep thermocline that lies at around 70 m. The physical separation of the bottom waters
leads to a strong chemical stratification. Oxygen is depleted below 120 m, and dissolved iron
accumulates below 180 m (Costa et al., 2015).
On a longer time scale, the Lake Towuti Drilling project provides useful information in
reconstructing its paleoenvironment (Russell et al., 2016). Paleoclimate proxies suggest that
Lake Towuti experienced a substantially lower water level between c.a. 29 and 16 kya (Vogel
et al., 2015, Goudge et al., 2017). This reduced hydrological connectivity with other Malili
12
Figure 1.3: Towuti lake level via satellite altimetry Sulistioadi et al. (2015)
lakes (Costa et al., 2015). These events correspond to the global Last Glacial Maximum (LGM)
period. Paleolimnological evidence, such as abundant red siderite clay (Russell et al., 2020)
and Mo isotope variations (Damanik et al., 2022) suggests that the lake was mixing regularly
during that period and had abundant oxygen throughout the water column. The variations in
the lake’s past primary production might also be made possible with changing in stratification
status. Ageli et al. (2022) suggests that Lake Towuti was once dominated by productive
planktonic communities. However, under modern conditions, the primary production is now
dominated by benthic organisms. This can be explained by the establishing of stratification
under modern conditions. This stratification blocks nutrient supply towards surface plankton
communities, thus limiting the primary production.
1.6 Research Questions
In the context of meromictic tropical lakes, some of the key questions revolve around the
quantitative controls on lake stratification, and the effect of them on chemical distributions.
1. What causes a tropical lake to become permanently stratified?
13
2. What magnitudes of changes in climate parameters might cause a meromictic lake to mix
episodically, or regularly?
3. In particular, how was the physical status of Lake Towuti different during the Last Glacial
Maximum?
4. How long does it take to deplete/replenish oxygen in bottom water after prolonged
stratification/mixing?
5. How does the strength and persistence of stratification or mixing affect the redox
chemistry in the monimolimnion?
While it is generally acknowledged that greater depths and smaller surface-to-depth ratios lead
to a greater probability of meromixis, lake-specific predictions generally have not yet been
possible, including stratification or destratification of a lake in the course of its geologic history.
Answers to these questions hold keys to the interpretations of paleolimnological records in these
lakes and linking them to the contemporary environmental conditions.
In order to answer these questions, we used a hydrodynamic model to perform an explorative
study, focusing it on the possible range of hydrodynamic behaviors in Lake Towuti. By
varying the climatic conditions on model input, we analyze the likely changes in Lake Towuti
stratification throughout the glacial-interglacial transition. The drastic changes may be verified
against the available sediment record. We then use the simulation results to infer the principal
drivers for the onset of meromixis in tropical lakes.
Importantly, obtaining a realistic hydrodynamic regime provides us with a better quantification
of the vertical mixing intensities (described in one-dimensional models by the turbulent eddy
diffusivity parameter, Kz). The magnitude and depth variation in this parameter determine
the rates of chemical transport throughout the water column, shape the vertical profiles of
chemical species, and regulate their concentrations. The obtained distribution Kz(z) is then
14
used as an input to a one-dimensional biogeochemical reaction-transport model. The simulated
chemical distributions are compared against observations. Exploration of that model then
allows us to conclude on the potential transitions between oxic, ferruginous (iron-rich) and
euxinic (sulfide-rich) conditions in the bottom waters of Lake Towuti, throughout its Quaternary
geologic history, and, by extrapolation, on the control factors for ferruginous-euxinic transitions
in other water bodies.
Chapter 2
Physical Model
This chapter uses a hydrodynamic model to investigate the stability of Lake Towuti and its
sensitivity to meteorological and climate conditions.
The primary questions to be answered are:
1. What causes a tropical lake to become permanently stratified?
2. What magnitudes of changes in climate parameters might cause a stably stratified lake to
mix episodically, or regularly?
3. How was the physical status of Lake Towuti different during the Last Glacial Maximum
(LGM), a period of substantially different climate for which paleolimnological
information is available?
15
16
2.1 Methodology
2.1.1 Modeling Framework
As observational data for Lake Towuti are scarce and discontinuous, we use modeling methods
to investigate the hydrodynamics and redox chemistry. The modeling framework of this project
is visualized in the flowchart Fig. 2.1. The simulation of modern Lake Towuti conditions is
termed as the ’nominal model’, in contrast to simulations investigating the influence of changed
climate and lake level conditions. The hydrodynamic model simulates the movement of water
and the dynamics of water properties, such as temperature and velocity. The hydrodynamic
model generates an important parameter - vertical eddy diffusivity, Kz. It is subsequently used
in the chemistry model as an approximation of how chemicals get transported throughout the
water column (see Chapter 3). The result of chemistry concentrations, when being compared to
currently measured chemical values, serve as a confirmation to the hydrodynamic model results.
We assume that chemical concentrations do not affect the physical properties of the water, thus,
there is no feedback constructed from the chemistry model to the physics model.
17
Figure 2.1: The modeling process. Rounded rectangles represent process / models; Document
symbols are used to indicate inputs and outputs
2.1.2 Hydrodynamic Model
The hydrodynamics of Lake Towuti is simulated with a one-dimensional k− ε model Simstrat,
version 3.0 (Goudsmit et al., 2002, Gaudard et al., 2019). The model was shown (Perroud et al.,
2009) to be a good choice for simulating vertical temperature profiles and thermocline positions
in deep stratified lakes. It was widely used to simulate lakes in Switzerland with rigorous
calibrations (Gaudard et al., 2019). It has also been applied to meromictic Lake Kivu (Thiery
et al., 2014, Ba¨renbold et al., 2022) with success.
18
Main Equations
Simstrat uses the following equations to describe temperature T [°C], mean horizontal velocities
u and v [ms−1], turbulent kinetic energy (TKE) per unit mass k [Jkg−1] and TKE dissipation
rate ε [Wkg−1]. Depth z [m] is defined as 0 at the surface and negative downwards Goudsmit
et al. (2002).
∂T
∂ t
=
1
A
∂
∂ z
(A(ν ′t +ν
′)
∂T
∂ z
)+
1
ρ0cp
∂Hsol
∂ z
+
dA
dz
Hgeo
Aρ0cp
(2.1)
∂u
∂ t
=
1
A
∂
∂ z
(A(νt +ν)
∂u
∂ z
)+ f v (2.2)
∂v
∂ t
=
1
A
∂
∂ z
(A(νt +ν)
∂v
∂ z
)− f u (2.3)
∂k
∂ t
=
1
A
∂
∂ z
(Aνk
∂k
∂ z
)+P+Pseiche+B− ε (2.4)
∂ε
∂ t
=
1
A
∂
∂ z
(Aνε
∂ε
∂ z
)+
ε
k
(cε1(P+Pseiche+ cε3B− cε2ε)) (2.5)
Here, ρ0 [kgm−3] and cp [Jkg−1 K−1] are reference density and specific heat of water. A [m2]
is the surface area of lake at depth z [m], Hsol [Wm−2] is the shortwave solar radiation flux into
the lake, Hgeo [Wm−2] is the geothermal heat flux, and f [s−1] is the Coriolis parameter. ν and
νt are molecular and turbulent viscosity for momentum, ν ′ and ν ′t are molecular and turbulent
diffusivity of temperature, and νε and νk are the turbulent diffusivity of energy dissipation and
TKE. These viscosity / diffusivity terms all have units of [m2 s−1]. P and B [Wkg−1] represent
19
the kinetic energy generated from shear stress production and buoyancy production:
P = νt((
∂u
∂ z
)2+(
∂v
∂ z
)2) (2.6)
B =−ν ′t N2 (2.7)
where N [s−1] is the buoyancy frequency defined as:
N =− g
ρo
∂ p
∂ z
. (2.8)
The model also simulates the effect of seiches on turbulent kinetic energy by adding a
production term Pseiche in equation 2.4 and 2.5. This term is an empirical fit of seiche energy as
a portion from the surface and bottom stress, discussed in detail by Goudsmit et al. (2002).
The molecular viscosity for momentum is:
ν = 1.5×10−6 m2 s−1 (2.9)
The molecular diffusivity for temperature is:
ν ′ = 1.5×10−7 m2 s−1 (2.10)
Finally, the turbulent viscosity and diffusivity for momentum, temperature, k, and ε can be
calculated with Kolmogorov-Prandtl relation. They are all linearly related to k
2
ε . Turbulent
viscosity for momentum is:
νt = cµ
k2
ε
(2.11)
20
Turbulent diffusivity for temperature is:
ν ′t = c
′
µ
k2
ε
(2.12)
Turbulent diffusivity for k is:
νk =
cµ
σk
k2
ε
(2.13)
Turbulent diffusivity for ε is:
νk =
cµ
σε
k2
ε
(2.14)
The values of the constants, cµ , c′µ , σk, and σε are given in Table B.1, along with other model
parameters.
In summary, this model provides a 1-D approximation of mixing based on viscosity and
diffusivity of momentum and temperature, which is scaled based on turbulent kinetic energy
and its dissipation rate. Convection of water and the associated convective mixing is represented
with greater diffusivity terms.
Simulation Configuration and Input
Inputs to the hydrodynamic model include lake bathymetry, initial distributions of temperature
and salinity, light attenuation, hydrological inputs, and meteorological variables.
Lake bathymetry is defined by a hypsographic curve, the areal extent as a function of depth,
at 5 meter intervals, interpolated from Fig. 1.1. The depression in the southern basin below
150 meter was excluded from the areal extent calculation because it is effectively isolated
from the main deep basin. Several model runs suggest that the simulated stratification was
not significantly affected by the details of lake morphology or by excluding the separate deep
basin.
21
Initial conditions required by Simstrat includes water column temperature t [°C], salinity S [‰],
horizontal velocity defined by u and v components [ms−1], TKE per unit mass k [Jkg−1], and
TKE dissipation rate ε [Wkg−1]. Initial horizontal velocities are set to 0 [ms−1], initial k is set
to 3×10−6 [Jkg−1] and initial ε is set to 5×10−10 [Wkg−1].
The initial temperature and salinity are obtained from CTD measurements in 2001. CTD
profiles for multiple years were obtained from Drs. Sean Crowe and Doug Haffner (personal
communication). These CTD profiles consist of temperature and conductivity measurements.
The density of water was calculated via the equation of state (Chen and Millero, 1986):
ρ0 = 0.9998395+6.7914×10−5t−9.0894×10−6t2+1.0171×10−7t3−1.2846∗10−4t4
+1.1592×10−11t5−5.0125×10−14t6+(8.181×10−4−3.85×10−6t+4.96×10−8t2)S
(2.15)
where ρ0 is the density of water at the surface with zero water pressure, aka potential density.
t and S are temperature and salinity at each depth. Salinity is approximated from the measured
conductivity using a simple scale function from Katsev et al. (2010) for the neighboring Lake
Matano:
Salinity [g kg−1] = 0.0010×Conductivity [µS cm−1] (2.16)
Here, conductivity is in [µScm−1], and the resultant salinity is in [gkg−1].
For the initial conditions, the temperature and salinity were extrapolated from the closest
available measured profiles.
The coefficient of light attenuation in the water column was estimated from the profile
of the photosynthetically active radiation (PAR) measured in 2015 (Sean Crowe, personal
22
communication). This was done by linear fitting the log-transformed profile between 17 and
110 m depth. Lake Towuti’s photic zone is very deep. Light attenuates to 1 % at approximately
60 m. This is a result of Lake Towuti being an oligotrophic lake, with few particles to scatter
light. In the simulations, light attenuation was assumed to be constant throughout seasons. The
magnitude of light attenuation strongly affects the depth of the simulated thermocline and is
potentially important for the agreement between model and observations.
River inflow and outflow data are unavailable, but since Lake Towuti has very small catchment-
to-lake area ratio, the hydrodynamics of the lake is not highly sensitive to river discharges.
Several trials were conducted with different orders of magnitude of discharges (20, 7×107
m3 s−1). The mixing pattern was not significantly affected. Thus, the default discharge values
were fixed at 20 m3 s−1 for both river inflow and outflow.
Surface forcing parameters required by Simstrat include zonal and meridional wind speeds,
surface air temperature, surface solar radiation (after atmospheric absorption/reflection), vapor
pressure, and cloud cover. These variables were obtained from ERA5 reanalysis dataset
(Hersbach et al., 2020) with hourly resolution. Care was taken to only select the available land
grids around Lake Towuti (121.4 °E, 2.75 °S) to avoid the nearby ocean climate. The time
series of average ERA climate variables are presented in Fig. 2.2.
The monthly average of ERA5 reanalysis dataset is compared with observational data in Fig.
2.3. Meteorological observations were available for Kendari (166 km away) and Palu (300 km
away) city. ERA5 reanalysis has a similar air temperature range. The wind speed is half or a
quarter of that being observed in the two cities. Relative humidity from ERA5 is lower than
observations during July - October. The surface shortwave radiation from the reanalysis data
has a similar trend with the sunshine minutes observed in Kendari. The cloud cover is higher
for ERA5 than observed. The total precipitation is within the range between Kendari and Palu.
The reanalysis data is the best choice for modeling Lake Towuti, as there is no direct
23
measurement of meteorology over the lake that can satisfy the time resolution requirement
for model simulation. Although discrepancies exist between observational and reanalysis
meteorology, there is little confidence in deciding which one best represents the real condition
over the lake surface. Both of the observed weather stations are far away, and are located
in different geographical environments (both Kendari and Palu are coastal cities at lower
elevations).
Adequate results were obtained from the uncorrected reanalysis data, except the absorption of
shortwave and longwave radiation are downscaled by a constant factor (App. B, see psw water
and plw). Other constants used in the model were set so that their values were deemed to
best represent tropical lakes (App. B). Otherwise no additional calibration and correction were
conducted due to the lack of long-term continuous monitoring data.
24
Figure 2.2: ERA5 reanalysis climate variables from 1989 to 2016, marked with ENSO events
indicated by the MEI.v2 index: El Nin˜o (MEI > 0.5) are marked with red bands; La Nin˜a (MEI
< -0.5) are marked with blue bands.
25
Figure 2.3: Comparison of the annual climate data around Lake Towuti from different
sources. ERA5 reanalysis data are presented as mean monthly averages derived for 2000-2016.
Observational data are from weather stations in Kendari and Palu city (https://www.weather-
atlas.com/, year range unknown).
26
2.1.3 Stratification Strength Indicators
We define several indicators to quantify the strength of stratification of the entire water column
of the lake. We use temperature difference indicator (TDI), heat content indicator (HCI), and
Schmidt Stability index (SSI) as characteristics of the entire water column at any given time.
The weakest temperature difference (WTD) and maximum buoyancy frequency (MBF) are used
to characterize the weakest and the strongest instances of stratification over the entire simulated
period.
The temperature difference indicator (TDI) is calculated as the difference between the maximum
temperature and minimum temperature in a water column. This is typically the difference
between lake surface temperature and bottom temperature. A TDI around zero indicates that
the water column has uniform temperature, indicating a whole water column mixing.
T DI(t) = max
z
(T (z, t))−min
z
(T (z, t)) (2.17)
The Weakest Temperature Difference (WT D) is defined as the smallest temperature difference
indicator (TDI) in the water column throughout the whole simulated time series:
WT D = min
t
(T DI) (2.18)
A WTD of zero would indicate the existence of a homogeneous temperature distribution
throughout the water column at least once during the simulation, indicative of the whole-water
column mixing.
The Maximum Buoyancy Frequency (MBF) is the maximum Brunt-Va¨isa¨la¨ frequency squared
(N2) anywhere in the water column. MBF characterizes the stability at the thermocline. Higher
27
MBF indicates a more stable thermocline being achieved during the simulation run:
MBF = max
t,z
(N2(z, t)) (2.19)
The Schmidt Stability (SS) characterizes the stability of the water column based on the density
gradient. It is defined as (Idso, 1973) :
SS =
g
A0
∫ z=0
z=zd
(z− zg) Az (ρ− ρ¯) dz (2.20)
It represents the energy needed to mix the water column from the surface down to a depth zd .
The Schmidt Stability index (SSI) is the SS taken to the maximum depth. Here A0 is the lake
surface area, Az is the areal extent at depth z, and zg is the depth of the center of the mass,
calculated by:
zg =
∫ z=0
z=zd ρ Az z dz∫ z=0
z=zd ρ Az dz
(2.21)
and ρ¯ is the average density of the water column, calculated by:
ρ¯ =
∫ z=0
z=zd ρ Az dz∫ z=0
z=zd Az dz
(2.22)
The Heat Content (HC) characterizes the susceptibility of the water column to mixing under
changes in heat fluxes. It is defined as the water heat content in the water column above depth
zd :
HC =
1
A0
∫ z=0
z=zd
ρ cw (T −T0) Az dz (2.23)
28
where T is the in situ temperature. T0 is an arbitrarily chosen reference temperature (28 °C
if not state otherwise), cw is the specific heat of water, 4200 [Jkg−1 °C−1]. The Heat content
indicator (HCI) is the heat content of the entire water column.
To eliminate the effects of the initial spin up period for the model, the first year of the simulation
was usually excluded from the indicator analyses.
2.1.4 Sensitivity Test
A Monte Carlo method was used to test the sensitivity of the model simulated stability to
climate variables. Multiple climate variables were varied randomly at the same time, and the
same model simulations were repeated for 100 randomized parameter sets. The parameters
were drawn from normal distributions, which covered the realistic ranges of the respected
variables. The standard deviations (σ ), chosen to represent relevant environmental ranges, are
listed in Table 2.1. For the wind speeds u and v, the normally-distributed randomization factors
were used multiplicatively, rather than additively. The average values of the meteorological
variables are calculated and compared with stratification indicators (WTD and MBF) calculated
for each model run. As the climate variables may be coupled among themselves in reality,
the distributions that characterize their real variability may deviate from the independently
chosen normal distributions used here. The Monte Carlo analysis serves as an exploration for
understanding the effects of climate variables on hydrodynamic circulation in tropical lakes.
Variables distribution standard deviation unit
Wind speed multiply Gaussian 0.4 m s−1
Air temperature Gaussian 1.5 °C
Vapor pressure Gaussian 1 mbar
Cloud cover Gaussian 0.1 %
Table 2.1: Variations in meteorological variables used in the Monte Carlo analysis.
29
2.1.5 Correlation with the ENSO Climate Index
Located in the tropical west Pacific, Lake Towuti might be susceptible to the variations of
large climate cycles, such as El Nin˜o-Southern Oscillation (ENSO). The correlation of the
ENSO index, MEI.v2 (https://psl.noaa.gov/enso/mei/) and the model simulated SSI has been
investigated by cross correlation. This is conducted by MATLAB normalized cross correlation
function, with inputs of the MEI.v2 index and SSI(t) from the nominal Simstrat simulation run
(2-hour resolution). The normalized cross correlation generates Pearson correlation coefficient
r as a function of lag. Positive lag indicates MEI.v2 prior to SSI. Positive correlation coefficient
indicates the El Nin˜o corresponds to strong SSI or La Nin˜a corresponds to weak SSI. A negative
correlation coefficient indicates the opposite effect, that is, El Nin˜o corresponds to weak SSI or
La Nin˜a corresponds to strong SSI.
2.1.6 Lake Towuti in the Last Glacial Maximum
Stratification of Lake Towuti during the Last Glacial Maximum was simulated to establish
correspondence with existing sediment records (Russell et al., 2020). During the Last Glacial
Maximum (LGM), the average air temperature was estimated to be 4 °C cooler, relative
humidity 15 % lower, and wind speed 2 – 3 times higher than modern conditions (estimated by
Dr. James Russell, personal communication).
Simulation Climate Lake Level
1 Modern 100
2 Modern 150
3 Modern 190
4 LGM 100
5 LGM 150
6 LGM 190
Table 2.2: Simulation runs for investigating the effect of LGM conditions.
30
2.2 Results
2.2.1 Stratification of Lake Towuti
The simulated temperature distributions reproduce the essential features of the observational
data, despite the input data being approximated from reanalysis with minimal calibration.
Simulated water temperatures range between 28 °C at the bottom and 30 °C at the surface,
matching the observational ranges. The seasonal thermocline deepens during the wet season,
from 40 to 60 m deep. The model somewhat underestimates surface water temperatures
and overestimates the depth of the thermocline compared to the mooring observations from
2012-2015 (Fig. 2.5). The overall differences between the mooring data and the simulation,
however, do not exceed 1 °C (Fig. 2.6).
Figure 2.4: Temperature distributions and evolution in Lake Towuti simulated using the
unmodified ERA5 reanalysis meteorological data as model input.
The long term simulation (Fig. 2.4) showed some temperature variability in the deep
monimolimnion that was not captured by shorter term mooring observations (Fig. 2.5) but was
31
Figure 2.5: Comparison of simulations with the temperature data from a moored instrument
array (Dr. James Russell, pers. comm.; observational data after June 2014 have less confidence
due to instrumentation problems).
recorded in the discrete CTD profiles taken over 2000-2009 (Fig. 2.7). Lake Towuti’s bottom
layer was previously assumed to have been isolated from surface waters, with heat exchange
occurring only slowly, via diffusion across a stable thermocline. However, model outputs
suggest one or more potential mixing events. The model predicts mixing in 2001 and 2002,
with temperature at the bottom of the lake cooling down as a result. From the 2002 dry season
onward, stratification became established and the bottom water slowly warmed up, as heat
diffused downward from the epilimnion. In the 2007/2008 dry season, the water column mixed
again, lowering the bottom lake temperature by almost half a degree, from the original 28 °C
to 27.5 °C. The simulated temperature profiles from these cooling and warming events seem to
match the available CTD profiles from 2000 to 2009 (Fig. 2.7). As the CTD profiles represent
32
Figure 2.6: Residuals from the model (model minus mooring temperature). Warm color
indicates warmer simulated temperature, colder color indicates colder simulated temperature
compared to the observation
discrete moments of time and were not taken in the same seasons, inferring changes from
such observations alone is difficult. The simulations illustrate the corresponding dynamics that
happened over multiple seasons. In particular, they show that the cooler water temperature in
the deep water observed in the 2009 CTD profile compared to the 2006 and earlier profiles is
most likely the result of a deep mixing event in the dry season of 2008.
The simulation provides information on the dynamics of the heat content and the stability of the
lake (Fig. 2.8). A complete mixing of the lake would homogenize the temperature in a water
column, corresponding to zero TDI and zero SSI. A complete mixing should be associated with
negative surface heat flux, which cools down the entire water column. Thus, the average HCI
would also decrease. After the 2008 mixing event, TDI and SSI rebounded and their average
values increased compared to the pre-2008 state. This indicates that the thermal gradient in the
lake increased after mixing. The more stable stratification lasted until the end of the simulation
33
Figure 2.7: Comparison between observational CTD temperature profiles and modeled profiles
at the corresponding time. Vertical dashed lines represent 28.5 °C. It is plotted to provide an
easy comparison between different CTD profiles.
time frame. Cooling of the bottom water thus made the lake more difficult to be fully mixed
again.
Simstrat also provides information on the lake surface heat flux. Fig. 2.9 shows the heat flux
in a typical stratified year (2012). Shortwave and longwave radiation from the sky provide the
positive inflow of heat. Longwave radiation out from the lake, sensible heat flux, and latent heat
flux make up for the negative outflow of heat. Overall, the seasonal variations of shortwave
and downward longwave radiation are the most significant components of the heat flux balance,
followed by the variations in latent heat flux.
The surface heat fluxes are affected by the meteorological variables. Surface shortwave
radiation, which is the incoming radiation from the sun, depends on the time of year and the
reflection and absorption of light by the atmosphere. Thus, when the cloud cover or relative
humidity are higher, the amount of shortwave radiation reaching the lake surface is reduced.
This leads to lower shortwave radiation inputs from April to July, and higher radiation from
August to November. Longwave radiation from the sky (or longwave in) is linked to air
temperature and humidity. Higher air temperature and humidity are usually observed during
the wet season. Thus, incoming longwave radiation is greater during the wet season, between
34
Figure 2.8: The simulated time series for the overall stratification indicators: temperature
difference indicator (TDI), maximum buoyancy frequency (MBF), heat content indicator (HCI),
and Schmidt Stability index (SSI).
November and June of the following year. Outbound longwave radiation emitted by the
lake follows Planck’s law (black-body radiation), and is only related to the surface water
temperature. Since the surface water temperature is warmer from November to April, that is
when the longwave radiation out from the lake is greatest (more negative). Sensible and latent
heat are determined by wind speed, air temperature, water temperature, and humidity (Imboden
and Wu¨est, 1995). As water temperature is always warmer than air temperature, these two
fluxes are always negative. Greater wind speed and lower humidity usually occur during the
dry season, leading to more negative sensible and latent heat fluxes. However, during the dry
season, the air temperature and water temperature are both cooler, leading to an opposite effect.
35
On balance, air temperature and water temperature exert a more significant influence on the
sensible and latent heat fluxes.
Figure 2.9: Heat fluxes across the lake surface calculated by Simstrat in 2012. From top to
bottom: shortwave radiation, longwave radiation reflected from sky into the lake (longwave in),
longwave radiation from lake into the atmosphere (longwave out), sensible heat, latent heat,
and the net heat flux. Red shaded regions mark a typical Sulawesi wet season, and blue shades
indicate a typical dry season.
Fig. 2.10 presents the correlation between the monthly averages of the climate variables from
36
the ERA5 reanalysis (same as the model inputs) and the monthly averages of the SSI. Each data
point represents a monthly sample from the nominal simulation. The corresponding correlation
coefficients and p-values are listed in Table 2.3 (note that the data points, being from the same
simulation, are not statistically independent, hence these statistical measures are provided only
for illustration). The Schmidt Stability index correlates with most of the climate variables.
The correlation is stronger with air temperature, vapor pressure, and cloud cover, and weaker
with wind speed. The correlation with solar radiation is not significant. This suggests that the
incidents of weaker stratification (low SSI) are associated mostly with cooler air temperatures
and drier conditions. Higher winds also occur at the times of weaker stratification. Note that
the correlation here does not necessarily indicate that a particular climate variable controls the
stability (causality), as they co-vary among themselves. In the Monte Carlo simulations in Fig.
2.13, where values of the climate variables were varied independently, vapor pressure and cloud
cover by themselves did not significantly affect lake stability. Lower vapor pressure and lower
cloud cover are, however, associated with the period’s lower air temperatures (dry season) in
the Lake Towuti region.
SSI correlated with r p
Wind Speed -0.245* 0.00001
Air Temperature 0.696* 0.00000
Solar Radiation -0.072 0.19087
Vapor Pressure 0.826* 0.00000
Cloud Cover 0.653* 0.00000
Table 2.3: Correlation between ERA5 climate variables and the SSI of Simstrat simulation run
with them as inputs.
We also inspected the connection between the strength of the lake stratification and the regional
climate variations. The response of lake stratification to large-scale variations, specifically, El
Nin˜o-Southern Oscillation (ENSO), is not straightforward. (Fig. 2.11). Several SSI troughs do
follow the periods of negative MEI.v2 index, such as from year 2000 to year 2008. However,
37
Figure 2.10: Scatter plots of SSI against model input climate variables (from ERA5 reanalysis).
Because the (monthly) data points correspond to monthly averages, a brief mixing of the water
column corresponds to low values of SSI but not necessarily zero.
anti-correlations between SSI and MEI.v2 index are also observed, such as from 1990 to 1995.
Similarly, periods of high stability (year 1997, Year 2003, Year 2010, and year 2016) correspond
with positive MEI.v2 index (i.e. El Nin˜o periods.), though some (shorter) periods of high
MEI.v2 did not apparently lead to greater stability (years 2005).
Fig. 2.12 shows the cross-correlation function between the ENSO index and the simulated
38
Figure 2.11: The MEI.v2 index of ENSO (> 0.5 - El Nin˜o, <−0.5 - La Nin˜a) and the moving
(12 months) mean of the Schmidt Stability (SSI) from the simulation.
lake stability, represented by SSI. The horizontal axis is the time lag in months. Negative
time lag corresponds to the ENSO indicators leading the SSI. Positive correlation coefficient
would indicate the El Nin˜o events being associated with higher lake stability and/or La Nin˜a
events being associated with lower stability. Negative correlation would indicate the opposite
relationship. The cross-correlation shows that the higher (lower) values of the MEI.v2 index
corresponds to increased (decreased) stability about four months later. However, the correlation
is not strong, for reasons discussed later (see Discussion 2.3.5). The correlations between
MEI.v2 index and the individual climate variables are included in Fig. D.4 for reference.
39
Figure 2.12: Normalized cross correlation between the MEI.v2 index and unfiltered Schmidt
Stability index (SSI) (all detrended). This figure shows that the high MEI.v2 indexes, i.e. El
Nin˜o events, are correlated with weaker stability simulated in the model, and vice versa for La
Nin˜a events. The maximum positive correlation occurs with a -4 month lag with the correlation
coefficient of 0.197. The maximum negative correlation occurs with a 14 months lag with the
correlation coefficient of 0.134. These are not strong correlations. If ENSO climate variation is
a precursor for the lake status change, the effect occurs approximately 4 months later.
2.2.2 Sensitivity Test Results
The Monte Carlo sensitivity analysis illustrates the effects of individual meteorological
variables on lake stratification (Fig. 2.13, 2.14). Wind speed and air temperature have the
most significant effects on Lake Towuti’s stratification strength. As expected, average wind
speed correlates negatively with the stratification strength while air temperature has a positive
correlation. In contrast to the single-simulation results in Fig. 2.10, vapor pressure and
cloud cover individually do not have strong effects on stratification, when their effects are
investigated independently in Monte-Carlo runs of the model. The sensitivity test suggests
the typical magnitudes of climate variables that would be required for long-term stratification
or, alternatively, might lead to whole-water column mixing (as indicated by WTD). For Lake
Towuti, the lake is likely to experience complete mixing if the average wind speed is above 2
40
ms−1 (doubled from ERA5 averages), or when the average air temperature is below 23 °C (1–2
°C cooler than ERA5 averages). On the other hand, if the average air temperature exceeds 27
°C, the stratification of the lake is likely to persist for the whole simulated period. Higher air
temperatures also lead to stronger thermocline during the stratified wet season (as indicated by
MBF). Wind has a weak negative effect on the strength of the thermocline. In this project, wind
speed is used to correlate with lake stability. However, in the future, it will be more sensible to
use wind power (related to cubed wind speed) for the purpose of this comparison since wind
power more linearly controls energy fluxes.
Figure 2.13: Correlation between the MBF indicator for the maximum strength of the
thermocline (dots), the WTD indicator for the minimum temperature difference (circles), and
the mean values of climate variables, from 100 runs of Monte Carlo simulations (2001 - 2008
period). The wind speed and air temperature have the most significant impact on the stability
of stratification.
Figures 2.13 and 2.13 compare the stability measures before and after the putative mixing event
of 2008. The second figure excludes the extreme results from mixing in 2008. Although 68 out
of 100 simulations that were run for the time period between 2001 and 2016 displayed mixing
41
Figure 2.14: Correlation between MBF (dots), WTD (circles), and the mean values of climate
variables, from 100 runs of Monte Carlo simulations (2009 - 2016 period). Notice that a lot of
simulations’ WTD between 2009 - 2016 period is non-zero.
(zero temperature gradient at some point during the simulation), almost all of these instances
occurred prior to 2008. Only three out of 68 simulation runs had mixing after 2009 (Fig. 2.13,
2.14, 2.15). Prior to the strong climatic event in 2008 that promoted most vigorous mixing, the
stability of the lake was weaker (Fig. 2.8). After that, the stability of Lake Towuti increased,
which made subsequent mixing events more difficult due to increased density gradient resulting
from cooler bottom water similar to the results seen in Fig. 2.2 and 2.8.
Fig. 2.15 illustrates the differences between these two different time ranges (before and
after 2008), as reflected in each of our two stability measures. The WTD values plot above
the 1:1 line, which indicates weaker temperature differences during the dry season (more
prone to mixing) before 2008 than after 2009. The MBF values generally plot below the 1:1
line, indicating that the wet seasons before 2008 were somehow more conducive to stronger
thermoclines.
42
Figure 2.15: WTD and MBF indicator comparison between pre 2008 and post 2009 from
simulation. With criteria of WTD < 0.4 °C indicating mixing, there are 68 simulations
simulated mixing before 2008, but only 3 of them also simulated mixing between 2009 - 2016
2.2.3 Lake Towuti during the Last Glacial Maximum
Simulations of Lake Towuti stratification under the hypothesized climatic conditions of the Last
Glacial Maximum (LGM) show that lowering of the lake level, combined with colder and drier
conditions, changes the stratification dynamics (Fig. 2.16). In a shallower lake, turbulence from
the interactions of currents with the lake bottom may increase the mixing rates in the interior
of the water column. Simulations suggest, however, that under the modern climate conditions
lowering the lake level alone would not significantly influence the whole water temperature
and the thermocline depth. As the lake level varied from 200 m to 100 m in the model, under
the modern climate conditions, the temperature at the surface and the thermocline depth did
not change significantly (Fig. 2.16 left column). The surface processes were still isolated
from the bottom water. However, under the LGM climate conditions, which produced stronger
negative heat fluxes at the surface, lowering lake level to 150 m resulted in a whole water column
mixing almost every year (Fig. 2.16 right column). Lowering lake level also increased surface
43
temperatures and maximum thermocline depths during wet seasons under the LGM climate
condition. The difference between the air temperature and the surface water temperature is kept
at four °C under both the modern and the LGM climate.
Figure 2.16: Simstrat simulation exploring conditions during the Last Glacial Maximum (LGM)
for different water levels (100 m, 150 m, and 190 m) and for different climate conditions
(left column : current ERA5 reanalysis climate conditions, right column : LGM climate
conditions described in subsection 2.1.6). The simulated processes are similar to the nominal
hydrodynamic simulation for the modern lake, but the outputs are cropped so that it represents
a relatively steady stage of establishing LGM hydrodynamics.
44
2.3 Discussion
2.3.1 Multi-year Stability and Potential Oligomixis of Lake Towuti
Contrary to current beliefs about the meromixis in the lake, Lake Towuti may experience
whole-water column mixing, i.e. it may be oligomictic rather than meromictic. The mixing in
Lake Towuti, as in other similar tropical lakes, is mostly induced by strong, negative surface
heat fluxes during the dry season. These heat fluxes continuously lower the mixolimnion
temperature to that comparable to monimolimnion, until they induce mixing. Simulations,
including sensitivity tests, show that the mixing regime of Lake Towuti may be borderline
between meromixis and oligomixis. If Lake Towuti is truly oligomictic, mixing can be expected
to happen about once a decade, at the times when the thermal gradient is weak and weather
conditions are especially cold and dry. Mixing events are unlikely to be very frequent, as cooling
of the bottom waters during mixing establishes a stronger temperature gradient in the following
season, while diffusion of heat from the mixolimnion to reheat the monimolimnion takes several
years (see below).
Our simulation (Fig. 2.4) suggests a possibility of a series of brief mixing events between 2000
and 2008. After a major mixing event in 2008, the bottom water cooled down significantly,
setting up for increased stability thereafter. The lake was more stably stratified from 2009
onward.
Evidence for deep mixing also can be seen in the observational data. CTD profiles from 2000
to 2009 (Fig. 2.7) show more variability in deep water temperatures than could be explained
if Lake Towuti were permanently stratified. In particular, temperature in the bottom water
had decreased in 2009 compared to 2006 (Fig. 2.7). In the absence of intrusions of cold
groundwater, this could only happen if the surface water cooled sufficiently to sink down into
the monimolimnion.
45
After deep mixing, stratification reestablishes after receiving enough surface warming during
the wet season. After restratification, the lake surface again becomes isolated from the
lake bottom. This isolation, however, does not inhibit the gradual diffusion of heat into the
monimolimnion. The overall picture of lakes’ mixing-stratification dynamics is discussed
below.
2.3.2 Mechanisms of Stratification and Mixing
Strength of Stratification
The strength of stratification in Lake Towuti can be characterized by the cumulative Schmidt
Stability (SS) and cumulative heat content (HC). Fig. 2.17 illustrates the typical SS and HC
profiles that were calculated based on available CTD and mooring data. During the wet season,
cumulative SS increases from 0 at the surface to 5 kJm−2 at around 100 m, to more than 10
kJm−2 at the maximum mooring depth 150 m. During the dry season, the magnitude of SS
is much smaller, at around 1 kJm−2 at 150 m. The magnitude of cumulative SS indicates the
mechanical energy required to homogenize the water column down to a given depth. During the
wet season, cumulative HC is 0 at the surface, 4×105 kJm−2 at the 100 m and 4.5×105 kJm−2
at 150 m. The magnitude of cumulative HC indicates the negative thermal energy required to
remove the heat and cool down all the water above the given depth.
We can compare these numbers, respectively, with the typical wind power and surface heat
fluxes at Lake Towuti. The wind power per unit area being transferred into a lake (Pkin,W ) can
be estimated using the approximation of Imboden and Wu¨est (1995):
Pkin,W = χρairC10W 310 (2.24)
Here, χ is the fraction of wind energy at 10 m above the surface being transferred into the
46
lake, typically 1–2%. ρair is the density of air, 1.2 kgm−3. C10 is the drag coefficient, which
depends on the lake surface roughness. It is typically 0.001 for wind speed smaller than 7 ms−1
(Imboden and Wu¨est, 1995). Thus, for an average wind speed (10 m above surface, W10) around
0.85 ms−1, the typical wind power transferred into the lake is around 1×10−5 Wm−2. For a
three-month season, the total wind energy is on the order of 0.1 kJm−2. By comparing this
value to a typical SS profile (Fig. 2.17), the mixing power induced by pure wind drag can only
mix the water column down to about 30 meters in a three month season.
Lake Towuti’s seasonal heat fluxes simulated by Simstrat are shown in Fig. 2.9. The net
surface heat flux varies between -150 Wm−2 from June to September to +150 Wm−2 from
September to January. The negative (positive) heat fluxes during the dry (wet) season is on the
order of 1×105–1×106 kJm−2. The annual heat imbalance is on the order of 1×105kJm−2.
Comparison with the HC in Fig. 2.17 shows that such heat fluxes are capable of extending
mixing substantially downward into the monimolimnion during the dry season, potentially even
to the lake bottom in more extreme conditions.
Thus, variations in heat fluxes are the major mechanism for stabilizing and destabilizing the
water column. Mixing induced by direct wind turbulent energy is limited in comparison, as
is commonly recognized (Imboden and Wu¨est, 1995). While these numbers serve as a simple
illustration of the relative importance of mechanical and thermal energies (Imboden and Wu¨est,
1995), the processes that distribute both mechanical and thermal energies throughout the lake
and their associated time scales are also important. They include seiche activities, turbulent
energy dissipation, and non-linear variations in temperature being an active tracer.
Heat Diffusion
Even a strongly stratified water column would not stop the exchanges of heat and substances
between the surface and bottom layers. The dominant transport process is diffusion. Molecular
47
(a) Cumulative stability (SS) and heat content (HC) calculated based on CTD data.
(b) Cumulative SS and HC calculated from mooring data
Figure 2.17: Schmidt Stability and Heat Content (relative to the water temperature at
the calculated depth). The magnitude of SS indicates the mechanical energy required to
homogenize the water column to the desired depth while the magnitude of HC indicates the
negative thermal energy required to cool down all the water above that depth.
diffusivity (commonly 1×10−9 m2 s−1) is usually too slow to significantly affect the molecule
concentration profiles in stagnant water columns. The diffusivity of heat has a greater magnitude
(commonly 1×10−6 m2 s−1). The turbulent diffusion, even during the stratified season in Lake
Towuti (1×10−6 m2 s−1), in contrast, increases the magnitude of vertical diffusivity for both
48
chemicals and heat.
The depth and time scales for the diffusion process in the monimolimnion can be estimated by
a back-of-the-envelope calculation using the Einstein’s equation for Brownian motion:
x¯2 = 2Dt (2.25)
The turbulent diffusivity under stratified conditions in Lake Towuti is D≈ 1×10−6 m2 s−1 (Fig.
3.1). Under this rate of diffusion, heat can be transported 60 meter downwards within a year:
xmax =
√
2∗10−6m2s−1 ∗31536000s≈ 60m (2.26)
This process can gradually warm up the monimolimnion temperature, and reduce the
temperature difference between the lake’s surface and bottom. The time scale for this
diffusional process to restore pre-mixing temperature gradient can be estimated:
∆T
∆t
≈ Kz∆T∆z
∆t ≈ 0.5°C/(1×10−6m2 s−1 1°C
(20 m)2
)
≈ 1×108 seconds≈ 3 years
(2.27)
This diffusional process transports heat into the monimolimnion. However, because the surface
water temperature is balanced with the surface atmospheric conditions at a much faster rate,
the state of the mixolimnion does not experience drastic changes in its average temperature.
Once the mixolimnion supplies heat to warm up the monimolimnion, the balance between
mixolimnion and surface weather condition will restore the surface water temperature, via
49
decreased negative fluxes of blackbody radiation, sensible and latent heat. These processes keep
mixolimnion warm as usual, without being significantly disturbed by the cold monimolimnion.
Increasing global average air temperature might further warm up the surface temperature, whose
effect will be discussed in detail later.
The Effect of Precipitation
The momentum and heat effect of the precipitation is not simulated in Simstrat. However, a
back-of-envelope calculation suggests that the effect of surface precipitation is small. The heat
content (per unit area) of the precipitation (reference to 28 °C) in one hour, 5 mm/h rainfall
event is:
HC/Area = ρw× rain f allDepth× cp×∆T
= 996 kgm−3×5×10−3 m×4189 J°C−1 kg−1× (21−28) °C
=−1.46×105 Jm−2
=−1.46×102 kJm−2
(2.28)
Compared to the value magnitude of heat content plot in Fig. 2.17, the cooling power of the one-
hour precipitation is not strong enough to destroy a stable thermocline unless the precipitation
is 100 to 1000 fold greater than 5 mm/hr or it lasts for a much longer time (100 to 1000 hours).
Stratification Mechanism
The stability of Lake Towuti is determined by the temperature gradient, as the salinity is very
low (conductivity of 140 - 195 µScm−1 is equivalent to salinity 1.4×10−5–1.9×10−5 gkg−1)
and has only a weak gradient. In the absence of a salinity gradient, the water temperature
gradient plays a key role in the stratification and mixing of Lake Towuti.
50
Seasonal climate patterns are the main cause of the variations in water temperature, and thus,
lake stability. Wet seasons come with higher air temperature, lower wind speed, higher vapor
pressure, and higher cloud cover, resulting in positive surface heat flux. This increases the heat
content in the lake’s epilimnion and establishes stronger stratification. The opposite happens
in the dry season, when air temperature decreases, wind speed increases, and humidity and
cloud cover decreases, leading to negative surface heat flux. This process cools down the lake
surface, reduces the overall temperature gradient, induces mixolimnion mixing, and destabilizes
the lake.
Mixolimnion mixing almost always starts in the dry season, when naturally the negative heat
flux cools down the surface. During extreme climate events, mixing can extend to the whole
water column. Over the multi-year time scale, air temperature influences lake stability the most.
Average air temperatures below 24 °C would result in mixing during the simulated period (2000
- 2016). Wind speed has the potential to destabilize the water column as well, but at more
extreme values. In Lake Towuti, if the average wind speed is greater than 2 ms−1 (double of
current wind speed), the lake is likely to mix. Wind affects the stability by both mechanical and
thermal effects, as higher wind decreases heat fluxes and induces instability.
Besides immediate climate variations over the lake, lake history also plays an important role in
how the lake reacts to changes in climate. Climate variables influence the heat flux at the surface
of the lake. A negative heat flux would only lead to mixing when it is sufficient to cool down the
whole mixolimnion temperature to that being comparable to monimolimnion. If mixolimnion
is deeper and warmer, or when monimolimnion is much cooler, it will need more sustained
and stronger negative heat flux to cool down the lake’s surface and induce whole water column
mixing.
Interaction between meteorological variations and lake history create a negative feedback loop
(Fig. 2.18). Full water column mixing reduces the temperature both at the lake’s surface and
51
bottom. The overall heat content and Schmidt stability decrease. After the event, bottom water
continues to stay cold, while the surface temperature increases. The surface heat flux plays an
important role in regulating the surface temperature so that it reaches an equilibrium dictated
by surface air temperature and other weather conditions. The reestablished surface temperature
accompanied by cooled bottom water leads to an increase in surface and bottom temperature
difference. The greater difference leads to an increase in Schmidt Stability and makes the
lake more resistant to mixing. Thus, even with similar magnitudes of climate variability, Lake
Towuti is less likely to be mixed again after 2008, when a major mixing and cooling happens.
The strong stratification, however, will not prohibit heat exchange between the surface and the
bottom part of the lake. This process is accomplished by heat diffusion (≈ 60 m/yr). This
divisional process warms up lake bottom temperature, and reduces the overall stratification
strength gradually. Overtime, the lake goes back to its original state where surface and bottom
temperature difference and lake stability are weak. The overall process resembles a negative
feedback cycle.
Figure 2.18: Conceptual diagram of the mixing-stratification dynamics in Lake Towuti. Starting
from a weak thermal gradient, when there is a strong negative surface heat flux, surface water
cools down to a point where the whole water column mixes. If the negative heat flux is
strong, the whole lake cools down by a significant amount. After this event, the lake gradually
reestablishes its surface temperature, while the bottom water stays cold. Lake stability after
restratification gets stronger, and thus the lake is less vulnerable to mixing. With enough time,
the gradual diffusional process will blend out this strong stratification, and make it weak again,
to its original state.
52
2.3.3 Lake Towuti During the Last Glacial Maximum
The cooler and drier condition during the LGM affects the stratification of Lake Towuti by
reducing the lake level, and reducing the surface heat flux. Hydrodynamic simulation applying
LGM climate and various water depths indicates that the lake level of Towuti being 150 meter
or less can lead to frequent whole water column mixing. This ties to the oxidized water column
conditions proposed by Russell et al. (2020).
An interesting phenomenon is observed with the LGM simulations. Varying lake level did not
result in significant difference in lake temperature structure under modern climate conditions.
However, under LGM climate conditions, the lowering of lake level increases thermocline
depth, and increases overall heat content. This is because shallower lake during LGM had
higher turbulent diffusivity, significantly increasing the rate of heat diffusion, allowing faster
penetration of surface heat towards the bottom. The overall surface heat flux will adjust to this
increase, and leads to the increase in overall heat content.
2.3.4 The Impact of Global Climate Changes
The effect of climate on stratification of tropical lakes is often invoked in discussions but rarely
quantified in mechanistic terms. In the context of climate warming, stratification of lakes
is believed to strengthen globally (Woolway and Merchant, 2019) as warming of the surface
increases the temperature difference from the deeper waters.
The global average air temperature has increased, on average, by 0.6 °C over the 20th century
(Albritton et al., 2001). The effect of the warming trend has been seen on global lakes. The
increase in surface water temperature on Lake Kivu has been estimated to be 0.12 °C per decade
(Katsev et al., 2014). The warming rate of Lake Superior summer temperature is 0.35 °C per
decade Austin and Colman (2008). From two-decades of simulations conducted in Lake Towuti,
however, the warming trend is not yet significant and it is difficult to be separated from the rest
53
of yearly variabilities.
From previous discussions, air temperature is recognized as the major control on surface lake
temperature and lake stability. Based on this conclusion, it is foreseeable that Lake Towuti, and
similar tropical lakes will experience increasing lake surface temperature and increase stability
under warming air temperature. The negative feedback proposed in Fig. 2.18 might still be in
effect if the rate of increasing air temperature does not exceed the rate of restoring temperature
gradient via heat diffusion. This feedback, however, will not prevent the warming of a lake’s
overall temperature. Detailed investigation is still needed to distinguish the patterns of climate
change during dry season versus during wet season, as the timing of the changes is important in
dictating if a lake will be fully mixed.
2.3.5 Correlation with ENSO
Located in the center of the Indo-Pacific warm pool (IPWP), Sulawesi Island is expected to be
influenced by decadal climate variations, including El Nin˜o–Southern Oscillation (ENSO). In
the tropical western Pacific region, El Nin˜o years often correspond to a drier climate with lower
precipitation (Xie et al., 2016). The common conception is that, the cooler and drier climate
conditions during El Nin˜o years will favor the mixing of the lake.
The correlation, however, is not typically strong concluded from the simulation. The mixing
event in 2002-2003 coincided with an El Nin˜o year while the 2008-2009 mixing event coincided
with a La Nin˜a year (Fig. 2.11). This could be because the correlation between ENSO index
and Indonesia local climate is not straightforward as well. Most climate studies have focused
on the effect of ENSO events on rainfall over Indonesia, and obtained somewhat satisfactory
correlations with regional variations (As-syakur et al., 2014, Supari et al., 2018, Firmansyah
et al., 2022). However, there are limited investigations on air temperature and wind speed. The
correlation between ENSO MEI.v2 index and ERA5 reanalysis meteorology over Lake Towuti
54
is looked at in more detail in appendix (Fig. D.4). Generally, the correlation is weak.
This weak correlation revealed the complexity of the ENSO system, since it was developed to
characterize the global multi-decadal climate cycle. Investigations on sea surface temperature
with the regional climate system, such as the Indo-Pacific Warm Pool (IPWP), might yield
stronger correlations and provide relevant knowledge about the lake stability under the influence
of climatic variabilities.
Besides the weak correlation between the ENSO index and local climate, the lack of
correspondence could be due to the fact that the stratification is also affected by the history of
the lake. Stability and heat content from previous seasons and years contribute to how the lake
responds to immediate weather conditions.
2.3.6 What Makes a Tropical Lake Meromictic
From the inspection on Lake Towuti dynamics, a general conclusion is given for tropical lakes.
Tropical lakes have very different climate cycles compared to mid or high latitude lakes. The
seasonality over the tropics is often classified into dry-wet seasons. However, compared to
mid-latitude lakes, the water temperature in tropical lakes rarely goes through drastic changes
annually. This provides the basis for tropical lakes to stay stably stratified.
For thermogenic tropical lakes, their stratification is mainly maintained by an often weak
temperature gradient between the surface and bottom water. They lack other stability
mechanisms since the salinity content is not pronounced. Thus, these lakes can be heavily
impacted by the changes of the surrounding climate. An episode of strong, negative heat
flux, induced by cooler and drier climate can disrupt this thermal stratification. The effect
of heat flux is related to the lake’s previous history. When the lake originally had a cold
monimolimnion, or a warmer and thicker mixolimnion, it takes greater effort for these negative
surface heat flux to cool down the lake’s surface and induce mixing.
55
Other properties of these lakes can affect how they respond to the changes in climate and surface
heat flux. A lake’s morphology plays a vital role. Deeper equatorial lakes with smaller surface
area protects the lake from variations in surface climate conditions. The maximum depth is
especially important (confirming the results observed from Kraemer et al. (2015)). As Katsev
et al. (2017) point out, some typical tropical thermogenic lake (Lake Malawi, Tanganyika, and
Matano) have seasonal thermocline depth at 30-40 m during wet season, and 80-100 m during
dry season. Based on this estimation, one can expect that tropical lakes that are less than 100
m deep are likely to be mixed annually. Tropical lakes that are 100 - 200 m deep, such as Lake
Towuti, are more stable, but still likely to be mixed occasionally. Deeper tropical lakes, such
as Lake Tanganyika, Malawi, and Matano, have greater probability of staying meromictic for
centuries. Lake Atitlan and Lake Toba, for example, for which long term monitoring data are
lacking, are likely to stay meromictic because their maximum depths are greater than 200 m.
The simple categorization based on lake depth is not definitive. Variations can be induced by
several other factors as well. Firstly, it depends on the details of the tropical climate where
the lake is based. Tropical lakes at higher latitude have stronger seasonal cycles that amplifies
the effect of negative surface heat flux during dry seasons. Additionally, lake trophic status
and light attenuation rate strongly control thermocline depth based on their effect on absorption
of solar radiation. Heiskanen et al. (2015) has pointed out the importance of water clarity on
their response to climate change. Eutrophic lakes often have stronger abilities to scatter light,
resulting in larger lake attenuation coefficient. Thus, the solar radiation absorption throughout
the water column will attenuate faster, leading to warmer and shallower epilimnion. These lakes
will be less susceptible to mixing.
56
2.3.7 Performance of the Model, Caveats of Interpretations
Overall, Simstrat forced by ERA5 reanalysis hourly climate matches the observational
temperature range and stratification patterns. The simulated temperature profiles are within 1
°C difference from all the observational profiles. Simulated profiles also match the pattern of
temporal progression of discrete CTD prior to 2009 and continuous mooring between 2012 and
2015.
However, there are few limitations from the hydrodynamic simulations in this project. (1)
Simstrat simulates mixing processes by imposing larger turbulence diffusion terms (ν terms,
related to turbulent kinetic energy and dissipation, proportional to
k2
ε
), rather than including
the physical advection process. (2) Simstrat assumes 1-D physics. Depending on where the
CTD and mooring data are measured in Lake Towuti, the observed temperature may indicate
a slightly different physical setting than what was modeled. The limitation of 1-D model also
includes the empirical simulations of seiche activities, which could be better modeled with
two or three-dimensional models. This was not helped by the fact that these empirical seiche
parameters are not calibrated in Lake Towuti due to its limited mooring data. (3) There are
no direct measurements of meteorology at the lake’s surface. Thus, the climate forcing inputs
are only available with reanalysis datasets. Available weather stations are far away from Lake
Towuti and are at lower altitudes. They cannot provide confident verification. (4) There are
other limitations in input data. Such as, there is no information on hydrological inputs, and no
information on light attenuation variability. Thus, they are assumed to be reasonable constants.
This might result in simulations deviating from real life values.
57
2.4 Conclusions
Using hydrodynamic model simulation, it is revealed that Lake Towuti has possibly been
completely mixed at least once in the 2000s. The mixing is induced largely by cooler air
during the dry season, when the lake has weak stability. The major mixing event cools down
the bottom water temperature, resulting in an enlarged temperature difference once the lake
restratifies again during wet season. It increases lake stability thereafter. This process can be
described by a negative feedback loop (Fig. 2.18).
Since the maximum seasonal thermocline depth for global major thermogenic meromictic
tropical lakes is around 80 m, we concluded a bathymetry division that determines if a tropical
lake (without salinity gradient) is going to be meromictic. The criteria is:
Maximum Depth Mixing status
< 100 m warm monomictic
100 m - 200 m oligomictic
> 200 m meromictic
Table 2.4: Concluded mixing criteria for tropical lakes without salinity gradient
Of course, this criteria fluctuates with the climate conditions and trophic levels. When the
climate is cooler, drier, and windier, the lakes are more susceptible to mixing. When a lake is
more eutrophic, resulting in faster light attenuation and shallower thermocline, it is less prone
to mixing.
The climate during the LGM was colder and drier. Thus, tropical lakes during that time were
prone to mixing. The cooler and drier conditions also resulted in reduced lake levels. For Lake
Towuti, we suggest it would frequently mix during the LGM if it were less than 150 meter deep,
setting the environmental conditions for frequent oxygenation.
Chapter 3
Biogeochemical Model
3.1 Introduction
This chapter explores the chemical dynamics in Lake Towuti, in particular the effects of
changes in physical stratification on the distributions of redox-sensitive chemical species:
oxygen, sulfate/sulfide, iron(II)/(III), and methane. It approaches this problem with the help of
a reaction-transport biogeochemical model.
The model simulates the rates at which the chemicals accumulate or become depleted in the
lake’s monimolimnion. The results should inform us about the timescales of stratification
or mixing that are required for the observed (Fig. 3.3 dot profiles, provided by Dr. Sean
Crowe) chemical concentrations. When stratification persists for a long time, oxygen eventually
becomes depleted. Anoxic environments favor reduced species (H2S, Fe
2+, and CH4). Mixing
of the entire water-column, in contrast, reintroduces oxygen into the deep waters and can
reset the concentrations of these reduced species. An agreement between the observations and
the distributions of O2, H2S and Fe2+ simulated under the conditions of a persistent stable
stratification would imply that no mixing happened in the recent past. Otherwise, deviations
58
59
in the concentrations of reduced species might indicate a possible non-steady state, e.g. after
an episodic mixing. Numerical simulations may then provide estimates for the time scales
involved.
To investigate whether the chemical dynamics in Lake Towuti may be used as an independent
way of verifying the findings from the previous chapter about the possibility of episodic mixing,
we are seeking answers to the following questions:
1. How long does it take to deplete oxygen in the lake bottom waters when the lake becomes
stably stratified after a complete mixing?
2. How quickly can oxygen be replenished, and to what concentrations, during the episodes
of weak mixing, such as during the potential seasonal events that are predicted by our
hydrodynamic model?
3. How do the dynamics of stratification and mixing affect the redox chemistry in the
monimolimnion? In particular, how do the concentrations of dissolved ferrous iron and
dissolved hydrogen sulfide vary with time after the oxygenation events, and following the
restoration of anoxic conditions?
3.2 Methodology
3.2.1 Chemistry Model
We used a one-dimensional reaction-transport biogeochemical model. The model considers
the kinetics of biogeochemical reactions together with the physical transport of the reactants
and products by turbulent diffusion. It features an explicit treatment of microbial biomasses
that catalyze the anaerobic processes of iron reduction, sulfate reduction, and methanogenesis
through the respective microbial metabolisms. The model was adopted from Sergei Katsev and
Itay Halevy (PNAS in prep) and its description is given in Appendix C. The model solves a
60
system of partial differential equations for the set of included chemical species. The dominant
equation describing the dynamics of chemical species is:
∂C j
∂ t
=
∂
∂ z
(
Kz(z)
∂C j
∂ z
)
− v j ∂C j∂ z +∑i
νiRi j (3.1)
Here, z is the depth below the lake surface, Kz is the depth-dependent vertical diffusion
coefficient, C j is the concentration of a chemical species, v j is the advection (or settling, for
solid particles) velocity, and Ri j are the rates of reactions related to the chemical species with
stoichiometric coefficients νi. The reactions involved in this model describe the transformations
of oxygen, carbon, sulfur, and iron (Table C.2). Table C.3 and C.4 describe the corresponding
kinetic formulations and parameter values. The model was supplied with parameters that were
the best estimates for the conditions in Lake Towuti.
3.2.2 Using the Kz from the Hydrodynamic Model
To describe the intensity of physical transport within the stratified water column, the chemical
model uses the Kz profile that is approximated based on the results from the hydrodynamic
model described in the previous chapter. As described above, the k−ε model Simstrat simulates
mixing as a diffusive process using molecular and turbulent eddy diffusivity. Advection or
convection of water, which are rare under stratified conditions, are not explicitly modeled
but rather approximated by roughly equivalent increases in eddy diffusivity. The turbulent
diffusivity of temperature from the model output is thought to be a good representation for
the transport (Kz) of chemical species in the chemistry model, as both temperature and solutes
are distributed by turbulent eddies, at rates typically above the molecular diffusion of heat and
well above the molecular diffusion of solutes. The v j term in Eq. 3.1 only applies to solid
chemical species, describing their downward settling.
61
The temporal evolution of Kz in Lake Towuti simulated by the hydrodynamic model is shown
in Fig. 3.1. Imposing this fully time-resolved evolution of Kz on the chemical model has
proved to be computationally challenging, however. Rather, to answer questions about the time
scales of the chemical dynamics formulated above, we considered two contrasting (simpler)
scenarios for the temporal variations in Kz. One (the ”stable stratification” scenario) describes
a typical, time-averaged, Kz profile during a prolonged stratification period. The values for this
scenario were extracted from the 2003 dry and wet seasons. The other (the ”mixing” scenario)
uses the Kz profile that characterizes the brief period of enhanced mixing, such as during the
2008 hypothesized mixing event. The values of Kz for this scenario were extracted from the
2008 dry season. As the sharp variations in Kz by several orders of magnitude across a narrow
thermocline generate excessive stiffness in the MATLAB pdepe solver, we further simplified
the profile in chemical simulations by using its analytical approximation, as illustrated in Fig.
3.2. The analytical profile was constructed using Eq.3.2:
Figure 3.1: Turbulent diffusivity of temperature, equivalent to Kz,
simulated by the hydrodynamic model discussed in Chapter 2
62
Kz(z) = K0z
(
fKz+
1− fKz
1+ e(z−HK)/hK
)
+3(
z
60
)4 (3.2)
Here K0z is the surface diffusivity, fKz is the fraction of the bottom diffusivity in relation to
the surface. Hk is the width of thermocline, and hk is the thermocline depth for mixing. The
parameters used to simulate a typical stratify or mixing Kz are listed in Table 3.1. The results
are visualized in Fig. 3.2.
stratify mix scenario 1 mix scenario 2 units
K0z 3×10−4 3×10−3 1×10−2 m2 s−1
fKz 1×10−3 1×10−1 1×10−1 m2 s−1
Hk 8 8 8 m
hk 40 40 40 m
Table 3.1: Parameters defining analytical Kz profiles
The analytical Kz was approximated based on the typical variations, tailored to the resolution
of the PDEs in MATLAB. While the analytical profile for the ”steady state” scenario has
a more diffused thermocline than the profiles obtained in the simulation (Fig. 3.2), the
wider ”effective” width of the thermocline does reflect the vertical range within which the
thermocline migrates seasonally in Lake Towuti. Key to the downward transport of oxygen
and the accumulation of reduced species is the minimum value of Kz within the thermocline,
as it is that minimum value that limits the vertical transport of substances between the surface
and deep waters. The minimum value in the analytical profile is somewhat higher than minima
in the individual simulated profiles (Fig. 3.2), but reflects the potential influences of the
brief periods of incomplete mixing during the dry season, when Kz in the thermocline region
increases by up to two orders of magnitude (Fig. 3.1). The lower values of Kz in the epilimnion
compared to the simulated values do not have a significant effect, as mixing in the epilimnion
is still vigorous enough to effectively transport substances across the surface mixed layer.
63
Figure 3.2: Typical stratification (left) and mixing (right) Kz profiles extracted from
the hydrodynamic model (square and circular markers) and their corresponding analytical
approximation (solid line) used in the biogeochemical model.
3.2.3 Simulation Procedure
The Kz profiles described above were treated within individual simulation runs as constant
profiles, without temporal dynamics. In order to mimic the dynamic changes in Kz, separate
simulation runs were performed, based on the following procedures:
1. Obtaining long-term stratification with bottom-water anoxia: Starting from oxic initial
conditions, the ”stable stratification” Kz profile was applied for 40 years, which was
sufficient for the bottom water to become anoxic for ten years and for reduced chemical
species to accumulate.
2. Simulating brief periods of enhanced mixing: Starting from anoxic initial conditions, the
mixing Kz profile was applied and the simulation was run for 2 months (which was the
characteristic mixing duration observed from the hydrodynamic model).
64
3. Simulating recovery of stratification after a brief mixing event: Starting from profiles
after the mixing event, the stratified Kz profile was applied again for seven years.
The obtained temporal evolution of the chemical profiles, recorded at regular time periods, were
then inspected, looking for the characteristic time scales over which the reduced chemicals
(particularly H2S and Fe2+) reached concentrations similar to their observed magnitudes.
65
3.3 Results
Starting from fully oxygenated condition, oxygen becomes depleted below 120 m after about
30 years of steady stratification (Fig. 3.3), causing accumulation of H2S, Fe2+ and CH4. Sulfate
(SO4) reduction decreases sulfate concentration below the oxycline but does not deplete it
completely, in agreement with observations and in contrast to Lake Matano (Crowe et al., 2011).
Sulfate is reduced into hydrogen sulfide, which could be oxicized by Fe(III) and precipitate out
as elemental sulfur. This process reduces overall sulfate + sulfide concentration within the water
column, especially after year 40, when highly reactive Fe(III) being produced by oxidizing
Fe2+ at oxycline reacts with sulfide and remove it from the water column. Hydrogen sulfide
accumulates to low micro molar levels, also similar to observations. The shape of the sulfide
profile achieved after 40 years of simulation time differs somewhat from the observations, with
the simulated profile accumulating higher concentration of sulfide around the sulfate reduction
zone. The H2S concentration reaches its maximum levels within a few years of the onset of
anoxia, whereas Fe2+ and CH4 continue to build up. Simulation had numerical problems to
continue after 40 years. However, extending the stratification beyond 40 years is expected to
further decrease the concentration of H2S and increase the Fe2+ concentrations to significantly
higher levels.
Two months of mixing at the intensities suggested by the hydrodynamic model transport
oxygen to deeper and previously anoxic waters. The resultant concentrations of oxygen in the
monimolimnion, nevertheless, remain low. For the weaker mixing scenario, when the Kz ranges
between 3×10−3 at the surface to 3×10−4 m2 s−1 at the bottom, the mixing season fails to
fully oxygenate the deeper part of the lake, allowing for some presence of reduced species
there (Fig. 3.4a). For the stronger mixing scenario, when the Kz ranges between 1×10−2 at the
surface to 1×10−3 m2 s−1 at the bottom, oxygen increases to tens of µmolL−1 at the bottom
of the lake, which is roughly a quarter of oxygen saturation concentration. This downward
66
Figure 3.3: Development of chemical stratification starting from fully
oxygenated initial condition. A typical ”stratification” Kz was applied over
the course of the entire simulation. Red circles are observational profiles
(Dr. Sean Crowe, pers. comm.)
propagation of oxygen eliminates the reduced species (H2S, Fe2+, and CH4) from the water
column after about one month (Fig. 3.4b).
Resumption of stratification after brief and incomplete circulation restores the anoxia
gradually. Following the weaker mixing scenario (Fig. 3.5a), H2S slowly re-accumulates.
The accumulation rate this time ( 2 µmolL−1 in seven years of anoxia) is slower compared to
when the anoxia was established after a prolonged stratification after fully oxic conditions ( 5
µmolL−1 in 4 years of anoxia, Fig. 3.3). Fe2+ gradually accumulates as well, returning to its
pre-mixing level at a rate comparable to the rate in Fig. 3.3. Following the stronger mixing
scenario, however, oxygen persists in the deeper waters for several years and the reduced
substances do not accumulate (Fig. 3.5b).
67
(a) The ”weak mixing” scenario. Mixing Kz ranges from 3×10−3 to 3×10−4 ms−2
(b) The ”strong mixing” scenario. Similar configuration as (a), but Kz ranges
between 1×10−2 and 1×10−3 ms−2.
Figure 3.4: Chemical profiles during a mixing event. Starting from the initial conditions
generated from the 40 years of stratification (Fig. 3.3), a typical mixing Kz was applied over the
course of the simulation (2 months).
68
(a) Restratifying. Start from conditions after weak mixing (Figure 3.4a)
(b) Similar configuration as (a), but start from strong mixing (Figure 3.4b
Figure 3.5: Chemical profiles during the restratification process. Starting from conditions after
mixing (Fig. 3.4), a typical stratify Kz was applied over the course of the simulation (seven
years)
69
3.4 Discussion
3.4.1 Gentle Mixing Reproduces Observed Reduced Water Column
The simulations suggest that, in Lake Towuti, a fully oxygenated water column needs about
thirty years of stratification to have oxygen depleted below the 120 m depth. Once oxygen is
depleted, hydrogen sulfide and dissolved ferrous iron build up in the anoxic monimolimnion.
When the H2S and Fe2+ concentrations are both tens of µmol l−1, the two aqueous species can
coexist with respect to mineral FeS solubility (Rickard, 2006). Longer simulation could not
be completed, for numerical reasons. However, interpreting from the trend of decreasing H2S
(by reactive Fe(III) oxidation) and continuously increasing Fe2+ (due to recycling of Fe2+ from
sediments), it is suggested that Lake Towuti eventually favors ferruginous conditions.
The relatively stable balance of H2S and Fe2+ in Lake Towuti is disturbed by introducing oxygen
through mixing. The simulations suggest that the rate of oxygen introduction may be slower
than one might envision in a turnover that would be analogous to the regular mixing of temperate
lakes. Under the weak mixing scenario, the bottom part of the lake is not fully oxygenated
within the two months of the mixing period. Even under a stronger mixing scenario where Kz is
on the order of 1×10−3 to 1×10−2 m2 s−1, the bottom water does not bring enough oxygen to
oxidize all the H2S and Fe2+ in a month. The introduction of oxygen via turbulent mixing in the
bottom water, however, generated reactive iron, which oxidized H2S. The produced elemental
sulfur S0 then precipitated out and left the reaction zone. Restratification after this mixing does
not reestablish the state observed when the anoxia was first established. The accumulation of
H2S is now slower, probably due to the suppression from reactive iron.
A ”strong mixing” event, even though it did not introduce oxygen to full saturation, was
sufficient to deplete all the reduced species after one month. A seven-year restratification after
the stronger mixing event was not sufficient to deplete oxygen and re-accumulate the reduced
70
species to the observed concentrations.
These results corroborate the findings from the hydrodynamic model and further challenge the
hypothesis that Lake Towuti is permanently meromictic. If the lake had been stably stratified for
more than 50 years without the introduction of oxygen via turbulent diffusivity, it is foreseeable
that the anoxic condition would build up Fe2+ concentrations to hundreds µmolL−1 level,
similar to the level recorded in the neighboring Lake Matano Crowe et al. (2011). Lake Matano
has Fe2+ concentrations around 140 µmolL−1 in its lower strata, and only 0.3 µmolL−1H2S
within a narrow layer near the chemocline. Compared to Lake Matano, Lake Towuti has Fe2+
and H2S concentrations of tens of µmolL−1. This suggests that the chemistry of Lake Towuti
is affected by periodic weak ventilation events, in contrast to Lake Matano that is permanently
stratified.
3.4.2 Limitations of the Chemical Model
The uncertainties in the model could come from several sources. First, some of the model
inputs and parameters are constrained within relatively wide ranges (Appendix B. For example,
the surface flux of FeIII (F0FeIII), the recycling efficiency of reduced species in the sediment
( frec), and the setting velocity of solid species (vsettle) are not measured directly for Lake
Towuti but estimated, based on their likely realistic ranges. There are also uncertainties for
chemical reaction rates that were established based on reaction kinetics typically observed
in other environments and thermodynamics. Further, the time-independent analytical profiles
for Kz, which differed from the hydrodynamic simulation outputs, could provide additional
uncertainty. This chapter explored the upper and lower limits of Kz and suggested possible
chemical dynamics in the lake. The model essentially simulates a one dimensional water
column without areal extent. It only establishes boundary conditions at the surface and bottom,
without considering the chemical exchanges that could happen at horizontal boundaries.
71
Adding some essential chemistry components can further support the interpretation of model
outputs, but this was not achieved in this study. If we have more information about conservative
irons (e.g. [Cl– ]), they can be used to calibrate or verify the transportation processes simulated
by the chemistry model. Additionally, if we can add simulations about nutrient concentrations
(e.g.[PO43 – ], [NH4+]), the model will provide more information regarding the effect of physical
stratification on nutrient cycling and primary production.
3.5 Conclusion
The observed profiles of O2, and µmolL−1 levels of H2S and Fe2+ coexisting in Lake Towuti
can be reproduced by the model if stratification and anoxia persist for up to a decade. Short and
weak mixing events could have happened between several months and several years prior to
the measurements. The chemical profiles that were mostly taken in 2013 are thus consistent
with the possibility of weak circulation and low-level oxygenation as a result of the 2008
weather anomaly. Anoxia longer than a decade would lead to further consumption of H2S
and accumulation of Fe2+ to levels above those observed, potentially to ≈ 100 uM as in Lake
Matano.
Chapter 4
Conclusion
This project simulates the physics and chemistry of Lake Towuti with hydrodynamic and
biogeochemical models. The results expand our knowledge of tropical lakes’ mixing and
stratification and the resultant chemical dynamics.
The hydrodynamic simulations suggest that, for thermogenically stratified tropical lakes, their
stratification status is determined by local climate conditions and lake morphology. Among the
climate variables, air temperature has the most significance. Air temperature strongly controls
the surface heat flux, determining whether the surface of the lake cools down sufficiently to
initiate thermal convection. Wind speed is the second most important factor. It affects both the
heat flux and the momentum flux in a lake. A change in surface heat flux that eliminates the
established thermal gradient (i.e. when lake surface cools down to bottom water temperature)
causes the whole water column to mix through thermal convection, but the process proceeds
slower than overturns in temperate lakes. The bathymetry of a tropical lake also influences its
mixing status. For lakes in the same climate, such as neighboring Lakes Towuti and Matano,
mixing status is mainly determined by their depths. For stable stratification, the depth of a lake
needs to significantly exceed the typical thermocline depth. Large thermogenically stratified
72
73
tropical lakes tend to have stable thermoclines at depths around 80 m (Katsev et al., 2017).
Thus, stratification in lakes shallower than 100 m is susceptible to annual variations in weather,
and may be mixed frequently. Lakes between 100 m to 200 m deep are likely to be oligomictic.
Lakes deeper than 200 m are likely to be meromictic. Changes in lake depth in response to
changes in climate thus can lead to changes in stratification regime. In particular, for Lake
Towuti, the cooler, drier, and windier conditions during the LGM, which likely decreased the
water level to around 150 m, must have induced frequent mixing and oxygenation of the bottom
waters. This agrees well with paleolimnological evidence (Vogel et al., 2015, Goudge et al.,
2017).
Lakes’ mixing regimes regulate their chemical properties. In lakes that sustain stratification
for many years, oxygen becomes consumed in the deeper layers, even under oligotrophic
conditions. For Lake Towuti, the transition from full oxygenation (with oxygen in deeper
waters approaching saturation levels) to deep-water anoxia is suggested by the simulations to
take about 30 years. The anoxia allows the accumulation of reduced species, such as hydrogen
sulfide, ferrous iron, and methane. Abundances of sulfur vs. iron minerals in the lake’s bedrock
and catchment determine whether sulfide or ferrous iron would dominate. For Malili Lake
system, the catchment bedrock is dominated by ultramafic species who provide abundant
iron sources into the lake, suggesting the tendency of Lake Matano and Lake Towuti being
ferruginous lakes.It remains an interesting question whether the balance between the iron-rich
conditions (ferruginous) and the sulfide-rich conditions (euxinic) is unstable, i.e. whether
accumulation of either substance tends to further promote its accumulation at the expense of
another. The abundance of iron suppresses the production of sulfide, by precipitating it out
with Fe2+ as FeS, by oxidizing it with reactive FeIII and precipitating it out as elemental sulfur
(S0), or by encouraging the abundance of iron reducers that compete with sulfate reducers for
the organic metabolic substrates. Our simulations suggest that coexistence of Fe2+ and H2S is
74
possible at low (µmolL−1) concentrations, as observed in Lake Towuti. This balance between
Fe2+ and H2S is made possible by episodic weak ventilation of the deep waters, which limits
the concentration of ferrous iron and resets the reduced species concentrations to low levels.
The simulation results of the chemistry model suggest that it takes at least one-month of
persisted mixing to reintroduce oxygen into the bottom layer of Lake Towuti, or similar
stratified tropical lakes. This rate is slower than what is typically observed in mid-latitude
lakes, where the convection process transports oxygen at a much quicker rate. Thus, it is
possible that lake Towuti might experience a periodical mixing event that homogenizes water
temperature without completely ventilating the water bottom. Although anoxia conditions are
observed from 2013 to 2015, it is possible for Lake Towuti to experience full water column
mixing immediately before this.
The length of stratification and anoxia have important implications on the primary productivity
of a lake. On one side, stratification can limit the availability of nutrients at the lake surface and
limit plankton’s primary production. On the other hand, the stratification can lead to oxygen
depletion and results in reduced environment at the lake bottom, releasing the phosphorus
trapped by iron oxides at the lake sediments, introducing nutrients to the water column to
enhance the water column primary production. The net effect is related to the initial phosphorus
concentration in the water column and sediment. Thus the variations in climate conditions and
changes in lake stability have profound influence on its primary production.
Chapter 5
Bibliography
Ageli, M. K., Hamilton, P. B., Bramburger, A. J., Weidman, R. P., Song, Z., Russell, J., Vogel,
H., Bijaksana, S., and Haffner, G. D. Benthic-Pelagic state changes in the primary trophic
level of an ancient tropical lake. Palaeogeography, Palaeoclimatology, Palaeoecology, 594
(March):110937, 2022. ISSN 00310182. doi: 10.1016/j.palaeo.2022.110937.
Ajayi, T. and Gupta, I. A review of reactive transport modeling in wellbore integrity
problems. Journal of Petroleum Science and Engineering, 175(January):785–803, 2019.
ISSN 09204105. doi: 10.1016/j.petrol.2018.12.079.
Albritton, D. L., Filho, L. G. M., Cubasch, U., Dai, X., Ding, Y., Griggs, D. J., Hewitson, B.,
Isaksen, I., Karl, T., McFarland, M., Meleshko, V. P., Mitchell, J. F. B., Noguer, M., Nyenzi,
B. S., Oppenheimer, M., Penner, J. E., Pollonais, S., Stocker, T., Trenberth, K. E., Allen,
M. R., Baede, A. P. M., Church, J. A., Ehhalt, D. H., Folland, C. K., Giorgi, F., Gregory,
J. M., Haywood, J. M., House, J. I., Hulme, M., Jaramillo, V. J., Jayaraman, A., Johnson,
C. A., Joussaume, S., Karoly, D. J., Kheshgi, H., Que´re´, C. L., Mata, L. J., McAvaney,
B. J., Mearns, L. O., Meehl, G. A., Iii, B. M., Mugara, R. K., Prather, M., Prentice, C.,
Ramaswamy, V., Raper, S. C. B., Salinger, M. J., Wang, M.-X., Watson, R. T., Yap, K.-
S., Joos, F., Ramirez-Rojas, A., Stone, J. M. R., and Zillman, J. 2001: Observed Climate
Variability and Change. In Climate Change 2001: The Scientific Basis. Contribution of
Working Group I to the Third Assessment Report of the Intergovernmental Panel on Climate
Change. Cambridge University Press, Cambridge, United Kingdom, 2001.
Alcocer, J. Mexican Meromictic Lakes: What We Know So Far. In Gulati, R. D., Zadereev,
E. S., and Degermendzhi, A. G., editors, Ecology of Meromictic Lakes, volume 228, pages
353–375. Springer International Publishing, Cham, 2017. ISBN 978-3-319-49141-7 978-3-
319-49143-1. doi: 10.1007/978-3-319-49143-1 12.
Alin, S. R. and Johnson, T. C. Carbon cycling in large lakes of the world: A synthesis of
75
76
production, burial, and lake-atmosphere exchange estimates. Global Biogeochemical Cycles,
21(3):1–12, 2007. ISSN 08866236. doi: 10.1029/2006GB002881.
Anderson, R. Y. Meromictic Lakes and Varved Lake Sediments in North America. U.S.
Government Printing Office, 1985.
As-syakur, A. R., Adnyana, I. W. S., Mahendra, M. S., Arthana, I. W., Merit, I. N., Kasa, I. W.,
Ekayanti, N. W., Nuarsa, I. W., and Sunarta, I. N. Observation of spatial patterns on the
rainfall response to ENSO and IOD over Indonesia using TRMM Multisatellite Precipitation
Analysis (TMPA). International Journal of Climatology, 34(15):3825–3839, 2014. ISSN
1097-0088. doi: 10.1002/joc.3939.
Austin, J. and Colman, S. A century of temperature variability in Lake Superior. Limnology and
Oceanography, 53(6):2724–2730, 2008. ISSN 1939-5590. doi: 10.4319/lo.2008.53.6.2724.
Ba¨renbold, F., Kipfer, R., and Schmid, M. Dynamic modelling provides new insights into
development and maintenance of Lake Kivu’s density stratification. Environmental Modelling
& Software, 147:105251, January 2022. ISSN 1364-8152. doi: 10.1016/j.envsoft.2021.
105251.
Berner, R. A. A New Geochemical Classification of Sedimentary Environments. SEPM
Journal of Sedimentary Research, Vol. 51, 1981. ISSN 1527-1404. doi: 10.1306/
212F7C7F-2B24-11D7-8648000102C1865D.
Boehrer, B., von Rohden, C., and Schultze, M. Physical Features of Meromictic Lakes:
Stratification and Circulation. In Ecology of Meromictic Lakes, pages 15–34. Springer, Cham,
2017. doi: 10.1007/978-3-319-49143-1 2.
Canfield, D. E. and Thamdrup, B. Towards a consistent classification scheme for geochemical
environments, or, why we wish the term ’suboxic’ would go away: Editorial. Geobiology, 7
(4):385–392, 2009. ISSN 14724677. doi: 10.1111/j.1472-4669.2009.00214.x.
Chen, C.-T. A. and Millero, F. J. Thermodynamic Properties for Natural Waters Covering
Only the Limnological Range. Limnology and Oceanography, 31(3):657–662, 1986. ISSN
19395590. doi: 10.4319/lo.1986.31.3.0657.
Chitamwebwa, D. B. R. Meromixis, stratification and internal waves in Kigoma waters of
Lake Tanganyika. In Lindqvist, O. V., Mo¨lsa¨, H., Salonen, K., and Sarvala, J., editors, From
Limnology to Fisheries: Lake Tanganyika and Other Large Lakes, pages 59–64. Springer
Netherlands, Dordrecht, 1999. ISBN 978-90-481-5339-8 978-94-017-1622-2. doi: 10.1007/
978-94-017-1622-2 5.
Costa, K. M., Russell, J. M., Vogel, H., and Bijaksana, S. Hydrological connectivity and mixing
of Lake Towuti, Indonesia in response to paleoclimatic changes over the last 60,000 years.
Palaeogeography, Palaeoclimatology, Palaeoecology, 417:467–475, 2015. ISSN 00310182.
doi: 10.1016/j.palaeo.2014.10.009.
77
Crowe, S. A., Katsev, S., Leslie, K., Sturm, A., Magen, C., Nomosatryo, S., Pack, M. A.,
Kessler, J. D., Reeburgh, W. S., Roberts, J. A., Gonza´lez, L., Douglas Haffner, G., Mucci, A.,
Sundby, B., and Fowle, D. A. The methane cycle in ferruginous Lake Matano. Geobiology,
9(1):61–78, January 2011. ISSN 14724677. doi: 10.1111/j.1472-4669.2010.00257.x.
Damanik, A., Wille, M., Ahmad, Q., Russell, J. M., Bijaksana, S., Cahyarini, S. Y., and Vogel,
H. Mo isotope variability records climate-driven water column redox changes in ferruginous
Lake Towuti , Indonesia over the past ˜ 30 kyrs. In EGU General Assembly 2022, Vienna,
Austria, 2022.
De Crop, W. and Verschuren, D. Mixing regimes in the equatorial crater lakes of western
Uganda. Limnologica, 90:125891, September 2021. ISSN 0075-9511. doi: 10.1016/j.limno.
2021.125891.
Engesgaard, P. and Christensen, T. H. A review of chemical solute transport models. Nordic
Hydrology, 19(3):183–216, 1988. ISSN 00291277. doi: 10.2166/nh.1988.0013.
Fang, X. and Stefan, H. G. Long-term lake water temperature and ice cover
simulations/measurements. Cold Regions Science and Technology, 24(3):289–304, July
1996. ISSN 0165-232X. doi: 10.1016/0165-232X(95)00019-8.
Firmansyah, A. J., Nurjani, E., and Sekaranom, A. B. Effects of the El Nin˜o-Southern
Oscillation (ENSO) on rainfall anomalies in Central Java, Indonesia. Arabian Journal
of Geosciences, 15(24):1746, December 2022. ISSN 1866-7538. doi: 10.1007/
s12517-022-11016-2.
Frey, D. G. A limnological reconnaissance of Lake Lanao. SIL Proceedings, 1922-2010, 17(2):
1090–1102, November 1969. ISSN 0368-0770. doi: 10.1080/03680770.1968.11895955.
Froelich, P. N., Klinkhammer, G. P., Bender, M. L., Luedtke, N. A., Heath, G. R., Cullen,
D., Dauphin, P., Hammond, D., Hartman, B., and Maynard, V. Early oxidation of
organic matter in pelagic sediments of the eastern equatorial Atlantic: Suboxic diagenesis.
Geochimica et Cosmochimica Acta, 43(7):1075–1090, July 1979. ISSN 0016-7037. doi:
10.1016/0016-7037(79)90095-4.
Fukushima, T., Matsushita, B., Subehi, L., Setiawan, F., and Wibowo, H. Will hypolimnetic
waters become anoxic in all deep tropical lakes? Scientific Reports, 7(1):45320, May 2017.
ISSN 2045-2322. doi: 10.1038/srep45320.
Fukushima, T., Setiawan, F., Subehi, L., Jiang, D., and Matsushita, B. Water temperature and
some water quality in Lake Toba, a tropical volcanic lake. Limnology, August 2022. ISSN
1439-8621, 1439-863X. doi: 10.1007/s10201-022-00703-4.
Gaudard, A., Ra˚man Vinna˚, L., Ba¨renbold, F., Schmid, M., and Bouffard, D. Toward an open
access to high-frequency lake modeling and statistics data for scientists and practitioners – the
78
case of Swiss lakes using Simstrat v2.1. Geosci. Model Dev., 12(9):3955–3974, September
2019. ISSN 1991-9603. doi: 10.5194/gmd-12-3955-2019.
Goudge, T. A., Russell, J. M., Mustard, J. F., Head, J. W., and Bijaksana, S. A 40,000 yr record
of clay mineralogy at Lake Towuti, Indonesia: Paleoclimate reconstruction from reflectance
spectroscopy and perspectives on paleolakes on Mars. Bulletin of the Geological Society of
America, 129(7-8):806–819, 2017. ISSN 19432674. doi: 10.1130/B31569.1.
Goudsmit, G. H., Burchard, H., Peeters, F., and Wu¨est, A. Application of k-ε turbulence models
to enclosed basins: The role of internal seiches. Journal of Geophysical Research C: Oceans,
107(12):23–1, 2002. ISSN 01480227. doi: 10.1029/2001jc000954.
Guilbaud, R., Poulton, S. W., Butterfield, N. J., Zhu, M., and Shields-Zhou, G. A. A global
transition to ferruginous conditions in the early Neoproterozoic oceans. Nature Geoscience,
8(6):466–470, June 2015. ISSN 1752-0908. doi: 10.1038/ngeo2434.
Hambright, K. D., Gophen, M., and Serruya, S. Influence of long-term climatic changes on the
stratification of a subtropical, warm monomictic lake. Limnology and Oceanography, 39(5):
1233–1242, 1994. ISSN 19395590. doi: 10.4319/lo.1994.39.5.1233.
Hasberg, A. K. M., Bijaksana, S., Held, P., Just, J., Melles, M., Morlock, M. A., Opitz, S.,
Russell, J. M., Vogel, H., and Wennrich, V. Modern sedimentation processes in Lake Towuti,
Indonesia, revealed by the composition of surface sediments. Sedimentology, 66(2):675–698,
February 2019. ISSN 00370746. doi: 10.1111/sed.12503.
Heiskanen, J. J., Mammarella, I., Ojala, A., Stepanenko, V., Erkkila¨, K.-M., Miettinen, H.,
Sandstro¨m, H., Eugster, W., Leppa¨ranta, M., Ja¨rvinen, H., Vesala, T., and Nordbo, A.
Effects of water clarity on lake stratification and lake-atmosphere heat exchange. Journal
of Geophysical Research: Atmospheres, 120(15):7412–7428, 2015. ISSN 2169-8996. doi:
10.1002/2014JD022938.
Hersbach, H., Bell, B., Berrisford, P., Hirahara, S., Hora´nyi, A., Mun˜oz-Sabater, J., Nicolas,
J., Peubey, C., Radu, R., Schepers, D., Simmons, A., Soci, C., Abdalla, S., Abellan, X.,
Balsamo, G., Bechtold, P., Biavati, G., Bidlot, J., Bonavita, M., De Chiara, G., Dahlgren,
P., Dee, D., Diamantakis, M., Dragani, R., Flemming, J., Forbes, R., Fuentes, M., Geer, A.,
Haimberger, L., Healy, S., Hogan, R. J., Ho´lm, E., Janiskova´, M., Keeley, S., Laloyaux, P.,
Lopez, P., Lupu, C., Radnoti, G., de Rosnay, P., Rozum, I., Vamborg, F., Villaume, S., and
The´paut, J.-N. The ERA5 global reanalysis. Quarterly Journal of the Royal Meteorological
Society, 146(730):1999–2049, 2020. ISSN 1477-870X. doi: 10.1002/qj.3803.
Hostetler, S. W. Simulation of Lake Evaporation with an Energy Balance-Eddy Diffusion Model
of Lake Temperature: Model Development and Validation, and Application to Lake-Level
Variations at Harney-Malheur Lake, Oregon. PhD thesis, University of Oregon, 1987.
Hutchinson, G. E. and Lo¨ffler, H. The Thermal Classification of Lakes. Proceedings of the
79
National Academy of Sciences, 42(2):84–86, 1956. ISSN 0027-8424. doi: 10.1073/pnas.42.
2.84.
Idso, S. B. On the concept of lake stability. Limnology and Oceanography, 18(4):681–683, July
1973. ISSN 0024-3590. doi: 10.4319/lo.1973.18.4.0681.
Imboden, D. M. and Wu¨est, A. Mixing Mechanisms in Lakes. In Lerman, A., Imboden,
D. M., and Gat, J. R., editors, Physics and Chemistry of Lakes, pages 83–138. Springer
Berlin Heidelberg, Berlin, Heidelberg, 1995. ISBN 978-3-642-85132-2. doi: 10.1007/
978-3-642-85132-2 4.
Jassby, A. and Powell, T. Vertical patterns of eddy diffusion during stratification in Castle Lake,
California1. Limnology and Oceanography, 20(4):530–543, 1975. ISSN 1939-5590. doi:
10.4319/lo.1975.20.4.0530.
Joehnk, K. and Umlauf, L. Modelling the metalimnetic oxygen minimum in a medium sized
alpine lake. Ecological Modelling, 136(1):67–80, January 2001. ISSN 03043800. doi:
10.1016/S0304-3800(00)00381-1.
Katsev, S. and Crowe, S. A. Organic carbon burial efficiencies in sediments: The power law
of mineralization revisited. Geology, 43(7):607–610, 2015. ISSN 19432682. doi: 10.1130/
G36626.1.
Katsev, S., Crowe, S. A., Mucci, A., Sundby, B., Nomosatryo, S., Douglas Haffner, G., and
Fowle, D. A. Mixing and its effects on biogeochemistry in the persistently stratified, deep,
tropical Lake Matano, Indonesia. Limnology and Oceanography, 55(2):763–776, March
2010. ISSN 00243590. doi: 10.4319/lo.2010.55.2.0763.
Katsev, S., Aaberg, A. A., Crowe, S. A., and Hecky, R. E. Recent Warming of Lake Kivu. PLoS
ONE, 9(10):e109084, October 2014. ISSN 1932-6203. doi: 10.1371/journal.pone.0109084.
Katsev, S., Verburg, P., Lliro´s, M., Minor, E. C., Kruger, B. R., and Li, J. Tropical
Meromictic Lakes: Specifics of Meromixis and Case Studies of Lakes Tanganyika, Malawi,
and Matano. In Ecology of Meromictic Lakes, pages 277–323. Springer, Cham, 2017. doi:
10.1007/978-3-319-49143-1 10.
Kling, G. W., Tuttle, M. L., and Evans, W. C. The evolution of thermal structure and water
chemistry in Lake Nyos. Journal of Volcanology and Geothermal Research, 39(2):151–165,
November 1989. ISSN 0377-0273. doi: 10.1016/0377-0273(89)90055-3.
Kraemer, B. M., Anneville, O., Chandra, S., Dix, M., Kuusisto, E., Livingstone, D. M.,
Rimmer, A., Schladow, S. G., Silow, E., Sitoki, L. M., Tamatamah, R., Vadeboncoeur, Y., and
McIntyre, P. B. Morphometry and average temperature affect lake stratification responses to
climate change. Geophysical Research Letters, 42(12):4981–4988, 2015. ISSN 1944-8007.
doi: 10.1002/2015GL064097.
Kranenburg, W., Tiessen, M., Veenstra, J., de Graaff, R., Uittenbogaard, R., Bouffard, D.,
80
Sakindi, G., Umutoni, A., Van de Walle, J., Thiery, W., and van Lipzig, N. 3D-modelling of
Lake Kivu: Horizontal and vertical flow and temperature structure under spatially variable
atmospheric forcing. Journal of Great Lakes Research, 46(4):947–960, 2020. ISSN
03801330. doi: 10.1016/j.jglr.2020.05.012.
Likongwe, P., Kihoro, J., Ngigi, M., and Jamu, D. Modeling spatial non-stationarity of Chambo
in South-East Arm of Lake Malawi. Asian Journal of Applied Science and Engineering, 4
(2):81–90, 2015.
MacIntyre, S. and Melack, J. M. Meromixis in an equatorial African soda lake1. Limnology and
Oceanography, 27(4):595–609, 1982. ISSN 1939-5590. doi: 10.4319/lo.1982.27.4.0595.
Mironov, D. V. Parameterization of lakes in numerical weather prediction: Description of a lake
model. Technical Report 11, German Weather Service, O enbacham Main, Germany, August
2008.
Oleksy, I. A. and Richardson, D. C. Climate Change and Teleconnections Amplify Lake
Stratification With Differential Local Controls of Surface Water Warming and Deep Water
Cooling. Geophysical Research Letters, 48(5):e2020GL090959, 2021. ISSN 1944-8007.
doi: 10.1029/2020GL090959.
Perroud, M., Goyette, S., Martynov, A., Beniston, M., and Anneville, O. Simulation of
multiannual thermal profiles in deep Lake Geneva: A comparison of one-dimensional lake
models. Limnology and Oceanography, 54(5):1574–1594, 2009. ISSN 00243590. doi:
10.4319/lo.2009.54.5.1574.
Poulton, S. W. and Canfield, D. E. Ferruginous conditions: A dominant feature of the ocean
through Earth’s history. Elements, 7(2):107–112, 2011. ISSN 18115217. doi: 10.2113/
gselements.7.2.107.
Poulton, S. W., Fralick, P. W., and Canfield, D. E. The transition to a sulphidic ocean ∼ 1.84
billion years ago. Nature, 431(7005):173–177, September 2004. ISSN 1476-4687. doi:
10.1038/nature02912.
Rejma´nkova´, E., Koma´rek, J., Dix, M., Koma´rkova´, J., and Giro´n, N. Cyanobacterial blooms in
Lake Atitlan, Guatemala. Limnologica, 41(4):296–302, December 2011. ISSN 0075-9511.
doi: 10.1016/j.limno.2010.12.003.
Richardson, D. C., Melles, S. J., Pilla, R. M., Hetherington, A. L., Knoll, L. B., Williamson,
C. E., Kraemer, B. M., Jackson, J. R., Long, E. C., Moore, K., Rudstam, L. G., Rusak,
J. A., Saros, J. E., Sharma, S., Strock, K. E., Weathers, K. C., and Wigdahl-Perry, C. R.
Transparency, Geomorphology and Mixing Regime Explain Variability in Trends in Lake
Temperature and Stratification across Northeastern North America (1975–2014). Water, 9
(6):442, June 2017. ISSN 2073-4441. doi: 10.3390/w9060442.
81
Rickard, D. The solubility of FeS. Geochimica et Cosmochimica Acta, 70(23):5779–5789,
December 2006. ISSN 00167037. doi: 10.1016/j.gca.2006.02.029.
Russell, J. M., Bijaksana, S., Vogel, H., Melles, M., Kallmeyer, J., Ariztegui, D., Crowe,
S., Fajar, S., Hafidz, A., Haffner, D., Hasberg, A., Ivory, S., Kelly, C., King, J., Kirana,
K., Morlock, M., Noren, A., O’Grady, R., Ordonez, L., Stevenson, J., von Rintelen, T.,
Vuillemin, A., Watkinson, I., Wattrus, N., Wicaksono, S., Wonik, T., Bauer, K., Deino, A.,
Friese, A., Henny, C., Imran, Marwoto, R., Ngkoimani, L. O., Nomosatryo, S., Safiuddin,
L. O., Simister, R., and Tamuntuan, G. The Towuti Drilling Project: Paleoenvironments,
biological evolution, and geomicrobiology of a tropical Pacific lake. Scientific Drilling, 21
(July):29–40, 2016. ISSN 18163459. doi: 10.5194/sd-21-29-2016.
Russell, J. M., Vogel, H., Bijaksana, S., Melles, M., Deino, A., Hafidz, A., Haffner, D., Hasberg,
A. K., Morlock, M., von Rintelen, T., Sheppard, R., Stelbrink, B., and Stevenson, J. The
late quaternary tectonic, biogeochemical, and environmental evolution of ferruginous Lake
Towuti, Indonesia. Palaeogeography, Palaeoclimatology, Palaeoecology, 556(August), 2020.
ISSN 00310182. doi: 10.1016/j.palaeo.2020.109905.
Sahoo, G. B., Forrest, A. L., Schladow, S. G., Reuter, J. E., Coats, R., and Dettinger, M. Climate
change impacts on lake thermal dynamics and ecosystem vulnerabilities. Limnology and
Oceanography, 61(2):496–507, 2016. ISSN 19395590. doi: 10.1002/lno.10228.
Shakeri Yekta, S., Svensson, B. H., Bjo¨rn, A., and Skyllberg, U. Thermodynamic modeling of
iron and trace metal solubility and speciation under sulfidic and ferruginous conditions in full
scale continuous stirred tank biogas reactors. Applied Geochemistry, 47:61–73, 2014. ISSN
18729134. doi: 10.1016/j.apgeochem.2014.05.001.
Sigurdsson, H., Devine, J. D., Tchua, F. M., Presser, F. M., Pringle, M. K. W., and Evans,
W. C. Origin of the lethal gas burst from Lake Monoun, Cameroun. Journal of Volcanology
and Geothermal Research, 31(1):1–16, March 1987. ISSN 0377-0273. doi: 10.1016/
0377-0273(87)90002-3.
Skeel, R. D. and Berzins, M. A Method for the Spatial Discretization of Parabolic Equations
in One Space Variable. SIAM Journal on Scientific and Statistical Computing, 11(1):1–32,
January 1990. ISSN 0196-5204, 2168-3417. doi: 10.1137/0911001.
Sommer, T., Schmid, M., and Wu¨est, A. The role of double diffusion for the heat and salt
balance in Lake Kivu. Limnology and Oceanography, 64(2):650–660, 2019. ISSN 19395590.
doi: 10.1002/lno.11066.
Stepanenko, V. and Lykossov, V. Numerical modeling of the heat and moisture transport in
a lake-soil system. Russian Journal for Meteorology and Hydrology, (3):95–104, January
2005.
Stepanenko, V. M., Goyette, S., Martynov, A., Perroud, M., Fang, X., and Mironov, D. First
82
steps of a Lake Model intercomparison project: LakeMIP. Boreal Environment Research, 15
(2):191–202, 2010. ISSN 12396095.
Stetler, J. T., Girdner, S., Mack, J., Winslow, L. A., Leach, T. H., and Rose, K. C. Atmospheric
stilling and warming air temperatures drive long-term changes in lake stratification in a large
oligotrophic lake. Limnology and Oceanography, 66(3):954–964, 2021. ISSN 19395590.
doi: 10.1002/lno.11654.
Subin, Z. M., Riley, W. J., and Mironov, D. An improved lake model for climate simulations:
Model structure, evaluation, and sensitivity analyses in CESM1. Journal of Advances in
Modeling Earth Systems, 4:M02001, February 2012. ISSN 1942-2466. doi: 10.1029/
2011MS000072.
Sulistioadi, Y. B., Tseng, K.-H., Shum, C. K., Hidayat, H., Sumaryono, M., Suhardiman, A.,
Setiawan, F., and Sunarso, S. Satellite radar altimetry for monitoring small rivers and lakes
in Indonesia. Hydrology and Earth System Sciences, 19(1):341–359, January 2015. ISSN
1027-5606. doi: 10.5194/hess-19-341-2015.
Supari, Tangang, F., Salimun, E., Aldrian, E., Sopaheluwakan, A., and Juneng, L. ENSO
modulation of seasonal rainfall and extremes in Indonesia. Climate Dynamics, 51(7):2559–
2580, October 2018. ISSN 1432-0894. doi: 10.1007/s00382-017-4028-8.
Swanner, E. D., Lambrecht, N., Wittkop, C., Harding, C., Katsev, S., Torgeson, J., and Poulton,
S. W. The biogeochemistry of ferruginous lakes and past ferruginous oceans. Earth-Science
Reviews, 211(June):103430, 2020. ISSN 00128252. doi: 10.1016/j.earscirev.2020.103430.
Taillefert, M. and Gaillard, J.-F. Reactive transport modeling of trace elements in the water
column of a stratified lake: Iron cycling and metal scavenging. Journal of Hydrology, 256:
16–34, September 2001. doi: 10.1016/S0022-1694(01)00524-8.
Thiery, W., Stepanenko, V. M., Fang, X., Jo¨hnk, K. D., Li, Z., Martynov, A., Perroud, M., Subin,
Z. M., Darchambeau, F., Mironov, D., and Van Lipzig, N. P. LakeMIP Kivu: Evaluating the
representation of a large, deep tropical lake by a set of one-dimensional lake models. Tellus,
Series A: Dynamic Meteorology and Oceanography, 66(1):1–18, 2014. ISSN 16000870. doi:
10.3402/tellusa.v66.21390.
Vaillant, J. J., Haffner, G. D., and Cristescu, M. E. The ancient lakes of indonesia: Towards
integrated research on speciation. In Integrative and Comparative Biology, volume 51, pages
634–643. Oxford Academic, October 2011. doi: 10.1093/icb/icr101.
Verburg, P. and Hecky, R. E. The physics of the warming of Lake Tanganyika by climate
change. Limnology and Oceanography, 54(6 PART 2):2418–2430, 2009. ISSN 00243590.
doi: 10.4319/lo.2009.54.6 part 2.2418.
Vogel, H., Russell, J. M., Cahyarini, S. Y., Bijaksana, S., Wattrus, N., Rethemeyer, J., and
Melles, M. Depositional modes and lake-level variability at Lake Towuti, Indonesia, during
83
the past ˜29 kyr BP. Journal of Paleolimnology, 54(4):359–377, September 2015. ISSN
15730417. doi: 10.1007/s10933-015-9857-z.
Vuillemin, A., Friese, A., Alawi, M., Henny, C., Nomosatryo, S., Wagner, D., Crowe, S. A.,
and Kallmeyer, J. Geomicrobiological features of ferruginous sediments from lake Towuti,
Indonesia. Frontiers in Microbiology, 7(JUN):1007, 2016. ISSN 1664302X. doi: 10.3389/
fmicb.2016.01007.
Walker, K. F. and Likens, G. E. Meromixis and a reconsidered typology of lake circulation
patterns: With 3 tables in the text. SIL Proceedings, 1922-2010, 19(1):442–458, October
1975. ISSN 0368-0770. doi: 10.1080/03680770.1974.11896084.
Weinberger, S. and Vetter, M. Using the hydrodynamic model DYRESM based on results of a
regional climate model to estimate water temperature changes at Lake Ammersee. Ecological
Modelling, 244:38–48, 2012. ISSN 03043800. doi: 10.1016/j.ecolmodel.2012.06.016.
Woolway, R. I. and Merchant, C. J. Worldwide alteration of lake mixing regimes in response
to climate change. Nature Geoscience, 12(4):271–276, April 2019. ISSN 17520908. doi:
10.1038/s41561-019-0322-x.
Woolway, R. I., Meinson, P., No˜ges, P., Jones, I. D., and Laas, A. Atmospheric stilling leads to
prolonged thermal stratification in a large shallow polymictic lake. Climatic Change, 141(4):
759–773, April 2017. ISSN 1573-1480. doi: 10.1007/s10584-017-1909-0.
Woolway, R. I., Sharma, S., Weyhenmeyer, G. A., Debolskiy, A., Golub, M., Mercado-Bettı´n,
D., Perroud, M., Stepanenko, V., Tan, Z., Grant, L., Ladwig, R., Mesman, J., Moore, T. N.,
Shatwell, T., Vanderkelen, I., Austin, J. A., DeGasperi, C. L., Dokulil, M., La Fuente, S.,
Mackay, E. B., Schladow, S. G., Watanabe, S., Marce´, R., Pierson, D. C., Thiery, W.,
and Jennings, E. Phenological shifts in lake stratification under climate change. Nature
Communications, 12(1):1–11, 2021. doi: 10.1038/s41467-021-22657-4.
Wuest, A., Piepke, G., and Halfman, J. Combined effects of dissolved solids and temperature
on the density stratification of Malawi. In Paleoclimatology of the Eastern African Lakes,
pages 183–202. Gordon and Breach, Toronto, 1996.
Xie, S. P., Kosaka, Y., Du, Y., Hu, K., Chowdary, J. S., and Huang, G. Indo-western
Pacific ocean capacitor and coherent climate anomalies in post-ENSO summer: A review.
Advances in Atmospheric Sciences, 33(4):411–432, 2016. ISSN 02561530. doi: 10.1007/
s00376-015-5192-6.
Appendix A
Glossary and Acronyms
Care has been taken in this thesis to minimize the use of jargon and acronyms, but this cannot
always be achieved. This appendix contains a table of acronyms and their meaning.
A.1 Acronyms
Table A.1: Acronyms
Acronym Meaning
CTD Conductivity, Temperature, and Depth
TD Temperature Difference
TDI Temperature Difference Indicator
HC Heat Content
HCI Heat Content Indicator
SS Schmidt Stability
SSI Schmidt Stability Index
Continued on next page
84
85
Table A.1 – continued from previous page
Acronym Meaning
WTD Weakest Temperature Difference
MBF Maximum Buoyancy Frequency (squared)
LGM Last Glacial Maximum
ENSO El Nin˜o-Southern Oscillation
IPWP Indo-Pacific Warm Pool
PDE Partial Differential Equation
SR Sulfate Reduction
IR Iron Reduction
Appendix B
Hydrodynamic Modeling Parameters
Table B.1: Physical parameters
Parameter / Method Unit Value Comment
Turbulence Model k− ε
Stability Function Quasi-equilibrium
Flux Condition No flux
Use Filtered Wind false
Seiche Normalization integral
Wind Drag Model Wu¨est and Lorke, 2003
Inflow Mode Inflow is density-driven
Pressure Gradients Bottom Friction
Ice Model off
Snow Model Off
Time Step seconds 60
Continued on next page
86
87
Table B.1 – continued from previous page
Parameter / Method Unit Value Comment
Start Day 2001-09-25
End Day 2016-12-31
Latitude ° -2
Air Pressure mbar 1010
αseiche 1 0.0011046 Fraction of wind energy which goes
into seiche energy
QN2 1 0.75 Fit parameter for distribution of seiche
energy
cd 1 0.002 Bottom drag coefficient
Hgeo Wm−2 0 Geothermal heat flux
psw water 1 0.8 Fit parameter for absorption of short
wave radiation from sky
plw 1 0.95 Fit parameter for absorption of
infrared radiation from sky
pwind f 1 1 Fit parameter for convective and latent
heat fluxes
βsol 1 0.3 Fraction of short-wave radiation
directly absorbed as heat by water
Initial Seiche Energy J 1×1010 Initial seiche energy
cµ 1 0.09 constant for calculating turbulent
viscosity for momentum
Continued on next page
88
Table B.1 – continued from previous page
Parameter / Method Unit Value Comment
c′µ 1 0.072 constant for calculating turbulent
diffusivity for temperature
σk 1 1 constant for calculating turbulent
diffusivity for k
σε 1 1.3 constant for calculating turbulent
diffusivity for ε
Appendix C
Chemsitry Model
C.1 Model Description
The chemistry model is developed by Sergei Katsev and Itay Halevy (PNAS in prep). The
following paragraphs explained the mechanism of this model.
The basis of this chemistry model is constructed by a system of partial differential equations
represented in Eq. C.1. The solution is obtained with MATLAB pdepe solver (reference to
Skeel and Berzins (1990)).
∂C j
∂ t
=
∂
∂ z
(
Kz(z)
∂C j
∂ z
)
− v j ∂C j∂ z +∑i
νiRi j (C.1)
Here, Kz(z) is the depth-dependent turbulent diffusivity, v j is advection velocity, but in this
case, only refers to the settling velocity. Ri j is the ith reaction rates for chemical species j. Its
stoichiometric coefficient for chemical j is νi. Species considered for constructing the systems
of equations includes POM, Ac, O2, SO
2–
4 , TS H2S + HS
–, FeOOH, Fe2+, S0, FeS, and the
biomass of the microbial populations that catalyze the corresponding reactions: XIR, XSR, XCH4,
89
90
XS–AOM, XFe–AOM. The reactions involved are listed in table C.2. DIC, ammonium, and pH are
assumed to be constant. Microbial reactions are described in the following bioenergetic model.
The bioenergetic model describes the energy of microbial metabolism as:
∆G = λcat∆Gcat +∆Gana (C.2)
where ∆Gcat is the Gibbs free energy of catabolism process, which generates energy for
microbial consumption, and ∆Gana describes the anabolism process that microbes consume
energy for building up biomass. ∆G is thus the energy being dissipated into the environment.
λcat is the number of times a catabolism process is required to provide the energy needed for
forming one C-mol of biomass with anabolism. Thus, the maximum growth yield by per mol
of reaction is:
Ymax =
1
λcat
=
∆Gcat
∆G−∆Gana (C.3)
Cell-specific cataboblic reaction rates (r) is a function of substrate (S) concentrations:
r =
1
X
dS
dt
= FKFT (C.4)
Where X is the biomass in C-mols. FK is the kinetic factor, and FT is the thermodynamic factor.
FK depends on maximum cell-specific rate Vmax and substrate concentrations:
FK =Vmax∏ SiKim+Si (C.5)
FT is formed based on Gibbs free energy of catabolic reaction and ATP energy release:
FT = 1− exp
(
∆Gcat +m∆GAT P
χRT
)
(C.6)
91
Table C.1: Chemical model boundary conditions
Parameter Value Unit Comment
top boundary
[O2] 200 µmol l−1 fixed concentration of oxygen
[SO4] 40 µmol l−1 fixed concentration of sulfate
F0FeIII 950 mmol/m
2/yr surface flux of FeIII
F0G 0.1 mmol/m
2/yr surface flux of OM
bottom boundary
FbottomFeII frec*F
bottom
FeIII mmol/m
2/yr bottom flux of FeII
*frec: fraction of FeIII being recycled, it is related to organic matter flux, Fe(III)
flux, oxygen concentration, and H2S concentration
Thus, variations in chemical species (S) and biomass (X) can be expressed as :
dSs
dt
= AsmrmbXb (C.7)
dXb
dt
=∑
n
YnmrmbXb−λ (mG)Xb (C.8)
Asm is the stoichiometric coefficient for substrate s in a catabolic reaction m. rmb is the rate of
this reaction. Xm is the active biomass involved in stimulating this reaction. Term ∑nYnmrmbXb
describes the yield of biomass Xb associated with anabolic reaction n and catabolic reaction m.
Term −λ (mG)Xb on the other hand, describes the death of microbes related to the maintenance
energy mG.
The important boundary conditions of the model are illustrated in Fig. C.1 and are defined
in Table C.1. Chemicals not mentioned here have surface boundary conditions of zero
concentration for dissolved species and zero flux for particulate species. The bottom boundary
conditions are zero flux for dissolved species. Solid species sink out when they reach the
sediment-water interface at their settling rate defined by vsettle in Table C.4. A fraction of
Fe(III) precipitated into the sediment is reduced and recycled as Fe2+.
92
Figure C.1: Illustration of the important boundary conditions and reactions of the
biogeochemistry model.
The primary production is assumed to be distributed in the upper water column according to an
analytical exponential function:
PP(z) = PPtotalke−kz (C.9)
where k = 0.0655, assuming the abundance of primary production decreases at the same rate
as light attenuation. PP0 is 8300 µmol/L/yr, the total primary production of the lake integrated
with depth. This is an assumption based on Lake Towuti being an oligotrophic lake. This
process produces POM, which remineralizes following a power law of organic age and available
93
oxygen level (Katsev and Crowe, 2015):
k(z, [O2]) = f kox+(1− f )kanox (C.10)
where the partitioning coefficient is f = [O2]/(KO2 + [O2]). The reactivities for the oxic and
anoxic mineralization are (Katsev and Crowe, 2015):
kox = 10−0.312τ−0.977 (C.11)
kanox = 10−1.1τ−0.857 (C.12)
The age τ of organic material (in years) was calculated at each depth as:
τ = τinit +
∫ z
0
dz
vdi f f + vsettleG
(C.13)
Where vdi f f is the effective downward velocity of organic matter, due to turbulent diffusion,
which can be approximated as:
vdi f f =
1
z
Kz
− ( z2Kz )2
∂Kz
∂ z
(C.14)
vsettleG is the linear settling velocity of organic matter, defined in Table C.4.
94
C.2 Biogeochemical Model Reactions
Table C.2: Chemical model reactions
Reaction ∆G0 Comment
POM−−→ CH3COO– hydrolysis
CH3COO– + 2O2 −−→ 2HCO3 – + H+ remineralisation
SO42 –
SulfateReducer−−−−−−−−→ HS– sulfate reduction
FeOOH IronReducer−−−−−−→ Fe2+ iron reduction
SO 2–4 + CH4 HS
– + HCO –3 + H2O -33 anaerobic oxidation of
methane (via sulfate)
8FeOOH + CH4 + 15H+ −−→ 8Fe2+ + HCO3 – + 13H2O -117 anaerobic oxidation
of methane (via iron
oxides)
H2S + 2O2 −−→ SO42 – + 2H+ -830 sulfide oxidation
4Fe2+ + O2 + 6H2O−−→ 4FeOOH + 8H+ -371 iron oxidation
2FeOOH + H2S + 4H+ −−→ S0 + 2Fe2+ + 4H2O -32 sulfide oxidation via
iron oxides
Fe2+ + HS– −−→ FeS + H+ precipitation
4S0 + 4H2O−−→ SO42 – + 3HS– + 5H+ +55
95
C.3 Biogeochemical Model Reaction Rates
Table C.3: Chemical reaction rates
Rate
RPPO2 = RPP
RG = k(z, [O2])[POM]
RO2 = RG
[O2]
[O2]+KO2
RHSox = kHSox[H2S][O2]
RFeox = kFeox[Fe2+][O2]
RHSFeIII = kHSFeIII[FeOOH][H2S]
RFeS = 5×102 ∗ (ΩFeS−1)
where ΩFeS =
[Fe2+][HS]
[H+]KFeS
HS =
[H2S]
1+[H+]/KHS
RCH4ox = kCH4ox[CH4][O2]
Rdisp = kdisp[S0]
96
C.4 Biogeochemical Model Parameters
Table C.4: Chemical model parameters
Parameter Value Unit Comment
τinit 0.5 y initial age of organic matter
PP 22 mmol/m2/d gross primary production
KAc 5 uM affinity for acetate for SR
15.5 uM affinity for acetate for IR
12 uM affinity for acetate for methanogens
KO2 2 uM affinity for oxygen
KO2inh 1 uM inhibition of anaerobes by O2
KSO4 5 uM affinity for sulfate
KCH4 30 uM affinity for methane (for AOM)
KFeIII 1000 uM affinity for goethite
KFeR 20 uM affinity for nanophase precipitates
KHS 9.3E-8 uM acid dissociation constant for H2S
KFeS 2.5E-3 uM affinity for iron sulfide
vsettle 1800 m/y settling velocity
vsettleG 600 m/y settling velocity for organic matter
frec 0.1 constant related to fraction of FeIII being
recycled from the sediment
kHSox 1 1/uM/y rate constant
kFeox 1 1/uM/y
kHSFeIII 2×10−1 1/uM/y
kFeS 5×10−1 1/uM/y
kCH4ox 1 1/uM/y
kdisp 2×10−2 1/y
T 29 oC Temperature
Continued on next page
97
Table C.4 – continued from previous page
Parameter Value Unit Comment
pH 7 pH
∆GAT P 30 kJ/mol cost of producing ATP (for FT )
χ 8,6,2 - stoich. number of Jin&Bethke
m 1,5/4,1/4 - stoich. number of Jin&Bethke
Vmax 1.0 molrx/molx/h microbial kinetics for SR
× 1/10 microbial kinetics for IR
× 1.0 microbial kinetics for meth.
× 1/25 microbial kinetics for S-AOM
× 3/10 microbial kinetics for Fe-AOM
× 1/50 microbial kinetics for FeR-AOM
λdeath 2E-4 1/h death constant
kdeath 2E-3 1/h death rate from energy starvation
Pmaint 4.56 kJ/molx/h min. maintenance power
Kanalim 1E-2 uM limiting conc. for anabolism
HCO3 – 3E-3 M
NH4+ 1E-5 M nutrient for anabolism
Appendix D
Extra Figures
Figure D.1: Root mean squared error (RMSE) between model and mooring (2012 - 2015)
temperature data as a function of time. The dashed line indicates an average RMSE of the
dataset. The solid vertical line represents June 2014. After this time, the confidence in mooring
data decreased. This figure illustrates that, RMSE significantly increased during 2012-2013 wet
season. However, this increasing of error did not accumulate in the proceeding simulation.
98
99
Figure D.2: RMSE between model and mooring (2012 - 2015) temperature data as a function
of depth. Dashed line indicates the overall RMSE of the dataset. RMSE is larger than average
at the surface, and is minimal at around 100 m. This indicates that Simstrat simulation errors
are mostly observed at the surface rather than deeper part of the lake and that the modeling of
surface heat flux over Lake Towuti results in discrepancy between simulation and observation.
Bottom temperature simulated by simstrat is more satisfactory.
100
Figure D.3: RMSE sensitivity with meteorological variables. This figure is generated with
the same data in 2.13, except the y-axis is now plotting the RMSE between the model and
mooring data. This figure illustrates that the current configuration of meteorological variables
is a satisfactory representation that allows simulation agreement with observation. No better
combinations of meteorology can significantly improve the performance of the model.
101
Figure D.4: Cross correlation between meteorological variables and MEI.v2 index. negative
lag indicates MEI.v2 is prior to the local meteorological variable. The correlations are mostly
weak (with Pearson’s r smaller than 0.5 for all the climate variables). MEI.v2 index represents
the air temperature over Lake Towuti several months later, but has delays in representing the
wind speed, solar radiation, vapor pressure, and cloud cover over Lake Towuti.