Two leaf canopy representation and solar info

[davidorme]

As we discussed on 16th May (@MarionBWeinzierl and @AmyOctoCat weren't there), Keith Bloomfield in Colin's group (@kjbloom) is interested in adding a new canopy representation to pyrealm. This is different from the canopy model that @AmyOctoCat is working on in #231 as part of the community development in #227 and elsewhere.

The key difference here affects what currently happens in PModel.estimate_productivity. That is currently implements what is usually called a "bigleaf" model: all of the incident sunlight is basically expected to be orthogonal to a single continuous leaf surface. The incoming radiation is µmol m2 s1, so the predictions are all rates based on a single big square leaf with overhead light.

Keith wants to implement what is known as a "two leaf, two stream" model. This represent the canopy as two partitions of sunlight and shade leaves and then also partitions direct and diffuse light. There is the added complexity that you can consider the beam azimuth of the light, which affects how the light interacts with the canopy. We could apply that to get proper diurnal conditions at subdaily timescales but also more accurately estimate GPP at a daily scale by using the correct azimuth angle for solar noon given the latitude.

I think we can replace the pmodel.estimate_productivity method with alternative canopy classes (e.g. BigleafProductivity, TwoLeafProductivity) and those should be useable with both the standard and subdaily PModel implementations.

The two leaf model would then also need users to specify beam angle as well as FAPAR and PPFD, so we’d need to provide methods to estimate those angles from latitude and time of day. There is the solar functionality in the upcoming v1.0.0 of pyrealm - it adds SPLASH - but that is currently focussed on daily solar fluxes. We might be able to re-factor some of it to produce the subdaily angles, but I’ve no problem using existing packages, as long as they support calculations for arrays of multiple sites.

Keith has provided a draft implementation using R:

#################################################################################
### kjb workings 2024, incorporating a SIMPLE canopy radiative transfer model ###
#################################################################################

## After several rounds of discussion, we try to follow the workings of the BESS model (Ryu et al. 2011) which draws heavily on de Pury & Farquhar 1997):
# TWO LEAF: sun and shade groups
# TWO STREAM: direct beam and diffuse light

### June 2023, Victor Flo Sierra has written an R function for irradiance partitioning based on the worked example provided in deP&F.

### To aid comparisons we draw on the same FLUXNET2015 subset as presented in Mengoli et al. (2021).  
# Our first step is to generate a sub-daily tibble for input to the P-model function.

library(tidyverse)
library(lubridate)


### We have the FLUXNET file from Giulia in csv format

df_Vie_hh <- read_csv(file = "Mengoli et al 2018_input-data/2014/FLX_BE-Vie_2014.csv",
                      show_col_types = FALSE) %>% 
  
  # render the time-stamp field more tractable
  mutate(step = ymd_hm(TIMESTAMP_START),
         date = as_date(step),
         month = lubridate::month(date),
         sitename = "BE-Vie") %>% 
  
  # set our desired date-range, here August
  filter(month == 8) %>% 
  
  # apply quality control filter to exclude medium to poor records:
  filter(NEE_CUT_REF_QC <= 1, VPD_F_QC <= 1, TA_F_QC <= 1, CO2_F_MDS_QC <= 1, PA_F_QC <= 1) %>% 
  
  # sort out our units for PA (here kPa) and VPD (here hPa); I'm going to leave PPFD as µmol Photons:
  mutate(vpd = VPD_F * 100,
         patm = PA_F * 1000) %>% 
  
  # select the variables we need for the P-model function
  select(sitename, date, step, PPFD_IN, vpd, TA_F, CO2_F_MDS, patm, GPP_DT_CUT_REF) %>% 
  
  # adopt the input names recognised by rpmodel:
  rename(tc = TA_F, ppfd = PPFD_IN, co2 = CO2_F_MDS, gpp_obs = GPP_DT_CUT_REF)


## We have daily fAPAR (MODIS) from Beni's workings so I can combine them here (ignoring sub-daily variations)

df_Vie_beni <- read_csv(file = "Mengoli et al 2018_input-data/2014/data_Beni2_BE-Vie_2014.csv", 
                        show_col_types = FALSE) %>% 
  mutate(date = lubridate::ymd(TIME)) %>% 
  select(c(date, fapar_spl))

## We also have half-hourly LAI (MODIS) estimates downloaded from Giulia's github repository

df_Vie_lai <- read_csv(file = "Mengoli et al 2018_input-data/2014/data_LAI_BE-Vie_2014.csv", 
                        show_col_types = FALSE) %>% 
  mutate(step = lubridate::ymd_hm(TIMESTAMP_START)) %>% 
  select(c(step, LAI))


# and we use the join argument to merge these separate tibbles

df_Vie_hh <- df_Vie_hh %>% 
  left_join(df_Vie_beni, by = "date") %>% 
  left_join(df_Vie_lai, by = "step") %>% 
  rename(fapar = fapar_spl)



# let's have a look at that - are the values reasonable?

summary(df_Vie_hh)

# no obvious problems there. 

######
## We'll start with the existing Big-Leaf representation of the P-model, but at the half-hourly time-step:
######

library(rpmodellt) # from Victor, based on the P-model with acclimation presented in Mengoli et al. 2021


### We use this run to generate big-leaf Vcmax estimates relating to the top of the canopy.  
## from Julia's workings these estimates are acclimated - as based on a weighted mean for the antecedent period (here 15 days) and optimised for noon conditions:


df_output <- as_tibble(rpmodel_subdaily(TIMESTAMP = df_Vie_hh$step, tc = df_Vie_hh$tc, 
                                          vpd = df_Vie_hh$vpd, co2 = df_Vie_hh$co2, ppfd = df_Vie_hh$ppfd, 
                                          patm =df_Vie_hh$patm, fapar = df_Vie_hh$fapar, 
                        # optional extras
                        # in particular notice the acclimation options (time-frame, averaging etc.)
                        LAI = NA,  
                        elv = 493, 
                        u=NA, # windspeed
                        ustar = NA,
                        canopy_height = NA, 
                        z = NA, leafwidth = NA,
                        netrad = NA, 
                        beta = 146.0, c_cost = 0.41,
                        do_leaftemp = FALSE,  gb_method = "simple", #Choudhury_1988"
                        do_acclimation = TRUE,
                        upscaling_method = "noon", hour_reference_T = 12,
                        acclim_days = 15, weighted_accl = TRUE, 
                        epsleaf = 0.96, #thermal absorptivity of the leaf
                        energy_params = list(
                          ste_bolz = 5.67e-8, #W m^-2 K^-4
                          cpm = 75.38, #J mol^-1 ºC-1
                          J_to_mol = 4.6, #Conversion factor from J m-2 s-1 (= W m-2) to umol (quanta) m-2 s-1
                          frac_PAR = 0.5, #Fraction of incoming solar irradiance that is PAR
                          fanir = 0.35 #Fraction of NIR absorbed
                          )
                        ))


## So we have a tibble with 19 variables generated by the P-model run.  


# and I combine that (in a very generic way) with the starting dataframe; and drop those simulated variables that we won't use subsequently

df_input <- df_Vie_hh %>% 
  cbind(df_output) %>% 
  select(-c(gpp, ca, chi, xi, mj, mc, iwue, gs, e, vpd_leaf, ns_star)) %>% 
  # I'll attach a suffix to the rpmodel generated variables to avoid any confusion
  # here I opt for 'assim' rather than 'gpp' for consistency of units
  rename(gpp_pmod = assim, gammastar_pmod = gammastar, kmm_pmod = kmm, ci_pmod = ci,
         vcmax_pmod = vcmax, jmax_pmod = jmax, vcmax25_pmod = vcmax25, rd_pmod =rd) 
  

## we now generate extinction coefficients to express the vertical gradient in photosynthetic capacity after the equation provided in Figure 10 of Lloyd et al. (2010):

df_input <- df_input %>% 
  mutate(kv_Lloyd = exp(0.00963 * vcmax_pmod - 2.43)) # Here I opt for vcmax rather than Vcmax25 


######
### Next the radiance partitioning after de Pury & Farquhar; 
## for now we ignore the clumping index term introduced in BESS
######

## we use a function provided by Victor Flo Sierra:  
# A number of parameters are provided from the deP&F paper (e.g. Tables 2 and 3):
# fa: scattering coefficient of PAR in the atmosphere
# sigma: leaf scattering coefficient of PAR (reflection and transmissivity)
# rho_cd: canopy reflection coefficient for diffuse PAR
# kd_prime: diffuse and scattered diffuse PAR extinction coefficient


library(solrad) # a CRAN package used for calculating solar radiation and related variables


irradiance_partition <- function(sun_shade = c("sun","shade"), TIMESTAMP, PPFD, LAI, PA, 
                                 PA0 = 101325, fa = 0.426, sigma = 0.15, 
                                 rho_cd = 0.036, kd_prime = 0.719, lat, long){
  
  stopifnot("TIMESTAMP must be a POSIXct, POSIXtl or a Date object" = all(lubridate::is.instant(TIMESTAMP)))
  if(any(PA <10000)) warning("Check if PA is in Pascal")
  if(any(PPFD > 3000)) warning("Check if PPFD is in micromol m-2 s-1")
  
  #Solar elevation angle
  DOY <- lubridate::yday(TIMESTAMP)
  hour <- lubridate::hour(TIMESTAMP)
  DOY_dec <- solrad::DayOfYear(TIMESTAMP)
  beta_angle <- solrad::Altitude(DOY = DOY_dec, 
                                 Lat = lat,
                                 Lon = long,
                                 SLon = long,
                                 DS = 0) * pi/180 # converting degrees to radians
  # beam extinction coefficient
  kb = 0.5/sin(beta_angle)
  # I overwrite Victor's step here, if beta is small (say 5°), then kb should be large
  #kb = ifelse(beta_angle<=0, 1e-10, kb)
  kb <- if_else(beta_angle * (180/pi) <= 1, 30, kb) # a beta of 1° would deliver a kb of >28
  # beam and scattered beam extinction coefficient
  kb_prime = 0.46/sin(beta_angle)
  #kb_prime = ifelse(beta_angle<=0, 1e-10,kb_prime )
  kb_prime <- if_else(beta_angle * (180/pi) <= 1, 27, kb_prime)
  # fraction of diffuse radiation
  m = (PA/PA0)/sin(beta_angle)
  fd = (1-0.72^m)/(1+0.72^m*(1/fa - 1))
  # beam irradiance horizontal leaves
  rho_h = (1-(1-sigma)^0.5)/(1+(1-sigma)^0.5)
  # beam irradiance uniform leaf-angle distribution
  rho_cb = (1-exp(-2*rho_h*kb/(1+kb)))
  # diffuse irradiance
  I_d = PPFD*fd
  I_d = ifelse(I_d<0,0,I_d)
  # beam irradiance
  I_b = PPFD*(1-fd)
  # scattered beam irradiance
  I_bs = I_b*((1-rho_cb)*kb_prime*exp(-kb_prime*LAI)-(1-sigma)*kb*exp(-kb*LAI))
  # Irradiance sun exposed
  I_c = (1-rho_cb)*I_b*(1-exp(-kb_prime*LAI))+(1-rho_cd)*I_d*(1-exp(-kd_prime*LAI))
  
  a = I_b*(1-sigma)*(1-exp(-kb*LAI))
  b = I_d*(1-rho_cd)*(1-exp(-(kd_prime+kb)*LAI))*kd_prime/(kd_prime+kb)
  c = I_b*(((1-rho_cb)*(1-exp(-(kb_prime+kb)*LAI))*(kb_prime/(kb_prime+kb)))-
             (1-sigma)*((1-exp(-2*kb*LAI)))/2)
  
  Isun = a+b+c
  
  # I introduce a clause here to exclude hours of obscurity (here a beta_angle > 1°, expressed in Radians)
  Ishade = if_else(beta_angle > 0.02,
                   I_c - Isun,
                   0)
  
  # don't quite understand this bit about sun_shade, replace for now with a list of the two fractions
  #if(sun_shade == "sun"){return(Isun)}else
    #if(sun_shade == "shade"){return(Ishade)}else{"sun_shade variable must be either sun or shade"}
  
  output <- list(
    Icsun = Isun,
    Icshade = Ishade,
    sEa = beta_angle, # solar elevation angle (aka solar altitude angle)
    kb = kb # beam irradiance extinction coefficient
    )
  
  return(output)
  
  }


## and we can run that partitioning function supplying inputs from our mini BE-Vie dataset


df_dePF <- as_tibble(irradiance_partition(TIMESTAMP = df_input$step, PPFD = df_input$ppfd, 
                                          PA = df_input$patm, LAI = df_input$LAI, 
                                          lat = 50.30493, long = 5.99812))


# and again, I combine that with the starting dataframe

df_input <- df_input %>% 
  cbind(df_dePF) 


#######
### Canopy photosynthesis
### So the two-leaf, two-stream model has separate calculations for the two fractions: sun and shade
######

## we derive separate photosynthetic parameters for the two groups, following BESS (or Mengoli) formulations:

df_input <- df_input %>% 
  
  # start with carboxylation
  mutate(Vmax25_canopy = LAI * vcmax25_pmod * ((1 - exp(-kv_Lloyd)) / kv_Lloyd),
         Vmax25_sun = LAI * vcmax25_pmod * ((1 - exp(-kv_Lloyd - kb * LAI))/ (kv_Lloyd + kb * LAI)),
         Vmax25_shade = Vmax25_canopy - Vmax25_sun,
         
         # next we convert those sun and shade values for ambient temperature using an Arrhenius function
         Vmax_sun = Vmax25_sun*exp(64800*(tc-25)/(298*8.314*(tc+273))),
         Vmax_shade = Vmax25_shade*exp(64800*(tc-25)/(298*8.314*(tc+273))),
         
         # now the photosynthetic estimates
         Av_sun = Vmax_sun * (ci_pmod - gammastar_pmod)/(ci_pmod + kmm_pmod),
         Av_shade = Vmax_shade * (ci_pmod - gammastar_pmod)/(ci_pmod + kmm_pmod), 
         
         ## and now separate Jmax estimates for sun and shade;
         
         Jmax25_sun = 29.1 + 1.64 * Vmax25_sun, # Eqn 31, after Wullschleger
         Jmax25_shade = 29.1 + 1.64 * Vmax25_shade,
         
         #phi0 = 0.352/8 + 0.022*tc - 0.00034*tc^2,
         
         #Jmax25_sun = (4 * phi0 * Icsun) / 
           #sqrt((1 / (1 - (0.41*(ci_pmod + 2*gammastar_pmod)/ (ci_pmod - gammastar_pmod))^(2/3))) - 1),
         #Jmax25_shade = (4 * phi0 * Icshade) / 
           #sqrt((1 / (1 - (0.41*(ci_pmod + 2*gammastar_pmod)/ (ci_pmod - gammastar_pmod))^(2/3))) - 1),
         
         
         # temperature correction (Mengoli 2021 Eqn 3b); relevant temperatures given in Kelvin
         Jmax_sun = Jmax25_sun * exp((43990/8.314) * (1/298 - 1/(tc+273))),
         Jmax_shade = Jmax25_shade * exp((43990/8.314) * (1/298 - 1/(tc+273))),
         
         # and now to calculate J and Aj for each group
         
         #J_sun = (4 * phi0 * Icsun) / sqrt(1 + ((4 * phi0 * Icsun) / Jmax_sun)^2),
         #J_shade = (4 * phi0 * Icshade) / sqrt(1 + ((4 * phi0 * Icshade) / Jmax_shade)^2),
         
         J_sun = (Jmax_sun * Icsun * (1 - 0.15) / (Icsun + 2.2 * Jmax_sun)),
         J_shade = (Jmax_shade * Icshade * (1 - 0.15) / (Icshade + 2.2 * Jmax_shade)),
         
         
         Aj_sun = (J_sun / 4) * (ci_pmod - gammastar_pmod) / (ci_pmod + 2 * gammastar_pmod),
         Aj_shade = (J_shade / 4) * (ci_pmod - gammastar_pmod) / (ci_pmod + 2 * gammastar_pmod),

         
         # gross assimilation for each class of leaves is the minimum of the two rates
         # again restricted to daylight hours (solar Elevation angle > 1°, here in Radians)
         Acanopy_sun = if_else(sEa > 0.02, 
                               if_else(Av_sun > Aj_sun, Aj_sun, Av_sun),
                               0),
         Acanopy_shade = if_else(sEa > 0.02, 
                                 if_else(Av_shade > Aj_shade, Aj_shade, Av_shade),
                                 0),
         
         # and we account for canopy respiration
         gpp_canopy = if_else(sEa > 0.02,
                              Acanopy_sun + Acanopy_shade - (rd_pmod * LAI),
                              0))

[davidorme]

@j-emberton I'm thinking that we move to a new Productivity abstract base class, with then derived class signatures something like:

class TwoLeafProductivity(ProductivityABC):

    def __init__(self, pmodel: PModel | SubdailyPModel, ppfd: NDArray, fapar: NDArray, beam_angle: NDArray):
        ...

[j-emberton]

I think this is reasonable. @AmyOctoCat and I will need to join up to discuss how to define the Productivity ABC so that it makes sense for the two productivity implementations. When I've used ABCs in the past, it's been to tie down specific methods, return types, that kind of thing to make sure all the derived classes are consistent. I assume we'll do something similar here?

[davidorme]

That's my first thought yes, but I need to talk to @kjbloom to try and figure out the details.

One critical difference at the moment is that the __init__ method for the PModel does a lot of the calculations and then there is that separate step estimate_productivity to scale up the predictions to the incoming radiation. The SubdailyModel currently does that all in one step: the behaviour of the model involves a lagged process that requires the radiation.

We may have to implement the TwoLeaf canopy separately for the two implementations (using shared internal functions). I think there may be a way to make the PModel and SubdailyModel have more similar implementations, but we'd need to look in detail at the internals.

[davidorme]

Or... it might be considerably easier to leave things as they are and have a separate post-processing class. Take an existing bigleaf model (of either PModel or SubdailyModel) and then update the predictions to a TwoLeafModel. That doesn't feel quite as elegant, but it is probably much simpler and avoids a lot of refactoring.

My feeling right now is that a lot of the ABC design will be about trying to reconcile these two specific implementations (BigLeaf/TwoLeaf). There may never be another implementation, so if (right now) we can simply make an addon converter, rather than engineering a new ABC, that might be a better place to start.

[j-emberton]

I get what you're saying and for the time being let's approach it that way.

For the future though, we might want to have a think about whether the PModel class is structured in the right way to be flexible enough for future modifications. You're in a much better position to understand that than me, but if we decide the class needs a refactor I'm happy to do it.

[davidorme]

In complete agreement on the flexibility longer term structure - and aligning PModel and SubdailyModel to make them more extensible would be an added bonus. Normally, I'm all over future proofing code (to an unhealthy degree) but my instinct here is that we don't know enough about that bigger picture and should go for the quicker win and establish some benchmarks.

[davidorme]

Just to unpack that last paragraph, we need something to make it easy for users to get solar azimuth and altitude.

• We've got a bunch of solar calculations in pyrealm.splash.solar - I don't know how much (if any) of that can be repurposed and extended but it is currently daily estimates of radiation. The internals of pyrealm.core.solar.DailySolarFluxes are a big lump of code at the moment, so we could move more into pyrealm.core.solar (a single function at the moment) and extend that core functionality. The advantage here is that all uses numpy already but may be missing too much.

• Other packages - there are a lot of options but they seem to be falling into single location solar energy implementations (https://github.com/pingswept/pysolar - also doesn't look like it has long term support) or insanely precise galactic and solar physics packages (https://www.sunpy.org/ - worried about conda install and other non-python packages; https://github.com/astropy/astropy - would allow us to apply pyrealm on any astronomical object, opening up exciting astrobiology applications, but possibly a bit over the top in terms of precision!).

[davidorme]

The pvlib.irradiance options look interesting, but the library does seem to be Location and PV Array based. I guess there is a lot that looks interesting there, but the key things for this particular feature are the zenith and azimuth - we currently typically rely on datasets providing shortwave downwelling radiation (e.g. WFD / WFDE5) which I think account for cloud effects. But it would be neat if the solution we pick here could also provide PPFD.

[j-emberton]

So, looking at pyrealm.splash.solar.DailySolarFluxes, it does two useful things. It calculates the PPDF which is useful for us, and it calculates the daily solar rad.

I propose something along these lines (based on what you've been hinting at):

• break out the physics calcs from DailySolarFluxes into standalone functions. These can then be called as required by any Class/Dataclass.

• meanwhile we create another dataclass that does basically the same stuff as pyrealm.splash.solar.DailySolarFluxes but on an arbitrary timescale i.e. for any two points in time it can calculate the solar irradiance and PPFD for that time period. We can do this by, e.g. for hourly calcs, calculating the solar hour angle then getting the solar zenith angle. It's then simple to get the solar rad for the hour. Might be simpler to start with 1hr blocks.

It wouldn't be difficult to add astronomy options in the future, and it would be very topical given the interest in space colonisation etc.

[davidorme]

I think that makes sense. That allows us to support:

• Calculating the daily estimates from the SPLASH model (pyrealm.splash.splash.SplashModel) as we currently need. The pyrealm.splash.solar.DailySolarFluxes calculates attributes that aren't used in the SPLASH model, so we could probably simplify those calculations.
• Calculate solar azimuth, zenith and fluxes for use in two leaf model (and possibly extend to the current big leaf model).
• Provide an alternative source of PPFD for use in fitting the P Model, which is a nice additional feature.

[davidorme]

I should add that the original SPLASH model implementation - which is included in the pyrealm_build_data.splash module - was written by Tyler Davis (https://bitbucket.org/labprentice/splash/src/master/). I've fairly radically refactored that code and numpy-ized it but I didn't really touch the internals of DailySolarFluxes.__init__ at all. I wanted to make it more modular but didn't fully understand where the natural function breaks were, so that would be an added bonus.

[davidorme]

I don't know the solar physics at all here, but if there are a wide range of equations for estimating downwelling radiation and solar azimuth and altitude, I think we'd be happy to sacrifice some precision for speed, at least initially. We can always add method arguments later to call out to more precise implementations.