## Marine Biogeochemical Modelling and Data Assimilation

Richard J. Matear and E. Jones

Abstract The inclusion of biogeochemistry into the Global Ocean Data Assimilation Experiment systems represents an exciting opportunity that involves significant challenges. To help articulate these challenges we review marine biogeochemical modeling and the existing applications of biogeochemical data assimilation. The challenges of biogeochemical data assimilation stem from the large model errors associated with biogeochemical models, the computational demands of the global data assimilation systems, and the strong non-linearity between biogeochemical state variables. We use the ocean state estimation problem to outline an approach to adding biogeochemical data assimilation to the Global Ocean Data Assimilation Experiment systems. Our approach allows the biogeochemical model parameters to be spatially and temporally varying to enable the data assimilation system to track the observed biogeochemical fields. The approach is based on addressing the challenges of biogeochemical data assimilation to improve both the state estimation of the biogeochemical fields and the underlying biogeochemical model.

### 12.1 Introduction

Marine biogeochemical (BGC) modelling is a key approach to helping us understand the biochemical processes responsible for the transfer of nutrients and carbon between inorganic and organic pools. Quantifying these biochemical processes is essential to understanding carbon cycling in the ocean, the air-sea exchange of carbon, the impact of climate variability and change on marine ecosystems, and the link between ocean physics and ocean biology.

CSIRO Wealth from Oceans National Research Flagship, CSIRO Marine and Atmospheric Research, Hobart, Australia e-mail: [email protected]

A. Schiller, G. B. Brassington (eds.), Operational Oceanography in the 21st Century, 295

DOI 10.1007/978-94-007-0332-2_12, © Springer Science+Business Media B.V. 2011

Data assimilation represents a new and exciting tool to advance ocean biogeo-chemical studies. The coupling of physical and biogeochemical data assimilation is a natural evolution of the Global Ocean Data Assimilation Experiment (GODAE). Data assimilation fuses models with a diverse set of observations to provide a more consistent view of the physical and biological state of the ocean. Extending the GODAE effort to include BGC data assimilation presents an exciting opportunity to expand the researchers interested in using ocean data assimilation products and Brasseur et al. (2009) provides a summary of BGC applications of data assimilation. However, delivering these new BGC data assimilation products is not a trivial task. Here we will focus on the challenges of utilizing the existing physical ocean data assimilation systems to include BGC data assimilation. To help set the context of this discussion we will first briefly review the basic components of the biogeochemical models that could be incorporated into GODAE data assimilation systems. Second, we will briefly discuss previous BGC data assimilation effort. Finally, we discuss the challenges of extending the existing GODAE data assimilation systems to include biogeochemical models to deliver biogeochemical data products. These tasks will be computationally demanding since the GODAE data assimilation systems already required a substantial computing resources and adding BGC to these systems will only further increase these requirements. Therefore, computational constraints will also be a factor to achieving BGC data assimilation.

### 12.2 Biogeochemical Modelling

Marine biota play an important role in nutrient and carbon cycling in the ocean. Our interest in prognostic marine biogeochemical models arises from the need to better understand, quantify, and eventually predict the oceans role in the global carbon cycle and biological processes involved in the cycling of carbon in the oceans. Here, I will focus on BGC models that involve carbon and nutrient cycling, with a link to the lowest trophic levels of the marine ecosystem. This clearly is not the only application of BGC modelling and other applications include the prediction of harmful algal blooms (Franks 1997), the quantitative understanding of oceanic food webs up to fish (Brown et al. 2010), and the possible impact of marine sulfur emissions on the formation of cloud condensation nuclei (Gabric et al. 1998). To articulate what BGC modelling represents, we have chosen to discuss a simple model that describes the nitrogen cycling and its links to the lowest trophic levels of the marine ecosystem. The key processes to capture in the model is the biological uptake of nitrogen as well as the subsequent remineralization of the organic matter back into to an inorganic form. The example BGC model is based on the model that was applied to the macro-nutrient replete region south of Tasmania to investigate phytoplankton dynamics (Kidston et al. 2010). The BGC model could easily be extended to carbon and other nutrients but even in its simple form it can deliver useful BGC data products like primary productivity, which underpins the whole marine ecosystem (Brown et al. 2010).

The BGC model naturally sub-divides vertically into two domains: the photic zone where light is sufficient to allow phytoplankton photosynthesis and the deeper aphotic zone where there is no photosynthesis. For the example BGC model, the focus is on the photic zone where the evolution of our 4 BGC state variables nitrate (N), phytoplankton (P), zooplankton (Z) and detritus (D) are described by the following equations. For completeness the equations are presented here and for more detail refer to Kidston et al. (2010). The following model equations describe the evolution of BGC state variables in the mixed layer in mmol N/m3 where M is the mixed layer depth and Table 12.1 summarizes the BGC model parameters.

dt M

dD 2

dz M

— = ixdD + Y2Z + ixpP - J (M, t, N)P + --(NO - N) (12.4)

dt M

Parameter |
Symbol |
Value |
Units | |

Phytoplankton model parameters | ||||

Initial slope of P-I curve |
a |
0.256 |
day-1/(W m-2) | |

Photosynthetically active radiation |
PAR |
0.43 |
- | |

Light attenuation due to water |
k w |
0.04 |
m-1 | |

a |
0.27 |
day-1 | ||

Maximum growth rate parameters |
b |
1.066 |
- | |

c |
1.0 |
C-1 | ||

Half saturation constant for N uptake |
k |
0.7 |
mmol N/m-3 | |

Phytoplankton mortality |
VP |
0.01 |
day-1 | |

Zooplankton model parameters | ||||

Assimilation efficiency |
Yi |
0.925 |
- | |

Maximum grazing rate |
g |
1.575 |
day-1 | |

Prey capture rate |
s |
1.6 |
(mmol N/m-2)- |
1day-1 |

Quadratic mortality |
VZ |
0.34 |
(mmol N/m-3)- |
1day-1 |

Excretion |
Y2 |
0.01 bcT |
day-1 | |

Detritus model parameters | ||||

Remineralisation rate |
Vd |
0.048 bcT |
day-1 | |

Sinking velocity |
WD |
18.0 |
m/day |

The phytoplankton Eq. 12.1 includes phytoplankton growth (J(M, t, N)P), loss due to zooplankton grazing (G(P, Z)) and phytoplankton mortality (^pP), and changes due to ocean physics (last term in Eq. 12.1). Ocean physics includes the dilution of P as the mixed layer depth (M) deepens with time (h+(t) = min(0, dM)) while shallowing of M has no impact on P because no new water is added. The m in the ocean physics terms represents the vertical diffusive loss of phytoplankton from the mixed layer.

The key environmental drivers of the phytoplankton growth are temperature (T), light levels (I) and nitrate concentrations (N) and the growth rate is given by the following equations.

24 M

The light levels are influenced by the daily variation in the incident solar radiation at the surface of the ocean (I(0, t)), the fraction of the solar radiation available for photosynthesis (PAR), the depth of the mixed layer, and the extinc-

— kz tion of light with depth (e w ). In these equations Jmx is the maximum phytoplankton growth with unlimited light and N, which reduces under low nitrate concentrations.

The zooplankton Eq. 12.2 represents the balance between growth due to phytoplankton grazing (G(P, Z)) and losses due to zooplankton excretions (y2Z) and mortality (^Z2) and ocean physics (last term). The grazing of phytoplankton by zooplankton is given by the following equation s e P2

Only a fraction of the grazed phytoplankton (yt) goes directly into zooplankton growth while the remainder is transferred to detritus.

For detritus (Eq. 12.3), there is an input from phytoplankton mortality, and zooplankton grazing and mortality, while detritus decomposition (pPp) and sinking (wDdDD) provide losses. The last term in the P equation represents ocean physics which is identical to the Eq. 12.1. The detrital sinking term "exports" nitrogen from the model domain. The sinking of detritus into the ocean and subsequent decompo sition back into it inorganic form (nitrate) creates a vertical gradient in nitrate that increases with depth.

For nitrate (Eq. 12.4), phytoplankton uptake provides the loss, detritus reminer-alization provides a source with an additional external supply due to ocean physics term (I^+L+M [No - N]), by resupplying the nitrate back into the domain by specifying the nitrate concentration below the mixed layer (No). The physical processes supplying sub-surface water to the upper ocean are important to providing nitrate for phytoplankton growth and they are an important link between the physical and biological system.

The simple BGC model can be expanded to include multiple nutrients (e.g. carbon, alkalinity, iron, phosphate), phytoplankton (e.g. size and functional dependent), zooplankton and detritus state variables (e.g. Vichi et al. 2007). The additional BGC state variables add complexity to the model and introduce additional model parameters that must be specified. Hence, there is a compromise between model complexity and information available to constrain the BGC model (Matear 1995). Further the parameterizations in the BGC model often rely on empirical relationships which differs from the ocean physical model who's governing equations have a strong theoretical connection (e.g. Navier-Stokes equation). Therefore, the trade-off between using a simple versus complex biological model comes down to the questions one is trying to address.

The simple BGC model only represents processes occurring in the mixed layer where the fields are well mixed and the system can be configured as a 0D system. The only physical processes included in the BGC model are the vertical supply of nitrate due mixed layer deepening and the vertical mixing between the mixed layer and the deep water. If one wanted to include this BGC model into the GODAE system it would require explicitly resolving the vertical dimension (e.g. Schartau and Oschlies 2003a) and coupling the BGC state variables with the ocean circulation (e.g. Oschlies and Schartau 2005). In such a formation, the evolution of the BGC state variables would be influenced by the biological processes described above and the ocean dynamics which would transport the BGC state variables around the ocean.

Although the 3D BGC model would couple the physical and biochemical processes to evolve the BGC state variables it is important to note that within the photic zone generally the biological transformations are much greater than physical influences (e.g. phytoplankton doubling occurs over days where as ocean advection is more like weeks in global ocean models). Therefore within the photic zone over daily time-scales there is the potential to focus on biochemical processes. There are many examples where the BGC model is only applied in a vertical dimension just focusing on processes occurring in the mixed layer. The dominance of the biological processes over the physics fails in high flow regions where horizontal advection can transport the biological pools and this is clearly evident in the ocean colour imagery of Chlorophyll a where the spirally nature of the Chlorophyll a attest to the importance of the ocean dynamics (Fig. 12.1). Below the photic zone, in the aphotic zone the physics and biological processes are both important to the evolution of the BGC state variables and generally it is not possible to neglect either process (Fig. 12.2).

04 May 2000

Fig. 12.1 A 1-day SeaWiFS ocean colour image of surface Chlorophyll a (mg Chla/m3) from western Australia on May 4, 2000. The contour lines denote the sea surface height anomaly (SSHA) field from the same period. Note the two anti-cyclonic eddies (negative SSHA closed contours) connect with the shelf and transporting high Chlorophyll a water off-shore. (Moore et al. 2007)

Fig. 12.1 A 1-day SeaWiFS ocean colour image of surface Chlorophyll a (mg Chla/m3) from western Australia on May 4, 2000. The contour lines denote the sea surface height anomaly (SSHA) field from the same period. Note the two anti-cyclonic eddies (negative SSHA closed contours) connect with the shelf and transporting high Chlorophyll a water off-shore. (Moore et al. 2007)

13 12

13 12

150 200 Days since Jan 1

Phytoplankton Zooplankton Detritus

50 100 150 200

Days since Jan 1

Fig. 12.2 The seasonal evolution of the 4 state variables of the 0D BGC model for the SAZ-Sense P1 site for the model parameter values given in Table 12.1. Top: Nitrate concentrations, Bottom: Phytoplankton, zooplankton and detritus concentrations all in mmol N/m3

12.3 Biogeochemical Data Assimilation

The application of data assimilation separates into two types of problems: 1. parameter estimation and 2. state estimation. The two approaches reflected different philosophies on how to fuse the BGC model with the observations with both delivering useful but different information. In both approaches, data is used to constrain the evolution of the state variables but with parameter estimation it is the model parameters that are modified to fit the constraints while in the state estimation the state variables are modified to fit the observations. The state estimation approach generally forms the foundation of GODAE physical data assimilation with such effort delivering either ocean forecast or reanalysis products of the time evolving 3-D ocean physical state. However, BGC data assimilation studies have tended to focus on data assimilation for parameter estimation. Although a complete description of the two approaches is beyond the scope of this paper, the following is a brief summary of the application of these two data assimilation approaches.

### 12.3.1 Parameter Estimation

An obvious feature of the simple BGC model discussed in Sect. 12.2 is the large number of model parameters that must be specified to simulated the BGC. The setting of these parameters can be partially accomplished from observations but for many parameters no direct observation estimate is available. Further, even for the more observable parameters like the parameters that control phytoplankton growth (e.g. a and Jmax) much uncertainty still exists because either the parameter values change with time or the values determined from the individual phytoplankton species may not be applicable to the entire ecosystem. The inability to directly specify all the model parameters forces one to determine the parameters values by tuning the model to reproduce the observations. This is a tedious and time consuming approach. The attraction of data assimilation is that it provides a means to generate a set of model parameters that reflects the observations, determines the values of the poorly known parameters, and provides insight into which parameters are constrained by the model (Kidston et al. 2010; Schartau and Oschlies 2003b).

There is now a long list of studies which have used data assimilation methods to estimate model parameters of ecosystem models (see Gregg et al. 2009 for a list of these studies). There is even a webpage setup to explore parameter estimation of a suite of different BGC models at a number of different sites (http://www.ccpo.odu. edu/marjy/Testbed/Workshop1.html).

A recent model comparison study highlights some of the issues using BGC model (Friedrichs et al. 2007). In their study, they assessed the ability of 12 different BGC models of varying complexity to fit observations from two different sites.

They showed that all models performed equally well when calibrated at one site but only models with multiple plankton state variables were able to have the same model parameters for both sites. The study emphasizes the need for spatial variability in the model parameters to account for the different ecological regimes. Consequently, no one unique set of model parameters exists for the included ecosystem components such as phytoplankton for the entire ocean and flexibility in the model parameters needs to be incorporated in the model application.

The uncertainty in the biological model is not just in the model parameters but extends to the choice of equations used to describe the biological system (Franks 2009). Both parameter and model formulation uncertainty introduce large model error in BGC modelling, which lends itself to parameter estimation studies which explore the parameter space, model complexity and model formulation (e.g. Matear 1995). The ability to address all three of these issues demonstrates the value in the parameter estimation approach for BGC data assimilation. In addition, parameter estimation can provide insight into what observations are critical to building and constraining more realistic models, and identifying the critical model parameters required to reproduce the observations (Kidston et al. 2010). The latter result provides a convenient way to identify a subset of critical model parameters, which capture the key observed dynamics of the biological system, and then explore how these parameters affect the dynamics of the system (e.g. Friedrichs et al. 2006).

One important aspect to note with the parameter estimation approach is that unrealistic parameter values may be estimated because important processes are excluded from the model formulation (these are call structural errors in the model). Therefore, the estimated model parameter values must be ecologically assessed and deemed plausible otherwise the formulated model has structural errors.

### 12.3.2 Ocean State Estimation

The attraction of state estimation data assimilation is that it provides a way to incorporate both physical and biological observations into the numerical models to obtain the evolution of BGC fields that are dynamically consistent with the observations and provide a tool to extend the observations in both space and time (Lee et al. 2009). The application of state estimation is to overcome limitations in the model by correcting the ocean state to produce a more realistic evolution of the ocean state (Natvik and Evensen 2003a). The approach provides a way of limiting the impact of model errors (parameter, formulation, initialization and forcing) to better hindcast and forecast the ocean state.

The study of Gregg (2008) provides a nice review of sequential BGC data assimilation studies which I briefly summarize here. Not surprising the focus of these studies is on assimilating ocean colour surface Chlorophyll a concentrations into their BGC models since this data product provides the best spatial and temporal data of the ocean biological system.

The first example of sequential data assimilation directly inserted CZCS chlorophyll into a 3-dimensional model of the southeast US coast (Ishizaka 1990). They produced immediate improvements in their chlorophyll simulation but these improvements did not last more than a couple of days before the model simulation diverge from the observations. The divergence of the model simulation from the observations reflected a bias in the biological model to overestimate Chlorophyll a. Correcting such a bias would be crucial to obtaining better and longer lasting state estimation results. More recently, the Ensemble Kalman Filter (EnKF) has been used to assimilate SeaWiFS ocean color Chlorophyll a data into a 3D North Atlantic model (Natvik and Evensen 2003 a, b). They showed the EnKF updated ocean state was consistent with both the observed phytoplankton and nitrate concentrations. However, this study did not show any comparison to unassimilated data to enable a quantitative assessment of the impact BGC state estimation (Gregg et al. 2009). Subsequently, Gregg (2008) used the Conditional Relaxation Analysis Method (CRAM) to sequential assimilate multi-year SeaWiFS data into a 3D BGC model. Not surprisingly they improved surface Chlorophyll a estimates since this was the state variable being assimilated. A more independent assessment was made by looking at the impact of data assimilation on the simulated depth-integrated primary production which showed a much more modest improvement with data assimilation compared to the improvement in Chlorophyll a.

Recently, Hemmings et al. (2008) presented a nitrogen balance scheme with the aim of assimilating ocean color Chlorophyll a to improve estimates of the seawa-ter pCO2. They used 1D simulations at two sites in the North Atlantic (30N and 50N) to assess the performance of their scheme. The scheme exploits the covari-ance between Chlorophyll a and the other biological state variables at a fraction of the computational cost of the multivariate EnKF scheme. To do this Hemmings et al. (2008) use 1D model simulations with varying model parameters to extract the relationship between simulated Chlorophyll a and the other biological variables. This information is then used to project the assimilated Chlorophyll a data on to all biological state variables. They show the nitrogen balancing approach improves the ocean pCO2 simulation over the case where only phytoplankton and DIC are updated. At 30N, RMS error of the surface pCO2 was reduced by more than 50% (4 ^ atm to less than 2 ^ atm). While at 50N, the surface pCO2 RMS error showed little improvement but the bias in the field was reduced from -2 ^ atm to nearly zero.

At present, the application of state estimation has been confined to utilizing remotely sensed ocean colour Chlorophyll a. To extend this information from the surface into the ocean interior will required coupled biological-physical models. Further, many of the colour images are corrupted by clouds and filling these data gaps is another obvious outcome of state estimation. Finally, it quite possible that the fields that we are most interested in are not sufficiently observed to provide the spatial and temporal coverage desired (e.g. pCO2). For these cases, data assimilation provides tool to exploit the existing observations to generated the fields of interest.

12.4 Challenges of Adding BGC Data Assimilation to GODAE Systems

### 12.4.1 Background

The Global Ocean Data Assimilation Experiment (GODAE) was conceived to provide a 3D depiction of global ocean circulation at eddy resolution through data assimilation that was consistent with physical fields and dynamical constraints (Lee et al. 2009). The expansion of the GODAE state estimation to BGC fields is a natural evolution that exploits the GODAE effort to provide BGC state estimates consistent with both the physical and biological information(Brasseur et al. 2009).

GODAE effort involves three streams: (1) mesoscale ocean analysis and forecast, (2) initialization of seasonal-interannual prediction, and (3) state estimation (reanalysis) products (Lee et al. 2009). In discussing the issues of including BGC into the GODAE effort the focus is on how to modify the state estimation problem. A range of assimilation methods have been applied by various groups to produce the ocean state estimation, ranging from adjoint to sequential methods.

The adjoint method (e.g. MOVE from Japan and the ECCO group (www.ecco-group.org)) (Lee et al. 2009) is analogous to the parameter estimation approach discussed earlier but implemented such that the ocean state is obtained by optimizing not only the model parameters but also control variables such as the initial state of the ocean and the surface forcing. Although the adjoint-based estimation products are characterized by consistency with the physical model equations they required the formulation of an "adjoint model" to do the data assimilation. The development of the adjoint of the couple physical biological model is possible (e.g. Matear and Holloway 1995; Schlitzer 2002) but it is not a trivial task. For this reason, it is more realistic to expect the first GODAE systems to include BGC to be based on sequential data assimilation methods and the focus of the following discussion is on how to get BGC into this type of data assimilation system.

The sequential methods as implemented by various ocean reanalysis groups (e.g. Bluelink from Australia (Oke et al. 2008); FOAM from the United Kingdom (Martin et al. 2007); Mercator from France (Brasseur et al. 2005)) are typically computationally more efficient than the adjoint method. The sequential approach allows the estimated state to deviate from an exact solution of the underlying physical model by applying statistical corrections to the ocean state. Such corrections act as internal sources/sinks of heat, salt, and momentum. Interested readers should read Zaron (this issue) for more a thorough presentation of the sequential data assimilation. To apply sequential data assimilation one needs information on how observations of a state variable projects onto the other state variables at that time, that includes all state variables for all model grid points. This information is referred to as the multi-variate Background Error Covariances (BECs) and is calculated from the anomalies in the state variable evolution in an unassimilating model. In practice the BECs are estimated from either an ensemble of simulations with different initial values of the state variables (Brasseur et al. 2005) or computed from a multi-year run of just one simulation (Oke et al. 2008). The benefits to using BECs are: They reflect the length-scales and the anisotropy of the ocean circulation in different regions. They provide the information on the covariances between different state variables in a dynamically consistent way (we use the term dynamically consistent to describe an ocean state that can be generated by the model). Finally they are easily generalized to assimilate different observation types in a single step. Therefore, one might expect the approach could easily be generalized to include BGC information. There have been a number of attempts to do this (Natvik and Evensen 2003a, b).

We will now explore these issues using our simple BGC model a pedagogical tool to evaluate whether BECs provide a suitable avenue for BGC state estimation data assimilation.

### 12.4.2 Potential Issues and Solutions

A key data set for biogeochemical state estimation is ocean colour Chlorophyll a hence the following discussion will focus on how constraints on phytoplankton concentrations affect the BGC state variables. Using BECs appears an obvious way to implement BGC data assimilation, however the BGC model uncertainty and biases can have a huge impact on the estimated BECs and the resulting data assimilation state estimation. For 3D BGC modelling, the calculation of the BECs poses a challenge since the biochemical processes controlling the flows between state variables are uncertain. This uncertainty is reflected in model parameter uncertainty, model parameters evolving with time, model formulation uncertainty and errors in the model structure (i.e. model simplification and missing processes). Capturing this uncertainty with an ensemble approach (e.g. EnKF) will be computationally demanding for a 3D BGC model because it will be difficult to include a sufficient number of ensemble members to account for all the BGC model uncertainties. The EnOI method (Brassington et al. 2007) provides a more efficient calculation of the BECs for the physical state variables by using statistics generated from a 9-year run of the unassimilating ocean model. Although computationally doable with the 3D BGC model, the BECs calculated from this method are not time-varying, which is probably not appropriate for BGC state variables since their present state, e.g. whether phytoplankton are growing or not, alters the BECs (Hemmings et al. 2008).

BECs generation within the 3D ocean circulation model at present is not computationally feasible for GODAE systems, however if we are mainly interested in the BGC in the photic zone there is the potential to treat the 3D BGC model as a set of 0D mixed layer representation for each ocean surface grid point of the domain. The evolution of the BGC fields in the mixed layer is generally controlled by BGC processes and applying this domain decomposition is often a reasonable assumption. Examples of such an approach include most of the parameter estimation studies referred to in the previous section. From the oceans dynamics perspective the key physical processes driving the BGC fields are the vertical supply of nutrient to the photic zone and the availability of light in the photic zone. Both these processes were incorporated into the simple BGC model presented earlier.

The key physical information needed for the 0D BGC model is the vertical velocity, the vertical mixing rate, and the temporal evolution of the mixed layer depth and temperature. These are already standard state variables generated by GODAE data assimilation systems therefore they are available to be used by our 0D BGC model. For the BGC data assimilation additional assessments of these physical products are needed given their importance to the BGC model behaviour but this is best pursued through improvements in the physical data assimilation system. With the rapidly growing number of Argo temperature and salinity profiles the observational information to assess these fields and improve them physical data assimilation systems provides a clear way forward.

By focusing the calculation of the BECs on just the 0D BGC model of the photic zone, generating the BECs from an ensemble of simulations becomes computationally doable. But what should the BECs look like? Given the non-linearity of the BGC model we expect a complex relationship between the BGC state variables. The Hemmings et al. (2008) nitrogen balancing approach was based on developing mechanistic links between the different BGC state variables. Although they discuss the BECs complexity and time dependent nature it is instructive to review this issue with our simple 0-D BGC model.

### 12.4.3 Pilot BGC Data Assimilation Using BECs

The BGC model described in Sect. 12.3 was used to explore to explore the relationships between the BGC state variables. To compute the BECs, we perform 100 random perturbations of the initial conditions of these state variables on January 1. The perturbations to the initial state were gaussian with a standard deviations of 0.01 mmol N/m3 and correct to have a mean of zero to ensure nitrogen conservation. The simulations, which occur during the growth phase of the phytoplankton, reveals several important features of the BGC model and its BECs. The random perturbations of the initial BGC state variables do not cause a long term change in the BGC state variable trajectories and all perturbations return back to the original trajectory of the unperturbed model with an e-folding time-scale of about 25 days (Fig. 12.3). This reflects a bias for the BGC state variables to follow the same trajectory as the unassimilated model. The bias in the model behaviour could possibly be reduced if we first used data assimilation to estimate the model parameters from the observations prior to tackling the sequential data assimilation. However, we would never expect to account for all the model biases everywhere in the model domain and further we want the data assimilation to help correct these deficiency to produce a more realistic model behaviour. The decay in perturbations BGC state variables involves a damped oscillation which reflects the unbalanced nature of the initial perturbations. For a non-linear model like our 0D BGC model, it is common to generate an unbalanced state when the state variables are updated, which gener-

b Days since Jan 1

Fig. 12.3 For the 0D BGC model at the SAZ-Sense P1 site: a the standard deviation in the BGC state variable anomalies generated from an ensemble of 100 simulations obtained by randomly perturbing the initial values on January 1 by ±0.01 mmol N/m3. b Evolution the BGC state variable anomalies for just one of the ensemble member. The anomalies in the BGC state variables are define as the difference in the state variables from the solution given in Fig. 12.2

b Days since Jan 1

Fig. 12.3 For the 0D BGC model at the SAZ-Sense P1 site: a the standard deviation in the BGC state variable anomalies generated from an ensemble of 100 simulations obtained by randomly perturbing the initial values on January 1 by ±0.01 mmol N/m3. b Evolution the BGC state variable anomalies for just one of the ensemble member. The anomalies in the BGC state variables are define as the difference in the state variables from the solution given in Fig. 12.2

ates an undesirable response in the model and this condition should be avoided. The simulated relationship between phytoplankton (P) and other state variables is also complex (Fig. 12.4) with clear negative relationship with nitrate (N) but no obvious one with zooplankton (Z) and detritus (D) (not shown). The negative phytoplankton-nitrate relationship is the phytoplankton growth dependent response where more phytoplankton equates to greater phytoplankton growth and increased nitrate uptake and reduced nitrate concentrations and vice versa. However, the phy-toplankton-zooplankton relationship does not appear to reveal any pattern. The link between P and Z becomes clearer when we plot the phase diagram of the 100 per-

0.012- |
! È |
0.0080.004 -0.000 |
" "^"if „ rf-„fr. |
- | ||||

!| o |
i * »< |
- | ||||||

CL |
- 0.008- 0.012 |
x* x |
- | |||||

- . |
1 1 1 1 1 1 1 1 008 - .004 0.000 0.004 |
0.008 | ||||||

a |
1 1 1 1 1 1 1 1 1 |
1 | ||||||

0.012 |
- | |||||||

to £ |
0.008 |
X X * x X * |
- | |||||

Se |
0.004 |
xX xx " x * |
- 0.008 - 0.012 |
* x * * * * * *x * x Xx x * xX" Xx«" X « |
x |
— |

- 0.02 - 0.01 0.00 0.01 0.02 b Zooplankton anomaly (mmol N / m3)

- 0.02 - 0.01 0.00 0.01 0.02 b Zooplankton anomaly (mmol N / m3)

Fig. 12.4 On January 7 from the 100 member ensemble simulations, the relationship between the anomaly in phytoplankton concentrations with the anomaly in nitrate (a), and zooplankton (b) concentrations. The anomalies in the BGC state variables are define as the difference in the state variables from the solution given in Fig. 12.2

turbations (Fig. 12.5). Because of the lag between Z response to P perturbation, we miss the strong connection between the two state variables. The lag between P and Z creates the unbalanced response of the BGC model to random perturbations, which causes limit cycles in our model simulations. Capturing this time-evolving connection between Z to P is not trivial and further it would be dependent on the choice of model parameter values and the BGC state prior to the perturbation.

For all the issues discussed above the use of the BECs for the BGC data assimilation is complicated and not an attractive way to tackle BGC state estimation. How should we do it?

By focusing on one simulation where we increase the P concentration by 0.01 mmol N/m3 and observe how the system evolves can help us develop an ap-

- 0.02 - 0.01 0.00 0.01 0.02 Zooplankton anomaly (mmol / m3)

Fig. 12.5 The phase relationship between phytoplankton and zooplankton anomaly concentrations from the ensemble of simulations for 50 days after initialization. Each asterisk is a day and line is an ensemble member. The anomalies in the BGC state variables are define as the difference in the state variables from the solution given in Fig. 12.2

- 0.02 - 0.01 0.00 0.01 0.02 Zooplankton anomaly (mmol / m3)

Fig. 12.5 The phase relationship between phytoplankton and zooplankton anomaly concentrations from the ensemble of simulations for 50 days after initialization. Each asterisk is a day and line is an ensemble member. The anomalies in the BGC state variables are define as the difference in the state variables from the solution given in Fig. 12.2

proach to exploiting the P observations more effectively. This simulation is analogous to the situation where the BGC model simulation under-estimates P and we want the model to behave more like the observed value. Under such a state perturbation within a few days the perturbed phytoplankton concentration rapidly returns to the pre-nudge value and little is gained from the assimilation with the exception we have generated an unbalance perturbation which generates an oscillation in the BGC state variables that persists for many more days (Fig. 12.6a). The unbalanced nature of the state update is most apparent in the fields not directly assimilated, for example Z shows much greater variability than P. Without continual data assimilation we would not maintain a realistic P and one would continue to produce spurious behaviour in the other BGC state variables. However, for this site where nitrate is always in excess (Fig. 12.2) the more direct way to control P is through its loss to grazing. By reducing grazing by 10% (Fig. 12.6b) a much more sustained phytoplankton response is achieved and the response of the other BGC state variables lack the spurious variability because the change in the BGC state is balanced. Such a modification of the BGC model can be justified on the grounds that the grazing parameter is uncertain.

The modification to grazing is correcting the model's under-estimate of P by altering the underlying equations rather than just changing P to reflect the observations. Changing the grazing rate alters the trajectory of the model simulation leading to a persistent increase in P, Z, D and a persistent decrease in N (Fig. 12.6b). The persistent nature of the changes in the BGC state variables is a desirable feature since it demonstrates that "assimilating" P can have a lasting impact on the model simulation and thus avoid the situation where data assimilation only provides a short term improvement to the BGC state evolution (e.g. Ishizaka 1990).

0.01

0.008

0.006

0.004

0.002

0.002

0.004

0.06 | |

0.05 | |

1 | |

z |
0.04 |

Ö | |

E | |

£ | |

0.03 | |

n | |

o | |

ra | |

■u |
0.02 |

n | |

e | |

c | |

n | |

o o |
0.01 |

>. | |

ra | |

E o |
0 |

n | |

< | |

- 0.01 | |

- 0.02 |

10 15 20 25 Days since Jan 1

10 15 20 25

Days since Jan 1

Fig. 12.6 a 0D BGC model evolution of the state variable anomalies for 0.01 mmol N/m3 additional to the phytoplankton concentration. b State variable anomalies when the zooplankton grazing rate on phytoplankton is decreased by 10%. In the two figures, the BGC state variable anomalies are denoted by nitrate (dark blue), phytoplankton (cyan), zooplankton (yellow), and detritus (red) all in units of mmol N/m3. The anomalies in the BGC state variables are define as the difference in the state variables from the solution given in Fig. 12.2

Now I accept this is a contrived example where the model has an obvious bias but this is probably always expected given model parameter and formulation uncertainties. In previous twin experiments where the sequential data assimilation was applied to the BGC state variables and judged to be a success it must be emphasized that model generated observations are used in the data assimilation therefore they are unbiased with respect to the model used to do the data assimilation (Eknes and Evensen 2002). Therefore the impact of model bias is omitted in their assessment of the success of the data assimilation, which would not be the case in a real world application.

The simple example highlights the benefits of perturbing the model parameters and Dowd (2007); Mattern et al. (2010); Jones et al. (2010) all provide more extensive examples of how to perturb the underlying model parameters to reproduce the observed state evolution of the BGC fields. These balanced methods (i.e. only the model parameters are changed) deliver estimates of parameters concurrently with the model state without the assumption that model parameters are constant in time. Using the time evolving parameter values to alter the BGC state variables to fit the observations has ecological meaning since we expect biological processes controlling these parameter values to be spatially and temporally varying. Ideally, for the GODAE BGC data assimilation we would like to carry the ensemble of different model parameters like in the 1D studies by Dowd (2007); Mattern et al. (2010); Jones et al. (2010) into the 3D BGC model but the computational overhead of running multiple BGC models is not possible. Hence, one will need to reduce the potential parameters to just one parameter set. Perhaps we can make some ecological defensible choices in how the model parameters evolve in time and space in a way that is consistent with model parameter uncertainty and uncertainty in the observations.

12.4.4 Proposed Sequential Data Assimilation for BGC State Estimation

In the previous section, I discussed several issues in applying state estimation for BGC fields. I would now like to propose an approach to incorporating the BGC into the state estimation of a GODAE system. The following is the proposed sequence of steps to do sequential BGC state estimation, which is also outlined in Fig. 12.7.

1. convert the observed remotely sensed ocean colour Chlorophyll a concentration into an estimate of surface phytoplankton concentrations to provide the observational constraint for the BGC data assimilation. Note, that an alternative would be to included Chlorophyll a as a state variable in the BGC model and included in the BGC model a formulation for the N:Chlorophyll a ratio of phytoplankton (e.g. Hemmings et al. (2008)).

2. using the seasonal climatology of derived phytoplankton observations apply parameter estimation to the 0D BGC representation of the 3D BGC model at

PD1 (observation) Pi (4)

## Guide to Alternative Fuels

Your Alternative Fuel Solution for Saving Money, Reducing Oil Dependency, and Helping the Planet. Ethanol is an alternative to gasoline. The use of ethanol has been demonstrated to reduce greenhouse emissions slightly as compared to gasoline. Through this ebook, you are going to learn what you will need to know why choosing an alternative fuel may benefit you and your future.

## Post a comment