tests/testthat/test-phytoplankton-fluxes.R

In StrathE2E2: End-to-End Marine Food Web Model

test_that("check phytoplankton nutrient uptake flux calculations", {

#-----------------------------------------------------------------------------------------------------------------
#
# This test uses a 'Testbed' version of the model in which all of the external drivers are set to constant values.
# We expect the model outputs to converge to a steady state under these conditions. Then we take the masses of the
# food web variables from the steady state model output, and manually derive the uptake fluxes between them.
# These manually derived fluxes are then compared with the fluxes generated in the 'flux_matrix' outputs from the model.
# So the test is checking that the flux calculations in the C-code of the package are correctly evaluated.
#
# The Testbed model is stored in /StrathE2E2/tests/testdata/models/
#
# The model has been run for >1000 years to ensure a steady state, and the end-state saved as initial conditions.
# So the model should be at steady state from time 0. In fact, the input csv file containing the initial conditions
# does not hold sufficient precision (number of decimal places) to ensure perfect steady state from the onset of a run.
# This is manifest as some tiny oscillations in the first year or so of a run before everything settles to its steady state.
# Hence we run the Testbed model for 10 years and then base the test on the final year of the run.
# 
# For convenience, the Testbed model uses the fitted parameters from the North Sea implementation of StrathE2E2.
# 
# An important points to note:
# - The temperature driving data for all three spatial compartmets of the Testbed
# model are set to 10C, wich matches the Q10 reference temperature in the fixed parameter inputs. This means that when manually
# calculating the uptake fluxes there is no need to compute any temperature corrections for the maximum uptake rate
# parameters. However, the C-code in the model IS doing the Q10 calculation, so this is an important aspect of the testing.
# As a consequence, it is very important not to chnage the temperature driving values in the Testbed model - they must be
# equal to the Q10 reference temperature value.
#
# This test applies the above approach to uptake by nitrate and ammonia by phytoplankton.
#
#-----------------------------------------------------------------------------------------------------------------
#
# The following non-exported functions in the package are needed for the test:
# source("internal.R")
# source("read_fitted_parameters.R")
# source("read_fixed_parameters.R")
# source("read_physics_drivers.R")
# source("read_physical_parameters.R")
#
#-----------------------------------------------------------------------------------------------------------------

#Some additional functions required for assembling the tests

#######################
#Function to pull the final year annual average mass of a variable out of the results
pull_final_year_mean<-function(results,vtoget){
vgot<-(elt(results$output, vtoget))
vmean<-mean(vgot[(length(vgot)-360):(length(vgot))])
vmean
}
######################
#Function to pull flux data out of the flow matrix generated by the model
extract_model_flux<-function(matrixname,preyname,predname){
predi<-which(colnames(matrixname)==predname)
preyi<-which(colnames(matrixname)==preyname)
model_flux<-(matrixname[preyi,predi]/360)
model_flux
}
#######################
#Autotrophic uptake function
f2<-function(a,b,c,k1,k2,k3)
#          nut,phyt,light,umax,hs,Lsat
{
  x<-c/k3
  minlit<-min(1,x)
  return ((minlit * (b*k1*a))/(a+k2))
}
#######################




#-----------------------------------------------------------------------------------------------------------------

#Run the Testbed model and pull out data from the results object

#Run the testbed model - assuming R home dir is "/GitLab/StrathE2E2/tests/testthat"
model<-e2e_read(model.name="Testbed",
                  model.variant="Const",
                  model.ident="base",
                  models.path="../testdata/models")
results <- e2e_run(model,nyears=20,csv.output=FALSE)
#e2e_plot_ts(model,results,selection="ECO")
#e2e_extract_start(model, results)

#Pull out the inshore and offshore final year data for phytoplankton and zoopkankton groups

#names(results$output[1:78]) # Lists the names of state varibales in the header

phyt_so<-pull_final_year_mean(results,"phyt_so")
phyt_si<-pull_final_year_mean(results,"phyt_si")

nitrate_so<-pull_final_year_mean(results,"nitrate_so")
nitrate_si<-pull_final_year_mean(results,"nitrate_si")

ammonia_so<-pull_final_year_mean(results,"ammonia_so")
ammonia_si<-pull_final_year_mean(results,"ammonia_si")


#Pull out the flow matrix from the results
flow_matrix<-elt(results$final.year.outputs,"flow_matrix_all_fluxes")
flow_matrix_ns<-elt(results$final.year.outputs,"flow_matrix_excl_spawn_recruit")

#colnames(flow_matrix)   #Lists the column names in the flow matrix
#rownames(flow_matrix)   #Lists the row names in the flow matrix

#-----------------------------------------------------------------------------------------------------------------




#Pull data out of the model object, ie data that have been assembled from the input csv files

#Pull out the model configuration to get areas and volumes

model.path<-model$setup$model.path
	read.model.setup(model.path)			# Models/Model/Variant/MODEL_SETUP.csv
	physical.parameters	<- read_physical_parameters(model.path)
	so_depth<-elt(physical.parameters,"so_depth")
	si_depth<-elt(physical.parameters,"si_depth")
	d_depth<-elt(physical.parameters,"d_depth")
        area_s0<-elt(physical.parameters,"x_area_s0")
        area_s1<-elt(physical.parameters,"x_area_s1")
        area_s2<-elt(physical.parameters,"x_area_s2")
        area_s3<-elt(physical.parameters,"x_area_s3")
        area_d0<-elt(physical.parameters,"x_area_d0")
        area_d1<-elt(physical.parameters,"x_area_d1")
        area_d2<-elt(physical.parameters,"x_area_d2")
        area_d3<-elt(physical.parameters,"x_area_d3")
shallowprop<-(area_s0+area_s1+area_s2+area_s3)
volume_so<-so_depth*(1-shallowprop)
volume_d<-d_depth*(1-shallowprop)
volume_si<-si_depth*(shallowprop)


# Pull out the uptake function parameters needed to manually calculate the fluxes

model.path<-model$setup$model.path
	read.model.setup(model.path)			# Models/Model/Variant/MODEL_SETUP.csv
	fitted.parameters	<- read_fitted_parameters(model.path)
	fixed.parameters	<- read_fixed_parameters(model.path)
        physics.drivers         <- read_physics_drivers(model.path)

satlight_phyt <- elt(fixed.parameters,"satlight_phyt")

sslight<-physics.drivers$sslight[1]
so_logespm<-physics.drivers$so_logespm[1]
si_logespm<-physics.drivers$si_logespm[1]

u_phyt<-elt(fitted.parameters,"u_phyt")
h_phyt<-elt(fitted.parameters,"h_phyt")

PREF_NIT_phyt<-elt(fitted.parameters,"PREF_NIT_phyt")
PREF_AMM_phyt<-elt(fitted.parameters,"PREF_AMM_phyt")

inshore_phyt_depth_prop <- elt(physical.parameters,"xinshore_phyt_prop_depth")
lightSPM_intercept      <- elt(physical.parameters,"xlightSPM_intercept")
lightSPM_slope          <- elt(physical.parameters,"xlightSPM_slope")

#Light calculations

 kvertattn_o <-  lightSPM_intercept + lightSPM_slope * exp(so_logespm)     # Parameterised relationship between light vertical attenuation coefficient and SPM
 kvertattn_i <-  lightSPM_intercept + lightSPM_slope * exp(si_logespm)     # Parameterised relationship between light vertical attenuation coefficient and SPM

 phyt_propsl_o <-        (((1/kvertattn_o)*exp(-kvertattn_o*0)) - ((1/kvertattn_o)*exp(-kvertattn_o*so_depth)))/(so_depth)
 phyt_propsl_i <-        (((1/kvertattn_i)*exp(-kvertattn_i*0)) - ((1/kvertattn_i)*exp(-kvertattn_i*(si_depth*inshore_phyt_depth_prop))))/(si_depth*inshore_phyt_depth_prop)

 phyt_S_layer_light_o <- phyt_propsl_o*sslight
 phyt_S_layer_light_i <- phyt_propsl_i*sslight


#-----------------------------------------------------------------------------------------------------------------

#Now calculate the fluxes

calc_flux_nit_phyt_o <- f2(nitrate_so,phyt_so,phyt_S_layer_light_o,u_phyt*PREF_NIT_phyt,(h_phyt*volume_so),satlight_phyt)
calc_flux_nit_phyt_i <- f2(nitrate_si,phyt_si,phyt_S_layer_light_i,u_phyt*PREF_NIT_phyt,(h_phyt*volume_si),satlight_phyt)

calc_flux_nit_phyt_w <- calc_flux_nit_phyt_o + calc_flux_nit_phyt_i

model_flux_nit_phyt<-extract_model_flux(flow_matrix,"wcnitrate","phyt")

#calc_flux_nit_phyt_w
#model_flux_nit_phyt

#............................

calc_flux_amm_phyt_o <- f2(ammonia_so,phyt_so,phyt_S_layer_light_o,u_phyt*PREF_AMM_phyt,(h_phyt*volume_so),satlight_phyt)
calc_flux_amm_phyt_i <- f2(ammonia_si,phyt_si,phyt_S_layer_light_i,u_phyt*PREF_AMM_phyt,(h_phyt*volume_si),satlight_phyt)

calc_flux_amm_phyt_w <- calc_flux_amm_phyt_o + calc_flux_amm_phyt_i

model_flux_amm_phyt<-extract_model_flux(flow_matrix,"wcammonia","phyt")

#calc_flux_amm_phyt_w
#model_flux_amm_phyt

#............................
#............................

#Implement testthat checks...

expect_equal(model_flux_nit_phyt,calc_flux_nit_phyt_w , 1e-7)
expect_equal(model_flux_amm_phyt,calc_flux_amm_phyt_w , 1e-7)

})