J. J. Vallino Ecosystems Center

Improving marine ecosystem models:

use of data assimilation and mesocosm experiments

J. J. Vallino

Marine Biological Laboratory

Woods Hole, MA 02543

Submitted

12 November, 1998

Revised:

Keywords: data assimilation, parameter estimation, mesocosm, food web model, adjoint method, optimization.

Abstract

Our inability to accurately model marine food webs severely limits the prognostic capabilities of current generation marine biogeochemistry models. To address this problem we examine the use of data assimilation and mesocosm experiments to facilitate the development of food web models. The components of the data assimilation demonstrated include the construction of measurement models, the adjoint technique to obtain gradient information on the objective function, the use of parameter constraints, incorporation of discrete measurements and assessing parameter observability. We also examine the effectiveness of classic and contemporary optimization routines used in data assimilation.

A standard compartment-type food web model is employed with an emphasis on organic matter production and consumption. Mesocosm experiments designed to examine the interaction of inorganic nitrogen with organic matter provide the data used to constrain the model. Although we are able to obtain reasonable fits between the mesocosm data and food web model, the model lacks the robustness to be applicable across trophic gradients, such as those occurring in coastal environments. The robustness problem is due to inherent structural problems that render the model extremely sensitive to parameter values. Furthermore, parameters governing actual ecosystems are not constants, but rather vary as a function of environmental conditions and species abundance, which increases the sensitivity problem. We conclude by briefly discussing possible improvements in food web models and the need for rigorous comparisons between models and data (a modeling workbench) so that performance of competing models can be assessed. Such a workbench should facilitate systematic improvements in prognostic marine food web models.

1. Introduction

Aquatic food web modeling has primarily focused on conceptual models where the objective is to capture, qualitatively, the behavior of aquatic ecosystems. For example, models have been designed to examine the development of phytoplankton blooms associated with the shoaling of the mixed layer (Evans and Parslow, 1985; Sverdrup, 1953)§§, the importance of the microbial loop (Taylor and Joint, 1990)§, the chaotic behavior of food webs (Beckers and Nihoul, 1995)§, and the “paradox of the phytoplankton” (Hutchinson, 1961; Stone, 1990)§§. Although conceptual models have proved to be valuable for elucidating principles governing aquatic ecosystems, we must consider them qualitative because they are often only compared to limited experimental data sets, observations, or to no data at all. These models often lack the breadth to fully capture ecosystem biogeochemistry. Consequently, there has been relatively little improvement in the quantitative, predictive capabilities of aquatic food web models, even though these models have been under development for at least 50 years (Totterdell, 1993)§. With improvements in transport and ocean circulation models and the growing interests to assess impacts of global change, there is an increasing need to develop quantitative ecosystem models where the focus is not on concepts but rather on prognostic capabilities, efficiency and accuracy (Sarmiento et al., 1993)§.

The difficulties in quantitative modeling of aquatic foods are well appreciated (Evans and Fasham, 1993a; Platt et al., 1981)§§, and can be categorized by three main challenges: 1) understanding and modeling of 4D complex multidisciplinary dynamics; 2) 4D data acquisition; and 3) rigorous model-data comparisons and data assimilation. Although complete ecosystem models must meet all three challenges, it is possible to decouple the modeling effort into independent components in order to facilitate overall model development. As will be demonstrated, mesocosm experiments (i.e., enclosed, experimental, aquatic ecosystems, see Grice et al. (1982)§) can facilitate this decoupling.

Marine ecosystems are complex systems that are governed by nonlinear growth behavior of the constituent organisms, by highly dynamic predator-prey interactions between organisms, and by hydrodynamic and other external drivers that can exhibit complex patterns themselves. Hence, developing models that can accurately predict the concentrations of organisms and the rates of the biogeochemical transformations they mediate over large scales in open systems is a non-trivial task (Challenge 1). However, transport processes can be removed and decoupled form the biogeochemical processes by the use of mesocosm experiments because these systems can be designed to be well-mixed. By reducing the system dimensions from 4D (i.e., space and time) to 1D (time only), greater effort can be placed on improving the ecosystem model, and this reduction also lowers computational overhead.

The large effort required to obtain comprehensive data sets is unavoidable in biological oceanography (Challenge 2). This challenge is even greater in systems with complex circulation, because it is difficult to close elemental balances without either detailed flux measurements across system boundaries or a well-calibrated transport model. Indeed, often both are necessary. Although there is no substitute for field observations, again mesocosm experiments can greatly facilitate problem solution, as these systems can be intensely sampled, and only boundary fluxes across the air-water interface need be measured or accounted for. Mesocosms also allow manipulation of experimental conditions to investigate system dynamics in various regions of state space. Interestingly, while mesocosms have been used to experimentally study ecosystem processes and impact of pollutants (Grice and Reeve, 1982)§, only a few cases exist where they have been utilized for model development (Baretta-Bekker et al., 1998)§. In this manuscript, we will rely on mesocosm experiments as the source for observations. Hence, we will investigate only 1D data assimilation.

Even when observations are available, it is difficult to employ them for model improvement. As a result, aquatic food web models go largely untested and uncalibrated. Challenge 3, which addresses this difficulty, can be broken into two sub-areas: a) integration of data and models (data assimilation) to improve model state estimates, model parameter estimates, or both; b) adaptive improvements of observation strategies (which data to sample, which sensors and platforms to use, and when and where to sample) based on data and on integrated data/model simulation.

A general definition of data assimilation is the integration of models with data to improve the estimation of a system's state. The emphasis of data assimilation can be either data centric, or model centric (McLaughlin, 1995)§. In the fields of oceanography and meteorology, the former is often the case, while the latter is often called parameter estimation. In oceanography and meteorology, data are often abundant and the models are accurate as they are derived from first principles (conservation of mass, energy and momentum). However, because of sensitivity to initial conditions (i.e., chaos) as well as inaccurate model parameterizations, these models require frequent updates to prevent model divergence from the true state. In this case, data which are more or less randomly distributed in time and space, are interpolated to a model grid, and compared to a model forecast. By appropriate weighting of these two state assessments, an improved estimate of the true state can be generated and used as an initial condition for the next model forecast (Bergamasco et al., 1993; Daley, 1991; Evensen, 1994; Tziperman and Thacker, 1989; Wiggins, 1972)§§§§§. The Kalman filter is a classic example of model-data integration in this manner (Burger et al., 1998)§, but other methods exist, including adjoint method, inverse analysis, and optimal interpolation (Courtier et al., 1993; Robinson et al., 1998)§§. In model-centric data assimilation, data are used to improve estimates of model parameters instead of the model state (Bennett, 1992; Evensen et al., 1998; Gunson et al., 1999; Robinson et al., 1998)§§§§, although it is possible to achieve both objectives of state and parameter estimation simultaneously (Malanotte-Rizzoli and Tziperman, 1996)§. This manuscript addresses the parameter estimation component of data assimilation for marine food web models.

Minimizing parameter uncertainty is particularly important in food web models. Because marine food web models are quite sensitive to parameter values, model dynamics can change substantially with slightly different parameter values. Consequently, a poor fit between model output and observations can result from either poor model structure due to inappropriate physics or growth models, or from poor selection of model parameters. This ambiguity between structural uncertainties and parameter uncertainties makes it difficult to distinguish superiority between two or more competing models, and partially explains the plethora of models describing similar, if not the same, phenomena, such as phytoplankton growth (Geider and Osborne, 1992, pg. 159)§, (Fasham et al., 1990; Taylor and Joint, 1990)§§ zooplankton growth (Fasham et al., 1990; Franks et al., 1986; Gunson et al., 1999)§§§ or mortality closure terms (Steele and Henderson, 1992)§. If parameter uncertainty can be removed by assimilation of experimental data, then models with superior functional structure can be identified, and selected for improvement systematically. Of course, we need to define a standardized experimental data set that models can be tested against for this approach to be successful (a modeling workbench (Evans and Fasham, 1993b)§).

In addition to recursive model development, integration of measurements with models must also be considered (Challenge 3b). State variables used in models often do not correspond to quantities that are directly measurable in the field, which makes it difficult to use all available observations and to compare models to observations. Some measurements may encompass several state variables, while others may only represent a portion of a state variable. Particulate organic carbon (POC) is a frequently and easily obtained measurement. Yet, POC is typically not a state variable in most models because it represents carbon from both living organisms and detritus in various states of degradation. This type of model-observation mismatch is not often a significant problem, as we will show that it is usually possible to construct measurement models that relate state variables to measured variables (Robinson et al., 1998)§. However, mapping state variables to partial observations of that variable can be problematic. For example, biomass measurements may be available for a particular species, while the model may use an aggregated compartment representing a large number of species. In this case, no obvious mapping exists between the aggregated state variable and the observed quantity without increasing the complexity of the model. These types of problems are best resolved by better integration of experimental design with model development.

In this manuscript, we use data assimilation to find the optimum set of parameter values that minimizes error between model output and observations. For observations we make use of experimental data collected from marine mesocosms. The manuscript illustrates 1) the use of mesocosm experiments to generate data for model development, 2) data assimilation using several different optimization algorithms with mesocosm data and an organic-matter-based food web model, and 3) problems with the structure of current generation food web models. Although there have been several studies focused on parameter estimation in marine food web models, the majority have used only simulated data (Crispi and Mosetti, 1993; Gunson et al., 1999; Ishizaka, 1993; Lawson et al., 1995; Lawson et al., 1996; Marsili-Libelli, 1992)§§§§§§ or very limited data (Fasham and Evans, 1995; Marcos and Payre, 1988; Matear, 1995; Prunet et al., 1996)§§§§ to constrain the model. This manuscript is unique in that it is a rigorous model-data comparison. The Appendix contains the food web model equations and descriptions of the optimization routines.
2. Data Assimilation

A standard state space model is used to describe the biogeochemistry and food web dynamics of the marine ecosystem. The state variables represent those quantities of the system that are dynamic and best approximate the system in the context of the modeling objectives. The state space model in vector form is given by,

A nonlinear measurement model (Robinson et al., 1998)§,

Integration of the state model (Eq. 1) generates a prediction of the state variables over time for a given set of initial conditions, x0, and model parameters, p. The deterministic measurement model (Eq. 3) translates the predicted state variables to observation space so that the model output can be compared to observations, as given by the residual vector, µ §. The objective of the present data assimilation is to find the set of initial conditions and model parameters that minimizes some aspect of the residual vector, typically the sum of the squared residuals (i.e., least squares), but other objective functions are possible (Janssen and Heuberger, 1995)§. Formally, we seek here the set of parameters, k, that minimize the objective (or cost) function,

An optimum, k*, of a scalar function occurs where the gradient of the objective function in parameter space vanishes, µ § (or µ §), and the Hessian matrix of J(k*) at this point, µ §, is either positive definite (a minimum) or negative definite (a maximum). However, because the hyper-surface defined by the nonlinear function J(k) need not be quadratic, it is possible to have more than one optimum parameter set. Considering only minima of J(k), a minimum is said to be a global minimum if no other lesser or equal minima exist in or on the domain of interest; otherwise, it is considered a local minimum. Clearly, it is desirable to find the global minimum of J(k), but this is a non-trivial problem because the global minimum is only a local property of the objective function. Consequently, it is often necessary to exhaustively search all parameter space to locate the global minimum, which is usually quite costly. Consequently, development of efficient global optimization algorithms is an area of intense research (Barhen et al., 1997)§.

Optimization routines fall into two categories: those that search for local optima and those that search for global optima. The principal numerical techniques for local optimization are the methods of steepest descent, simplex, and conjugate gradient (Press et al., 1986)§ (see the Appendix). The majority of local optimization methods for smooth objective functions employ variations of the conjugate gradient technique (Fletcher and Reeves, 1964)§. The basic idea of this technique is to choose a set of conjugate search directions such that a minimization conducted along one search direction does not corrupt a minimization conducted along a previous search direction. Corruption of the previous iterates is why methods like steepest descent are inefficient (Press et al., 1986)§. For a quadratic objective function,
µ §, (5)
the optimum is found in n conjugate line-searches for an n-dimension system. Of course, the objective function defined by Eq. 4 is not quadratic. However, in a neighborhood sufficiently close to an optimum, µ §, the objective function can be locally approximated, via a Taylor series expansion, by the quadratic function given by Eq. 5, where the matrix A and vector b are the Hessian and gradient of J(k), and the constant c is the objective function, all evaluated around the current search point, µ §, and with µ § replaced by µ §. The details of most conjugate gradient methods differ in how the Hessian and the gradient of the objective function are constructed and stored. Powel's method attempts to construct the conjugate search directions solely from evaluations of the objective function. Fletcher-Reeves and quasi-Newton methods construct conjugate search directions by evaluating both the objective function and its gradient. Others require information on the Hessian at any given point (Chow et al., 1994)§.

Since knowledge of µ § can greatly facilitate location of optima, it is often desirable to calculate µ § for those optimization routines that can use it. However, since an analytical solution for Eq. 1 typically does not exist, the objective function given by Eq. 4 only implicitly depends on k, so that µ § must be determined numerically. One approach to calculate µ § is known as the sensitivity method (see the Appendix), but this approach was found to be too computationally intensive.

A more elegant technique to compute µ §, based on variational calculus and extensively employed in optimal control theory (Kirk, 1970)§, is the adjoint method (Courtier et al., 1993; Crispi and Mosetti, 1993; Marcos and Payre, 1988)§§§. In this method, a set of m adjoint variables, a(t) (also know as costate or Lagrange multipliers), are introduced, which allow the state space model constraints (Eq. 1) to be directly incorporated into an augmented objective function,µ §
µ § (6)
The solution of the augmented objective function can then be found by variational calculus, taking derivatives with respect to p, x0, and a. Using Eq. 6 and conveniently defining the evolution of the adjoint variables by the following adjoint equations and boundary conditions,

µ § (7)
the gradient of µ § is given by:

4. Implementation

In this section we describe application of nonlinear optimization for estimating parameters in a fairly complex aquatic food web model based on experimental data obtained from a mesocosm experiment.
a. Experimental setup and data

Estuaries and coastal zones receive a significant input of organic material exported from terrestrial ecosystems (Peterson et al., 1995; Smith and Hollibaugh, 1993; Smith and Mackenzie, 1987)§§§. In order to examine how marine food web communities process and use this dissolved organic matter (DOM), a mesocosm experiment was conducted in Woods Hole, MA, USA in which 7 m3 seawater enclosures were augmented with combinations of dissolved inorganic nitrogen (DIN) and DOM. The enclosures were polyethylene bags equipped with flotation collars and structural hoops which were deployed in Great Harbor. The DOM in these experiments was prepared by leaching leaf litter in seawater. The mesocosm experiment consisted of four treatments: Bag A, control (no additions); Bag B, one-time addition of DOM at the start of the experiment resulting in an increase of 300 ƒÝM DOC (dissolved organic carbon); Bag C, daily additions of nitrate, phosphate, and silica equivalent to 5 ƒÝM N, 0.5 ƒÝM P, and 7 ƒÝM Si, respectively; Bag D, treatments B and C combined. All treatments also received NaH13CO3 as a C-tracer.

The experiment was conducted outdoors under ambient lighting (Fig. 1) from 10 Sep. to 30 Sep. 1994. Since the focus of this paper is on data assimilation techniques and aquatic food web models, experimental interpretation will not be presented here and only the following subset of the full suite of measurements will be used for data assimilation (Eq. 3): dissolved organic carbon (DOC, ƒÝM C); particulate organic carbon (POC, ƒÝM C); particulate organic nitrogen (PON, ƒÝM N); dissolved inorganic nitrogen (DIN, ƒÝM N); chlorophyll a (ƒÝg l-1); net primary productivity at a specified depth (NPP, ƒÝM C d-1); bacterial productivity (BP, ƒÝM C d-1); light extinction coefficient (K, m 1). Measurements were taken daily or every other day.

b. Food web model

The food web model (Eq. 1) used to capture the dynamics of the mesocosm experiment (Fig. 2, and see the Appendix) consists of the following ten state variables: autotrophs (A(t), ƒÝM C), heterotrophs (Z(t), ƒÝM C), dissolved inorganic nitrogen (N(t), ƒÝM N), dissolved labile organic C (OCL(t), ƒÝM C), dissolved labile organic N (ONL(t), ƒÝM N), dissolved refractory organic C (OCR(t), ƒÝM C), dissolved refractory organic N (ONR(t), ƒÝM N), detrital C (DC(t), ƒÝM C), detrital N (DN(t), ƒÝM N), and bacteria (B(t), ƒÝM C). In this model C and N are linked where the biotic compartments (i.e., A(t), Z(t), and B(t)) are assumed to have fixed C:N ratios, while the C:N ratios of the non-living organic compartments are free to vary. For processes, bacteria utilize labile organic matter pools and can either remineralize organic N or immobilize DIN depending on the C:N ratio of the labile DOM. The aggregated heterotrophs pool graze both the bacteria and autotrophs and remineralize N. The autotrophs compete with bacteria for DIN, and excrete both labile and refractory DOM. Decomposition of detrital material into DOM and decomposition of refractory DOM into labile DOM are governed by first order kinetics. The overall model is similar to other aggregated, Monod-based growth models (Evans and Parslow, 1985; Fasham et al., 1990; Kremer and Nixon, 1978; Moloney and Field, 1991; Moran et al., 1988; Pace et al., 1984)§§§§§§, except that a significant emphasis has been placed on modeling organic matter production and consumption, which is consistent with the goals of the mesocosm experiment (see the Appendix).