metab.mle {LakeMetabolizer}

Metabolism calculated from the maximum likelihood estimates of the parameters in a standard linear regression model

Description

Process-error-only model with parameters fitted via maximum likelihood estimation (MLE). This function runs the maximum likelihood metabolism model on the supplied gas concentration and other supporting data.

Usage

metab.mle(do.obs, do.sat, k.gas, z.mix, irr, wtr, error.type = "OE", ...)

Arguments

do.obs

Vector of dissolved 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

error.type

Option specifying if model should assume pure Process Error 'PE' or Observation Error 'OE'. Defaults to observation error 'OE'.

...

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

Details

The model has the three parameters, c_1, c_2, epsilon, and has the form

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}

The maximum likelihood estimates of model parameters can be accessed via attributes(metab.mle(...))[["params"]]

Note

Currently, missing values in any arguments will result in an error, so freq must always equal nobs.

Author(s)

Luke A Winslow, Ryan Batt, GLEON Fellows

References

Hanson, PC, SR Carpenter, N Kimura, C Wu, SP Cornelius, TK Kratz. 2008 Evaluation of metabolism models for free-water dissolved oxygen in lakes. Limnology and Oceanography: Methods 6: 454:465.

Solomon CT, DA Bruesewitz, DC Richardson, KC Rose, MC Van de Bogert, PC Hanson, TK Kratz, B Larget, R Adrian, B Leroux Babin, CY Chiu, DP Hamilton, EE Gaiser, S Hendricks, V Istvanovics, A Laas, DM O'Donnell, ML Pace, E Ryder, PA Staehr, T Torgersen, MJ Vanni, KC Weathers, G Zhuw. 2013. Ecosystem Respiration: Drivers of Daily Variability and Background Respiration in Lakes around the Globe. Limnology and Oceanography 58 (3): 849:866. doi:10.4319/lo.2013.58.3.0849.

See Also

metab, metab.bookkeep, metab.ols, metab.kalman, 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
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, ]
z.mix = ts.thermo.depth(wtr)

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.mle(doobs[,2], do.sat, k.gas, z.mix[,2], irr[,2], wtr[,3])