This document is compiled by “knitting” the Rmarkdown file knit_calib_pmodel.Rmd. All source code is available in the repository calib_pmodel on bitbucket, group labprentice (login required).

1. Collect data

Read forcing data and P-model output into a common dataframe. To collect all daily forcing data and flux measurement data from the FLUXNET 2015 dataset, fAPAR forcing data, and SOFUN output data into a single data frame, use get_modobs.R. First adjust its header definitions:

##--------------------------------------------------------------------
## MANUAL SETTINGS
##--------------------------------------------------------------------
myhome               = "~/"
simsuite             = "fluxnet2015"  
outputset            = c( "s14" )
add_swcvars          = TRUE
add_swcvars_etbucket = FALSE
overwrite            = TRUE
overwrite_dosites    = TRUE 
outdir               = "~/data/fluxnet_sofun/"
##--------------------------------------------------------------------

Then run in R:

source("get_modobs.R")

This performs a data “cleaning” (see clean_fluxnet.R), calculates additional variables (e.g. soilm_mean, soilm_obs_mean) and writes all data as an R data frame (tibble) that we’re using for the calibration into the file df_modobs_fluxnet2015_*.Rdata.

1. P-model calibration

After writing df_modobs_fluxnet2015_*.Rdata by the previous step, we simply have to read in the data, remove cold and dry days and calibrate the apparent quantum yield efficiency parameter kphio_app using GPP from the FLUXNET 2015 data (GPP_NT_VUT_REF) and the linear regression function lm() using least-squares minimisation. After running SOFUN with an arbitrary parameter value kphio_app, the calibrate_pmodel.R calculates the optimised kphio_app to best match observations during relatively moist and warm days. SOFUN writes parameter values as meta information into NetCDF outputs. calibrate_pmodel.R reads this and and calculates by how much this (arbitrary) value has to be scaled to best match observations. These are the steps:

  1. Load the calibration function
source("calib_pmodel.R")
  1. Load the data collected by ‘get_modobs.R’
load( "df_modobs_fluxnet2015_s14_with_SWC_v4.Rdata" )  # loads df_fluxnet
ndays_0 <- nrow( df_fluxnet )
print( paste( "total number of days, unfiltered: ", ndays_0 ) )
## [1] "total number of days, unfiltered:  449620"
  1. Aggregate to monthly data
require(dplyr)
require(lubridate)
df_fluxnet <- df_fluxnet %>% group_by( mysitename, year(date), month(date) ) %>%
                                                         summarise_all( mean, na.rm=TRUE )
  1. Read GPP NetCDF output to get parameter value of quantum yield efficiency that was used in simulations
## simulation set
iset <- "s14"

## a random site
runname <- "FR-Pue"

## output directory
dirn <- paste( myhome, "sofun/output_nc_fluxnet2015_sofun/", iset, "/", sep="" )

## a random variable
filn <- paste0( runname, ".d.gpp.nc" )

## read meta information
require(ncdf4)
path         <- paste0( dirn, filn )
nc           <- nc_open( path )
kphio_used   <- ncatt_get( nc, varid=0, attname="param_kphio_GrC3" )$value %>% as.numeric()
fapar_source <- ncatt_get( nc, varid=0, attname="fapar_source" )$value
nc_close(nc)
  1. Evaluate optimised parameter for a range of soil moisture cutoffs
soilm_cutoff_range <- seq(0.1, 0.9, 0.05)
list_sens_soilm <- lapply( soilm_cutoff_range, function(x) calib_pmodel( df_fluxnet, 5.0, x, kphio_used ) )
kphio_opt  <- lapply( list_sens_soilm, function(x) x$kphio_opt ) %>% unlist()
kphio_corr <- lapply( list_sens_soilm, function(x) x$kphio_corr ) %>% unlist()
df_sens_soilm_cutoff <- tibble( soilm_cutoff=soilm_cutoff_range, kphio_opt=kphio_opt, kphio_corr=kphio_corr )
with( df_sens_soilm_cutoff, plot( soilm_cutoff, kphio_opt, pch=16, main=paste0("Sensitivity to soil moisture cutoff (temp. cutoff = 5 deg C, fAPAR data: ", fapar_source, ")"), font.main=1, cex.main=0.9 ))

  1. Evaluate optimised parameter for a range of temperature cutoffs
temp_cutoff_range <- seq(0, 20, 1)
list_sens_temp <- lapply( temp_cutoff_range, function(x) calib_pmodel( df_fluxnet, x, 0.4, kphio_used ) )
kphio_opt  <- lapply( list_sens_temp, function(x) x$kphio_opt ) %>% unlist()
kphio_corr <- lapply( list_sens_temp, function(x) x$kphio_corr ) %>% unlist()
df_sens_temp_cutoff <- tibble( temp_cutoff=temp_cutoff_range, kphio_opt=kphio_opt, kphio_corr=kphio_corr )
with( df_sens_temp_cutoff, plot( temp_cutoff, kphio_opt, pch=16, main=paste0("Sensitivity to temperature cutoff (soilm. cutoff = 0.4, fAPAR data: ", fapar_source, ")"), font.main=1, cex.main=0.9 ))

  1. Plot modelled vs. observed for given soil moisture (0.4) and temperature (5 deg C) cutoff ==> correction factor is 0.7785415. Optimal kphio is 0.05792348.
source("analyse_modobs.R")
stats <- with( filter( df_fluxnet, temp>5 & soilm_mean>0.4), analyse_modobs( GPP_NT_VUT_REF, 0.7785415 * gpp, yintersect0=FALSE, xlab=expression( paste( "modelled GPP (gC m"^{-1}, " d"^{-1}, ")")), ylab=expression( paste( "observed GPP (gC m"^{-1}, " d"^{-1}, ")")), do.plot=TRUE ) )