metab.kalman {LakeMetabolizer}

Metabolism calculated from parameters estimated using a Kalman filter

Description

A state space model accounting for process and observation error, with the maximum likelihood of parameteres estimated using a Kalman filter. Also provides a smoothed time series of oxygen concentration.

Usage

metab.kalman(do.obs, do.sat, k.gas, z.mix, irr, wtr, ...)

Arguments

do.obs

Vector of dissovled oxygen concentration observations, mg O[2] L^{-1}

do.sat

Vector of dissolved oxygen saturation values based on water temperature. Calculate using o2.at.sat

k.gas

Vector of kGAS values calculated from any of the gas flux models (e.g., k.cole) and converted to kGAS using k600.2.kGAS

z.mix

Vector of mixed-layer depths in meters. To calculate, see ts.meta.depths

irr

Vector of photosynthetically active radiation in \mu mol\ m^{-2} s^{-1}

wtr

Vector of water temperatures in ^{\circ}C. Used in scaling respiration with temperature

...

additional arguments; currently "datetime" is the only recognized argument passed through ...

Details

The model has four parameters, c_1, c_2, Q, H, and consists of equations involving the prediction of upcoming state conditional on information of the previous state (a_{t|t-1}, P_{t|t-1}), as well as updates of those predictions that are conditional upon information of the current state (a_{t|t}, P_{t|t}). a is the

v=k.gas/z.mix

a_t = c_1*irr_{t-1} + c_2*log_e(wtr_{t-1}) + v_{t-1}*do.sat_{t-1}

beta = e^{-v}

do.obs_t = a_t/v_{t-1} + -e^{-v_{t-1}}*a_t/v_{t-1} + beta_{t-1}*do.obs_{t-1} + epsilon_t

The above model is used during model fitting, but if gas flux is not integrated between time steps, those equations simplify to the following:

F_{t-1} = k.gas_{t-1}*(do.sat_{t-1} - do.obs_{t-1})/z.mix_{t-1}

do.obs_t=do.obs_{t-1}+c_1*irr_{t-1}+c_2*log_e(wtr_{t-1}) + F_{t-1} + epsilon_t

The parameters are fit using maximum likelihood, and the optimization (minimization of the negative log likelihood function) is performed by optim using default settings.

GPP is then calculated as mean(c1*irr, na.rm=TRUE)*freq, where freq is the number of observations per day, as estimated from the typical size between time steps. Thus, generally freq==length(do.obs).

Similarly, R is calculated as mean(c2*log(wtr), na.rm=TRUE)*freq.

NEP is the sum of GPP and R.

Value

A data.frame with columns corresponding to components of metabolism

GPP

numeric estimate of Gross Primary Production, mg O_2 L^{-1} d^{-1}

R

numeric estimate of Respiration, mg O_2 L^{-1} d^{-1}

NEP

numeric estimate of Net Ecosystem production, mg O_2 L^{-1} d^{-1}

Use attributes to access more model output:

smoothDO

smoothed time series of oxygen concentration (mg O[2] L^{-1}), from Kalman smoother

params

parameters estimated by the Kalman filter (c_1, c_2, Q, H)

Note

If observation error is substantial, consider applying a Kalman filter to the water temperature time series by supplying wtr as the output from temp.kalman

Author(s)

Ryan Batt, Luke A. Winslow

References

Batt, Ryan D. and Stephen R. Carpenter. 2012. Free-water lake metabolism: addressing noisy time series with a Kalman filter. Limnology and Oceanography: Methods 10: 20-30. doi: 10.4319/lom.2012.10.20

See Also

temp.kalman, watts.in, metab, metab.bookkeep, metab.ols, metab.mle, metab.bayesian

Examples

library(rLakeAnalyzer)
doobs <- load.ts(system.file('extdata',
                            'sparkling.doobs', package="LakeMetabolizer"))
wtr <- load.ts(system.file('extdata',
                          'sparkling.wtr', package="LakeMetabolizer"))
wnd <- load.ts(system.file('extdata',
                          'sparkling.wnd', package="LakeMetabolizer"))
irr <- load.ts(system.file('extdata',
                          'sparkling.par', package="LakeMetabolizer"))

#Subset a day
Sys.setenv(TZ='GMT')
mod.date <- as.POSIXct('2009-07-08', 'GMT')
doobs <- doobs[trunc(doobs$datetime, 'day') == mod.date, ]
wtr <- wtr[trunc(wtr$datetime, 'day') == mod.date, ]
wnd <- wnd[trunc(wnd$datetime, 'day') == mod.date, ]
irr <- irr[trunc(irr$datetime, 'day') == mod.date, ]

k600 <- k.cole.base(wnd[,2])
k.gas <- k600.2.kGAS.base(k600, wtr[,3], 'O2')
do.sat <- o2.at.sat.base(wtr[,3], altitude=300)

metab.kalman(irr=irr[,2], z.mix=rep(1, length(k.gas)),
            do.sat=do.sat, wtr=wtr[,2],
            k.gas=k.gas, do.obs=doobs[,2])