Data Assimilation for Marine Ecosystem Models

Data Assimilation for

Marine Ecosystem Models

Data Assimilation for

Marine Ecosystem Models

Proefschrift

ter verkrijging van de graad van doctor

aan de Technische Universiteit Delft,

op gezag van de Rector Magnificus prof. ir. K.C.A.M. Luyben;

voorzitter van het College voor Promoties

in het openbaar te verdedigen op dinsdag 16 juli 2013 om 10.00 uur

door

Joanna Sylwia PELC

Master of Science in Applied Mathematics

geboren te Nowa S´ol, Polen.

Samenstelling promotiecommissie:

Rector Magnificus voorzitter

Prof. dr. ir. A. W. Heemink Technische Universiteit Delft, promotor Prof. dr. A. M. Moore University of California, Santa Cruz Prof. dr. ir. A. E. Mynett UNESCO-IHE

Prof. dr. R. W. P. M. Laane University of Amsterdam Prof. dr. ir. H. X. Lin Technische Universiteit Delft Prof. dr. ir. H. H. G. Savenije Technische Universiteit Delft Dr. ir. G. Y. El Serafy Deltares

Dit onderzoek kwam tot stand met steun van Deltares

Dr. ir. Ghada Y. El Serafy heeft als begeleider in belangrijke mate aan de totstandkoming van het proefschrift bijgedragen

isbn 978-94-6203-388-7

Copyright c° 2013 by Joanna S. Pelc

All rights reserved. No part of the material protected by this copyright notice may be reproduced or utilized in any form or by any means, electronic or mechanical, including photocopying, recording or by any information storage and retrieval system, without the prior permission of the author.

Typeset by the author with the LA_{TEX Documentation System.}

Acknowledgments

I owe my deepest gratitude to all the people I worked with during my PhD re-search. This PhD thesis is a result of support and guidance from many people.

Firstly, I would like to thank my supervisors Arnold Heemink and Ghada El Serafy for giving me this opportunity, believing in me, and having so much patience for my stubborn perfectionism. Ghada, thank you for your warm sup-port and freedom you gave me during my research. Arnold, thank you for your heartwarming guidance, and for your positive attitude in every situation.

I want to thank my work colleagues in Deltares and Delft University of Tech-nology for a friendly and collaborative atmosphere. A special thanks to Hans Los, thank you for your enthusiastic lessons and discussions about ecosystem models, and for your support and warm encouragements throughout my research. Many thanks to Jan van Beek, thank you for your assistance with the model, and solving many of my Delft3D emergencies. Also, I would like to thank Sandra Gaytan, thank you for your help and assistance with the satellite data. I would also like to thank Martin Verlaan, thank you for translating the abstract of this thesis into Dutch.

During my PhD work I had an opportunity to collaborate with great people at Nansen Environmental and Remote Sensing Center in Bergen and at University of California in Santa Cruz. I am very grateful to Ehouarn Simon, Laurent Bertino, Jerome Fiechter and Andrew Moore, thank you for having me there and for the very fruitful collaborations, I have learned a lot while working with you all.

I would like to thank the University of Zielona G´ora for giving me the solid background in mathematics. My sincere gratitude to Jolanta Misiewicz and Roger Cook for creating a bridge between the University of Zielona G´ora and Delft University of Technology, which gave me the opportunity to complete my master of science in applied mathematics at Delft University of Technology. Without them I would not be where I am today.

Important thanks goes to two special friends – Ivan Garcia and Remus Hanea, without whom my PhD times would not have been the same. Guys, thank you for all the good ’musketeer’ times we shared. Ivan, we have started the PhD research

and especially to have someone to discuss with, and sometimes also to disagree with. I particularly enjoyed our engineer vs mathematician discussions, which often showed me an ”orthogonal” perspective. Remus, there is no words I could use to express my gratitude to you. Thank you for your vivid and motivating guidance throughout all my time in Delft, as my master supervisor, and as my friend. My scientific path has been in a great deal ”data-assimilated” by you. Without you, I would not be walking this path. Thank you!

During all this time, I was pursuing my passion for dancing, by training and performing with salsa show teams. I would like to thank to both of my teams, SoSalsa! and Drankenwoede, for all the great times we shared, and for being like a family to me. I would also like to thank all the party animals and salsa addicts for all dances, parties and craziness, thank you all for keeping me sane through all this time.

A special mention goes to all the little fairies, who were cheering me up throu-ghout my PhD times by surprising me with chocolates, cakes, roasted almonds, and happy notes – thank you Dana Stuparu, Jos´e Reinders, Dick Vethaak, and Aditya Kamalapurkar. Sometimes the smallest things bring the biggest smile.

I would like to thank all my friends, who have supported me and made me feel home during all these years in Holland. Especially, I want to thank Kasia KrugÃla and my ’Love Board’: Aurinke, Hendrik, Hiske, Jos´e, and Roosmarijn, guys thank you for always being there for me.

Finally, I want to thank my family, for all their support, encouragements and good vibes which always keep me going. A special thank to my grandma Alfreda KÃl¸ebek for ingraining in me the endless craving to learn, and my mom Krystyna Czubik for inspiring in me ambition, strength and faith, that I can always achieve anything I want.

Joanna Pelc Delft, June 2013

Summary

In the presence of environmental issues caused by climate change and eutrophi-cation, ecosystem models capable of providing accurate predictions are of consi-derable interest. Even the most advanced ecological models are not capable of reproducing reality perfectly, therefore it is essential to integrate them with the available observations. Data assimilation techniques serve as tools to calibrate and improve model accuracy by combining it with the data from measurements. Although ecosystem models have evolved to be considered as relatively pro-ficient predictive tools, they are still missing a unified set of rules to govern the whole system. This is compensated by adding various biogeochemical processes associated with a substantial number of uncertain parameters, which are an im-portant source of uncertainty in ecosystem models. To date, these have been successfully established through manual calibrations, with help of available mea-surements. However, the manual process is often time consuming, and becomes more complicated as the number of parameters to be estimated increases. Moreo-ver, the parameter values are usually dependent on the region for which they are calibrated, therefore the process often needs to be repeated when a new ecosystem region is considered. Data assimilation methods provide a semi-automatic tool for estimating many parameters simultaneously, which when combined with expert knowledge can be a significant improvement over manual calibration.

Since the parameters are a primary source of uncertainty, their estimation is the main focus of this work. Additionally, as the initial conditions are not always well known, these are also included and updated together with the parameters. Based on available data assimilation methods, variational techniques were selected as the most suitable. These were first introduced in meteorology for initial condi-tion estimacondi-tion, and further it also proved to be a proficient tool for parameter estimation. The same method has already been used successfully for ecosystem calibration. However variational techniques require the implementation of the model derivatives, which is a challenge for the ecological models, since they are renowned for their high degree of nonlinearity and issues of non-differentiability. Moreover, as ecosystem models become more advanced and sophisticated, it is

important to search for alternative solutions for obtaining their adjoint models. To address this issue, a model reduced four dimensional variational data as-similation (4D-Var) is used. The method was designed in such a way that the implementation of the adjoint of the tangent linear approximation of the original model is not required. Using an ensemble of simulations from the original model, a proper orthogonal decomposition is used to build a reduced model. Then the adjoint of the original model is approximated by the adjoint of the tangent linear approximation of the reduced model. The method is easily adapted to tackle initial condition estimation, which is proficiently estimated in the reduced space. To gain a better understanding of the model reduced 4D-Var method and its behavior in a marine ecosystem application, the technique is first applied to a simple 1D ecological model. The method is demonstrated in a twin experiment framework, where synthetically generated surface phytoplankton data are used. Three parameters are calibrated in a combination with an estimation of the initial conditions. Different control strategies are explored, showing that updating the initial condition is essential to obtain accurate parameter calibrations, which moti-vates their combined estimation in the further applications. Although the system used to investigate the method is small, it allows us to address and identify impor-tant issues when applying data assimilation methods in ecosystem models before integration in real complex configurations. Based on the relatively good results obtained in the 1D ecological model, the next step is to apply the model reduced 4D-Var method and assess its performances in a realistic large scale ecosystem.

Next, the technique is implemented in the 2D North Sea coupled physical-ecological model (BLOOM/GEM). The model consists of detailed hydrodynamics, suspended sediment and river loads. Among many other processes, it takes into account primary production, nutrient cycle, phytoplankton species competition and their adaptation to different limiting conditions. Justified by generally well mixed conditions in the North Sea a vertically mixed 2D version of the model is used. It is a result of many years of experience, and it has been repeatedly calibrated and validated.

The main goal is to improve the model capability to predict the chlorophyll-a concentrchlorophyll-ation, which is the mchlorophyll-ain indicchlorophyll-ator of the wchlorophyll-ater quchlorophyll-ality. Followed by sensitivity analysis, supported by expert opinion, a number of parameters were selected as the most significant in the BLOOM/GEM model. A number of twin experiments are performed aimed at the calibration of two leading parameters from the sensitivity list. Similar control strategies as in the simple ecosystem case were applied, confirming that estimation of the initial condition enhances the parameter calibration in terms of accuracy and computational efficiency. Pa-rameter estimation resulted in 80%-90% of improvement with respect to the prior error (100% corresponds to the error of the prior parameters), whereas the initial condition updated in the same experiments obtained a 10%-25% of improvement with respect to the prior error. Moreover, the results show that the method is capable of significantly improving the prediction of chlorophyll-a, reaching up to 35% of total improvement over two years. The relatively good performance of the

model reduced 4D-Var method in the 2D North Sea BLOOM/GEM model de-monstrates its potential as a calibration tool for enhancing the model predictions using real chlorophyll-a measurements derived from the remotely sensed MERIS data. Therefore, further experiments are performed, this time with the real data, following the same control vector strategies as for the twin experiments, where the parameters and initial conditions are estimated. The predictions of chlorophyll-a concentration resulting from data assimilation are validated against remotely sen-sed MERIS measurements during a nearly two year period after the assimilation. The performance of the assimilated model is enhanced with respect to the original model, showing that for some control strategies up to 10% improvement of the model occurs, which verifies the model reduced approach to be a useful tool for improving model predictions of chlorophyll-a.

Samenvatting

Data Assimilatie voor Mariene Ecosysteem Modellen

Ecosysteemmodellen, die in staat zijn om nauwkeurige voorspellingen te leveren, zijn van groot belang voor milieuproblemen veroorzaakt door klimaatverande-ring of eutrofi¨eklimaatverande-ring. Zelfs de meest geavanceerde ecologische modellen zijn niet in staat om de werkelijkheid perfect na te bootsen, daarom is het gebruik van metingen essentieel. Met behulp van data-assimilatie methoden kunnen modellen worden gekalibreerd en kan de nauwkeurigheid worden verbeterd door metingen in het model te integreren. Hoewel ecosysteemmodellen zich ontwikkeld hebben tot een relatief betrouwbaar voorspellend gereedschap, is er nog geen overkoepe-lende theorie die het hele systeem omvat. Hiervoor wordt gecompenseerd door een groot aantal onzekere parameters te introduceren voor de verschillende bioche-mische processen, wat een aanzienlijke onzekerheid in deze modellen veroorzaakt. Tot op heden heeft men deze parameters met succes gekalibreerd met behulp van de beschikbare metingen. Dit handmatige werkproces is vaak tijdrovend en wordt gecompliceerder naarmate het aantal parameters toeneemt. Bovendien zijn de waarden van de parameters afhankelijk van het gebied waarvoor zij zijn geka-libreerd, waardoor de kalibratie moet worden herhaald als men een andere regio wil beschouwen. Data-assimilatie methoden bieden de mogelijkheid om op een semiautomatische manier een groot aantal parameters tegelijkertijd te schatten. Gecombineerd met kennis van experts kan dit een belangrijke verbetering opleve-ren in vergelijking met handmatige kalibratie.

Omdat de parameters de grootste bron van onzekerheid vormen, is het schat-ten ervan het hoofddoel van dit proefschrift. Omdat ook de beginvoorwaarden niet altijd bekend zijn worden deze toegevoegd en samen met de parameters ges-chat. Van de beschikbare data-assimilatiemethoden zijn de variationele metho-den geselecteerd als het meest geschikt. Variationele methometho-den wermetho-den het eerst ge¨ıntroduceerd in de meteorologie voor het schatten van beginvoorwaarden, maar ze bleken daarnaast ook heel geschikt voor het schatten van parameters. Hoewel deze methoden zijn ook met succes toegepast voor de kalibratie van

modellen, blijft het bepalen van de benodigde modelafgeleiden een uidaging omdat ecosysteemmodellen sterk niet lineair zijn en er bovendien vaak problemen zijn met de differentieerbaarheid. Tenslotte worden ecosysteemmodellen steeds gea-vanceerder en complexer, waardoor het belangrijk wordt om te zoeken naar een alternatieve manier om de geadjungeerde van het model te bepalen.

Als alternatief wordt hier een modelreductie methode, model reduced four dimensional variational data-assimilation (4D-Var), gebruikt. De methode is zo ontworpen dat het implementeren van de geadjungeerde van het gelineariseerde model niet nodig is. Met een ensemble van simulaties wordt eerst de dominante deelruimte bepaald door berekening van de ”proper orthogonal decomposition”, waarmee vervolgens het model wordt gereduceerd tot een veel kleinere omvang. De geadjungeerde van het oorspronkelijke model wordt benaderd met de geadjun-geerde van het gereduceerde model, die eenvoudig is te bepalen. Om een beter inzicht te krijgen in het gedrag van model reduced 4D-Var voor mariene ecosys-teemtoepassingen, wordt de techniek eerst toegepast op een simpel 1D ecologisch model. De methode wordt gedemonstreerd aan de hand van een tweeling ex-periment, waarbij kunstmatig gegenereerde metingen van fytoplankton worden geassimileerd. Er worden tegelijkertijd zowel drie parameters als de toestand van het model geschat. Er worden diverse varianten verkend, die laten zien dat het schatten van de beginconditie essentieel is voor het nauwkeurig schatten van de parameters, wat een goede reden is om beide tegelijk te schatten in toekomstige toepassingen. Hoewel dit model vrij klein is stelt het ons toch in staat om de belangrijkste factoren te identificeren voor de overstap naar een meer complex model. Uitgaande van de goede resultaten voor dit kleinere 1D ecologische mo-del is toepassing van momo-del reduced 4D-Var voor een toepassing van realistische omvang een logische volgende stap.

Vervolgens is de methode ge¨ımplementeerd in het gekoppelde fysisch-ecologisch model (BLOOM/GEM) voor de Noordzee. Het model bevat een gedetailleerde hydrodynamica, zwevend sediment en riviertoevoer. Naast vele andere processen, word rekening gehouden met primaire productie, de nutri¨enten kringloop, compe-titie tussen soorten fytoplankton en de aanpassing aan verschillende limiterende omstandigheden. Door de goed vertikaal gemengde omstandigheden in de Noord-zee is het gebruik van een 2D model goed te rechtvaardigen. Dit model is het resultaat van vele jaren ervaring en herhaaldelijk kalibreren en valideren.

Het belangrijkste doel van deze studie is het verbeteren van de mogelijkheden van het model om de chlorofyl-a concentratie te voorspellen, wat de belangrijkste indicator voor waterkwaliteit is. Na een gevoeligheidsstudie, is met expertkennis een selectie gemaakt met de meest invloedrijke parameters in het BLOOM/GEM model. Eerst zijn er tweelingexperimenten met de belangrijkste 2 parameters uit-gevoerd, die met een vergelijkbare aanpak als voor het eenvoudiger 1D model, laten zien dat het gecombineerd schatten van parameters en beginvoorwaarden zowel de nauwkeurigheid als rekentijd verbetert. Het gecombineerd schatten van parameters en de begintoestand resulteerde in een relatieve verbetering van 80%-90% ten opzichte van de nauwkeurigheid (100%) zonder data-assimilatie, daarbij

xi werd de nauwkeurigheid van de begintoestand 10%-25% hoger. Bovendien wordt ook de nauwkeurigheid van de voorspelling van chlorofyl-a beduidend beter, tot een verbetering van gemiddeld 35% over een periode van twee jaar. De goede pres-taties van model reduced 4D-Var voor het BLOOM/GEM model voor de Noord-zee toont de kracht van de methode om modelvoorspellingen van chlorofyl-a te verbeteren met metingen van de MERIS satelliet. Daarom zijn er tenslotte expe-rimenten uitgevoerd met echte metingen, waarbij dezelfde strategie werd gevolgd als voor de eerdere experimenten, waarbij parameters en en beginvoorwaarden tegelijk werden geschat. De modelvoorspellingen van de chlorofyl-a zijn gevali-deerd met MERIS gegevens gedurende een periode van bijna 2 jaar volgend op de data-assimilatie. De prestaties van het model met data-assimilatie zijn duidelijk beter dan het originele model, waarbij in een aantal gevallen de verbetering tot 10% oploopt. Dit laat zien dat model reduced 4D-Var een nuttige aanpak is voor het voorspellen van chlorofyl-a.

Contents

Chapter 1

Introduction

1.1 Marine ecosystem background . . . 1

1.2 Motivation . . . 3

1.3 Main objective . . . 5

1.4 Overview of the thesis . . . 6

Chapter 2

Primary production model

/

2.1.1 Competition and adaptation . . . 12

2.1.2 Phytoplankton processes . . . 14

2.1.3 Nutrient cycling processes . . . 14

2.1.4 Oxygen processes . . . 16

2.1.5 Light extinction . . . 16

2.1.6 Mineralization . . . 17

2.1.7 Sediment processes . . . 17

2.1.8 Grazing . . . 17

2.2 North Sea application . . . 18

2.2.1 Schematization . . . 18

2.2.2 Hydrodynamics . . . 19

2.2.3 North Sea specific processes . . . 20

Chapter 3

Application of model reduced 4D-Var to a 1D ecosystem

3.2 1D Ecological Model . . . 24

3.3 Variational data assimilation . . . 25

3.3.1 Incremental 4D-Var . . . 27

3.3.2 Model Reduced 4D-Var . . . 30

3.4 Framework of the experiments . . . 31 xiii

3.4.1 Twin experiment setup . . . 32

3.4.2 Observations . . . 33

3.4.3 Background error covariance matrices . . . 33

3.4.4 Minimization . . . 34

3.4.5 Reduced model setup . . . 34

3.5 Results and Discussion . . . 36

3.6 Conclusions . . . 46

Chapter 4

Model reduced 4D-Var for

/

: twin experiments

4.2 Primary production model BLOOM/GEM . . . 53

4.2.1 Water quality and ecological processes . . . 53

4.2.2 North Sea application . . . 55

4.3 Model reduced variational data assimilation . . . 57

4.3.1 Model Reduced 4D-Var . . . 58

4.4 Framework of the experiments . . . 59

4.4.1 Flow chart of the algorithm . . . 60

4.4.2 Minimization constraints . . . 61

4.4.3 Twin experiment setup . . . 61

4.4.4 Observations . . . 62

4.4.5 Background error covariance matrices . . . 62

4.5 Results and Discussion . . . 63

4.5.1 Energy scenarios . . . 63

4.5.2 Short assimilation window . . . 67

4.5.3 Long assimilation window . . . 76

4.5.4 Comparison . . . 80

4.5.5 Computational cost . . . 81

4.6 Conclusions . . . 83

Chapter 5

Model reduced 4D-Var for

/

:

5.2 MERIS derived chlorophyll-a measurements . . . 87

5.3 Primary production model BLOOM/GEM . . . 91

5.3.1 Ecosystem components . . . 92

5.3.2 Model configuration . . . 93

5.3.3 Boundaries and river discharges . . . 93

5.4 Model reduced variational data assimilation . . . 93

5.5 Framework of the experiments . . . 95

5.5.1 Control vector strategies . . . 96

5.5.2 Assimilation window . . . 96

5.5.3 Constraints . . . 96

Contents xv

5.6 Results and discussion . . . 99

5.6.1 Two biological parameters . . . 99

5.6.2 Five biological parameters . . . 106

5.6.3 Comparison . . . 109

5.6.4 In situ comparison . . . 111

5.7 Conclusions . . . 112

Chapter 6

Conclusions and recommendations

Chapter A

Projection matrix P

Chapter B

Parameter transformation

Chapter 1

Introduction

”How inappropriate to call this planet Earth when it is quite clearly Ocean.

”

— Arthur C. Clarke

1.1 Marine ecosystem background

Environmental decay was triggered by human activity already many years ago. With its dramatic impact on the ocean, it has altered the biological content of the waters, which have become more amenable to harmful algal blooms. It threatens the natural wildlife and water quality. Moreover, the whole process is enhanced by climate change and the constant growth of the human population.

Followed by the rapid growth of the population, an increased demand for food has resulted in global over fishing of the marine waters (Watson and Pauly, 2001b; Pauly et al., 2002). While it has been taking place already for many cen-turies ago, it has accelerated dramatically since 1950s due to the development of new technologies (Watson and Pauly, 2001a), resulting in over fishing and major ecological extinctions (Jackson et al., 2001). Humans have become the top preda-tor in coastal food webs, which has changed their structure dramatically, causing a decline in many species which used to be responsible for the natural filtering of the water. The plankton grazers were so efficient that the symptoms of eutro-phication were limited. However, over fishing has caused nearly their extinction, which resulted in an increase of phytoplankton, jellyfish and enhanced microbial population (Jackson et al., 2001).

The forests, especially those on the coast, are very important as they are natural filters for the nutrient run off from lands. Due to the increased demand for food and housing, forests are being replaced by agriculture and buildings, resulting in increased use of fertilizers, deforestation and sewage pollution. These factors have led to an excess of nutrients in estuaries, known as eutrophication, leading to enormous algal blooms, which often result in oxygen depletion during

decay. This creates areas low in dissolved oxygen, sometimes even anoxic (dead zones), and when the concentration falls below 30%, the water conditions become lethal for most of the sea life. Such environment is suitable for bacteria and jellyfish, which have a negative influence on other species.

In order to mitigate the harmful effects of eutrophication, many countries have initiated programs to remove or reduce nutrient accumulation in the waters. Agricultural runoff enriches the water mainly with nitrates and phosphates, ho-wever the debate over which nutrients should be controlled has resulted in much controversy (Howarth and Marino, 2006). Many scientists argue that the key nutrient responsible for eutrophication is nitrogen (Brush, 2008). On the other hand, Schindler et al. (2008) carried out an experiment for 37 years, to validate the theory that controlling nitrogen inputs can control eutrophication, and found it was phosphorus not nitrogen that is the critical nutrient that should be control-led. Their study showed, that in some cases removal of nitrogen even enhanced primary production. Even though opinions were divided, the control of phosphate was initiated in freshwaters and estuaries in 1960s, whereas the control of nitrates did not begin until 2001 (Howarth and Marino, 2006).

Biological modeling received even more attention when the algal bloom fre-quency, intensity and range expansions began to change enhanced by the climate change (Hallegraeff, 1993, 2010). According to Paerl and Huisman (2008) global warming amplifies eutrophication, and moreover, warmer waters are more ame-nable to blooms. A particular threat is due to the harmful algal blooms. The dominant toxic species, cyanobacteria (blue-green algae), known as the first plant on Earth to produce oxygen, now threatens human health and ecosystem life (Paerl and Huisman, 2008). Hallegraeff (2009) also draws attention to an increa-sing harmful algal blooms due to the climate change, and highlights the fact that different species of algae react quite differently to increases in ocean temperature. As a result we should expect range expansions of warm-water species, that may out compete cold-water species. This may result in migrations and significant expan-sions in the range of some species, followed by unexpected harmful algal blooms in the areas which are currently poorly monitored (Hallegraeff, 2010). Particu-larly vulnerable are the areas of agricultural nutrient run off, increasing water temperatures, hurricanes and many other factors, where current ocean conditions tend to favor toxin producing species (Hallegraeff, 2010; Peperzak, 2003).

Marine ecologists, Jeremy Jackson warns that we are experiencing environ-mental decay, which drastically alters the ocean. This process started many years ago, mainly due to human activity, but has gone unnoticed. As a consequence, ocean conditions are becoming more favorable for bacteria and jellyfish, just as it was a half-billion years ago (Kenneth R. Weiss ”A Primeval Tide of Toxins”, July 30, 2006, Los Angeles Times).

Castle and Rodgers Jr. (2009) carried out research based on the ancient algal observations from the fossil records, and found that extremely abundant cyano-bacteria blooms coincided with the major mass extinction. They concluded that cyanobacteria was an important factor during the major past mass extinctions,

Introduction 3 and also stated that the climate conditions during these events were similar to current climate conditions.

On the other hand, Boyce et al. (2010) studied archives of data available from 1899 to present, and found that ocean warming is correlated with a long-term de-clining trend in phytoplankton biomass. This is yet another debatable statement, as opposed to the scientists warnings about the rapid incline of the algal blooms. Despite which of these opinions is true, they demonstrate that marine ecosystems are complex, and we have a lot to learn.

Although there is still a lot to understand about the marine ecosystem dyna-mics, it is commonly acknowledged that environmental conditions are decaying, and that enhanced monitoring and management strategies are required to improve the description of the biological content of the waters given by ecological models.

1.2 Motivation

In the presence of environmental issues caused by climate change and eutrophi-cation, ecosystem models capable of providing accurate predictions are of consi-derate interest. Despite of many modeling challenges (Doney, 1999; Jørgensen, 2008), ecological models have evolved to be considered as relevant predictive tools. However, most sophisticated systems are not able to reproduce the reality com-pletely, therefore it is essential to integrate them with the available observations. Data assimilation techniques serve as a tool to calibrate and improve model ac-curacy by combining them with the data given by the measurements.

Ecological models, however, while often very advanced and complicated are still missing a unified set of rules, that govern the whole system (Jørgensen, 2008; Los, 2009; Doron et al., 2011). This is often compensated for by adding extra processes, which come with a large number of extra parameters, and these of-ten introduce significant uncertainty (Doron et al., 2011). That is the reason why an important source of uncertainties in ecosystem models is assigned to the poorly known parameters. With a given observation set, these could be estimated through laboratory experiments. However, when a large number of parameters is considered, it is much harder to manage the calibration problem manually. Mo-reover, the calibration needs to be repeated every time a new ecosystem region is considered. Data assimilation methods provide tools to estimate many parameters simultaneously, with the possibility of accounting for their correlations.

There are two distinct classes of data assimilation methods, both capable of accounting for imperfect parameters. One is the class of variational (inverse) techniques, which search for an optimal set of control variables, such that a cost function which measures the distance between the model and observations and the model and a prior is minimized. The second class is represented by forward methods, which assimilate data sequentially in time. Hence, they are often refer-red to as sequential methods. Their main objective is to correct the estimated variables at every time an observation becomes available. Despite the differences

between the two classes of methods, both methodologies are equivalent for linear systems (Lorenc, 1986).

Variational techniques have been successfully applied to improve ecosystem predictions. Several approaches have been used, such as ecosystem state estima-tion (Natvik et al., 2001), updating the input of a biological model by improving its coupled hydrodynamical model (Fiechter et al., 2011), and the calibration of the ecological parameters. The latter can be classified in terms of the optimization technique used, among which are simulated annealing (Matear, 1995), gradient steepest descend (Fennel et al., 2001), Green’s function method (Miller et al., 2000), and the most popular the adjoint method (Friedrichs, 2002; Zhao et al., 2005). Sequential methods were also applied to perform data assimilation for ecosystems, mainly for the state estimation. One of the first sequential methods used in ecosystem applications was the Extended Kalman Filter (Carmillet et al., 2001), soon after very common became the Ensemble Kalman Filter (Allen et al., 2002; Eknes and Evensen, 2002; Natvik and Evensen, 2003; Simon and Bertino, 2009). For an overview of sequential methods used in oceanography and ecology see Bertino et al. (2003). Gaussian anamorphosis extensions of ensemble-based Kalman filters have been suggested by Bertino et al. (2003) to tackle the pro-blems of sub-optimality of the filter arising from the non-Gaussian distributions of most of the biological variables and parameters. These approaches can be ea-sily applied in realistic configurations (Simon and Bertino, 2009) and have been proved to be efficient tools for calibration poorly known parameters (Doron et al., 2011; Simon and Bertino, 2012) For a detailed overview of data assimilation me-thods in application to biological models see Gregg (2008).

In this work the aim is the simultaneous estimation of constant parameters and the initial conditions, therefore a variational technique is a suitable choice. The method used in this work is four dimensional variational data assimilation (4D-Var), an adjoint method, first introduced in meteorology for the initial condition estimation problem (Le Dimet and Talagrand, 1986; Talagrand and Courtier, 1987). However, it has proven to be a powerful tool for ecosystem calibration, especially as these models become more and more sophisticated. The number of their bio-logical components, as well as their parameters is increasing, and even simple ecosystem models can have strong nonlinear behavior. Moreover biogeochemi-cal dynamics often introduce nondifferentiability into the system, and for such a challenging environment obtaining an adjoint becomes nontrivial. Also the model resolutions are much finer than in the past, which results in a large dimension of the model states. This introduces a limitation for using the finite difference gra-dient approximations, since these are not suitable for large problems. A number of methods have been proposed to deal with such challenging configurations by ob-taining the adjoint in a reduced space (Vermeulen and Heemink, 2006; Cao et al., 2007; Fang et al., 2009). Although the methods differ in their approach, all are ba-sed on the proper orthogonal decomposition (POD), also known as the Karhune-Lo´eve transform (KLT), principal component analysis (PCA) or the method of empirical orthogonal functions (EOF) (Pearson, 1901; Shlens, 2009). None of the

Introduction 5 model reduced schemes have yet been used in ecological applications.

Based on a number of simulations of the original model, proper orthogonal decomposition is used to obtain a reduced model. The model-reduced 4D-Var is performed in the reduced space, therefore, the implementation of the adjoint of the tangent linear approximation of the original model is not required. Instead, it is approximated by the adjoint of the tangent linear approximation of the reduced model. The method is easily extended to the control of the initial condition, hence the parameter calibration is coupled together with the initial condition estimation. Due to the limitations of the adjoint model development for ecosystem models, such approaches may serve as an attractive tool for these applications. There-fore, the aim of this work is to evaluate the feasibility of a reduced adjoint ap-proach in calibrating ecosystem models. The model-reduced 4D-Var proposed by Vermeulen and Heemink (2006) was effectively used in several applications, such as groundwater flow (Vermeulen et al., 2005), shallow-water flow (Altaf et al., 2009, 2010), oil reservoir optimization (Kaleta et al., 2011), and morphodynamics (Garcia et al., 2013). An advantage was shown especially for models characteri-zed by periodic behavior (Altaf et al., 2009), since for these type of models, the number of required model simulations to obtain the reduced model is relatively small. Due to seasonality, ecosystem models show periodic behavior, therefore this method also serves as a potential tool for ecological applications.

1.3 Main objective

The main objective of this work is to explore and evaluate the feasibility of the model reduced four dimensional variational data assimilation method (4D-Var) for application to ecosystem models. It can be further divided into the following challenges which are addressed in the thesis:

Adjoint approximation. Typically, ecosystem models have strong nonli-near behavior, and consist of nondifferentiable biogeochemical dynamics. Such systems are extremely challenging in terms of obtaining their adjoints, whereas the standard finite difference technique used for the gradient approximations is no longer an efficient solution for the currently very large systems. The model re-duced 4D-Var method is used to build an approximation of the ecosystem adjoint in the reduced space, so as to maintain the accuracy similar to the full space, as well as keep the computational cost feasible.

Parameter calibration. The available data for the ecological models mostly consist of chlorophyll-a measurements. In the model chlorophyll-a concen-tration is an output resulting from the total performance of biological substances, which are governed by many parameters. The objective is to estimate the para-meter values based on the chlorophyll-a observations. However, there are many parameters that have a similar impact at the green pigment concentration. Even if the parameters have different functions in the model, they still can result in nearly the same contribution to the chlorophyll-a value. Such parameters are

very difficult to estimate when only chlorophyll-a measurements are given. Initial condition estimation. The initial conditions for the biological concentrations are of the same size as the model states, which especially in large scale ecosystems may become extremely large. In this work the size of the initial condition of the large scale ecosystem is of the order 105_{. This increases the}

size of the control vector dramatically, and therefore makes its estimation more demanding. An additional challenge of the initial state control is maintaining a physically plausible structure by using appropriate background error covariance strategies. In this work the initial condition is estimated in the reduced space, which decreases the size of the control vector substantially, while maintaining high accuracy of the updates. Moreover, operating in the reduced space utilizes the already available adjoint approximation for the initial condition, and the design of the method maintains the physical structure of the initial condition.

1.4 Overview of the thesis

The thesis is organized as a collection of published work preceded by a general introduction and an extended description of the main ecological model used in this work, and followed with general conclusions. Therefore, sometimes there might be a slight overlap of the presented material, since the goal was to keep the publications as close to their original form as possible.

The primary production model targeted in this research is the 2D North Sea coupled physical-ecological model BLOOM/GEM, therefore a separate chapter is devoted for an extended description of this model (see Chapter 2). To get a better understanding of the model reduced 4D-Var method and its behavior in a marine ecosystem application, the method is first described and applied to a simple 1D ecological model. The method is applied in a twin experiment frame-work, where synthetically generated surface phytoplankton data is used. Three parameters are calibrated in combination with estimation of the initial conditions. Different control strategies are explored, showing that updating the initial condi-tion is essential to obtain accurate parameter calibracondi-tions, which motivates the combined estimation approach in the further applications. The results are presen-ted in Chapter 3. Although the system dimension used to investigate the method is small, it is representative of yearly cycles of a phytoplankton bloom, it iden-tifies important issues that can arise, when applying data assimilation methods in ecosystem models, before integration in real complex configurations. Based on the relatively good results obtained in the 1D ecological model, the next step is to apply the model reduced 4D-Var method and assess its performance in a realistic large scale ecosystem. The technique is implemented for the 2D North Sea cou-pled physical-ecological model and presented in Chapter 4. The main goal is to improve the model predictions, with the focus on the chlorophyll-a concentration, which is achieved by biological parameter calibration. Based on sensitivity ana-lysis, supported by expert opinions, a number of parameters were selected as the

Introduction 7 most significant in the BLOOM/GEM model. A number of twin experiments were performed aiming at the calibration of two leading parameters from the sensiti-vity list. Similar control strategies as in the simple ecosystem case were applied, confirming that initial condition estimation enhances the parameter calibration in terms of accuracy and computational efficiency. Moreover, the results show that the method is capable of significantly improving the chlorophyll-a descrip-tion. The relatively good performance of the model reduced 4D-Var method in the application to the 2D North Sea BLOOM/GEM model demonstrates its po-tential as a calibration tool to enhance model predictions using real chlorophyll-a measurements derived from the remotely sensed MERIS data, which is examined in Chapter 5. Following the same control vector strategies as for the twin ex-periments, the parameters and initial conditions are estimated. The predictions of chlorophyll-a concentration resulting from the assimilated model are validated against remotely sensed MERIS measurements during a nearly two year period after the assimilation. The model performance is enhanced with respect to the un-assimilated model, showing for some control strategies up to 10% improvement of the model compared to the validation data, which demonstrates that the model reduced approach is a useful tool for improving the chlorophyll-a model predic-tions. Finally, the conclusions and recommendations are drawn in Chapter 6.

Chapter 2

Primary production model

BLOOM/GEM

Phytoplankton blooms are important factors in ecological and water quality mo-deling. In order to model algal blooms, algal biomass is represented by the concen-tration of chlorophyll-a. To model chlorophyll-a concenconcen-tration, primary produc-tion, the nutrient cycle, as well as the phytoplankton species composiproduc-tion, the generic ecological model (BLOOM/GEM) (Blauw et al., 2009; Los et al., 2008) is used. The model was developed such that it can be applied to any aquatic system, and consists of detailed underlying hydrodynamics, suspended sediment and river loads, which are required for ecological modeling. A 2D application of this model to the southern North Sea is used in this work.

BLOOM/GEM is a part of the Delft3D modeling system of Deltares, which consists of modules targeted at different processes. The main module is the hydro-dynamical model, which calculates the advective-dispersive transport (Delft3D-FLOW as described in Lesser et al. (2004)). Possible additional modules are for waves, morphology, and suspended sediments, although these are optional and are not used in the model setup in this work.

The transport of the substances is commonly described by the advection-dispersion-reaction equations, i.e.

∂C ∂t = −u ∂C ∂x − v ∂C ∂y − w ∂C ∂z + ∂ ∂x µ Dx∂C ∂x ¶ + ∂ ∂y µ Dy∂C ∂y ¶ + ∂ ∂z µ Dz∂C ∂z ¶ + S(x, y, z) + P (x, y, z) (2.1) where C denotes the concentration (kg m−3_{), u, v, w are the velocity components}

(m s−1_{), D}

x, Dy, Dzare the components of the dispersion tensor (m2s−1), x, y, z

are the coordinates of the three spatial dimensions (m), and t denotes time. S denotes source or sink of mass due to waste loads and boundaries (kg m−3 _s−1_),

and P denotes source or sink of mass due to physical, biochemical and biological processes (kg m−3 _s−1_).

Equation (2.1) is a basis for substance transport in the BLOOM/GEM model. Hence, the concentration change in time is caused by: the advective and dispersive transport, the mass change due to waste loads and boundaries (S), and due to physical, chemical and biological reactions, which convert one substance into another and/or make it appear or disappear (P ). The latter are also referred to as the water quality and ecological processes, and are described in the following section.

2.1 Water quality and ecological processes

Part of equations (2.1) are the physical, chemical and biological processes (P ), responsible for the reactions between the model substances. They may increase or decrease the concentration of a substance, or they may transform it into another substance. The main reactions affect algae growth and mortality,

mineraliza-Algae P N C N NH4-N NO3-N P PO4-P Detritus P N C settling settling respiration photosynthesis Nutrient mineralisation mineralisation metabolism mortality DO production consumption reaeration Detritus in Sediment C N P Si Si Si N2 denitrification

mineralisation & nitrification

autolysis Si consumption nitrification Grazers grazing grazing oxygen consumption biodeposition AIP adsorption Microphytobenthos C N P Si AIP in sediment settling mortality photosynthesis

Figure 2.1: Schematic overview of all state variables and processes in GEM model.

The variables in grey and the processed indicated by dashed lines have not been used in the North Sea application. DO denotes dissolved oxygen, and AIP stands for adsorbed inorganic phosphorus. (Los et al., 2008; Blauw et al., 2009)

Primary production model bloom/gem 11 tion of organic matter, nutrient uptake and release, and oxygen production and consumption. All the possible processes in BLOOM/GEM are illustrated sche-matically in Figure 2.1.

There are three nutrients considered in the nutrient cycle in the BLOOM/GEM model: nitrogen (N ), phosphorus (P ) and silicate (Si). They appear in the mo-del as dissolved inorganic concentrations of ammonium (N H₄+), nitrate (N O₃−), ortho-phosphate (P O4) and dissolved silicate (Si). The carbon (C) cycle is

mo-deled only partially by BLOOM/GEM. The model controls the mass balance of organic carbon, whereas the dissolved inorganic carbon is assumed to be unlimi-ted in the system. See Figure 2.2 for an example of model simulation presenting the annual mean of the nutrients.

Ammonium: mean 2003 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 (a) Ammonium Nitrate: mean 2003 0 0.1 0.2 0.3 0.4 0.5 0.6 (b) Nitrate Phosphate: mean 2003 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 (c) Phosphate Silica: mean 2003 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 (d) Silica

Figure 2.2: Example of the 2D BLOOM/GEM model simulation showing the annual

mean for the nutrient concentrations.

.

Dead particulate organic matter is modeled in BLOOM/GEM as detritus sepa-rately for water and sediment, where fractions for different nutrients are available (DetC, DetN, DetP, DetSi, DetSedC, DetSedN, DetSedP, DetSedSi).

The formulation of the model allows for modeling many different algae species (including different species for green and blue-green algae (cyanobacteria), as well as macroalgae). In this work the application of BLOOM/GEM includes four phytoplankton species: diatoms, flagellates, dinoflagellates, and P haeocystis. See

Figure 2.3 and 2.5 for a distribution of species within a year 2003. J F M A M J J A S O N D J 0 1 2 3 4 5 Time Chlorophyll−a in algae

Mean over the North Sea domain

Diatoms Phaeocystis Flagelattes Dinoflagellates

Figure 2.3: Example of model simulation showing the chlorophyll-a concentration per

alga species averaged over the North Sea domain and shown within year 2003 [mg-Chlfa/m3].

2.1.1 Competition and adaptation

The core of the model formulation is the competition between algae species, as well as the modeling of the adaptation of phytoplankton to different nutrient and light conditions (Los et al., 1984; Los and Brinkman, 1988).

Most phytoplankton species adapt rapidly to the changing environmental conditions (Harris, 1986). In order to model the adaptation of species to dif-ferent conditions in BLOOM/GEM, each species is divided into three phenotype groups:

• energy type (E-phenotype): species are modeled as E-phenotype, when they

are subjected only to the light limitation. Hence, these are the fastest growing types (high growth rates, low mortality rates), and they have high N/C and P/C ratios, since they appear when none of the nutrients are limiting.

• nitrogen type (N-phenotype): species are modeled as N-phenotypes when

nitrogen becomes a limiting factor, hence these types have lower N/C ratios. Due to limiting nutrient condition they have lower maximum growth rates, and higher mortality rates. N-phenotypes are also characterized by higher settling velocities, and lower chlorophyll-a content.

Primary production model bloom/gem 13

• phosphorus type (P-phenotype): species are modeled as P-phenotypes when

phosphorus becomes a limiting factor, hence these types have lower P/C ratios. Due to limiting nutrient conditions they have lower maximum growth rates, and higher mortality rates. P-phenotypes are also characterized by higher settling velocities, and lower chlorophyll-a content.

Each of the phenotypes is modeled as a separate state variable. When the light and/or nutrient availability change, a species can be instantaneously switched from one phenotype to another within the same alga species group, such that it is always optimally adapted to the current conditions. The model allows for coexistence of different phenotypes. In Figure 2.4 three phenotypes of diatoms are shown as an example.

Diatoms E−phenotype: MAR−APR−MAY 2003

0 0.01 0.02 0.03 0.04 0.05 (a) E-phenotype

Diatoms N−phenotype: MAR−APR−MAY 2003

0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 (b) N-phenotype

Diatoms P−phenotype: MAR−APR−MAY 2003

0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 (c) P-phenotype Figure 2.4: Example of the 2D BLOOM/GEM model simulation showing seasonal

mean concentration for different diatom phenotypes; the mean was calculated for a period March–May 2003.

Competition between the phytoplankton species is the main process in the model (Los et al., 1984; Los and Brinkman, 1988; Los and Wijsman, 2007). The change in species composition resulting from competition is slower from the adap-tation processes, and it depends on algal growth and mortality rates. Different species groups coexist together, which involves changing dominance of the spe-cies. BLOOM/GEM assumes that fast growing species become dominant when resources (light and nutrients) are abundant, whereas the slow growing (efficient) species dominate when the resources are limited. This is based on the theory of K-and r- strategies (Harris, 1986; Reynolds et al., 1983), according to which the spe-cific algae requirements for nutrient uptake together with their growth rates are essential in determining the species composition resulting from competition. For example shade favoring species, like cyanobacteria will often dominate in shaded areas, whereas in the areas with higher light intensity there will be more light favo-ring species such as diatoms and green algae (Los, 2009). By this principle each phytoplankton species maximizes its own benefit, and it is equivalent to maxi-mizing the total net production of the all phytoplankton groups (Blauw et al., 2009). Therefore, an optimization technique can be used to calculate the algal

biomass by maximizing the total net production. The BLOOM/GEM model uses computationally efficient linear programming methods to do this, as described in Danzig (1963). Next, the model finds a phytoplankton species composition by maximizing the total net production with specific constraints on their growth, mortality, light extinction, and nutrient availability (Blauw et al., 2009).

2.1.2 Phytoplankton processes

The processes which relate to phytoplankton in BLOOM/GEM, next to the com-petition and adaptation, are primary production, respiration and mortality. Al-gae capture energy from the sunlight and through the process of photosynthesis convert dissolved organic matter into organic compounds and dissolved oxygen. Nitrogen uptake is preferably done through ammonium (N H+

4), and in the case

when this is depleted then via nitrate (N O−

3). Certain blue-green algae species

are able to fix nitrogen in its elementary form (N2) from the atmosphere, and

later they release it as ammonium into the water (WL | Delft Hydraulics, 2007). However, in this work no cyanobacteria species are included, so nitrogen fixation is not a part of the current setup of BLOOM/GEM.

Mortality process controls nutrients recycling when primary producers die. Namely, after algae die, a fraction of the phytoplankton biomass is released as dissolved inorganic matter ready for the uptake for other algae, which is done through the process of autolysis. The remaining fraction of the biomass is released into the water as detritus. Grazing is an additional loss process, which may be modeled as a separate module, included as a grazing function, or as a part of the total mortality (Los and Wijsman, 2007), as it is done in this work.

Settling process is modeled with settling velocity which is constant over time, independent from turbulence and bottom shear stress. Different settling velocities are specified for different algae species and phenotypes (Blauw et al., 2009).

Due to the time scale of a cell division, which is in the order of one day, typically the simulation time step for the phytoplankton processes in the BLOOM/GEM model is chosen to be 24 hours. The transport and a few of processes are simulated using a short time step, typically 30 minutes.

2.1.3 Nutrient cycling processes

In general the nutrient cycle processes model the uptake of nutrients by algae, which are further released back into the system in the form of dissolved inorganic nutrients and dead organic matter, referred to as detritus. The detritus in the water is reprocessed into the dissolved inorganic nutrients through the microbial decomposition (mineralization). Similarly the dead organic matter in the sedi-ment is reprocessed into inorganic nutrients, which are further released from the sediment into the water column (several formulations are available).

Ammonium is converted into nitrate through the process of nitrification

N H+

Primary production model bloom/gem 15

Chlf−a in Diatoms: MAR−APR−MAY 2003

0 0.5 1 1.5 2 2.5 3

(a) Diatoms MAM

Chlf−a in Diatoms: JUN−JUL−AUG 2003

0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 (b) Diatoms JJA

Chlf−a in Diatoms: SEP−OCT−NOV 2003

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 (c) Diatoms SON

Chlf−a in Phaeocystis: MAR−APR−MAY 2003

0 1 2 3 4 5 6 7 8 (d) Phaeocystis MAM

Chlf−a in Phaeocystis: JUN−JUL−AUG 2003

0 0.2 0.4 0.6 0.8 1

(e) Phaeocystis JJA

Chlf−a in Phaeocystis: SEP−OCT−NOV 2003

0 0.1 0.2 0.3 0.4 0.5 (f) Phaeocystis SON

Chlf−a in Flagellates: MAR−APR−MAY 2003

0 0.2 0.4 0.6 0.8 1 1.2 (g) Flagellates MAM

Chlf−a in Flagellates: JUN−JUL−AUG 2003

0 0.5 1 1.5 2 2.5 3 3.5 4 (h) Flagellates JJA

Chlf−a in Flagellates: SEP−OCT−NOV 2003

0 0.5 1 1.5 2 2.5 3

(i) Flagellates SON

Chlf−a in Dinoflagellates: MAR−APR−MAY 2003

0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 (j) Dinoflag. MAM

Chlf−a in Dinoflagellates: JUN−JUL−AUG 2003

1 2 3 4 5 6 (k) Dinoflag. JJA

Chlf−a in Dinoflagellates: SEP−OCT−NOV 2003

0 0.5 1 1.5 2 2.5 3 (l) Dinoflag. SON Figure 2.5: Example of the 2D BLOOM/GEM model simulation showing seasonal

means for chlorophyll-a concentration per alga species [mgChlfa/m3]

This usually happens in the sediment and deep water layers, as these are richer in dissolved oxygen. Whereas, in the areas where the concentration of dissolved oxygen drops below a specified level the rates for nitrification are reduced. In anaerobic regions, nitrate is subjected to denitirification, which processes it into nitrogen

N O−

3 + H+→ 0.5 N2−+ 1.25 O2+ 0.5 H2O (2.3)

Additionally, the elementary nitrogen may be fixed back into ammonium by blue algae or some specific bacteria species, and part of it is released into the atmos-phere in form of nitrogen gas.

Similarly, there are processes available for phosphorus, such as adsorption and desorption of ortho-phosphate (P O4), and sedimentation of adsorbed inorganic

phosphorus (AIP). The adsorption of phosphate is a physical-chemical process which depends on the water pH. Primary production may increase the pH, which causes desorption of the phosphate from the suspended sediment, which results in an increased availability of that nutrient for the primary producers.

Silicate is a nutrient used only by diatoms, and it is only available in the form of dissolved silicate. It is replenished by the slow dissolution of opal silicate, which comes from the residue of the silicate skeletons of diatoms.

2.1.4 Oxygen processes

There are several processes in the BLOOM/GEM model which relate to the dis-solved oxygen concentration. A large fraction of the oxygen is produced by the phytoplankton through the primary production processes. There are several pro-cesses which uptake the oxygen, such as algal respiration, mineralization of detri-tus (water and sediment), and nitrification (Los et al., 2008).

Depending on the saturation of the dissolved oxygen concentration in the wa-ter column, oxygen may enwa-ter from or escape to the atmosphere. Replenishment of the oxygen from the atmosphere is called reaertion. The atmosphere is an im-portant source of dissolved oxygen especially in the areas where the organic matter discharges from rivers (WL | Delft Hydraulics, 2007). These areas are particularly prone to oxygen depletion, since its large amount is consumed during the minera-lization process. The rate of the oxygen uptake through decomposition is higher than that of the reaertion from the atmosphere, and therefore the dissolved oxy-gen concentration will drop until the organic matter is sufficiently reduced. This phenomena creates the so called dissolved oxygen sags (WL | Delft Hydraulics, 2007).

2.1.5 Light extinction

Light availability is an essential factor for algal growth, therefore the correct evaluation of light conditions in the water is very important. The extinction of light is modeled as an exponential decrease of light intensity with depth based

Primary production model bloom/gem 17 on the Lambert-Beer formula. It is calculated in the model from the given light irradiation (photosynthetic active radiation - PAR) and the light extinction due to substances in the water:

Iz= Ioe−kz (2.4)

where Izis the underwater light intensity at depth z, Io is the surface irradiance,

k is the total light extinction coefficient in the water column, which is calculated

as a sum of the background extinction (due to water itself and non-modeled sub-stances) and extinction caused by suspended inorganic matter, humic substances (determined by salinity), detritus and self shading of phytoplankton. Each of the substances in the model have been assigned a specific extinction coefficient (Los and Bokhorst, 1997; Los et al., 2008).

2.1.6 Mineralization

Mineralization (decomposition) of particulate organic matter is the reverse pro-cess to photosynthesis, in which organic compounds and oxygen are transfor-med into carbon dioxide and water. There are two types of mineralization, one which involves decomposition of detritus resulting from algal mortality, and second which decomposes the organic matter originating from the waste loads (WL | Delft Hydraulics, 2007).

2.1.7 Sediment processes

Sediment is an important part of the water system. It is a living environment, which can be a sink or a source of particulate or dissolved matter resulting from the interaction between the water column and the sediment. The organic and inorganic particulate matter are stored in the sediment as a result of sedimen-tation. Through the process of mineralization in the bottom the organic matter is recycled into nutrients and released back into the water column. Part of the detritus is resuspended back in the water column due to the erosion processes included in the model. Whereas, the part of the substances are moved from the active sediment layer into the deeper sediment layers through the burial process (WL | Delft Hydraulics, 2007; Los et al., 2008). As the deeper sediment layers are not modeled in BLOOM/GEM they are a sink for the nutrient cycle in the model. The burial rate is not generally known for marine ecosystems, and it is assumed to also include other unknown source/sinks in the mass balance (Blauw et al., 2009).

2.1.8 Grazing

A couple of formulations are available to model dynamically the zooplankton and zoobenthos in BLOOM/GEM. However, their simulation is a very complex process, and still does not provide better results than grazing introduced through a forcing function (WL | Delft Hydraulics, 2007). Therefore, typically the grazer

biomass is introduced to the model as a forcing, which is derived from a field data. Different types of grazer biomass can be included, which are further used to simulate the grazing effects by the model (Blauw et al., 2009). Alternatively, grazing may be also controlled by enhancing the mortality rates, as it is done for the North Sea applications of the model (Los et al., 2008). As the mortality rates for the E-phenotypes are lower than those of the N- and P-phenotypes, the mortality effect is enhanced at the peak of the spring bloom, when the phenotypes switch from the energy limited to the efficient phenotypes. This happens at the same time when the zooplankton grazing is expected to be dominant as well, and therefore, enhancing the mortality formulation accounts for grazing in the model (Los et al., 2008).

2.2 North Sea application

Development of a model to describe and predict the ecosystem of the Dutch coastal zone began in the early 1980s (Van Pagee et al., 1988; Moll and Radach, 2003; Radach and Moll, 2006; Lacroix et al., 2007). The work was ongoing at various groups and institutes, which focused on different geographical regions and modeling objectives. The BLOOM/GEM North Sea application is one of several existing large-scale 3D ecosystem models applied to the North Sea area. It is a result of many years of experience, and it has been repeatedly calibrated and validated. It is a result of many years of experience, and it has been repeatedly calibrated and validated. Many in situ observations obtained from the monitoring stations, with the majority from the Dutch coastal area, were used. It included utilizing measurements of many different substances from many years, aiming to globally optimize the model performance in as far as monitoring stations are representative (Los et al., 2008).

2.2.1 Schematization

The considered region of interest starts at the English Channel and continues through the southern part of the North Sea, which stretches around 400 km north of the Wadden Islands reaching up to 57◦_{N, including regions such as Dogger}

Bank, Oyster Ground, Southern and German Bight, see Figure 2.6.

The schematization used in the BLOOM/GEM North Sea application in this work is called ”coarse ZUNO-grid”, see Figure 2.6. It uses a curve-linear model grid in order to represent the Dutch rounded coastline, and to enable the transport along the shoreline, which is typical for this area. To assure high accuracy in the coastal areas of interest, a relatively high resolution is used in that area, i.e. approximately 1 × 1 km in the Dutch coastal zone. Whereas, the northern part of the domain is much coarser, i.e. approximately 20 × 20 km. This results in a horizontal dimension of 65 × 134, which is used for hydrodynamic simulations, whereas 4350 active horizontal grid cells are used for ecological simulations. There

Primary production model bloom/gem 19

Figure 2.6: Model grid and the map of the southern North Sea used in the hydrodynamic

and BLOOM/GEM simulations. Left: the map of the southern North Sea including its bathymetry. Right: model grid known as ”coarse ZUNO grid”.

are 10 sigma layers in the vertical direction, which are finer near the bed and near the surface (around 4% of the local water depth) and coarser in the middle of the water column (around 20% of the local water depth). In this work 2-dimensional (vertically mixed) version of the model is used for the ecological application, which can be justified by generally well mixed conditions in the Dutch coastal zone, i.e. the main area of interest. The hydrodynamics are calculated using full 3D schematization, which are later aggregated to 2D and used as off-line input for the 2D version of the ecological model (Los and Blaas, 2010).

2.2.2 Hydrodynamics

The hydrodynamical model (Delft3D-Flow) used to calculate the advective dis-persive transport was extensively calibrated and validated before it was coupled with the BLOOM/GEM module (Lesser et al., 2004; Los et al., 2008).

As the water temperature varies considerably over the modeled area, a 3D temperature model (part of Delft3D-Flow) output is used as a temperature forcing for the ecological simulations. The wind conditions are also quite important in the area of interest, as these have a large influence on the long-shore and cross-shore currents in the coastal areas. Therefore, the wind as well as other meteorological conditions are taken into account by forcing the hydrodynamical model with the wind and pressure fields (Proctor et al., 2002; Delhez et al., 2004; Los et al., 2008). Wind speed, day length, and solar irradiance are adopted from historic measurements from a near shore station of Royal Dutch Meteorological Institute, and their values are extrapolated to be adequate to the areas of interest of the North Sea.

The exchange of water masses through the open boundaries of the English Channel and the northern boundary at 57◦_{N are calculated using a large scale}

hy-drodynamical model covering the entire Continental Shelf (Gerritsen and Bijlsma, 1988). The boundary conditions for the substance concentrations are assumed constant in time along the northern boundary, and the input through the English Channel is given on a monthly basis (based on the study of Laane et al. (1993)). There are multiple fresh water and nutrient sources in the model given by the river discharges. The main rivers of the Netherlands, Germany, France and UK are included. Light limitation is an important factor for primary production mo-deling, and it is highly affected by the inorganic suspended matter. Therefore, it is essential to include it in the model. Results from a separate module simulation of the inorganic suspended matter give a base for creating input for the ecological model (Los et al., 2008).

Chlf BACK MAR−APR−MAY 0 1 2 3 4 5 6 7 8 9 10 (a) Mar-Apr-May 2003 Chlf BACK JUN−JUL−AUG 2 4 6 8 10 12 (b) Jun-Jul-Aug 2003 Chlf BACK SEP−OCT−NOV 1 2 3 4 5 6 7 (c) Sep-Oct-Nov 2003 Figure 2.7: Example of the 2D BLOOM/GEM model simulation showing seasonal

means of chlorophyll-a concentration.

2.2.3 North Sea specific processes

In the main areas of interest of the North Sea most of the mineralization of organic matter happens in the water column. Therefore, a simple benthic formulation of the mineralization is used in this application. In the shallow areas, such as the Wadden Sea, the mineralization rates for the sediment are larger, and in some areas may even exceed those in the water. As these shallow areas are not entirely accounted for in the simplified formulation, it is partly the reason why the model results at these areas are not expected to be as accurate as those in the areas of interest.

Chapter 3

Application of model reduced

4D-Var to a 1D ecosystem model

∗

The model reduced 4D-Var (Vermeulen and Heemink, 2006) is investigated to test its feasibility in ecosystem application. Ecological models are known for their high nonlinearity and issues with non-differentiability. However, since the method is performed in the reduced space, the implementation of the adjoint of the tangent linear approximation of the original model is not required.

Twin experiments are conducted in a 1D ecological model. Surface phytoplank-ton data are used, with 30% log-normally distributed measurement error. Three parameters are chosen for calibration. The method performs very well in the setup where a perfect initial condition is used, as well as for the combined parameter and initial condition estimation. A relatively well calibrated initial condition contri-butes to accurate parameter estimations. Not accounting for the wrongly assigned initial condition leads to incorrectly calibrated parameters.

3.1 Introduction

In the presence of environmental issues caused by climate change and eutrophi-cation (Nixon, 1995; Hallegraeff, 1993, 2009, 2010; Pauly et al., 2002; Gregg et al., 2003, 2005; Peperzak, 2003; Paerl and Huisman, 2008; Brush, 2008; Schindler et al., 2008), ecosystem models capable to provide accurate predictions are of high inter-est. Despite many modeling challenges (Doney, 1999; Jørgensen, 2008), ecological models have evolved to be considered as relevant predictive tools. However, most sophisticated systems are not able to reproduce reality completely, therefore it is essential to integrate them with the available observations. Data assimilation

∗_{Appeared in revised form in Ocean Modelling: Pelc, J. S., Simon, E., Bertino, L., El Serafy,}

G., and Heemink, A. W. (2012). Application of model reduced 4D-Var to a 1D ecosystem model.

Ocean Modelling, 57-58:43-58.